3.30 Real-Time and Imaginary-Time TDDFT

Real-Time TDDFT - as the name suggests - calculates the response of a Kohn-Sham system in real time. The time-dependent Kohn-Sham equation

$$i\frac{\partial }{\partial t}{\psi }^{\mathrm{KS}}(𝐫,t)={\mathscr{H}}^{\mathrm{KS}}[\rho (t),t]{\psi }^{\mathrm{KS}}(𝐫,t)$$

(3.110)

is solved and describes the time-evolution of the electronic wave function ${\psi }^{\mathrm{KS}}(𝐫,t)$. The Hamiltonian has an explicit time-dependence in this case, e.g. by an external electric field and/or a time-dependent external ionic potential.
In contrast to the linear response approach, real-time TDDFT is able to capture the whole nonlinear characteristic of the given system in real time. Possible applications include absorption spectrum calculations based on the dynamical dipole response, high-harmonic generation simulations or ion bombardment simulations based on non-adiabatically coupled electron-ion dynamics (Ehrenfest dynamics). Please see [212] for a comprehensive review.

What is actually done is the time propagation of the atomic basis function coefficients $\{c_{in}(t)\}$ from the LCAO ansatz

$${\psi }_n^{\mathrm{KS}}(𝐫,t)=\sum _i^{N_{\mathrm{basis}}}{c}_{in}(t){\varphi }_i(𝐫-{𝐑}_{I(i)}),n\in N^{\mathrm{occ}}.$$

The coefficients are expressed as a matrix by considering all (occupied) states and basis indices: $𝐂\in {ℂ}^{N_{\mathrm{basis}}\times N_{\mathrm{occ}}}$. The time-dependent KS matrix equation to be solved is then

$$\frac{d}{dt}𝐂(t)=-i{𝐒}^{-1}𝐇(t)𝐂(t)$$

where $𝐒\leftrightarrow ⟨{\varphi }_i|{\varphi }_j⟩$ is the overlap matrix and $𝐇\leftrightarrow ⟨{\varphi }_i|{\mathscr{H}}^{\mathrm{KS}}|{\varphi }_j⟩$ is the Hamiltonian matrix. The efficient and accurate solution of this equation is the key mechanic in a real-time TDDFT implementation. Note that very large basis sets will lead to an ill-conditioned overlap matrix, causing ${𝐒}^{-1}$ to be a problematic object, yielding entirely wrong results or blow-up.
The time-dependence of the Hamiltonian is implicit via the time-dependent density and explicit by the possible dependence on an external field, e.g. a laser or a dynamical ionic potential. An external electric field can here be incorporated via

	$\mathrm{Length}\mathrm{gauge}:$	$H_{ij}(t)=H_{ij}^{\mathrm{KS}}[\rho (t)]+𝐄(t)\cdot ⟨{\varphi }_i|𝐫|{\varphi }_j⟩$
	$\mathrm{Velocity}\mathrm{gauge}:$	$H_{ij}(t)=H_{ij}^{\mathrm{KS}}[\rho (t)]-i𝐀(t)\cdot ⟨{\varphi }_i|\nabla |{\varphi }_j⟩+{\displaystyle \frac{1}{2}}{𝐀}^2(t)S_{ij}$

where the dipole approximation – neglecting any spatial dependence of the electric field – is used. Electric field $𝐄$ and vector potential $𝐀$ are connected by the relation $𝐄=-{\partial }_t𝐀$. Both gauges are physically equivalent but have different technical implications. In this implementation, only the velocity gauge can be applied in case of periodic systems. We note here that the (atomic ZORA) relativistic kinetic operator is currently not considered in the expression for the velocity gauge but that is on the to-do list (it can nevertheless be used). Scaled ZORA can not be applied.

In Ehrenfest dynamics, non-adiabatically coupled motion between the electronic and ionic subsystems is simulated. Mobile atoms and thus mobile basis functions lead to a modified differential equation for the electrons,

	${\displaystyle \frac{d}{dt}}𝐂(t)$	$=-i{𝐒}^{-1}(𝐇+𝐆)𝐂(t),$
	$G_{ij}$	$=i⟨{\varphi }_i|{\dot{𝐑}}_{I(j)}\cdot \nabla |{\varphi }_j⟩,$

where the matrix $𝐆$ describes the time-dependence of the atomic basis functions and effectively conserves the norm of the time-propagated wavefunctions.
Furthermore, the forces on the nuclei are complemented by specific non-adiabatic contributions,

$${𝐅}_I={𝐅}_I^{\mathrm{HF}}+{𝐅}_I^{\mathrm{MP}}+{𝐅}_I^{\mathrm{𝐗𝐂}}+{𝐅}_I^{\mathrm{DBC}}+{𝐅}_I^{\mathrm{NC}},$$

where ${𝐅}_I^{\mathrm{HF}}$ are the standard Hellmann-Feynman forces, ${𝐅}_I^{\mathrm{MP}}$ are multipole force contributions and ${𝐅}_I^{\mathrm{XC}}$ are XC related forces, e.g. GGA. The force term

$${𝐅}_I^{\mathrm{DBC}}=-\sum _n^{N_{\mathrm{occ}}}\sum _{ij}^{N_{\mathrm{basis}}}{f}_n{c}_{in}^{\ast }{c}_{jn}\left[⟨{\nabla }_I{\varphi }_i|{\mathscr{H}}^{\mathrm{KS}}|{\varphi }_j⟩+⟨{\varphi }_i|{\mathscr{H}}^{\mathrm{KS}}|{\nabla }_I{\varphi }_j⟩\right]$$

stands for "Dynamical Basis Correction" and can be seen as the time-dependent analogue to Pulay forces. The last term is denoted as the nonadiabatic coupling force, possibly consisting of two terms:

$${𝐅}_I^{\mathrm{NC}}={𝐅}_I^{\mathrm{EC}}+{𝐅}_I^{\mathrm{𝐌𝐂}}.$$

Including the full above-mentioned expression guarantees full energy- and momentum conservation, but usually only incorporating ${𝐅}_I^{\mathrm{EC}}$ is sufficient and noticeably cheaper. The energy-conserving force term can be written as

	${𝐅}_I^{\mathrm{EC}}$	$={\displaystyle \sum _n^{N_{\mathrm{occ}}}}{f}_n{𝐜}_n^{†}\left[{\mathrm{𝐇𝐒}}^{-1}{𝐁}_I+{\left({\mathrm{𝐇𝐒}}^{-1}{𝐁}_I\right)}^{†}\right]{𝐜}_n,$
	${𝐁}_{ij,I}$	$=⟨{\varphi }_i|{\nabla }_I|{\varphi }_j⟩.$

The remaining force term restoring momentum conservation is given as

	${𝐅}_I^{\mathrm{MC}}$	$=i{\displaystyle \sum _n^{N_{\mathrm{occ}}}}{f}_n{𝐜}_n^{†}\left[{𝐖}_I^{†}-{𝐖}_I+{𝐃}^{†}{𝐒}^{-1}{𝐁}_I-{\left({𝐃}^{†}{𝐒}^{-1}{𝐁}_I\right)}^{†}\right]{𝐜}_n,$
	${𝐖}_{ij,I}$	$={\dot{𝐑}}_I⟨{\nabla }_I{\varphi }_i|{\nabla }_I{\varphi }_j⟩,$
	$D_{ij,I}$	$={\dot{𝐑}}_I\cdot ⟨{\varphi }_i|{\nabla }_I|{\varphi }_j⟩.$

The integration of the nuclear equations of motion is performed via a Velocity Verlet based algorithm, i.e.,

	${\dot{𝐑}}_I(t+\mathrm{\Delta }t/2)$	$=\dot{𝐑}(t)+{\displaystyle \frac{\mathrm{\Delta }t}{2M_I}}{𝐅}_I(t),$
	${𝐑}_I(t+\mathrm{\Delta }t)$	$={𝐑}_I(t)+\mathrm{\Delta }t{\dot{𝐑}}_I(t+\mathrm{\Delta }t/2),$

where $\mathrm{\Delta }t$ is the corresponding Ehrenfest time stepping. Based on each coordinate change, a geometry update is performed. Please see [183, 235] for further details.

Depending on what task has to be performed, different physical observables are of interest in a RT-TDDFT simulation. The computation of the most important ones is briefly explained here. These are also computed and sent to output by default.

• •

  Electronic energy: The time-dependent electronic energy for the current time step is computed via

  	$E_{\mathrm{el}}^{\mathrm{KS}}(t)=$	${\displaystyle \sum _{n=1}^{{\mathrm{N}}_{\mathrm{occ}}}}{f}_n⟨{\psi }_n(t)|{\mathscr{H}}^{\mathrm{KS}}(t)|{\psi }_n(t)⟩-{\displaystyle \int {\mathrm{d}}^3𝐫\rho (𝐫,t)v_{\mathrm{xc}}^{\mathrm{ad}}(𝐫)}$
  		$+E_{\mathrm{xc}}^{\mathrm{ad}}[\rho ](t)-{\displaystyle \frac{1}{2}}{\displaystyle \int {\mathrm{d}}^3𝐫\rho (𝐫,t)v_{\mathrm{H}}(𝐫,t)}+E_{\mathrm{nuc}}(t).$

  where the adiabatic approximation is used, i.e., the instantaneous density $\rho (𝐫,t)$ is employed to evaluate the functionals (this is indicated by the ad superscript). The electronic energy is not necessarily a conserved quantity, e.g. when energy is absorbed from an external field.

• •

  Total electronic and nuclear kinetic energy: In case when Ehrenfest dynamics is performed, the nuclear kinetic energy contribution to the total energy has to be taken into account, i.e.

  	$E_{\mathrm{tot}}(t)$	$=E_{\mathrm{el}}^{\mathrm{KS}}(t)+E_{\mathrm{nuc},\mathrm{kin}}(t)$
  		$=E_{\mathrm{el}}^{\mathrm{KS}}(t)+{\displaystyle \frac{1}{2}}{\displaystyle \sum _I}{M}_I{\dot{𝐑}}_I^2(t).$

  In the absence of external fields, this quantity must be conserved.

• •

  Electronic dipole moment: Important quantities that can be computed from the dynamical electronic dipole response are the polarisability $\alpha (\omega )$ and absorption strength $S(\omega )$:

  $$𝝁(t)=\int d^{3}r𝐫\rho (t)\to {\alpha }_{ij}(\omega )=\frac{ℱ[{\mu }_i](\omega )}{ℱ[E_j](\omega )}\to S(\omega )=\frac{2\omega }{3\pi }\mathrm{Im}\left\{\mathrm{Tr}\left[\alpha (\omega )\right]\right\}$$

  The electronic dipole moment is by default computed for cluster systems. One can also compute tranisition dipole moments which are obtained by projection onto the ground state Kohn-Sham orbitals [46].

• •

  Electronic current: From the current density $\text{j}(\text{r},t)$, one can compute the conductivity $\sigma (\omega )$ and the dielectric function $ϵ(\omega )$ via the total current $\text{I}(t)$:

  	$\text{j}(\text{r},t)$	$={\displaystyle \frac{1}{2}}{\displaystyle \sum _{nij}}{f}_n[c_{in}^{\ast }(t)c_{jn}(t)⟨{\varphi }_i|-i\nabla +\text{A}(t)|{\varphi }_j⟩+\mathrm{c}.\mathrm{c}.]\to \text{I}(t)=-{\displaystyle \frac{1}{V}}{\displaystyle \int }{d}^{3}r\text{j}(\text{r},t)$
  	$\to$	${\sigma }_{ij}(\omega )={\displaystyle \frac{ℱ[I_i](\omega )}{ℱ[E_j](\omega )}}\to ϵ(\omega )=1+i4\pi {\displaystyle \frac{\sigma (\omega )}{\omega }}$

• •

  Trajectories: In case of Ehrenfest dynamics, trajectories are written to .xyz file/s.

Please note that a python-based post-processing tool named eval_tddft.py is located in the utilities/rt-tddft folder which is made for RT-TDDFT output files. It can currently compute, visualize and write the absorption strength function/polarisability (by input electronic dipole and external field) and also the conductivity/dielectric function (by input electronic current and external field) with possibly more functionality to come.

Some general important remarks regarding different topics are given here. Please also pay attention to this.

• •

  Exchange-Correlation: Due to technical and/or formal limitations, not all available XC functionals can be used. We tested for LDA- and GGA-based functionals which should work fine. Other choices were not tested or are in development process and the code will stop for unreasonable settings. Hybrid functionals like B3LYP can be used for finite systems. This also applies for the long-range corrected LC-$\omega$PBE functional.

• •

  Compatibility: We tested our software for Intel-based parallelization together with the MKL numerical linear algebra library, and for GNU OpenMPI-based parallelization together with the numerical linear algebra libraries ScaLAPACK/OpenBLAS (versions 2019 / 2020). It seems like there exists a bug in the recent MKL versions (Intel Parallel Studio XE 2019 / 2020) that affects the linear equation solver subroutines pzgetrs and pzgesv, causing memory leakage at every call, resulting in excessive memory usage by the code. This was not observed with the non-Intel framework. These operations are conducted when the Crank-Nicolson propagator is used or when the scaling and squaring method for matrix exponentials is employed. We implemented a fix that will result in the code using explicit inversion instead of solving a linear system, i.e.:

  	$\mathrm{𝐋𝐢𝐧𝐞𝐚𝐫}\mathrm{𝐬𝐨𝐥𝐯𝐞𝐫}:\mathrm{solve}\mathrm{𝐀𝐗}=𝐁\mathrm{for}𝐗$		(3.111)
  	$\mathrm{𝐈𝐧𝐯𝐞𝐫𝐬𝐢𝐨𝐧}:\mathrm{compute}{𝐀}^{-1}\mathrm{for}𝐗={𝐀}^{-1}𝐁$		(3.112)

  For this, the user has to provide the flag

  • RT_TDDFT_crank_nicolson_solve_inv .true.

  in the control.in file.

Quite ’typical’ values for time steppings are 0.1 a.u. $=0.0024$ fs or for electric field amplitudes 0.01 a.u. $=51.42$ eV/Å  (in the perturbative regime).
This keyword must be set by the user or the code will not run which is intended to prevent unit confusions.

The choice of the time step heavily affects the accuracy and stability of the time-propagation and is restricted by internal or external degrees of freedom, e.g., the spectral range of the Hamiltonian or the frequency of an external field oscillation. Generally, a small time step yields better results but significantly increases the computation time. The possible time step also depends on the propagation scheme which is explained later. A conservative choice is dt_prop $\le 0.1\mathrm{a}.\mathrm{u}.(=0.0024\mathrm{fs}$) but $0.2$ a.u. to $0.5$ a.u. might also be possible. Remember using consistent units set by the RT_TDDFT_input_units keyword. Note that the time step is also a convergence parameter for spectral information.

While this makes no difference formally, it has some technical implications. When doing periodic simulations, only the velocity gauge can be used. In the cluster case, the length gauge should be used because it is computationally cheaper. However, the formal gauge invariance is broken in any numerical approximation and it is likely possible that basis set and/or real-space grid settings converge differently for length and velocity gauge.
Note that the gauge is a universal setting for each individual RT-TDDFT simulation and applies to all fields specified.

The choice of the external field depends on the type of calculation that one wants to perform. To obtain an absorption spectrum via the dynamical electronic dipole response, a delta-kick field can be applied by, e.g., setting type = 2, center, width and Ez to some values which will apply a Gaussian pulse of given amplitude at center time/width and along the z-axis. width should then be chosen quite small, e.g., to 0.05 fs. Alternatively, one can provide a constant field for the first time step, which saves computation time and provides a sharp pulse naturally. This would be type = 1 and t_end = dt.
To apply a field pulse with a defined wavelength, one can, e.g., set type = 3, freq, cycle, Ez which will result in a modulated Gaussian pulse of given wavelength around center time.
Multiple instances of this keyword with different settings can be given, e.g., to apply multiple different pulses over time. Internally, all defined fields are simply added up, so pay attention to temporal localization (if intended).
For all field types, the length and velocity gauge implementation are related by $𝐄(t)=-\frac{d𝐀(t)}{dt}$.
Note that the field amplitudes are always given as the amplitudes of the electric field, even in the case when the velocity gauge is used (see RT_TDDFT_td_field_gauge).
Currently, the following field types are implemented:

• •

  0: No field at all. This is the default when RT_TDDFT_td_field is not specified.

• •

  1: Constant Electric field.

  	$𝐄(t)$	$={𝐄}_0$
  	$𝐀(t)$	$=-{𝐄}_0(t-t_0)$

• •

  2: Delta kick pulse via Gaussian:

  	$𝐄(t)$	$={𝐄}_0\exp \left({\displaystyle \frac{t-t_c}{t_w}}\right){\left(1+\exp \left({\displaystyle \frac{t-t_c}{t_w}}\right)\right)}^{-2}$
  	$𝐀(t)$	$=-{𝐄}_0{t}_w\left(1-{\left(\exp \left({\displaystyle \frac{t-t_c}{t_w}}\right)\right)}^{-1}\right)$

• •

  3: Localized field oscillation via modulated Gaussian:

  	$𝐄(t)$	$={𝐄}_0\exp \left(-{\displaystyle \frac{{(t-t_c)}^2}{2w^2}}\right)\sin \left(\omega (t-t_c)\right)$
  	$𝐀(t)$	$={\displaystyle \frac{{𝐄}_0}{\omega }}\exp \left(-{\displaystyle \frac{{(t-t_c)}^2}{2w^2}}\right)\cos \left(\omega (t-t_c)\right)$

  with $w=t_w/(2\sqrt{2\log (2)})$, $t_w=\mathrm{FWHM}$.

• •

  4: Sinusoidal field with reference frequency:

  	$𝐄(t)$	$={𝐄}_0\sin \left({\displaystyle \frac{\omega }{c}}(t-t_c)\right)$
  	$𝐀(t)$	$={𝐄}_0{\displaystyle \frac{c}{\omega }}\cos \left({\displaystyle \frac{\omega }{c}}(t-t_c)\right)$

• •

  5: Sin² envelope pulse:

  	$𝐄(t)$	$={𝐄}_0(\sin {\left({\displaystyle \frac{\pi }{c}}(t-t_c)\right)}^2\sin (\omega (t-t_c))$
  		$-\left({\displaystyle \frac{\pi }{c\omega }}\sin \left({\displaystyle \frac{2\pi }{c}}(t-t_c)\right)\cos (\omega (t-t_c))\right))$
  	$𝐀(t)$	$={\displaystyle \frac{{𝐄}_0}{\omega }}\sin {\left({\displaystyle \frac{\pi }{c}}(t-t_c)\right)}^2\cos (\omega (t-t_c))$

• •

  6: Ramp Electric field:

  	$𝐄(t)$	$={𝐄}_0(t-t_0)$
  	$𝐀(t)$	$=-{\displaystyle \frac{{𝐄}_0}{2}}{(t-t_0)}^2$

• •

  7: Delta kick Electric field via Lorentzian:

  	$𝐄(t)$	$={𝐄}_0{\displaystyle \frac{\frac{t_w}{2\pi }}{{(t-t_c)}^2+{\left(\frac{1}{2}{t}_w\right)}^2}}$
  	$𝐀(t)$	$={𝐄}_0{\displaystyle \frac{\arctan \left(2(\frac{(t-t_c)}{t_w})+\frac{\pi }{2}\right)}{\pi }}$

• •

  8: Constant Vector Potenital:

  	$𝐄(t)$	$=\{\begin{array}{cc}{𝐄}_0,t=t_0\hfill & \\ -{𝐄}_0,t=t_{\mathrm{end}}\hfill & \\ 0,\text{otherwise}\hfill & \end{array}$
  	$𝐀(t)$	$={𝐄}_0$

where $t_0=$ t_start, $\omega =2\pi f$ and $f=$ freq, $c=$ cycle, $t_c=$ center, $t_w=$ width, $t_{e}nd=$ t_end, and ${𝐄}_0={(E_x,E_y,E_z)}^T$. Parameters not occurring in a function definition are not referenced and can be set to any value.
Note that user-defined fields can easily be included in the corresponding subroutine src/rt-tddft/src/rt_tddft_ext_field.f90.
Remember using consistent units set by the RT_TDDFT_input_units keyword.

If this setting is specified for a collinear calculation (see keyword spin in the FHI-aims manual), only the spin $=$ spin time-dependent eigenvectors will be excited by the specific external field. This may, e.g., be used to model triplet excitations.

The possible integration schemes are (see, e.g., [111, 53, 17]) :

• •

  exponential_midpoint (EM):

  $$𝐔(t_k+\mathrm{\Delta }t,t_k)=\exp \left(-i\mathrm{\Delta }t{𝐒}^{-1}𝐇(t_k+\mathrm{\Delta }t/2)\right)$$

• •

  crank_nicolson (CN): The order ${\mathrm{N}}_{\mathrm{CN}}$ can be set via the keyword
  RT_TDDFT_crank_nicolson_order order.

  $$𝐔(t_k+\mathrm{\Delta }t,t_k)=\frac{{\sum }_{m=0}^{{\mathrm{N}}_{\mathrm{CN}}}{\left(-i\mathrm{\Delta }t{𝐒}^{-1}𝐇(t_k+\mathrm{\Delta }t/2)\right)}^m}{{\sum }_{m=0}^{{\mathrm{N}}_{\mathrm{CN}}}{\left(i\mathrm{\Delta }t{𝐒}^{-1}𝐇(t_k+\mathrm{\Delta }t/2)\right)}^m}$$

• •

  etrs (Enforced Time-Reversal Symmetry):

  $$𝐔(t_k+\mathrm{\Delta }t,t_k)=\exp \left(-i\frac{\mathrm{\Delta }t}{2}{𝐒}^{-1}𝐇(t_k+\mathrm{\Delta }t)\right)\exp \left(-i\frac{\mathrm{\Delta }t}{2}{𝐒}^{-1}𝐇(t_k)\right)$$

• •

  cfm4 (Commutator-Free Magnus Expansion 4):

  	$𝐔(t_k+\mathrm{\Delta }t,t_k)=\exp \left(-i\mathrm{\Delta }t{𝐒}^{-1}\left({\alpha }_1{𝐇}_1+{\alpha }_2{𝐇}_2\right)\right)\exp \left(-i\mathrm{\Delta }t{𝐒}^{-1}\left({\alpha }_1{𝐇}_2+{\alpha }_1{𝐇}_2\right)\right)$
  	${\alpha }_{1/2}={\displaystyle \frac{1}{4}}\mp {\displaystyle \frac{\sqrt{3}}{6}},{𝐇}_{1/2}=𝐇\left(t_k+\left({\displaystyle \frac{1}{2}}\mp {\displaystyle \frac{\sqrt{3}}{6}}\right)\mathrm{\Delta }t\right)$

• •

  m4 (Magnus Expansion 4):

  	$𝐔(t_k+\mathrm{\Delta }t,t_k)$	$=\left({\displaystyle \frac{\mathrm{\Delta }t}{2}}{𝐒}^{-1}\left({𝐇}_1+{𝐇}_2\right)+\mathrm{\Delta }{t}^2{\displaystyle \frac{\sqrt{3}}{12}}[{𝐒}^{-1}{𝐇}_2,{𝐒}^{-1}{𝐇}_1]\right)$
  	${𝐇}_{1/2}$	$=𝐇\left(t_k+\left({\displaystyle \frac{1}{2}}\mp {\displaystyle \frac{\sqrt{3}}{6}}\right)\mathrm{\Delta }t\right)$

• •

  runge_kutta_4: Standard explicit Runge-Kutta 4 scheme. Works only for very small time steps but is a good choice for doing accurate reference simulations because its properties are well-known.

• •

  parallel_transport_cn: Crank-Nicolson integration scheme within the parallel transport approach [158]. This requires the accelerated Anderson solver, see RT_TDDFT_propagator_solver. [EXPERIMENTAL]

  	$𝐂(t_k+\mathrm{\Delta }t)=$	$𝐂(t_k)-i{\displaystyle \frac{\mathrm{\Delta }t}{2}}({𝐒}^{-1}𝐇(t_k)𝐂(t_k)-𝐂(t_k)\left({𝐂}^{†}(t_k)𝐇(t_k)𝐂(t_k)\right)$
  		$-i{\displaystyle \frac{\mathrm{\Delta }t}{2}}({𝐒}^{-1}𝐇(t_k+\mathrm{\Delta }t)𝐂(t_k+\mathrm{\Delta }t)$
  		$-𝐂(t_k+\mathrm{\Delta }t)\left({𝐂}^{†}(t_k+\mathrm{\Delta }t)𝐇(t_k+\mathrm{\Delta }t)𝐂(t_k+\mathrm{\Delta }t)\right)$

Note that all above mentioned implicit schemes are defined as implicit, as they require evaluation of the Hamiltonian matrix at some future time. Currently, there are two ways implemented on how to solve this type of equation which can be set by the RT_TDDFT_propagator_solver keyword.
It should also be noted here that, when using a predictor-corrector scheme (which is the default solver for all implicit propagators), accuracy can be enhanced by using extrapolation, see the keyword RT_TDDFT_extrapolate_predictor.
 
Some additional schemes are implemented for testing purposes, namely the M6 (m6), CFM4:3 (cfm4_3) and CFM6:5 (cfm6_5) schemes [17].
 
For general purpose, we recommend to use the crank_nicolson propagation scheme.

This keyword must usually not be set explicitly since every RT_TDDFT_propagator has its own default setting for the solver. Anyway, one can set the solver manually, e.g., to use one of the implicit exponential propagators with extrapolation, which is much cheaper since the corrector integration of the Hamiltonian matrix is not applied (however, this may be much less stable, not justifying the setting). Some combinations are not possible and a warning will be given that default settings will be used. Please use this option with caution.

Performing additional corrector steps until a convergence criterion is reached will formally increase accuracy. It should nevertheless be noted that the biggest error is probably introduced already by the predictor step and that significant improvements by doing additional corrector steps are not to be expected. Every corrector update involves another integration of the Hamiltonian matrix which is the most expensive part in RT-TDDFT usually. Before doing more corrector steps, it probably makes more sense to choose a smaller time step – if accuracy is a problem. If one wants to use this functionality anyway, values of ${10}^{-5}$-${10}^{-9}$ would make sense usually. Note that a maximum number of 10 corrector steps is hardcoded to avoid unnecessarily many corrector updates (one should check settings when this happens).

Please see RT_TDDFT_use_precor_tol for further remarks. The default value is usually a good choice.

The Hamiltonian (matrix) used in the predictor step will be extrapolated by a given method which can be set via the RT_TDDFT_ham_extrapolation keyword. This also means that additional storage is needed. The extrapolated Hamiltonian (matrix) is applied in the predictor step. Tests indicate that extrapolation can noticeably improve accuracy. However, one should check this setting when unexpected instabilities occur.

The default setting indicating linear extrapolation works well in tests, whereas, in contrast, the Lagrange extrapolation approach sometimes produces odd results and should be used only with caution. Some of our tests indicate that higher-order extrapolation ($\ge 4$) can yield instabilities, thus one should check this.

This only makes sense if an explicit scheme like runge_kutta_4 is used as a predictor, or, and this is the usual setting here, if any applicable implicit scheme is used in combination with extrapolation. A well-functioning scheme consists of an extrapolated exponential_midpoint predictor in combination with the cfm4 propagator for the corrector. This keyword is ignored if no predictor-corrector scheme is in use.

Currently, three different approaches to compute matrix exponentials needed for the exponential-type propagators can be used. The diagonalization-based approach is usually working perfectly but in the current implementation, it cannot be used for non-Hermitian matrices – this can be the case for effective Hamiltonian matrices occurring in Ehrenfest dynamics. In this case, the scaling and squaring - method [142] will be applied. Anyway, this setting will then be adjusted automatically. One can modify this setting mainly for benchmarking or testing. Note that using the taylor method requires to choose an order of the expansion, see
RT_TDDFT_exponential_taylor_order.

The Taylor method applies to general matrices and can be an economic choice. However, the order must be set with caution, as we found that the here employed all-electron approach likely requires higher orders. It may not even work when heavy elements are involved. This setting also highly depends on the time step. We found values of 10-20 to work well in certain cases.

Usually, 1st order seems to suffice for most applications. However, it can be beneficial to use higher order, e.g., 2/3, when instabilities are observed, especially possibly for heavier elements.

This only makes sense when the software bug explained above occurs. Setting the flag to .true. should resolve the issue. It is possible that the performance is impacted by this setting.

This can be used to prevent a simulation of ’blowing up’, i.e., becoming unstable, as measured by the degree of energy conservation. The start time should be chosen such that the reference value makes sense: for example, when applying a finite resonant excitation, the total energy will increase, but should stay constant thereafter; thus, the reference value should be set a while after the excitation ends.

Imaginary time propagation can be used as an alternative approach to the standard SCF procedure to converge a ground-state solution [93]. This does only require a time step as control parameter and may be used to converge situations where SCF fails. Note that we recommend the Crank-Nicolson propagator for this functionality (higher orders may also be beneficial in order to enable higher time steps).
After convergence is reached, the eigenvectors are written back into the native data structures so that postprocessing methods should be usable.

These parameters are explained in detail by Walker et al. [311]. We refer to algorithm ’AA’ in this reference. The parameter $m$ = order, $k_{\max }=$ iter as in the reference. The tolerance parameter depends on the setting RT_TDDFT_anderson_cond_type and should generally be chosen as a very small number. The most important parameter here is iter, which can be crucial to obtain a converged fixed-point iteration, but may also increase the possible numerical demand.

This defines which quantity is used as convergence measure in the fixed-point iteration. ev_maximum defines that $\max \{|{\psi }_n^{(i)}-{\psi }_n^{(i+1)}|\forall n\in {\mathrm{N}}_{\mathrm{occ}}\}$ must be below tol as defined via RT_TDDFT_anderson_order_iter. ev_individual means that all solutions must be below the threshold tol. delta_rho only requires that the average change in charge density $\mathrm{\Delta }{\rho }^i$ must be below the tolerance.

The time step for Ehrenfest dynamics must obviously be at least the same value as for electron propagation. Usually, it can be chosen around 5- to 20-fold while maintaining stable propagation. Since updating the geometry and calculating forces is expensive, the Ehrenfest time step should be as large as possible. Additionally, the specific interplay of time propagation, force update and geometry update works very good when dt_geo $=N\times$ dt with $N=2^n$. A reasonable choice for doing Ehrenfest dynamics should be dt $=0.1$ a.u. in combination with dt_geo $=0.4-0.64$ a.u. Monitoring energy conservation and norm conservation of the eigenvectors is advised when trying out settings.
Remember using consistent units set by the RT_TDDFT_input_units keyword.

Usually, Ehrenfest dynamics start with the begin of a RT-TDDFT simulation when requested, but employing this key can modify the start time in case specific tasks require this.
Remember using consistent units set by the RT_TDDFT_input_units keyword.

Usually, it seems that using the ’normal’ non-adiabatic Ehrenfest forces conserving the energy (but not necessarily the momentum) is sufficient, but by using this keyword, the complete expression is evaluated – this might be important in specific situations and is of course more general. Anyway it is more expensive.

All output files have the format PREFIX.rt-tddft.OBSERVABLE.SUFFIX or
PREFIX.rt-tddft-ehrenfest.OBSERVABLE.SUFFIX where PREFIX can be set by the user to specify simulation settings. For example, energies are by default written to a file named output.rt-tddft.energies.dat. The corresponding quantities are always specified (format, units, etc.) in comment headers inside the files. Note that for any file output, existing files will not be overwritten. Instead, new files named PREFIX.N.rt-tddft.OBSERVABLE.SUFFIX with N $=1,2,..$ will be created.

The external field and all corresponding parameters are written to a file before a real-time simulation in this case. This information can, e.g., be used for postprocessing or visualization. Note that our postprocessing tool eval_tddft.py requires this file.

Default cube output is deactivated until this flag is given and the output frequency is specified. All options that can be used for cube output via the output cube keyword can be applied here and must be set in addition to this keyword.
Remember using consistent units set by the RT_TDDFT_input_units keyword.

Since RT-TDDFT simulations are usually performed for a significant number of time-integration/density update/force calculation/etc. operations, the associated output will result in very large amounts of mostly unnecessary data which can be avoided by setting this keyword. The default should yield sufficient information on what is happening, but setting it to 2 will give some more information about numerical properties which can be important to follow.
Please note that important observables, e.g., energies, dipole moment, geometry, etc. are written to separate files which can be controlled further and which eases post-processing (see RT_TDDFT_write_file_prefix).

In the RT-TDDFT module, several routines from the parent code are called which may produce their own output. By default, the output level for these calls is set to minimal, i.e., not output will be sent to the console. However, this can be overridden by this keyword. We refer to the keywork output_level in the FHI-aims manual for specific choices.

Energy should always be observed in order to validate a properly running simulation, i.e., total energy conservation. Note that when values are written to a file too, the corresponding file is named PREFIX.rt-tddft.energies.dat where PREFIX can be changed via the RT_TDDFT_write_file_prefix keyword.

The electronic dipole is usually an observable of major interest in many RT-TDDFT applications. Note that when values are written to a file too, the corresponding file is named PREFIX.rt-tddft.dipole.dat where PREFIX can be changed via the
RT_TDDFT_write_file_prefix keyword.

This keyword can be used to save unnecessary computation time for components that are not of interest. However, the computational effort for observables is generally low, as compared to the remaining functionality.

The computation of transition dipole moments can be requested here. The lower and upper index define the set of occupied/unoccupied Kohn-Sham states for which time-dependent transition dipole moments will be computed. It is necessary that the native keyword empty_states is set accordingly high (e.g., to the number of basis functions or, e.g., 10000), else it is possible that not enough unoccupied states are available in the corresponding arrays. This is one of the rather computationally expensive observables.

This keyword can be used to save unnecessary computation time for components that are not of interest.

If given, the dipole expectation values for all occupied (and thus, time-propagated) states are computed and sent to output. Note that this computation requires the explicit evaluation of the dipole matrix $⟨{\varphi }_i|𝐫|{\varphi }_j⟩$ and could cause additional computational cost.

The magnetic moment is usually an observable of major interest in RT-TDDFT applications for circular dichroism. Note that when values are written to a file too, the corresponding file is named
PREFIX.rt-tddft.magmom.dat where PREFIX can be changed via the
RT_TDDFT_write_file_prefix keyword.

This keyword can be used to save unnecessary computation time for components that are not of interest. However, the computational effort for observables is generally low, as compared to the remaining functionality.

The electronic current is often the main observable of interest in periodic RT-TDDFT simulations (analog to the dipole in finite systems). Note that when values are written to a file too, the corresponding file is named PREFIX.rt-tddft.current.dat where PREFIX can be changed via the RT_TDDFT_write_file_prefix keyword.

Note that when values written are to a file too, the corresponding file is named PREFIX.rt-tddft-ehrenfest.trajectory.xyz where PREFIX can be changed via the
RT_TDDFT_write_file_prefix keyword.

This quantity is computed as

$$N_e^{\mathrm{ex}}(t)=\sum _{n=1}^{{\mathrm{N}}_{\mathrm{occ}}}{f}_n(N_{electrons}-{|⟨{\psi }_n(0)|{\psi }_n(t)⟩|}^2),$$

(3.113)

i.e., via projection of time-dependent orbitals on ground-state orbitals [238]. This refers to the number of electrons that are possibly excited within the introduced dynamics.

All necessary information to perform a RT-TDDFT restart is written to file/s in this case. In detail, all the eigencoefficients, coordinates, settings, etc. Multiple entries of this keyword for different times can be given. Note that this can potentially produce huge files, especially when periodic boundary conditions are in use (as for every $𝐤$-point, eigenvectors are written).

See RT_TDDFT_restart_write for details. When doing very demanding simulations and/or very long simulations, one should set a sensible value here to avoid needing to restart from the beginning after a job was canceled. Writing a restart file automatically when a kill signal is received is on the to-do list.

Besides the physical quantities, the restart file contains all relevant real-time simulation parameters which are compared to what is given in control.in. The settings must be the same or the simulation will abort. This also works for Ehrenfest dynamics. Note that when performing a restart, FHI-aims will run a SCF cycle, but the resulting eigenvectors, etc. will be overwritten from what is read in the restart file - based on this, the real-time simulation will be re-initialized.