# Mass-Zero constrained dynamics for simulations based on orbital-free density
functional theory

A. Coretti  Faculty of Physics, University of Vienna, Kolingasse 14-16, 1090
Vienna, Austria  Department of Mathematical Sciences, Politecnico di Torino,
I-10129 Torino, Italy  Centre Européen de Calcul Atomique et Moléculaire
(CECAM), Ecole Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland
T. Baird  Centre Européen de Calcul Atomique et Moléculaire (CECAM), Ecole
Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland  R. Vuilleumier
PASTEUR, Département de chimie, École normale supérieure, PSL University,
Sorbonne Université. CNRS, 75005 Paris, France  S. Bonella
<EMAIL_ADDRESS>Centre Européen de Calcul Atomique et Moléculaire
(CECAM), Ecole Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland

###### Abstract

A new algorithm for efficient and fully time-reversible integration of first-
principles molecular dynamics based on orbital-free density functional theory
(OFDFT) is presented. The algorithm adapts to this nontrivial case the
recently introduced Mass-Zero (MaZe) constrained dynamics. The formalism
ensures that full adiabatic separation is enforced between nuclear and
electronic degrees of freedom and, consequently, that the exact Born-
Oppenheimer probability for the nuclei is sampled. Numerical integration of
the MaZe dynamics combines standard molecular dynamics algorithms, e.g. Verlet
or velocity Verlet, with the SHAKE method to impose the minimum conditions on
the electronic degrees of freedom as a set of constraints. The developments
presented in this work, that include a bespoke adaptation of the standard
SHAKE algorithm, ensure that the quasilinear scaling of OFDFT is preserved by
the new method for a broad range of kinetic and exchange-correlation
functionals, including nonlocal ones. The efficiency and accuracy of the
approach is demonstrated via calculations of static and dynamic properties of
liquid sodium in the constant energy and constant temperature ensembles.

First-principles Molecular Dynamics, Orbital-free Density Functional Theory,
Constrained Molecular Dynamics

††preprint: AIP/123-QED

##  I Introduction

In this paper, we detail the use of the Mass-Zero (MaZe) constrained dynamics
for the interesting case of first-principles molecular dynamics based on
Orbital-Free Density Functional Theory (OFDFT).

MaZe is a general method for simulating the evolution of a set of physical
degrees of freedom (dofs) whose dynamics depends on parameters subject to
given conditions.coretti:2018b ; bonella:2020  In MaZe these parameters are
treated as degrees of freedom of zero inertia in an extended dynamical system
and the conditions are imposed as constraints. The method is formulated in the
Lagrangian framework, enabling the key properties of the approach, namely full
time reversibility and exact sampling of the probability density of the
physical dofs, to emerge naturally from a fully consistent dynamical and
statistical viewpoint. An important application domain for MaZe is adiabatic
dynamics. In the following we shall focus, in particular, on the first-
principles adiabatic evolution of a set of coupled ionic (slow) and electronic
(fast) dofs. The energy of this system is a function of the ionic positions
and a functional of the electronic density. In OFDFT-based simulations, the
density is represented in a basis, e.g. plane waves, and the expansion
coefficients play the role of the electronic dofs. The well-separated
timescales of the motion of the slow and fast dofs enable one to adopt the
Born-Oppenheimer approximation for the dynamics. In this approximation,
classical evolution equations govern the ionic motion, but the forces acting
on these dofs depend parametrically on the values of the fast variables. These
values, in turn, are determined by the condition that the energy is at a
minimum with respect to the fast dofs for each ionic configuration along the
trajectory.

Typical numerical methods for first-principles MD can be grouped into two main
families. The first directly implements Born-Oppenheimer dynamics: standard
classical MD for the ions is combined with minimization of the energy with
respect to the electronic dofs, for example via nonlinear conjugate gradient
aguado:2003a ; jahn:2004  or via ad hoc evolutions as in the always stable
predictor-corrector method.kolafa:2004 ; genzer:2004  The second implements
Car-Parrinello MD:car:1985  the fast dofs are treated as dynamical variables
(of finite mass) by defining an appropriate Lagrangian and the minimum is
approximately tracked by the overall dynamics.pastore:1991  A hybrid method
that introduces new dynamical variables, in addition to those directly
associated with the electronic dofs, and uses them to provide good starting
values for energy minimizers such as the conjugate-gradient method, was also
proposed in the extended Lagrangian Born-Oppenheimer (XL-BOMD)
framework.niklasson:2006 ; niklasson:2020 ; niklasson:2021  These methods are
routinely used and have produced remarkable results, but they suffer from
problems that motivate further developments. In particular, full convergence
of minimization based approaches can be too expensive for large-scale
simulations of condensed-phase systems.pounds:2009  Incomplete convergence has
serious consequences on the evolution, including instabilities,pacaud:2018
violation of adiabaticity remler:1990 ; pulay:2004  and therefore energy
drifts and breaking of time reversibility of the propagation. Violation of
adiabaticity affects also Car-Parrinello MD aguado:2003a  and forces either
the use of thermostats (which, however, lead to questions on the statistical
ensemble sampled) or of a very small mass associated to the auxiliary
variables (which, however, imposes a very small timestep for the numerical
propagation). Full time reversibility or ease of portability of the algorithm
is problematic for the predictor-corrector and XL-BOMD methods, that have
system dependent parameters or forces that must be tuned for each
calculation.kolafa:2004 ; genzer:2004  Indeed, the violation of strict
adiabaticity is intrinsic to extended-Lagrangian methods as long as the mass
of the auxiliary variables is finite.payne:1986 ; pastore:1991 ; payne:1992 ;
marx:2012-book

MaZe provides an alternative approach that eliminates or mitigates the
problems mentioned above. It employs the framework of constrained dynamics
goldstein:2002-book ; ciccotti:1986  to ensure fully adiabatic propagation
and, consequently, correct sampling of the target probability for the physical
dofs. Using an adapted SHAKE algorithm,ryckaert:1977  it enables efficient
satisfaction of the minimum condition with a level of precision close to the
numerical limit.coretti:2018b ; bonella:2020 ; coretti:2020a  The SHAKE-based
integration algorithm is symplectic leimkuhler:2015-book  and fully time
reversible, leading to long-time stability of the simulation with a timestep
determined by the typical frequencies of the motion of the slow dofs. Recent
applications of MaZe that illustrate these properties cover the simulation of
classical models for polarizable systems,coretti:2018b ; coretti:2020a ;
marin-lafleche:2020  including the nontrivial case of classical polarizable
systems in external magnetic field,girardier:2021  and proof-of-principle
OFDFT calculations that employed simple energy functionals.bonella:2020

In the following, we considerably extend our previous exploratory work for
MaZe in OFDFT simulations by tackling two important issues. Firstly, as
discussed in detail in Sec. II, the MaZe OFDFT constrained variables are the
complex coefficients in a Fourier-series expansion of the electronic density.
This implies that the Lagrange multipliers associated to our problem are also
complex. The generalization of the SHAKE algorithm to this case is presented
in this work. The second, critical problem that we address stems from the
specific nature of the minimum constraint. As a rule, SHAKE has been applied
to relatively simple constraints, involving only small subsets of degrees of
freedom, e.g. fixing bond lengths in molecules where the functional form of
the constraint is quadratic in the coordinates of the few atoms involved in
each bond. This local structure is key for the rapid convergence of the
iterative algorithm (almost) universally adopted for solving the constrained
evolution of the system.ryckaert:1977  The minimum constraint appearing in
MaZe, on the other hand, is typically a global and highly nonlinear function
of many auxiliary degrees of freedom that, for general energy functionals,
cannot be solved efficiently with the standard algorithm. In this work then,
we introduce an alternative procedure, based on work of Weinbach and
Elber,weinbach:2005  to solve the constraint equations and prove its
remarkable efficacy. This development is essential in view of MaZe deployment
to large-scale simulations based on OFDFT with functionals of arbitrary
complexity, including nonlocal kinetic and exchange-correlation
functionals.wang:1992 ; wang:1999

The paper is organized as follows. We start by discussing in detail the MaZe
formalism for the specific case of OFDFT first-principles MD at constant
number of particles, volume and energy. We then show how the constrained
dynamics is adapted to the case of complex Lagrange multipliers. In Sec. III,
we discuss first the solution of the constraints via the standard algorithm
(see Sec. III.1) and highlight its high numerical cost for generic kinetic and
exchange-correlation functionals. In Sec. III.2 the new, numerically efficient
and general algorithm for the solution of the constraint equations is then
introduced. The accuracy and efficiency enabled by these developments are
illustrated in the Results section for constant energy and constant
temperature systems. Here, we first validate the algorithm via benchmark
calculations on liquid sodium and then compare the numerical performance (both
in terms of number of operations per timestep and in clocking of the timestep)
of the new and standard SHAKE algorithms. In particular, we demonstrate that,
for commonly adopted energy functionals at the GGA level, the algorithm scales
quasilinearly with system size.

##  II MaZe dynamics for orbital-free DFT

OFDFT was proposedpearson:1993  as a computationally convenient alternative to
DFT dynamics based on the use of Kohn-Sham orbitals. marx:2012-book  The main
advantage of OFDFT is that, in contrast to Kohn-Sham-based DFT, it enables a
representation of the total energy that does not require auxiliary orthonormal
orbitals. This entails that the overall cost of the calculation per timestep
is quasilinear in the system’s volume $\Omega$,
($\mathcal{O}(\Omega\log\Omega)$) and memory requirements conveniently scale
linearly as $\mathcal{O}(\Omega)$. Unfortunately, developing reliable and
transferable kinetic-energy functionals for general systems (and for different
phases of matter) is an open challenge ligneres:2005-inbook  that prevents
universal use of the method. In spite of this well-known issue, OFDFT remains
a useful approach, often enabling the simulation of system sizes out of reach
of alternative methods.ligneres:2005-inbook ; witt:2018  The method is also
implemented in different high performance codes for the simulation of
materials, such as GPAW mortensen:2005 ; Enkovaara:2010  or PROFESS.chen:2015
Very recently, the use of machine learning methods to obtain flexible and
accurate kinetic-energy functionals meyer:2020 ; ghasemi:2021  has also led to
a renewed interest in OFDFT.

In OFDFT, the ground state electronic energy is written as a function of the
$N$ ionic positions, $\bm{R}\equiv\\{\bm{R}_{1},...,\bm{R}_{N}\\}$, and as a
functional of the electronic density $n$

$\begin{array}[]{cc}E[\bm{R},n]&\equiv
E_{\text{KE}}[n]+E_{\text{Hart}}[n]+E_{\text{xc}}[n]\\\
&+E_{\text{ext}}[\bm{R},n]+E_{\text{Coul}}(\bm{R}),\end{array}$ (1)

where $E_{\text{KE}}[n]$ is the kinetic energy of a noninteracting electron
gas at the density of the physical system, and $E_{\text{Hart}}[n]$ and
$E_{\text{xc}}[n]$ are the Hartree and exchange-correlation energies,
respectively. $E_{\text{Coul}}(\bm{R})$ is the classical Coulomb energy of the
ions and $E_{\text{ext}}[\bm{R},n]$ is the interaction energy between the
electrons and the ions. In MD simulations, the density is typically
represented in the plane-wave basis

$n(\bm{r})=\sum_{\bm{G}}\tilde{n}_{\bm{G}}e^{i\bm{G}\cdot\bm{r}},$ (2)

where $\bm{G}=\frac{2\pi}{L}(i,j,k)$, with $L$ the size of the cubic box and
$i$, $j$, $k$ integers. This expansion is truncated with
$\tilde{n}_{\bm{G}}=0$ if any component of the wavevector is larger than a
cutoff, $|G^{l}|>G_{\text{c}}$, $l=x,y\text{ or }z$. Since the density in real
space $n(\bm{r})$ is, by definition, a real function, the Fourier coefficients
$\tilde{n}_{\bm{G}}$ (typically complex quantities), satisfy the relationship
$\tilde{n}_{-\bm{G}}=\tilde{n}^{\star}_{\bm{G}}$. Once the Fourier
representation is adopted, the ground state energy becomes a function of the
ionic positions and the expansion coefficients.

First-principles OFDFT molecular dynamics relies on the Born-Oppenheimer
approximation: the forces on the ions are computed at each timestep as the
gradient of the ground state energy with respect to the ionic positions, with
the energy computed at an electronic density that minimizes its value. In the
plane wave representation of the density, this minimum condition is equivalent
to imposing that the derivatives of the ground state energy with respect to
the coefficients of the Fourier expansion are null. In taking these
derivatives, we have the choice to adopt as independent variables the real and
imaginary parts of $\tilde{n}_{\bm{G}}$ or to consider instead
$\tilde{n}_{\bm{G}}$ and $\tilde{n}^{\star}_{\bm{G}}$. Indeed, the definition
of the complex derivatives (also known as Wirtinger derivatives)

$\displaystyle\frac{\partial E}{\partial\tilde{n}^{\star}_{\bm{G}}}$
$\displaystyle=\frac{1}{2}\left(\frac{\partial
E}{\partial\text{Re}[\tilde{n}_{\bm{G}}]}+\mathrm{i}\frac{\partial
E}{\partial\text{Im}[\tilde{n}_{\bm{G}}]}\right),$ (3)
$\displaystyle\frac{\partial E}{\partial\tilde{n}_{\bm{G}}}$
$\displaystyle=\frac{1}{2}\left(\frac{\partial
E}{\partial\text{Re}[\tilde{n}_{\bm{G}}]}-\mathrm{i}\frac{\partial
E}{\partial\text{Im}[\tilde{n}_{\bm{G}}]}\right),$

(where $\text{Re}[\tilde{n}_{\bm{G}}]$ and $\text{Im}[\tilde{n}_{\bm{G}}]$ are
the real and imaginary parts of $\tilde{n}_{\bm{G}}$) ensures the equivalence
of these choices. In the following, we adopt $\tilde{n}_{\bm{G}}$ and
$\tilde{n}^{\star}_{\bm{G}}$ as independent variables. Furthermore, from now
on all Fourier coefficients are taken in the interval
$\bm{G}\in[\bm{0},\bm{G}_{\text{max}}]$, a shorthand notation to indicate that
wavevectors span one-half space (we denote $N_{G}$ the number of such
vectors). With these choices, the Born-Oppenheimer conditions become

$\displaystyle\sigma_{\bm{G}_{\alpha}}(\bm{R},\tilde{\bm{n}},\tilde{n}_{\bm{G}_{0}})=\frac{\partial
E(\bm{R},\tilde{\bm{n}},\tilde{n}_{\bm{G}_{0}})}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}=0\
\ \text{for}\ \ \alpha\in[1,N_{G}];$ (4)
$\displaystyle\sigma^{\star}_{\bm{G}_{\alpha}}(\bm{R},\tilde{\bm{n}},\tilde{n}_{\bm{G}_{0}})=\frac{\partial
E(\bm{R},\tilde{\bm{n}},\tilde{n}_{\bm{G}_{0}})}{\partial\tilde{n}_{\bm{G}_{\alpha}}}=0\
\ \text{for}\ \ \alpha\in[0,N_{G}].$

In the equation above, we have introduced, for future convenience, the vector

$\tilde{\bm{n}}=\\{\tilde{n}_{\bm{G}_{1}},\ldots,\tilde{n}_{\bm{G}_{N_{G}}},\tilde{n}^{\star}_{\bm{G}_{1}},\ldots,\tilde{n}^{\star}_{\bm{G}_{N_{G}}}\\},$
(5)

that contains the $2N_{G}$ Fourier coefficients other than the one
corresponding to $\bm{G}_{0}=(0,0,0)$. This latter coefficient is associated
to the additional condition that the total number of electrons
$N_{\text{el}}$, must be conserved along the dynamics. In the plane wave
basis, this condition is expressed as

$\tilde{n}_{\bm{G}_{0}}\Omega=N_{\text{el}},$ (6)

where $\Omega$ is the volume of the simulation box.

The MaZe evolution for OFDFT first-principles dynamics is obtained by
promoting the coefficients $\tilde{n}_{\bm{G}_{\alpha}}$ and
$\tilde{n}^{\star}_{\bm{G}_{\alpha}}$ to the role of dynamical variables. The
evolution equations, obtained from a Lagrangian formulation of the dynamics of
the ions and these additional variables, are then complemented with the
conditions in Eqs. (4) and (6), which are interpreted as constraints. In the
following, we impose the conservation of the number of electrons by enforcing
Eq. (6) and the minimum condition on $\tilde{n}_{\bm{G}_{0}}$ as an initial
condition. We then exclude this quantity from the set of auxiliary dynamical
variables and consider it as a constant in the simulations. To simplify the
notation, we also drop the explicit indication of the dependence of the ground
state energy and of the constraints on $\tilde{n}_{\bm{G}_{0}}$.

To proceed, we assign (temporarily) a finite (scalar) inertia, $\mu$, to the
additional variables $\tilde{\bm{n}}$, and define the following Lagrangian for
the system111We define the scalar product of two complex vectors,
$\bm{a},\bm{b}\in\mathbb{C}^{N}$ as
$\bm{a}^{\star}\cdot\bm{b}=\sum\limits_{i=1}^{N}a_{i}^{\star}b_{i}$.

$L=\frac{1}{2}\sum_{I}^{N}M_{I}\dot{\bm{R}}^{2}_{I}+\frac{\mu}{2}\dot{\bm{\tilde{n}}}^{\star}\cdot\dot{\bm{\tilde{n}}}-E(\bm{R},\tilde{\bm{n}}).$
(7)

We then introduce the vector of the $2N_{G}$ constraints

$\bm{\sigma}=\\{\sigma_{\bm{G}_{1}},\ldots,\sigma_{\bm{G}_{N_{G}}},\sigma^{\star}_{\bm{G}_{1}},\ldots,\sigma^{\star}_{\bm{G}_{N_{G}}}\\},$
(8)

and the $2N_{G}$ Lagrange multipliers

$\bm{\lambda}=\\{\lambda_{\bm{G}_{1}},\ldots\lambda_{\bm{G}_{N_{G}}},\lambda^{\star}_{\bm{G}_{1}},\ldots,\lambda^{\star}_{\bm{G}_{N_{G}}}\\},$
(9)

in order to form the constrained Lagrangian

$\mathcal{L}=L+\bm{\sigma}^{\star}\cdot\bm{\lambda}.$ (10)

Note that, at variance with typical situations in molecular dynamics, the
constraints and the Lagrange multipliers are, in general, complex quantities,
making it necessary to generalize standard algorithms for their treatment. The
prescribed form for the Lagrangian respects two important properties. The
constraint term, in fact, can be rewritten as

$\bm{\sigma}^{\star}\cdot\bm{\lambda}=\text{Re}\left(\frac{\partial
E}{\partial\tilde{\bm{n}}}\right)\cdot\text{Re}\left(\bm{\lambda}\right)+\text{Im}\left(\frac{\partial
E}{\partial\tilde{\bm{n}}}\right)\cdot\text{Im}(\bm{\lambda}),$ (11)

where we have introduced the vector of derivatives:

$\displaystyle\frac{\partial E}{\partial\tilde{\bm{n}}}=\left\\{\frac{\partial
E}{\partial\tilde{n}_{\bm{G}_{1}}},\ldots,\frac{\partial
E}{\partial\tilde{n}_{\bm{G}_{N_{G}}}},\frac{\partial
E}{\partial\tilde{n}^{\star}_{\bm{G}_{1}}},\ldots,\frac{\partial
E}{\partial\tilde{n}^{\star}_{\bm{G}_{N_{G}}}}\right\\}.$ (12)

Eq. (11) shows that the scalar product of the constraints and Lagrange
multipliers is real and, therefore, so is the constrained Lagrangian, Eq. (7).
Furthermore, the expression verifies, once more, the equivalence of choosing
$\tilde{n}_{\bm{G}_{\alpha}}$ and $\tilde{n}^{\star}_{\bm{G}_{\alpha}}$ or the
real and imaginary parts of the density coefficients as independent variables.
Indeed, with the latter choice, the equation indicates that, as expected and
consistently, the constraints

$\frac{\partial
E}{\partial\text{Re}[\tilde{n}_{\bm{G}_{\alpha}}]}=0;\qquad\frac{\partial
E}{\partial\text{Im}[\tilde{n}_{\bm{G}_{\alpha}}]}=0,$ (13)

would be satisfied, via the — real — Lagrange multipliers
$\text{Re}\left(\lambda_{\bm{G}_{\alpha}}\right)$ and
$\text{Im}\left({\lambda}_{\bm{G}_{\alpha}}\right)$.

The constrained evolution equations for the MaZe OFDFT system are then given
by

$\displaystyle M_{I}\ddot{\bm{R}}_{I}$
$\displaystyle=-\nabla_{\bm{R}_{I}}E(\bm{R},\tilde{\bm{n}})+(\nabla_{\bm{R}_{I}}\bm{\sigma}^{\star})\cdot\bm{\lambda};$
(14) $\displaystyle\mu\ddot{{\tilde{\bm{n}}}}$ $\displaystyle=-\frac{\partial
E(\bm{R},\tilde{\bm{n}})}{\partial\tilde{\bm{n}}^{\star}}+\frac{\partial\bm{\sigma}^{\star}}{\partial\tilde{\bm{n}}^{\star}}\cdot\bm{\lambda}$
$\displaystyle\equiv\frac{\partial
E(\bm{R},\tilde{\bm{n}})}{\partial\tilde{\bm{n}}^{\star}}+\Sigma^{\dagger}\cdot\bm{\lambda},$

where we have introduced the matrix
$\Sigma^{\dagger}\equiv\frac{\partial\bm{\sigma}^{\star}}{\partial\tilde{\bm{n}}^{\star}}$,
which is the Hermitian conjugate of

$\displaystyle\\{\Sigma\\}_{\alpha,\beta}=\\{\Sigma\\}_{N_{G}+\alpha,N_{G}+\beta}^{\star}$
$\displaystyle=\frac{\partial\sigma_{\bm{G}_{\alpha}}}{\partial\tilde{n}_{\bm{G}_{\beta}}};$
(15)
$\displaystyle\\{\Sigma\\}_{\alpha,N_{G}+\beta}=\\{\Sigma\\}_{N_{G}+\alpha,\beta}^{\star}$
$\displaystyle=\frac{\partial\sigma_{\bm{G}_{\alpha}}}{\partial\tilde{n}^{\star}_{\bm{G}_{\beta}}}.$

The evolution equations for the additional variables can be rearranged by
first observing that the vector derivative of the energy with respect to the
expansion coefficients is the null vector since it corresponds to the
constraints. Dividing the evolution equation for the additional variables by
the inertia we obtain

$\displaystyle M_{I}\ddot{\bm{R}}_{I}$
$\displaystyle=-\nabla_{\bm{R}_{I}}E(\bm{R},\tilde{\bm{n}})+(\nabla_{\bm{R}_{I}}\bm{\sigma}^{\star})\cdot\bm{\lambda};$
(16) $\displaystyle\ddot{\tilde{\bm{n}}}$
$\displaystyle=\frac{1}{\mu}\Sigma^{\dagger}\cdot\bm{\lambda}.$

The key step in the MaZe approach follows by considering the limit $\mu\to 0$
for the second equation above and imposing that this limit is taken keeping
the accelerations of $\tilde{\bm{n}}$ finite. This implies that, for every
component of the vector of the Lagrange multipliers, the ratio
$\gamma_{\bm{G}_{\beta}}=\lim_{\mu\to 0}\frac{\lambda_{\bm{G}_{\beta}}}{\mu}$
must remain finite, entailing that the Lagrange multipliers are proportional
to $\mu$. In the limit of zero mass for the auxiliary variables, then, the
constraint forces acting on the ionic degrees of freedom vanish and the system
becomes

$\displaystyle M_{I}\ddot{\bm{R}}_{I}$
$\displaystyle=-\nabla_{\bm{R}_{I}}E(\bm{R},\tilde{\bm{n}});$ (17)
$\displaystyle\ddot{\tilde{\bm{n}}}$
$\displaystyle=\Sigma^{\dagger}\cdot\bm{\gamma}.$

with $\bm{\gamma}=\lim_{\mu\to 0}\frac{\bm{\lambda}}{\mu}$.

Eq. (17) defines the MaZe dynamics for OFDFT. Some important properties must
be noted. Firstly, since the Lagrange multipliers $\bm{\lambda}$ go to zero
with the auxiliary mass $\mu$, the evolution of the physical variables is not
affected by the constraint forces. Secondly, the evolution of
$\tilde{\bm{n}}$, controlled only by the constraint forces, satisfies by
construction the minimum conditions at all times. These conditions are thus
automatically fulfilled also in the first of Eq. (17) which is then fully
equivalent to Born-Oppenheimer dynamics. This implies that, by rigorously
enforcing the mass-zero limit for the auxiliary variables, MaZe provides an
evolution that leads to the exact adiabatic dynamics for the physical
variables. This has the important consequence that the probability density
sampled from the dynamical system for the ions is exactly the one typically
associated with Born-Oppenheimer dynamics.bonella:2020 ; coretti:2020a  MaZe
is then free of sampling errors typically encountered in schemes with finite
values of the auxiliary variables’ mass. From the point of view of the
algorithm, the first equation can be integrated via any standard MD algorithm
(e.g. Verlet) with a timestep determined only by the force acting on the ions.
At each timestep, the Lagrange multipliers $\bm{\gamma}$, that appear as
unknown, time-dependent parameters, must also be determined. This is done
adopting the SHAKE strategy. ryckaert:1977 ; ciccotti:1986  In the next
section, details of the algorithm for OFDFT MaZe dynamics are provided.

##  III MaZe algorithm for orbital-free DFT

The numerical integration of the dynamical system in Eq. (17) proceeds as
follows. The evolution of the ionic degrees of freedom is performed via
standard classical MD algorithms. Here we discuss the case of Verlet
propagation but there is no difficulty in adopting, for example, the velocity
version of the scheme. Thus

$\bm{R}_{I}(t+\delta t)=2\bm{R}_{I}(t)-\bm{R}_{I}(t-\delta t)-\frac{\delta
t^{2}}{M_{I}}\nabla_{\bm{R}_{I}}E\bigl{(}\bm{R}(t),\tilde{\bm{n}}(t)\bigr{)},$
(18)

where $\delta t$ is the timestep of the simulation. This algorithm requires,
at each step, the computation of the forces acting on the ions after
advancement of their positions. These forces

$\bm{F}_{I}(t+\delta t)=-\nabla_{\bm{R}_{I}}E\bigl{(}\bm{R}(t+\delta
t),\tilde{\bm{n}}(t+\delta t)\bigr{)},$ (19)

depend both on the advanced ionic coordinates and on the values of the
$\tilde{\bm{n}}(t+\delta t)$. The latter can be formally obtained as

$\displaystyle\tilde{\bm{n}}(t+\delta t;\bm{\gamma})=$ $\displaystyle
2\tilde{\bm{n}}(t)-\tilde{\bm{n}}(t-\delta t)+\delta
t^{2}\Sigma(t)^{\dagger}\cdot\bm{\gamma}$ (20) $\displaystyle=$
$\displaystyle\tilde{\bm{n}}^{P}+\Sigma(t)^{\dagger}\cdot\tilde{\bm{\gamma}}.$

In the second line above, we have implicitly defined
$\tilde{\bm{n}}^{P}=2\tilde{\bm{n}}(t)-\tilde{\bm{n}}(t-\delta t)$, to be
called later the “provisional” value of the coefficient, and we have scaled
the Lagrange multipliers — as usually done in SHAKE — by the square of the
timestep, so $\tilde{\bm{\gamma}}=\delta t^{2}\bm{\gamma}$. The updated
density coefficients depend on the set of Lagrange multipliers that are, at
this stage, unknown. Following the SHAKE strategy, the multipliers are
determined by imposing that the constraints

$\displaystyle\bm{\sigma}\bigl{(}\bm{R}(t+\delta t),\tilde{\bm{n}}(t+\delta
t;\tilde{\bm{\gamma}})\bigr{)}$ $\displaystyle=\bm{0},$ (21)

are satisfied at the positions predicted by the MD algorithm. Eqs. (21) are,
in general, a system of nonlinear equations to be solved for the Lagrange
multipliers. The solution of the nonlinear system is obtained via an adapted
iterative Newton-Raphson scheme in which the values of the Lagrange
multipliers at iteration ${\kappa}+1$ are given by

${\tilde{\bm{\gamma}}^{{\kappa}+1}}={\tilde{\bm{\gamma}}^{\kappa}}-\omega\mathbb{J}^{-1}(\tilde{\bm{n}}^{\kappa})\bm{\sigma}\bigl{(}\bm{R}(t+\delta
t),\tilde{\bm{n}}^{\kappa}\bigr{)}.$ (22)

In the equation above, we have introduced
$\mathbb{J}^{-1}(\tilde{\bm{n}}^{\kappa})$, the inverse of the Jacobian
$\mathbb{J}(\tilde{\bm{n}}^{\kappa})\equiv\frac{\partial\bm{\sigma}(\tilde{\bm{n}}^{\kappa})}{\partial\tilde{\bm{\gamma}}}$
whose matrix elements are given by

$\displaystyle\\{\mathbb{J}(\tilde{\bm{n}}^{\kappa})\\}_{\alpha,\beta}=\\{\mathbb{J}(\tilde{\bm{n}}^{\kappa})\\}_{N_{G}+\alpha,N_{G}+\beta}^{\star}$
$\displaystyle=\frac{\partial\sigma_{\bm{G}_{\alpha}}(\tilde{\bm{n}}^{\kappa})}{\partial\tilde{\gamma}_{\bm{G}_{\beta}}};$
(23)
$\displaystyle\\{\mathbb{J}(\tilde{\bm{n}}^{\kappa})\\}_{\alpha,N_{G}+\beta}=\\{\mathbb{J}(\tilde{\bm{n}}^{\kappa})\\}_{N_{G}+\alpha,\beta}^{\star}$
$\displaystyle=\frac{\partial\sigma_{\bm{G}_{\alpha}}(\tilde{\bm{n}}^{\kappa})}{\partial\tilde{\gamma}^{\star}_{\bm{G}_{\beta}}},$

with $\alpha,\beta=1,...,N_{G}$, and $\tilde{\bm{n}}^{\kappa}$ is the vector
of the electronic density coefficients at iteration $\kappa$. The parameter
$\omega$ in Eq. (22) is a scaling factor that, following common practice in
minimization algorithms, nocedal:2006-book  modulates the update of the
Lagrange multipliers. Previous calculations barth:1995  show that using such a
scaling improves convergence in standard implementations of SHAKE. At each
iteration in the cycle, the Fourier coefficients are updated as

$\displaystyle{\tilde{\bm{n}}^{\kappa+1}}$
$\displaystyle={\tilde{\bm{n}}}^{\kappa}+\Sigma(t)^{\dagger}\cdot\tilde{\bm{\gamma}}^{\kappa},$
(24)

and used in Eq. (22) to update the vector of constraints. The iteration cycle
is typically initialized using null values for the Lagrange multipliers or
assigning them as the converged values of the previous timestep. The process
is stopped when $\max_{\bm{G}_{\alpha}}|{\sigma_{\bm{G}_{\alpha}}}|$ becomes
smaller than a predefined tolerance.

The convergence of the Newton-Raphson algorithm is normally
quadratic.nocedal:2006-book  The numerical efficiency of the method therefore
hinges on the cost of the different quantities and operations necessary to
perform the iterations. Below, we summarize the main characteristics of these
quantities by showing — when necessary — their explicit component-wise
expression, reporting the numerical cost of the key operations in the
algorithm. Detailed assessment of these numerical costs is given in the
Appendices to avoid interrupting the main text with lengthy calculations.

Let us begin with the evaluation of the constraint vector
$\bm{\sigma}\bigl{(}\bm{R}(t+\delta t),\tilde{\bm{n}}^{\kappa}\bigr{)}$. From
Eq. (21), we have to compute, for example,

$\sigma_{\bm{G}_{\alpha}}\bigl{(}\bm{R}(t+\delta
t),\tilde{\bm{n}}^{\kappa}\bigr{)}=\frac{\partial E\bigl{(}\bm{R}(t+\delta
t),\tilde{\bm{n}}^{\kappa}\bigr{)}}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}},$
(25)

where the ground state energy, $E$ is expressed as the sum in Eq. (1). The
$\sigma^{\star}_{\bm{G}_{\alpha}}\bigl{(}\bm{R}(t+\delta
t),\tilde{\bm{n}}^{\kappa}\bigr{)}$ can be obtained via complex conjugation.
Eq. (25) highlights that the calculation of the constraints is identical to
that of the force on the electronic coefficients in Car-Parrinello OFDFT
pearson:1993  or of the derivatives employed in conjugate-gradient
minimization in Born-Oppenheimer OFDFT chen:2015 . Direct evaluation shows
that the derivatives are obtained at an $\Omega\log\Omega$ cost for any form
of the energy functional, and in particular of the kinetic energy functional.
We do not report the explicit expressions for these quantities since they are
standard and available in the literature.wang:1999 ; chen:2015

Secondly, the evolution equations (20) and the update of the density
coefficients (24) require a matrix-vector product
$\Sigma^{\dagger}\cdot\bm{\gamma}$. In Appendix A we show that this can be
performed with $\Omega\log\Omega$ scaling for GGA functionals and in Appendix
B that the same scaling can be achieved for nonlocal functionals. Having a
procedure to compute the matrix-vector product between $\Sigma^{\dagger}$ and
a generic vector, also offers a means to compute the matrix-vector product of
a generic vector with $\mathbb{J}$, exploiting the fact that the Jacobian can
be expressed as
$\mathbb{J}=\Sigma(\tilde{\bm{n}}^{\kappa})\cdot\Sigma(t)^{\dagger}$. The
latter relationship can be proven by considering the definition of the matrix
elements in Eq. (23). For example,

$\displaystyle\\{\mathbb{J}$
$\displaystyle(\tilde{\bm{n}}^{\kappa})\\}_{\alpha,\beta}=\frac{\partial\sigma_{\bm{G}_{\alpha}}}{\partial\tilde{\gamma}_{\bm{G}_{\beta}}}$
$\displaystyle=\sum^{N_{G}}_{\lambda=1}\biggl{[}\frac{\partial\sigma_{\bm{G}_{\alpha}}}{\partial\tilde{n}_{\bm{G}_{\lambda}}}\frac{\partial\tilde{n}_{\bm{G}_{\lambda}}}{\partial\tilde{\gamma}_{\bm{G}_{\beta}}}+\frac{\partial\sigma_{\bm{G}_{\alpha}}}{\partial\tilde{n}^{\star}_{\bm{G}_{\lambda}}}\frac{\partial\tilde{n}^{\star}_{\bm{G}_{\lambda}}}{\partial\tilde{\gamma}_{\bm{G}_{\beta}}}\biggr{]}$
$\displaystyle=\sum^{N_{G}}_{\lambda=1}\biggl{[}\frac{\partial\sigma_{\bm{G}_{\alpha}}(\tilde{\bm{n}}^{\kappa})}{\partial\tilde{n}_{\bm{G}_{\lambda}}}\frac{\partial\sigma^{\star}_{\bm{G}_{\beta}}(t)}{\partial\tilde{n}^{\star}_{\bm{G}_{\lambda}}}+\frac{\partial\sigma_{\bm{G}_{\alpha}}(\tilde{\bm{n}}^{\kappa})}{\partial\tilde{n}^{\star}_{\bm{G}_{\lambda}}}\frac{\partial\sigma^{\star}_{\bm{G}_{\beta}}(t)}{\partial\tilde{n}_{\bm{G}_{\lambda}}}\biggr{]}$
$\displaystyle=\left\\{\Sigma(\tilde{\bm{n}}^{\kappa})\cdot\Sigma(t)^{\dagger}\right\\}_{\alpha,\beta}.$
(26)

In going from the first to the second line in the equations above, the
derivative of the density coefficients with respect to the Lagrange
multipliers was computed from Eq. (20) and the notation indicates that
derivatives of the constraints are computed either at the provisional values
of the coefficients, $\tilde{\bm{n}}^{\kappa}$, or with variables computed at
time $t$. Similar expressions hold for the other blocks of the matrix proving
the result. Finally, the solution of the Newton-Raphson system, Eq. (22),
requires evaluation and inversion of the Jacobian, Eq. (23). These steps are
characteristic of the MaZe approach and nontrivial. The number of elements in
the Jacobian matrix, in fact, scales quadratically with the number of
electronic coefficients and brute force inversion of the matrix would take the
cost of MaZe dynamics to cubic order in CPU operations and to quadratic order
in memory requirements. Indeed, it is well known even in standard applications
of classical constrained MD, that the dimensions of the Jacobian are typically
too large to allow for direct inversion. The commonly adopted strategy is then
to consider the diagonal approximation of the Jacobian (i.e., all off-diagonal
elements are set to zero) enabling inversion at linear cost in the matrix
dimensions. In the next subsection, we discuss in detail this consolidated
approach, that we shall indicate as standard SHAKE and that has been employed
also in previous work with MaZe. For energy functionals beyond the GGA class,
such as those including the more recent and accurate nonlocal kinetic
energies, however, linear scaling of this algorithm fails due to the intrinsic
cost of computing the diagonal elements of the Jacobian. In the following
subsection, we then describe an effective and completely general alternative
implementation of SHAKE that preserves the $\Omega\log\Omega$ scaling of OFDFT
simulations also for nonlocal functionals. Both implementations of SHAKE
discussed in the following avoid propagation of the (numerical) error in the
enforcement of the constraint at different timesteps. Furthermore, they are
symplectic and time reversible,leimkuhler:1994 ; leimkuhler:2004-book  leading
to stable propagation for very long times.

Before moving to the details of the SHAKE implementations, we note that the
algorithm described so far for a system of ions evolving at constant energy
(NVE ensemble) can be applied with very limited changes also to ionic
evolution in other ensembles. In particular, ionic motion at constant
temperature will be considered in the Results section. In the full adiabatic
limit, the overall system is characterized in this case by a steady
nonequilibrium regime in which the electronic degrees of freedom remain
strictly at zero temperature and respect the minimum conditions, while the
ions move at an assigned temperature, $T$. The coupled dynamical system is
then obtained by considering thermostatted evolution equations (e.g. Langevin)
for the ionic degrees of freedom and constraining the electronic coefficients
at each configuration along the ionic trajectory as described above.

###  III.1 Standard SHAKE (STD-SHAKE)

The diagonal approximation of the Jacobian matrix defined in Eq. (23) is given
by

$\displaystyle\\{\mathbb{J}_{D}\\}_{\alpha,\beta}$
$\displaystyle=\frac{\partial\sigma_{\bm{G}_{\alpha}}}{\partial\tilde{\gamma}_{\bm{G}_{\beta}}}\delta_{\alpha,\beta};$
(27) $\displaystyle\\{\mathbb{J}_{D}\\}_{\alpha,N_{G}+\beta}$
$\displaystyle=0;$ $\displaystyle\\{\mathbb{J}_{D}\\}_{N_{G}+\alpha,\beta}$
$\displaystyle=0;$
$\displaystyle\\{\mathbb{J}_{D}\\}_{N_{G}+\alpha,N_{G}+\beta}$
$\displaystyle=\frac{\partial\sigma^{\star}_{\bm{G}_{\alpha}}}{\partial\tilde{\gamma}^{\star}_{\bm{G}_{\beta}}}\delta_{\alpha,\beta},$

leading to approximate Eq. (22) as

${\tilde{\bm{\gamma}}^{{\kappa}+1}}={\tilde{\bm{\gamma}}^{\kappa}}-\omega\mathbb{J}_{D}^{-1}(\tilde{\bm{n}}^{\kappa})\bm{\sigma}\bigl{(}\bm{R}(t+\delta
t),\tilde{\bm{n}}^{\kappa}\bigr{)}.$ (28)

This diagonal approximation trivially decouples the first $N_{G}$ components
of the vectors $\bm{\sigma}$ and $\bm{\gamma}$ from the second $N_{G}$
components. Furthermore, as the second $N_{G}$ equations in the diagonal
system corresponds to the complex conjugate of the first $N_{G}$, we only need
to solve for the $N_{G}$-dimensional vector
$\tilde{\gamma}_{\bm{G}_{\alpha}}$. Eq. (28) is then restricted to $N_{G}$
components with the $N_{G}\times N_{G}$ diagonal matrix
$\left\\{\mathbb{J}_{D}\right\\}_{\alpha,\alpha}=\frac{\partial\sigma_{\bm{G}_{\alpha}}}{\partial\tilde{\gamma}_{\bm{G}_{\beta}}}\delta_{\alpha,\beta}$.
Previous experience shows that the use of such diagonal approximation of the
Jacobian does not hinder the speed of convergence of the Newton-Raphson
process, in particular when — as usual in applications of SHAKE in classical
MD — this matrix is diagonally dominant. Furthermore, in typical classical MD
applications, the evaluation of the diagonal elements of the Jacobian does not
pose a numerical challenge because the constraints have simple forms (e.g.
fixed bond distances in molecules), and involve small subsets of degrees of
freedom (e.g. neighboring atoms in the molecules). This is, unfortunately, not
the case for the constraints considered in this work. To appreciate the effort
required to compute the matrix elements of the Jacobian, let us consider the
different terms in Eq. (III). From the definition of the constraints, we have,
for example,

$\displaystyle\frac{\partial\sigma_{\bm{G}_{\alpha}}(\tilde{\bm{n}}^{\kappa})}{\partial\tilde{n}_{\bm{G}_{\lambda}}}$
$\displaystyle=\frac{\partial^{2}E[\tilde{\bm{n}}^{\kappa}]}{\partial\tilde{n}_{\bm{G}_{\lambda}}\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}};$
(29)
$\displaystyle\frac{\partial\sigma_{\bm{G}_{\alpha}}(\tilde{\bm{n}}^{\kappa})}{\partial\tilde{n}^{\star}_{\bm{G}_{\lambda}}}$
$\displaystyle=\frac{\partial^{2}E[\tilde{\bm{n}}^{\kappa}]}{\partial\tilde{n}^{\star}_{\bm{G}_{\lambda}}\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}},$

with analogous expressions for the other terms. The second derivatives of the
energy with respect to the expansion coefficients, in particular for the
kinetic and exchange-correlation terms, are typically available via libraries
of energy functionals for DFT calculations. In the following, we refer to the
lib_xc lehtola:2018  library that has been used in our implementation of MaZe
OFDFT. Explicit expressions for the second derivatives are more easily
discussed by considering separately the terms in Eq. (1). Let us begin by
observing that the Coulomb interaction between the ions (that does not depend
on the electronic density) and the interaction energy between electrons and
ions, $E_{\rm{ext}}[\bm{R},n]$ (that depends linearly on the electronic
density) do not contribute to the second derivatives. The Hartree term, on the
other hand, is quadratic in the electronic density. Using its expression in
Fourier space, it is immediate to show that

$\frac{\partial^{2}E_{\rm{Hart}}[n]}{\partial\tilde{n}_{\bm{G}_{\beta}}\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}=\frac{4\pi\Omega}{|\bm{G}_{\alpha}|^{2}}\delta_{\alpha,\beta},$
(30)

indicating that this term is diagonal in reciprocal space. The evaluation of
the second derivatives of the exchange-correlation (xc) and kinetic energy
functionals is less straightforward. In the following, we shall consider the
GGA class of approximation for these quantities, so they will depend, in
general, on the electronic density and its gradient. GGA xc functionals are
often used in condensed-phase applications both in orbital-free and Kohn-Sham
DFT,perdew:1992  indicating that the algorithm discussed here has a broad
interest. Within OFDFT, GGA kinetic energy functionals include the Thomas-
Fermi-von Weizsäcker vonweizsacker:1935  (TFvW) and Perrot perrot:1994
approximations discussed in the Results section. The lib_xc library processes
the kinetic and xc functionals in the same way, and in the following we shall
discuss directly the sum of these two terms, i.e. $E_{\text{KXC}}[n,\nabla
n]=E_{\text{KE}}[n,\nabla n]+E_{\text{xc}}[n,\nabla n]$. A description of the
manipulations leading from the lib_xc output to a form of the second
derivatives suitable for our purposes is given in Appendix A. Here we report
the final result (see eq. (53))

$\displaystyle\frac{\partial^{2}E_{\text{KXC}}}{\partial\tilde{n}_{\bm{G}_{\beta}}\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}=\frac{\Omega}{N_{D}}\sum_{\bm{r}}\exp[\mathrm{i}\bm{G}_{\beta}\cdot\bm{r}]\exp[-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}]$
(31)
$\displaystyle\times\biggl{[}\Xi(\bm{r})+\mathrm{i}(G^{l}_{\beta}-G^{l}_{\alpha})\Phi^{l}(\bm{r})+G^{l}_{\beta}G^{m}_{\alpha}\Psi^{lm}(\bm{r})\biggr{]},$

where $N_{D}$ is the number of points used to discretize the density in real
space, $l$ and $m$ identify Cartesian components of the three-dimensional
vectors and repetition of these indexes indicates a sum over them. The
functions $\Xi(\bm{r})$, $\Phi^{l}(\bm{r})$ and $\Psi^{lm}(\bm{r})$ are
defined in Appendix A. Here we note for future convenience that they are real
periodic functions in real space. The second derivatives with respect to
different combination of coefficients and complex conjugates of coefficients
that are needed in the Jacobian matrix can be obtained from the expression
above using the rules of complex conjugation. Note the Hessian matrix above is
dominated by its diagonal for which $\bm{G}_{\beta}-\bm{G}_{\alpha}=\bm{0}$
and all terms in the sum over $\bm{r}$ interfere constructively. The same is
true for the term
$\frac{\partial^{2}E_{\text{KXC}}}{\partial\tilde{n}^{\star}_{\bm{G}_{\beta}}\partial\tilde{n}_{\bm{G}_{\alpha}}}$.
However, the two terms
$\frac{\partial^{2}E_{\text{KXC}}}{\partial\tilde{n}_{\bm{G}_{\beta}}\partial\tilde{n}_{\bm{G}_{\alpha}}}$
and
$\frac{\partial^{2}E_{\text{KXC}}}{\partial\tilde{n}^{\star}_{\bm{G}_{\beta}}\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}$,
entering in the expressions of
$\frac{\partial\sigma^{\star}_{\bm{G}_{\alpha}}}{\partial\tilde{n}_{\bm{G}_{\beta}}}$
and
$\frac{\partial\sigma_{\bm{G}_{\alpha}}}{\partial\tilde{n}^{\star}_{\bm{G}_{\beta}}}$,
respectively, are smaller in magnitude since they all carry a phase with
wavevector $\pm(\bm{G}_{\beta}+\bm{G}_{\alpha})$.

The structure of the Hartree and KXC terms motivates our choice of the
$2N_{G}$ complex vectors defined in Eqs. (8) and (9), and to solve the
constraint in Fourier space. This choice, in fact, ensures that the Jacobian
matrix is diagonally dominant facilitating the convergence of the standard
SHAKE algorithm.

Substituting Eq. (30) and the second derivatives of the KXC term, Eq. (31), in
Eq. (III), we obtain

$\displaystyle\left[\mathbb{J}_{D}\right]_{{\bm{G}_{\alpha}},\bm{G}_{\alpha}}=$
$\displaystyle\left(\frac{4\pi\Omega}{|\bm{G}_{\alpha}|^{2}}\right)^{2}+\frac{4\pi\Omega}{|\bm{G}_{\alpha}|^{2}}\frac{\Omega}{N_{D}}\sum_{\bm{r}}\Bigl{[}\Xi_{\text{P}}(\bm{r})+G^{l}_{\alpha}G^{m}_{\alpha}\Psi_{\text{P}}^{lm}(\bm{r})\Bigr{]}+\frac{4\pi\Omega}{|\bm{G}_{\alpha}|^{2}}\frac{\Omega}{N_{D}}\sum_{\bm{r}}\Bigl{[}\Xi_{t}(\bm{r})+G^{l}_{\alpha}G^{m}_{\alpha}\Psi_{t}^{lm}(\bm{r})\Bigr{]}$
(32)
$\displaystyle+\biggl{(}\frac{\Omega}{N_{D}}\biggr{)}^{2}\sum_{\bm{r},\bm{r}^{\prime}}e^{-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}}e^{\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}^{\prime}}\sum^{N_{G}}_{\lambda=1}e^{\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}}e^{-\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}^{\prime}}$
$\displaystyle\qquad\times\Bigl{[}\Xi_{\text{P}}(\bm{r})+\mathrm{i}(G^{l}_{\lambda}-G^{l}_{\alpha})\Phi^{l}_{\text{P}}(\bm{r})+G^{l}_{\lambda}G^{m}_{\alpha}\Psi_{\text{P}}^{lm}(\bm{r})\Bigr{]}\Bigl{[}\Xi_{t}(\bm{r}^{\prime})-\mathrm{i}(G^{l}_{\lambda}-G^{l}_{\alpha})\Phi_{t}^{l}(\bm{r}^{\prime})+G^{l}_{\lambda}G^{m}_{\alpha}\Psi_{t}^{lm}(\bm{r}^{\prime})\Bigr{]}$
$\displaystyle+\biggl{(}\frac{\Omega}{N_{D}}\biggr{)}^{2}\sum_{\bm{r},\bm{r}^{\prime}}e^{-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}}e^{\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}^{\prime}}\sum^{N_{G}}_{\lambda=1}e^{-\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}}e^{\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}^{\prime}}$
$\displaystyle\qquad\times\Bigl{[}\Xi_{\text{P}}(\bm{r})-\mathrm{i}(G^{l}_{\lambda}+G^{l}_{\alpha})\Phi_{\text{P}}^{l}(\bm{r})-G^{l}_{\lambda}G^{m}_{\alpha}\Psi_{\text{P}}^{lm}(\bm{r})\Bigr{]}\times\Bigl{[}\Xi_{t}(\bm{r}^{\prime})+\mathrm{i}(G^{l}_{\lambda}+G^{l}_{\alpha})\Phi_{t}^{l}(\bm{r}^{\prime})-G^{l}_{\lambda}G^{m}_{\alpha}\Psi_{t}^{lm}(\bm{r}^{\prime})\Bigr{]},$

where the first term in the RHS arises from the Hartree contribution to the
energy functional, while the others stem from the kinetic and xc functionals
or mixed terms. The Hartree contribution and the mixed terms can be easily
seen to scale linearly with respect to the number of grid points. The
difficult contributions in terms of numerical complexity are the last two
sums, which we now discuss.

The three sums appearing in these kinetic and xc terms indicate that the brute
force calculation of these contributions scales cubically with the number of
grid points in direct and reciprocal space, a prohibitive bottleneck for
realistic calculations. In Appendix C, we show how this cost can be reduced,
for GGA functionals, to $\Omega\log\Omega$ scaling by exploiting periodicity
in direct space. This result will be used for comparison purposes in the
Results section.

###  III.2 SHAKE with linear conjugate gradient (LCG-SHAKE)

An alternative way of performing the Newton-Raphson step, Eq. (22), can be
derived by building on an approach first suggested by Weinbach and
Elber.weinbach:2005  These authors proposed to obtain the value of the
constraints at each iteration $\kappa$ by considering an approximate and,
importantly, symmetric form of the Jacobian matrix which enables to solve the
Newton-Raphson system using linear conjugate gradient. This algorithm was
conceived to facilitate parallelization of classical constrained MD and
applied to the simulation of biological systems.weinbach:2005  In the
following we show how the approach can be adapted to the solution of the
Newton-Raphson step for MaZe OFDFT and demonstrate that it overcomes the
limitations of standard SHAKE by preserving linear scaling for general forms
of the energy functional, including those that employ nonlocal kinetic
energies.

Let us consider again Eq. (22): defining
$\delta\tilde{\bm{\gamma}}=\tilde{\bm{\gamma}}^{\kappa+1}-\tilde{\bm{\gamma}}^{\kappa}$,
this equation can be trivially rewritten as the system of linear equations

$\mathbb{J}(\tilde{\bm{n}}^{\kappa})\delta\tilde{\bm{\gamma}}=-\omega\bm{\sigma}\bigl{(}\bm{R}(t+\delta
t),\tilde{\bm{n}}^{\kappa}\bigr{)}.$ (33)

Our goal is to solve this system as efficiently as possible using linear
conjugate gradient. In analogy with Ref. 26, that considered the case of a
system of real equations, we observe that linear conjugate gradient cannot be
immediately applied because the Jacobian matrix is not hermitian (symmetric in
Ref. 26) due to the fact that the explicit expressions for its matrix elements
(see for example Eq. (III)) contain derivatives of the constraints evaluated
with the density computed at time $t$ and at the provisional value
$\tilde{\bm{n}}^{\kappa}$. Similar to Ref. 26, we then approximate the matrix
elements of the Jacobian as

$\displaystyle\\{\mathbb{J}(t)\\}_{\alpha,\beta}$
$\displaystyle=\frac{\partial\sigma_{\bm{G}_{\alpha}}}{\partial\tilde{\gamma}_{\bm{G}_{\beta}}}$
(34)
$\displaystyle\approx\sum^{N_{G}}_{\lambda=1}\biggl{[}\frac{\partial\sigma_{\bm{G}_{\alpha}}(t)}{\partial\tilde{n}_{\bm{G}_{\lambda}}}\frac{\partial\sigma^{\star}_{\bm{G}_{\beta}}(t)}{\partial\tilde{n}^{\star}_{\bm{G}_{\lambda}}}+\frac{\partial\sigma_{\bm{G}_{\alpha}}(t)}{\partial\tilde{n}^{\star}_{\bm{G}_{\lambda}}}\frac{\partial\sigma^{\star}_{\bm{G}_{\beta}}(t)}{\partial\tilde{n}_{\bm{G}_{\lambda}}}\biggr{]}$
$\displaystyle=\left\\{\Sigma(t)\cdot\Sigma(t)^{\dagger}\right\\}_{\alpha,\beta},$

and similarly for the other blocks of the matrix. Compared to the exact SHAKE
expression, all derivatives in the equations above are evaluated using the
density at time $t$, with similar expressions for the other blocks in Eq.
(23). We denote this approximate Jacobian as
$\mathbb{J}\bigl{(}\tilde{\bm{n}}(t)\bigr{)}\equiv\mathbb{J}(t)$. Direct
inspection of the approximate Jacobian shows that it is hermitian.
Furthermore, we have

${\bm{\zeta}}^{\star}\mathbb{J}(t){\bm{\zeta}}={\bm{\zeta}}^{\star}\Sigma(t)\Sigma^{\dagger}(t){\bm{\zeta}}>0,$
(35)

for all ${\bm{\zeta}}\neq\bm{0}$, proving that the matrix is positive
definite. The solution of

$\mathbb{J}(t)\delta\tilde{\bm{\gamma}}=-\omega\bm{\sigma}\bigl{(}\bm{R}(t+\delta
t),\tilde{\bm{n}}^{\kappa}\bigr{)},$ (36)

(that we take as an approximant for Eq. (33)) can then be obtained by direct
application of linear conjugate gradient with complex variables, with or
without preconditioning.joly:1993  (Note that using linear conjugate gradient,
requires the application of the matrix $\mathbb{J}(t)$ on complex vectors of
size $2N_{G}$. This operation can be performed with $\Omega\log\Omega$ scaling
for either GGA (see Appendix A) and nonlocal (see Appendix B) energy
functionals.) In Ref. 26 the authors showed that the solution of this system
converges to the same set of Lagrange multipliers of the standard SHAKE. In
the Results section, we also show that this method is more efficient than
standard SHAKE when applied to OFDFT dynamics.

##  IV Benchmarks and Results

The algorithms described in the previous sections are tested and compared in
the orbital-free DFT framework on constant energy (NVE) and constant
temperature (NVT) simulations of a system of particles of sodium in different
thermodynamic states and with different kinetic energy functionals. We start
by detailing the relative efficiency of the two methods and then, having shown
its merits, we move to compute static and dynamical statistical properties of
the system with the LCG-SHAKE. In particular, the reliability of the algorithm
for NVT calculations (never tested before), is demonstrated. Finally, we
discuss preliminary comparisons the performance of MaZe to that of standard
nonlinear conjugate gradient (NLCG) based Born-Oppenheimer dynamics by
considering energy conservation and time reversibility for a realistic test
system.

For all calculations discussed in this section, periodic boundary conditions
are enforced in all directions and the electrostatic interactions are treated
via the Ewald summation method with a cutoff radius of
$R_{\text{c}}=1.06\mathrm{\,\mbox{\AA}}$. For these tests, Slater dirac:1930 ;
bloch:1929  and Perdew & Zunger perdew:1981  functionals are used for the
exchange and correlation contributions $E_{\text{xc}}$ to the total energy,
respectively. Functionals adopted for the kinetic energy are detailed in the
following. The interaction of the atomic nucleus with the electrons up to the
last closed shell is modeled using the pseudopotential for sodium described by
Topp and Hopfield. topp:1973  The cutoff for the plane-wave expansion is set
to $G_{\text{c}}=10\mathrm{\,\mbox{$E_{\mathrm{h}}$}}$ and the timestep is set
to $\delta t=1\mathrm{\,fs}$ for all simulations. The NVE evolution of the
physical degrees of freedom is performed using the velocity Verlet algorithm.
As for the evolution of the coefficients of the density, the provisional
values used as the initial guess for the Newton–Raphson iterative procedure in
SHAKE is computed through the Verlet algorithm. The first two values of the
coefficients — necessary for the initialization of the Verlet algorithm — are
obtained through two minimizations of the energy functional using a NLCG
algorithm and the full MaZe dynamics for the system is started from the second
simulation step. The NVT ensemble is sampled via the use of a Langevin
thermostat where the evolution of the physical degrees of freedom is
integrated using the algorithm proposed by Vanden-Eijnden and Ciccotti.vanden-
eijnden:2006  The same approach is also used in the NVE runs — with the
friction parameter of the thermostat set to $\nu=10^{-5}\mathrm{\,au}$ — to
equilibrate the systems to the target temperatures (see below), via runs of
duration equal to $10\mathrm{\,ps}$. The number of electrons is kept constant
along the simulation by excluding the $\bm{G}=(0,0,0)$ component of the
electronic density expansion in reciprocal space from the dynamical system as
detailed in the previous sections. Tests are focused on two different
thermodynamic states of the system,222The melting temperature of Sodium at
room pressure is $T_{m}=371\mathrm{\,K}$. solid ($T_{s}\approx
25\mathrm{\,K}$, $\rho_{s}=969.5\mathrm{\,kg}\mathrm{\,m^{-3}}$) and liquid
($T_{l}\approx 430\mathrm{\,K}$,
$\rho_{l}=926.8\mathrm{\,kg}\mathrm{\,m^{-3}}$).

To set the stage for the production runs, we recall that, as discussed in Sec.
III.1, when using the standard SHAKE algorithm, the calculation of the matrix
$\mathbb{J}_{D}$ defined in Eq. (27) is the most expensive ingredient in the
solution of the constraints. To explore this calculation, we implemented and
computed it, as per its $\Omega\log\Omega$ scaling expression in Eq. (70), for
a system of $N=128$ sodium atoms in the solid (box size of
$L_{s}=17.1\mathrm{\,\mbox{\AA}}$) and liquid (box size
$L_{l}=17.4\mathrm{\,\mbox{\AA}}$). We considered two different kinetic energy
functionals: Thomas-Fermi-von Weizsäckervonweizsacker:1935  (TFvW) and
Perrot.perrot:1994  Two simulations of 1000 steps for the solid and liquid
systems were performed, computing $\mathbb{J}_{D}$ for configurations obtained
every 100 steps. Results for liquid Sodium are reported as a function of the
modulus of the wave-vector, as red (TFvW) and blue (Perrot) squares in Figure
1. Calculations for the solid (not reported) show very similar behavior. Note
that the results for the 10 configurations selected in these test runs are
essentially superimposed for both curves. These plots show then that the shape
of the denominator does not change significantly along the simulation. This
behavior is an additional confirmation of the fact that the Hessian of the
Hartree term in the energy functional — independent of the electronic density
and thus constant along the simulation — has a dominant role in the
computation of the denominator, therefore corroborating the discussion after
Eq (30) and Eq. (31).

In Figure 1, the numerical values of the SHAKE denominator are also compared
to the analytical approximation:

$\left[\mathbb{J}_{D}\right]^{\text{analytical}}_{{\bm{G}_{\alpha}},\bm{G}_{\alpha}}=\biggl{(}\frac{\Omega}{4\rho}\biggr{)}^{2}|\bm{G}_{\alpha}|^{4}+\bigl{(}4\pi\Omega\bigr{)}^{2}|\bm{G}_{\alpha}|^{-4},$
(37)

where $\Omega$ is the volume of the simulation box and
$\rho=N_{\text{el}}/\Omega$ is the uniform electron density of the system.
This analytical approximation is obtained as the sum of the diagonal Hartree
contribution, which is dominant at small $|\bm{G}_{\alpha}|$, and of the von
Weizsäcker contribution to the kinetic energy functional, dominant at large
$|\bm{G}_{\alpha}|$. The latter is given by

$E_{\mathrm{vW}}[n]=\frac{1}{8}\int_{\Omega}\mathrm{d}^{3}\bm{r}\frac{|\nabla
n(\bm{r})|^{2}}{n(\bm{r})}.$ (38)

Assuming $n(\bm{r})=\rho$ constant at the denominator in the integral, we have

$\displaystyle
E_{\mathrm{vW}}[n]\approx\frac{1}{8\rho}\int_{\Omega}\mathrm{d}^{3}\bm{r}|\nabla
n(\bm{r})|^{2}\approx\frac{\Omega}{4\rho}\sum_{\alpha=1}^{N_{G}}|\bm{G}_{\alpha}|^{2}\tilde{n}^{\star}_{\bm{G}_{\alpha}}\tilde{n}_{\bm{G}_{\alpha}}.$
(39)

Figure 1: Main plot: Diagonal of the SHAKE matrix (SHAKE denominator) for
liquid sodium for different kinetic energy functionals (red points TFvW, blue
points Perrot) compared with the analytic approximation obtained from Eq. (37)
(yellow curve). Inset: Ratio between the SHAKE denominator and the analytic
approximation.

In Figure 1 the analytical approximation, Eq. (37), is shown as the yellow
curve and turns out to be remarkably accurate, at least for the system under
investigation. (Similar results were obtained for solid Sodium, not shown.)
Adoption of this expression for the SHAKE denominator might be used to reduce
the computational cost of our calculations. In the following, however, we
employ it only as an effective preconditioner in the LCG-SHAKE algorithm
described in Section III.2.

We now move to the comparison of the efficiency of the standard SHAKE and LCG-
SHAKE implementations of MaZe. The tests below have been performed using the
same set of configurations employed in the discussion above. We begin by
showing, in Figure 2, the path to convergence of the Newton-Raphson solution
of the constraint equations for the two different algorithms. The figure
reports, as a function of the total number of iterations needed to converge,
the modulus of the largest component of the constraint vector below a given
threshold. In red, results for the standard SHAKE are reported, showing the
characteristic monotonic decrease of
$\mathrm{max}_{\bm{G}_{\alpha}}|\sigma_{\bm{G}_{\alpha}}|$ with number of
iterations already observed in previous MaZe implementations.bonella:2020  The
blue curves, on the other hand, show the path to convergence of the LCG-SHAKE.
We recall that for the latter, the linear system, Eq. (33), needs to be solved
via linear conjugate gradient, implying an additional iterative procedure.
This linear system is considered solved when the square-root of the norm of
the residuals is less than $5\cdot 10^{-14}\mathrm{\,au}$. The iterations
needed for this solution are represented in Figure 2 as the horizontal
segments. Both for standard SHAKE and for the LCG-SHAKE, we set the
convergence threshold of the Newton-Raphson algorithm to
$10^{-10}\mathrm{\,au}$. A successive over-relaxation parameter $\omega=0.35$
and $\omega=0.21$ is set in the standard SHAKE algorithm when the TFvW and
Perrot kinetic energy functionals are used, respectively. For the LCG-SHAKE
algorithm, $\omega=1$ proved the best choice in all tests performed. The
panels in the figure show that, for the functional considered and irrespective
of the solid (upper panel) or liquid (lower) state of the system, both
algorithms converge in a relatively small number of iterations. Furthermore,
LCG-SHAKE clearly outperforms standard SHAKE when the Perrot functional is
adopted. The apparent advantage of standard SHAKE for the TFvW functional is
however limited to the number of iterations. In fact, comparing the efficiency
of the two algorithms in terms of average wall time to solution, Figure 3, we
see that the algorithms (both very efficient) have relatively similar
performances for the TFvW calculations, while LCG-SHAKE is faster than
standard SHAKE for Perrot, most notably in the case of the liquid simulations.

Figure 2: Total number of iterations to reach convergence for Newton-Raphson
solution of the equation of constraints via standard SHAKE and LCG-SHAKE
algorithms. In the LCG-SHAKE case the horizontal (parallel to the $x$-axis)
segments represent the number of linear CG iterations used to solve the
linearized equation of the constraints, as the values on the $y$-axis refers
to the Newton-Raphson iterations. Figure 3: Average wall time per timestep
comparison between Standard SHAKE and LCG-SHAKE algorithms for different
kinetic energy functionals and different thermodynamic states.

We conclude this survey of the numerical properties of MaZe by demonstrating,
for the more efficient LCG-SHAKE, the linear scaling of the computational cost
with respect to the system size. For this particular benchmark, 10 simulation
steps of a system of solid sodium modelled with the TFvW kinetic energy
functional have been performed for a number of particles ranging from $N=2000$
to $N=500094$, and the average wall time per timestep has been reported in
Figure 4 as a function of the number of particles in the simulation, showing a
very good linear scaling.

Figure 4: Performance scaling of a system of solid sodium with TFvW kinetic
energy functional with respect to the number of particles for the LCG-SHAKE
algorithm.

Accuracy and reliability of the LCG-SHAKE MaZe method in performing physically
realistic simulations are now shown on a system of liquid sodium modelled
through the Perrot kinetic energy functional. Two sets of simulations are
presented in the following, in the NVE and NVT ensemble. The temperature of
the system is the same for both sets of simulations, $T=434\mathrm{\,K}$. In
all cases, the calculation of the coefficients of the density in reciprocal
space is performed using MaZe together with the LCG-SHAKE algorithm, where the
provisional value of the additional dofs is computed through a standard Verlet
algorithm.

We begin by considering the radial distribution function of the system, shown
in Figure 5. The red curve shows the result obtained for an NVE run of
duration $100\mathrm{\,ps}$. Position and height of the two peaks are in
agreement with previous numerical studies qian:1990  performed at the same
temperature and are compatible with experimental results greenfield:1971 ;
lee:1978  obtained at $T_{\text{exp}}=373\mathrm{\,K}$. The green and purple
curves, on the other hand, report NVT results for two different values of the
friction coefficient $\nu$ obtained from $50\mathrm{\,ps}$ long simulations.
As expected, within noise, this static property is independent of the choice
of the friction and identical to the NVE output. Figure 6 shows the velocity
autocorrelation function obtained from the constant energy run. Both the
characteristics of this curve and the diffusion coefficient obtained
integrating the result over time (shown as the red curve in Figure 7) are in
very good agreement with previous numerical results.shimojo:1994  In
particular, our computed value of the NVE diffusion is equal to $5.37\pm
0.16\mathrm{\,cm^{2}}\mathrm{\,s^{-1}}$ which agrees with the numerical
results reported in Ref. 51 and with the experimental results reported in Ref.
55 for the temperature of $T=434\mathrm{\,K}$ at which the simulations are
performed. Figure 7 shows also the value of the diffusion coefficients
computed in NVT simulations with different values of the friction
coefficients. In this case, the result naturally depends on the choice of
$\nu$, but we observe that, again as expected, it approaches the reference NVE
for small friction further illustrating the reliability of our canonical MaZe.

Figure 5: Comparison between radial distribution functions from NVT runs
performed using the algorithm in Ref. 50 for a Langevin thermostat for
different values of the friction parameters $\nu$. The NVE result is obtained
from a 100$\mathrm{\,ps}$ simulation of $N=128$ particles of sodium at
$T=434$$\mathrm{\,K}$. The other values of $\nu$ are obtained from
50$\mathrm{\,ps}$ NVT simulations of $N=128$ particles of sodium at
$T^{*}=434$$\mathrm{\,K}$. Figure 6: Normalized velocity autocorrelation
function for a 100$\mathrm{\,ps}$ NVE simulation of $N=128$ particles of
liquid sodium at $T=434$$\mathrm{\,K}$. Figure 7: Comparison between diffusion
coefficients from NVT runs for different values of the friction parameter
$\nu$. The NVE result is obtained from a 100$\mathrm{\,ps}$ simulation of
$N=128$ particles of Na at $T=434$$\mathrm{\,K}$. The other values of $\nu$
are obtained from 50$\mathrm{\,ps}$ NVT simulations of $N=128$ particles of Na
at $T^{*}=434$$\mathrm{\,K}$.

In the NVT ensemble, the number of iterations to converge the calculation of
the electronic density is also affected by the choice of $\nu$. In Figure 8
the path to convergence for a typical run is shown for two different values of
the friction coefficient: $\nu=10^{-3}$$\mathrm{\,au}$, in purple, and
$\nu=10^{-5}$$\mathrm{\,au}$, in green, and compared with the NVE result, in
red.

Figure 8: Total number of iterations to reach convergence for NVT simulation
with LCG-SHAKE for different values of the friction parameter $\nu$. The
horizontal (parallel to the $x$-axis) segments represent the number of linear
CG iterations used to solve the linearized equation of the constraints, as the
values on the $y$-axis refers to the Newton-Raphson iterations. Results are
for a typical minimization path along a simulation.

The trends in Figure 8 indicate that the increase in the number of iterations
necessary to solve the equation of the constraints with increasing $\nu$ is
mainly due to a difference in the value of the maximum component of the
constraints at the start of the iterative loop. Since the thermostat only
affects the dynamics of the physical degrees of freedom, this depends on the
fact that the provisional value of the electronic coefficients used to
initialize the SHAKE routine at each timestep, is a less accurate first guess
for higher values of the ionic friction. Indeed, once the initial penalty
associated to the violation of the constraint at zero iteration is absorbed,
typically in one Newton-Raphson step, the speed of convergence becomes
comparable with that of NVE. Furthermore, even for large values of $\nu$, the
total number of iterations to convergence remains relatively small.

We conclude this section by benchmarking LCG-SHAKE MaZe against standard Born-
Oppenheimer dynamics, based on nonlinear conjugate-gradient (NLCG)
minimization.press:1992-book  Tests are particularly focused on the time-
reversibility properties of the algorithms and relative efficiency of the two
approaches. Note however that for Born-Oppenheimer dynamics we have used a
simple textbook implementation of the NLCG algorithm. The initial density for
NLCG is taken at each timestep to be the final density of the last timestep,
thus no extrapolation was performed. Similarly, no preconditioning of the NLCG
was used. The comparison is performed on a set of simulations of $N=16$
particles of liquid sodium ($L=8.68\mathrm{\,\mbox{\AA}}$) using the Perrot
kinetic energy functional.pearson:1993 ; perrot:1994  The low number of
particles has been chosen to ensure that possible numerical errors arising in
the simulations are mainly due to the propagation algorithm and not to the
computation of the forces. Time reversibility of the algorithms is tested as
follows: Starting from the same initial configuration, a
$\mathcal{T}=10\mathrm{\,ps}$ simulation is performed using either nonlinear
conjugate-gradient Born-Oppenheimer or MaZe to enforce the minimum energy
condition (forward trajectory). The velocities of the final configuration of
the forward runs are then reversed and a new 10$\mathrm{\,ps}$ simulation is
performed (backward trajectory). Along the forward and backward trajectories,
positions of the particles are recorded every 10$\mathrm{\,fs}$. Average root
mean square displacements normalized to the length of the simulation box

$\mathrm{RMSD}(t)=\frac{1}{NL}\sum_{i=1}^{N}\sqrt{\bigl{(}\bm{r}_{i}^{+}(\mathcal{T}-t)-\bm{r}_{i}^{-}(t)\bigr{)}^{2}},$
(40)

are computed for corresponding configurations along backward and forward
trajectories for the two algorithms. The convergence threshold of the minimum
condition is set to a very stringent value (residual lower than
$10^{-10}$$\mathrm{\,au}$ for the modulus of the maximum component of the
derivative of the energy, compare with Tables 1 and 2).

Figure 9: Test of time-reversibility properties of MaZe (LCG-SHAKE, curves
with red symbols) with respect to standard Born-Oppenheimer alternatives,
where the minimization of the energy functional is performed with standard
nonlinear conjugate-gradient (NLCG, curves with blue symbols) techniques with
different convergence thresholds. The test is based on comparing
configurations along a forward- and a backward-in-time trajectory and
comparing the average RMSD for corresponding configurations (see Eq. (40) and
relative discussion). The upper limit on the $y$-axis represents an estimate
of the maximum allowed value of the RMSD considering the presence of the PBCs
(and therefore the fact that two particles are never more than
$\frac{L\sqrt{3}}{2}$ distant from each other). An intermediate threshold of
$10^{-5}$ for MaZe is also presented (red line open triangles) to demonstrate
the levels of time reversibility and stability attainable with intermediate
tolerances.

Results for LCG-SHAKE MaZe and NLCG Born-Oppenheimer with this threshold are
shown as the red and blue open squares symbols, respectively, in Figure 9. The
figure shows that results for the two algorithms are very similar with this
strict tolerance. The slow violation of time reversibility observed can be
attributed to the numerical error associated to the finite timestep and round
off that are amplified by the Lyapunov instability.

With this strong convergence requirement, however, the average time needed by
NLCG to solve the minimization problem is ten times larger than the one needed
by MaZe to solve the equation of the constraints. Table 1 reports timings for
LCG-SHAKE for different tolerances. For each user-chosen tolerance, we report
not only the observed timing but also the actual residual achieved. The final
residual can be significantly lower than the target tolerance as only one
Newton-Raphson (NR) iteration can lead to high convergence. The average number
of NR iterations is also reported in this table.

Table 1: Performance assessment for MaZe using LCG-SHAKE algorithm for the solution of the constraint equation. The user-chosen tolerance of the iterative procedure is reported in the second column, while in the third column the actual residual on the constraints after the self-consistent cycle is shown. Wall timings for a single MD step on a single core (AMD Ryzen 7 5800X CPU) and number of Newton-Raphson iterations to reach convergence are reported in the fourth and fifth column respectively. Method | Tol ($\mathrm{\,au}$) | Res ($\mathrm{\,au}$) | Time (s) | NR Iters
---|---|---|---|---
LCG-SHAKE | $10^{-1}$ | $1.4\times 10^{-4}$ | $3.3\times 10^{-2}$ | 1
$10^{-3}$ | $2.0\times 10^{-6}$ | $5.8\times 10^{-2}$ | 1
$10^{-5}$ | $1.3\times 10^{-6}$ | $8.0\times 10^{-2}$ | 1
$10^{-7}$ | $5.0\times 10^{-9}$ | $1.5\times 10^{-1}$ | 2
$10^{-10}$ | $2.2\times 10^{-11}$ | $2.7\times 10^{-1}$ | 3

For comparison, Table 2 reports the timings for our textbook implementation of
the NLCG method. In order to allow for a one to one correspondence with the
LCG-SHAKE case, the NLCG tolerance is set to the observed residual in Table 1.

Table 2: Performance assessment of standard Born-Oppenheimer MD using nonlinear conjugate gradient algorithm. Comparison between the methods is carried out by setting the same tolerance (second column) for the nonlinear conjugate gradient method as the residual on the constraints shown in Table 1. Wall timings for a single MD step for the same system used for tests reported in Table 1 are reported in the last column. An additional high tolerance is shown in the first line where the performances in terms of wall time are comparable between the two methods. Note that MaZe algorithm is faster than NLCG even when the precision is nine order of magnitudes higher than standard nonlinear conjugate gradient. Method | Tol ($\mathrm{\,au}$) | Time (s)
---|---|---
NLCG | $8.1\times 10^{-2}$ | $3.8\times 10^{-1}$
$8.8\times 10^{-5}$ | $9.9\times 10^{-1}$
$1.3\times 10^{-6}$ | 2.3
$4.4\times 10^{-9}$ | 3.7
$1.6\times 10^{-11}$ | 4.6

To obtain comparable average times per timestep along the simulation, the
tolerance of the nonlinear conjugate gradient has to be relaxed by nine orders
of magnitude compared the lowest tolerance achieved with LCG-SHAKE. NLCG runs
at low convergence, however, lead to dramatic effects and, in particular, the
algorithm suffers an essentially complete loss of time reversibility as shown
in Figure 9, where the blue open circles shows that the average distance
between time-reversed configurations quickly reaches the box dimensions. We
also report in Figure 9, the RMSD for an LCG-SHAKE run with intermediate
tolerance, corresponding to only one NR iteration per timestep. It can be seen
that such a trajectory is stable, with an RMSD that parallels the RMSD with
the lowest tolerance.

Essentially identical results are obtained for energy conservation in the
different runs. Figure 10 shows the MaZe (red curve with open squares) and
NLCG (blue curve with open squares) energies obtained with the same timestep.
Also in this case, when the same stringent convergence threshold is imposed,
the results for the two algorithms are superimposed. However, a significant
drift in the ionic energy appears as shown (blue curve with open circles) when
the convergence in NLCG is relaxed to match the timings of MaZe.

Figure 10: Test of energy-conservation properties of MaZe (LCG-SHAKE, curves
with red symbols) with respect to standard Born-Oppenheimer alternatives,
where the minimization of the energy functional is performed with standard
nonlinear conjugate gradient (NLCG, curves with blue symbols) techniques with
different convergence thresholds. An intermediate threshold of $10^{-5}$ for
MaZe is also presented (red line open triangles) to prove the stability and
energy conservation of the algorithm for intermediate tolerances. The total
energy on the $y$-axis is shifted by the value at the beginning of the
simulation. Results are for the forward trajectories only.

##  V Conclusions

This work generalizes previous results for OFDFT first-principles dynamics
using the MaZe approach in two important ways. Firstly, a detailed derivation
of the evolution equations has been presented, highlighting, in particular,
the consequences of the complex nature of the Lagrange multipliers associated
with the evolution of the coefficients of the expansion of the electronic
density in the plane-wave basis. The OFDFT MaZe dynamical system can be
integrated combining standard MD algorithms for the evolution of the ions with
the SHAKE method to enforce the minimum condition on the electronic degrees of
freedom. Efficient implementation of the constraints for generic forms of the
kinetic and exchange-correlation functionals, however, is highly nontrivial.
We have shown how to achieve linear scaling in the standard implementation of
SHAKE for GGA functionals, and introduced an original approach that enables
$\Omega\log\Omega$ scaling of the method for arbitrary forms of the kinetic
and exchange-correlation functional, including nonlocal expressions where
standard SHAKE fails. The accuracy and efficiency of both algorithms has been
demonstrated via simulations of a challenging benchmark system, solid and
liquid sodium. Use of the MaZe approach beyond constant energy simulations has
also been illustrated, presenting results for NVT simulations of the same
system. An in-depth comparison of OFDFT MaZe with alternative state-of-the-art
approaches, in particular XL-BOMD and advanced implementations of NLCG Born-
Oppenheimer dynamics including preconditioning and smart initialization, will
be performed in future studies. Here we make some general preliminary
comments. Numerically, when compared to the original Car-Parrinello dynamics,
MaZe has the advantage of enabling the use of a timestep consistent with the
ion dynamics instead of the shorter timesteps typically used in Car-Parrinello
MD. It also enforces full adiabatic separation via the limit of zero mass for
the auxiliary variables. Together with the fact that SHAKE efficiently
converges the minimum condition essentially to numerical precision, this
prevents energy drifts and avoids the need for thermostats to stabilize the
dynamics. Even though the overall approach is different, MaZe shares some
similarities the recent XL-BOMD method, such as the presence of a limit of the
fictitious mass going to zero. However, in XL-BOMD an extra dissipation of the
dynamics of the electronic dofs is needed in practice and the accuracy of
propagating the electronic dofs still remains dependent on the integration
timestep $\delta t$ in XL-BOMD, going as $\delta t^{-2}$.niklasson:2014  MaZe
also avoids the numerical cost of computing the kernel involved in XL-BOMD.
The advantage over Car-Parrinello and XL-BOMD comes at the non-trivial cost of
the computation of the Lagrange multipliers. The calculations presented in
this work, however, show that this numerical cost can be controlled and
remains quasilinear with respect to the system size. The developments
presented in this paper, and the remarkable efficiency of the algorithms
introduced, then open the way to applications of MaZe to more realistic OFDFT
models of condensed-phase systems and set the stage to the future extension of
the approach to Kohn-Sham based DFT.

## Acknowledgments

The authors would like to thank Giovanni Ciccotti, Benjamin Rotenberg and
Mathieu Salanne for fruitful discussions. Part of this work was supported
under the European Union’s Horizon 2020 research and innovation programme
under Grant Agreement No. 676531 (project E-CAM). R.V. would also like to
acknowledge hospitality at CECAM, where this work started.

## Data Availability Statement

The data that support the findings of this study together with the source code
used for the numerical simulations are available from the corresponding author
upon reasonable request.

##  Appendix A Expressions of forces and Hessians and $\Omega\log\Omega$
scaling for matrix-vector multiplications of standard GGA functionals

The calculations in this paper require the value of the energy functional and
its first and second derivatives. The manipulations presented in the main text
are for the most part general but in our specific implementation we exploit
the lib_xc library to obtain these quantities. lib_xc takes as input the
electronic density and its derivatives in real space and returns the value of
$\Epsilon$, the energy per unit particle, and its derivatives in real space.
The library employs the electronic density $n(\bm{r})$ and the modulus square
of its gradient $\varsigma(\bm{r})=|\nabla n(\bm{r})|^{2}$ as independent
“variables”. This is the general formulation for a GGA functional.

The Hartree contribution to first and second derivatives is trivial, the most
problematic term being the kinetic-exchange-correlation contribution that we
adress first. The kinetic-exchange-correlation energy is returned as

$E_{\text{KXC}}[n(\bm{r}),\varsigma(\bm{r})]=\frac{\Omega}{N_{\text{D}}}\sum_{\bm{r}}\varepsilon_{\text{KXC}}\bigl{(}n(\bm{r}),\varsigma(\bm{r})\bigr{)},$
(41)

where $N_{D}$ is the number of points in the direct space grid used to
represent the discretized density and
$\varepsilon_{\text{KXC}}\bigl{(}n(\bm{r}),\varsigma(\bm{r})\bigr{)}$ the
local kinetic-exchange-energy per unit volume. The library also outputs the
derivatives of the functional with respect to the density and to the modulus
square of the gradients:

$\displaystyle\upsilon_{\text{KXC}}^{n}\bigl{(}n(\bm{r}),\varsigma(\bm{r})\bigr{)}$
$\displaystyle=\frac{\partial\varepsilon_{\text{KXC}}\bigl{(}n(\bm{r}),\varsigma(\bm{r})\bigr{)}}{\partial
n(\bm{r})};$ (42)
$\displaystyle\upsilon_{\text{KXC}}^{\varsigma}\bigl{(}n(\bm{r}),\varsigma(\bm{r})\bigr{)}$
$\displaystyle=\frac{\partial\varepsilon_{\text{KXC}}\bigl{(}n(\bm{r}),\varsigma(\bm{r})\bigr{)}}{\partial\varsigma(\bm{r})}.$

Finally, lib_xc provides also the following derivatives, associated to the
Hessian of the functional

$\displaystyle f_{\text{KXC}}^{nn}\bigl{(}n(\bm{r}),\varsigma(\bm{r})\bigr{)}$
$\displaystyle=\frac{\partial^{2}\varepsilon_{\text{KXC}}\bigl{(}n(\bm{r}),\varsigma(\bm{r})\bigr{)}}{\partial^{2}n(\bm{r})};$
(43) $\displaystyle
f_{\text{KXC}}^{n\varsigma}\bigl{(}n(\bm{r}),\varsigma(\bm{r})\bigr{)}$
$\displaystyle=\frac{\partial^{2}\varepsilon_{\text{KXC}}\bigl{(}n(\bm{r}),\varsigma(\bm{r})\bigr{)}}{\partial
n(\bm{r})\partial\varsigma(\bm{r})};$ $\displaystyle
f_{\text{KXC}}^{\varsigma\varsigma}\bigl{(}n(\bm{r}),\varsigma(\bm{r})\bigr{)}$
$\displaystyle=\frac{\partial^{2}\varepsilon_{\text{KXC}}\bigl{(}n(\bm{r}),\varsigma(\bm{r})\bigr{)}}{\partial^{2}\varsigma(\bm{r})}.$

All derivatives above are defined per unit volume. The lib_xc outputs can be
easily combined to obtain the first and second derivatives of the kinetic and
exchange correlation energy functional with respect to the Fourier
coefficients of the electronic densities (see sections III.1 and III.2). As an
example, we describe in some detail the calculation of the first derivative of
the kinetic and exchange functional. We have:

$\displaystyle\frac{\partial
E_{\text{KXC}}}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}$
$\displaystyle=\frac{\Omega}{N_{D}}\sum_{\bm{r}}\frac{\partial\varepsilon_{\text{KXC}}(\bm{r})}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}$
(44)
$\displaystyle=\frac{\Omega}{N_{D}}\sum_{\bm{r}}\left[\frac{\partial\varepsilon_{\text{KXC}}(\bm{r})}{\partial
n(\bm{r})}\frac{\partial
n(\bm{r})}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}+\frac{\partial\varepsilon_{\text{KXC}}(\bm{r})}{\partial\varsigma(\bm{r})}\frac{\partial\varsigma(\bm{r})}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}\right]$
$\displaystyle=\frac{\Omega}{N_{D}}\sum_{\bm{r}}\left[\frac{\partial
n(\bm{r})}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}\upsilon_{\text{KXC}}^{n}(\bm{r})+\frac{\partial\varsigma(\bm{r})}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}\upsilon_{\text{KXC}}^{\varsigma}(\bm{r})\right],$

where the derivatives of the kinetic-exchange-correlation kernel have been
defined in Eqs. (42). The derivatives of $n(\bm{r})$ and $\varsigma(\bm{r})$
with respect to the coefficients of the electronic density in reciprocal space
can be easily computed. From the definition of the Fourier transform of the
electronic density, we have

$\frac{\partial
n(\bm{r})}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}=\exp[-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}].$
(45)

The derivative of $\varsigma(\bm{r})$ can be computed using the chain rule

$\frac{\partial\varsigma(\bm{r})}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}=\frac{\partial\varsigma(\bm{r})}{\partial\bm{\nabla}n(\bm{r})}\cdot\frac{\partial\bm{\nabla}n(\bm{r})}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}.$
(46)

Above, the gradient of the density is given by

$\displaystyle\bm{\nabla}n(\bm{r})$
$\displaystyle=\sum^{N_{G}}_{\beta=1}\Bigl{[}\tilde{n}_{\bm{G}_{\beta}}\exp[\mathrm{i}\bm{G}_{\beta}\cdot\bm{r}]\mathrm{i}\bm{G}_{\beta}-\tilde{n}^{\star}_{\bm{G}_{\beta}}\exp[-\mathrm{i}\bm{G}_{\beta}\cdot\bm{r}]\mathrm{i}\bm{G}_{\beta}\Bigr{]},$
(47)

taking the derivative of this quantity with respect to
$\tilde{n}^{\star}_{\bm{G}_{\alpha}}$, and recalling that
$\varsigma(\bm{r})=|\bm{\nabla}n(\bm{r})|^{2}$, Eq. (46) can then be recast as

$\frac{\partial\varsigma(\bm{r})}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}=-2\bm{\nabla}n(\bm{r})\cdot\mathrm{i}\bm{G}_{\alpha}\exp[-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}].$
(48)

Substituting Eqs. (45) and (48) in Eq. (44) we get

$\displaystyle\frac{\partial
E_{\text{KXC}}}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}=\frac{\Omega}{N_{D}}\sum_{\bm{r}}\upsilon_{\text{KXC}}^{n}(\bm{r})\exp[-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}]-\frac{\Omega}{N_{D}}\biggl{[}2\sum_{\bm{r}}\upsilon_{\text{KXC}}^{\varsigma}(\bm{r})\bm{\nabla}n(\bm{r})\exp[-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}]\biggr{]}\cdot\mathrm{i}\bm{G}_{\alpha}.$
(49)

The expression above can be conveniently rewritten as the sum of two Fourier
transforms:

$\displaystyle\frac{\partial
E_{\text{KXC}}}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}$
$\displaystyle=\frac{\Omega}{N_{D}}\biggl{[}\mathrm{FT}\Bigl{[}\upsilon_{\text{KXC}}^{n}(\bm{r})\Bigr{]}-2\mathrm{FT}\Bigl{[}\upsilon_{\text{KXC}}^{\varsigma}(\bm{r})\bm{\nabla}n(\bm{r})\Bigr{]}\cdot\mathrm{i}\bm{G}_{\alpha}\biggr{]},$
(50)

and computed numerically via discrete Fourier transforms
($\mathrm{FT}[\bullet]$ in the equation above) with an $\Omega\log\Omega$
cost.

The second derivatives of the kinetic-exchange-correlation functional are
obtained by differentiation of the result above and using, where appropriate,
the definitions in Eqs. (43) and the chain rule for the derivatives. For
example, we have

$\displaystyle\frac{\partial^{2}E_{\text{KXC}}}{\partial\tilde{n}_{\bm{G}_{\beta}}\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}$
$\displaystyle=\frac{\Omega}{N_{D}}\sum_{\bm{r}}\exp[\mathrm{i}\bm{G}_{\beta}\cdot\bm{r}]\exp[-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}]$
(51)
$\displaystyle\times\biggl{[}f^{nn}(\bm{r})-2\mathrm{i}\partial_{l}n(\bm{r})(-G_{\beta}^{l}+G_{\alpha}^{l})f^{n\varsigma}(\bm{r})+2G_{\beta}^{l}G_{\alpha}^{m}\bigl{(}\upsilon^{\varsigma}(\bm{r})\delta_{lm}+2\partial_{l}n(\bm{r})\partial_{m}n(\bm{r})f^{\varsigma\varsigma}(\bm{r})\bigr{)}\biggr{]}.$

In the equation above, $\partial_{l}n(\bm{r})$ indicates the derivative of the
density with respect to the Cartesian component $l$ of the position $\bm{r}$
and repeated Latin indexes are summed over. For future convenience, we define
the functions

$\displaystyle\Xi(\bm{r})$ $\displaystyle\equiv f^{nn}(\bm{r});$ (52)
$\displaystyle\Phi^{l}(\bm{r})$ $\displaystyle\equiv
2\partial_{l}n(\bm{r})f^{n\varsigma}(\bm{r});$
$\displaystyle\Psi^{lm}(\bm{r})$ $\displaystyle\equiv
2\bigl{(}\upsilon^{\varsigma}(\bm{r})\delta_{lm}+2\partial_{l}n(\bm{r})\partial_{m}n(\bm{r})f^{\varsigma\varsigma}(\bm{r})\bigr{)},$

to rewrite Eq. (51) as

$\displaystyle\frac{\partial^{2}E_{\text{KXC}}}{\partial\tilde{n}_{\bm{G}_{\beta}}\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}=\frac{\Omega}{N_{D}}\sum_{\bm{r}}\exp[\mathrm{i}\bm{G}_{\beta}\cdot\bm{r}]\exp[-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}]\times\biggl{[}\Xi(\bm{r})+\mathrm{i}(G_{\beta}^{l}-G_{\alpha}^{l})\Phi^{l}(\bm{r})+G_{\beta}^{l}G_{\alpha}^{m}\Psi^{lm}(\bm{r})\biggr{]}.$
(53)

The second derivatives with respect to different combination of coefficients
and complex conjugates of coefficients that are needed in the calculation of
the matrix
$\Sigma^{\dagger}=\Sigma^{\dagger}=\Sigma^{\dagger}_{\mathrm{Hart}}+\Sigma^{\dagger}_{\mathrm{KXC}}$
can be obtained from the expression above using the rules of complex
conjugation.

While the Hartree term is diagonal in the density coefficients when using
plane waves

$\left\\{\Sigma^{\dagger}_{\mathrm{Hart}}\right\\}_{\alpha,\beta}=\left\\{\Sigma^{\dagger}_{\mathrm{Hart}}\right\\}_{\alpha,N_{G}+\beta}=\frac{4\pi\Omega}{|\bm{G}_{\alpha}|^{2}}\delta_{\alpha,\beta},$
(54)

and so the matrix-vector product is trivial, the final result for the matrix
$\Sigma^{\dagger}_{\mathrm{KXC}}$ is given by

$\displaystyle\left\\{\Sigma^{\dagger}_{\mathrm{KXC}}\right\\}_{\alpha,\beta}$
$\displaystyle=\frac{\Omega}{N_{D}}\sum_{\bm{r}}\exp[\mathrm{i}\bm{G}_{\beta}\cdot\bm{r}]\exp[-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}]\times\biggl{[}\Xi(\bm{r})+\mathrm{i}(G^{l}_{\beta}-G^{l}_{\alpha})\Phi^{l}(\bm{r})+G^{l}_{\beta}G^{m}_{\alpha}\Psi^{lm}(\bm{r})\biggr{]};$
(55)
$\displaystyle\left\\{\Sigma^{\dagger}_{\mathrm{KXC}}\right\\}_{\alpha,N_{G}+\beta}$
$\displaystyle=\frac{\Omega}{N_{D}}\sum_{\bm{r}}\exp[-\mathrm{i}\bm{G}_{\beta}\cdot\bm{r}]\exp[-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}]\times\biggl{[}\Xi(\bm{r})-\mathrm{i}(G^{l}_{\beta}+G^{l}_{\alpha})\Phi^{l}(\bm{r}))-G^{l}_{\beta}G^{m}_{\alpha}\Psi^{lm}(\bm{r})\biggr{]}.$

which is more complicated to treat. Considering the three terms in the square
brackets separately and starting from the first one, the matrix-vector product
can however be written

$\displaystyle\frac{\Omega}{N_{D}}\sum_{\beta=1}^{N_{G}}\sum_{\bm{r}}\exp[\mathrm{i}\bm{G}_{\beta}\cdot\bm{r}]\exp[-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}]\Xi(\bm{r})\gamma_{\bm{G}_{\beta}}\frac{\Omega}{N_{D}}\sum_{\beta=1}^{N_{G}}\sum_{\bm{r}}\exp[-\mathrm{i}\bm{G}_{\beta}\cdot\bm{r}]\exp[-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}]\Xi(\bm{r})\gamma^{\star}_{\bm{G}_{\beta}}$
(56)
$\displaystyle=2\Omega\sum_{\bm{r}}\exp[-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}]\Xi(\bm{r})\mathrm{IFT}[\bm{\gamma}](\bm{r})=2\Omega\mathrm{FT}\bigr{[}\Xi(\bm{r})\mathrm{IFT}[\bm{\gamma}](\bm{r})\bigl{]}(\bm{G}_{\alpha}),$

where $\mathrm{FT}[\bullet]$ denotes a discrete Fourier transform and
$\mathrm{IFT}[\bullet]$ an inverse discrete Fourier transform. Note that the
discrete inverse Fourier transform is defined including the term
$\bm{G}=\bm{G}_{0}=(0,0,0)$, which is excluded from the constrained degrees of
freedom in our approach. In the numerical implementation of the Fourier
transform above, we thus set $\gamma_{\bm{G}_{0}}=0$ to include this term in
the sum over the vectors in reciprocal space.

Using Fast Fourier Transforms these operations can thus be done with
$\Omega\log\Omega$ scaling. Similar manipulations show that the same is true
also for the other two terms.

##  Appendix B $\Omega\log\Omega$ scaling for matrix-vector multiplications of
nonlocal functionals

For kinetic energy functionals with a nonlocal kernel, we consider the form
discussed by Wang et al. for density-independent kernels wang:1999 , which we
write in the more general form

$T_{\mathrm{K}}^{f,g}=C_{\mathrm{K}}\sum_{\bm{r}}\sum_{\bm{r}^{\prime}}f\bigl{(}n(\bm{r})\bigr{)}w(\bm{r}-\bm{r}^{\prime})g\bigl{(}n(\bm{r}^{\prime})\bigr{)},$
(57)

where $C_{K}$, is a constant, $f$ and $g$ are generic functions of the density
(Wang et al. considered only powers, compare with Eq. (46) of Ref. 28) and
$w(\bm{r}-\bm{r}^{\prime})$ is the nonlocal kernel. Products of this form are
conveniently expressed in Fourier space as

$T_{\mathrm{K}}^{f,g}=C_{\mathrm{K}}\sum_{\bm{G}}\widetilde{F}^{\star}(\bm{G})\widetilde{W}(\bm{G}),\widetilde{G}(\bm{G}),$
(58)

where $\widetilde{F}(\bm{G})$, $\widetilde{G}(\bm{G})$ and
$\widetilde{W}(\bm{G})$ indicates the Fourier transforms of the functions
$f\bigl{(}n(\bm{r})\bigr{)}$, $g\bigl{(}n(\bm{r})\bigr{)}$ and $w(\bm{r})$,
respectively. Note that these Fourier transforms include the term
$\bm{G}=\bm{G}_{0}=(0,0,0)$ and the sum extends from $-\bm{G}_{\mathrm{max}}$
to $\bm{G}_{\mathrm{max}}$.

Let us begin by considering the first derivative of a kinetic energy
functional of this nonlocal form. We have:

$\displaystyle\frac{\partial
T_{\mathrm{K}}^{f,g}}{\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}$
$\displaystyle=C_{\mathrm{K}}\sum_{\bm{r}}\sum_{\bm{r}^{\prime}}\frac{\partial
f}{\partial
n}\bigl{(}n(\bm{r})\bigr{)}e^{-\mathrm{i}\bm{G}_{\alpha}\bm{r}}w(\bm{r}-\bm{r}^{\prime})g\bigl{(}n(\bm{r}^{\prime})\bigr{)}+C_{\mathrm{K}}\sum_{\bm{r}}\sum_{\bm{r}^{\prime}}f\bigl{(}n(\bm{r})\bigr{)}w(\bm{r}-\bm{r}^{\prime})\frac{\partial
g}{\partial
n}\bigl{(}n(\bm{r}^{\prime})\bigr{)}e^{-\mathrm{i}\bm{G}_{\alpha}\bm{r}^{\prime}}$
(59) $\displaystyle=C_{\mathrm{K}}\sum_{\bm{r}}\frac{\partial f}{\partial
n}\bigl{(}n(\bm{r})\bigr{)}e^{-\mathrm{i}\bm{G}_{\alpha}\bm{r}}\mathrm{IFT}[\widetilde{W}(\bm{G})\widetilde{G}(\bm{G})](\bm{r})+C_{\mathrm{K}}\sum_{\bm{r}^{\prime}}\mathrm{IFT}[\widetilde{W}^{\star}(\bm{G})\widetilde{F}(\bm{G})]^{\star}(\bm{r}^{\prime})\frac{\partial
g}{\partial
n}\bigl{(}n(\bm{r}^{\prime})\bigr{)}e^{-\mathrm{i}\bm{G}_{\alpha}\bm{r}^{\prime}}$
$\displaystyle=C_{\mathrm{K}}\mathrm{FT}\biggl{[}\frac{\partial f}{\partial
n}\bigl{(}n(\bm{r})\bigr{)}\mathrm{IFT}[\widetilde{W}(\bm{G})\widetilde{G}(\bm{G})](\bm{r})\biggr{]}(\bm{G}_{\alpha})+C_{\mathrm{K}}\mathrm{FT}\biggl{[}\mathrm{IFT}[\widetilde{W}^{\star}(\bm{G})\widetilde{F}(\bm{G})]^{\star}(\bm{r}^{\prime})\frac{\partial
g}{\partial n}\bigl{(}n(\bm{r}^{\prime})\bigr{)}\biggr{]}(\bm{G}_{\alpha})$

As observed by Wang et al., the evaluation of this term — which only involves
direct and inverse Fourier transforms — can be performed numerically with
$\Omega\log\Omega$ effort.

We now show that also the product of the matrix of second derivatives of the
nonlocal kinetic energy with the $\bm{\gamma}$ needed for the solution of our
constraints (see also first part of this Appendix for GGA functionals) can be
performed with the same numerical scaling. We start from

$\displaystyle\sum_{\beta=1}^{N_{G}}\biggl{\\{}\Sigma^{\dagger}_{T}\biggr{\\}}_{\alpha,\beta}\gamma_{\bm{G}_{\beta}}=\sum_{\beta=1}^{N_{G}}\frac{\partial^{2}T_{\mathrm{K}}^{f,g}}{\partial\tilde{n}_{\bm{G}_{\beta}}\partial\tilde{n}^{\star}_{\bm{G}_{\alpha}}}\gamma_{\bm{G}_{\beta}}=$
$\displaystyle
C_{\mathrm{K}}\sum_{\bm{r}}\sum_{\bm{r}^{\prime}}\frac{\partial^{2}f}{\partial
n^{2}}\bigl{(}n(\bm{r})\bigr{)}e^{-\mathrm{i}\bm{G}_{\alpha}\bm{r}}\biggl{(}\sum_{\beta=1}^{N_{G}}e^{\mathrm{i}\bm{G}_{\beta}\bm{r}}\gamma_{\bm{G}_{\beta}}\biggr{)}w(\bm{r}-\bm{r}^{\prime})g\bigl{(}n(\bm{r}^{\prime})\bigr{)}$
(60) $\displaystyle+$ $\displaystyle
C_{\mathrm{K}}\sum_{\bm{r}}\sum_{\bm{r}^{\prime}}\frac{\partial f}{\partial
n}\bigl{(}n(\bm{r})\bigr{)}e^{-\mathrm{i}\bm{G}_{\alpha}\bm{r}}w(\bm{r}-\bm{r}^{\prime})\frac{\partial
g}{\partial
n}\bigl{(}n(\bm{r}^{\prime})\bigr{)}\biggl{(}\sum_{\beta=1}^{N_{G}}e^{\mathrm{i}\bm{G}_{\beta}\bm{r}^{\prime}}\gamma_{\bm{G}_{\beta}}\biggr{)}$
$\displaystyle+$ $\displaystyle
C_{\mathrm{K}}\sum_{\bm{r}}\sum_{\bm{r}^{\prime}}\frac{\partial f}{\partial
n}\bigl{(}n(\bm{r})\bigr{)}\biggl{(}\sum_{\beta=1}^{N_{G}}e^{\mathrm{i}\bm{G}_{\beta}\bm{r}}\gamma_{\bm{G}_{\beta}}\biggr{)}w(\bm{r}-\bm{r}^{\prime})\frac{\partial
g}{\partial
n}\bigl{(}n(\bm{r}^{\prime})\bigr{)}e^{-\mathrm{i}\bm{G}_{\alpha}\bm{r}^{\prime}}$
$\displaystyle+$ $\displaystyle
C_{\mathrm{K}}\sum_{\bm{r}}\sum_{\bm{r}^{\prime}}f\bigl{(}n(\bm{r})\bigr{)}w(\bm{r}-\bm{r}^{\prime})\frac{\partial^{2}g}{\partial
n^{2}}\bigl{(}n(\bm{r}^{\prime})\bigr{)}\biggl{(}\sum_{\beta=1}^{N_{G}}e^{\mathrm{i}\bm{G}_{\beta}\bm{r}^{\prime}}\gamma_{\bm{G}_{\beta}}\biggr{)}e^{-\mathrm{i}\bm{G}_{\alpha}\bm{r}^{\prime}}$

where we have used the derivative of Eq. (59) with respect to
$\tilde{n}_{\bm{G}_{\beta}}$. Adding to this the second term of the matrix-
vector product
$\sum_{\beta=1}^{N_{G}}\left\\{\Sigma^{\dagger}_{T}\right\\}_{\alpha,N_{G}+\beta}\gamma^{\star}_{\bm{G}_{\beta}}$
and defining
$\gamma_{D}(\bm{r})=\sum_{\beta=1}^{N_{G}}e^{\mathrm{i}\bm{G}_{\beta}\bm{r}}\gamma_{\bm{G}_{\beta}}+\sum_{\beta=1}^{N_{G}}e^{-\mathrm{i}\bm{G}_{\beta}\bm{r}}\gamma^{\star}_{\bm{G}_{\beta}}=\mathrm{FT}[\gamma](\bm{r})$
(where again, in the numerical implementation of the Fourier transform FT, we
set $\gamma_{\bm{G}_{0}}=0$), we can rewrite the matrix-vector product as

$\displaystyle\sum_{\beta=1}^{N_{G}}\left\\{\Sigma^{\dagger}_{T}\right\\}_{\alpha,\beta}\gamma_{\bm{G}_{\beta}}+\sum_{\beta=1}^{N_{G}}\left\\{\Sigma^{\dagger}_{T}\right\\}_{\alpha,N_{G}+\beta}\gamma^{\star}_{\bm{G}_{\beta}}$
$\displaystyle=C_{\mathrm{K}}\sum_{\bm{r}}\sum_{\bm{r}^{\prime}}\frac{\partial^{2}f}{\partial
n^{2}}\bigl{(}n(\bm{r})\bigr{)}e^{-\mathrm{i}\bm{G}_{\alpha}\bm{r}}\gamma_{D}(\bm{r})w(\bm{r}-\bm{r}^{\prime})g\bigl{(}n(\bm{r}^{\prime})\bigr{)}$
(61) $\displaystyle+$ $\displaystyle
C_{\mathrm{K}}\sum_{\bm{r}}\sum_{\bm{r}^{\prime}}\frac{\partial f}{\partial
n}\bigl{(}n(\bm{r})\bigr{)}e^{-\mathrm{i}\bm{G}_{\alpha}\bm{r}}w(\bm{r}-\bm{r}^{\prime})\frac{\partial
g}{\partial n}\bigl{(}n(\bm{r}^{\prime})\bigr{)}\gamma_{D}(\bm{r}^{\prime})$
$\displaystyle+$ $\displaystyle
C_{\mathrm{K}}\sum_{\bm{r}}\sum_{\bm{r}^{\prime}}\frac{\partial f}{\partial
n}\bigl{(}n(\bm{r})\bigr{)}\gamma_{D}(\bm{r})w(\bm{r}-\bm{r}^{\prime})\frac{\partial
g}{\partial
n}\bigl{(}n(\bm{r}^{\prime})\bigr{)}e^{-\mathrm{i}\bm{G}_{\alpha}\bm{r}^{\prime}}$
$\displaystyle+$ $\displaystyle
C_{\mathrm{K}}\sum_{\bm{r}}\sum_{\bm{r}^{\prime}}f\bigl{(}n(\bm{r})\bigr{)}w(\bm{r}-\bm{r}^{\prime})\frac{\partial^{2}g}{\partial
n^{2}}\bigl{(}n(\bm{r}^{\prime})\bigr{)}\gamma_{D}(\bm{r}^{\prime})e^{-\mathrm{i}\bm{G}_{\alpha}\bm{r}^{\prime}}.$

The structure of the four terms in the expression above is the same as that of
the two terms involved in the expression of the first derivative of the
kinetic energy functional, Eq. (59). Using the same manipulations, each one of
them can then be expressed as a sequence of direct and inverse Fourier
transforms, thus leading to an $\Omega\log\Omega$ numerical effort in the
calculations.

##  Appendix C $\Omega\log\Omega$ scaling for computation of the SHAKE matrix
for GGA functionals

In the following we show how the cost of computing the diagonal elements of
the SHAKE matrix (see eq. (32)) can be reduced to $\Omega\log\Omega$ scaling
by exploiting periodicity in direct space. We report below the expression of
these diagonal elements for convenience:

$\displaystyle\left[\mathbb{J}_{D}\right]_{{\bm{G}_{\alpha}},\bm{G}_{\alpha}}=$
$\displaystyle\left(\frac{4\pi\Omega}{|\bm{G}_{\alpha}|^{2}}\right)^{2}+\frac{4\pi\Omega}{|\bm{G}_{\alpha}|^{2}}\frac{\Omega}{N_{D}}\sum_{\bm{r}}\Bigl{[}\Xi_{\text{P}}(\bm{r})+G^{l}_{\alpha}G^{m}_{\alpha}\Psi_{\text{P}}^{lm}(\bm{r})\Bigr{]}+\frac{4\pi\Omega}{|\bm{G}_{\alpha}|^{2}}\frac{\Omega}{N_{D}}\sum_{\bm{r}}\Bigl{[}\Xi_{t}(\bm{r})+G^{l}_{\alpha}G^{m}_{\alpha}\Psi_{t}^{lm}(\bm{r})\Bigr{]}$
$\displaystyle+\biggl{(}\frac{\Omega}{N_{D}}\biggr{)}^{2}\sum_{\bm{r},\bm{r}^{\prime}}e^{-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}}e^{\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}^{\prime}}\sum^{N_{G}}_{\lambda=1}e^{\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}}e^{-\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}^{\prime}}$
$\displaystyle\qquad\times\Bigl{[}\Xi_{\text{P}}(\bm{r})+\mathrm{i}(G^{l}_{\lambda}-G^{l}_{\alpha})\Phi^{l}_{\text{P}}(\bm{r})+G^{l}_{\lambda}G^{m}_{\alpha}\Psi_{\text{P}}^{lm}(\bm{r})\Bigr{]}\Bigl{[}\Xi_{t}(\bm{r}^{\prime})-\mathrm{i}(G^{l}_{\lambda}-G^{l}_{\alpha})\Phi_{t}^{l}(\bm{r}^{\prime})+G^{l}_{\lambda}G^{m}_{\alpha}\Psi_{t}^{lm}(\bm{r}^{\prime})\Bigr{]}$
$\displaystyle+\biggl{(}\frac{\Omega}{N_{D}}\biggr{)}^{2}\sum_{\bm{r},\bm{r}^{\prime}}e^{-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}}e^{\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}^{\prime}}\sum^{N_{G}}_{\lambda=1}e^{-\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}}e^{\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}^{\prime}}$
$\displaystyle\qquad\times\Bigl{[}\Xi_{\text{P}}(\bm{r})-\mathrm{i}(G^{l}_{\lambda}+G^{l}_{\alpha})\Phi_{\text{P}}^{l}(\bm{r})-G^{l}_{\lambda}G^{m}_{\alpha}\Psi_{\text{P}}^{lm}(\bm{r})\Bigr{]}\times\Bigl{[}\Xi_{t}(\bm{r}^{\prime})+\mathrm{i}(G^{l}_{\lambda}+G^{l}_{\alpha})\Phi_{t}^{l}(\bm{r}^{\prime})-G^{l}_{\lambda}G^{m}_{\alpha}\Psi_{t}^{lm}(\bm{r}^{\prime})\Bigr{]},$

which we rewrite as

$\displaystyle\left[\mathbb{J}_{D}\right]_{{\bm{G}_{\alpha}},\bm{G}_{\alpha}}=\left(\frac{4\pi\Omega}{|\bm{G}_{\alpha}|^{2}}\right)^{2}+\frac{4\pi\Omega}{|\bm{G}_{\alpha}|^{2}}\frac{\Omega}{N_{D}}\biggl{(}\sum_{\bm{r}}\Bigl{[}\Xi_{\text{P}}(\bm{r})+G^{l}_{\alpha}G^{m}_{\alpha}\Psi_{\text{P}}^{lm}(\bm{r})\Bigr{]}+\sum_{\bm{r}}\Bigl{[}\Xi_{t}(\bm{r})+G^{l}_{\alpha}G^{m}_{\alpha}\Psi_{t}^{lm}(\bm{r})\Bigr{]}\biggr{)}+\left[\mathbb{J}^{\text{KXC}}_{D}\right]_{{\bm{G}_{\alpha}},\bm{G}_{\alpha}}$
(62)

The first term involving the Hartree contribution, as well as the mixed terms,
have a numerical cost that trivially scales linearly with system size. We will
discuss in this appendix the pure kinetic and exchange-correlation
contribution
$\left[\mathbb{J}^{\text{KXC}}_{D}\right]_{{\bm{G}_{\alpha}},\bm{G}_{\alpha}}$.
Let us begin by observing that performing the products in the triple sum above
results in a total of eight terms, all of the generic form

$\displaystyle h^{ab}_{1}(\bm{G}_{\alpha})$
$\displaystyle=\sum_{\bm{r}^{\prime}=0}^{\bm{L}}\sum_{\bm{r}=0}^{\bm{L}}\sum_{\lambda=1}^{N_{G}}\biggl{[}e^{-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}}e^{\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}}e^{\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}^{\prime}}e^{-\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}^{\prime}}\times(\mathrm{i}G^{l}_{\lambda})^{a}(\mathrm{i}G^{m}_{\lambda})^{b}\phi(\bm{r})\psi(\bm{r}^{\prime})\biggr{]}$
(63)

where we have explicitly indicated the range of the sums in direct space,
$a,b=0,1$ and the functions $\phi$ and $\psi$ are proxies for the $\Xi$ and
$\Phi$ computed at iteration $\kappa$ or at time $t$ values of the electronic
coefficients. Focusing at first on the sum over $\bm{r}$, note that, since
both $\Phi$ and $\Xi$ are functions of the (periodic) energy of the system,
the function
$\phi(\bm{r})e^{-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}}e^{\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}}$
is periodic in $\bm{r}$ with period $\bm{L}$. The sums in the equation above
can then be reorganized as follows

$\displaystyle h^{ab}_{1}(\bm{G}_{\alpha})$
$\displaystyle=\sum_{\bm{r}^{\prime}=0}^{\bm{L}}\sum_{\bm{r}=\bm{r}^{\prime}}^{\bm{L}+\bm{r}^{\prime}}\sum_{\lambda=1}^{N_{G}}\biggl{[}e^{-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}}e^{\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}}e^{\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}^{\prime}}e^{-\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}^{\prime}}\times(\mathrm{i}G^{l}_{\lambda})^{a}(\mathrm{i}G^{m}_{\lambda})^{b}\phi(\bm{r})\psi(\bm{r}^{\prime})\biggr{]}$
(64)

from which, using the transformation
$\bm{r}=\bm{r}^{\prime}+\bm{r}^{\prime\prime}$,

$\displaystyle h^{ab}_{1}(\bm{G}_{\alpha})$
$\displaystyle=\sum_{\bm{r}^{\prime}=0}^{\bm{L}}\sum_{\bm{r}^{\prime\prime}=0}^{\bm{L}}\sum_{\lambda=1}^{N_{G}}\biggl{[}e^{-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}^{\prime\prime}}e^{\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}^{\prime\prime}}\times(\mathrm{i}G^{l}_{\lambda})^{a}(\mathrm{i}G^{m}_{\lambda})^{b}\phi(\bm{r}^{\prime}+\bm{r}^{\prime\prime})\psi(\bm{r}^{\prime})\biggr{]}$
(65)

The terms in Eq. (65) can now be rearranged to obtain

$\displaystyle h^{ab}_{1}(\bm{G}_{\alpha})$
$\displaystyle=\sum_{\bm{r}^{\prime\prime}=0}^{\bm{L}}e^{-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}^{\prime\prime}}\biggl{[}\biggl{(}\sum_{\lambda=1}^{N_{G}}e^{\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}^{\prime\prime}}(\mathrm{i}G^{l}_{\lambda})^{a}(\mathrm{i}G^{m}_{\lambda})^{b}\biggr{)}\biggl{(}\sum_{\bm{r}^{\prime}=0}^{\bm{L}}\phi(\bm{r}^{\prime}+\bm{r}^{\prime\prime})\psi(\bm{r}^{\prime})\biggr{)}\biggr{]}$
(66)

The sum over $\bm{r}^{\prime}$ in the equation above is, by definition, the
spatial correlation function
$\left(\phi\star\psi\right)(\bm{r}^{\prime\prime})=\sum_{\bm{r}^{\prime}=0}^{\bm{L}}\phi(\bm{r}^{\prime}+\bm{r}^{\prime\prime})\psi(\bm{r}^{\prime})$.

The second triple sum in Eq. (32) can be treated in the same way to obtain

$\displaystyle h^{ab}_{2}(\bm{G}_{\alpha})$
$\displaystyle=\sum_{\bm{r}^{\prime\prime}=0}^{\bm{L}}e^{-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}^{\prime\prime}}\biggl{[}\biggl{(}\sum_{\lambda=1}^{N_{G}}e^{-\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}^{\prime\prime}}(-\mathrm{i}G^{l}_{\lambda})^{a}(-\mathrm{i}G^{m}_{\lambda})^{b}\biggr{)}\biggl{(}\sum_{\bm{r}^{\prime}=0}^{\bm{L}}\phi(\bm{r}^{\prime}+\bm{r}^{\prime\prime})\psi(\bm{r}^{\prime})\biggr{)}\biggr{]}$
(67)

Combining these results, we have

$h^{ab}_{1}(\bm{G}_{\alpha})+h^{ab}_{2}(\bm{G}_{\alpha})=\sum_{\bm{r}^{\prime\prime}=0}^{\bm{L}}\left(\phi\star\psi\right)(\bm{r}^{\prime\prime})e^{-\mathrm{i}\bm{G}_{\alpha}\cdot\bm{r}^{\prime\prime}}\sum_{\lambda=1}^{N_{G}}\biggl{[}e^{\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}^{\prime\prime}}(\mathrm{i}G^{l}_{\lambda})^{a}(\mathrm{i}G^{m}_{\lambda})^{b}+e^{-\mathrm{i}\bm{G}_{\lambda}\cdot\bm{r}^{\prime\prime}}(-\mathrm{i}G^{l}_{\lambda})^{a}(-\mathrm{i}G^{m}_{\lambda})^{b}\biggr{]}$
(68)

The expression above is particularly convenient for our purposes because the
remaining sums over direct and reciprocal space are all, directly or via
trivial manipulations, direct or inverse Fourier transforms that can be
computed at $\Omega\log\Omega$ (or equivalently $N_{D}\log N_{D}$) cost using
libraries implementing the fast Fourier transform (FFT) algorithm. More in
detail, the sum over $\bm{G}_{\lambda}$ is, by definition, the inverse
(discrete) Fourier transform of the product
$(iG^{l}_{\lambda})^{a}(iG^{m}_{\lambda})^{b}$ multiplied by $N_{D}$, while
the correlation function $\left(\phi\star\psi\right)$ can be computed as

$\left(\phi\star\psi\right)(\bm{r}^{\prime\prime})=\sum_{\bm{r}^{\prime}}\phi(\bm{r}^{\prime}+\bm{r}^{\prime\prime})\psi(\bm{r}^{\prime})=N_{\text{D}}\mathrm{IFT}\bigl{[}\mathrm{FT}[\phi]^{\star}\mathrm{FT}[\psi]\bigr{]}$
(69)

via the convolution theorem. Thus,

$\displaystyle h^{ab}_{1}(\bm{G}_{\alpha})+h^{ab}_{2}(\bm{G}_{\alpha})$
$\displaystyle=N^{2}_{D}\sum_{a,b}\operatorname{FT}\biggl{[}\biggl{(}\operatorname{IFT}[(\mathrm{i}G^{l}_{\lambda})^{a}(\mathrm{i}G^{m}_{\lambda})^{b}]\biggr{)}\biggl{(}\operatorname{IFT}\bigl{[}\operatorname{FT}[\phi]^{\star}\operatorname{FT}[\psi]\bigr{]}\biggr{)}\biggr{]}$
(70)

Note that the FTs in Eq. (70) return a vector of length $N_{D}$ whose
components are the $N_{D}$ nonzero elements of $h^{ab}_{1}+h^{ab}_{2}$. This
implies that, together with the other terms in Eq. (62), the calculation of
the SHAKE denominator is reduced to a number of $\mathcal{O}(N_{\text{D}}\log
N_{\text{D}})$ operations equal to the number of terms of the form of the LHS
of Eq. (70). Finally, the diagonal part of the KXC contribution to the
hessian, $\mathbb{J}^{\text{KXC}}_{D}$, is obtained as a linear combination of
these terms $h^{ab}_{1}+h^{ab}_{2}$, with some coefficients involving a
product with $G^{l}_{\alpha}$ or $G^{m}_{\alpha}$. This however comes at
nearly no extra cost and is in any case linear in the system size so that the
overall cost proportional to $\mathcal{O}(N_{\text{D}}\log N_{\text{D}})$ is
achieved.

## References

  * (1) Alessandro Coretti, Sara Bonella, and Giovanni Ciccotti.  Communication: Constrained molecular dynamics for polarizable models.  The Journal of Chemical Physics, 149(19):191102, 2018.
  * (2) Sara Bonella, Alessandro Coretti, Rodolphe Vuilleumier, and Giovanni Ciccotti.  Adiabatic motion and statistical mechanics via mass-zero constrained dynamics.  Physical Chemistry Chemical Physics, page 10.1039.D0CP00163E, 2020\.
  * (3) Andrés Aguado, Leonardo Bernasconi, and Paul A. Madden.  Interionic potentials from _ab initio_ molecular dynamics: The alkaline earth oxides CaO, SrO, and BaO.  The Journal of Chemical Physics, 118(13):5704–5717, 2003.
  * (4) Sandro Jahn, Paul A. Madden, and Mark Wilson.  Dynamic simulation of pressure-driven phase transformations in crystalline ${\mathrm{al}}_{2}{\mathrm{o}}_{3}$.  Physical Review B, 69:020106, Jan 2004.
  * (5) J. Kolafa.  Time-reversible always stable predictor-corrector method for molecular dynamics of polarizable molecules.  Journal of Computational Chemistry, 25(3):335–342, 2004.
  * (6) Jiří Genzer and Jiří Kolafa.  Molecular dynamics of potential models with polarizability: comparison of methods.  Journal of Molecular Liquids, 109(2):63 – 72, 2004.  Sixth Liblice Conference on the Statistical Mechanics of Liquids.
  * (7) R. Car and M. Parrinello.  Unified approach for molecular dynamics and density-functional theory.  Physical Review Letters, 55:2471–2474, Nov 1985.
  * (8) G. Pastore, E. Smargiassi, and F. Buda.  Theory of _ab initio_ molecular-dynamics calculations.  Physical Review A, 44:6334–6347, Nov 1991.
  * (9) Anders M. N. Niklasson, C. J. Tymczak, and Matt Challacombe.  Time-reversible Born–Oppenheimer molecular dynamics.  Physical Review Letters, 97:123001, Sep 2006.
  * (10) Anders M. N. Niklasson.  Extended Lagrangian Born–Oppenheimer molecular dynamics using a Krylov subspace approximation.  The Journal of Chemical Physics, 152(10):104103, 2020.
  * (11) Anders MN Niklasson.  Extended lagrangian born–oppenheimer molecular dynamics for orbital-free density-functional theory and polarizable charge equilibration models.  The Journal of Chemical Physics, 154(5):054101, 2021.
  * (12) M Pounds, S Tazi, M Salanne, and P A Madden.  Ion adsorption at a metallic electrode: an _ab initio_ based simulation study.  Journal of Physics: Condensed Matter, 21(42):424109, sep 2009.
  * (13) Fabien Pacaud, Mathieu Salanne, Thibault Charpentier, Laurent Cormier, and Jean-Marc Delaye.  Structural study of Na2O-B2O3-SiO2-La2O3 glasses from molecular simulations using a polarizable force field.  Journal of Non-Crystalline Solids, 499:371 – 379, 2018.
  * (14) Dahlia K Remler and Paul A Madden.  Molecular dynamics without effective potentials via the car-parrinello approach.  Molecular Physics, 70(6):921–966, 1990.
  * (15) Peter Pulay and Géza Fogarasi.  Fock matrix dynamics.  Chemical physics letters, 386(4-6):272–278, 2004.
  * (16) MC Payne, JD Joannopoulos, DC Allan, MP Teter, and David H Vanderbilt.  Molecular dynamics and ab initio total energy calculations.  Physical review letters, 56(24):2656, 1986.
  * (17) Mike C. Payne, Michael P. Teter, Douglas C. Allan, TA Arias, and ad JD Joannopoulos.  Iterative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradients.  Reviews of modern physics, 64(4):1045, 1992.
  * (18) Dominik Marx and J urg Hutter.  Ab initio molecular dynamics: Basic theory and advanced methods.  Cambridge University Press, Cambridge, 2012.
  * (19) Herbert Goldstein.  Classical mechanics.  Addison Wesley, San Francisco, 2002.
  * (20) G. Ciccotti and J.-P. Ryckaert.  Molecular dynamics simulation of rigid molecules.  Computer Physics Reports, 4(6):346 – 392, 1986.
  * (21) Jean-Paul Ryckaert, Giovanni Ciccotti, and Herman J.C Berendsen.  Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes.  Journal of Computational Physics, 23(3):327 – 341, 1977.
  * (22) A. Coretti, L. Scalfi, C. Bacon, B. Rotenberg, R. Vuilleumier, G. Ciccotti, M. Salanne, and S. Bonella.  Mass-zero constrained molecular dynamics for electrode charges in simulations of electrochemical systems.  The Journal of Chemical Physics, 152(19):194701, May 2020.
  * (23) B. Leimkuhler and C. Matthews.  Molecular dynamics : with deterministic and stochastic numerical methods.  Springer, Cham, 2015.
  * (24) Abel Marin-Laflèche, Matthieu Haefele, Laura Scalfi, Alessandro Coretti, Thomas Dufils, Guillaume Jeanmairet, Stewart K. Reed, Alessandra Serva, Roxanne Berthin, Camille Bacon, Sara Bonella, Benjamin Rotenberg, Paul A. Madden, and Mathieu Salanne.  Metalwalls: A classical molecular dynamics software dedicated to the simulation of electrochemical systems.  Journal of Open Source Software, 5(53):2373, 2020.
  * (25) D. D. Girardier, A. Coretti, G. Ciccotti, and S. Bonella.  Mass-zero constrained dynamics and statistics for the shell model in magnetic field.  The European Physical Journal B, 94(8):158, 2021.
  * (26) Yael Weinbach and Ron Elber.  Revisiting and parallelizing SHAKE.  Journal of Computational Physics, 209(1):193 – 206, 2005.
  * (27) Lin-Wang Wang and Michael P. Teter.  Kinetic-energy functional of the electron density.  Physical Review B, 45:13196–13220, Jun 1992.
  * (28) Yan Alexander Wang, Niranjan Govind, and Emily A. 1.  Orbital-free kinetic-energy density functionals with a density-dependent kernel.  Physical Review B, 60(24):16350–16358, December 1999.  Publisher: American Physical Society.
  * (29) M Pearson, E Smargiassi, and P A Madden.  _Ab initio_ molecular dynamics with an orbital-free density functional.  Journal of Physics: Condensed Matter, 5(19):3221–3240, may 1993\.
  * (30) Vincent L. Lignères and Emily A. Carter.  An Introduction to Orbital-Free Density Functional Theory, pages 137–148.  Springer Netherlands, Dordrecht, 2005.
  * (31) William C. Witt, Beatriz G. del Rio, Johannes M. Dieterich, and Emily A. Carter.  Orbital-free density functional theory for materials research.  Journal of Materials Research, 33(7):777–795, 2018.
  * (32) J. J. Mortensen, L. B. Hansen, and K. W. Jacobsen.  Real-space grid implementation of the projector augmented wave method.  Phys. Rev. B, 71:035109, Jan 2005.
  * (33) J Enkovaara, C Rostgaard, J J Mortensen, J Chen, M Dułak, L Ferrighi, J Gavnholt, C Glinsvad, V Haikola, H A Hansen, H H Kristoffersen, M Kuisma, A H Larsen, L Lehtovaara, M Ljungberg, O Lopez-Acevedo, P G Moses, J Ojanen, T Olsen, V Petzold, N A Romero, J Stausholm-Møller, M Strange, G A Tritsaris, M Vanin, M Walter, B Hammer, H Häkkinen, G K H Madsen, R M Nieminen, J K Nørskov, M Puska, T T Rantala, J Schiøtz, K S Thygesen, and K W Jacobsen.  Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method.  Journal of Physics: Condensed Matter, 22(25):253202, jun 2010.
  * (34) Mohan Chen, Junchao Xia, Chen Huang, Johannes M Dieterich, Linda Hung, Ilgyou Shin, and Emily A Carter.  Introducing PROFESS 3.0: An advanced program for orbital-free density functional theory molecular dynamics simulations.  Computer Physics Communications, 190:228–230, 2015.
  * (35) Ralf Meyer, Manuel Weichselbaum, and Andreas W. Hauser.  Machine learning approaches toward orbital-free density functional theory: Simultaneous training on the kinetic energy density functional and its functional derivative.  Journal of Chemical Theory and Computation, 16(9):5685–5694, 2020\.  PMID: 32786898.
  * (36) S. Alireza Ghasemi and Thomas D. Kühne.  Artificial neural networks for the kinetic energy functional of non-interacting fermions.  The Journal of Chemical Physics, 154(7):074107, February 2021.  Publisher: American Institute of Physics.
  * (37) Jorge Nocedal and Stephen J. Wright.  Numerical optimization.  Springer, New York, 2006.
  * (38) Eric Barth, Krzysztof Kuczera, Benedict Leimkuhler, and Robert D. Skeel.  Algorithms for constrained molecular dynamics.  Journal of Computational Chemistry, 16(10):1192–1209, October 1995\.
  * (39) Benedict J. Leimkuhler and Robert D. Skeel.  Symplectic numerical integrators in constrained Hamiltonian systems.  Journal of Computational Physics, 112(1):117 – 125, 1994.
  * (40) B. Leimkuhler and S. Reich.  Simulating Hamiltonian dynamics.  Cambridge University Press, Cambridge, UK New York, 2004.
  * (41) Susi Lehtola, Conrad Steigemann, Micael J.T. Oliveira, and Miguel A.L. Marques.  Recent developments in libxc — A comprehensive library of functionals for density functional theory.  SoftwareX, 7:1–5, January 2018.
  * (42) John P. Perdew, J. A. Chevary, S. H. Vosko, Koblar A. Jackson, Mark R. Pederson, D. J. Singh, and Carlos Fiolhais.  Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation.  Phys. Rev. B, 46:6671–6687, Sep 1992.
  * (43) von Weizsäcker, C. F.  Zur theorie der kernmassen.  Zeitschrift für Physik, 96(7-8):431–458, 1935.
  * (44) F Perrot.  Hydrogen-hydrogen interaction in an electron gas.  Journal of Physics: Condensed Matter, 6(2):431–446, January 1994\.
  * (45) Pascal Joly and Gèrard Meurant.  Complex conjugate gradient methods.  Numerical Algorithms, 4(3):379–406, October 1993.
  * (46) P. A. M. Dirac.  Note on exchange phenomena in the thomas atom.  Mathematical Proceedings of the Cambridge Philosophical Society, 26(3):376–385, 1930.
  * (47) F. Bloch.  Bemerkung zur elektronentheorie des ferromagnetismus und der elektrischen leitfähigkeit.  Zeitschrift für Physik, 57(7):545–555, 1929.
  * (48) J. P. Perdew and Alex Zunger.  Self-interaction correction to density-functional approximations for many-electron systems.  Phys. Rev. B, 23:5048–5079, May 1981.
  * (49) William C. Topp and John J. Hopfield.  Chemically motivated pseudopotential for sodium.  Physical Review B, 7:1295–1303, Feb 1973.
  * (50) Eric Vanden-Eijnden and Giovanni Ciccotti.  Second-order integrators for Langevin equations with holonomic constraints.  Chemical Physics Letters, 429(1-3):310–316, September 2006.
  * (51) Guo-Xin Qian, M. Weinert, G. W. Fernando, and J. W. Davenport.  First-principles calculation of the activation energy for diffusion in liquid sodium.  Physical Review Letters, 64(10):1146–1149, March 1990.
  * (52) A. J. Greenfield, J. Wellendorf, and Nathan Wiser.  X-ray determination of the static structure factor of liquid na and k.  Phys. Rev. A, 4:1607–1616, Oct 1971.
  * (53) T. Lee, J. Bisschop, W. van der Lugt, and W.F. Van Gunsteren.  Radial distribution functions of liquid na and cs.  Physica B+C, 93(1):59–62, 1978.
  * (54) Fuyuki Shimojo, Kozo Hoshino, and Mitsuo Watabe.  Dynamical correlation functions and memory functions of liquid sodium—a molecular dynamics simulation.  Journal of the Physical Society of Japan, 63(1):141–155, 1994.
  * (55) M. W. Ozelton and R. A. Swalin.  Self-diffusion in liquid sodium at constant volume and constant pressure.  Philosophical Magazine, 18(153):441–451, September 1968.
  * (56) William H. Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling.  Numerical recipes in C: The art of scientific computing.  Cambridge University Press, Cambridge Cambridgeshire New York, 1992.
  * (57) Anders M. N. Niklasson and Marc J. Cawkwell.  Generalized extended Lagrangian Born-Oppenheimer molecular dynamics.  The Journal of Chemical Physics, 141(16):164123, October 2014.