I'm working on a Real-Time TDDFT implementation and I am currently comparing different propagation schemes for the propagation of the Kohn-Sham wave function, $$ \phi(t+\Delta t) = \hat{\mathcal{U}}\phi(t). $$ An external field is incorporated by a vector potential and the kinetic energy matrix element is modified such that the hamilton matrix gets a complex part: $$ \mathcal{H}_{ij} = \mathcal{H}_{ij}^{(R)} + i \mathcal{Z}_{ij} $$ where $\mathcal{H}^{(R)}$ denotes all the real parts in the hamilton matrix and $\mathcal{Z}$ is a matrix that incorporates the imaginary part of the terms arising due to the vector potential. I think details are not that important here. Among other things, I use the instantaneous total energy $E_\text{tot}(t)$ to compare different propagators$^*$.

Most implemented integration schemes are exponential-based, i.e. $$ \mathcal{U}=\exp\left(f(\mathcal{H})\right), $$ e.g. Exponential Midpoint with $f(\mathcal{H})=-i\mathcal{S}^{-1}\mathcal{H}(t+\Delta t/2)$ or Enforced Time-Reversal Symmetry with $f(\mathcal{H})=-i\Delta t/2(\mathcal{S}^{-1}\mathcal{H}(t+\Delta t)+\mathcal{S}^{-1}\mathcal{H}(t))$ ($\mathcal{S}$ is the overlap matrix). These methods mostly compute the same results but when I use the Crank-Nicolson propagator, $$ \mathcal{U}=\frac{\mathcal{S}-i\mathcal{H}(t)\Delta t/2}{\mathcal{S}+i\mathcal{H}(t)\Delta t/2}, $$ it only reproduces the results when $\mathcal{Z}=0$, else I get (a while after the pulse ends, constant) amplitude oscillations that even persist after a pulse is switched off. Reducing the step size does not change anything.

enter image description here

So my question is what the reason for this could be. I don't know too much about details of applying these propagation schemes to the given kind of problem.

$^*$ $E_\text{tot}(t) = \sum_i f_i \epsilon_i(t) - F[\rho(t)] + E_{nn}$ where $\epsilon_i$ are the eigenvalues of the time-dependent hamiltonian, $F$ denotes some functionals incorporating XC and electrostatic corrections and $E_{nn}$ is the internuclear repulsion energy.

  • 2
    $\begingroup$ This is likely due to L-stability. People think of the Crank-Nicholson scheme as "good for stiff" because it's A-stable, but that's a very weak form of stability. In reality, the production-quality ODE solvers for stiff equations all have L-stability, which is the dampening of oscillatory modes (at infinity). This might be an example of where L-stability is required instead of just A-stability. $\endgroup$ – Chris Rackauckas Nov 24 '17 at 14:14
  • $\begingroup$ Thanks alot for your hint! I will have a look at that. Could you maybe recommend literature for this topic? $\endgroup$ – Lukk Dec 6 '17 at 9:06
  • $\begingroup$ Hairer's second book on ODE methods (on stiff solvers) has a very good treatment of this topic. $\endgroup$ – Chris Rackauckas Dec 6 '17 at 15:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.