yambo-code / yambo

This is the official GPL repository of the yambo code
http://www.yambo-code.eu/
GNU General Public License v2.0
99 stars 39 forks source link

INV and RWA integrators notes #95

Closed andreamarini closed 1 month ago

andreamarini commented 5 months ago

Again on the RT_integrators.

Can you please share here the math underlying the INV and RWA integrators?

I can't get the INV math from the suborutines.

The RWA puzzles me as it should not be an approximation but a rotating frame representation (RWF) that is exact as it eliminates the independent particle oscillations.

Actually it is based on the popular solution of a system of first order linear ODEs using diagonalization of the linear coefficient.

https://www.math.uh.edu/~dblecher/EMo.pdf

In the case of G^< the matrix to diagonalize is already diagonal.

sangallidavide commented 5 months ago

Can you please share here the math underlying the INV and RWA integrators?

Andrea, from what I remember, you coded the RWA. :-)

For the INV it comes from yambo_nl. I also attach here the notes sent via email. integrators.pdf

The RWA puzzles me as it should not be an approximation but a rotating frame representation (RWF) that is exact as it eliminates the independent particle oscillations.

Yeap, fully agree. I already commented in the past that it should not be called RWA, since it is not an approximation

andreamarini commented 2 months ago

I am working on the time integrators.

I read in integrators.pdf that the linarization (at the basis of the PERT approach) is exact in HF

image

Can't undestand what you mean with linearization. In what? Time or $$\rho(t)-\rho(t_{ref})$$?

sangallidavide commented 2 months ago

Linearization in the dependency of h on rho: h[rho] = M*rho. This is what we do when we define the collisions M. It also implies that h[rho(t)]-h[rho^eq] = M*(rho(t)-rho^eq). This makes the implementation "numerically sensitive"

andreamarini commented 2 months ago

Based on the math below the INV/EXP are approximations valid ONLY in the slow time evolution regime.

The key passage is the splitting of the evolution operator.

If my notes are correct the validity of these approximations (that assume Eq (10)) are not based on linear response but on the assumption of small (metallic like) frequencies dictating the density matrix oscillations.

This the only physical argument it looks ok for me to motivate

image

scan.pdf

sangallidavide commented 2 months ago

I think we already discussed this.

Yeah, the EXP integrator can be proved to be good in case the time evolution of the hamiltonian (not the dynamics in general) is smooth. In yambo_rt the fast dynamics is dictated by the coherent part, which is dictated by the hamiltonian (Not sure what is "I" in the attached notes). Notice that this is exact at the IP level, and indeed the EXP reduces to what is called "RWA" since the hamiltonian is time independent.

Within other approximations, such as HF or HSEX, the hamiltonian acquires fast oscillations. However, this does not prove that in such cases the EXP integrator is worse then EULER.

Indeed the EXP and INV integrators have another important property, e.g. they are unitary. The INV (combined with RK2, e.g. INV+RK2, it becomes the popular CRANK-NICKOLSON, thera are more reasons for which this is popular) has been shown to be numerically more stable that the EULER (or EULER+RK2). To test this, just take a td-hsex simulation in the test-suite and try to perform it with different time-steps with both EULER and INV. You will see that INV remains stable and closer to the correct solution for a larger time-step.

andreamarini commented 2 months ago

Yeah, the EXP integrator can be proved to be good in case the time evolution of the hamiltonian (not the dynamics in general) is smooth.

In my notes I demonstrate that the commutator above is a necessary condition.

Please provide alternative math if you don't agree.

(Not sure what is "I" in the attached notes).

image

Notice that this is exact at the IP level, and indeed the EXP reduces to what is called "RWA" since the hamiltonian is time independent.

In the time-independent case the commutator is trivially zero.

Attached complete math. Again if you don't agree please provide math. THX

scan.pdf

sangallidavide commented 2 months ago

Yeah, at the IP level and when the external field is off, the I_k(t) is time independent.

What you prove is a "sufficient condition" for EXP being better than EULER.

When I_k(t) is time oscillating we do not know from the math which is better between EULER and EXP(/INV).

andreamarini commented 2 months ago

In the meeting we said that all coding should be supported by math.

At the moment my math demonstrates that it is possible to split the evolution operator only and only if the above commutator is zero.

Without further math this is what we have and what the code should be based on

andreamarini commented 2 months ago

yambo-tests/TESTS/MAIN/hBN/RT

Script attached

./do.tcsh -e 100 -i RK2 -l none
./do.tcsh -e 100 -i EULER -l EXP
./do.tcsh -e 100 -i EULER -l EXPACC

do.tcsh.txt

Results (y axis polarization)

image

sangallidavide commented 2 months ago

EXP and INV are alternative to EULER. All three can be combined with "RKn" or "HEUN"

The combinations EULER+EXP or EULER+EXPACC does not make sense. Not sure what was selected in input.

Instead, "EULER+RK2" (which is selected by default just specifying RK" in input) is likely more precised than simple "EXP", being a multi step approach. One should compare

See the comments here:

I remember part of the discussion was already done here: https://github.com/yambo-code/yambo/issues/93

andreamarini commented 2 months ago

I just run the master where the input file Integrator combinations are fixed in RT_Integrators_init regardless of what the user defines.

The output above is what a normal user would get.

Let me write it again: my math (2 pages of simple derivations) demonstrates that the EXP/INV are based on the slow dynamics ansatz.

I maybe wrong but we should base our discussions on math.

Base.in (needed to run the tests above) attached.

base.in.txt