# Shadow Energy Functionals and Potentials in Born-Oppenheimer Molecular
Dynamics

Anders M. N. Niklasson<EMAIL_ADDRESS>Christian F. A. Negre  Theoretical
Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545

###### Abstract

In Born-Oppenheimer molecular dynamics (BOMD) simulations based on density
functional theory (DFT), the potential energy and the interatomic forces are
calculated from an electronic ground state density that is determined by an
iterative self-consistent field optimization procedure, which in practice
never is fully converged. The calculated energies and the forces are therefore
only approximate, which may lead to an unphysical energy drift and
instabilities. Here we discuss an alternative shadow BOMD approach that is
based on a backward error analysis. Instead of calculating approximate
solutions for an underlying exact regular Born-Oppenheimer potential, we do
the opposite. Instead, we calculate the exact electron density, energies and
forces, but for an underlying approximate shadow BO potential energy surface.
In this way the calculated forces are conservative with respect to the
approximate shadow potential and generate accurate molecular trajectories with
long-term energy stability. We show how such shadow BO potentials can be
constructed at different levels of accuracy as a function of the integration
time step, $\delta t$, from the constrained minimization of a sequence of
systematically improvable, but approximate, shadow energy density functionals.
For each energy functional there is a corresponding ground state BO potential.
These pairs of shadow energy functionals and potentials are higher-level
generalizations of the original “$0$th-level” shadow energy functionals and
potentials used in extended Lagrangian BOMD [Eur. Phys. J. B 94, 164 (2021)].
The proposed shadow energy functionals and potentials are useful only within
this extended dynamical framework, where also the electronic degrees of
freedom are propagated as dynamical field variables together with the atomic
positions and velocities. The theory is quite general and can be applied to MD
simulations using approximate DFT, Hartree-Fock or semi-empirical methods, as
well as to coarse-grained flexible charge models.

first principles theory, electronic structure theory, molecular dynamics,
extended Lagrangian, self-consistent field, minimization, non-linear
optimization, Broyden, quasi-Newton method, Anderson mixing, Pulay mixing,
DIIS

††preprint: LA-UR-22-29595

##  I Introduction

The general notion of a shadow molecular dynamics provides a highly powerful
concept that helps us understand and design accurate and computationally
efficient simulation schemes Yoshida (1990); Grebogi  _et al._ (1990);
Toxvaerd (1994); Gans and Shalloway (2000); Bond and Leimkuhler (2007);
Toxvaerd  _et al._ (2012); Hammonds and Heyes (2020). The idea behind shadow
molecular dynamics is based on a backward error analysis. Instead of
calculating approximate forces and energies for an underlying exact potential
energy surface, it is often easier to calculate exact forces and energies, but
for an underlying approximate shadow potential (or shadow Hamiltonian). In
this way important physical properties of the simulated shadow dynamics such
as time-reversibility, the conservation of the total energy and the phase-
space area, can be fulfilled, because the forces of the shadow dynamics can be
generated exactly. In practice, shadow dynamics simulation methods are
therefore often both more accurate and computationally more efficient compared
to alternative techniques. In particular, their long-term accuracy and
stability are often superior.

The shadow dynamics terminology was originally introduced in the analysis and
explanation of the accuracy and long-term stability of symplectic or geometric
integration schemes such as the velocity Verlet algorithm in terms of a shadow
Hamiltonian Yoshida (1990); Grebogi  _et al._ (1990); Toxvaerd (1994). Here we
use the notion of a shadow molecular dynamics in the slightly more general
form that is associated with a backward error analysis. A shadow dynamics is
then generated, for example, when rapid changes or discontinuities from
cutoffs in the exact interatomic potential are smoothed out with an
approximate shadow potential for which we can calculate the exact forces and
use longer integration time steps Hammonds and Heyes (2021); Müser (2022).

Shadow molecular dynamics was originally introduced in the context of
classical molecular mechanics. More recently, the concept of a shadow dynamics
has been applied also to non-linear self-consistent field (SCF) theory in
quantum-mechanical Born-Oppenheimer molecular dynamics (QMD) simulations based
on extended Lagrangian Born-Oppenheimer molecular dynamics (XL-BOMD) Niklasson
_et al._ (2007); Niklasson (2008); Cawkwell and Niklasson (2012); Hutter
(2012); Lin  _et al._ (2014); Souvatzis and Niklasson (2014); Niklasson (2017,
2021a). The idea of a shadow molecular dynamics has been applied also to the
non-linear time-dependent dynamics of superfluidity Patrick Henning (2021), as
well as to flexible charge equilibration models Niklasson (2021b, a).

In this article we will revisit the construction of the approximate shadow
energy functionals and potentials used in XL-BOMD simulations and show how
their accuracy can be systematically improved to higher-orders as a function
of the integration time step, $\delta t$. It is important to note that these
shadow energy functionals and potentials are designed and useful only as parts
of molecular dynamics simulations within the framework of XL-BOMD, where also
the electronic degrees of freedom are propagated as extended dynamical
variables together with the atomic positions and velocities. The interatomic
forces calculated from the gradients of the shadow Born-Oppenheimer potential
are exact only in this dynamical setting. For static, non-dynamical systems,
the corresponding interatomic forces are only approximate, and in general not
even very accurate.

In regular QMD simulations Marx and Hutter (2000); Tuckerman (2010); Niklasson
(2021a) the Born-Oppenheimer potential and the interatomic forces are
calculated on-the-fly from the ground-state electronic structure, which is
determined from an iterative SCF optimization of some constrained non-linear
energy functional, that is given, for example, from Hartree-Fock or DFT
Roothaan (1951); McWeeny (1959); Mermin (1963); Hohenberg and Kohn (1964);
Kohn and Sham (1965); Mermin (1965); Parr and Yang (1989); Dreizler and Gross
(1990); Engel and Dreizler (2011). In practice the iterative SCF optimization
is never fully converged and always approximate. This may create small errors
in the ground state electron density, but these small errors can break time-
reversibility and lead to non-conservative forces. Accumulated over time, the
small errors from the approximate SCF optimization will therefore become
significant. Often the errors appear as an unphysical systematic drift in the
total energy, where the incompletely converged electronic structure behaves as
an artificial heat source or sink Remler and Madden (1990); Pulay and Fogarasi
(2004); Herbert and Head-Gordon (2005); Niklasson  _et al._ (2006); Kühne  _et
al._ (2007), which invalidates the QMD simulations. In the more recent
formulations of XL-BOMD the shadow Born-Oppenheimer potential is designed to
avoid the computational overhead and convergence errors in the iterative SCF
optimization.

In XL-BOMD the iterative SCF optimization procedure is avoided by including
the electronic degrees of freedom as extended dynamical variables, in the
spirit of Car-Parrinello molecular dynamics Car and Parrinello (1985);
Niklasson (2021a), in addition to the atomic positions and velocities.
However, in contrast to Car-Parrinello molecular dynamics, a constrained
optimization is still required to calculate the exact electronic ground state,
but the optimization is performed for an approximate shadow energy functional.
This optimization can be performed exactly in a single step and no iterative
process is needed. The ground state energy then defines the shadow Born-
Oppenheimer potential and the corresponding conservative forces. The ability
of XL-BOMD to avoid an iterative optimization and still generate exact
conservative forces is thus of great practical interest, both by reducing the
computational cost and by improving the accuracy and long-term stability of
the molecular dynamics simulations.

The shadow potential approximates the exact fully-converged regular Born-
Oppenheimer potential. In the original shadow potential formulation of XL-
BOMD, the error in the forces and the potential energies scale with the size
of the integration time step, $\delta t$, to the second, ${\cal O}(\delta
t^{2})$, and fourth order, ${\cal O}(\delta t^{4})$, respectively. However,
there seems to be no way to improve the order of the scaling. The only way to
boost the accuracy is to reduce the size of the integration time step. Here we
will show how higher-levels of accuracy in the forces and shadow Born-
Oppenheimer potentials can be achieved from the constrained minimization of a
sequence of systematically improvable, but approximate, shadow energy
functionals. For each energy functional there is a corresponding ground state
Born-Oppenheimer potential. The accuracies of these pairs of shadow
functionals and potentials are determined by the size of the integration time
step, $\delta t$. The increased level of accuracy is thus meaningful only in
the context of molecular dynamics simulations.

A higher-order accuracy in the shadow Born-Oppenheimer potential will often
improve the long-term stability of a QMD simulation. This is of particular
interest in QMD simulations of chemical systems that may have unsteady charge
solutions or chemical reactions, for example, where the electronic energy gap
between the Highest Occupied Molecular Orbital (HOMO) and the Lowest
Unoccupied Molecular Orbital (LUMO) is opening and closing along the molecular
trajectories.

The theory will be explained in terms of general Hohenberg-Kohn DFT Hohenberg
and Kohn (1964); Kohn and Sham (1965); Mermin (1965); Parr and Yang (1989);
Dreizler and Gross (1990); Engel and Dreizler (2011) and is applicable to a
broad range of methods, including Hartree-Fock and Kohn-Sham DFT Roothaan
(1951); McWeeny (1959); Mermin (1963); Hohenberg and Kohn (1964); Kohn and
Sham (1965); Mermin (1965); Parr and Yang (1989); Dreizler and Gross (1990);
Engel and Dreizler (2011), approximate DFT and semi-empirical methods Elstner
_et al._ (1998); Finnis  _et al._ (1998); et al. (2020); Dewar and Thiel
(1977); Dewar  _et al._ (1985); Stewart (2013); Bannwarth  _et al._ (2018);
Dral  _et al._ (2019); Malone  _et al._ (2020); Zhou  _et al._ (2020);
Bannwarth  _et al._ (2020), as well as to various coarse-grained polarizable
charge equilibration models Vesely (1977); Sprik and Klein (1988); Mortier
_et al._ (1986); Rappe and III (1991); Lamoureux and Roux (2003); Verstraelen
_et al._ (2013); Naserifar  _et al._ (2017); York and Yang (1996); Tabacchi
_et al._ (2002); Niklasson (2021b).

The article is outlined as follows. First we review the construction of the
“0th-level” shadow energy functional and Born-Oppenheimer potential used in
the original shadow potential formulation of XL-BOMD. We then describe how an
improved “1st-level” pair of shadow energy functional and Born-Oppenheimer
potential can be constructed. We then derive the equations of motion in an
adiabatic limit, where we assume that the extended electronic motion is rapid
compared to the slower moving nuclei. This is consistent with the underlying
Born-Oppenheimer approximation. Thereafter we discuss generalizations to
higher $m$th-level pairs of shadow energy functionals and potentials. The
integration of the equations of motion for the electronic degrees of freedom
is then explained, where we use a low-rank preconditioned Krylov subspace
approximation. To better understand the shadow functionals and potentials we
consider the relationship to the Harris-Foulkes functional Harris (1985);
Foulkes and Haydock (1989) in the static non-dynamical case for Kohn-Sham DFT.
We also apply our theory for the $1$st-level pairs of shadow energy
functionals and potentials to a simple flexible charge equilibration model
that corresponds to an orbital-free coarse-grained DFT. Thereafter, to
summarize the results, we present a pseudocode for XL-BOMD simulations using a
1st-level shadow energy functional and potential. We demonstrate the improved
scaling and ability to treat unstable chemical systems using the 1st-level
shadow energy functional and Born-Oppenheimer potential based on self-
consistent charge density functional tight-binding (SCC-DFTB) theory Harrison
(1980); Foulkes and Haydock (1989); Porezag  _et al._ (1995); Elstner  _et
al._ (1998); Finnis  _et al._ (1998); Frauenheim  _et al._ (2000); Koskinen
and Mäkinen (2009); Gaus  _et al._ (2011); Aradi  _et al._ (2015); et al.
(2020). At the end we give a brief summary and our conclusions.

##  II Generalized Shadow Functionals and Potentials

To present the pairs of energy functionals and Born-Oppenheimer potentials we
will use Hohenberg-Kohn density functional theory Hohenberg and Kohn (1964).
The corresponding Kohn-Sham expressions are generated by replacing the
universal energy functional with its orbital-dependent Kohn-Sham energy
functional Kohn and Sham (1965). Generalization to Hartree-Fock theory and
semi-empirical methods, as well as to coarse-grained orbital-free flexible
charge models, should be straightforward Niklasson (2021a).

###  II.1 Born-Oppenheimer Potential

In Hohenberg-Kohn DFT Hohenberg and Kohn (1964); Parr and Yang (1989);
Dreizler and Gross (1990); Engel and Dreizler (2011), the relaxed ground state
electron density, $\rho_{\rm min}({\bf r})$, is given from a constrained
minimization of an energy density functional, $E[{\bf R},\rho]$, over all
physically relevant electron densities, $\rho$, 111All physically relevant
electron densities determined by anti-symmetric electron wavefunctions
Dreizler and Gross (1990); Engel and Dreizler (2011) that integrates to the
total number of electrons, $N_{e}$, i.e.

$\rho_{\rm min}({\bf r})=\arg\min_{\rho}\left\\{E[{\bf
R},\rho]\left|\int\rho({\bf r})d{\bf r}=N_{e}\right.\right\\}.$ (1)

The DFT energy functional,

$E[\rho]\equiv E[{\bf R},\rho]=F[\rho]+\int V_{\rm ext}({\bf R,r})\rho({\bf
r})d{\bf r},$ (2)

includes a system-independent, non-linear, universal electron functional,
$F[\rho]$, and an energy term with an external potential, $V_{\rm ext}({\bf
R,r})$, which we here assume is from ions at the atomic positions, ${\bf
R}=\\{{\bf R}_{I}\\}$. The universal energy functional, $F[\rho]$, includes
all the electron-electron interactions and the kinetic energy term. To keep it
general, we may also assume ensemble generalizations where $F[\rho]$ accounts
for thermal effects, including the entropy contribution at finite electronic
temperatures Mermin (1963, 1965); Parr and Yang (1989); Engel and Dreizler
(2011); Pittalis  _et al._ (2011); Pribram-Jones  _et al._ (2014); Niklasson
(2021a). In the corresponding Kohn-Sham DFT the thermal effects introduces
fractional occupation numbers of the Kohn-Sham orbitals Mermin (1965); Parr
and Yang (1989), which is important to be able to describe, for example,
metallic systems at finite temperatures and to stabilize calculations of
systems with a small or vanishing electronic energy gap.

In the Born-Oppenheimer approximation Heitler and London (1927); Born and
Oppenheimer (1927); Marx and Hutter (2000); Tuckerman (2002) the Born-
Oppenheimer potential energy surface, $U({\bf R})$, is determined for the
fully relaxed electronic ground state, i.e.

$U({\bf R})=E[\rho_{\rm min}]+V_{\rm nn}({\bf R}),$ (3)

which includes the additional ion-ion repulsion energy term, $V_{\rm nn}({\bf
R})$. The motion of the atoms can then be generated by integrating Newton’s
equation of motion,

$M_{I}{\bf\ddot{R}}_{I}=-\nabla_{I}U({\bf R}),$ (4)

where $\\{M_{I}\\}$ are the atomic masses, one for each atom $I$, and the dots
denote the time derivatives.

In general, the calculation of the ground state density, $\rho_{\rm min}({\bf
r})$, requires some form of iterative optimization procedure or SCF approach,
because of the non-linearity of the universal energy functional, $F[\rho]$.
For example, in Kohn-Sham DFT the SCF optimization requires repeated
diagonalizations of the effective single-particle Kohn-Sham Hamiltonians. This
can cause a significant computational overhead and in practice the solution is
never fully converged and only approximate. Force terms that in general are
very difficult, if not impossible to calculate in practice, like

$\displaystyle\int\left(\delta E[\rho]\big{/}\delta\rho({\bf
r})\right)\left(\partial\rho({\bf r})\big{/}\partial{\bf
R}_{I}\right)\big{|}_{\rho\approx\rho_{\rm min}}d{\bf r},$ (5)

are therefore not vanishing exactly, because $\left(\delta
E[\rho]\big{/}\delta\rho({\bf r})\right)$ is vanishing only if $\rho({\bf
r})=\rho_{\rm min}({\bf r})$ 222Notice that exact ground state, $\rho_{\rm
min}({\bf r})$, and its approximate ground state solution,
$\rho\approx\rho_{\rm min}$, depend on ${\bf R}$, i.e. $\rho_{\rm min}({\bf
r})\equiv\rho_{\rm min}({\bf R},{\bf r})$, but we have dropped the explicit
${\bf R}$-dependencies in our simplified notation. Insufficiently converged
solutions for the electronic ground state density, and where the non-vanishing
force term in Eq. (5) is ignored, therefore lead to non-conservative forces
that may invalidate a molecular dynamics simulation Remler and Madden (1990);
Pulay and Fogarasi (2004); Herbert and Head-Gordon (2005); Niklasson  _et al._
(2006); Kühne  _et al._ (2007). Recent formulations of XL-BOMD were developed
to overcome these shortcomings Niklasson (2021a).

###  II.2 Zeroth-Level Shadow Functional and Born-Oppenheimer Potential

In the more recent formulations of XL-BOMD Niklasson (2021a), the energy
functional, $E[\rho]$ in Eq. (2), is approximated by a linearized shadow
energy functional,

${\cal E}^{(0)}[\rho,n^{(0)}]=E[n^{(0)}]+\int\frac{\delta
E[\rho]}{\delta\rho({\bf r})}\Big{|}_{n^{(0)}}\left(\rho({\bf r})-n^{(0)}({\bf
r})\right)d{\bf r},$ (6)

which is given by a linearization of $E[\rho]$ around some approximate $0$th-
level ground state density, $n^{(0)}({\bf r})\approx\rho_{\rm min}({\bf r})$.
More generally, we can create a $0$th-level shadow energy functional Niklasson
(2021b, a) by some approximation, where

${\cal E}^{(0)}[\rho,n^{(0)}]=E[n^{(0)}]+{\cal O}(|\rho-n^{(0)}|^{2}).$ (7)

This generalization is of particular interest in formulations of orbital-free
flexible-charge equilibration models. It allows more freedom in the
construction of the shadow energy functional, e.g. where parts of $E[\rho]$
are expanded to second order in $\rho$ to guarantee a unique ground state
solutionNiklasson (2021b, a). The corresponding $n^{(0)}$-dependent ground
state electron density, $\rho_{\rm min}[n^{(0)}]$, is then given by the
constrained minimization as in Eq. (1), where

$\rho_{\rm min}[n^{(0)}]({\bf r})=\arg\min_{\rho}\left\\{{\cal
E}^{(0)}[\rho,n^{(0)}]\left|\int\rho({\bf r})d{\bf r}=N_{e}\right.\right\\}.$
(8)

With the minimization we here mean the lowest stationary solution over all
physically relevant electron densities with $N_{e}$ number of electrons. The
relaxed ground state density then defines the approximate,
$n^{(0)}$-dependent, shadow Born-Oppenheimer potential,

${\cal U}^{(0)}({\bf R},n^{(0)})={\cal E}^{(0)}\left[\rho_{\rm
min}[n^{(0)}],n^{(0)}\right]+V_{\rm nn}({\bf R}).$ (9)

The advantage with this $0$th-level shadow energy functional, ${\cal
E}^{(0)}[\rho,n^{(0)}]$ in Eq. (6), is that the ground state density,
$\rho_{\rm min}[n^{(0)}]({\bf r})$, can be calculated without requiring any
iterative optimization procedure to find a SCF solution – at least if we have
found some appropriate shadow energy functional, ${\cal
E}^{(0)}[\rho,n^{(0)}]$, consistent with Eq. (7). Instead, the exact ground
state electron density can be calculated directly in a single step, because
all the non-linearities in $E[\rho]$ with respect to $\rho$ that would require
an iterative solution have been removed in ${\cal E}^{(0)}[\rho,n^{(0)}]$. In
Kohn-Sham density functional theory, the exact minimization is reached in a
single construction and diagonalization of the Kohn-Sham Hamiltonian, and in
the corresponding coarse-grained charge equilibration models Niklasson (2021b,
a) the relaxed ground state is given from the solution of a quasi-diagonal
system of linear equations, which has a simple direct analytical solution. In
this way, any possible convergence problems and associated inconsistencies
between the calculated ground state density, $\rho_{\rm min}[n^{(0)}]$, and
the shadow Born-Oppenheimer potential, ${\cal U}^{(0)}({\bf R},n^{(0)})$, are
avoided.

Because of the linearization in the energy functional the error in the shadow
Born-Oppenheimer potential is of second order in the residual function,
$f[n^{(0)}]({\bf r})=\rho_{\rm min}[n^{(0)}]({\bf r})-n^{(0)}({\bf r})$, i.e.

$\left|{\cal U}^{(0)}-U\right|\propto\left|\rho_{\rm
min}[n^{(0)}]-n^{(0)}\right|^{2}.$ (10)

The approximate density, $n^{(0)}$, therefore needs to be close to the relaxed
ground state density, $\rho_{\rm min}[n^{(0)}]$ or $\rho_{\rm min}$, to ensure
that the error in the approximate shadow potential is small 333If
$n^{(0)}=\rho_{\rm min}$ then $\rho_{\rm min}[n^{(0)}]=n^{(0)}$ and ${\cal
U}^{(0)}=U$.. Below we will show how this is achieved in QMD simulations by
propagating the approximate ground state density, $n^{(0)}$, as a dynamical
field variable within an extended Lagrangian formulation, where $n^{(0)}\equiv
n^{(0)}({\bf r},t)$ is propagated by a harmonic oscillator that is centered
around the optimized ground state density, $\rho_{\rm min}[n^{(0)}]({\bf r})$,
along the molecular trajectories. But before we present the extended
Lagrangian molecular dynamics scheme we will show how the $0$th-level shadow
energy functional and Born-Oppenheimer potential can be improved in accuracy.

###  II.3 First-Level Shadow Functional and Born-Oppenheimer Potential

The accuracy of the approximate $0$th-level shadow energy functional, ${\cal
E}^{(0)}[\rho,n^{(0)}]$ in Eq. (6) or Eq. (7), can be improved. However, a
straightforward expansion of $E[\rho]$ to higher orders in $\rho$ would not
help, because this would require some iterative solution to the constrained
minimization problem of a non-linear energy functional. Instead, we have to
improve the accuracy of the approximate $0$th-level energy functional without
loosing the linearity in $\rho$. We can achieve this by improving the estimate
of $n^{(0)}$ to be even closer to the exact ground state density, $\rho_{\rm
min}$, in Eq. (1). This can be accomplished with an updated and more accurate
density, $n^{(1)}({\bf r})$, which is given by a single Newton optimization
step (for multiple steps see III.3),

$\begin{array}[]{l}{\displaystyle n^{(1)}({\bf r})\equiv n^{(1)}[n^{(0)}]({\bf
r})=n^{(0)}({\bf r})}\\\ {}{}\\\ {\displaystyle~{}~{}~{}~{}~{}-\int
K^{(0)}({\bf r,r^{\prime}})\left(\rho_{\rm min}[n^{(0)}]({\bf
r^{\prime}})-n^{(0)}({\bf r^{\prime}})\right)d{\bf r^{\prime}}},\end{array}$
(11)

where the kernel $K^{(0)}({\bf r,r^{\prime}})$ is the inverse Jacobian of the
residual function, $f[n^{(0)}]({\bf r})=\rho_{\rm min}[n^{(0)}]({\bf
r})-n^{(0)}({\bf r})$. This means that

$\int K^{(0)}({\bf r,r^{\prime}})\frac{\delta\left(\rho_{\rm
min}[n^{(0)}]({\bf r^{\prime}})-n^{(0)}({\bf r^{\prime}})\right)}{{\delta
n^{(0)}({\bf r^{\prime\prime}})}}d{\bf r^{\prime}}=\delta({\bf
r-r^{\prime\prime}}).$ (12)

The Newton step in Eq. (11) (under reasonable conditions) is quadratically
convergent such that

$|\rho_{\rm min}-n^{(1)}|\propto|\rho_{\rm min}-n^{(0)}|^{2}\propto|\rho_{\rm
min}[n^{(0)}]-n^{(0)}|^{2}.$ (13)

We here assume that the functional is sufficiently well-behaved and that
$n^{(0)}$ is close enough to the exact ground state density, $\rho_{\rm min}$,
to achieve the quadratic convergence.

The shadow energy functional can now be improved in accuarcy by using the
updated density, $n^{(1)}$, instead of $n^{(0)}$ in the linearization of the
energy functional. This updated and improved approximate $1$st-level shadow
energy functional is then given by

$\begin{array}[]{l}{\displaystyle{\cal
E}^{(1)}\left[\rho,n^{(1)}\right]=E\left[n^{(1)}\right]}\\\ {}{}\\\
{\displaystyle~{}~{}~{}~{}+\int\frac{\delta
E\left[\rho]\right]}{\delta\rho({\bf r})}\Big{|}_{n^{(1)}}\left(\rho({\bf
r})-n^{(1)}({\bf r})\right)d{\bf r}},\\\ \end{array}$ (14)

or more generally as an approximation where

$\displaystyle E[\rho]$ $\displaystyle={\cal
E}^{(1)}\left[\rho,n^{(1)}\right]+{\cal O}(|\rho-n^{(1)}|^{2}).$ (15)

The updated $n^{(1)}$-dependent ground state density is then given from the
constrained minimization,

$\rho_{\rm min}[n^{(1)}]({\bf r})=\arg\min_{\rho}\left\\{{\cal
E}^{(1)}\left[\rho,n^{(1)}\right]\left|\int\rho({\bf
r})=N_{e}\right.\right\\},$ (16)

with respect to variationally stationary solutions. This ground state density
defines our $1$st-level shadow Born-Oppenheimer potential,

$\begin{array}[]{l}{\displaystyle{\cal U}^{(1)}({\bf R},n^{(0)})\equiv{\cal
U}^{(1)}({\bf R},n^{(1)}[n^{(0)}])}\\\ {}{}\\\ {\displaystyle={\cal
E}^{(1)}\left[\rho_{\rm min}[n^{(1)}],n^{(1)}\right]+V_{\rm nn}({\bf
R})}.\end{array}$ (17)

It is important to have the $1$st-level shadow potential, ${\cal U}^{(1)}({\bf
R},n^{(0)})$, expressed as a function of $n^{(0)}$ and not of $n^{(1)}$. We
can do so because $n^{(1)}$ is determined from $n^{(0)}$ in Eq. (11), where
$n^{(1)}\equiv n^{(1)}[n^{(0)}]$. We will take advantage of this relation in
the next section, where $n^{(0)}$ is propagated as a dynamical field variable,
$n^{(0)}({\bf r},t)$.

The constrained minimization in Eq. (16) can be achieved, in general, in a
single step without requiring any iterative optimization procedure, thanks to
the linear dependency of $\rho$ in the shadow energy functional, ${\cal
E}^{(1)}\left[\rho,n^{(1)}\right]$. No iterative self-consistent optimization
procedure is needed.

The error in the $1$st-level shadow potential scales as

$|{\cal U}^{(1)}-U|\propto|\rho_{\rm
min}[n^{(1)}]-n^{(1)}|^{2}\propto|\rho_{\rm min}[n^{(0)}]-n^{(0)}|^{4},$ (18)

thanks to the quadratic convergence of the Newton update of $n^{(1)}$ in Eq.
(11), where the size of the residual function, $f[n^{(0)}]({\bf r})=\rho_{\rm
min}[n^{(0)}]({\bf r})-n^{(0)}({\bf r})$, decays quadratically in a single
Newton step.

The Newton step in Eq. (11) is similar to an SCF iteration step. However, in
Kohn-Sham DFT, the Newton update does not require any additiontal Hamiltonian
diagonalization. In QMD simulations the Newton step can be performed using a
preconditioned Krylov subspace expansion Niklasson (2020a, b, 2021a); Das and
Gavini (2022), where each Krylov subspace vector can be determined from
response calculations using quantum perturbation theory. The preconditioned
Krylov subspace expansion used to approximate the kernel, $K^{(0)}$, acting on
the residual function Niklasson (2020a, b, 2021a); Das and Gavini (2022);
Negre  _et al._ (2022) is described in more detail in Sec. III.5.

In practice the preconditioned Krylov subspace expansion of the kernel,
$K^{(0)}$, is truncated and only approximate. The density update in Eq. (11)
is then given by a quasi-Newton step, which in general has a slower, non-
quadratic convergence.

###  II.4 Pairs of Shadow Functionals and Potentials

A key concept in our presentation are pairs of energy functionals and Born-
Oppenheimer potentials. The potential is always given from a constrained
minimization over the electron density of an energy functional, where the
initial pair of electronic energy functional and Born-Oppenheimer potential,
corresponding to regular DFT, is given by

$\Big{\\{}E[\rho],U({\bf R})\Big{\\}}.$ (19)

This pair in then replaced, first by the $n^{(0)}$-dependent $0$th-level
shadow energy functional and potential,

$\left\\{{\cal E}^{(0)}[\rho,n^{(0)}],{\cal U}^{(0)}({\bf
R},n^{(0)})\right\\},$ (20)

and then by the $n^{(0)}$-dependent $1$st-level shadow energy functional and
potential,

$\left\\{{\cal E}^{(1)}[\rho,n^{(0)}],{\cal U}^{(1)}({\bf
R},n^{(0)})\right\\}.$ (21)

The regular functional-potential pair in Eq. (19) is in practice difficult to
represent exactly, because the calculated Born-Oppenheimer potential, $U({\bf
R})$, at least in practice, is never given by the exact ground state of the
energy functional, $E[\rho]$. An accurate match between $E[\rho]$ and $U({\bf
R})$ can only be achieved by an expensive iterative optimization procedure,
because of the non-liniarity of $E[\rho]$. This is in contrast to the $0$th-
level shadow functional-potential pair in Eq. (20), which easily are matched
at only a modest cost, because no iterative optimization is required. What we
have presented so far is how we can construct an updated 1st-level pair of
shadow energy functionals and potentials in Eq. (21) that also can be matched
exactly. This higher-level generalization has an improved level of accuracy.

##  III Extended Lagrangian Born-Oppenheimer Molecular Dynamics

In a QMD simulation the initial approximate ground state density,
$n^{(0)}({\bf r})$, around which the linearization is performed for the
construction of the shadow energy functional will get further and further away
from the corresponding exact ground state density, $\rho_{\rm min}({\bf r})$,
as the atoms are moving away from the initial configuration. The accuracy of
the shadow energy functional and the corresponding shadow Born-Oppenheimer
potential will then get successively worse. The density, $n^{(0)}({\bf r})$,
therefore needs to be updated. One way is to update the density as a function
of the atomic positions, for example, where $n^{(0)}({\bf r})\equiv
n^{(0)}({\bf R,r})=\sum_{I}n^{\rm atom}_{I}({\bf r}-{\bf R}_{I})$, is the
superposition of separate neutral atomic electron denisities, $\\{n^{\rm
atom}_{I}({\bf r}-{\bf R}_{I})\\}$, centered around the atomic positions,
${\bf R}=\\{{\bf R}_{I}\\}$. However, this would lead to difficulties
calculating forces, as in Eq. (5), including all the density-dependent energy
terms – one for each atom. If $n^{(0)}({\bf r})$ would be the variational
ground state these terms would all vanish, but this is only true for the
optimized densities, $\rho_{\rm min}[n^{(0)}]({\bf r})$ or $\rho_{\rm
min}[n^{(1)}]({\bf r})$, with respect to the shadow potentials, ${\cal
U}^{(0)}({\bf R},n^{(0)})$ or ${\cal U}^{(1)}({\bf R},n^{(0)})$. Even in this
case partial derivatives, $\partial n^{(0)}/\partial{\bf R}_{I}$, would need
to be calculated. A solution to these problems is offered by XL-BOMD, where
$n^{(0)}$ is propagated as a dynamical field variable, $n^{(0)}({\bf r},t)$
Niklasson (2021a).

###  III.1 Extended Lagrangian

In XL-BOMD, we include the approximate ground state density, $n^{(0)}({\bf
r})$, and its time derivative as additional dynamical field variables,
$n^{(0)}({\bf r},t)$ and ${\dot{n}}^{(0)}({\bf r},t)$, in an extend Lagrangian
formalism, beside the nuclear positions and their velocities, ${\bf R}(t)$ and
${\bf\dot{R}}(t)$. The dynamics of $n^{(0)}({\bf r},t)$ is generated by an
extended harmonic oscillator that is centered around the optimized ground
state of the shadow potential, $\rho_{\rm min}[n^{(0)}]$, along the molecular
trajectories. In this way $n^{(0)}({\bf r},t)$ closely follows the ground
state such that the error in the shadow potential does not increase along the
trajectory.

In the Euler-Lagrange equations of motion, the partial derivatives only appear
with respect to each single dynamical variable, with all the other dynamical
variables being constant. The calculations of $n^{(0)}$-dependent force terms,
e.g.

$\displaystyle\int\frac{\delta{\cal U}^{(m)}}{\delta n^{(0)}({\bf
r})}\frac{\partial n^{(0)}({\bf r})}{\partial{\bf R}_{I}}d{\bf
r},~{}~{}(m=0~{}{\rm or}~{}1),$ (22)

can therefore be avoided. Additional force terms, such as

$\displaystyle\int\left(\frac{\delta{\cal
E}^{(m)}[\rho,n^{(m)}]}{\delta\rho({\bf
r})}\right)\left(\frac{\partial\rho({\bf r})}{\partial{\bf
R}_{I}}\right)\Big{|}_{\rho=\rho_{\rm min}[n^{(m)}]}d{\bf r},$ (23)

can also be ignored, because $\rho_{\rm min}[n^{(m)}]$ is determined from the
condition that

$\displaystyle\frac{\delta{\cal
E}^{(m)}[\rho,n^{(m)}]}{\delta\rho}\Big{|}_{\rho=\rho_{\rm
min}[n^{(m)}]}=0,~{}~{}(m=0~{}{\rm or}~{}1).$ (24)

This reduces not only the computational cost, but also makes it possible to
calculate “exact” conservative forces that generate stable long-term molecular
trajectories.

We can now define the $1$st-level extended Lagrangian, ${\cal L}^{(1)}$, in
XL-BOMD, using the $1$st-level shadow Born-Oppenheimer potential, where

$\begin{array}[]{l}{\displaystyle{\cal L}^{(1)}({\bf
R,\dot{R}},n^{(0)},{\dot{n}^{(0)}})=\frac{1}{2}\sum_{I}M_{I}|{\bf\dot{R}}_{I}|^{2}-{\cal
U}^{(1)}({\bf R},n^{(0)})}\\\ {}{}\\\
{\displaystyle+\frac{1}{2}\mu\int|{\dot{n}}^{(0)}({\bf r})|^{2}d{\bf
r}-\frac{1}{2}\mu\omega^{2}\iint\left(\rho_{\rm min}[n^{(0)}]({\bf
r})-n^{(0)}({\bf r})\right)}\\\ {}{}\\\ {\displaystyle\times T^{(0)}({\bf
r,r^{\prime}})\left(\rho_{\rm min}[n^{(0)}]({\bf r^{\prime}})-n^{(0)}({\bf
r^{\prime}})\right)d{\bf r}d{\bf r^{\prime}}}.\end{array}$ (25)

Here $n^{(0)}({\bf r},t)$ is treated as a dynamical field variable with its
time derivative, ${\dot{n}}^{(0)}({\bf r},t)$, and some chosen mass parameter,
$\mu$. This is in addition to the regular dynamical variables of the atomic
motion, ${\bf R}$ and ${\bf\dot{R}}$. The atomic masses are given by
$\\{M_{I}\\}$. The frequency of the extended harmonic oscillator is set by
$\omega$ and the harmonic well is centered around $\rho_{\rm
min}[n^{(0)}]({\bf r})$. $T^{(0)}({\bf r,r^{\prime}})$ is a symmetric positive
definite metric tensor given by the square of a kernel, $K^{(0)}({\bf
r,r^{\prime}})$, where

$T^{(0)}({\bf r,r^{\prime}})=\int\left(K^{(0)}({\bf
r,r^{\prime\prime}})\right)^{\dagger}K^{(0)}({\bf
r^{\prime\prime},r^{\prime}})d{\bf r^{\prime\prime}}.$ (26)

We define the kernel, $K^{(0)}({\bf r,r^{\prime}})$, as the inverse Jacobian
of the residual function, $f[n^{(0)}]({\bf r})=\rho_{\rm min}[n^{(0)}]({\bf
r})-n^{(0)}({\bf r})$. This means that the kernel, $K^{(0)}({\bf
r,r^{\prime}})$, is the same as in Eqs. (11) and (12). In this way the
dynamical density variable, $n^{(0)}({\bf r},t)$, evolves as if it would
oscillate around the much closer approximation to the exact ground state, i.e.
$n^{(1)}({\bf r})$ from the Newton update, compared to the more approximate,
$\rho_{\rm min}[n^{(0)}]({\bf r})$. This definition of the kernel simplifies
the equations of motion that we will derive below at the same time as it also
improves the accuracy of the shadow Born-Oppenheimer potential by evolving
$n^{(0)}({\bf r},t)$ around a closer approximation to the exact ground state,
$\rho_{\rm min}({\bf r})$.

The only difference to the original formulation of XL-BOMD is that the
Lagrangian, ${\cal L}^{(1)}({\bf R,\dot{R}},n^{(0)},{\dot{n}^{(0)}})$, in Eq.
(25) uses the $1$st-level shadow Born-Oppenheimer potential, ${\cal
U}^{(1)}({\bf R},n^{(0)})$, instead of the $0$-th level, ${\cal U}^{(0)}({\bf
R},n^{(0)})$.

###  III.2 Equations of Motion

The Euler-Lagranges equations for ${\cal L}^{(1)}({\bf
R,\dot{R}},n^{(0)},{\dot{n}^{(0)}})$,

$\frac{d}{dt}\left(\frac{\partial{\cal L}^{(1)}({\bf
R,\dot{R}},n^{(0)},{\dot{n}^{(0)}}))}{\partial{\bf\dot{R}}_{I}}\right)=\frac{\partial{\cal
L}^{(1)}({\bf R,\dot{R}},n^{(0)},{\dot{n}^{(0)}}))}{\partial{\bf R}_{I}}$ (27)

and

$\frac{d}{dt}\left(\frac{\delta{\cal L}^{(1)}({\bf
R,\dot{R}},n^{(0)},{\dot{n}^{(0)}}))}{\delta{\dot{n}}^{(0)}({\bf
r})}\right)=\frac{\delta{\cal L}^{(1)}({\bf
R,\dot{R}},n^{(0)},{\dot{n}^{(0)}}))}{\delta n^{(0)}({\bf r})}$ (28)

give us the equations of motion,

$\begin{array}[]{l}{\displaystyle M_{I}{\bf\ddot{R}}_{I}=-\frac{\partial{\cal
U}^{(1)}({\bf R},n^{(0)})}{\partial{\bf R}_{I}}\Big{|}_{n^{(0)},n^{(1)}}}\\\
{}{}\\\ {\displaystyle~{}~{}-\int\frac{\delta{\cal U}^{(1)}({\bf
R},n^{(0)})}{\delta n^{(1)}({\bf r})}\frac{\partial n^{(1)}({\bf
r})}{\partial{\bf R}_{I}}\Big{|}_{n^{(0)}}d{\bf r}}\\\ {}{}\\\
{\displaystyle~{}~{}~{}-\frac{1}{2}\mu\omega^{2}\frac{\partial}{\partial{\bf
R}_{I}}\iint\left(\rho_{\rm min}[n^{(0)}]({\bf r})-n^{(0)}({\bf
r})\right)T^{(0)}({\bf r,r^{\prime}})}\\\ {}{}\\\
{\displaystyle~{}~{}~{}~{}\times\left(\rho_{\rm min}[n^{(0)}]({\bf
r^{\prime}})-n^{(0)}({\bf r^{\prime}})\right)\Big{|}_{n^{(0)}}d{\bf r}d{\bf
r^{\prime}}}\\\ \end{array}$ (29)

and

$\begin{array}[]{l}{\displaystyle\mu{\ddot{n}}^{(0)}({\bf
r})=-\frac{\delta{\cal U}^{(1)}({\bf R},n^{(0)})}{\delta n^{(0)}({\bf r})}}\\\
{}{}\\\ {\displaystyle-\frac{1}{2}\mu\omega^{2}\frac{\delta}{\delta
n^{(0)}({\bf r})}\iint\left(\rho_{\rm min}[n^{(0)}]({\bf
r^{\prime}})-n^{(0)}({\bf r^{\prime}})\right)}\\\ {}{}\\\
{\displaystyle~{}\times T^{(0)}({\bf
r^{\prime},r^{\prime\prime}})\left(\rho_{\rm min}[n^{(0)}]({\bf
r^{\prime\prime}})-n^{(0)}({\bf r^{\prime\prime}})\right)d{\bf
r^{\prime}}d{\bf r^{\prime\prime}}}\end{array}$ (30)

These equations are far from trivial to use in a QMD simulation. However, the
equations of motion are simplified if we impose an adiabatic limit, in the
same way as for the original Born-Oppeheimer approximation, where we assume
that the electronic degrees of freedom are fast compared to the slower nuclear
motion. To derive the equations of motion in this adiabatic limit we first
assert the following frequency dependencies in the residual functions,

$\left|\rho_{\rm min}[n^{(0)}]-n^{(0)}\right|\propto\omega^{-2},$ (31)

and

$\left|\rho_{\rm min}[n^{(1)}]-n^{(1)}\right|\propto\omega^{-4},$ (32)

which are assumed to be valid in the limit of $\omega\rightarrow\infty$. These
adiabatic relations are difficult to prove a priori, but they can be shown to
hold a posteriori by integrating the equations of motions that have been
derived under the assumptions of Eqs. (31) and (32). This will demonstrate
below in Fig. 1.

Using the asserted adiabatic scaling relations in Eqs. (31) and (32), we find
(under reasonable conditions) from the definition of $n^{(1)}[n^{(0)}]({\bf
r})$ in Eq. (11) that

$\begin{array}[]{l}{\displaystyle\frac{\delta n^{(1)}[n^{(0)}]({\bf
r})}{\delta n^{(0)}({\bf r^{\prime\prime}})}=\frac{\delta n^{(0)}({\bf
r})}{\delta n^{(0)}({\bf r^{\prime\prime}})}}\\\ {}{}\\\ {\displaystyle-\int
K^{(0)}({\bf r,r^{\prime}})\frac{\delta\left(\rho_{\rm min}[n^{(0)}]({\bf
r^{\prime}})-n^{(0)}({\bf r^{\prime}})\right)}{{\delta n^{(0)}({\bf
r^{\prime\prime}})}}d{\bf r^{\prime}}}\\\ {}{}\\\
{\displaystyle-\int\frac{\delta K^{(0)}({\bf r,r^{\prime}})}{\delta n({\bf
r^{\prime\prime}})}\left(\rho_{\rm min}[n^{(0)}]({\bf
r^{\prime}})-n^{(0)}({\bf r^{\prime}})\right)d{\bf r^{\prime}}}\\\ {}{}\\\
{\displaystyle=\delta({\bf r-r^{\prime\prime}})-\delta({\bf
r-r^{\prime\prime}})}\\\ {}{}\\\ {\displaystyle-\int\frac{\delta K^{(0)}({\bf
r,r^{\prime}})}{\delta n({\bf r^{\prime\prime}})}\left(\rho_{\rm
min}[n^{(0)}]({\bf r^{\prime}})-n^{(0)}({\bf r^{\prime}})\right)d{\bf
r^{\prime}}}\\\ {}{}\\\ {\displaystyle\propto\left|\rho_{\rm
min}[n^{(0)}]({\bf r^{\prime}})-n^{(0)}({\bf
r^{\prime}})\right|\propto\omega^{-2}}.\end{array}$ (33)

Using the same assertions we also find that

$\left|\frac{\delta{\cal U}^{(1)}}{\delta
n^{(1)}}\right|\propto\left|\rho_{\rm
min}[n^{(1)}]-n^{(1)}\right|\propto\omega^{-4}.$ (34)

This gives us

$\displaystyle\left|\frac{\delta{\cal U}^{(1)}}{\delta n^{(0)}}\right|$
$\displaystyle=\left|\frac{\delta{\cal U}^{(1)}}{\delta n^{(1)}}\frac{\delta
n^{(1)}}{\delta n^{(0)}}\right|$ (35) $\displaystyle\propto\left|\rho_{\rm
min}[n^{(1)}]-n^{(1)}\right|\times\left|\rho_{\rm
min}[n^{(0)}]-n^{(0)}\right|$ (36)
$\displaystyle\propto\omega^{-4}\times\omega^{-2}.$ (37)

The scaling relation in Eq. (31) also mean that the last gradient term in Eq.
(29) becomes proportional to $\mu$. The asserted scaling relations above
inserted in the equations of motion in Eqs. (29) and (30) then give us,

$\begin{array}[]{l}{\displaystyle M_{I}{\bf\ddot{R}}_{I}=-\frac{\partial{\cal
U}^{(1)}({\bf R},n^{(0)})}{\partial{\bf R}_{I}}\Big{|}_{n^{(0)}}+{\cal
O}\left(\omega^{-4}\right)+{\cal O}\left(\mu\right)},\\\ \end{array}$ (38)

and

$\begin{array}[]{l}{\displaystyle{\ddot{n}}^{(0)}({\bf r})={\cal
O}(\mu^{-1}\omega^{-6})+{\cal O}(\omega^{-2})}\\\ {}{}\\\
{\displaystyle~{}~{}~{}~{}-\omega^{2}\int K^{(0)}({\bf
r,r^{\prime}})\left(\rho_{\rm min}[n^{(0)}]({\bf r^{\prime}})-n^{(0)}({\bf
r^{\prime}})\right)d{\bf r^{\prime}}},\end{array}$ (39)

where we have assumed that $\delta T^{(0)}/\delta n^{(0)}$ is bounded and
$\omega$-independent as $\omega\rightarrow\infty$. We can then derive the
equations of motion in the adiabatic limit, where $\omega\rightarrow\infty$
combined with the mass-zero limit $\mu\rightarrow 0$, which here is chosen
such that $\mu\omega^{4}\rightarrow{\rm constant}$. This is a classical
analogue to the Born-Oppenheimer approximation, where we simply stick with the
original Born-Oppenheimer assumption that the electronic degrees of freedom is
evolving on a much faster time scale compared to a slower nuclear motion. In
this adiabatic limit we get the final equations of motion for XL-BOMD with the
$1$st-level updated shadow Born-Oppenheimer potential,

$\displaystyle M_{I}{\bf\ddot{R}}_{I}=$
$\displaystyle-\frac{\partial}{\partial{\bf R}_{I}}{\cal U}^{(1)}({\bf
R},n^{(0)})\Big{|}_{n^{(0)}},$ (40) $\displaystyle{\ddot{n}}^{(0)}({\bf r})=$
$\displaystyle-\omega^{2}\int K^{(0)}({\bf r,r^{\prime}})\left(\rho_{\rm
min}[n^{(0)}]({\bf r^{\prime}})-n^{(0)}({\bf r^{\prime}})\right)d{\bf
r^{\prime}}.$ (41)

Because $n^{(0)}\equiv n^{(0)}({\bf r},t)$ is a dynamical field variable in
XL-BOMD, the partial derivatives in Eq. (40) with respect to the nuclear
coordinates are evaluated under a constant electron density, $n^{(0)}$. Thus,
even if $n^{(0)}$ is not the variationally optimized ground state density, we
can still calculate the exact forces in the adiabatic equations of motion for
XL-BOMD. However, this does not work for static calculations. It only works in
the context of XL-BOMD, where the electronic degrees of freedom are propagated
dynamically. The 1st-level updated shadow potential is thus mainly useful only
in this dynamical setting.

The equations of motion, Eqs. (40) and (41), are almost identical to the
original equations of motion for XL-BOMD using the $0$th-level Born-
Oppenheimer potential Niklasson (2017, 2020a, 2021a). The only difference is
that we now have the $1$st-level Born-Oppenheimer potential, ${\cal
U}^{(1)}({\bf R},n^{(0)})$, instead of the original $0$th-level ${\cal
U}^{(0)}({\bf R},n^{(0)})$. The error terms neglected in the adiabatic limit,
where $\mu\propto\omega^{-4}$, indicates that the error in the interatomic
force term should scale as $\omega^{-4}$. This scaling will also be
demonstrated below in Fig. 2.

It is important to note that even if we only would use some approximation of
the kernel, $K({\bf r,r^{\prime}})$, the same equations of motion, in Eqs.
(40) and (41), can be derived in an adiabatic limit. The only difference is
that the adiabatic limit has to be modified such that
$\mu\omega^{m}\rightarrow{\rm constant}$ for some value $m\in[1,4]$, and with
a modified assertion, where $|\rho_{\rm
min}[n^{(1)}]-n^{(1)}|\propto\omega^{-m}$ for some value of $m\in[2,4]$. The
scaling of the errors in the forces and the potential energy will then be
different and less favorable. Of critical importance is only that we calculate
the forces from the shadow potential, ${\cal U}^{(1)}$, defined by the
optimized ground state of a shadow energy functional that has been linearized
around some updated $n^{(0)}$-dependent density, $n^{(1)}({\bf r})\equiv
n^{(1)}[n^{(0)}]({\bf r})$, and where $\big{|}\delta{\cal
U}^{(1)}\big{/}\delta n^{(1)}\big{|}\propto|\rho_{\rm min}[n^{(1)}]-n^{(1)}|$.
Replacing the Newton update of the electron density in Eq. (11) with an
approximate quasi-Newton scheme or any other SCF-like iteration update, should
therefore also work under the same conditions. This observation may also help
explain why some earlier versions of XL-BOMD Niklasson  _et al._ (2006);
Niklasson (2008); Steneteg  _et al._ (2010); Arita  _et al._ (2014); Aradi
_et al._ (2015); Peters  _et al._ (2017) often works quite well, but where a
few SCF cycles often were required in each time step prior to the force
evaluations, while the extended electronic degrees of freedom was propagated
dynamically. Our analysis here shows us why and when we can expect these
initial versions of XL-BOMD to work or fail. This insight appears analogous to
how, for example, solving a system of non-linear equations with some simple ad
hoc mixed iterations (which often works), can be replaced by a more
transparent and efficient conjugate gradient or Newton-based method. Once we
understand the theoretically more rigorous alternative we also understand why
and when the ad hoc method works and how it can be improved.

The equations of motion in Eqs. (40) and (41), in combination with the
definition of the $1$st-level shadow energy functional, Eq. (14), and the
Born-Oppenheimer potential, Eq. (17), are some of the key results of this
article.

###  III.3 Higher-Level Generalizations

Higher $m$th-level generalizations of the pairs of shadow energy functionals
and potentials can also be designed, where the approximate higher-level
density approximations to the exact ground state are updated with repeated
Newton steps,

$\begin{array}[]{l}{\displaystyle n^{(m)}({\bf r})\equiv n^{(m)}[n^{(0)}]({\bf
r})}\\\ {}{}\\\ {\displaystyle\equiv
n^{(m)}\left[n^{(m-1)}[\ldots[n^{(0)}]]\right]({\bf r})=n^{(m-1)}({\bf r})}\\\
{}{}\\\ {\displaystyle-\int K^{(m-1)}({\bf r,r^{\prime}})\left(\rho_{\rm
min}\left[n^{(m-1)}\right]({\bf r^{\prime}})-n^{(m-1)}({\bf
r^{\prime}})\right)d{\bf r^{\prime}}}.\end{array}$ (42)

The corresponding linearized $m$th-level shadow energy density functionals are
then given by,

$\begin{array}[]{l}{\displaystyle{\cal
E}^{(m)}\left[\rho,n^{(m)}\right]=E\left[n^{(m)}\right]}\\\ {}{}\\\
{\displaystyle+\int\frac{\delta E\left[\rho\right]}{\delta\rho({\bf
r})}\Big{|}_{n^{(m)}}\left(\rho({\bf r})-n^{(m)}({\bf r})\right)d{\bf
r}}.\end{array}$ (43)

The constrained electronic ground state optimization then gives us the ground
state density,

$\begin{array}[]{l}\rho_{\rm min}\left[n^{(m)}\right]({\bf r})\\\ {}{}\\\
~{}~{}=\arg\min_{\rho}\left\\{{\cal
E}^{(m)}\left[\rho,n^{(m)}\right]\left|\int\rho({\bf r})d{\bf
r}=N_{e}\right.\right\\},\end{array}$ (44)

which defines the $m$th-level shadow Born-Oppenheimer potentials,

$\begin{array}[]{l}{\cal U}^{(m)}({\bf R},n^{(0)})={\cal
E}^{(m)}\left[\rho_{\rm min}\left[n^{(m)}\right],n^{(m)}\right]+V_{\rm
nn}({\bf R}).\end{array}$ (45)

The adiabatic equations of motion from an $m$th-level extended Lagrangian,
${\cal L}^{(m)}$, follows in the same way as above, where

$\displaystyle M_{I}{\bf\ddot{R}}_{I}=$
$\displaystyle-\frac{\partial}{\partial{\bf R}_{I}}{\cal U}^{(m)}({\bf
R},n^{(0)})\Big{|}_{n^{(0)}},$ (46) $\displaystyle{\ddot{n}}^{(0)}({\bf r})=$
$\displaystyle-\omega^{2}\int K^{(0)}({\bf r,r^{\prime}})\left(\rho_{\rm
min}[n^{(0)}]({\bf r^{\prime}})-n^{(0)}({\bf r^{\prime}})\right)d{\bf
r^{\prime}}.$ (47)

As for the $1$st-level approximation, we have used the nested dependencies of
$n^{(m)}$ on $n^{(0)}$ and let the shadow potential be a functional of
$n^{(0)}$. While the above higher-order generalization is straightforward, we
have found it of little value in practice, because the accuracy is, in
general, already very high at the $0$th-level and virtually exact at the
$1$st-level. For example, we tried to show numerically that the error in the
shadow potential energy surface, which scales as $\omega^{-4}$ for ${\cal
U}^{(0)}({\bf R},n^{(0)})$ Niklasson (2017), scales as $\omega^{-8}$ for
${\cal U}^{(1)}({\bf R},n^{(1)})$. However, in practical simulations this
scaling was not possible to observe, because the error in the $1$st-level
shadow potential for any normal integration time steps was already at machine
precision and no relevant scaling could be demonstrated. Instead, it has to be
demonstrated indirectly from the scaling of $|\rho_{\rm
min}[n^{(1)}]-n^{(1)}|$ from which we get the scaling of $|{\cal
U}^{(1))}-U|\propto|\rho_{\rm min}[n^{(1)}]-n^{(1)}|^{2}$. In the following we
will therefore ignore any higher-level generalizations beyond the $1$st-level.

###  III.4 Integrating the electronic equation of motion

To integrate the equations of motion for the nuclear degrees of freedom in Eq.
(40) we can use a leapfrog velocity Verlet scheme, whereas the integration of
the harmonic oscillator equation of motion in Eq. (41) for the extended
electronic degrees of freedom requires some care. In principle, the same
Verlet integration scheme could be used also for the electronic propagation.
However, typically we need to include some weak form of dissipation that keeps
$n^{(0)}({\bf r})$ synchronized with the trajectories of the atomic positions
and the exact Born-Oppenheimer ground state Niklasson  _et al._ (2009);
Steneteg  _et al._ (2010); Zheng  _et al._ (2011); Odell  _et al._ (2009,
2011); Niklasson (2021a). This modified Verlet integration scheme has the
following form,

${\displaystyle{\bf n}^{(0)}_{j+1}=2{\bf n}^{(0)}_{j}-{\bf
n}^{(0)}_{j-1}+\delta t^{2}{\bf\ddot{n}}^{(0)}_{j}+\alpha\sum_{l=0}^{l_{\rm
max}}c_{l}{\bf n}^{(0)}_{j-l}},$ (48)

where we use a convenient vector notation, with ${\bf n}^{(0)}_{j}\equiv{\bf
n}^{(0)}(t_{0}+j\delta t)\in{\boldmath R}^{N}$, $~{}j=0,1,2,\ldots$ . The
first three terms on the right-hand side of Eq. (48) are the regular Verlet
terms, whereas the last term is an additional weak dissipative ad hoc damping
force. An optimized set of coefficients of $\alpha$, $\\{c_{l}\\}$, and the
dimensionless constant $\kappa=\delta t^{2}\omega^{2}$, for various orders of
$l_{\rm max}$ can be found in Ref. Niklasson  _et al._ (2009). As an
alternative to such modified Verlet integration schemes, we may connect the
electronic degrees of freedom to a thermostat, i.e. a stochastic Langevin-like
dynamics or a chained Nose-Hoover thermostat, which also keeps the electronic
degrees of freedom synchronized with the ground state solution determined by
the nuclear coordinates Leven and Head-Gordon (2019); An  _et al._ (2020).

In the initial time step we can set all densities $\\{{\bf n}_{j}\\}$ equal to
the optimized regular Born-Oppenheimer ground state density,
${\boldsymbol{\rho}}_{\rm min}$.

It is important to note that we always use a constant for the product $\delta
t^{2}\omega^{2}=\kappa$ in our simulations. This means that $\delta
t\propto\omega^{-1}$, as long as we use the same Verlet integration scheme
with a constant size of the integration time step, $\delta t$. This controls
the way we can understand the scaling and the order of the accuracy, for
example, of the forces (See Fig. 2), as a function of the chosen size of the
integration time step, $\delta t$, or the inverse frequency, $\omega^{-1}$.

###  III.5 Approximating the kernel with preconditioned Krylov subspace

In addition to the modified Verlet integration, we also need to approximate
the kernel, $K({\bf r,r^{\prime}})$, both in the integration of the electronic
degrees of freedom in Eq. (41) and for the Newton update of the density,
$n^{(0)}$, to $n^{(1)}$ in Eq. (11). The kernel is the same and it is acting
on the same residual, apart from a trivial constant factor $\omega^{2}$. The
approximation of the kernel acting on the residual therefore only needs to be
performed once every integration time step.

In the more convenient matrix-vector notation, Eq. (41) or Eq. (47) is given
by

${\bf\ddot{n}}^{(0)}=-\omega^{2}{\bf K}\left({\boldsymbol{\rho}}^{(0)}_{\rm
min}[{\bf n}^{(0)}]-{\bf n}^{(0)}\right),$ (49)

where ${\bf K}\in{\boldmath R}^{N\times N}$, ${\bf K}={\bf J}^{-1}$,
${\boldsymbol{\rho}}^{(0)}_{\rm min}[{\bf n}^{(0)}]\in{\boldmath R}^{N}$, and
${\bf n}^{(0)}\in{\boldmath R}^{N}$. We can rewrite this equation of motion in
an equivalent preconditioned form,

${\bf\ddot{n}}^{(0)}=-\omega^{2}\left({\bf K}_{0}{\bf J}\right)^{-1}{\bf
K}_{0}\left({\boldsymbol{\rho}}^{(0)}_{\rm min}[{\bf n}^{(0)}]-{\bf
n}^{(0)}\right),$ (50)

where we have introduced a preconditioner, ${\bf K}_{0}\approx{\bf J}^{-1}$.
${\bf J}$ is the Jacobian of the residual function,

$\displaystyle{\bf f}({\bf n}^{(0)})=$
$\displaystyle{\boldsymbol{\rho}}^{(0)}_{\rm min}[{\bf n}^{(0)}]-{\bf
n}^{(0)}.$ (51)

If we use the notation,

$\displaystyle{\bf f}_{{\bf v}_{k}}({\bf n}^{(0)})\equiv$ $\displaystyle{\bf
K}_{0}\frac{\partial{\bf f}({\bf n}^{(0)}+\lambda{\bf
v}_{k})}{\partial\lambda}\Big{|}_{\lambda=0}={\bf K}_{0}{\bf J}{\bf v}_{k},$
(52)

it is possible to show that the preconditioned Jacobian, ${\bf K}_{0}{\bf J}$,
can be approximated by a low-rank (rank-$m$) approximation,

$\displaystyle{\bf K}_{0}{\bf J}\approx\sum_{kl}^{m}{\bf f}_{{\bf
v}_{k}}L_{kl}{\bf v}_{l}^{\rm T},$ (53)

for some set of vectors $\\{{\bf v}_{k}\\}$, and with ${\bf L}={\bf O}^{-1}$,
where $O_{ij}={\bf v}_{i}^{T}{\bf v}_{j}$ and $m<N$ Niklasson (2020a). The
directional derivatives of ${\bf f}({\bf n})$ in the direction of ${\bf
v}_{k}$ (or Gateaux derivatives) in Eq. (52) can be calculated using quantum
perturbation theory Niklasson (2020a); Niklasson  _et al._ (2015); Nishimoto
(2017).

The low-rank inverse of the preconditioned Jacobian, ${\bf K}_{0}{\bf J}$, is
then given by a pseudoinverse,

$\displaystyle\left({\bf K}_{0}{\bf J}\right)^{-1}\approx\sum_{kl}^{m}{\bf
v}_{k}{M_{kl}}{{\bf f}^{T}_{{\bf v}_{l}}},$ (54)

with ${\bf M}={\bf O}^{-1}$, where $O_{ij}={\bf f}_{{\bf v}_{i}}^{T}{\bf
f}_{{\bf v}_{j}}$. By chosing the vectors, $\\{{\bf v}_{k}\\}$, from an
orthogonalized preconditioned Krylov subspace Niklasson (2020a),

$\displaystyle\left\\{{\bf v}_{k}\right\\}\in{\rm span}^{\perp}\left\\{{\bf
K}_{0}{{\bf f}}({\bf n}^{(0)}),({\bf K}_{0}{\bf J})^{1}{\bf K}_{0}{{\bf
f}}({\bf n}^{(0)}),\right.$ (55) $\displaystyle\left.({\bf K}_{0}{\bf
J})^{2}{\bf K}_{0}{{\bf f}}({\bf n}^{(0)}),({\bf K}_{0}{\bf J})^{3}{\bf
K}_{0}{{\bf f}}({\bf n}^{(0)}),\ldots\right\\},$ (56)

we can rapidly reach a well-converged approximation of the kernel, ${\bf K}$,
acting on the residual function. The advantage with the preconditioner, ${\bf
K}_{0}$, is that it typically reduces the number of Krylov subspace vectors
(or low-rank updates) necessary to reach convergence. However, in principle
the preconditioner is not needed and the low-rank Krylov subspace
approximation works well also without preconditioning Das and Gavini (2022).

If we let $\Delta{\bf n}^{(0)}$ denote the result of the kernel acting on the
residual, i.e.

$\displaystyle\Delta{\bf n}^{(0)}=$ $\displaystyle\left({\bf K}_{0}{\bf
J}\right)^{-1}{\bf K}_{0}\left({\boldsymbol{\rho}}^{(0)}_{\rm min}[{\bf
n}^{(0)}]-{\bf n}^{(0)}\right)$ (57) $\displaystyle\approx$
$\displaystyle\left(\sum_{kl}{\bf v}_{k}{M_{kl}}{{\bf f}^{T}_{{\bf
v}_{l}}}\right){\bf K}_{0}{\bf f}({\bf n}^{(0)}),$ (58)

we find that the electronic equation of motion in Eq. (41) and the Newton step
in Eq. (11) are given by

$\displaystyle{\bf\ddot{n}}^{(0)}=-\omega^{2}\Delta{\bf n}^{(0)},$ (59)
$\displaystyle{\bf n}^{(1)}={\bf n}^{(0)}-\Delta{\bf n}^{(0)}.$ (60)

This clearly shows how the approximation of $\Delta{\bf n}^{(0)}$ only needs
to be performed once every time step for the $1$st-level generalized shadow
XL-BOMD in Eqs. (40) and (41). This simplification is another of our key
results.

The cost of calculating a preconditioner, ${\bf K}_{0}$, can be expensive, but
in QMD simulations the preconditioner can be reused, often over thousands of
integration time steps (or the whole simulation) before an updated
preconditioner is needed. In practice the overhead is therefore quite small.
In simulations of regular stable molecular systems a scaled delta function
typically works perfectly well as a preconditioner. The main cost of the
preconditioned subspace expansion required to approximate $\Delta{\bf
n}^{(0)}$ in Eq. (58) is therefore the construction of the residual response
vectors, $\\{{\bf f}_{{\bf v}_{k}}\\}$, from the directional perturbations in
$\\{{\bf v}_{k}\\}$. In Kohn-Sham DFT these response vectors can be calculated
from quantum perturbation theory Niklasson  _et al._ (2015); Nishimoto (2017).
If we assume that we have already performed a diagonalization of the
unperturbed Kohn-Sham Hamiltonian, ${\bf H}[n^{(0)}]$, to find $\rho_{\rm
min}[n^{(0)}]$ in Eq. (8), these response vectors are fairly easy to calculate
as no additional diagonalizations are needed Das and Gavini (2022); Negre  _et
al._ (2022). Nevertheless, the calculation of the residual response vectors,
$\\{{\bf f}_{{\bf v}_{k}}\\}$, is the main bottleneck of the Krylov subspace
expansion. In Kohn-Sham DFT, using an atomic-orbital basis, the cost is
typically dominated by the transformations back and forth between the non-
orthogonal atomic-orbital basis and the molecular-orbital eigenbasis, in which
the response summations are performed Niklasson (2020a, b).

###  III.6 The Harris-Foulkes functional

The $0$th and $1$st-level pairs of shadow energy functionals and Born-
Oppenheimer potentials presented in this article, e.g. as in Eqs. (6)-(9), are
quite general and easy to apply in different applications, e.g. to orbital-
based Kohn-Sham DFT, density matrix methods in Hartree-Fock theory, to
orbital-free polarizable charge equilibration models, or in a slightly
different form even to time-dependent models for superfluidity Niklasson
(2021a); Patrick Henning (2021). These linearized shadow energy functionals
and their constrained optimization may, of course, appear somewhat trivial. It
is only in combination with XL-BOMD that the shadow energy functionals and
Born-Oppenheimer potentials become useful and consequential. The proposed
$0$th or $1$st-level shadow energy functionals and potentials are only
meaningful within the dynamical simulation framework of XL-BOMD, where the
electronic degrees of freedom appear as dynamical field variables in addition
to the nuclear positions and velocities. Only then can we calculate the exact
forces from the shadow Born-Oppenheimer potential. For a static, non-dynamical
problem, the corresponding forces are only approximate, and would, in general,
require a well-converged iterative SCF optimization procedure to achieve any
reasonable accuracy.

For the static, non-dynamical problem, and for the particular case of Kohn-
Sham DFT, the optimized ground-state shadow Born-Oppenheimer potential,
although conceptually different, is interchangeable with the Harris-Foulkes
energy density functional Harris (1985); Foulkes and Haydock (1989); Niklasson
and Cawkwell (2014). The Harris-Foulkes (HF) functional, $E_{\rm
HF}[\rho_{0}]$, is an approximate energy expression for the electronic ground-
state energy in Kohn-Sham DFT that depends on some input density,
$\rho_{0}({\bf r})$, where

$\displaystyle E_{\rm HF}[\rho_{0}]=$
$\displaystyle\sum_{i}f_{i}\varepsilon_{i}-\frac{1}{2}\iint\frac{\rho_{0}({\bf
r})\rho_{0}({\bf r^{\prime}})}{|{\bf r-r^{\prime}}|}d{\bf r}d{\bf r^{\prime}}$
(61) $\displaystyle+E_{\rm xc}[\rho_{0}]-\int V_{\rm xc}[\rho_{0}]({\bf
r})\rho_{0}({\bf r})d{\bf r}.$ (62)

Here $\\{\varepsilon_{i}\\}$ are the eigenvalues of the Kohn-Sham Hamiltonian,
$H_{\rm KS}[\rho_{0}]$, calculated for the input density, $\rho_{0}({\bf r})$,
$\\{f_{i}\\}$ are the occupation numbers, $E_{\rm xc}[\rho_{0}]$ is the
exchange-correlation energy functional with the corresponding exchange-
correlation potential, $V_{\rm xc}[\rho_{0}]({\bf r})=\delta E_{\rm
xc}[\rho]\big{/}\delta\rho({\bf r})\big{|}_{\rho_{0}}$.

Apart from the nuclear-nuclear repulsion term (and possibly an additional
electronic entropy contribution), $E_{\rm HF}[\rho_{0}]$, has the same form as
the $0$th-level shadow potential, ${\cal U}^{(0)}({\bf R},n^{(0)})$, with
$\rho_{0}=n^{(0)}$. The difference is that ${\cal U}^{(0)}({\bf R},n^{(0)})$,
as it appears in XL-BOMD, represents an exact ground-state shadow Born-
Oppenheimer potential, which is determined from a variationally optimized
shadow energy functional, ${\cal E}^{(0)}[\rho,n^{(0)}]$, with some external
and electrostatic potentials that are given by the nuclear positions, ${\bf
R}(t)$, and a separate dynamical variable density, $n^{(0)}({\bf r},t)$.
Because $n^{(0)}({\bf r},t)$ is a dynamical field variable of the extended
Lagrangian in Eq. (25), forces in the Euler-Lagrange’s equations of motion can
easily be calculated from the partial derivatives of ${\cal U}^{(m)}({\bf
R},n^{(0)})$ with respect to a constant density, $n^{(0)}({\bf r},t)$. This is
in contrast to the Harris-Foulkes functional, which is an approximate energy
density functional expression for the (static) Kohn-Sham ground state energy,
where the density $\rho_{0}({\bf r})$ represents, either overlapping ${\bf
R}$-dependent atomic charge densities, or some iteratively and partially SCF
updated (and thus ${\bf R}$-dependent) input density. The Harris-Foulkes
energy functional is thus best used for estimating the electronic ground state
energy for approximate densities. Accurate calculations of the interatomic
forces would still require a regular iterative SCF optimization, or the
additional calculation of the gradients of the electron density with respect
to the atomic positions.

In XL-BOMD the pairs of shadow energy functionals and potential energy
surfaces therefore play a different role, and allows for computationally
simple and accurate calculations of conservative interatomic forces in
molecular dynamics simulations, without relying on the Hellmann-Feynman
theorem Marx and Hutter (2000); Feynman (1939). However, the linearized shadow
energy functionals and optimized Born-Oppenheimer potentials presented and
derived here, as in Eqs. (6)-(9), provide an alternative and probably more
transparent and straightforward approach to derive and understand the Harris-
Foulkes functional in Kohn-Sham DFT. The procedure in Eqs. (6)-(9) is also
easy to generalize and apply to a broad variety of other energy expressions
besides the Kohn-Sham energy functional Niklasson (2021a). As an example, in
the section below, we will use the approach in Eqs. (6)-(9) to the design a
shadow energy functional and potential for a coarse-grained flexible charge
equilibration model.

###  III.7 Coarse-grained flexible charge model

Flexible charge models can be derived from an atomic coarse-graining of DFT
Vesely (1977); Sprik and Klein (1988); Mortier  _et al._ (1986); Rappe and III
(1991); Lamoureux and Roux (2003); Naserifar  _et al._ (2017); York and Yang
(1996); Tabacchi  _et al._ (2002); Verstraelen  _et al._ (2013); Niklasson
(2021b, a). They often serve as simplified or conceptual versions of DFT and
can also be used to illustrate our shadow energy functionals in Born-
Oppenheimer molecular dynamics.

In the simplest form of flexible charge models the electronic energy
functional in DFT is approximated by the energy function

$\displaystyle E({\bf
R,q})=\sum_{I}\chi_{I}q_{I}+\frac{1}{2}\sum_{I}U_{I}q_{I}^{2}+\frac{1}{2}\sum_{IJ}^{I\neq
J}q_{I}\gamma_{IJ}q_{J},$ (63)

where ${\bf q}=\\{q_{I}\\}$ is the coarse-grained charge density, represented
by net partial charges (or electron occupations) of each atom $I$, $\chi_{I}$
are the estimated atomic electronegativities, $U_{I}$ the chemical hardness or
Hubbard-U parameters, and $\gamma_{IJ}$ describe the Coulomb interactions
between penetrating spherical atom-centered charge densities centered at atom
$I$ and $J$. At large interatomic distances these interactions decay as
$\gamma_{IJ}\rightarrow|{\bf R}_{I}-{\bf R}_{J}|^{-1}$ and at short-range
distances the onsite limit, $\gamma_{IJ}\rightarrow U_{I}$, is reached as
$|{\bf R}_{I}-{\bf R}_{J}|\rightarrow 0$.

The electronic ground state is given from the constrained minimization, where

$\displaystyle{\bf q}_{\rm min}=\min_{\bf q}\left\\{E({\bf
R,q})\left|\sum_{I}q_{I}=0\right.\right\\}.$ (64)

This minimization requires the solution of a full system of linear equations,
which is the main computational bottleneck. If an iterative solver is used the
optimized solutions need to be well-converged to provide accurate conservative
forces in a molecular dynamics simulation. The optimized ground state charges
then gives us the Born-Oppenheimer potential,

$\displaystyle U({\bf R})=E({\bf R},{\bf q}_{\rm min})+V({\bf R}).$ (65)

The molecular trajectories can then be generated from the integration of
Newton’s equations of motion,

$\displaystyle M_{I}{\bf\ddot{R}}_{I}=-\nabla_{I}U({\bf R}).$ (66)

Following the approach in Eqs. (6)-(9), a $0$th-level shadow energy function,
${\cal E}^{(0)}({\bf R,q},{\bf n}^{(0)})\approx E({\bf R,q})$, can be
constructed from a partial linearization of $E({\bf R,q})$ around some
approximate ground state solution, ${\bf n}^{(0)}\approx{\bf q}_{\rm min}$,
where

$\displaystyle{\cal E}^{(0)}({\bf R,q},{\bf n}^{(0)})=$
$\displaystyle\sum_{I}\chi_{i}q_{i}+\frac{1}{2}\sum_{I}U_{I}q_{I}^{2}$ (67)
$\displaystyle+\frac{1}{2}\sum_{I\neq
J}(2q_{I}-n_{I}^{(0)})\gamma_{IJ}n_{J}^{(0)}.$ (68)

The constrained minimization (the lowest stationary solution) of this shadow
energy function gives us the ${\bf n}^{(0)}$-dependent ground state density,

$\displaystyle{\bf q}_{\rm min}[{\bf n}^{(0)}]=\arg\min_{\bf q}\left\\{{\cal
E}^{(0)}({\bf R,q},{\bf n}^{(0)})\left|\sum_{I}q_{I}=0\right.\right\\}$ (69)

and the corresponding $0$th-level shadow Born-Oppenheimer potential,

$\displaystyle{\cal U}^{(0)}({\bf R},{\bf n}^{(0)})={\cal E}^{(0)}({\bf
R},{\bf q}_{\rm min}[{\bf n}^{(0)}],{\bf n}^{(0)})+V({\bf R}).$ (70)

The shadow energy function, ${\cal E}^{(0)}({\bf R,q},{\bf n}^{(0)})$, is
constructed such that ${\bf q}_{\rm min}[{\bf n}^{(0)}]$ is determined by a
quasi-diagonal system of linear equations that has a trivial analytical
solution Niklasson (2021b, a).

To introduce the 1st-level update we can improve the ground-state estimate of
${\bf n}^{(0)}$ with a Newton step,

$\displaystyle{\bf n}^{(1)}\equiv{\bf n}^{(1)}[{\bf n}^{(0)}]={\bf
n}^{(0)}-\Delta{\bf n}^{(0)},$ (71)

where

$\displaystyle\Delta{\bf n}^{(0)}=$ $\displaystyle\left({\bf K}_{0}{\bf
J}\right)^{-1}{\bf K}_{0}\left({\boldsymbol{\rho}}^{(0)}_{\rm min}[{\bf
n}^{(0)}]-{\bf n}^{(0)}\right),$ (72)

which can be approximated, for example, by the preconditioned low-rank Newton
step as in Eq. (58). Notice that this updated approximate charge vector is
${\bf n}^{(0)}$-dependent, i.e.

$\displaystyle{\bf n}^{(1)}\equiv{\bf n}^{(1)}[{\bf n}^{(0)}]={\bf
n}^{(0)}-\Delta{\bf n}^{(0)}.$ (73)

The updated $1$st-level energy function is now given by

$\displaystyle{\cal E}^{(1)}({\bf R,q},{\bf n}^{(1)})=$
$\displaystyle\sum_{I}\chi_{i}q_{i}+\frac{1}{2}\sum_{I}U_{I}q_{I}^{2}$ (74)
$\displaystyle+\frac{1}{2}\sum_{I\neq
J}(2q_{I}-n_{I}^{(1)})\gamma_{IJ}n_{J}^{(1)}.$ (75)

The optimized ground state density is then given from the constrained
minimization, where

$\displaystyle{\bf q}_{\rm min}[{\bf n}^{(1)}]=\arg\min_{\bf q}\left\\{{\cal
E}^{(1)}({\bf R,q},{\bf n}^{(1)})\left|\sum_{I}q_{I}=0\right.\right\\}.$ (76)

The shadow energy function ${\cal E}^{(1)}({\bf R,q},{\bf n}^{(1)})$ is
constructed in the same way as ${\cal E}^{(0)}({\bf R,q},{\bf n}^{(0)})$ such
that ${\bf q}_{\rm min}[{\bf n}^{(1)}]$ also is determined by a quasi-diagonal
system of linear equations that has a trivial analytical solution Niklasson
(2021b, a). This gives us the corresponding $1$st-level shadow Born-
Oppenheimer potential,

$\displaystyle{\cal U}^{(1)}({\bf R},{\bf n}^{(0)})={\cal E}^{(1)}({\bf
R},{\bf q}_{\rm min}[{\bf n}^{(1)}],{\bf n}^{(1)})+V({\bf R}),$ (77)

where ${\bf n}^{(1)}\equiv{\bf n}^{(1)}[{\bf n}^{(0)}]$. The shadow Born-
Oppenheimer potential can then be used in an extended Lagrangian formulation
Niklasson (2021b, a), which in an adiabatic limit gives us the equations of
motion,

$\displaystyle M_{I}{\bf\ddot{R}}_{I}=-\nabla_{I}{\cal U}^{(1)}({\bf R},{\bf
n}^{(0)})\big{|}_{{\bf n}^{(0)}},$ (78)
$\displaystyle{\bf\ddot{n}}^{(0)}=-\omega^{2}\Delta{\bf n}^{(0)}.$ (79)

The nuclear coordinates and velocities can then be integrated using a standard
velocity Verlet integration scheme and for the evolution of the atomic partial
charges, ${\bf n}(t)$, we can use the modified Verlet integation scheme
including some additional weak dissipative damping forces as in Eq. (48).

This example with a coarse-grained flexible charge equilibration model
demonstrates the general applicability of our shadow molecular dynamics
approach and how it can be used to construct pairs of shadow energy
functionals and potentials for XL-BOMD simulations at different levels of
accuracy.

##  IV Pseudocode

The easiest way to summarize the generalized $1$st-level update of the shadow
energy functional and Born-Oppenheimer potential in XL-BOMD is to describe the
method in a step-by-step procedure using a pseudocode. Algorithm 1 gives a
schematic picture of what an XL-BOMD simulation using the $1$st-level shadow
potential, ${\cal U}^{(1)}$, would look like for an orbital-dependent Kohn-
Sham like electronic structure theory. It is expressed in a matrix-vector
notation that is well-suited, for example, for SCC-DFTB simulations. All
$0$th-level superscript, (0), as in ${\bf n}^{(0)}$, have been dropped to
simplify the notation. Here ${\bf S}$ is a basis-set overlap matrix and ${\bf
H}$ is the effective single-particle (Kohn-Sham) Hamiltonian. Of critical
importance is the construction of $\Delta{\bf n}\equiv\Delta{\bf n}^{(0)}$
with a low-rank preconditioned Krylov subspace approximation using quantum
response calculations. In contrast to the most recent XL-BOMD schemes, we now
need two diagonalizations per time step, instead of only one. Algorithm 1
provides a compact summary of the most important results of this article.

Algorithm 1 Pseudocode for the XL-BOMD scheme using the $1$st-level updated
shadow Born-Oppenheimer potential, ${\cal U}^{(1)}({\bf R},n^{(0)})$. Matrix-
vector notation is used and the $0$th-level ${(0)}$-superscripts, i.e. as in
${\bf n}^{(0)}$ or $\Delta{\bf n}^{(0)}$, have been dropped for brevity. One
rank-$m$ approximation of $\Delta{\bf n}$ and two Hamiltonian diagonalizations
are required in each time step.

$\mbox{Atomic masses and positions,~{}}{\bf M}=\\{M_{I}\\},~{}{\bf R}=\\{{\bf
R}_{I}\\}$

$\mbox{Get ground state,~{}}{\bf q}_{\rm min},\mbox{~{}with regular SCF}$

${\bf q}_{\rm min}\Rightarrow\mbox{``exact''}~{}{U}({\bf R})~{}\mbox{and
forces},{\bf F}=\\{{\bf F}_{I}\\}$

$\mbox{Initialize charges,}~{}{\bf n}_{j}={\bf q}_{\rm min},~{}j=1,2,\ldots,k$

$\mbox{Initialize velocities,}~{}{\bf V}=\\{{\bf V}_{I}\\}$

$\mbox{Estimate preconditioner,}~{}{\bf K}_{0}={\bf J}^{-1}$

$\mbox{Initial}~{}\Delta{\bf n}=({\bf K}_{0}{\bf J})^{-1}{\bf K}_{0}({\bf
q}_{\rm min}[{\bf n}]-{\bf n})={\bf 0}$

$t=t_{0}$

while $t<t_{\rm max}$ do

${\bf V}_{I}={\bf V}_{I}+(\delta t/2){\bf F}_{I}/M_{I}$

${\bf n}_{0}=2{\bf n}_{1}-{\bf n}_{2}-\delta t^{2}\omega^{2}\Delta{\bf
n}+\alpha\sum_{l=0}^{k}c_{l}{\bf n}_{1-l}$

${\bf n}_{k}={\bf n}_{k-1},~{}\ldots,~{}{\bf n}_{2}={\bf n}_{1},~{}{\bf
n}_{1}={\bf n}_{0},~{}{\bf n}={\bf n}_{0}$

${\bf R}_{I}={\bf R}_{I}+\delta t{\bf V}_{I}$

${\bf H}[{\bf n}]={\bf H}[{\bf R},{\bf n}],~{}{\bf S}={\bf S}[{\bf R}],~{}{\bf
Z}={\bf S}^{-1/2}$

${\bf q}_{\rm min}[{\bf n}]\Leftarrow\mbox{from diagonalized}~{}{\bf
Z}^{T}{\bf H}[{\bf n}]{\bf Z}$

$\Delta{\bf n}=({\bf K}_{0}{\bf J})^{-1}{\bf K}_{0}({\bf q}_{\rm min}[{\bf
n}]-{\bf n})~{}\mbox{with~{}rank-}m~{}\mbox{approx.}$

${\bf n}^{(1)}={\bf n}-\Delta{\bf n},~{}\mbox{approximate Newton step}$

${\bf q}_{\rm min}[{\bf n}^{(1)}]\Leftarrow\mbox{from diagonalized}~{}{\bf
Z}^{T}{\bf H}[{\bf n}^{(1)}]{\bf Z}$

${\bf q}_{\rm min}[{\bf n}^{(1)}]\Rightarrow\mbox{shadow}~{}{\cal
U}^{(1)}({\bf R},{\bf n})~{}\mbox{and forces},~{}{\bf F}$

${\bf V}_{I}={\bf V}_{I}+(\delta t/2){\bf F}_{I}/M_{I}$

$t=t+\delta t$

end while

##  V Examples

We will demonstrate the accuracy and performance of the shadow energy
functionals and Born-Oppenheimer potentials in XL-BOMD simulations using SCC-
DFTB theory Harrison (1980); Foulkes and Haydock (1989); Porezag  _et al._
(1995); Elstner  _et al._ (1998); Finnis  _et al._ (1998); Frauenheim  _et
al._ (2000); Koskinen and Mäkinen (2009); Gaus  _et al._ (2011); Aradi  _et
al._ (2015); et al. (2020). SCC-DFTB theory can be seen as a framework for
different levels of approximations of density functional theory. Here we will
use the scheme given by a second-order expansion in the charge density
fluctuations around a reference density of overlapping neutral atomic charge
distributions, where the atomic net Mulliken partial charges are used to
describe the long-range electrostatic interactions. In this way the continuous
charge density, $\rho({\bf r})$, of regular DFT becomes vectorized with one
net partial charge per atom, ${\bf q}=\\{q_{I}\\}$. The fluctuating partial
charges are optimized self-consistently to account for interatomic charge
transfer and the response to the long-range electrostatic interactions. In a
general SCC-DFTB scheme this requires a repeated set of constructions of an
approximate effective single-particle Kohn-Sham Hamiltonian, diagonalizations,
charge calculations from the eigenfunctions, and Coulomb potential summations,
until a self-consistent charge convergence is reached. SCC-DFTB theory
therefore follows the same iterative SCF procedure as a regular first-
principles Kohn-Sham DFT calculation. Here we will also use a thermal DFTB
theory, where we assume fractional occupation numbers of the molecular
orbitals determined by the Fermi function at some given electronic
temperature, $T_{e}$, including an electronic entropy term Mermin (1963,
1965); Parr and Yang (1989); Engel and Dreizler (2011); Pittalis  _et al._
(2011); Pribram-Jones  _et al._ (2014); Niklasson (2021a). The fractional
occupation numbers are important to better stabilize the electronic structure
calculations when the electronic HOMO-LUMO energy gap is small or vanishing.
This also affects how we perform the response calculations of $\\{{\bf
f}_{{\bf v}_{i}}\\}$ in the Krylov subspace approximation in Eq. (54) of the
preconditioned kernel Niklasson  _et al._ (2015); Nishimoto (2017); Niklasson
(2017, 2021a); Das and Gavini (2022); Negre  _et al._ (2022).

For our implementation and XL-BOMD simulations we use a developers version of
the LATTE software package Cawkwell and et al. (2010); Cawkwell and Niklasson
(2012); Krishnapryian  _et al._ (2017) that closely follows Alg. 1. As
preconditioner we use an exact calculation of the kernel in the first time
step, and we use a sufficient number of low-rank updates to achieve an
approximate quadratic convergence in the Newton updates. The maximum number,
$m$, of Krylov subspace vectors, i.e. in the rank-$m$ approximation, never
exceeds 6.

First we will look at the asserted scaling expressed in Eqs. (31) and (32)
that were assumed in the derivation of the equations of motion, in Eq. (40)
and Eq. (41). Thereafter, we will demonstrate the advantage of the $1$st-level
update of the shadow energy functional and Born-Oppenheimer potential, ${\cal
U}^{(1)}$, compared to the original $0$th-level approximation for XL-BOMD
simulations of an unstable, charge-sensitive, chemical system.

Figure 1:  Scaling of the residual error terms as a function of time step,
$\delta t$, or harmonic oscillator frequency, $\omega$, for a system of
amorphous carbon with 55 atoms using periodic boundary conditions. The
simulations are performed with a constant dimensionless constant,
$\kappa=\delta t^{2}\omega^{2}$, which means that $\delta
t\propto\omega^{-1}$. XL-BOMD based on the enhanced $1$st-level shadow Born-
Oppenheimer potential, ${\cal U}^{(1)}({\bf R},n^{(0)})$, was used following
Alg. 1. The root-mean square errors (RMSE) are given by the root-mean-square
of the residuals, ${\bf q}_{\rm min}[{\bf n}^{(0)}]-{\bf n}^{(0)}$ and $q_{\rm
min}[{\bf n}^{(1)}]-{\bf n}^{(1)}$, and are averaged over snapshots of 100
integration time steps. The dashed lines indicates the exact $\delta
t^{2}\sim\omega^{-2}$ and $\delta t^{4}\sim\omega^{-4}$ scalings.

###  V.1 Scaling

Figure 1 shows the approximate scaling of the root mean square errors (RMSE)
given by the root mean square of the residuals, ${\bf q}_{\rm min}[{\bf
n}^{(0)}]-{\bf n}^{(0)}$ and ${\bf q}_{\rm min}[{\bf n}^{(1)}]-{\bf n}^{(1)}$
for simulations of amorphous carbon. The results of the simulations confirm
the assumed scaling orders, where $|{\bf q}_{\rm min}[{\bf n}^{(0)}]-{\bf
n}^{(0)}|\propto\omega^{-2}$ and $|{\bf q}_{\rm min}[{\bf n}^{(1)}]-{\bf
n}^{(1)}|\propto\omega^{-4}$. These scalings were asserted a priori in the
derivation of the equations of motion in an adiabatic limit as
$\omega\rightarrow\infty$. The scaling of the RMSE extracted from the XL-BOMD
simulations shown in Fig. 1 confirms these assumption. Notice that the
$\omega^{-1}\propto\delta t$, because our integration scheme, Eq. (48), has
been chosen such that $\delta t^{2}\omega^{2}$ is a dimensionless constant,
${\kappa}=\delta t^{2}\omega^{2}$.

Figure 2:  Scaling of the fractional error in the interatomic forces for water
and amorphous carbon as a function of the integation time step $\delta t$, or
harmonic oscillator frequency, $\omega$. XL-BOMD based on the enhanced $1$st-
level shadow potential, ${\cal U}^{(1)}({\bf R},n^{(0)})$, was used. The
simulations are performed with a constant dimensionless constant,
$\kappa=\delta t^{2}\omega^{2}$, which means that $\delta
t\propto\omega^{-1}$. The fractional error was estimated from an on-the-fly
comparison with the “exact” fully converged Born-Oppenheimer forces, where the
error was averaged over all the atoms and force components over a snapshot of
100 integration time steps. The dashed lines indicates the exact $\sim\delta
t^{4}$ or $\sim\omega^{-4}$ scalings.

The error in the $0$th-level shadow Born-Oppenheimer potential scales as
$|{\cal U}^{(0)}-U|\propto|{\bf q}_{\rm min}[{\bf n}^{(0)}]-{\bf
n}^{(0)}|^{2}$. This means that the error in the sampling of the $0$th-level
shadow Born-Oppenheimer potential, ${\cal U}^{(0})$, scales as $\delta t^{4}$
with the integration time step, which has been confirmed previously, e.g.
Niklasson (2017). The new $1$st-level updated shadow Born-Oppenheimer
potential, ${\cal U}^{(1)}$, has the same form for the error, where $|{\cal
U}^{(1)}-U|\propto|{\bf q}_{\rm min}[{\bf n}^{(1)}]-{\bf n}^{(1)}|^{2}$. This
means, from the scaling demonstrated in Fig. 1, that the error in the sampling
of the shadow Born-Oppenheimer potential ${\cal U}^{(1)}$ scales at $\delta
t^{8}$. It is hard to demonstrate this scaling of the error in ${\cal
U}^{(1)}$ directly, because the error converges to quickly and saturates at a
level set by the available numerical precision. Here we therefore only show
this $\delta t^{8}$-scaling indirectly, from the $\delta t^{4}$ or
$\omega^{-4}$-scaling of $|{\bf q}_{\rm min}[{\bf n}^{(1)}]-{\bf n}^{(1)}|$ in
Fig. 1.

From the derivation of the equations of motion with the $1$st-level shadow
potential in Eqs. (40) and (41), we made the estimate that the equations of
motion for the atomic positions should have an error that scales as
$\propto\omega^{-4}$. In Fig. 2 we show the results of simulations of an
amorphous Carbon and a water system, were we find that the fractional error in
the evaluated forces for the $1$st-level ${\cal U}^{(1)}$ shadow potential
scale at $\propto\delta t^{4}$. This confirms the previously estimated
scaling. This is in contrast to the original $0$th-level shadow Hamiltonian
formulation of XL-BOMD using ${\cal U}^{(0)}$ with an error in the forces that
is only of second order, $\propto\delta t^{2}$ Niklasson (2017).

The dramatic improvement in the scaling of the error as a function of the
integration time step may seem impressive. Nevertheless, often the improved
behavior only has a minor effect on the accuracy and stability of XL-BOMD
simulations. It is only for highly unstable systems, where the improved
scaling and accuracy from the $1th$-level update of the shadow energy
functional and Born-Oppenheimer potential play a role. For such problems we
find that stable molecular trajectories often can be achieved with a slightly
longer integration time step than what otherwise would be possible with the
original $0$th-level shadow energy functional and Born-Oppenheimer potential.

Another important observation is that the higher-degree of accuracy in the
force evaluations may be useful if higher-order symplectic integration schemes
are used. In previous studies, using earlier versions of XL-BOMD, we found
that we needed a fairly tight SCF convergence prior to the force evaluations
for the higher-order symplectic integration schemes in order to take full
advantage of their improved accuracy Odell  _et al._ (2009, 2011). The $1$st-
level shadow energy functional and Born-Oppenheimer potential should therefore
be well-suited in combination with various 4th-order symplectic integrations
schemes Odell  _et al._ (2009).

Figure 3:  XL-BOMD simulations based on SCC-DFTB theory of an artificial
highly reactive randomized mixture of liquid nitromethane (49 atoms with
periodic boundary conditions). The upper panel a) shows the statistical
temperature, the middle panel b) shows the fluctuations in the total energy
per atom, and the lower panel c) shows the HOMO-LUMO electronic energy gap. An
integration time step of $\delta t=0.2$ fs was used in combination with a
fractional occupation number corresponding to an electronic temperature,
$T_{e}=1,500$ K. The $1$st-level updated shadow potential, ${\cal U}^{(1)}$,
(blue lines) shows a more stable dynamics without the more pronounced
fluctuations in the total energy fluctuations of the $0$th-level shadow
potential, ${\cal U}^{(0)}$, (red lines).

###  V.2 Unstable mixture of nitromethane

To demonstrate the advantage of the $1$st-level shadow energy functional and
Born-Oppenheimer potential compared to the original $0$th-level approach, we
will look at a chemically unstable system, with a small or vanishing HOMO-LUMO
energy gap. Such systems are often difficult to study, in particular with
regular direct quantum-mechanical Born-Oppenheimer molecular dynamics methods.
As an example we have chosen an artificial mixture of liquid nitromethane,
(CH3NO2)7, where a handful randomly chosen atoms have switched places. This
artificial testbed system is highly unstable and exothermic reactions occurs
within a few hundred femtoseconds. This is illustrated in Fig. 3. We find a
significantly improved stability in the simulation with the $1$st-level
updated shadow potential, ${\cal U}^{(1)}$, (blue solid lines) compared to the
original $0$th-level shadow potential, ${\cal U}^{(0)}$, (red dashed lines) as
indicated by the fluctuations in the total energy shown in the mid panel b).
Only by reducing the integration time step, $\delta t$, or possibly by
increasing the electronic temperature, is it possible to stabilize the XL-BOMD
simulation using the original $0$th-level shadow potential.

##  VI Summary and discussion

In this article we have introduced a generalization of the shadow energy
functionals and Born-Oppenheimer potentials used in XL-BOMD. The original 0th-
level shadow energy functional generates a Born-Oppenheimer potential that has
an error in the fourth-order, ${\cal O}(\delta t^{4})$, of the integration
time step, $\delta t$, and with an error in the interatomic forces that is of
second-order, ${\cal O}(\delta t^{2})$. With the 1st-level update the error in
the potential energy can be reduced to scale as ${\cal O}(\delta t^{8})$,
where the error in the calculated interatomic forces scales as ${\cal
O}(\delta t^{4})$. The main additional cost using the $1$st-level instead of
the $0$th-level shadow potential is the cost of an extra Hamiltonian
diagonalization. We showed how this improved level of accuracy helps stabilize
the integration of the molecular trajectories, which can be of particular
importance for unstable, charge-sensitive, reactive systems with a small or
vanishing electronic HOMO-LUMO energy gap. The improved scaling in the error
of the potential and forces may also be of interest in the application of
higher-order symplectic integration schemes Niklasson (2008); Odell  _et al._
(2009, 2011). These higher-order schemes are of no use unless they can be
matched by force evaluations with a comparable or higher level of accuracy.

The ability to systematically improve the accuracy of the Born-Oppenheimer
potential has many similarities with earlier versions of XL-BOMD Niklasson
_et al._ (2006); Niklasson (2008); Steneteg  _et al._ (2010), where often a
few SCF steps were needed prior to each force evaluation. However, with the
detailed analysis supported by the concept of a shadow dynamics or a backward
error analysis, we now have a more transparent description of why and when
this is the case and how we can optimize the efficiency of our XL-BOMD
simulations. The key idea is the construction of pairs of shadow energy
functionals and potentials, where the shadow potential is given from an exact,
yet computationally cheap, ground-state optimization of a linearized shadow
energy functional. In combination with XL-BOMD, where the electronic degrees
of freedom is propagated dynamically, the shadow Born-Oppenheimer potentials
can then be used to calculate conservative interatomic forces that generates
accurate molecular trajectories with long-term energy stability.

The generalized shadow energy functionals and Born-Oppenheimer potentials were
demonstrated using Kohn-Sham based SCC-DFTB theory. However, the underlying
theory was derived in a general form that also applies to other electronic
structure theories, inlcuding Hartee-Fock and orbital-free DFT. As an example
we also discussed an extension to flexible charge equilibration models, which
can be derived as coarse-grained versions of Hohenberg-Kohn DFT. The higher-
level generalization of the shadow energy functionals and Born-Oppenheimer
potentials presented here are therefore applicable to a broad variety of
electronic structure methods and flexible charge models within the framework
of XL-BOMD.

##  VII Acknowledgements

This work is supported by the U.S. Department of Energy Office of Basic Energy
Sciences (FWP LANLE8AN,“Next generation quantum-based molecular dynamic”) and
by the U.S. Department of Energy through the Los Alamos National Laboratory.
Discussions with Joshua Finkelstein are gratefully acknowledged. Los Alamos
National Laboratory is operated by Triad National Security, LLC, for the
National Nuclear Security Administration of the U.S. Department of Energy
Contract No. 892333218NCA000001.

## References

  * Yoshida (1990) H. Yoshida, Phys. Lett. A 150, 262 (1990).
  * Grebogi  _et al._ (1990) C. Grebogi, S. M. Hammel, J. A. Yorke, and T. Saur, Phys. Rev. Lett. 65, 1527 (1990).
  * Toxvaerd (1994) S. Toxvaerd, Phys. Rev. E 50, 2271 (1994).
  * Gans and Shalloway (2000) J. Gans and D. Shalloway, Phys. Rev. E 61, 4587 (2000).
  * Bond and Leimkuhler (2007) S. D. Bond and B. J. Leimkuhler,  _Molecular dynamics and the accuracy of numerically computed averages_ (Cambride University Press, United Kingdom, 2007).
  * Toxvaerd  _et al._ (2012) S. Toxvaerd, O. J. Heilmann, and J. C. Dyre, J. Chem. Phys. 136, 224106 (2012).
  * Hammonds and Heyes (2020) K. D. Hammonds and D. M. Heyes, J. Chem. Phys. 152, 024114 (2020).
  * Hammonds and Heyes (2021) K. D. Hammonds and D. M. Heyes, J. Chem. Phys. 154, 174102 (2021).
  * Müser (2022) M. H. Müser, Molecular Simulation 48, 1393 (2022), https://doi.org/10.1080/08927022.2022.2094430 .
  * Niklasson  _et al._ (2007) A. M. N. Niklasson, C. J. Tymczak, and M. Challacombe, J. Chem. Phys. 126, 144103 (2007).
  * Niklasson (2008) A. M. N. Niklasson, Phys. Rev. Lett. 100, 123004 (2008).
  * Cawkwell and Niklasson (2012) M. J. Cawkwell and A. M. N. Niklasson, J. Chem. Phys. 137, 134105 (2012).
  * Hutter (2012) J. Hutter, WIREs Comput. Mol. Sci. 2, 604 (2012).
  * Lin  _et al._ (2014) L. Lin, J. Lu, and S. Shao, Entropy 16, 110 (2014).
  * Souvatzis and Niklasson (2014) P. Souvatzis and A. M. N. Niklasson, J. Chem. Phys. 140, 044117 (2014).
  * Niklasson (2017) A. M. N. Niklasson, J. Chem. Phys. 147, 054103 (2017).
  * Niklasson (2021a) A. M. N. Niklasson, Eur. Phys. J. B 94, 164 (2021a).
  * Patrick Henning (2021) A. M. N. N. Patrick Henning, Kinetic and Related Models 14, 303 (2021).
  * Niklasson (2021b) A. M. N. Niklasson, J. Chem. Phys. 154, 0000 (2021b).
  * Marx and Hutter (2000) D. Marx and J. Hutter, “Modern methods and algorithms of quantum chemistry,” (ed. J. Grotendorst, John von Neumann Institute for Computing, Jülich, Germany, 2000) 2nd ed.
  * Tuckerman (2010) M. E. Tuckerman,  _Statistical Mechanics: Theory and Molecular Simulation_ (Oxford University Press, New York, 2010).
  * Roothaan (1951) C. C. J. Roothaan, Rev. Mod. Phys. 23, 69 (1951).
  * McWeeny (1959) R. McWeeny, Phys. Rev. 114, 1528 (1959).
  * Mermin (1963) N. D. Mermin, Annals of Physics 21, 99 (1963).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B:864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, 1133 (1965).
  * Mermin (1965) N. D. Mermin, Phys. Rev. B 137, A1441 (1965).
  * Parr and Yang (1989) R. G. Parr and W. Yang,  _Density-functional theory of atoms and molecules_ (Oxford University Press, Oxford, 1989).
  * Dreizler and Gross (1990) R. Dreizler and K. Gross,  _Density-functional theory_ (Springer Verlag, Berlin Heidelberg, 1990).
  * Engel and Dreizler (2011) E. Engel and R. Dreizler,  _Density-functional theory_ (Springer Verlag, Berlin Heidelberg, 2011).
  * Remler and Madden (1990) D. K. Remler and P. A. Madden, Mol. Phys. 70, 921 (1990).
  * Pulay and Fogarasi (2004) P. Pulay and G. Fogarasi, Chem. Phys. Lett. 386, 272 (2004).
  * Herbert and Head-Gordon (2005) J. Herbert and M. Head-Gordon, Phys. Chem. Chem. Phys. 7, 3269 (2005).
  * Niklasson  _et al._ (2006) A. M. N. Niklasson, C. J. Tymczak, and M. Challacombe, Phys. Rev. Lett. 97, 123001 (2006).
  * Kühne  _et al._ (2007) T. D. Kühne, M. Krack, F. R. Mohamed, and M. Parrinello, Phys. Rev. Lett. 98, 066401 (2007).
  * Car and Parrinello (1985) R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
  * Elstner  _et al._ (1998) M. Elstner, D. Poresag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert, Phys. Rev. B 58, 7260 (1998).
  * Finnis  _et al._ (1998) M. W. Finnis, A. T. Paxton, M. Methfessel, and M. van Schilfgarde, Phys. Rev. Lett. 81, 5149 (1998).
  * et al. (2020) B. H. et al., J. Chem. Phys. 152, 124101 (2020).
  * Dewar and Thiel (1977) M. J. S. Dewar and W. Thiel, Theoret.Chim. Acta 46, 89 (1977).
  * Dewar  _et al._ (1985) M. J. S. Dewar, E. G. Zoebisch, E. F. Healy, and J. J. P. Stewart, J. Am. Chem. Soc. 107, 3902 (1985).
  * Stewart (2013) J. J. P. Stewart, J. Mol. Model. 19, 1 (2013).
  * Bannwarth  _et al._ (2018) C. Bannwarth, S. Ehlert, and S. Grimme, J. Chem. Theory Comput. 15, 1652 (2018).
  * Dral  _et al._ (2019) P. O. Dral, X. Wu, and W. Thiel, J. Chem. Theory Comput. 15, 1743 (2019).
  * Malone  _et al._ (2020) W. Malone, B. Nebgen, A. White, Y. Zhang, H. Song, J. A. Bjorgaard, A. E. Sifain, B. Rodriguez-Hernandez, V. M. Freixas, S. Fernandez-Alberti, A. E. Roitberg, T. R. Nelson, and S. Tretiak, Journal of Chemical Theory and Computation 16, 5771 (2020), pMID: 32635739, https://doi.org/10.1021/acs.jctc.0c00248 .
  * Zhou  _et al._ (2020) G. Zhou, B. Nebgen, N. Lubbers, W. Malone, A. M. N. Niklasson, and S. Tretiak, Journal of Chemical Theory and Computation 16, 4951 (2020), pMID: 32609513, https://doi.org/10.1021/acs.jctc.0c00243 .
  * Bannwarth  _et al._ (2020) C. Bannwarth, E. Caldeweyher, S. Ehlert, A. H. ans P. Pracht, J. Seibert, S. Spicher, and S. Grimme, WIREs Comput. Lol. Sci. 11, 1 (2020).
  * Vesely (1977) F. J. Vesely, J. Comput. Phys. 24, 361 (1977).
  * Sprik and Klein (1988) M. Sprik and M. L. Klein, The Journal of Chemical Physics 89, 7556 (1988), https://doi.org/10.1063/1.455722 .
  * Mortier  _et al._ (1986) W. J. Mortier, S. K. Ghosh, and S. Shankar, Journal of the American Chemical Society 108, 4315 (1986), https://doi.org/10.1021/ja00275a013 .
  * Rappe and III (1991) A. K. Rappe and W. A. G. III, J. Phys. Chem 95, 3358 (1991).
  * Lamoureux and Roux (2003) G. Lamoureux and B. T. Roux, J. Chem. Phys. 119, 3025 (2003).
  * Verstraelen  _et al._ (2013) T. Verstraelen, P. W. Ayers, V. Van Speybroeck, and M. Waroquier, The Journal of Chemical Physics 138, 074108 (2013), https://doi.org/10.1063/1.4791569 .
  * Naserifar  _et al._ (2017) S. Naserifar, D. J. Brooks, W. A. Goddard, and V. Cvicek, The Journal of Chemical Physics 146, 124117 (2017), https://doi.org/10.1063/1.4978891 .
  * York and Yang (1996) D. M. York and W. Yang, The Journal of Chemical Physics 104, 159 (1996), https://doi.org/10.1063/1.470886 .
  * Tabacchi  _et al._ (2002) G. Tabacchi, C. J. Mundy, J. Hutter, and M. Parrinello, J. Chem. Phys. 117, 1416 (2002).
  * Harris (1985) J. Harris, Phys. Rev. B 31, 1770 (1985).
  * Foulkes and Haydock (1989) W. M. C. Foulkes and R. Haydock, Phys. Rev. B 39, 12520 (1989).
  * Harrison (1980) W. A. Harrison,  _Electronic structure and the properties of solids: the physics of the chemical bond_ (Dover, New York, 1980).
  * Porezag  _et al._ (1995) D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, and R. Kaschner, Phys. Rev. B 51, 12947 (1995).
  * Frauenheim  _et al._ (2000) T. Frauenheim, G. Seifert, M. Elstner, Z. Hajnal, G. Jungnickel, D. Poresag, S. Suhai, and R. Scholz, Phys. Stat. sol. 217, 41 (2000).
  * Koskinen and Mäkinen (2009) P. Koskinen and V. Mäkinen, Computational Materials Science 47, 237 (2009).
  * Gaus  _et al._ (2011) M. Gaus, Q. Cui, and M. Elstner, J, Chem. Theory Comput. 7, 931 (2011).
  * Aradi  _et al._ (2015) B. Aradi, A. M. N. Niklasson, and T. Frauenheim, J. Chem. Theory Comput. 11, 3357 (2015).
  * Note (1) All physically relevant electron densities determined by anti-symmetric electron wavefunctions Dreizler and Gross (1990); Engel and Dreizler (2011).
  * Pittalis  _et al._ (2011) S. Pittalis, C. R. Proetto, A. Floris, A. Sanna, C. Bersier, K. Burke, and E. K. U. Gross, Phys. Rev. Lett. 107, 163001 (2011).
  * Pribram-Jones  _et al._ (2014) A. Pribram-Jones, S. Pittalis, E. K. U. Gross, and K. Burke, in  _Frontiers and Challenges in Warm Dense Matter_ , edited by F. Graziani, M. P. Desjarlais, R. Redmer, and S. B. Trickey (Springer International Publishing, Cham, 2014) pp. 25–60.
  * Heitler and London (1927) W. Heitler and F. London, Z. Phys. 44, 455 (1927).
  * Born and Oppenheimer (1927) M. Born and R. Oppenheimer, Ann. Phys. 389, 475 (1927).
  * Tuckerman (2002) M. E. Tuckerman, J. Phys.: Conden. Matter 14, 1297 (2002).
  * Note (2) Notice that exact ground state, $\rho_{\rm min}({\bf r})$, and its approximate ground state solution, $\rho\approx\rho_{\rm min}$, depend on ${\bf R}$, i.e. $\rho_{\rm min}({\bf r})\equiv\rho_{\rm min}({\bf R},{\bf r})$, but we have dropped the explicit ${\bf R}$-dependencies in our simplified notation.
  * Note (3) If $n^{(0)}=\rho_{\rm min}$ then $\rho_{\rm min}[n^{(0)}]=n^{(0)}$ and ${\cal U}^{(0)}=U$.
  * Niklasson (2020a) A. M. N. Niklasson, J. Chem. Phys. 152, 104103 (2020a).
  * Niklasson (2020b) A. M. N. Niklasson, J. Chem. Theory Comput. 16, 3628 (2020b).
  * Das and Gavini (2022) S. Das and V. Gavini, “Accelerating self-consistent field iterations in kohn-sham density functional theory using a low rank approximation of the dielectric matrix,”  (2022).
  * Negre  _et al._ (2022) C. F. A. Negre, M. E. Wall, and A. M. N. Niklasson, “Graph-based quantum response theory and shadow born-oppenheimer molecular dynamics,”  (2022).
  * Steneteg  _et al._ (2010) P. Steneteg, I. A. Abrikosov, V. Weber, and A. M. N. Niklasson, Phys. Rev. B 82, 075110 (2010).
  * Arita  _et al._ (2014) M. Arita, D. R. Bowler, and T. Miyazaki, J. Chem. Theory Comput. 10, 5419 (2014).
  * Peters  _et al._ (2017) L. D. M. Peters, J. Kussmann, and C. Ochsenfeld, J. Chem. Theory Comput. 13, 5479 (2017).
  * Niklasson  _et al._ (2009) A. M. N. Niklasson, P. Steneteg, A. Odell, N. Bock, M. Challacombe, C. J. Tymczak, E. Holmstrom, G. Zheng, and V. Weber, J. Chem. Phys. 130, 214109 (2009).
  * Zheng  _et al._ (2011) G. Zheng, A. M. N. Niklasson, and M. Karplus, J. Chem. Phys. 135, 044122 (2011).
  * Odell  _et al._ (2009) A. Odell, A. Delin, B. Johansson, N. Bock, M. Challacombe, and A. M. N. Niklasson, J. Chem. Phys. 131, 244106 (2009).
  * Odell  _et al._ (2011) A. Odell, A. Delin, B. Johansson, M. J. Cawkwell, and A. M. N. Niklasson, J. Chem. Phys. 135, 224105 (2011).
  * Leven and Head-Gordon (2019) I. Leven and T. Head-Gordon, Phys. Chem. Chem. Phys. 21, 18652 (2019).
  * An  _et al._ (2020) D. An, S. Y. Cheng, T. Head-Gordon, L. Lin, and J. Lu, “Convergence of stochastic-extended lagrangian molecular dynamics method for polarizable force field simulation,” (2020), arXiv:1904.12082 [math.NA] .
  * Niklasson  _et al._ (2015) A. M. N. Niklasson, M. J. Cawkwell, E. H. Rubensson, and E. Rudberg, Phys. Rev. E 92, 063301 (2015).
  * Nishimoto (2017) Y. Nishimoto, J. Chem. Phys. 146, 084101 (2017).
  * Niklasson and Cawkwell (2014) A. M. N. Niklasson and M. Cawkwell, J. Chem. Phys. 141, 164123 (2014).
  * Feynman (1939) R. P. Feynman, Phys. Rev. 56, 340 (1939).
  * Cawkwell and et al. (2010) M. J. Cawkwell and et al., “LATTE,” (2010), Los Alamos National Laboratory (LA- CC-10004), http://www.github.com/lanl/latte.
  * Krishnapryian  _et al._ (2017) A. Krishnapryian, P. Yang, A. M. N. Niklasson, and M. J. Cawkwell, J. Chem. Theory Comput. 13, 6191 (2017).