Generalized Extended Lagrangian Born-Oppenheimer Molecular Dynamics

Extended Lagrangian Born-Oppenheimer molecular dynamics based on Kohn-Sham density functional theory is generalized in the limit of vanishing self-consistent field optimization prior to the force evaluations. The equations of motion are derived directly from the extended Lagrangian under the condition of an adiabatic separation between the nuclear and the electronic degrees of freedom. We show how this separation is automatically fulfilled and system independent. The generalized equations of motion require only one diagonalization per time step and are applicable to a broader range of materials with improved accuracy and stability compared to previous formulations.


I. INTRODUCTION
Born-Oppenheimer molecular dynamics simulations, where classical molecular trajectories are propagated by forces that are calculated on-the-fly from the relaxed electronic ground state in each time step, provide a powerful tool in materials science, chemistry and biology [1,2]. However, the computational cost is high. In regular Born-Oppenheimer molecular dynamics simulations based on Kohn-Sham density functional theory [3][4][5][6], the cost scales with the cube of the number of atoms, which is multiplied by a significant prefactor given by the number of self-consistent field iterations that are required to find the relaxed electronic ground state prior to the force evaluation in each time step. If a sufficient degree of convergence is not achieved, the electronic system behaves like a heat sink or source, gradually draining or adding energy to the atomic system due to a broken time-reversibility in the propagation of the underlying electronic degrees of freedom [7][8][9][10][11]. Recently, an extended Lagrangian formulation of Born-Oppenheimer molecular dynamics was proposed that overcomes this problem [12][13][14][15][16][17]. By including an auxiliary electron density as a dynamical variable, the equations of motion can be integrated without breaking time reversibility, even under approximate selfconsistent field convergence. This framework enables Born-Oppenheimer molecular dynamics simulations [18,19] with accurate long-term energy conservation using fast linear scaling solvers [20][21][22] that can treat much larger systems than previously possible.
The extended Lagrangian formulation was recently investigated in the limit of vanishing self-consistent field convergence prior to the force calculations. This fast approach requires only one diagonalization per time step [17,23,24]. In this limit, the iterative ground state optimization, which is necessary in regular Born-Oppenheimer molecular dynamics, is replaced by a dynamical optimization acting over time [25,26]. However, this fast self-optimizing dynamics relies on the ability to keep the electronic degrees of freedom close to the exact self-consistent ground state, which is approximated using a simple linear mixing. For many materials systems exhibiting, for example, bond breaking or charge sloshing, this linear mixing is insufficient and leads to instabilities. In this article we revisit extended Lagrangian Born-Oppenheimer molecular dynamics in the limit of vanishing self-consistent field optimization and propose a more general formulation that is applicable to a broader range of materials. In contrast to previous formulations we also derive the equations of motion directly from the extended Lagrangian without intermediate steps. This is possible in an adiabatic limit under the condition of a frequency separation between the nuclear and the electronic degrees of freedom, which we show is automatically fulfilled.
We present our generalized extended Lagrangian Born-Oppenheimer molecular dynamics in terms of Kohn-Sham density functional theory [3][4][5][6] and we use the electron density as the variable for the electronic degrees of freedom. However, the framework should be generally applicable also to other electronic structure formulations, for example, density matrices in Hartree-Fock calculations and wavefunctions or Green's functions in methods for strongly correlated electrons.

A. Generalized extended Born-Oppenheimer Lagrangian
Our generalized extended Born-Oppenheimer Lagrangian is defined by The first term is the kinetic energy for the nuclear degrees of freedom, i.e. the atomic coordinates R = {R I } and the velocitiesṘ = {Ṙ I }, where the dot denotes the time derivative. The second term is a modified Kohn-Sham Born-Oppenheimer potential energy, U (R, n), including nuclear-nuclear repulsion. A more precise definition of U (R, n) is given in Sec. II B below. The last two terms form a harmonic oscillator for the extended electronic degrees of freedom, i.e. the electron density n(r) and its time derivativeṅ(r). Here µ is an electron mass parameter and ω is the frequency determining the curvature of the harmonic well centered around the optimized density ρ(r). The electron density ρ(r) is a fully relaxed and optimized ground state density, but for an approximate linearized electron-electron interaction. Details about the definition of ρ(r) are given below in Sec. II B. In contrast to the previous formulation of extended Lagrangian Born-Oppenheimer molecular dynamics [12], the harmonic potential includes a non-local kernel K(r, r ′ ). This kernel is defined by the relation A more detailed explanation for this particular choice, which forces the dynamical variable n(r) to evolve close to the exact fully optimized ground state density, is given in Sec. II C below. However, a constant density-independent kernel can also be used.

B. Modified Born-Oppenheimer potential energy surface
To construct the modified potential energy surface U (R, n) in Eq. (1), we start with the Kohn-Sham Born-Oppenheimer potential energy surface U (R) [4][5][6], which is given by a constrained functional minimization, where The first term in Eq. (4) is the kinetic Kohn-Sham energy for an electron density ρ(r) that is given from an occupied (occ) set of orthonormalized orbitals {ψ i }. The second term contains the external or pseudopotential potential, V ext (R, r), i.e. the electron-ion interaction, and the third term is the electron-electron interaction energy, consisting of the Hartree term and the exchange-correlation energy. The V nn (R) term in Eq. (3) is the ion-ion repulsion energy. The constrained minimization in Eq. (3) is performed over all normalized electron densities, representable by orthonormal orbitals. The minimization defines a ρ-independent potential energy surface U (R), which can be calculated through a non-linear iterative self-consistent field optimization procedure that involves solving a sequence of Kohn-Sham eigenvalue equations. The exact self-consistent (sc) ground state density, ρ sc (r), is the density for which the minimum U (R) is obtained in Eq. (3).
The modified Born-Oppenheimer potential energy surface, U(R, n) in Eq. (1), is constructed as U (R) if we approximate the electron-electron interaction term E ee [ ρ] with a functional linearization around the electron density n(r), Our modified Born-Oppenheimer potential energy surface is then given by the constrained minimization, where the electron-electron interaction in Eq. (5) has been replaced by the linearized expression in Eq. (6), i.e.
The density n(r) appears in U(R, n) as a dynamical variable in the same way as R. A given set of dynamical variables, n(r) and R, determine the relaxed optimized ground state density ρ(r), which is attained at the minimum that defines U(R, n) in Eq. (7), in the same way as the non-linear iterative optimization of regular Kohn-Sham density functional theory determines the exact self-consistent ground state density ρ sc (r) in Eq. (3), i.e.
The main difference is that the constrained minimum of the linearized form in Eq. (7) or (9) To highlight this direct relationship between the optimized electron density ρ(r) from Eq. (9) and the dynamical variable density n(r) we may use the equivalent notation ρ[n](r) ≡ ρ(r), as above. The leading difference between the modified potential energy surface U(R, n) and the regular Kohn-Sham potential U (R), which follows from the linearization, is of second order with respect to the error in the density. As long as the dynamical variable density n(r) evolves close to the exact ground state density ρ sc (r), the difference between the modified potential energy surface U(R, n) and the exact Kohn-Sham potential U (R) will be small. This motivates our specific choice of the generalized kernel in Eq. (2) as will be shown in the next section.

C. Ground state optimization kernel
Our particular choice of the generalized kernel K(r, r ′ ) in Eq. (2) forces the dynamical variable density n(r) to evolve close to the exact self-consistent ground state density ρ sc (r), such that the modified Born-Oppenheimer potential energy surface U(R, n) is an accurate approximation to the regular self-consistent Kohn-Sham potential U (R).
Let the functional F (n) be the difference between ρ[n](r) and n(r), i.e.
The equation F (n) = 0 is fulfilled for the exact self-consistent ground state density, i.e. when n(r) = ρ sc (r). If n(r) is sufficiently close to the exact ground state density ρ sc (r), we may approximate the ground state density through a Newton minimization step [27], which schematically is given by where J is the Jacobian, δF/δn. This means that or in a more explicit form where the non-local operator K(r, r ′ ) is the inverse Jacobian of F (n) = ρ[n] − n defined by which corresponds to our definition of the kernel in Eq. (2). The error in the approximation in Eq. (14), at least under reasonable conditions, is of second order, i.e. O (ρ − n) 2 , which can be understood from the quadratic convergence of the Newton scheme for sufficiently close initial values of n and assuming F (n) is smooth [27]. Our particular choice of the non-local kernel in the generalized extended Lagrangian, Eq. (1), is thus derived from an optimization of the extended harmonic potential to accurately mimic the oscillation of n(r) around the exact ground state density ρ sc (r). This improves the accuracy of the modified potential energy surface U(R, n) and helps stabilize the dynamics.

D. Equations of motion
Previous derivations of the equations of motion in the limit of vanishing self-consistent field convergence [17,24] started from the decoupled equations of motion of extended Lagrangian Born-Oppenheimer molecular dynamics in the limit of full self-consistent field convergence. Then, a simplified set of equations of motion were derived by relaxing the condition of the self-consistent field optimization. Such an approach can be difficult to interpret. Here we will follow a more direct path in the derivation, without assuming an intermediate step of full self-consistent field convergence. Instead, we will apply an adiabatic condition that assumes a separation in the frequencies between the nuclear and the electronic degree of freedom.
The Euler-Lagrange equations of motion for the generalized extended Lagrangian, L(R,Ṙ, n,ṅ) in Eq. (1), describe a coupled dynamical system of nuclear degrees of freedom with some unknown highest characteristic frequency Ω and an electronic degrees of freedom with the frequency ω. For normal integration time steps used in regular Born-Oppenheimer molecular dynamics we will show in Sec. II G that ω > Ω and that we have a natural system-independent adiabatic separation between the two frequencies. This natural adiabatic separation between slower nuclear and faster electronic degrees of freedom motivates our derivation of the equation of motion in the limit when ω 2 /Ω 2 → ∞. In this limit the difference between ρ(r) and n(r) is small and scales as the inverse of the square of the electron frequency ω 2 , i.e. (ρ[n] − n) = O(Ω 2 /ω 2 ). This assumed scaling of (ρ[n] − n) as a function of ω 2 /Ω 2 → ∞ is difficult to prove rigorously, but it is straightforward to demonstrate that it holds a posteriori (See Fig. 3). In this adiabatic limit we further choose to reduce the electron mass parameter µ → 0 such that lim µω → constant. The Euler-Lagrange equations of motion in this limit are where the forces in the first equation are given from the partial derivative of U(R, n) under the condition of a constant density n(r), since it is a dynamical variable, and where K(r, r ′ ) is defined by Eq. (2). A more detailed derivation of the equations of motion is given in the appendix. The modified Born-Oppenheimer potential energy surface U(R, n) is calculated through a single diagonalization, which yields the optimized density ρ(r), Eq. (10). No additional diagonalization of the Kohn-Sham Hamiltonian is necessary to integrate the equations of motion above -only one Hamiltonian diagonalization per molecular dynamics time step is necessary. However, the additional calculation of the kernel K(r, r ′ ) may introduce a significant overhead.
The key results of this paper are the equations of motion, Eq. (16), and their derivation from the generalized extended Lagrangian, Eq. (1), under the condition of an adiabatic separation between the frequencies of the nuclear and the electronic degrees of freedom. This allows stable molecular dynamics simulations for a broader range of materials systems than previously possible. The extended range of systems that can be treated with the new generalized formalism can be understood from the stability conditions of our previous formulation of extended Lagrangian Born-Oppenheimer molecular dynamics, as will be discussed in the next section.

E. Stability
It is easy to see that previous formulations of extended Lagrangian Born-Oppenheimer molecular dynamics in the limit of vanishing self-consistent field optimization [17,23,24] are recovered with a scaled δ-function approximation of the kernel, i.e. for K(r, r ′ ) = −c δ(r − r ′ ) in Eq. (16), where c ∈ [0, 1]. This corresponds to approximating the exact ground state density ρ sc (r) by a simple linear mixing, i.e. where so that The equation of motion for the electronic degrees of freedom with the scaled δ-function approximation of the kernel in Eq. (16),n therefore corresponds to a dynamics where the auxiliary density n(r) evolves around a self-consistent ground state density that is approximated by a linear mixing between ρ[n] and n. This linear mixing provides a reasonable approximation of the self-consistent ground state density if the electronic energy functional E(R, ρ) has a simple convex form, such that This allows a stable integration of the electronic degrees of freedom in Eq. (16) using the local scaled δ-function approximation of the kernel and requires only one Kohn-Sham diagonalization per time step. However, instabilities may occur, for example, during bond breaking or in large metallic systems exhibiting charge sloshing. In this case ρ(r) has to be further updated through an iterative self-consistent field optimization. Thus, while previous formulations, in the limit requiring only one diagonalization per time step, rely on stability and convergence using a constant linear mixing, which is possibly guaranteed only for convex functionals [28], the generalized formulation with a nonlocal kernel will have a broader range of applications. We can expect that the difference in applicability should be related to the range of systems for which a simple constant linear mixing between densities is sufficient to reach the self-consistent ground state solution in a regular Kohn-Sham calculation [28], in comparison to problems where Newton-like iterations, such a Broyden or Kerker mixing are necessary [29][30][31][32].

F. Constant of motion
The constrained optimization of ρ(r) in Eq. (9) that defines U(R, n) in Eq. (7) and thus the nuclear forces in Eq. (16), corresponds to a ground state relaxation of the electronic degrees of freedom with respect to the regular Kohn-Sham energy functional as in Eq. (3), but with a linearized expression for the electron-electron interaction, Eq. (6). In this way U(R, n) as defined by Eq. (7) can be understood in terms of an optimized and thus ρ-independent Born-Oppenheimer potential energy surface with a relaxed electron density, but which, thanks to the linearized electron-electron interaction term, can be obtained through a single diagonalization, Eq. (10). We may therefore interpret our one-diagonalization-per-time-step extended Lagrangian Born-Oppenheimer molecular dynamics as an exact Kohn-Sham Born-Oppenheimer molecular dynamics simulation, but with a constant of motion of a modified shadow Hamiltonian, which is valid in the same adiabatic limit as the equations on motion, Eq. (16). Thus, instead of integrating the exact Kohn-Sham Born-Oppenheimer equations of motion using an approximate self-consistent field optimization that leads to a systematic drift in the energy [7,8], we perform the integration based on an exact, fully optimized ground state density but for an approximate shadow Hamiltonian. In this way the energy drift of regular Born-Oppenheimer molecular dynamics can be avoided. The argument is conceptually very similar to the backward analysis for symplectic integration schemes [47]. As long as ρ(r) and n(r) are close to the self-consistent ground state density ρ sc (r), E Tot and the potential energy surface U(R, n) will be very close to exact Kohn-Sham Born-Oppenheimer molecular dynamics. This condition is automatically fulfilled when there is an adiabatic separation between ρ(r) and n(r). As will be shown in the next section, this separation is naturally fulfilled and material independent.

G. Adiabatic separation
Under the condition of an adiabatic separation between a slower nuclear motion characterized by a frequency Ω, and a faster auxiliary electronic motion governed by the frequency of the extended harmonic oscillator ω, we may expect, based on the classical adiabatic theorem, that the dynamical variable density n(r), which is controlled by the harmonic oscillator frequency ω, closely should follow the exact ground state density ρ sc (r), which is governed by a slower nuclear frequency Ω. As we will show here, this condition of adiabatic separation between Ω and ω is automatically fulfilled. We can make a straightforward estimate of this separation assuming a standard Verlet integration of n(r) in the equations of motion in Eq. (16), It can be shown that the highest possible value of the dimensionless variable κ = δt 2 ω 2 that provides a stable Verlet integration is for κ = 2 [12,13,33]. For higher order symplectic integration schemes, significantly larger values of κ are possible [34], though these schemes require multiple force evaluations per time step. We further assume that the shortest period of the nuclear degrees of freedom, T n (with the corresponding frequency Ω), is integrated in about 20 integration time steps such that δt = T n /20, which is a reasonable choice both in classical and regular Born-Oppenheimer molecular dynamics simulations. In this case we have that i.e. ω ≈ 4.5 × Ω. Thus, for a normal choice of integration time step as used in classical or regular Born-Oppenheimer molecular dynamics, we have a natural and system independent separation between the faster extended harmonic oscillator frequency ω and the slower nuclear frequency Ω. It is this natural adiabatic separation that motivates the derivation of the equations of motion in the adiabatic limit when ω 2 /Ω 2 → ∞. The ratio that governs the error, i.e. Ω 2 /ω 2 ≈ 1/20, may seem large to motivate the full adiabatic limit. However, as we will demonstrate in the examples below, this natural and material independent separation is sufficient to provide highly accurate molecular dynamics simulations.

III. EXAMPLES
We have performed an implementation of our generalized extended Lagrangian Born-Oppenheimer dynamics in the self-consistent charge density functional based tight-binding code LATTE [18,23]. It has a particularly simple form that allows a direct calculation of the exact non-local kernel. In density functional based tight-binding theory [35][36][37][38] the continuous charge density n(r) is represented by the atomic net Mulliken charges n = {n i } for each atom i. In the same way, the optimized ground state density ρ(r) is represented by net Mulliken charges q = {q i }. In this formalism the generalized electronic equations of motion in Eq. (16) are given bÿ where the kernel is calculated from with The derivatives δq j /δn k have been calculated using density matrix perturbation theory [39] and the nuclear degrees of freedom have been integrated using the velocity Verlet algorithm. For the electronic degrees of freedom we used a modified Verlet integration scheme [13,14,33], where (2), which is calculated at t = 0.5 ps and then kept constant for t > 0.5 ps, provides a very accurate dynamics, which is very close to an exact regular Born-Oppenheimer molecular dynamics.
This integration scheme is neither symplectic nor exactly time-reversible when α > 0. However, without the last dissipative force term in the integration scheme, accumulation of numerical noise may lead to instabilities.
The equations of motion in previous formulations [17,23,24] are recovered by a scaled δ-function approximation of the δq/δn derivative, where δq j /δn k ∼ δ jk . In this case the equations of motion are simplified tö As was previously discussed in Sec. II E, this corresponds to a harmonic oscillator evolving around a ground state density that is approximated by a linear mixing. As long as linear mixing gives an improved ground state estimate over n(r), this is a simple, efficient and stable approach. For many systems there is therefore no visible effect if we switch from the scaled δ-function approximation to the exact calculation of the non-local kernel in Eq. (2). However, for material systems where a simple constant linear mixing is not sufficient to reach a ground state solution the previous approach will fail. For such systems (or if the scaling constant c has been set too large or too small) instabilities occur in the time integration of n(r) and the simulations may diverge. In this case multiple self-consistent field iterations using some more advanced optimization scheme are necessary to improve the approximation of ρ(r), which significantly increases the computational cost. Figure 1 depicts some results of using various approximations to the kernel K(r, r ′ ) in Eq. (2) for simulations of liquid acrylonitrile. Only one Hamiltonian diagonalization (or density matrix construction) per time step was used. Liquid acrylonitrile (C 3 H 3 N) was chosen because it shows significant instabilities with a dynamics that is hard to stabilize with a simple scaled δ-function approximation of the kernel K(r, r ′ ). The two solid curves show what happens when we switch from a simple scaled δ-function approximation to the exact non-local expression in Eq. (2), which thereafter is kept constant after its calculation at t = 0.5 ps. The two dashed lines show the fluctuations in the total energy if we continue using the scaled δ-function approximation. As is clearly seen, the exact non-local kernel provides a very efficient stabilization of the dynamics even if it is kept fixed as a constant preconditioner after its initial calculation. In fact, recalculating the kernel for each new time step (not shown) makes very little difference  10), using a linear scaling density matrix expansion scheme [18,21].
to the simulation in Fig. 1. This indicates that it is possible to use more efficient, approximate calculations of the kernel in the first place, possibly Broyden based schemes [29][30][31], of Kerker mixing techniques [32], as well as other more recent methods [40,41] with similar results. These optimization schemes that are used in a broad range of self-consistent electronic structure methods are based on various non-local approximations of the exact kernel and reduces the computational cost significantly [27]. The upper panel of Fig. 2 shows the fluctuations in the total energy for a simulation of liquid methane. The effect on the total energy fluctuations from the exact non-local kernel, which is included at t = 0.25 ps, is fairly small. However, the error in the charge distribution (lower panel) shows a significant shift, but the error is small already at t < 0.25 ps when the local scaled δ-function approximation is used for the kernel. The δ-function approximation was scaled with an artificially small prefactor, c = 0.25, to enhance the effect of the switch to the non-local exact kernel. If a prefactor in the range of c ∈ [0.5, 0.7] had been chosen, the difference would have been much smaller.
The computational cost of calculating the exact kernel in Eq. (2) increases with the number of atoms such that a system with 100 atoms introduces an additional effort in a single calculation of K(r, r ′ ) corresponding to about 100 normal molecular dynamics time steps without the calculation of K(r, r ′ ). Fortunately, the calculation is trivially parallelizable and may thus add very little overhead in the total wall-clock time of a simulation, in particular, if the kernel only needs to be updated every few thousand time steps. However, for simulations of very large systems, which may require a more frequent update of the kernel, the full exact calculation of K(r, r ′ ) has to be avoided.
In Fig. 3 we demonstrate the scaling of the difference between ρ[n](r) and n(r) as a function of ω. The exact nonlocal kernel was calculated once in the beginning of the simulations and was thereafter kept constant. The root mean square error between ρ[n](r) and n(r) was calculated for 2,000 time steps. Since δt 2 ω 2 = κ, which occurs as a constant dimensionless factor in the modified Verlet integration, Eq. (22) or (27), we have that ω = √ κ/δt. Different values of ω are thus given by modifying the integration time step δt. The figure clearly demonstrates how the difference between ρ[n](r) and n(r) scales as ω −2 , which is consistent with the necessary assumption in the derivation of the equations of motion in the adiabatic limit as ω → ∞. This also illustrates the tunable accuracy of the simulation. Since ρ[n](r)− n(r) = F (n) in Eq. (11), the dynamics converges quadratically toward exact self-consistent Kohn-Sham Born-Oppenheimer molecular dynamics as ω → ∞. However, already for a normal integration time step when δt = 0.5 fs, the root mean square error is less than 0.005. The sensitivity to the integration time step and the value of ω can also be investigated for the vibrational spectra. Figure 4 shows the (shifted) vibrational density of states (DOS) of liquid Methane given from the Fourier transform of the velocity autocorrelation function [42]. After an initial equilibration of 50,000 time steps the velocity autocorrelation was sampled over 100,000 time steps every femtosecond for simulations using four different integration time steps δt. The lower curves show the results for Born-Oppenheimer molecular dynamics simulations using multiple self-consistent field iterations in the optimization of the ground state density before each force evaluation ("exact" BOMD). The four upper curves show the same results using our fast generalized extended Lagrangian formulation that requires only one diagonalization per MD time step with the optimized non-local kernel. The curves are virtually all on top of each other without any sensitivity to ω in for range of integration time steps. The longest integration time step, δt = 0.6 ps, is about 1/20 of the shortest period corresponding to the highest frequency of the system. For time steps δt = 0.75 ps we started to se an instability with a small but systematic drift in the energy both for the optimized Born-Oppenheimer simulation and the extended Lagrangian scheme with the non-local kernel.

IV. DISCUSSION
The generalized extended Lagrangian formulation of Born-Oppenheimer molecular dynamics allows a more accurate and stable dynamics for a broader class of systems. Calculations involving bond breaking or charge sloshing, which may occur in simulations of chemical reactions or large metallic systems, are often not possible to converge using a simple constant linear mixing of densities. The generalized non-local kernel acting on the electronic degrees of freedom can avoid such shortcomings and thus significantly increase the application range of extended Lagrangian Born-Oppenheimer molecular dynamics simulations in the fast limit that only requires one diagonalization per time step.
The derivation of the generalized equations of motion based on the adiabatic condition has several connections to the analysis of the equations of motion of extended Born-Oppenheimer molecular dynamics in its original formulation that recently was performed by Lin et al. [16] as well as the more generalized framework described by Hutter [15]. Their analysis also includes several interesting comparisons and relations to Car-Parrinello molecular dynamics [1,26,43], which is based on an extended Lagrangian approach to first principles molecular dynamics simulations that was introduced almost 30 years ago. The concept of a modified shadow Hamiltonian plays an important role in classical molecular dynamics [44][45][46][47]. Instead of an approximate numerical integration of an underlying exact Hamiltonian dynamics, a class of symplectic integration schemes corresponds to an exact numerical integration of an underlying approximate shadow Hamiltonian. In this way properties of the dynamical flow of a Hamiltonian dynamics can be rigorously conserved. Here we used the concept of a shadow Hamiltonian in the context of an incomplete self-consistent field optimization to describe and understand the accuracy in the long-term conservation of the total energy.
Although conceptually different and presented in a form that is easy to generalize to other functionals, the modified Born-Oppenheimer potential energy surface U(R, n) with the linearized electron-electron interaction in Eq. (7) has the same form as the Harris-Foulkes functional [48,49], with a second-order error term The main difference is that U(R, n) is represented in terms of a variationally optimized energy functional for a given external and electrostatic potential that are determined by the nuclear positions R and a dynamical variable density n(r). Since n(r) occurs as a dynamical variable in the extended Lagrangian, forces that are calculated from the Euler-Lagrange equations of motion are given by partial derivatives of U(R, n) with respect to a constant density n(r) in Eq. (16). This is not possible in a Harris-Foulkes functional, where n(r) represents either overlapping R-dependent atomic charge densities or the successive (non-converged and thus R-dependent) input densities in a regular iterative Kohn-Sham optimization procedure. Thus, our modified Born-Oppenheimer potential energy surface U(R, n) plays a quite different role and allows computationally simple and accurate calculations of the forces without relying on the Hellmann-Feynman theorem [1]. The important role of n(r) as a dynamical variable is sometimes misunderstood. It is only thanks to the underlying dynamical framework that our molecular dynamics simulations work. Extended Lagrangian Born-Oppenheimer molecular dynamics is sometimes referred to simply as a time-reversible extrapolation scheme for the electronic degrees of freedom from one step to the next. That picture is misleading in the same way as if we would view extended Lagrangian Car-Parrinello molecular dynamics simply as an ad hoc numerical interpolation technique for wavefunctions.
The calculations of the exact non-local kernel K(r, r ′ ) increases the computational cost of extended Lagrangian Born-Oppenheimer molecular dynamics. However, for many systems a local scaled δ-function approximation is sufficient. Other approximate expressions based on, for example, Broydens method or Kerker mixing could also be applied. In fact, it is straightforward to invent numerous flavors of local or semi-local approximations of the Kernel that require different levels of computational complexity. However, each update of an approximate Kernel, for example based on the Broyden scheme, would introduce a computational overhead similar to a full self-consistent field optimization. An understanding of the efficiency of various existing non-local kernel approximations can be estimated from their previous applications to accelerate the self-consistent field optimization in regular Kohn-Sham density functional theory with a fixed external potential [28-31, 40, 41]. However, there is exist no perfect, fully automatic, black box solution to the self-consistent field optimization problem and the same limitation can therefore be expected also for application in the generalized form of extended Lagrangian Born-Oppenheimer molecular dynamics.

V. SUMMARY AND CONCLUSIONS
We have proposed a generalized framework for extended Lagrangian Born-Oppenheimer molecular dynamics, Eqs.
(1) and (16), that improves the stability and extends the range of applications in the limit of fast quantum based molecular dynamics simulations that require only one diagonalization per time step. We showed how the equations of motion, Eq. (16), can be derived directly from the Lagrangian, Eq. (1), under the condition of an adiabatic separation between the nuclear and electronic degrees of freedom. We also showed how this adiabatic separation is automatically fulfilled for normal choices of the integration time step.
We believe that the general ideas presented in this paper can be directly applied to a broad class of Born-Oppenheimer-like molecular dynamics simulation techniques based on, for example, equilibrated or constrained charge relaxation for reactive or polarizable force field calculations, orbital-free density functional theory, as well as density matrix, Green's or correlated wavefunction based methods beyond the effective Kohn-Sham single particle formulation of density functional theory.

VII. APPENDIX
Previous derivations of the equations of motion in the limit of vanishing self-consistent field convergence followed a different path [17,24] from what will be given here. In earlier derivations a set of Euler-Lagrange equations of motion for extended Lagrangian Born-Oppenheimer molecular dynamics were first derived in the limit of full self-consistent field convergence as µ → 0. Then, a new set of equations of motion were derived by relaxing the condition of the self-consistent field optimization. That approach is not fully transparent. Below we will make a derivation of the equations of motion directly from the generalized extended Lagrangian without the intermediate step of assuming a full self-consistent field convergence. Instead, we will assume an adiabatic separation between the electronic and nuclear degrees of freedom and derive the equations of motion in the limit as ω/Ω → ∞ and µ → 0. This approach is motivated by the natural adiabatic separation between the nuclear and the electronic degrees of freedom that was analyzed in Sec. II G.