# Real-Time Propagation and Ehrenfest MD

In real time time-dependent DFT, instead of solving the static Schroedinger equation (SE), the aim is to solve the time dependent SE (TDSE)

Contrary to the time independent case, there is no variational principle for the total energy, and the total energy has to be replaced by the quantum mechanical action,

for which it is valid that the function \(\Phi(t)\) that makes the action stationary will be the solution. The resulting time-dependent Kohn-Sham (TD-KS) equations read as

with \(\phi_i\) the KS orbitals and the time dependent density

\(\textbf{F}(t)\) is a time dependent external field, if it has been specified. By applying an external electric field, we intend to simulate the interaction between the electron and an external radiation. The light is then treated classically using, either the length gauge for non periodic systems, or the velocity gauge for periodic systems.

A rather general derivation for Ehrenfest (non-adiabatic quantum-classical molecular) dynamics can be obtained starting from the action of a system. For a situation in which the the electrons are treated quantum mechanically while the nuclei are treated classically the total action can be written as the sum of the two environments

and the quantum mechanical action as above defined. The equations of motion are then derived by making the action stationary

Evaluating these expressions in the framework of TDDFT, the equations of motion become

for the nuclear motion, wile for the electrons the time dependent SE as given above is used. These equations are valid for the exact wavefunctions. In the Kohn-Sham approach, the wavefunctions are replaced by a linear combination of basis functions. For plane waves the equations remain the same, since they do not depend on the nuclear coordinates and therefore are independent of the change of the ionic positions during MD. The case of a Gaussian basis set requires a more detailed analysis. Using Gaussian basis sets for the expansion of the wavefunctions, an implicit dependence of the wavefunctions on the nuclear positions is introduced. Since in an MD approach the nuclear positions are function of time, the time derivative in the quantum mechanical action has to be replaced by the total time derivative

Due to the introduction of the basis, the independent variables for making the action constant become the expansion coefficients \(C_{j\alpha}(t)\) of the molecular orbital \(j\) in the contracted basis functions \(\phi_\alpha({\bf r},t)\)

where

Hence, in Ehrenfest dynamics an additional contribution due to the derivative of the basis functions becomes part of the Hamiltonian used in the exponential operator. Instead of being purely imaginary, the matrix in the exponential of the propagator becomes fully complex.

## Running real time time dependend DFT in CP2K

To run RTP or Ehrenfest MD with CP2K, the corresponding RUN_TYPE in
the [GLOBAL] (#CP2K_INPUT.GLOBAL) section has to be selected, and the MD
section has to be present to specify the time step length `TIMESTEP`

and the desired number of steps
(`STEPS`

). It is crucial to set an appropriate time step to propagate the electronic degrees of
freedom, i.e., in the order of atto-seconds. All other input parameters related to the RT-TDDFT run
are specified from the section
REAL_TIME_PROPAGATION. The simulation needs to
start from the electronic density at \(t_0\) (`INITIAL_WFN`

). This can be the ground state obtained by
means of an initial SCF optimization (`SCF_WFN`

) or providing the restart file of a previously
computed state (`RESTART_WFN`

). The propagation can also be restarted setting the initial
wavefunction as `RT_RESTART`

and providing the correct restart file, which has a different format
with respect to the SCF restart file. The path to the restart file is in both cases provided by the
usual WFN_RESTART_FILE_NAME key.

Three different propagators are available in CP2K, the enforced time-reversible symmetry propagator (ETRS), the exponential midpoint (EM) propagator, and the Crank-Nicholson propagator which can be seen as a first order Pad’e approximation of the EM propagator. The ETRS approach starts with an exponential approximation to the evolution operator \(\hat{U}(t,0)=\exp{-it\hat{H}(0)}\), to then compute the final time-reversible and unitary propagator self-consistently. In the real time propagation scheme (fixed ionic positions) the self consistent solution only involves the calculation of the new Kohn-Sham matrix for the propagated coefficients. For Ehrenfest MD, the iterative procedure is embedded into the integrator of the nuclear EOM. Hence, each iteration step for Ehrenfest dynamics involves a real time propagation step and a complete evaluation of the forces. Due to this, more iterations are needed to reach self consistency for the propagator. The convergence criterion implemented in CP2K is defined as

with \(\Delta C\) as the difference of coefficient matrices in two successive steps, and \(\epsilon\)
given by `EPS_ITER`

.

In terms of computational cost, the most expensive part is the evaluation of the matrix exponential
in the propagator. Four different methods among those listed in have been implemented in CP2K, i.e.,
the Taylor expansion, he diagonalization, the Pad’e approximation, and the Arnoldi subspace
iteration, to be selected via the `MAT_ESP`

input key. The Arnoldi method often provides a superior
performance. Comparing the theoretical scaling, the Arnoldi method is expected to be about 5 times
as fast as Pad’e or Taylor. However, the Pade approximation can be sometime the faster and more
stable choice than the Arnoldi method (e.g., large time step). Since propagation schemes require an
iterative procedure, it is convenient to apply an extrapolation scheme in order to speed up
convergence. For Ehrenfest dynamics, an extrapolation based on the product of the density and the
overlap matrix turns out to be the a good choice. Nevertheless, the quality of the different
extrapolation strongly depends on the system applied to and the settings of the simulation.
Therefore, which method to use and the extrapolation order needs to be tested on the system of
interest. For very large systems, the density-based method can be used to achieve linear scaling
performance (see the
DENSITY_PROPAGATION
keyword).

## The external electric field

One of the most relevant application domains for RT-TDDFT is the study of light?matter interactions, e.g., in the field of spectroscopy, excited state dynamics and radiation damage. To mimic these phenomena, at any time during the propagation, it is possible to apply a time dependent electric field \({\bf E}(t)\). The applied field is in general modulated by an envelope function, \(E_{\text{env}}(t)\) and is defined as:

Where \(\textbf{P}\) is the field polarization, \(\omega_{0}\) is the carrying frequency at which the field oscillates, and \(\phi\) its initial phase. The characteristic of the applied field as well as its time extension are provided through the section EFIELD. The time dependent electric field defined within this section can only be used in combination with RT-TDDFT. By default the coupling between the electric field and the electronic degrees of freedom is described within the length gauge, by adding to the Hamiltonian the dipole coupling term \( e \textbf{E}(t) \cdot {\bf r} \). This approach is only valid for isolated molecular systems, and not when periodic boundary conditions are applied. The velocity-gauge form of the equations suitable for periodic systems is obtained through a gauge transformation involving the vector potential

In the time dependent KS equations, the vector potential \({\bf A}(t)\) appears in the kinetic energy
term and, in the case where non-local pseudopotentials are used, the gauge field also transforms the
electron?ion interaction. To use this representation the `VELOCITY_GAUGE`

input keyword has to be
activated. The total energy varies with time as the applied external field interacts with the
system. To monitor the time evolution of the field as well as of the various terms contributing to
the total electronic energy, the corresponding print key can be activated from the
REAL_TIME_PROPAGATION section.

with \(\Delta C\) as the difference of coefficient matrices in two successive steps, and \(\epsilon\)
given by `EPS_ITER`

.

## Resonant Excitation

The input file provided below starts a RT–TDDFT simulation for an isolated carbon monoxide molecule perturbed by a resonant X-ray field. Such real-time calculation aims to promote electrons from some core occupied state to empty. What states are going to be addressed depends on the frequency of the field, which has to be resonant to some core state transition. The electronic dynamics are propagated as the field is applied, leading to an electronic excitation increasing with time, if the field frequency matches the aimed electronic transition energy.

This equation is propagated for each \(i\) electron with a time step \(\delta t\): the MOs evolved
collectively over time in a **continuous manner**. Especially, the energy evolution is continuous:
the electrons cannot absorb exactly one photon at a given frequency to instantaneously go from one
electronic state to another.

If one applies a field resonant with a given electronic transition,
the MOs involved in this transition will gradually be transformed from their initial state to the
corresponding excited state. This will happen in real-time without defining any specific electronic
transition before starting the simulation. Some preliminary calculations to determine the spectral
features of the system under study are in general useful to better define the desired properties of
the perturbing field. For core state excitations these information can be obtained by XAS
simulations.

Once the RTP has started, the evolution of the electronic structure can be monitored by means of several descriptors, like the time dependent dipole moment, current density, total electronic density, as well as the projection of the time dependent orbitals onto some preselected reference states. This latter is particularly useful to verify the variation in population of specific states, while monitoring density differences is usually helpful to detect charge transfer processes.

By projecting the time-dependent molecular orbitals onto some reference states (e.g., the initial MOs) means computing the overlap between the propagated orbital \(\psi_i({\bf r}, t) = \sum_\alpha C_{i\alpha}(t)\phi({\bf r})\) and any reference orbital \(\psi^{\text{ref}}_m\). For instance, considering as reference orbitals the static unoccupied ones, the quantity

is an estimate of the number of excited electrons.

### CP2K input

RTP.inp is the input file used for such simulation:

```
@set EXC_STATE_1 LR-xasat2_1s_singlet_idx1-1.wfn
@set EXC_STATE_2 LR-xasat2_1s_singlet_idx2-1.wfn
&GLOBAL
PROJECT RTP
RUN_TYPE RT_PROPAGATION
PRINT_LEVEL MEDIUM
&END GLOBAL
&MOTION
&MD
ENSEMBLE NVE
STEPS 5000
TIMESTEP [fs] 0.00078
TEMPERATURE [K] 0.0
&END MD
&END MOTION
&FORCE_EVAL
METHOD QS
&DFT
&REAL_TIME_PROPAGATION
MAX_ITER 100
MAT_EXP ARNOLDI
EPS_ITER 1.0E-11
INITIAL_WFN SCF_WFN
&PRINT
&FIELD
FILENAME =applied_field
&END FIELD
&PROJECTION_MO
REFERENCE_TYPE SCF
REF_MO_FILE_NAME RTP-RESTART.wfn
REF_MO_INDEX -1
SUM_ON_ALL_REF .FALSE.
TD_MO_INDEX -1
SUM_ON_ALL_TD .FALSE.
&PRINT
&EACH
MD 1
&END EACH
&END PRINT
&END PROJECTION_MO
&PROJECTION_MO
REFERENCE_TYPE XAS_TDP
REF_MO_FILE_NAME ${EXC_STATE_1}
TD_MO_INDEX -1
SUM_ON_ALL_TD .FALSE.
&PRINT
&EACH
MD 1
&END EACH
&END PRINT
&END PROJECTION_MO
&PROJECTION_MO
REFERENCE_TYPE XAS_TDP
REF_MO_FILE_NAME ${EXC_STATE_2}
TD_MO_INDEX -1
SUM_ON_ALL_TD .FALSE.
&PRINT
&EACH
MD 1
&END EACH
&END PRINT
&END PROJECTION_MO
&END PRINT
&END REAL_TIME_PROPAGATION
&EFIELD
ENVELOP GAUSSIAN
! gaussian env in fs units
&GAUSSIAN_ENV
SIGMA [fs] 0.3073
T0 [fs] 1.3190
&END GAUSSIAN_ENV
! in W.cm-2
INTENSITY 4.08E+13
PHASE 0.0
POLARISATION 1 0 0
! this is 529 eV:
WAVELENGTH 2.34374655955
&END EFIELD
BASIS_SET_FILE_NAME BASIS_PCSEG2
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 1000
NGRIDS 5
REL_CUTOFF 60
&END MGRID
&QS
METHOD GAPW
EPS_FIT 1.0E-6
&END QS
&SCF
MAX_SCF 500
EPS_SCF 1.0E-8
&END SCF
&POISSON
POISSON_SOLVER WAVELET
PERIODIC NONE
&END POISSON
&XC
&XC_FUNCTIONAL PBE
&PBE
SCALE_X 0.55
&END
&END XC_FUNCTIONAL
&HF
FRACTION 0.45
&INTERACTION_POTENTIAL
POTENTIAL_TYPE TRUNCATED
CUTOFF_RADIUS 7.0
&END INTERACTION_POTENTIAL
&END HF
&END XC
&PRINT
&MULLIKEN OFF
&END MULLIKEN
&HIRSHFELD OFF
&END HIRSHFELD
&MOMENTS
PERIODIC .FALSE.
FILENAME =dipole
COMMON_ITERATION_LEVELS 100000
&EACH
MD 1
&END EACH
&END MOMENTS
&END PRINT
&END DFT
&SUBSYS
&CELL
ABC 10 10 10
ALPHA_BETA_GAMMA 90 90 90
PERIODIC NONE
&END CELL
&TOPOLOGY
COORD_FILE_NAME carbon-monoxide_opt.xyz
COORD_FILE_FORMAT XYZ
&END TOPOLOGY
&KIND C
BASIS_SET pcseg-2
POTENTIAL ALL
&END KIND
&KIND O
BASIS_SET pcseg-2
POTENTIAL ALL
&END KIND
&END SUBSYS
&END FORCE_EVAL
```

In the following the chosen parameters for real-time propagation, electric field, and the time-dependent projection are discussed.

### Real-Time Propagation parameters

In this simulation, we start from an SCF-optimized state which corresponds to the ground state by
setting INITIAL_WFN to `SCF_WFN`

.
Before the Real-Time Propagation starts, a ground state calculation will thus be performed. We use a
Molecular Orbital-based description of the wave function for the propagation along with the Arnoldi
approach to compute the exponential by setting
MAT_EXP to `ARNOLDI`

. For each time
step, the ERTS algorithm is used with a convergence threshold defined by
EPS_ITER, which we set to `1.0E-11`

.
Note that the smaller this threshold, the more iterations per time step will be needed to converge.
Hence, we have set the maximal iteration number
MAX_ITER at quite high value of 100.
The time step used is rather small since we have to describe a core-hole excitation process that
takes place with a typical frequency of 529 eV. Following the rule of thumb to set the time step to
about 10 times the field frequency, we set TIMESTEP to
`[fs] 0.00078`

### Field parameters

The field is defined by its envelope, its intensity, its polarization along the laboratory x, y, and z-axis, its wavelength, and the original phase. Several types of field envelopes can be used. Here we use a Gaussian one with a width of \(\sigma=0.3073\) fs and centered at \(T0=1.3190\) fs. The intensity used is 4.08E+13 W.cm\(^{-2}\). Along with a carrying frequency of 529 eV (approximately 2.34374655955 nm), it should promote about \(10^{-3}\) from the Oxygen 1s to the first available excited state.

```
&EFIELD
ENVELOP GAUSSIAN
&GAUSSIAN_ENV
SIGMA [fs] 0.3073
T0 [fs] 1.3190
&END GAUSSIAN_ENV
INTENSITY 4.08E+13
PHASE 0.0
POLARISATION 1 0 0
WAVELENGTH 2.34374655955
&END EFIELD
```

The total number of STEPS is set to `5000`

so that the field envelope
fits in the period of time simulated.

### Time-Dependent Projection parameters

Time-dependent projections can be defined in the REAL_TIME_PROPAGATION/PRINT section. One can define several PROJECTION_MO sections which will be done consequitively and the results stored in different files.

First, let us have a look at the projection toward the ground state:

```
&PROJECTION_MO
REFERENCE_TYPE SCF
REF_MO_FILE_NAME RTP-RESTART.wfn
REF_MO_INDEX -1
SUM_ON_ALL_REF .FALSE.
TD_MO_INDEX -1
SUM_ON_ALL_TD .FALSE.
&PRINT
&EACH
MD 1
&END EACH
&END PRINT
&END PROJECTION_MO
```

All the time-dependent Molecular Orbitals are projected by setting
TD_MO_INDEX to
`-1`

and these projections are stored separately by setting
SUM_ON_ALL_TD
to `.FALSE.`

.

This calculation is spin-independent so one does not have to define the spin of the MO to project.
The reference to projected to is loaded from the file `RTP-RESTART.wfn`

, which is the ground state
obtained after the SCF cycles. Note that you can define a wave-function that is not the ground state
as long as the wave-function description (basis set, number of electrons…) is the same as the one
used for the real-time propagation. The
REFERENCE_TYPE
defaults to `SCF`

.

All the molecular orbitals available in this reference wave-function will be used as a reference for
the projection by setting
REF_MO_INDEX to
`-1`

, and stored in separate files by setting
SUM_ON_ALL_REF
to `.FALSE.`

. The projection will be computed at each time step.

There are \(N_e/2 = 7\) molecular orbitals for carbon monoxide. There are thus 7 time-dependent MOs (the \(i\)s) and 7 reference MO (the \(m\)s), so that there will be \(7x7=49\) projection per time step:

where \(S_{\alpha \beta}\) is the overlap matrix between basis set orbitals.

Now, let us focus on the second and third projections which can be viewed as projections toward excited-states:

```
&PROJECTION_MO
REFERENCE_TYPE XAS_TDP
REF_MO_FILE_NAME ${EXC_STATE_1}
TD_MO_INDEX -1
SUM_ON_ALL_TD .FALSE.
&PRINT
&EACH
MD 1
&END EACH
&END PRINT
&END PROJECTION_MO
```

In this case, all the time-dependent MO are involved in the projection and stored separately. This
time, the reference wave function is supposed to be from an XAS_TDP calculation, see Linear Response
input file in the
tutorial archive. Running with
the proper parameters, this XAS_TDP run saves one .wfn file per excited state. Using
`REFERENCE_TYPE XAS_TDP`

, the projection uses automatically the state saved in the produced wave
function file as the reference:

where \(C_{\omega\alpha}\) is the \(a^{\text{th}}\) atomic coefficient of the excited state found in the XAS_TDP module. There will be 7 projections per time step in this case: one for each time-dependent MO. The excited state population associated with this excited state is:

The carbon monoxide molecule has a rotational symmetry along its CO bond. If one notes this axis \(z\), then the \(x\) and \(y\)-axis are equivalent by symmetry. It happens that the first available excited state for the Oxygen 1s is degenerate: there are two available excited states orthogonal in the \(xy\)-plane. That is why a third projection is requested which uses the other excited state proposed by the XAD_TDP module. Therefore, the excited state should be understood as the sum over the two equivalent excited states:

where \(\omega'\) stands for the other equivalent excited state, with \(\omega=\omega'=529\) eV.