# Fast method for quantum mechanical molecular dynamics

Anders M. N. Niklasson111Corresponding Author Email<EMAIL_ADDRESS>Marc J.
Cawkwell  Theoretical Division, Los Alamos National Laboratory, Los Alamos,
New Mexico 87545

###### Abstract

As the processing power available for scientific computing grows, first
principles Born-Oppenheimer molecular dynamics simulations are becoming
increasingly popular for the study of a wide range of problems in materials
science, chemistry and biology. Nevertheless, the computational cost of Born-
Oppenheimer molecular dynamics still remains prohibitively large for many
potential applications. Here we show how to avoid a major computational
bottleneck: the self-consistent-field optimization prior to the force
calculations. The optimization-free quantum mechanical molecular dynamics
method gives trajectories that are almost indistinguishable from an “exact”
microcanonical Born-Oppenheimer molecular dynamics simulation even when low
pre-factor linear scaling sparse matrix algebra is used. Our findings show
that the computational gap between classical and quantum mechanical molecular
dynamics simulations can be significantly reduced.

electronic structure theory, molecular dynamics, Born-Oppenheimer molecular
dynamics, tight-binding theory, self-consistent tight binding theory, self-
consistent-charge density functional tight-binding theory, density matrix,
linear scaling electronic structure theory, Car-Parrinello molecular dynamics,
self-consistent field, extended Lagrangian molecular dynamics

††preprint: LA-UR 12-23991

##  I Introduction

The past three decades have witnessed a dramatic increase in the use of the
molecular dynamics simulation method Karplus and McCammon (2002); Marx and
Hutter (2000). While it is unquestionably a powerful and widely used tool, its
ability to calculate physical properties is limited by the quality and the
computational complexity of the interatomic potentials. Among computationally
tractable models, the most accurate are explicitly quantum mechanical with
interatomic forces calculated on-the-fly using a nuclear potential energy
surface that is determined by the electronic ground state within the Born-
Oppenheimer approximation Wang and Karplus (1973); Leforestier (1978); Marx
and Hutter (2000). In Hartree-Fock Roothaan (1951); McWeeny (1960) or density
functional theory Hohenberg and Kohn (1964); Kohn and Sham (1965); Parr and
Yang (1989); Dreizler and Gross (1990), the electronic ground-state density is
given through a self-consistent-field (SCF) optimization procedure, which
involves iterative mixed solutions of the single-particle eigenvalue equations
and accounts for details in the charge distribution. Since the interatomic
forces are sensitive to the electrostatic potential Feynman (1939), molecular
dynamics simulations are often of poor quality without a high degree of self-
consistent-field convergence. This is unfortunate since the iterative self-
consistent-field procedure is computationally expensive and in practice always
approximate.

Recently there have been efforts to reduce the computational cost of the self-
consistent-field optimization without causing any significant deviation from
“exact” Born-Oppenheimer molecular dynamics simulations Pulay and Fogarasi
(2004); Niklasson et al. (2006); Kühne et al. (2006). In this article we go
one step further, and in analogy to time-dependent techniques such as
Ehrenfest molecular dynamics Ehrenfest (1927); Alonso (2008); Jakowski (2009)
or the Car-Parrinello method Car and Parrinello (1985); Marx and Hutter
(2000); Tuckerman (2002); Hartke (1992); Schlegel (1992); Herbert (1992);
Kirchner (2012), we show how the electronic ground state optimization can be
circumvented fully without any noticeable reduction in accuracy in comparison
to “exact” Born-Oppenheimer molecular dynamics.

Our optimization-free dynamics is based on a reformulation of extended
Lagrangian Born-Oppenheimer molecular dynamics Niklasson (2008) in the limit
of vanishing self-consistent-field optimization. The method is presented
within a general free energy formulation that is valid also at finite
electronic temperatures and should be applicable to a broad class of
materials. In addition to the removal of the costly self-consistent-field
optimization we also demonstrate compatibility with low pre-factor linear
scaling electronic structure theory Goedecker (1999); Bowler and Miyazaki
(2012); Niklasson (2002). The combined scheme provides a very efficient,
energy conserving, low-complexity method for performing accurate quantum
molecular dynamics simulations.

##  II Fast Quantum Mechanical Molecular Dynamics

###  II.1 Born-Oppenheimer molecular dynamics

Born-Oppenheimer molecular dynamics based on density functional theory can be
described by the Lagrangian

${\cal L}^{\rm BO}({\bf
R},{\bf\dot{R}})=\frac{1}{2}\sum_{k}M_{k}{\dot{R}}_{k}^{2}-U[{\bf R};\rho],$
(1)

where the potential energy,

$\begin{array}[]{l}{\displaystyle U[{\bf R};\rho]=2\sum_{i\in\rm
occ}\varepsilon_{i}-\frac{1}{2}\iint\frac{\rho({\bf r})\rho({\bf
r^{\prime}})}{|{\bf r^{\prime}-r}|}d{\bf r^{\prime}}d{\bf r}}\\\
~{}~{}~{}~{}~{}{\displaystyle-\int V_{\rm xc}[\rho]\rho({\bf r})d{\bf
r}+E_{\rm xc}[\rho]+E_{\rm zz}[{\bf R}],}\end{array}$ (2)

is calculated at the self-consistent electronic ground state density,
$\rho({\bf r})$, for the nuclear configuration ${\bf R}=\\{R_{k}\\}$ Parr and
Yang (1989); Dreizler and Gross (1990). Here, $\varepsilon_{i}$ are the
(doubly) occupied eigenvalues of the effective single-particle Kohn-Sham
Hamiltonian,

$H[\rho]=-\frac{1}{2}\nabla^{2}+V_{\rm n}({\bf R},{\bf r})+\int\frac{\rho({\bf
r^{\prime}})}{|{\bf r^{\prime}-r}|}d{\bf r^{\prime}}+V_{\rm xc}[\rho],$ (3)

where $V_{\rm xc}[\rho]$ is the exchange correlation potential, $V_{\rm
n}({\bf R},{\bf r})$ the external (nuclear) potential, and
$-\frac{1}{2}\nabla^{2}$ the kinetic energy operator. $E_{\rm zz}[{\bf R}]$ is
the electrostatic ion-ion repulsion and $E_{\rm xc}[\rho]$ the exchange
correlation energy.

If the electron density deviates from the ground state density $\rho$ by some
small amount $\delta\rho$, the error in the potential energy is essentially of
the order $\delta\rho^{2}$, depending on the particular formulation used for
calculating $U[{\bf R};\rho+\delta\rho]$ Harris (1985); Sutton et al. (1988);
Foulkes and Haydock (1989). However, since the Hellmann-Feynman theorem is
valid only at the ground state density, we do not have a simple expression for
the forces that avoids calculating derivatives of the electronic density,
$\partial(\rho+\delta\rho)/\partial R_{k}$. In practical calculations, the
accuracy of the potential energy can therefore not be expected to hold also
for the forces and a high degree of self-consistent-field convergence is
therefore typically required.

###  II.2 Extended Lagrangian molecular dynamics

Here we outline how we can circumvent the self-consistent-field procedure in
Born-Oppenheimer molecular dynamics. Instead of recalculating the ground state
density before each force evaluation with an iterative optimization procedure,
the idea here is to use an auxiliary density $n({\bf r})$, as in extended
Lagrangian Born-Oppenheimer molecular dynamics Niklasson (2008); Niklasson et
al. (2009); Steneteg et al. (2010); Zheng et al. (2011); Niklasson et al.
(2011), which evolves through a harmonic oscillator centered around the ground
state density $\rho({\bf r})$. Based on a general free energy formulation of
extended Lagrangian Born-Oppenheimer molecular dynamics Niklasson et al.
(2011) in the limit of vanishing self-consistent-field optimization, we define
the extended Lagrangian:

$\begin{array}[]{l}{\displaystyle{\cal L}({\bf
R},{\bf\dot{R}},n,{\dot{n}})=\frac{1}{2}\sum_{k}M_{k}{\dot{R}}_{k}^{2}-{\cal
U}[{\bf R};n]+T_{e}{\cal S}[{\bf R};n]}\\\
{\displaystyle+\frac{1}{2}\mu\int{\dot{n}}({\bf r})^{2}d{\bf
r}-\frac{1}{2}\mu\omega^{2}\int\left(\rho({\bf r})-n({\bf r})\right)^{2}d{\bf
r}}.\end{array}$ (4)

While the potential and entropy terms, ${\cal U}$ and ${\cal S}$, are well
defined at the ground state density Parr and Yang (1989), i.e. when $n=\rho$,
there are several different options when $n$ deviates from $\rho$, e.g. the
Harris-Foulkes functional Harris (1985); Sutton et al. (1988); Foulkes and
Haydock (1989). In a more general case, the potential energy and entropy term
may therefore also be determined by $n({\bf r})$ implicitly through an
additional function $\sigma[n({\bf r})]$, i.e. ${\cal U}[{\bf R};n]\equiv{\cal
U}[{\bf R};n,\sigma[n]]$ and ${\cal S}[{\bf R};n]\equiv{\cal S}[{\bf
R};n,\sigma[n]]$. Here $\sigma[n({\bf r})]$ is a temperature dependent density
given from the diagonal part of the real-space representation of the (doubly
occupied) density matrix, which is given through a Fermi-operator expansion
Parr and Yang (1989) of the effective single-particle Hamiltonian, $H[n]$,
i.e.

$\sigma({\bf r})\equiv\sigma[n({\bf
r})]=2\left.{\left(e^{\beta(H[n]-\mu_{0}I)}+1\right)^{-1}}\right\rvert_{\bf
r=r^{\prime}}.$ (5)

At zero electronic temperature the Fermi-operator expansion corresponds to a
step function with the step formed at the chemical potential, $\mu_{0}$. In
our Lagrangian above, $\mu$ and $\omega$ are fictitious mass and frequency
parameters of the harmonic oscillator and $\beta$ is the inverse electronic
temperature, i. e. $\beta=1/(k_{\rm B}T_{e})$. The purpose of the entropy-like
term ${\cal S}[{\bf R};n]$ is here to make the derived forces of our dynamics
variationally correct for a given entropy-independent density, $n({\bf r})$,
at any electronic temperature. This approach is different from the regular
formulation where the density is determined by the entropy through the
minimization of the electronic free energy functional Parr and Yang (1989);
Weinert and Davenport (1992); Wentzcovitch et al. (1992); Niklasson (2008b).

####  II.2.1 Equations of motion

The molecular trajectories corresponding to the extended free energy
Lagrangian ${\cal L}$ in Eq. (4) are determined by the Euler-Lagrange
equations of motion,

$\begin{array}[]{l}{\displaystyle
M_{k}{\ddot{R}}_{k}=\left.{-\frac{\partial{\cal U}[{\bf R};n]}{\partial
R_{k}}}\right\rvert_{n}+T_{e}\left.{\frac{\partial{\cal S}[{\bf
R};n]}{\partial R_{k}}}\right\rvert_{n}}\\\ {}{}\\\
{\displaystyle-\left.{\frac{\mu\omega^{2}}{2}\frac{\partial}{\partial
R_{k}}\int\left(\rho({\bf r})-n({\bf r})\right)^{2}d{\bf
r}}\right\rvert_{n}},\\\ \end{array}$ (6)

and

$\begin{array}[]{l}{\displaystyle\mu{\ddot{n}}({\bf
r})=\mu\omega^{2}\left(\rho({\bf r})-n({\bf r})\right)}\\\ {}{}\\\
{\displaystyle-\left.{\frac{\delta{\cal U}[{\bf R};n]}{\delta
n}}\right\rvert_{\bf R}+T_{e}\left.{\frac{\delta{\cal S}[{\bf R};n]}{\delta
n}}\right\rvert_{\bf R}},\end{array}$ (7)

where the partial derivatives are taken with respect to constant density, $n$,
or coordinates, ${\bf R}$. The limit $\mu\rightarrow 0$ gives us the equations
of motion of our extended Lagrangian dynamics,

$M_{k}{\ddot{R}}_{k}=-\left.{\frac{\partial{\cal U}[{\bf R};n]}{\partial
R_{k}}}\right\rvert_{n}+T_{e}\left.{\frac{\partial{\cal S}[{\bf
R};n]}{\partial R_{k}}}\right\rvert_{n}$ (8) ${\ddot{n}}({\bf
r})=\omega^{2}\big{(}\rho({\bf r})-n({\bf r})\big{)},$ (9)

where we have defined ${\cal S}[{\bf R};n]$ such that

$\left.{\frac{\delta{\cal U}[{\bf R};n]}{\delta n}}\right\rvert_{\bf
R}=T_{e}\left.{\frac{\delta{\cal S}[{\bf R};n]}{\delta n}}\right\rvert_{\bf
R}.$ (10)

As is shown in the Appendix (Sec. VII.1), the corresponding property for
$\partial{\cal S}/\partial R_{k}$ is also of importance for the calculation of
the Pulay force in Eq. (8). Notice that these equations still require a full
self-consistent field optimization, since the auxiliary density $n({\bf r})$
evolves around the ground state density $\rho({\bf r})$.

Since the nuclear degrees of freedom do not depend on the mass parameter $\mu$
in Eqs. (8) and (9), the total free energy,

$E_{\rm tot}=\frac{1}{2}\sum_{k}M_{k}{\dot{R}}_{k}^{2}+{\cal U}[{\bf
R};n]-T_{e}{\cal S}[{\bf R};n],$ (11)

is a constant of motion in the limit of vanishing $\mu$. Moreover, if $E_{\rm
tot}$ is close to the exact ground state free energy for approximate densities
$n({\bf r})$, we can also expect that the forces of the extended Lagrangian
dynamics should be accurate.

The forces in Eq. (8) are calculated at the approximate, unrelaxed, density
$n({\bf r})$ using a Hellmann-Feynman-like expression, where the partial
derivatives are taken with respect to a constant density $n({\bf r})$. This is
possible only because $n({\bf r})$ appears as an independent dynamical
variable. In general, as mentioned above, this can not be assumed, since the
Hellmann-Feynman force expression is formally applicable only at the ground
density. A more detailed derivation of explicit force expressions is given in
the Appendix.

####  II.2.2 Entropy contribution

Depending on the particular functional form chosen for the potential energy
term, ${\cal U}({\bf R};n)$, we may not have access to a simple explicit
expression of ${\cal S}[{\bf R};n]$ that fulfills Eq. (10). In this case an
approximate entropy term has to be used. This has no effect on the dynamics in
Eqs. (8) and (9), since the forces remain exact by definition. An
approximation of the entropy term therefore only affects the estimate of the
constant of motion, $E_{\rm tot}$, in Eq. (11).

We have found that the regular expression for the electronic entropy Parr and
Yang (1989),

${\cal S}[{\bf R};n]=-2k_{\rm
B}\sum_{i}\left\\{f_{i}\ln(f_{i})+(1-f_{i})\ln(1-f_{i})\right\\},$ (12)

which formally is defined only at the ground state density, i.e. when
$n=\rho$, typically provides a highly accurate approximation also for
approximate densities as will be illustrated in the examples below. Here
$f_{i}$ are the occupation numbers of the states, i.e. the eigenvalues of the
density matrix in Eq. (5). These are determined by the Fermi-Dirac
distribution of the single-particle eigenvalues $\varepsilon_{i}$ of the
Hamiltonian $H[n]$, i.e.

$f_{i}=\left[e^{\beta(\varepsilon_{i}-\mu)}+1\right]^{-1}.$ (13)

By comparing the calculation of $E_{\rm tot}$ in Eq. (11) using the
approximate entropy term, ${\cal S}[{\bf R};n]$, in Eq. (12) to “exact”, fully
optimized, Born-Oppenheimer molecular simulations, we can estimate the
accuracy of our dynamics.

###  II.3 Fast quantum mechanical molecular dynamics

As in extended Lagrangian Born-Oppenheimer molecular dynamics, the
irreversibility of regular Born-Oppenheimer molecular dynamics that is caused
by the self-consistent-field optimization, can be avoided, since the density
$n({\bf r})$ can be integrated using a reversible geometric integration
algorithm Leimkuhler and Reich (2004); Niklasson (2008); Niklasson et al.
(2009); Odell et al. (2009), e.g. the Verlet algorithm as in Eq. (16) below.
This prevents the unphysical drift in the energy and phase space of regular
Born-Oppenheimer molecular dynamics Pulay and Fogarasi (2004); Niklasson et
al. (2006); Kühne et al. (2006) and our dynamics will therefore exhibit long-
term stability of the free energy $E_{\rm tot}$ in Eq. (11).

A main problem so far is that we still need to calculate the self-consistent
ground state density $\rho({\bf r})$ in the integration of $n({\bf r})$ in Eq.
(9). Fortunately, various geometric integrations of the auxiliary density
$n({\bf r})$ in Eq. (9) are stable also for approximate ground state density
estimates of $\rho({\bf r})$, as long as the approximation of $\rho({\bf r})$
is at least infinitesimally closer to the exact ground state compared to
$n({\bf r})$. Using an integration time step of $\delta t$, this stability
holds if the value of the dimensionless variable $\kappa=\delta
t^{2}\omega^{2}$ is chosen to be appropriately small Niklasson et al. (2009);
Odell et al. (2009, 2011). For energy functionals that are convex in the
vicinity of the ground state density we may therefore replace $\rho({\bf r})$
in Eq. (9) by a linear combination $(1-c)n+c\sigma$ Dederichs and Zeller
(1983), which gives us the approximate equations of motion

$M_{k}{\ddot{R}}_{k}=-\left.{\frac{\partial{\cal U}[{\bf R};n]}{\partial
R_{k}}}\right\rvert_{n}+T_{e}\left.{\frac{\partial{\cal S}[{\bf
R};n]}{\partial R_{k}}}\right\rvert_{n},$ (14)

and

${\ddot{n}}({\bf r})=\omega^{2}\big{(}\sigma({\bf r})-n({\bf r})\big{)},$ (15)

where the constant $\omega^{2}$ has been rescaled by $c$. The Verlet
integration of Eq. (15), including a weak dissipation to avoid an accumulation
of numerical noise Niklasson et al. (2009); Steneteg et al. (2010),

${\displaystyle n_{t+\delta t}=2n_{t}-n_{t-\delta t}+\delta
t^{2}\omega^{2}\left(\sigma_{t}-n_{t}\right)+\alpha\sum_{k=0}^{K}c_{k}n_{t-k\delta
t},}$ (16)

is therefore stable if a sufficiently small positive value of $\kappa=\delta
t^{2}\omega^{2}$ is chosen Niklasson et al. (2009). Thus, without any self-
consistent-field optimization of $\rho({\bf r})$, the previously optimized
values of $\kappa$ in Ref. Niklasson et al. (2009); Odell et al. (2009, 2011)
should be rescaled by a positive factor $\leq 1$. Certain ill behaved (non-
convex) functionals with self-consistent-field instabilities Dederichs and
Zeller (1983) can not be treated in this framework.

The proposed molecular dynamics as given by Eqs. (14) and (15) is the central
result of this paper. The equations of motion do not involve any ground state
self-consistent-field optimization prior to the force evaluations and only one
single diagonalization or density matrix construction is required in each time
step. The frequency $\omega$ of the electronic density is well separated from
the nuclear vibrational oscillations. Using a value of $\delta
t\omega={\sqrt{\kappa}}\approx 1$ and an integration time step $\delta t$,
which is $\sim 1/15$ of the period of the nuclear motion, the frequencies
differ by a factor of 5. As will be demonstrated in the examples below, the
scheme is also fully compatible with linear scaling electronic structure
theory Goedecker (1999); Bowler and Miyazaki (2012). This compatibility is
crucial in order to simulate large systems. The removal of the costly ground
state optimization, in combination with low-complexity linear scaling solvers,
provide a computationally fast quantum mechanical molecular dynamics (Fast-
QMMD) that can match the fidelity and accuracy of regular Born-Oppenheimer
molecular dynamics.

There are several alternative approaches to derive or motivate the equations
of motion of the fast quantum mechanical molecular dynamics, Eqs. (14) and
(15), and details of the dynamics may vary depending on the choice of the
functional form of ${\cal U}({\bf R};n)$. However, the particular derivation
presented here is the most transparent and general approach that we have found
so far.

The equations of motion are given in terms of the electron density, but they
should be generally applicable to a large class of methods, such as Hartree-
Fock theory, which is analyzed in the Appendix (Sec. VII.1), or plane wave
pseudo-potential methods Steneteg et al. (2010). Here we will demonstrate our
fast quantum mechanical molecular dynamics scheme using self-consistent-charge
density functional tight-binding theory Sutton et al. (1988); Finnis et al.
(1998); Elstner et al. (1998); Finnis (2003), as implemented in the electronic
structure code latte Sanville et al. (2010), either with an orthogonal or a
non-orthogonal representation and both at zero and at finite electronic
temperatures. With this method we can easily reach the time and length scales
necessary to establish long-term energy conservation and linear scaling of the
computational cost. Details of the computational method and our particular
choices of ${\cal U}({\bf R};n)$ are given in the Appendix.

Figure 1: Total energy fluctuations, Eq. (11), using “exact” (4 SCF/step)
Born-Oppenheimer molecular dynamics (BOMD), and the fast quantum mechanical
molecular dynamics, Eqs. (14) and (15), (Fast-QMMD), with ($\tau>0$) or
without ($\tau=0$) thresholding applied in the low pre-factor linear scaling
solver Niklasson (2002). The simulations were performed with the molecular
dynamics program latte using self-consistent-charge density functional based
tight-binding theory in an orthogonal formulation at $T_{e}=0$, i.e. as in
Eqs. (57), (58) and (60).

Figure 2: Panel a) shows the x-plane phase space trajectory of a single carbon
atom (C) based on an “exact” (4 SCF/step) Born-Oppenheimer molecular dynamics
(BOMD, dashed line) and the fast quantum mechanical molecular dynamics (Fast-
QMMD, solid line). Panel b) shows the fluctuations in the net auxiliary charge
$n_{i}(t)$ and ground state charge $q_{i}(t)$ for the same carbon atom
($i$=C). The numerical threshold $\tau$ is applied in the linear scaling
solver Niklasson (2002). The simulations were performed with the program latte
using self-consistent-charge tight-binding theory in an orthogonal formulation
at zero electronic temperature, i.e. as in Eqs. (57), (58) and (60).

##  III Examples

###  III.1 Orthogonal representation

Table 1: Wall clock timings of the fast quantum mechanical molecular dynamics (Fast-QMMD) simulations in comparison to Born-Oppenheimer molecular dynamics (BOMD) (4 SCF/step), without ($\tau=0$) and with ($\tau>0$) a low pre-factor linear scaling solver for the density matrix Niklasson (2002) with threshold tolerance $\tau$. The program (latte in its orthogonal formulation at $T_{e}=0$) was executed on a single core of a 2.66 GHz Quad-Core Intel Xeon processor.  Polyethene chain C100H202 | Efficiency
---|---
BOMD ($\tau=0$) | 7.5 s/step
Fast-QMMD ($\tau=0$) | 1.5 s/step
Fast-QMMD ($\tau=10^{-5}$) | 0.61 s/step
Liquid Methane (CH4)100 | Efficiency
BOMD ($\tau=0$) | 12.5 s/step
Fast-QMMD ($\tau=0$) | 2.5 s/step
Fast-QMMD ($\tau=10^{-5}$) | 0.35 s/step

Figure 1 shows the fluctuations in the total energy (kinetic plus potential)
using the fast quantum mechanical molecular dynamics, Eqs. (14) and (15), as
implemented in Eqs. (57), (58) and (60), and an “exact” Born-Oppenheimer
molecular dynamics Niklasson (2008), for liquid methane (density = 0.422
g/cm3) at room temperature. The calculations were performed with the latte
molecular dynamics program using periodic boundary conditions and an
integration time step of $\delta t=0.5$ fs. Since the molecular system is
chaotic, any infinitesimally small deviation will eventually lead to a
divergence between different simulations. However, even after hundreds of time
steps and over 300 fs of simulation time the total energy curves are virtually
on top of each other as is seen in the inset. The same remarkable agreement is
seen in Fig. 2, which shows the projected phase space of an individual carbon
atom and the fluctuations of its net charge. In this case the C atom was
displaced compared to the simulation in Fig. 1 to further enhance the charge
fluctuations.

###  III.2 Linear scaling

The fast quantum mechanical molecular dynamics scheme is also stable in
combination with approximate linear scaling sparse matrix algebra Goedecker
(1999); Bowler and Miyazaki (2012). Using the recursive second order spectral
projection method for the construction of the density matrix Niklasson (2002)
with a numerical threshold, $\tau=10^{-5}$, below which all elements are set
to zero after each individual projection, we notice excellent accuracy and
stability in Fig. 1 without any systematic drift in the total energy.

Despite their high efficiency and low computational pre-factor compared to
alternative linear scaling electronic structure methods Rubensson and Rudberg
(2011), it has been argued that recursive purification algorithms are non-
variational and therefore incompatible with forces of a conservative system
Bowler and Miyazaki (2012), which is necessary for long-term energy
conservation. As is evident from Figs. 1 and 2, this is not a problem. The
graphs are practically indistinguishable from “exact” Born-Oppenheimer
molecular dynamics, without any signs of a systematic drift in the total
energy. The corresponding linear scaling compatibility with microcanonical
simulations was recently also demonstrated for self-consistent-field-optimized
extended Lagrangian Born-Oppenheimer molecular dynamics Cawkwell and Niklasson
(2012).

The gain in speed using the fast quantum mechanical molecular dynamics scheme
in comparison to Born-Oppenheimer molecular dynamics is illustrated by the
wall-clock timings shown in Tab. 1.

Figure 3: Total energy fluctuations, Eq. (11), using “exact” (4 SCF/step)
Born-Oppenheimer molecular dynamics (BOMD), and the fast quantum mechanical
molecular dynamics, Eqs. (14) and 15), (Fast-QMMD), with ($\tau=10^{-5}$) or
without ($\tau=0$) thresholding applied in the low pre-factor linear scaling
solver Niklasson (2002). The simulations were performed with the molecular
dynamics program latte in the non-orthogonal formulation at $k_{\rm B}T_{e}=0$
eV, i.e. as implemented in Eqs. (50), (51) and (55) with the entropy term
approximated by ${\cal S}=0$.

Figure 4: Total free energy fluctuations, Eq. (11), using “exact” (4 SCF/step)
Born-Oppenheimer molecular dynamics (BOMD), and the fast quantum mechanical
molecular dynamics, Eqs. (14) and (15), (Fast-QMMD). The simulations were
performed with the molecular dynamics program latte using the non-orthogonal
formulation at an electronic temperature of $k_{\rm B}T_{e}=0.5$ eV, i.e. as
implemented in Eqs. (50), (51) and (55) with the entropy term approximated by
Eq. (56).

###  III.3 Non-orthogonal representation

For non-orthogonal representations at finite electronic temperatures, a Pulay
force term and a finite approximate entropy contribution to the total free
energy have to be included. Figures 3 and 4 illustrate the total energy
fluctuations for the fast quantum mechanical molecular dynamics simulations of
a hydrocarbon chain as implemented in latte using Eqs. (50), (51) and (55),
with the approximate entropy term in Eq. (56). The electronic temperature of
the examples in Figure 3 is set to zero, $k_{\rm B}T_{e}=0$ eV, and for the
examples in Fig. 4, $k_{\rm B}T_{e}=2$ eV. In the first time step the initial
nuclear temperature, $T_{\rm init}$, was set to $300$ K using a Gaussian
distribution of the velocities. Despite the approximation of $\rho$ in Eq.
(15) and the approximate estimate of the entropy contribution to the free
energy there is virtually no difference seen between the fast quantum
mechanical and the Born-Oppenheimer molecular dynamics simulations.

As in the orthogonal case, the non-orthogonal formulation of our fast quantum
mechanical molecular dynamics is fully compatible with linear scaling
complexity in the construction of the density matrix at $T_{e}=0$ K. In Fig. 3
the reduced complexity simulation shows no significant deviation from “exact”
Born-Oppenheimer molecular dynamics. At finite electronic temperatures, the
linear scaling construction of the Fermi operator Niklasson (2003, 2008b) has
not yet been implemented.

###  III.4 Long-term stability and conservation of the total energy

To assess the long-term energy conservation and the stability we use a test
system comprised of 16 molecules of isocyanic acid, HNCO, at a density of 1.14
g cm-3. The system was first thermalized to a temperature of 300 K over a
simulation time of 12.5 ps by the rescaling of the nuclear velocities. The
simulations used an integration time step, $\delta t$, of 0.25 ps. The
simulations were performed using self-consistent tight-binding theory Sutton
et al. (1988); Finnis et al. (1998); Elstner et al. (1998); Finnis (2003) with
a non-orthogonal basis as implemented in latte, using Eqs. (50), (51) and (55)
with the entropy term approximated by Eq. (56).

Fast quantum mechanical molecular dynamics and “exact” Born-Oppenheimer
molecular dynamics simulations with 4 self-consistent field cycles per time
step were performed over 250,000 time steps (62.5 ps) with $T_{e}=0$ K and
$k_{\text{B}}T_{e}=0.5$ eV. The latter temperature is small with respect to
the HOMO-LUMO gap of HNCO, which is about 6.0 eV, yet the entropy term, Eq.
(12) or Eq. (56), contributes about 0.19 eV to the total energy owing to the
partial occupation of states in the vicinity of the chemical potential.
Trajectories computed at $T_{e}=0$ K with “exact” Born-Oppenheimer molecular
dynamics and the fast quantum mechanical molecular dynamics method without
($\tau=0$) and with ($\tau=10^{-5}$) linear scaling constructions of the
density matrix are presented in Fig. 5. The standard deviation of the
fluctuations of the total energy about its mean and an estimate of the level
of the systematic drift of the total energies are presented in Table 2. These
data show that the fast quantum mechanical molecular dynamics simulations
yield trajectories that are effectively indistinguishable from the “exact”
Born-Oppenheimer trajectories. Moreover, as was seen above, the fast quantum
mechanical molecular dynamics scheme appears to be fully compatible with
linear scaling construction of the density matrix and the resulting
approximate forces, since this trajectory differs from the “exact” Born-
Oppenheimer molecular dynamics trajectory only by a small-amplitude random-
walk of the total energy about its mean Cawkwell and Niklasson (2012). The
systematic drift in energy is several orders of magnitude smaller than in
previous attempts to combine linear scaling solvers with regular Born-
Oppenheimer molecular dynamics Mauri and Galli (1994); Horsfield et al.
(1996); Tushida (2008); Shimojo et al. (2008).

The trajectories computed with an electronic temperature corresponding to
$k_{\text{B}}T_{e}=0.5$ eV differ qualitatively from those computed with zero
electronic temperature. Figure 6 and Table 2 show that while the “exact” Born-
Oppenheimer trajectory conserves the free energy to an extremely high
tolerance over the duration of the simulation, the total free energy in the
fast quantum mechanical molecular dynamics simulation exhibit random-walk
behaviour about the mean value. Although the fast quantum mechanical molecular
dynamics simulations involve an approximate expression for the entropy, we
find that this alone cannot account for the level of fluctuations observed.
Instead, we have found that the rescaling of the $\kappa$ value in the
integration, Eq. (16), affects this random-walk. By changing the rescaling
factor to 3/4, instead of 1/2 as in all the other examples, the amplitude of
the random walk is significantly reduced. Nevertheless, the fast quantum
mechanical molecular dynamics trajectories at finite electronic temperature
exhibit systematic drifts in the total energy that are negligible and the
fluctuations of the total energy about the mean are of the same order as those
that arise from the application of the approximate linear scaling method at
$T_{e}=0$ K.

Figure 5: Total energy versus time for liquid isocyanic acid with a nuclear temperature of 300 K and $T_{e}=0$ K computed with “exact” Born-Oppenheimer MD and the Fast QMMD method with exact and approximate linear scaling density matrix constructions. The numerical threshold $\tau$ is applied in the linear scaling solver Niklasson (2002) below which all matrix elements are set to zero after each iteration. Figure 6: Total free energy versus time for liquid isocyanic acid with a nuclear temperature of 300 K and $k_{\text{B}}T_{e}=0.5$ eV computed with “exact” Born-Oppenheimer molecular dynamics (BOMD) and the fast quantum mechanical molecular dynamics (Fast-QMMD) method with $\kappa$ rescaled by 3/4 instead of 1/2. Table 2: Standard deviation, $\sigma$, of the total energy about its mean value and the upper bound of the systematic drift of the total energy, $E_{\text{drift}}$, computed from “exact” Born-Oppenheimer molecular dynamics (BOMD) and fast quantum mechanical molecular dynamics (Fast-QMMD) simulations of liquid isocyanic acid. The simulation were performed with the latte program in the non-orthogonal formulation, i.e. as implemented in Eqs. (50), (51) and (55) with the entropy term approximated by Eq. (56). $k_{\text{B}}T_{e}$ |  | $\sigma$ | $E_{\text{drift}}$
---|---|---|---
(eV) |  | ($\mu$eV) | ($\mu$eV/atom/ps)
| Fast-QMMD ($\tau=0$) | 0.315 | $5.10\times 10^{-3}$
0.0 | Fast-QMMD ($\tau=10^{-5}$) | 0.702 | 0.285
| BOMD (4 SCF/step) | 0.358 | $9.94\times 10^{-3}$
| Fast-QMMD $(1/2)\kappa$ | 2.47 | 1.43
0.5 | Fast-QMMD $(3/4)\kappa$ | 0.786 | $7.85\times 10^{-2}$
| BOMD (4 SCF/step) | 0.361 | $8.50\times 10^{-2}$

##  IV Convergence properties

Figure 7: The root mean square deviation (RMSD) between the fast quantum
mechanical molecular dynamics, Eqs. (14)-(15), and “exact” (4 SCF/step) Born-
Oppenheimer molecular dynamics, for the nuclear forces, the net Mulliken
charges and the total energy for a Naphthalene molecule at room temperature.
The simulation were performed with the latte molecular dynamics program using
self-consistent-charge tight-binding theory in an orthogonal formulation at
$T_{e}=0$, i.e. as implemented in Eqs. (57), (58) and (60) with ${\cal S}=0$..

The fast quantum mechanical molecular dynamics scheme, Eqs. (14) and (15), can
also be analyzed in terms of the convergence to “exact” Born-Oppenheimer
molecular dynamics as a function of the finite integration time step $\delta
t$. By comparing the deviation in forces, net Mulliken charges, and the total
energy, between the fast quantum mechanical molecular dynamics scheme and an
“exact” Born-Oppenheimer molecular dynamics as a function of $\delta t$ we can
study the consistency between the two methods. Figure 7 shows the difference
between a fully converged “exact” Born-Oppenheimer molecular dynamics
simulation and the fast quantum mechanical molecular dynamics scheme as
measured by the root mean square deviation over 200 fs of simulation time. We
find that the deviation of the nuclear forces, the charges $\\{q_{i}\\}$, as
well as the total energy difference are of the order $\delta t^{2}$ with a
small pre-factor. This convergence demonstrates a consistency between the fast
quantum mechanical scheme and Born-Oppenheimer molecular dynamics using Verlet
integration, where the optimization-free scheme behaves as a well controlled
and tunable approximation. As in “exact” Born-Oppenheimer molecular dynamics,
the dominating error is governed by the local truncation error arising from
the choice of finite integration time step $\delta t$, which is much larger
than any difference between the fast quantum mechanical molecular dynamics and
Born-Oppenheimer molecular dynamics.

##  V Summary and Conclusions

Based on a free energy formulation of extended Lagrangian Born-Oppenheimer
molecular dynamics in the limit of vanishing self-consistent-field
optimization, we have derived and demonstrated a fast quantum mechanical
molecular dynamics scheme, Eqs. (14) and (15), which with a high precision can
match the accuracy and fidelity of Born-Oppenheimer molecular dynamics. In
addition to the removal of the self-consistent-field optimization we have also
demonstrated compatibility with low pre-factor linear scaling solvers. The
combined scheme provides a very efficient, energy conserving, low-complexity
method to perform accurate quantum molecular dynamics simulations. Our
findings show how the computational gap between classical and quantum
mechanical molecular dynamics simulations can be reduced significantly.

##  VI Acknowledgment

We acknowledge support by the United States Department of Energy Office of
Basic Energy Sciences and the LANL Laboratory Directed Research and
Development Program. Discussions with E. Chisolm, J. Coe, T. Peery, S.
Niklasson, C. Ticknor, C.J. Tymczak, and G. Zheng, as well as stimulating
contributions at the T-Division Ten Bar Java group are gratefully
acknowledged. LANL is operated by Los Alamos National Security, LLC, for the
NNSA of the U.S. DOE under Contract No. DE-AC52-06NA25396.

##  VII Appendix

###  VII.1 Calculating the forces in Hartree-Fock theory

Here we present some details of the fast quantum mechanical molecular
dynamics, Eqs. (14) and (15), using a simple but general Hartree-Fock
formalism, which should be directly applicable to a broad class of hybrid and
semi-empirical electronic structure schemes. Instead of the auxiliary density
variable $n({\bf r})$ we will here use the more general density matrix $P$. In
this formalism the extended free-energy Lagrangian in Eq. (4) is given by

$\begin{array}[]{l}{\displaystyle{\cal L}({\bf
R},{\bf\dot{R}};P,{\dot{P}})=\frac{1}{2}\sum_{k}M_{k}{\dot{R}}_{k}^{2}-{\cal
U}[{\bf R};P]+T_{e}{\cal S}[{\bf R};P]}\\\ {\displaystyle+\frac{1}{2}\mu
Tr[{\dot{P}}^{2}]-\frac{1}{2}\mu\omega^{2}Tr[({\mathcal{D}}_{\rm
gs}-P)^{2}],}\\\ \end{array}$ (17)

with the potential energy chosen as

${\cal U}[{\bf R};P]=2Tr[hD(P)]+Tr\\{D(P)G[D(P)]\\},$ (18)

and ground state (gs) density matrix ${\mathcal{D}}_{\rm gs}$. ${\cal S}[{\bf
R};P]$ is an unspecified electronic entropy term, which will be determined by
the requirement to make the derived forces variationally correct, and $T_{e}$
is the electronic temperature. The Fockian (or the effective single-particle
Hamiltonian) is

$F[P]=h+G[P],$ (19)

with the short-hand notation, $G[P]=2J[P]-K[P]$, where $J[P]$ and $K[P]$ are
the conventional Coulomb and exchange matrices, and $h$ is the matrix of the
one-electron part Roothaan (1951); McWeeny (1960). The temperature dependent
density matrix

$D(P)=Z\left(e^{\beta(F^{\perp}[P]-\mu_{0}I)}+1\right)^{-1}Z^{T},$ (20)

which corresponds to $\sigma[n]$ in Eq. (5), is given as a Fermi function of
the orthogonalized Fockian,

$F^{\perp}[P]=Z^{T}F[P]Z.$ (21)

Here $Z$ and its transpose $Z^{T}$ are the inverse Löwdin or Cholesky-like
factors of the overlap matrix, $S$, determined by the relation

$Z^{T}SZ=I.$ (22)

At zero electronic temperature ($T_{e}=0$ K) the Fermi-operator expansion in
Eq. (20) is given by the Heaviside step function, with the step formed at the
chemical potential $\mu_{0}$, separating the occupied from the unoccupied
states.

The Euler-Lagrange equations of motion of ${\cal L}$ in Eq. (17) are given by

$\begin{array}[]{l}{\displaystyle
M_{k}{\ddot{R}}_{k}=-\left.{\frac{\partial{\cal U}}{\partial
R_{k}}}\right\rvert_{P}+T_{e}\left.{\frac{\partial{\cal S}}{\partial
R_{k}}}\right\rvert_{P}}\\\ {}{}\\\
{\displaystyle-\left.{\frac{1}{2}\mu\omega^{2}\frac{\partial}{\partial
R_{k}}Tr[({\mathcal{D}}_{\rm gs}-P)^{2}]}\right\rvert_{P},}\end{array}$ (23)

and

$\mu{\ddot{P}}=\mu\omega^{2}({\mathcal{D}}_{\rm
gs}-P)-\left.{\frac{\partial{\cal U}}{\partial P}}\right\rvert_{\bf
R}+T_{e}\left.{\frac{\partial{\cal S}}{\partial P}}\right\rvert_{\bf R}.$ (24)

A cumbersome but fairly straightforward derivation (see Ref. Niklasson (2008b)
for a closely related example), using the relation and notation
$Z_{R_{k}}=\partial Z/\partial R_{k}=-(1/2)S^{-1}S_{R_{k}}Z$, and defining the
${\cal S}$ term such that

$T_{e}\left.{\frac{\partial{\cal S}}{\partial P_{ij}}}\right\rvert_{\bf
R}=2Tr\left[F^{\perp}[D]D^{\perp}_{P_{ij}}\right]$ (25)

and

$T_{e}\left.{\frac{\partial{\cal S}}{\partial
R_{k}}}\right\rvert_{P}=2Tr\left[F^{\perp}[D]D^{\perp}_{R_{k}}\right],$ (26)

gives the equations of motion

$\begin{array}[]{l}{\displaystyle
M_{k}{\ddot{R}}_{k}=-2Tr[h_{R_{k}}D]-Tr[DG_{R_{k}}(D)]}\\\ {}{}\\\
{\displaystyle+Tr[(DF[D]S^{-1}+S^{-1}F[D]D)S_{R}]}\\\ {}{}\\\
{\displaystyle-\left.{\frac{1}{2}\mu\omega^{2}\frac{\partial}{\partial
R_{k}}Tr[({\mathcal{D}}_{\rm gs}-P)^{2}]}\right\rvert_{P}},\end{array}$ (27)

and

$\mu{\ddot{P}}=\mu\omega^{2}({\mathcal{D}}_{\rm gs}-P).$ (28)

Notice that because of matrix symmetry $P_{ij}$ is not independent form
$P_{ji}$. The partial derivatives of matrix elements $P_{ij}$ are therefore
both over $P_{ij}$ and $P_{ji}$. In the limit $\mu\rightarrow 0$, we get the
final equations of motion for the fast quantum mechanical molecular dynamics
scheme,

$\begin{array}[]{l}{\displaystyle
M_{k}{\ddot{R}}_{k}=-2Tr\left[h_{R_{k}}D\right]-Tr\left[DG_{R_{k}}(D)\right]}\\\
{}{}\\\
{\displaystyle+Tr\left[(DF[D]S^{-1}+S^{-1}F[D]D)S_{R_{k}}\right]},\end{array}$
(29) ${\ddot{P}}=\omega^{2}\left(D(P)-P\right),$ (30)

where we have included the substitution of ${\mathcal{D}}_{\rm gs}$ with
$D(P)$ in the same way as in Eq. (15), i.e. with $\omega^{2}$ rescaled by a
constant $c\leq 1$. The notation for the partial derivative of the two-
electron term is defined as $G_{R_{k}}(D)=\left.{\partial G(D)/\partial
R_{k}}\right\rvert_{D}$, i.e. under the condition of constant density matrix
$D$.

The last term in Eq. (29), which includes the basis-set dependence $S_{R_{k}}$
is the Pulay force term that here is given in a generalized form that is valid
also for non-idempotent density matrices at finite electronic tememperatures
Niklasson (2008b).

###  VII.2 Approximate Entropy contribution

The ${\cal S}[{\bf R};P]$ term is defined such that the two conditions in Eqs.
(25) and (26) are fulfilled. At the self-consistent ground state density, i.e.
when $P=D={\mathcal{D}}_{\rm gs}$, both these conditions are automatically
satisfied by the corresponding regular ground state (gs) electronic entropy
contribution to the free energy Parr and Yang (1989),

$\begin{array}[]{l}{\displaystyle{\cal S}_{\rm gs}[{\bf R};P]={\cal S}_{\rm
gs}[{\bf R};D^{\perp}(P)]}\\\ {}{}\\\ {\displaystyle=-2k_{\rm
B}Tr[D^{\perp}\ln(D^{\perp})+(I-D^{\perp})\ln(I-D^{\perp})]},\end{array}$ (31)

where the relation between $D^{\perp}$ and $D$ is given by the congruence
transformation

$D=ZD^{\perp}Z^{T}.$ (32)

A related derivation is given in Ref. Niklasson (2008b). Using the approximate
estimate ${\cal S}_{\rm gs}[{\bf R};P]$ in Eq. (31) when $P$ and $D$ deviate
from the ground state gives,

$\begin{array}[]{l}{\displaystyle T_{e}\left.{\frac{\partial{\cal S}}{\partial
P_{ij}}}\right\rvert_{\bf
R}=2Tr\left[F^{\perp}[P]D^{\perp}_{P_{ij}}\right]},{}{}\\\ {\displaystyle
T_{e}\left.{\frac{\partial{\cal S}}{\partial
R_{k}}}\right\rvert_{P}=2Tr\left[F^{\perp}[P]D^{\perp}_{R_{k}}\right]},\\\
\end{array}$ (33)

which only approximately fulfills the conditions in Eqs. (25) and (26). It is
possible to show that the error is linear in $\delta P=D-P$ by a linearization
of $F^{\perp}[D]$ around $P$. Since $D(P)$ and $P$ can be assumed to be close
to the ground state, $\delta P$ is small. From the scaling result illustrated
in Fig. 7 the error should therefore be quadratic in the integration time
step, i.e. $\sim\delta t^{2}$. We may therefore approximate the total free
energy using ${\cal S}_{\rm gs}[{\bf R};P]$, which is zero at $T_{e}=0$ K.
However, for the exact formulation and derivation of the equations of motion,
Eqs. (29) and (30), the entropy contribution, $T_{e}{\cal S}[{\bf R};P]$, is
unknown, both at finite and zero temperatures. As is seen in the equations of
motion, Eqs. (29) and (30), this does not affect the forces or the dynamics,
only the estimate of the constant of motion,

$E_{\rm tot}=\frac{1}{2}\sum_{k}M_{k}{\dot{R}}_{k}^{2}+{\cal U}[{\bf
R};P]-T_{e}{\cal S}[{\bf R};P],$ (34)

is approximated. By comparing the approximate $E_{\rm tot}$ to optimized
“exact” Born-Oppenheimer molecular dynamics simulations, the accuracy of the
dynamics can be estimated.

###  VII.3 Alternative potential energy forms

As an alternative to the potential energy, ${\cal U}({\bf R};P)$, in Eq. (18)
we may chose other functional forms that are equivalent at the ground state,
i.e. when $P=D={\cal D}_{gs}$. By using the Harris-Foulkes-like relation
Harris (1985); Foulkes and Haydock (1989),

${\displaystyle Tr[DG(D)]\approx Tr[(2D-P)G(P)]},$ (35)

which has an error of second order in $\delta P=D-P$, we may, for example,
choose

${\cal U}[{\bf R};P]=2Tr[hD(P)]+Tr\\{[2D(P)-P]G(P)\\},$ (36)

as our potential energy term. In this case, the equations of motion at
$T_{e}=0$ corresponding to Eqs. (29) and (30) become

$\begin{array}[]{l}{\displaystyle
M_{k}{\ddot{R}}_{k}=-2Tr\left[h_{R_{k}}D\right]-Tr\\{[2D-P]G_{R_{k}}(P)\\}}\\\
{}{}\\\
{\displaystyle+Tr\left[(DF[P]S^{-1}+S^{-1}F[P]D)S_{R_{k}}\right]},\end{array}$
(37)

and

${\ddot{P}}=\omega^{2}\left(D-P\right),$ (38)

with the constant of motion

$\begin{array}[]{l}{\displaystyle E_{\rm
tot}=\frac{1}{2}\sum_{k}M_{k}{\dot{R}}_{k}^{2}+2Tr[hD]}\\\ {}{}\\\
{\displaystyle+Tr\\{(2D-P)G(P)\\}-T_{e}{\cal S}({\bf R};P)}.\end{array}$ (39)

The entropy term that makes the nuclear forces variationally correct is here
fulfilled by the expression in Eq. (31). With this choice of potential our
dynamics only requires one Fockian or effective single particle Hamitonian
construction per time step. Unfortunately, the error in the Pulay force has
been found to be large compared to Eq. (29). The dynamics in Eqs. (37) and
(38) should therefore be used only for orthogonal representations, i.e. when
the overlap matrix $S=I$.

###  VII.4 Self-Consistent-Charge Density Functional Tight-Binding Theory

In self-consistent-charge density functional based tight-binding theory Sutton
et al. (1988); Finnis et al. (1998); Elstner et al. (1998); Finnis (2003) the
continuous electronic density, $\sigma({\bf n})$, or the density matrix,
$D(P)$, in Eq. (18) is replaced by the net Mulliken charges ${\bf q}[{\bf
n}]=\\{q_{i}\\}$ for each atom $i$, where ${\bf n}=\\{n_{i}\\}$ are the
dynamical variables corresponding to $P$. The potential energy functional
${\cal U}$ in Eq. (18) is then reduced to

${\cal U}[{\bf R};{\bf n}]=2\sum_{i\in\rm
occ}\varepsilon_{i}-\frac{1}{2}\sum_{i,j}q_{i}({\bf n})q_{j}({\bf
n})\gamma_{ij}+E_{\rm pair}[{\bf R}].$ (40)

Here $\varepsilon_{i}$ are the (doubly) occupied eigenvalues of the charge
dependent effective single-particle Hamiltonian

$\begin{array}[]{l}{\displaystyle H_{i\alpha,j\beta}[{\bf
n}]=h_{i\alpha,j\beta}}\\\ {}{}\\\
{\displaystyle+(1/2)\sum_{k\beta^{\prime}}\left(S_{i\alpha,k\beta^{\prime}}V^{ee}_{k\beta^{\prime},j\beta}+V^{ee}_{i\alpha,k\beta^{\prime}}S_{k\beta^{\prime},j\beta}\right)}\end{array}$
(41)

where

$V^{ee}_{j\beta,k\beta^{\prime}}=\sum_{l}q_{l}({\bf
n})\gamma_{jl}\delta_{jk}\delta_{\beta^{\prime}\beta},$ (42)

$h_{i\alpha,j\beta}$ is a parameterized Slater-Koster tight-binding
Hamiltonian, $S_{i\alpha,j\beta}$ the overlap matrix, $i$ and $j$ are atomic
indices and $\alpha$ and $\beta$ are orbital labels Sanville et al. (2010).
The net Mulliken charges are given by

$q_{i}[{\bf n}]\ =2\sum_{\alpha\in
i}\left(\varrho^{\perp}_{i\alpha,i\alpha}-{\varrho^{0}}^{\perp}_{i\alpha,i\alpha}\right),$
(43)

with the density matrix

$\varrho^{\perp}=\varrho^{\perp}[{\bf n}]=\left(e^{\beta(H^{\perp}[{\bf
n}]-\mu_{0}I)}+1\right)^{-1},$ (44)

using the orthogonalized Hamiltonian

$H^{\perp}[{\bf n}]=Z^{T}H[{\bf n}]Z.$ (45)

Here ${\varrho^{0}}$ is the density matrix of the corresponding separate non-
interacting atoms. The de-orthogonalized density matrix is

$\varrho=\varrho[{\bf n}]=Z\varrho^{\perp}[{\bf n}]Z^{T},$ (46)

and as above, the congruence transformation factors are defined through

$Z^{T}SZ=I,$ (47)

where $S$ is the basis set overlap matrix.

The electron-electron interaction in Eq. (40) is determined by $\gamma_{ij}$,
which decays like $1/R$ at large distances and equals the Hubbard repulsion
for the on-site interaction. $E_{\rm pair}[{\bf R}]$ is a sum of pair
potentials, $\phi(R)$, that provide short-range repulsion. The radial
dependence, $\zeta(R)$, of the Slater-Koster bond integrals, elements of the
overlap matrix, and the $\phi(R)$ are all represented analytically in latte by
the mathematically convenient form,

$\zeta(R)=A_{0}\prod_{i=1}^{4}\exp{(A_{i}R^{i})},$ (48)

where $A_{0}$ to $A_{4}$ are adjustable parameters that are fitted to the
results of quantum chemical calculations on small molecules. To ensure that
the off-diagonal elements of $h$ and $S$ and the $\phi(R)$ in our self-
consistent tight-binding implementation decay smoothly to zero at a specified
distance, $R_{\text{cut}}$, we replace the $\zeta(R)$ by cut-off tails of the
form,

$t(R)=B_{0}+\Delta R(B_{1}+\Delta R(B_{2}\\\ +\Delta R(B_{3}+\Delta
R(B_{4}+\Delta RB_{5}))))$ (49)

at $R=R_{1}$, where $\Delta R=R-R_{1}$ and $B_{0}$ to $B_{5}$ are adjustable
parameters. The adjustable parameters are parameterized to match the value and
first and second derivatives of $t(R)$ and $\zeta(R)$ at $R=R_{1}$ and to set
the value and first and second derivatives of $t(R)$ to zero at
$R=R_{\text{cut}}$.

####  VII.4.1 Non-orthogonal representation at $T_{e}\geq 0$

The fast quantum mechanical molecular dynamics scheme, Eqs. (29) and (30) or
Eqs. (8) and (9), using self-consistent tight-binding theory in its non-
orthogonal formulation is given by

$\begin{array}[]{l}{\displaystyle M_{k}{\ddot{R}}_{k}=-2Tr\left[\varrho
H_{R_{k}}\right]}\\\ {}{}\\\
{\displaystyle+\frac{1}{2}\sum_{i,j}q_{i}q_{j}\frac{\partial\gamma_{ij}}{\partial
R_{k}}+\sum_{i,j}q_{i}\gamma_{ij}\left.{\frac{\partial q_{j}}{\partial
R_{k}}}\right\rvert_{\varrho}}\\\ {}{}\\\ {\displaystyle+Tr[(S^{-1}H[{\bf
q}]\varrho+\varrho H[{\bf q}]S^{-1})S_{R_{k}}]-\frac{\partial E_{\rm
pair}[{\bf R}]}{\partial R_{k}},}\\\ \end{array}$ (50)

and

${\ddot{n}}_{i}=\omega^{2}\left(q_{i}-n_{i}\right),$ (51)

where

$H_{R_{k}}=\left.{\frac{\partial H}{\partial R_{k}}}\right\rvert_{\varrho}$
(52)

and

$S_{R_{k}}=\frac{\partial S}{\partial R_{k}}.$ (53)

The partial derivatives of $q_{j}$ and $H$ in Eqs. (50) and (52) are with
respect to a constant density matrix $\varrho$ in its non-orthogonal form,
i.e. including an $S$ dependence of $q_{j}$,

$\left.{\frac{\partial q_{j}}{\partial
R_{k}}}\right\rvert_{\varrho}=2\sum_{\alpha\in j}\left(\varrho
S_{R_{k}}\right)_{j\alpha,j\alpha}.$ (54)

The total energy is given by

$\begin{array}[]{l}{\displaystyle E_{\rm
tot}=\frac{1}{2}\sum_{k}M_{k}{\dot{R}}_{k}^{2}+2\sum_{i\in\rm
occ}\varepsilon_{i}}\\\ {}{}\\\
{\displaystyle-\frac{1}{2}\sum_{i,j}q_{i}q_{j}\gamma_{ij}+E_{\rm pair}[{\bf
R}]-T_{e}{\cal S}[{\bf R};{\bf n}]},\end{array}$ (55)

with the entropy contribution to the free energy approximated by

${\cal S}[{\bf R};n]\approx-2k_{\rm
B}\sum_{i}\left\\{f_{i}\ln(f_{i})+(1-f_{i})\ln(1-f_{i})\right\\}.$ (56)

Here $f_{i}=f_{i}[n]$ are the eigenstates of the Fermi operator expansion
$\varrho^{\perp}[n]$ of $H^{\perp}[{\bf n}]$ in Eq. (44).

####  VII.4.2 Orthogonal representation at $T_{e}=0$

For orthogonal formulations, i.e. when $S=I$, and at zero electronic
temperature, $T_{e}=0$, we will base our dynamics on the equations of motion
in Eqs. (37) and (38). In this case the fast quantum mechanical molecular
dynamics scheme, Eqs. (50)-(51), is given by

$\begin{array}[]{l}{\displaystyle M_{k}{\ddot{R}}_{k}=-2Tr\left[\varrho
H_{R_{k}}\right]+\frac{1}{2}\sum_{i,j}\left(n_{i}n_{j}\frac{\partial\gamma_{ij}}{\partial
R_{k}}\right)}\\\ {\displaystyle-\frac{\partial E_{\rm pair}[{\bf
R}]}{\partial R_{k}}},\end{array}$ (57)
${\ddot{n}}_{i}=\omega^{2}\left(q_{i}-n_{i}\right),$ (58)

where

$\\{H_{R_{k}}[{\bf n}]\\}_{i\alpha,j\beta}=\frac{\partial
h_{i\alpha,j\beta}}{\partial
R_{k}}+\sum_{l}n_{l}\frac{\partial\gamma_{il}}{\partial
R_{k}}\delta_{ij}\delta_{\alpha\beta}.$ (59)

The density matrix is given directly from the step function of the
Hamiltonian, $\varrho=\theta(\mu_{0}I-H[{\bf n}])$, without any de-
orthogonalization that requires the calculation of the inverse factorization
of the overlap matrix, Eq. (47). The constant of motion, $E_{\rm tot}$, is
approximate by

$\begin{array}[]{l}{\displaystyle E_{\rm
tot}=\frac{1}{2}\sum_{k}M_{k}{\dot{R}}_{k}^{2}+2\sum_{i\in\rm
occ}\varepsilon_{i}}\\\ {}{}\\\
{\displaystyle-\frac{1}{2}\sum_{i,j}(2n_{i}-q_{i})n_{j}\gamma_{ij}+E_{\rm
pair}[{\bf R}]}.\end{array}$ (60)

####  VII.4.3 General remarks

Apart from the first few initial molecular dynamics time steps, where we apply
a high degree of self-consistent-field convergence and set ${\bf n=q}$, no
ground state self-consistent-field optimization is required. The density
matrix, $\varrho$, and the Hamiltonian, $H$, necessary in the force
calculations (and for the total energy) are calculated only once per time step
in the orthogonal case with one additional construction of the Hamiltonian
required in non-orthogonal simulations. The numerical integration of the
equations of motion in Eq. (14) is performed with the velocity Verlet scheme
and in Eq. (15) with the modified Verlet scheme in Eq. (16) as described in
Ref. Niklasson et al. (2009). For the examples presented here we used the
modified Verlet scheme including dissipation ($\alpha>0$) with $K=5$ and the
constant $\kappa=\delta t^{2}\omega^{2}$ as given in Ref. Niklasson et al.
(2009) was rescaled by a factor $1/2$ in all examples except for one of the
test cases in Fig. 6.

## References

  * Karplus and McCammon (2002) M. Karplus and J. A. McCammon, Nat. Struct. Biol. 9, 646 (2002).
  * 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.
  * Wang and Karplus (1973) I. S. Y. Wang and M. Karplus, J. Am. Chem. Soc. 95, 8160 (1973).
  * Leforestier (1978) C. Leforestier, J. Chem. Phys. 68, 4406 (1978).
  * Roothaan (1951) C. C. J. Roothaan, Rev. Mod. Phys. 23, 69 (1951).
  * McWeeny (1960) R. McWeeny, Rev. Mod. Phys. 32, 335 (1960).
  * 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. B 140, A1133 (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. M. Dreizler and K. U. Gross, _Density-functional theory_ (Springer Verlag, Berlin Heidelberg, 1990).
  * Feynman (1939) R. P. Feynman, Phys. Rev. 56, 367 (1939).
  * Pulay and Fogarasi (2004) P. Pulay and G. Fogarasi, Chem. Phys. Lett. 386, 272 (2004).
  * Niklasson et al. (2006) A. M. N. Niklasson, C. J. Tymczak, and M. Challacombe, Phys. Rev. Lett. 97, 123001 (2006).
  * Kühne et al. (2006) T. D. Kühne, M. Krack, F. R. Mohamed, and M. Parrinello, Phys. Rev. Lett. 98, 066401 (2006).
  * Ehrenfest (1927) P. Ehrenfest, Z. Phys. 45, 455 (1927).
  * Alonso (2008) J. L. Alonso, X. Andrade, P. Echenique, F. Falceto, D. Prada-Garcia, A. Rubio, Phys. Rev. Lett. 101, 096403 (2008).
  * Jakowski (2009) J. Jakowski, and K. Morokuma, J. Chem. Phys. 130, 224106 (2009).
  * Car and Parrinello (1985) R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
  * Tuckerman (2002) M. Tuckerman, J. Phys.:Condens. Matter 50, 1297 (2002).
  * Hartke (1992) B. Hartke, and E. A. Carter, Chem. Phys. Lett. 189, 358 (1992).
  * Schlegel (1992) H. B. Schlegel, J. M. Millam, S. S. Iyengar, G. A. Voth, A. D. Daniels, G. Scusseria, and M. J. Frisch, J. Chem. Phys. 114, 9758 (2001).
  * Herbert (1992) J. Herbert, and M. Head-Gordon, J. Chem. Phys. 121, 11542 (2004).
  * Kirchner (2012) B. Kirchner J. di Dio Philipp, and J. Hutter, Top. Curr. Chem. 307, 109 Springer Verlag, Berlin Heidelberg, (2012).
  * Goedecker (1999) S. Goedecker, Rev. Mod. Phys. 71, 1085 (1999).
  * Bowler and Miyazaki (2012) D. R. Bowler and T. Miyazaki, Rep. Prog. Phys. 75, 036503 (2012).
  * Niklasson (2002) A. M. N. Niklasson, Phys. Rev. B 66, 155115 (2002).
  * Cawkwell and Niklasson (2012) M. J. Cawkwell and A. M. N. Niklasson, J. Chem. Phys. (accepted for publication).
  * Niklasson (2003) A. M. N. Niklasson, Phys. Rev. B 68, 233104 (2003).
  * Mauri and Galli (1994) F. Mauri, and G. Galli, Phys. Rev. B 50, 4316 (1994).
  * Horsfield et al. (1996) A. P. Horsfield, A. M. Bratkovsky, M. Fern, D. G. Pettifor, and M. Aoki, Phys. Rev. B 54, 12694 (1996).
  * Tushida (2008) E. Tushida, J. Phys.: Condens. Matter 20, 294212 (2008).
  * Shimojo et al. (2008) F. Shimojo, R. K. Kalia, A. Nakono, and P. Vashista, Phys. Rev. B 77, 085103 (2008).
  * 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).
  * Niklasson (2008) A. M. N. Niklasson, Phys. Rev. Lett. 100, 123004 (2008).
  * Niklasson et al. (2009) A. M. N. Niklasson, P. Steneteg, A. Odell, N. Bock, M. Challacombe, C. J. Tymczak, E. Holmström, G. Zheng, and V. Weber, J. Chem. Phys. 130, 214109 (2009).
  * Steneteg et al. (2010) P. Steneteg, I. A. Abrikosov, V. Weber, and A. M. N. Niklasson, Phys. Rev. B 82, 075110 (2010).
  * Zheng et al. (2011) G. Zheng, A. M. N. Niklasson, and M. Karplus, J. Chem. Phys. 135, 044122 (2011).
  * Niklasson et al. (2011) A. M. N. Niklasson, P. Steneteg, and N. Bock, J. Chem. Phys. 135, 164111 (2011).
  * Weinert and Davenport (1992) M. Weinert and J. W. Davenport, Phys. Rev. B 45, R13709 (1992).
  * Wentzcovitch et al. (1992) R. M. Wentzcovitch, J. L. Martins, and P. B. Allen, Phys. Rev. B 45, R11372 (1992).
  * Niklasson (2008b) A. M. N. Niklasson, J. Chem. Phys. 129, 244107 (2008b).
  * Dederichs and Zeller (1983) P. H. Dederichs and R. Zeller, Phys. Rev. B 28, 5262 (1983).
  * Leimkuhler and Reich (2004) B. Leimkuhler and S. Reich, _Simulating Hamiltonian Dynamics_ (Cambridge University Press, 2004).
  * Odell et al. (2009) A. Odell, A. Delin, B. Johansson, N. Bock, M. Challacombe, and A. M. N. Niklasson, 131, 244106 (2009), J. Chem. Phys.
  * Odell et al. (2011) A. Odell, A. Delin, B. Johansson, M. Cawkwell, and A. M. N. Niklasson, 135, 224105 (2011), J. Chem. Phys.
  * Sutton et al. (1988) A. P. Sutton, M. W. Finnis, D. G. Pettifor, and Y. Ohta, J. Phys. C: Solid State Phys. 21, 35 (1988).
  * Finnis et al. (1998) M. W. Finnis, A. T. Paxton, M. Methfessel, and M. van Schilfgaarde, Phys. Rev. Lett. 81, 5149 (1998).
  * Elstner et al. (1998) M. Elstner, D. Poresag, G. Jungnickel, J. Elstner, M. Haugk, T. Frauenheim, S. Suhai, and G. Seifert, Phys. Rev. B 58, 7260 (1998).
  * Finnis (2003) M. Finnis, _Interatomic forces in condensed matter_ (Oxford University Press, 2003).
  * Sanville et al. (2010) E. Sanville, N. Bock, W. M. Challacombe, A. M. N. Niklasson, M. J. Cawkwell, D. M. Dattelbaum, and S. Sheffield, _Proceedings of the Fourteenth International Detonation Symposium_ (Office of Naval Research, Arlington VA, ONR-351-10-185, 2010), pp. 91–101.
  * Rubensson and Rudberg (2011) E. H. Rubensson and E. Rudberg, J. Phys.: Condes. Matter 23, 075502 (2011).