Integration-scheme.md

+# Time dependent equations
+The `td-wslda` code solves equations of the form:
+```math
+i \dot{\varphi_n} = \hat{H}(\{\varphi_n\},t)\varphi_n
+``` 
+where the Hamiltonian is given by
+```math
+\hat{H} = 
+  \begin{pmatrix}
+    h_{\uparrow}(\bm{r},t) & \Delta(\bm{r},t) \\
+    \Delta^*(\bm{r},t)& -h_{\downarrow}^*(\bm{r},t) 
+  \end{pmatrix}
+```
+and $`\varphi_n`$ is shorthand notation for quasiparticle wavefunctions (`qpwfs`):
+```math
+\varphi_n(\bm{r},t)=
+\begin{pmatrix}
+    u_{n,\uparrow}(\bm{r},t) \\ 
+    v_{n,\downarrow}(\bm{r},t)
+  \end{pmatrix}
+```
+Note that all `qpwfs` fluctuate in time with a typical oscillating factor $`\exp(-i E_n t)`$, where $`E_n`$ is a quasi-particle energy. Since we evolve all states with quasi-particle energies from  
+the interval $`[-E_c,E_c]`$, the evolved states exhibit both very slow as well as very rapid oscillations in time.
+In order to improve accuracy of the time integration, in each integration step we subtract from the quasi-particle solution the typical frequency of oscillations, which in practice modifies 
+the equation to the form:  
+```math
+i \dot{\varphi_n} = (\hat{H}(\{\varphi_n\},t)-\langle H\rangle_n)\varphi_n
+```
+where $`\langle H\rangle_n`$ is ``instantaneous'' quasi-particle energy. Although this step introduces a significant numeric cost to the calculations, it greatly improves the accuracy of the time integration. 
+
+# ABM integration algorithm
+We convert time-dependent equations into from:
+```math
+ \dfrac{dy(t)}{dt} = f(y,t)
+```
+where:
+```math
+y(t) \equiv \varphi_n(\bm{r},t)
+```
+```math
+f(y,t) \equiv \dfrac{1}{i}(\hat{H}(\{\varphi_n\},t)-\langle H\rangle_n)\varphi_n(\bm{r},t)
+```
+
+Adams-Bashforth-Moulton (ABM) is designed as follow:
+* predictor of 4-th order:
+```math
+y_k^{(p)} = y_{k-1} + \dfrac{55}{24}\Delta t f_{k-1} - \dfrac{59}{24}\Delta t f_{k-2} + \dfrac{37}{24}\Delta t f_{k-3} - \dfrac{9}{24}\Delta t f_{k-4}
+``` 
+* corrector of 5-th order:
+```math
+y_k = y_{k-1} + \dfrac{251}{720}\Delta t f(y_k^{(p)},k\Delta t) +\dfrac{646}{720}\Delta t f_{k-1}-\dfrac{264}{720}\Delta t f_{k-2}+\dfrac{106}{720}\Delta t f_{k-3} - \dfrac{19}{720}\Delta t f_{k-4}
+```
+where $`f_{k}=f(y_k, k\Delta t)`$. The accuracy of the method is set by the corrector.  
+
+The algorithm requires 5 (wave function) buffers called: $`y_{k-1}`$, $`f_{k-1}`$, $`f_{k-2}`$, $`f_{k-3}`$, $`f_{k-4}`$.  
+
+The algorithm is as follow (operations listed within a single step have to be executed simultaneously):  
+
+1. 
+```math
+  \left\lbrace 
+  \begin{array}{lll}
+   y_{k-1} & \leftarrow & y_{k-1} + \dfrac{55}{24}\Delta t f_{k-1} - \dfrac{59}{24}\Delta t f_{k-2} + \dfrac{37}{24}\Delta t f_{k-3} - \dfrac{9}{24}\Delta t f_{k-4},\\
+   f_{k-4} & \leftarrow & y_{k-1} +\dfrac{646}{720}\Delta t f_{k-1}-\dfrac{264}{720}\Delta t f_{k-2}+\dfrac{106}{720}\Delta t f_{k-3} - \dfrac{19}{720}\Delta t f_{k-4}
+  \end{array}
+  \right. 
+```
+2. Compute densities from $`y_{k-1}`$ and formulate potentials (like $`U(\bm{r})`$ and $`\Delta(\bm{r})`$) needed to formulate Hamiltonian. This step requires MPI communication. For densities and potentials use additional buffers, note the are small in comparison to wave function buffers.
+
+3. $`y_{k-1} \leftarrow f(y_{k-1},k\Delta t)`$, this is first application of Hamiltonian, as Hamiltonian use potentials formulated in step 2. This step can be done is batched mode. This is a numerically intensive step.
+
+4. $`y_{k-1} \leftarrow \dfrac{251}{720}\Delta t y_{k-1} + f_{k-4}`$, this step can be merged with step 3. Now buffer $`y_{k-1}`$ holds wave functions for $`t=k\Delta t`$.
+
+5. Compute densities from $`y_{k-1}`$ and formulate potentials (like $`U(\bm{r})`$ and $`\Delta(\bm{r})`$) needed to formulate Hamiltonian. This step requires MPI communication. For densities and potentials use additional buffers, note the are small in comparison to wave function buffers.
+
+6.  Preparation of buffers for next step. Note, the last operation corresponds to second application of Hamiltonian, as Hamiltonian use potentials formulated in step 5. This step can be done is batched mode. This is numerically intensive step.
+```math
+  \left\lbrace 
+  \begin{array}{lll}
+   f_{k-4} & \leftarrow & f_{k-3},\\
+   f_{k-3} & \leftarrow & f_{k-2},\\
+   f_{k-2} & \leftarrow & f_{k-1},\\
+   f_{k-1} & \leftarrow & f(y_{k-1},k\Delta t).
+  \end{array}
+  \right. 
+```
+
+7. Increase $`k`$ by one and go to step 1.
\ No newline at end of file