Revision of the real-time Integrator methods

[andreamarini]

On the Yambo notes on overleaf i wrote a detailled analysis of what I know of the numerical approaches to the solution of the real-time dynamics, belonging to the Runge-Kutta family.

I added references and I detailed the math as much I could do.

Some observations follow:

As already stated in issue #95 I have demonstrated that:

• the EXP integrator (and the associated implementations in Yambo) is not part of the Runge-Kutta family and cannot be used as elemental step. From the book of Butchler (cited in the notes) the coeffiecients $a,b,c$ of the Runge-Kutta methods must respect specific internal constrains based on Taylor expantion and interpolation rules. So it cannot be replaced with something else.
• the EXP integrator could be used as standalone only if its correctness and stability is first demonstrated. Moreover I could not find references anywhere. If anyone knows a reference please comment.
• the EXP has nothing to do with Crank-Nicholson that is a second-order implicit RK method (see overleaf)

More observations about yambo implementation:

• naming seems very confusing to me. Heun of second-order is equivalent to RK2, for example. More integrators (including RK3) should be added
• the (\a) matrix in yambo is replaced by a vector. But this possible only for some RK approaches. RK3 for example requires a full matrix.
• in the iterative RK scheme the code does not update the explicit dependence of the coherent part and, indeed, the (c) factors are not used.
• as far as I can see the code assumes an iterative form of $y_n(t)$ in terms of $y_{n-1}(t)$ that is valid only in specific cases of super-simple $a$ matrices format and not, for example, in RK4 or RK3.

[sangallidavide]

Thanks Andrea for the notes.

I also point to this paper, and the Octopus code documentation
https://doi.org/10.1103/PhysRevE.96.063307
https://www.octopus-code.org/documentation/main/manual/calculations/time-dependent/

Here I just want to restrict the discussion to "RK2 + INV" of yambo (which is the integratore I would suggest as default) and "RK2 + EXP(n)"

The Crank-Nicolson (CN) method (equation 25 in the overleaf notes), is written as equation 21 (in the Phys.Rev. E, after a key approximation), and slightly differently in the Octopus documentation, e.g. with H evaluated at time-step T+dt/2 (following the "mid-point rule")

In both cases it is reduced to an explicit solver. However, in the Octopus notes, it requires to evaluate the hamiltonian at time t+dt/2, what you call "explicit-mid point) or RK22

To connect with yambo, this needs to be generalized to the case of the density matrix. This is done in equations (6) and (7) in my notes. This corresponds to the yambo implementation "RK2+INV", which, folllowing Octopus, and what done by Myrta and Claudio in yambo_nl, can be called CN (or approximated CN)
This is an approximation to the "implicit CN", since there is the approximations used in the Phys. Rev. E (right before eq. 21), or its "mid-point" refinement.

Similarly the mid-point or RK22 (simply RK2 in the yambo language) can be combined with the exponential integrator. This gives "RK2+EXP". Again, following Octopus, this is called Exponenital mid-point (EMP)

[andreamarini]

Dear Davide please point me to specific points in the overleaf notes you don't agree and/or you want to comment. And please define, coherently with the formalism introduced there, the formal derivation and definition of the RK2+EXP. I have tried to do it on a piece of paper but I cannot get it. Moreover to easily show it's application please consider the simple ODE I use at the end of the notes. There are also at least two points raised about the validity of the exp elemental evolution operator that I kindly ask you to consider.

…

On 30 August 2024 11:31:52 CEST, Davide Sangalli ***@***.***> wrote: Thanks Andrea for the notes. I also point to this paper, and the Octopus code documentation https://doi.org/10.1103/PhysRevE.96.063307 https://www.octopus-code.org/documentation/main/manual/calculations/time-dependent/ Here I just want to restrict the discussion to "RK2 + INV" of yambo (which is the integratore I would suggest as default) and "RK2 + EXP(n)" The Crank-Nicolson (CN) method (equation 25 in the overleaf notes), is written as equation 21 (in the Phys.Rev. E, after a key approximation), and slightly differently in the Octopus documentation, e.g. with H evaluated at time-step T+dt/2 (following the "mid-point rule")  In both cases it is reduced to an explicit solver. However, in the Octopus notes, it requires to evaluate the hamiltonian at time t+dt/2, what you call "explicit-mid point) or RK22 To connect with yambo, this needs to be generalized to the case of the density matrix. This is done in equations (6) and (7) in my notes. This corresponds to the yambo implementation "RK2+INV", which, folllowing Octopus, and what done by Myrta and Claudio in yambo_nl, can be called CN (or approximated CN) This is an approximation to the "implicit CN", since there is the approximations used in the Phys. Rev. E (right before eq. 21), or its "mid-point" refinement. Similarly the mid-point or RK22 (simply RK2 in the yambo language) can be combined with the exponential integrator. This gives "RK2+EXP". Again, following Octopus, this is called Exponenital mid-point (EMP)  -- Reply to this email directly or view it on GitHub: #126 (comment) You are receiving this because you were assigned. Message ID: ***@***.***>

[sangallidavide]

Dear Andrea, I do not have any specific point on which I do not agree. I just reported some extra info (a paper in particular) which extend the notes.

For the simplified model (equation (23) in your notes), this is too simple. It results in a time independent Hamiltonian, H=i*alpha where:

• the "EXP" integrator is exact, and identical to "RK2+EXP"
  Moreover, there is no difference between:
• CN, eq. 22 in your notes or eq. 20 in the Phys Rev. E
• the "simple approximated CN" or "INV" integrator, eq. 21 in the Phys. Rev. E
• the "mid point approximated CN" or "RK2+INV" integrator obtained replacing V(t) with V(t+dt/2) in eq. 21 in the Phys. Rev. E

About the exp integrator, and the time-ordering issue, check eqs. (7-9) of the Phys. Rev. E.

[andreamarini]

Dear Davide,

For the simplified model (equation (23) in your notes), this is too simple. It results in a time independent Hamiltonian, H=i*alpha where:

ok. Consider, then, the following ODE

$$\frac{d}{dt} y(t) =I\left[y(t)\right](t) $$

and write explicitly the math the leads of the RK2+EXP solution method pointing to the steps and approximations you used in Yambo.

Please use the overleaf using the notation I introduced.

Thanks

[andreamarini]

Hi Davide, I had a quick look at your overleaf notes and I have a basic comment.

If you agree about my notes this means that you agree that a general Runge-Kutta procedure is mathematically defined by equations (10) and (11).

The equations together with the definition of the corresponding Butcher's tables define inequitably the procedure.

This means that RK2, for example, corresponds to Eq.(15) and CN to Eq.(22).

So I don't understand how, if INV is Eq. (36) and thus corresponds to the Butcher's table defined in Eq. (22), you can combine it with RK2 that corresponds to a different table.

It would be helpful if you use, in your derivations, the formalism I introduced in Sec. IB (if you agree of course).

I have more comment but first I would like to define a common mathematical background.

[andreamarini]

The new implementation in phys-rt_integrators is ready.

I need still to finalize few aspects but the results are really promising.

I have added math and references on the overleaf notes. All can be reproduced with the h-BN test of the test-suite.

In short:

• now all RK (n<4) methods are implemented using a modular structure
• a simplifed Runge–Kutta–Fehlberg method is used to, given the field, estimate time-step and best method
• all time-evolution is done in the rotating frame that, now, it is not an approximation anymore.

Now a delicate aspect. From the tests I did the EXP solver does not perform well compared to the new implementation. On the basis of the fact that, to my understanding, lacks of a solid theoretical justification I propose to remove it.

The results I am referring to are on overleaf (here a copy)

I therefore propose to remove EXP integrators and related approximations (INV, LINEAR REGIME...). This will make the code much more simple and light.

[andreamarini]

The (simplified) Runge–Kutta–Fehlberg can be activated by using

Integrator= "best"           # [RT] Time-evolution RK integrator: "EULER/EMP/RK2/RK3/H3/RK42/RK41/best"
RTstep=1          as    # [RT] Real Time step: >0 absolute, <0 used as seed for the Runge-Kutta optimization procedure.
RKquestTresh= ERROR           # [RT][o/oo] Treshold used to identify the best dT/integrator

In the case above all integrators are tested starting from a base of dT=1as.

Also a specific integrator can be provided. The RKF method is activated when RTstep<0

Attached the input files and scripts I used in the h-BN case.

test.tar.gz

[andreamarini]

Updated figure using RK2+INV

[sangallidavide]

Are you comparing the following ?

• RK4 with 1 as time step
• EXP with 70 as time step
• RK2+INV with 70 as time step

[andreamarini]

Are you comparing the following ?

• RK4 with 1 as time step
• EXP with 70 as time step
• RK2+INV with 70 as time step

Yes.

[sangallidavide]

Can you put RK4 with 70 as as well ?

[andreamarini]

Can you put RK4 with 70 as as well ?

Why if the devel-rt-integrators EMP (N=2) with dT=70as works perfectly?

I am not comparing integrators in the phys-master but integrators of phy-master with phys-rt-integrators.

Moreover all math is on overleaf. Including the dT estimation.

[sangallidavide]

What is EMP ?

[andreamarini]

What is EMP ?

[sangallidavide]

So EMP is what we called RK2 so far.

So you are saying RK2 with 70 as is on top of RK4 with 1 as.
Question: is EMP combined with RWA, while the others are not?

[andreamarini]

So EMP is what we called RK2 so far.

No. On overleaf RK21 is RK2. RK22 is EMP.

In the new branch there is RK2 AND EMP.

So you are saying RK2 with 70 as is on top of RK4 with 1 as. Question: is EMP combined with RWA, while the others are not?

I am saying that in the new branch EMP (70 as) is on top of RK4 (1 as). Both run in the rotating frame.

RW is not an approximation in the new branch. It cannot be switched off.

[sangallidavide]

We should compare integrators in the same conditions.

If, in the new branch, the rotating frame is always used, than EMP means RK22 + rotating frame.
It should be compared with RK2+INV+RWA from the old branch. Otherwise we are comparing different things.

[andreamarini]

Davide. The situation I think it is really clear.

I say that RK2+INV is not a RK method and I also propose a new code that, to me, performs better.

You say that this is not clear and that RK2+INV is a formally well defined integrator.

Fine.

Math is on overleaf. RKn procedure is mathematically there, appendix A. We decided that code should be supported with math. So go ahead and support your RK2+INV with a math that demonstrates that it is a RK method. I asked you several times how can you use two Butcher tables at the same time. As far as I know you never answered.

Still you can say: RK2+INV in the rotating frame is better than EMP! Great. As I wrote several times I have no idea how to implement as I don't see how to define it mathematically. Also because if RK2+INV is ok you could also do RKn+INV. How?

Please. Let's not flood each other with tons of comments. At this moment we need math. Once we have it we discuss it and, eventually, move to the code.

thx

[sangallidavide]

So, we do not agree on the math. Fine.

Can we just do simple numerical tests under well defined conditions ?