# Ehrenfest dynamics implemented in the all-electron package exciting

Ronaldo Rodrigues Pela1,2,3, Claudia Draxl1,2 1 Physics Department and IRIS
Adlershof, Humboldt-Universität zu Berlin, Zum Großen Windkanal 2, 12489
Berlin  2 European Theoretical Spectroscopy Facility (ETSF)  3 Supercomputing
Department, Zuse Institute Berlin, Berlin, Germany<EMAIL_ADDRESS>
###### Abstract

Ehrenfest Dynamics combined with real-time time-dependent density functional
theory has proven to be a reliable tool to study non-adiabatic molecular
dynamics with a reasonable computational cost. Among other possibilities, it
allows for assessing in real time electronic excitations generated by ultra-
fast laser pulses, as e.g., in pump-probe spectroscopy, and their coupling to
the nuclear vibrations even beyond the linear regime. In this work, we present
its implementation in the all-electron full-potential package exciting. Three
cases are presented as examples: diamond and cubic BN relaxed after an initial
lattice distortion, and cubic BN exposed to a laser pulse. Comparison with the
Octopus code exhibits good agreement.

††: Electronic Structure

##  1 Introduction

Molecular dynamics (MD) is one of the most important techniques of
computational materials science. In one of the first reports on a computer MD,
Alder and Wainwright described it as tracking the movement of atoms which
interact by means of a given potential [1, 2]. In non-empirical MD, the
potential is not modeled _a priori_ via some parametrized expression, but
rather obtained “on the fly” through _ab initio_ calculations [2]. Quite often
then, the Born-Oppenheimer approximation (BOA) is employed. Thereby, each
nucleus is treated as a point charge moving under the Coulomb potential due to
the electrons and other nuclei, whereas the electrons are treated quantum
mechanically with a wavefunction that responds adiabatically to the nuclei
movement, remaining in an eigenstate in each time [3].

Despite many successful theoretical predictions, BOA is not appropriate for
several applications, e.g., those involving the interaction between
electromagnetic waves and matter. Processes in which the excited electrons
respond to a rapidly varying time-dependent field and then lose their energy
by coupling to lattice vibrations, are non-adiabatic and thus have to be
treated beyond the BOA [4, 5, 6, 7, 8]. Examples of such phenomena of growing
interest are pump-probe experiments, in which an initial laser pulse excites
the system, while a second one monitors its time-dependent relaxation towards
the ground state (or a stable excited-state configuration).

An efficient method to describe such non-adiabatic processes is the Ehrenfest
Dynamics [9, 10] that combines a classical treatment of the nuclei with a non-
adiabatic evolution of electronic wavefunctions – usually carried out in the
framework of real-time time-dependent density-functional theory (RT-TDDFT).
This approach has been successfully applied in many recent studies, as diverse
as electron injection into BN nanotubes with thousands of atoms [11],
collisions between oxygen and graphite clusters [12], excited carrier dynamics
in carbon nanotubes [13], ion collisions with carbon nanostructures [14],
photoinduced dynamics in semiconductor quantum dots [15], ultrafast electron
injection from a PbSe quantum dot into a TiO2 surface [16], electron and hole
transfer at a polymer-carbon nanotube heterojunction [17], or chlorine ion
colliding to a graphene nanoflake consisting of hundreds of atoms [18].

A non-negligible aspect in computational materials science, in particular
considering that we have entered the big-data era, is the validation of data.
In electronic- structure theory, there are different approaches to solve the
Kohn-Sham (KS) equations of DFT. Linearized augmented planewaves (LAPW)
combined with local-orbitals (lo) can be regarded as the basis to describe KS
wavefunctions with ultimate precision, thus serving as a reference for other
implementations or approximations [19, 20]. In this sense, the full-potential
LAPW+lo exciting code [21] can be considered as reference for high-precision
calculations. Moreover, it is open source and user-friendly, with more than 50
tutorials explaining its features and how to explore them [22]. In this
manuscript, we report how we combine Ehrenfest Dynamics with RT-TDDFT, another
most recent implementation [23]. To this extent, we briefly summarize the
theory and provide the formalism behind our implementation together with LAPW-
specific derivations. Finally, we show three different examples of
applications, before drawing our conclusions.

##  2 Theoretical framework

Within Ehrenfest dynamics, the nuclei behave as classical repelling point
charges that move being submitted to an external electric field and the one
imposed by themselves. The electrons dynamics, in turn, is typically modeled
by RT-TDDFT. Considering the external electric field given in the gauge
$\mathbf{E}(t)=-\frac{1}{c}\frac{d\mathbf{A}(t)}{dt}$, $\mathbf{A}(t)$ being
the vector potential and $c$ the speed of light, then the action integral for
the total system of nuclei and electrons can be written as [24, 25]

$\displaystyle\mathcal{A}$ $\displaystyle=$
$\displaystyle\int_{t_{1}}^{t_{2}}\mathrm{d}t\sum_{J}\left[\frac{M_{J}\dot{\mathbf{R}}_{J}^{2}}{2}+\frac{Z_{J}}{c}\dot{\mathbf{R}}_{J}\cdot\mathbf{A}(t)\right]+$
(1) $\displaystyle+$
$\displaystyle\int_{t_{1}}^{t_{2}}\mathrm{d}t\sum_{j}f_{j\mathbf{k}}w_{\mathbf{k}}\left\langle\psi_{j\mathbf{k}}(t)\left|\mathrm{i}\frac{\partial}{\partial
t}-\hat{H}(\mathbf{r},\mathbf{R},t)\right|\psi_{j\mathbf{k}}(t)\right\rangle.$

$M_{J}$ is the mass of a nucleus with index $J$, $Z_{J}$ its atomic number,
and $\mathbf{R}_{J}(t)$ its position. $|\psi_{j\mathbf{k}}(t)\rangle$ is the
single-particle KS wavefunction with wavevector $\mathbf{k}$, state index $j$,
and occupation $f_{j\mathbf{k}}$. $w_{\mathbf{k}}$ is the weight of the
k-point. The time-dependent KS hamiltonian is

$\hat{H}(\mathbf{r},\mathbf{R},t)=\frac{1}{2}\left(-\mathrm{i}\nabla+\frac{1}{c}\mathbf{A}(t)\right)^{2}+v_{H}(\mathbf{r},t)+v_{XC}(\mathbf{r},t)+v_{nucl}(\mathbf{r},\mathbf{R}),$
(2)

where $v_{H}$ and $v_{XC}$ are the time-dependent Hartree and exchange-
correlation potentials. $v_{nucl}(\mathbf{r},\mathbf{R})$ is the electrostatic
potential that comprises all the nuclei-nuclei and nuclei-electron
interactions

$v_{nucl}(\mathbf{r},\mathbf{R})=\displaystyle{\frac{1}{2}\sum_{I,J}^{I\not=J}\frac{Z_{I}Z_{J}}{|\mathbf{R}_{I}-\mathbf{R}_{J}|}-\sum_{J}\frac{Z_{J}}{|\mathbf{r}-\mathbf{R}_{J}|}}.$
(3)

###  2.1 Electronic equation of motion

exciting employs the LAPW+lo basis [21, 26]. Around the nuclei, non-
overlapping spheres (termed muffin-tin spheres or atomic spheres) are defined.
The space outside the muffin-tin (MT) spheres is called the interstitial
region. An LAPW function $|\phi_{\mathbf{G}+\mathbf{k}}\rangle$ with
wavevector $\mathbf{G}+\mathbf{k}$, where $\mathbf{G}$ is a reciprocal lattice
vector and $\mathbf{k}$ a crystal momentum from the first Brillouin zone (BZ),
are described by planewaves in the interstitial region, each of them being
smoothly augmented into atomic-like wavefunctions inside the atomic spheres.
lo’s, $|\phi_{\nu}\rangle$, here indexed by $\nu$, are functions defined only
inside one muffin-tin sphere and zero outside. In Ehrenfest dynamics, the
basis is implicitly time-dependent, since the muffin-tin spheres are centered
around the nuclei whose positions evolve in time. The KS wavefunctions
$|\psi_{j\mathbf{k}}(t)\rangle$ are expanded in terms of the basis as

$|\psi_{j\mathbf{k}}(t)\rangle=\sum_{\mathbf{G}}C_{j\mathbf{k}\mathbf{G}}(t)|\phi_{\mathbf{G}+\mathbf{k}}\rangle+\sum_{\nu}C_{j\mathbf{k}\nu}(t)|\phi_{\nu}\rangle=\sum_{\mu}C_{j\mathbf{k}\mu}(t)|\phi_{\mu\mathbf{k}}\rangle.$
(4)

Using the more compact notation indicated by the last summation in Eq. (4),
the time-evolution of the array $C_{j\mathbf{k}}(t)$, obtained by minimizing
the integral of action, follows

$S_{\mathbf{k}}(t)\dot{C}_{j\mathbf{k}}(t)=-\mathrm{i}H_{\mathbf{k}}(t)C_{j\mathbf{k}}(t)-B_{\mathbf{k}}(t)C_{j\mathbf{k}}(t),$
(5)

where
$[S_{\mathbf{k}}(t)]_{\mu\mu^{\prime}}=\langle\phi_{\mu\mathbf{k}}|\phi_{\mu^{\prime}\mathbf{k}}\rangle$
and
$[H_{\mathbf{k}}(t)]_{\mu\mu^{\prime}}=\langle\phi_{\mu\mathbf{k}}|\hat{H}(t)|\phi_{\mu^{\prime}\mathbf{k}}\rangle$
are the overlap and hamiltonian matrices, respectively, and
$B_{\mathbf{k}}(t)$ reflects the time evolution of the basis:

$B_{\mathbf{k}}(t)=\sum_{J}\dot{\mathbf{R}}_{J}\cdot\mathcal{B}_{J\mathbf{k}}(t),\qquad[\mathcal{B}_{J\mathbf{k}}(t)]_{\mu^{\prime}\mu}=\left\langle\phi_{\mu^{\prime}\mathbf{k}}\left|\frac{\partial\phi_{\mu\mathbf{k}}}{\partial\mathbf{R}_{J}}\right.\right\rangle.$
(6)

The auxiliary matrix $\mathcal{B}_{J\mathbf{k}}$ is obtained by considering
that the changes in the LAPW terms due to the movement of the nuclei are [27]:

$\left|\frac{\partial\phi_{\mathbf{k}+\mathbf{G}}}{\partial\mathbf{R}_{J}}\right\rangle=\mathrm{i}(\mathbf{k}+\mathbf{G})\left|\phi_{\mathbf{k}+\mathbf{G}}\right\rangle_{MT_{J}}-\left|\nabla\phi_{\mathbf{k}+\mathbf{G}}\right\rangle_{MT_{J}},$
(7)

where the subscript $MT_{J}$ indicates that the corresponding function is non-
zero only inside the MT sphere centered around the nucleus $J$. The analogous
expression for the lo’s reads [27]:

$\left|\frac{\partial\phi_{\nu}}{\partial\mathbf{R}_{J}}\right\rangle=-\left|\nabla\phi_{\nu}\right\rangle_{MT_{J}}.$
(8)

Inserting Eqs. (7) and (8) into Eq. (6), and considering that $\mu$ and
$\mu^{\prime}$ can be LAPWs or lo’s, leads to:

$[\mathcal{B}_{J\mathbf{k}}(t)]_{\mu^{\prime}\mu}=\left\\{\begin{array}[]{l}\phantom{-}\mathrm{i}(\mathbf{G}+\mathbf{k})\langle\phi_{\mathbf{G}^{\prime}+\mathbf{k}}|\phi_{\mathbf{G}+\mathbf{k}}\rangle_{MT_{J}}-\mathrm{i}\langle\phi_{\mathbf{G}^{\prime}+\mathbf{k}}|\mathbf{p}|\phi_{\mathbf{G}+\mathbf{k}}\rangle_{MT_{J}}\\\
-\mathrm{i}\langle\phi_{\mathbf{G}^{\prime}+\mathbf{k}}|\mathbf{p}|\phi_{\mathbf{G}+\mathbf{\nu}}\rangle_{MT_{J}}\\\
\phantom{-}\mathrm{i}(\mathbf{G}+\mathbf{k})\langle\phi_{\nu^{\prime}}|\phi_{\mathbf{G}+\mathbf{k}}\rangle_{MT_{J}}-\mathrm{i}\langle\phi_{\nu^{\prime}}|\mathbf{p}|\phi_{\mathbf{G}+\mathbf{k}}\rangle_{MT_{J}}\\\
-\mathrm{i}\langle\phi_{\nu^{\prime}}|\mathbf{p}|\phi_{\nu}\rangle_{MT_{J}}\end{array}\right.$
(9)

where $\mathbf{p}=-\mathrm{i}\nabla$ is the momentum operator. The inner
products with subscripts $MT_{J}$ indicate that the integrals which they
contain are evaluated only inside $MT_{J}$.

###  2.2 Equation of nuclear motion

Extremizing Eq. (1) with respect to $\mathbf{R}_{J}$ leads to the following
equation of motion for the nuclei:

$M_{J}\ddot{\mathbf{R}}_{J}(t)=\mathbf{F}_{ext,J}(\mathbf{R}_{J},t)+\mathbf{F}_{HF,J}(\mathbf{R}_{J},t)+\mathbf{F}_{corr,J}(\mathbf{R}_{J},t).$
(10)

The external force $\mathbf{F}_{ext,J}$ arises from the interaction with the
applied electric field:

$\mathbf{F}_{ext,J}(\mathbf{R}_{J},t)=Z_{J}\mathbf{E}(t)=-\frac{Z_{J}}{c}\frac{d\mathbf{A}}{dt}(t).$
(11)

The Hellmann-Feynman force, $\mathbf{F}_{HF,J}$, can be expressed as [28]

$\displaystyle\mathbf{F}_{HF,J}(\mathbf{R}_{J},t)$ $\displaystyle=$
$\displaystyle-\int\mathrm{d}\mathbf{r}\,n(\mathbf{r},t)\frac{\partial
v_{nucl}}{\partial\mathbf{R}_{J}}$ (12) $\displaystyle=$ $\displaystyle
Z_{J}\lim_{r_{J}\to
0}\frac{1}{r_{J}}\sum_{m=-1}^{1}\nabla[r_{J}v_{C,lm}(r_{J})Y_{lm}(\hat{r}_{J})]_{l=1},$

where $\mathbf{r}_{J}=\mathbf{r}-\mathbf{R}_{J}$, $v_{C,lm}(r_{J})$ is the
radial function that is multiplied by $Y_{lm}$ when $v_{C}=v_{nucl}+v_{H}$ is
expanded in spherical harmonics, and the gradient above is taken with respect
to $\mathbf{r}_{J}$.

While the Hellmann-Feynman force is common to implementations of any basis-set
type, so-called Pulay corrections need to be considered in the LAPW+lo
formalism [27, 28]. The additional term, $\mathbf{F}_{corr,J}$, arising due to
the basis incompleteness, is a sum of contributions from the valence
($\mathbf{F}_{val,J}$) and core ($\mathbf{F}_{core,J}$) states. Neglecting
second-order terms with respect to the displacement, one arrives at:

$\displaystyle\mathbf{F}_{core,J}(\mathbf{R}_{J},t)$ $\displaystyle=$
$\displaystyle-\int_{MT_{J}}\mathrm{d}\mathbf{r}\,n_{core}(\mathbf{r},t)\nabla
v_{KS}(\mathbf{r},t)$ (13) $\displaystyle=$
$\displaystyle\int_{MT_{J}}\mathrm{d}\mathbf{r}\,v_{KS}(\mathbf{r},t)\nabla
n_{core}(\mathbf{r},t),$ $\displaystyle\mathbf{F}_{val,J}(\mathbf{R}_{J},t)$
$\displaystyle=$
$\displaystyle-\sum_{j\mathbf{k}}w_{\mathbf{k}}f_{j\mathbf{k}}(C_{j\mathbf{k}})^{\dagger}(\mathcal{H}_{J\mathbf{k}}-\mathcal{S}_{J\mathbf{k}})C_{j\mathbf{k}}$
(14) $\displaystyle+$ $\displaystyle\int_{MT_{J}}\mathrm{d}\mathbf{r}(\nabla
n_{val}(\mathbf{r}))\left(v_{KS}(\mathbf{r})+\frac{\mathbf{A}^{2}}{2c^{2}}\right),$

$n_{core}$ and $n_{val}$ are the density of core and valence electrons,
respectively. The auxiliary matrices $\mathcal{S}_{J\mathbf{k}}$ and
$\mathcal{H}_{J\mathbf{k}}$ are defined as:

$\mathcal{S}_{J\mathbf{k}}(t)=H_{\mathbf{k}}(t)S_{\mathbf{k}}^{-1}(t)\mathcal{B}_{J\mathbf{k}}(t)+\mathcal{B}_{J\mathbf{k}}^{\dagger}(t)S_{\mathbf{k}}^{-1}(t)H_{\mathbf{k}}(t),$
(15)
$[\mathcal{H}_{J\mathbf{k}}(t)]_{\mu^{\prime}\mu}=\left\\{\begin{array}[]{l}\phantom{-}\mathrm{i}(\mathbf{G}-\mathbf{G}^{\prime})\left[\langle\phi_{\mathbf{G}^{\prime}+\mathbf{k}}|\hat{H}(t)|\phi_{\mathbf{G}+\mathbf{k}}\rangle_{MT_{J}}-\lambda_{J}\right]\\\
-\mathrm{i}(\mathbf{G}^{\prime}+\mathbf{k})\langle\phi_{\mathbf{G}^{\prime}+\mathbf{k}}|\hat{H}(t)|\phi_{\nu}\rangle_{MT_{J}}\\\
\phantom{-}\mathrm{i}(\mathbf{G}+\mathbf{k})\langle\phi_{\nu^{\prime}}|\hat{H}(t)|\phi_{\mathbf{G}+\mathbf{k}}\rangle_{MT_{J}}\\\
\phantom{-}0\mbox{ for the case lo-lo}\end{array}\right.$ (16)

The definition of the auxiliary function $\lambda_{J}$ is

$\lambda_{J}=\left(\frac{\mathbf{G}^{\prime}+\mathbf{k}}{2}+\frac{\mathbf{A}(t)}{c}\right)\cdot(\mathbf{G}+\mathbf{k})\frac{4\pi
R_{MT_{J}}^{3}}{\Omega}\frac{j_{1}(|\mathbf{G}-\mathbf{G}^{\prime}|R_{MT_{J}})}{|\mathbf{G}-\mathbf{G}^{\prime}|R_{MT_{J}}}\mathrm{e}^{\mathrm{i}\mathbf{(\mathbf{G}-\mathbf{G}^{\prime})}\cdot\mathbf{R}_{J}}$
(17)

where $\Omega$ is the unit cell volume.

###  2.3 Different time steps

As suggested in Refs. [10] and [29], we allow for two different time steps. A
smaller one, $\Delta t$, is taken to solve Eq. (5), which dictates the
dynamics of the electrons. The positions of the nuclei are kept fixed and are
only changed in periods of $\alpha\Delta t$, where $\alpha$ is a positive
integer. This is physically grounded on the larger masses of the nuclei, when
compared to electrons. At the same time, this procedure permits faster
calculations, since the evaluation of the forces on the ions, specially the
correction $\mathbf{F}_{val,J}(\mathbf{R}_{J},t)$, can be slow. In Section 3,
we will address the convergence behavior with respect to $\alpha$.

##  3 Examples of applications

###  3.1 BN in a distorted lattice

Our first example is BN in zinc-blend geometry with the experimental lattice
constant of $a=6.833$ bohr [30]. An $8\times 8\times 8$ k-grid and a basis-set
cutoff of rgkmax=6.0 are employed in the calculation. As initial condition for
the MD, nitrogen is displaced from $(0.25,0.25,0.25)$ to $(0.26,0.26,0.26)$
(in lattice coordinates), while boron is kept at the origin of the unit cell.
This gives rise to a restoring force, making the atoms oscillate around their
equilibrium positions. For the MD, we employ a time step of 0.01 a.u. for
electrons and 0.05 a.u. for the nuclei. The approximated enforced time-
reversal symmetry (AETRS) method together with a Taylor expansion up to the
4th order [23] are used as propagator for solving Eq. (5).

Figure 1: Oscillations of the distance $d$ between boron and nitrogen with
respect to the equilibrium bond length $d_{0}$ – comparison between exciting
and Octopus.

In Fig. 1, we depict in blue the evolution of the distance $d$ between boron
and nitrogen with respect to the equilibrium bond length $d_{0}$. In gray are
the results obtained with Octopus [31] employing a spacing of $a/40$ for the
real-space grid. There is very good agreement between the two codes, and the
difference in the period of oscillations, provided in Table 1, is less than
$0.37$%.

Table 1: Period $T$ (in a.u.) of oscillations in cubic BN (Fig. 1) and diamond (Fig. 3) obtained by MD. GS indicates the expected period from a separate set of ground-state calculations (see text).  |  $T$ [a.u.]
---|---
| BN | C
exciting \- MD | 1349 | 1239
Octopus | 1354 | 1242
exciting \- GS | 1294 | 1091

In this Table, we also provide the period predicted through
$T=2\pi\sqrt{\mu/k}$, where $\mu$ is the reduced mass of boron and nitrogen,
and $k$ is the spring constant obtained from a set of separate ground-state
(GS) calculations in which B and N are displaced from each other, and the
total energy with respect to the displacement is fitted with a parabola. The
period predicted from these GS calculations differs by approximately 4.0% from
the other two results. We attribute this difference to non-linear effects.
Even though the initial condition of placing N at $(0.26,0.26,0.26)$ may seem
a small perturbation, it is enough to reach the non-linear regime. This can be
seen in Fig. 1, since the oscillations are not symmetric with respect to the
equilibrium position.

###  3.2 Diamond in a distorted lattice

As a second example, we consider diamond with a distortion as initial
condition for the MD, similar to the previous subsection. We employ the
experimental lattice constant $a=6.7407$ bohr [32]. Like in the example of BN,
we displace one carbon atom to a non-equilibrium position, in this case
$(0.30,0.30,0.30)$, while keeping the other C atom at the origin. The
restoring force is expected to be larger than is the previous case and to
drive further in the non-linear regime. At first, we use the same time step of
0.05 a.u. for electrons and nuclei. The AETRS method together with a Taylor
expansion up to the 4th order [23] are used as propagator for solving Eq. (5).

Figure 2: Left: Impact of the basis set size, measured by the parameter
rgkmax, on the oscillations of the carbon-carbon distance $d$ relative to the
equilibrium bond length $d_{0}$. Right: Absolute value of the contributions to
the force acting on the displaced carbon atom over time (rgkmax=6.0).

Our calculations employ an $8\times 8\times 8$ k-grid, and rgkmax is varied
from 4.5 to 6.0, as shown in Fig. 2. We observe that starting from an rgkmax
of 5.0, the curves are almost indistinguishable. On the right side of Fig. 2,
for the case rgkmax=6.0, we compare the absolute value of the forces given in
Eq. (10): $\mathbf{F}_{HF,J}(\mathbf{R}_{J},t)$ and the correction
$\mathbf{F}_{corr,J}(\mathbf{R}_{J},t)$ split into
$\mathbf{F}_{core,J}(\mathbf{R}_{J},t)$ and
$\mathbf{F}_{val,J}(\mathbf{R}_{J},t)$. We see that the Hellmann-Feynman term
is much larger than the corrections. In this case, rgkmax=6.0 makes the LAPW
basis largely complete.

Figure 3 displays on the left the impact of the time-step multiplier $\alpha$
on the oscillations of the carbon atoms, when we consider time steps
$\alpha\Delta t$ and $\Delta t$, respectively, for updating the nuclei
positions and the KS wavefunctions. On the right, the speedup as a function of
$\alpha$ is shown. It is apparent that factors of 5 or 10 are possible without
losing significant precision.

Figure 3: Left: Influence of the multiplicative factor $\alpha$ on the MD of
diamond, if we update the atomic positions with a time step $\alpha\Delta t$,
$\Delta t$ being the time step for the electronic dynamics within RT-TDDFT.
Right: Speedup as a function of $\alpha$; the dotted line indicates the ideal
speedup.

A calculation performed with Octopus using a real-space grid with spacing of
$a/25$ is shown on the left side of Fig. 3 for comparison and highlighting
remarkable agreement for $\alpha=1$. In Table 1, we show the period obtained
with exciting (same time step for nuclei and electrons) and Octopus, the
difference being less than $0.25$%. Considering, as in the case of BN, a set
of GS calculations to predict the period from $T=2\pi\sqrt{\mu/k}$, we observe
a higher discrepancy of 12%, owing to larger non-linear effects that are also
seen in the asymmetry of the oscillations around the equilibrium position
(Figs. 2 and 3).

###  3.3 BN exposed to an intense electric field

In our third example, we start from the groundstate of zinc-blend BN and apply
an electric field $\mathbf{E}(t)=-\frac{1}{c}\frac{d\mathbf{A}}{dt}$ along the
[111] direction, where $\mathbf{A}(t)=\mathbf{A}_{0}f(t)\cos(\omega_{0}t).$ To
be able to observe a non-negligible effect of the electric field on the ions,
we choose an angular frequency $\omega_{0}$=0.5 a.u and set the magnitude of
$\mathbf{A}_{0}$ to 10.0 a.u. in each cartesian direction. For the envelope
function $f(t)$, we consider $\sin^{2}(\pi t/T_{p})$ until a time $T_{p}=1000$
a.u., and zero afterwards. This function mimics a gaussian envelope, quite
common in experiments.

Figure 4: Different contributions to the force felt by the boron atom in BN.
Only the $x$ component is shown, since symmetry implies that the other
components are equivalent.

In Fig. 4, we display the $x$ component of the different contributions to the
force felt by boron as given in Eqs. (10)-(14): $F_{HF}$, $F_{core}$,
$F_{val}$, and $F_{ext}$. The amplitudes of all contributions are small
compared to the example in Section 3.2 (see Fig. 2), where the diamond lattice
is artificially heavily distorted, a situation difficult to be reproduced by
applying a laser. The external force $F_{ext}$ assumes, as expected, the shape
of the applied field. Rapid oscillations are also observable in $F_{val}$ that
explicitly depends on the vector potential $\mathbf{A}(t)$. These oscillations
still remain after the field is no longer acting on the system, as a result of
$\mathbf{A}_{0}$ being high enough to excite to the non-linear regime.
$F_{HF}$ varies smoothly, since the interaction with the other nuclei is
expected to follow the timescale of the phonon modes, which is on the order of
1000 a.u. [33]. The same applies to $F_{core}$: Since we disregard the effect
of $\mathbf{A}(t)$ on the core electrons, $F_{core}$ should arise mainly from
displacements of the nuclei, which carry the core electrons along.

In Fig. 5, we compare the force on boron obtained with exciting and Octopus.
Since the force oscillates rapidly as shown in Fig. 4, we display its Fourier
transform for better comparison. Overall, we find good agreement.

Figure 5: Comparison of exciting and Octopus: Fourier transform of the force
acting on boron (for $\omega_{0}$=0.5 a.u.).

##  4 Conclusions

In this paper, we have reported the implementation of Ehrenfest dynamics
combined with RT-TDDFT in the all-electron full-potential package exciting
that implements the LAPW+lo basis. To illustrate its performance, we have
considered three examples: diamond and cubic BN with an initial lattice
distortion, and cubic BN excited by an oscillating electric field. Comparison
of our results with the Octopus code shows good agreement. The code will be
freely available in the next release of exciting (fluorine), and all results
underlying this publication are available on the NOMAD Laboratory [34, 35,
36].

This work was supported by the German Research Foundation within project Nr.
182087777 - SFB 951 (B11) and the priority program SPP2196 “Perovskite
Semiconductors,” project Nr. 424709454.

## References

## References

  * [1] Alder B J and Wainwright T E 1959 The Journal of Chemical Physics 31 459–466
  * [2] Andrade X, Castro A, Zueco D, Alonso J L, Echenique P, Falceto F and Rubio A 2009 Journal of Chemical Theory and Computation 5 728–742
  * [3] Scherrer A, Agostini F, Sebastiani D, Gross E and Vuilleumier R 2017 Physical Review X 7 031035 ISSN 2160-3308
  * [4] Nelson T R, White A J, Bjorgaard J A, Sifain A E, Zhang Y, Nebgen B, Fernandez-Alberti S, Mozyrsky D, Roitberg A E and Tretiak S 2020 Chemical Reviews 120 2215–2287
  * [5] Curchod B F E, Rothlisberger U and Tavernelli I 2013 ChemPhysChem 14 1314–1340 ISSN 1439-7641
  * [6] Crespo-Otero R and Barbatti M 2018 Chemical Reviews 118 7026–7068 ISSN 0009-2665
  * [7] Pradhan E, Sato K and Akimov A V 2018 Journal of Physics: Condensed Matter 30 484002 ISSN 0953-8984
  * [8] You P, Chen D, Lian C, Zhang C and Meng S 2021 Wiley Interdisciplinary Reviews: Computational Molecular Science 11 ISSN 1759-0876
  * [9] Ojanperä A, Havu V, Lehtovaara L and Puska M 2012 The Journal of Chemical Physics 136 144103
  * [10] Wang F, Yam C Y, Hu L and Chen G 2011 The Journal of Chemical Physics 135 044126
  * [11] Andermatt S, Cha J, Schiffmann F and VandeVondele J 2016 Journal of Chemical Theory and Computation 12 3214
  * [12] Isborn C M, Li X and Tully J C 2007 The Journal of Chemical Physics 126 134307
  * [13] Miyamoto Y, Rubio A and Tománek D 2006 Phys. Rev. Lett. 97(12) 126104 URL https://link.aps.org/doi/10.1103/PhysRevLett.97.126104
  * [14] Krasheninnikov A V, Miyamoto Y and Tománek D 2007 Phys. Rev. Lett. 99(1) 016104 URL https://link.aps.org/doi/10.1103/PhysRevLett.99.016104
  * [15] Prezhdo O V 2009 Accounts of Chemical Research 42 2005–2016 pMID: 19888715
  * [16] Long R and Prezhdo O V 2011 Journal of the American Chemical Society 133 19240–19249 pMID: 22007727
  * [17] Long R and Prezhdo O V 2014 Nano Letters 14 3335–3341
  * [18] Jia W, An D, Wang L W and Lin L 2018 Journal of Chemical Theory and Computation 14 5645–5652 ISSN 1549-9618
  * [19] Lejaeghere K, Bihlmayer G, Björkman T, Blaha P, Blügel S, Blum V, Caliste D, Castelli I E, Clark S J, Dal Corso A, de Gironcoli S, Deutsch T, Dewhurst J K, Di Marco I, Draxl C, Dułak M, Eriksson O, Flores-Livas J A, Garrity K F, Genovese L, Giannozzi P, Giantomassi M, Goedecker S, Gonze X, Grånäs O, Gross E K U, Gulans A, Gygi F, Hamann D R, Hasnip P J, Holzwarth N A W, Iuşan D, Jochym D B, Jollet F, Jones D, Kresse G, Koepernik K, Küçükbenli E, Kvashnin Y O, Locht I L M, Lubeck S, Marsman M, Marzari N, Nitzsche U, Nordström L, Ozaki T, Paulatto L, Pickard C J, Poelmans W, Probert M I J, Refson K, Richter M, Rignanese G M, Saha S, Scheffler M, Schlipf M, Schwarz K, Sharma S, Tavazza F, Thunström P, Tkatchenko A, Torrent M, Vanderbilt D, van Setten M J, Van Speybroeck V, Wills J M, Yates J R, Zhang G X and Cottenier S 2016 Science 351 ISSN 0036-8075
  * [20] Gulans A, Kozhevnikov A and Draxl C 2018 Phys. Rev. B 97(16) 161105
  * [21] Gulans A, Kontur S, Meisenbichler C, Nabok D, Pavone P, Rigamonti S, Sagmeister S, Werner U and Draxl C 2014 Journal of Physics: Condensed Matter 26 363202
  * [22] The exciting code http://exciting-code.org/
  * [23] Pela R R and Draxl C 2021 Electronic Structure 3 037001
  * [24] Saalmann U and Schmidt R 1996 Zeitschrift für Physik D Atoms, Molecules and Clusters 38 153–163
  * [25] Kunert T and Schmidt R 2003 The European Physical Journal D - Atomic, Molecular, Optical and Plasma Physics 25 15–24
  * [26] Singh D J and Nordström L 2006 Planewaves, Pseudopotentials, and the LAPW Method 2nd ed (New York: Springer) ISBN 0-387-28780-9
  * [27] Yu R, Singh D and Krakauer H 1991 Phys. Rev. B 43(8) 6411–6422 URL https://link.aps.org/doi/10.1103/PhysRevB.43.6411
  * [28] Kohler B, Wilke S, Scheffler M, Kouba R and Ambrosch-Draxl C 1996 Computer Physics Communications 94 31–48 ISSN 0010-4655 (Preprint mtrl-th/9511002)
  * [29] Li X, Tully J C, Schlegel H B and Frisch M J 2005 The Journal of Chemical Physics 123 084106
  * [30] Knittle E, Wentzcovitch R M, Jeanloz R and Cohen M L 1989 Nature 337 349
  * [31] Tancogne-Dejean N, Oliveira M J T, Andrade X, Appel H, Borca C H, Le Breton G, Buchholz F, Castro A, Corni S, Correa A A, De Giovannini U, Delgado A, Eich F G, Flick J, Gil G, Gomez A, Helbig N, Hübener H, Jestädt R, Jornet-Somoza J, Larsen A H, Lebedeva I V, Lüders M, Marques M A L, Ohlmann S T, Pipolo S, Rampp M, Rozzi C A, Strubbe D A, Sato S A, Schäfer C, Theophilou I, Welden A and Rubio A 2020 The Journal of Chemical Physics 152 124119
  * [32] Hom T, Kiszenik W and Post B 1975 Journal of Applied Crystallography 8 457–458
  * [33] Karch K and Bechstedt F 1997 Physical Review B 56 7404–7415 ISSN 1098-0121
  * [34] Draxl C and Scheffler M 2019 Journal of Physics: Materials 2 036001 URL https://doi.org/10.1088/2515-7639/ab13bb
  * [35] Draxl C and Scheffler M 2018 MRS Bulletin 43 676–682
  * [36] Nomad repository, dataset: Md-exciting-octopus. https://dx.doi.org/10.17172/NOMAD/2021.12.28-1