{'fraction_non_alphanumeric': 0.03776788851153531, 'fraction_numerical': 0.014111722684155054, 'mean_word_length': 5.123700250986016, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Extended Lagrangian Born-Oppenheimer molecular dynamics in the limit of vanishing self-consistent field optimization

We present an efficient general approach to first principles molecular dynamics simulations based on extended Lagrangian Born-Oppenheimer molecular dynamics in the limit of vanishing self-consistent field optimization. The reduction of the optimization requirement reduces the computational cost to a minimum, but without causing any significant loss of accuracy or longterm energy drift. The optimization-free first principles molecular dynamics requires only one single diagonalization per time step and yields trajectories at the same level of accuracy as"exact", fully converged, Born-Oppenheimer molecular dynamics simulations. The optimization-free limit of extended Lagrangian Born-Oppenheimer molecular dynamics therefore represents an ideal starting point for a robust and efficient formulation of a new generation first principles quantum mechanical molecular dynamics simulation schemes.


I. INTRODUCTION
With the rapid growth of available processing power, first principles molecular dynamics simulations, where the forces acting on the atoms are calculated on the fly using a quantum mechanical description of the electronic structure, are becoming an increasingly powerful tool in materials science, chemistry and biology [1]. While some early applications where performed already four decades ago [2,3], it was not until the development of efficient plane-wave pseudopotential methods [4][5][6][7][8][9] based on density functional theory [10,11] and the fast Fourier transform [12], that first principles molecular dynamics simulations became broadly applicable.
There are two major approaches to first principles molecular dynamics: a) Born-Oppenheimer molecular dynamics [1][2][3][7][8][9] and b) extended Lagrangian Car-Parrinello molecular dynamics [1,4,5,7,13,14,16,17,[17][18][19]. In Born-Oppenheimer molecular dynamics, the forces acting on the atoms are calculated at the relaxed electronic ground state in each time step, which provides a well defined and often very accurate approximation. A key problem, however, is that a straightforward implementation of Born-Oppenheimer molecular dynamics is unstable and does not conserve energy without a high degree of convergence in the electronic structure calculations. If this is not achieved, the electronic system behaves like a heat sink or source, gradually draining or adding energy to the atomic system [5,20]. Several techniques have therefore been developed that attempts to improve the efficiency of Born-Oppenheimer molecular dynamics and reduce the computational cost of the electronic optimization procedure [20][21][22][23]. In extended Lagrangian Car-Parrinello molecular dynamics, on the other hand, the computationally expensive ground state optimization is avoided. As in Ehrenfest based molecular dynamics [36][37][38], the electrons are instead treated as separate dynamical variables oscillating around the ground state. This approach permits a stable dynamics with a low computational cost per time step. Unfortunately, Car-Parrinello molecular dynamics simulations typically require shorter integration time steps and a system-dependent choice of electron mass parameters to yield reliable results in comparison on an "exact" Born-Oppnheimer molecular dynamics; although statistical averages are often in good agreement [1,19].
Recently, an extended Lagrangian formulation for a time-reversible Born-Oppenheimer molecular dynamics was proposed [24,25], which combines some of the best features of Car-Parrinello and regular Born-Oppenheimer molecular dynamics, while avoiding some of their most serious shortcomings. It has been argued that extended Lagrangian Born-Oppenheimer molecular dynamics can be seen as a general framework both for Born-Oppenheimer and Car-Parrinello molecular dynamics [19]. In this modern formalism of extended Lagrangian first principles molecular dynamics, Car-Parrinello molecular dynamics appears in the limit of vanishing self-consistent field optimization [19]. However, the optimization-free limit can be approached in different ways providing a variety of solutions. In this paper we show how extended Lagrangian Born-Oppenheimer molecular dynamics, in the limit of vanishing self-consistent field optimization, gives a first principles molecular dynamics at the same level of accuracy as "exact" Born-Oppenheimer molecular dynamics, but without requiring short integration time steps or a material dependent tuning of electron mass parameters as in Car-Parrinello molecular dynamics. The instability from the systematic energy drift associated with incomplete convergence of the electronic structure in regular Born-Oppenheimer molecular dynamics is also avoided. Our work here represents a generalization and first principles extension of recent work that was demonstrated for semiempirical self-consistent-charge tight-binding simulations [26].
The ability to achieve a high degree of accuracy in the limit of vanishing self-consistent field optimization serves two main purposes: 1) it simplifies the calculations with a reduction of the optimization cost to a minimum, and 2) it provides the ideal starting point for fully converged, i.e. "exact", Born-Oppenheimer molecular dynamics simulations when the requirement of accuracy is very high. The optimization-free limit of extended Lagrangian Born-Oppenheimer molecular dynamics therefore represents an efficient and robust framework for a new generation of first principles molecular dynamics simulations.

II. EXTENDED LAGRANGIAN BORN-OPPENHEIMER MOLECULAR DYNAMICS
Extended Lagrangian Born-Oppenheimer molecular dynamics [24] can be formulated in terms of a Lagrangian, where the regular Born-Oppenheimer Lagrangian defined at the electronic ground state density matrix D for a given set of nuclear coordinated, {R I } = R, has been extended with auxiliary dynamical variables for the electronic degrees of freedom, P 0 andṖ 0 , that evolve in a harmonic well centered around D. The potential energy U (R; D) is here the Hartree-Fock or Kohn-Sham energy functional including the ion-ion repulsion energy [27]. The parameter µ is a fictitious electron mass and ω is the frequency determining the curvature of the harmonic well. Euler-Lagrange equations, in the limit µ → 0 [24], gives the decoupled equations of motion: The partial derivative of U for the nuclear coordinate R I is taken with respect to a constant P 0 , since P 0 is an independent dynamical variable. The equations of motion can be integrated using a time-reversible symplectic scheme, both for the nuclear and electronic degrees of freedom [39,40]. By using a time-reversible P 0 as the initial guess of the iterative self-consistent field (SCF) optimization procedure, where the total Born-Oppenheimer energy, is stable without any long-term energy drift, even in the case of approximate convergence of P n [24][25][26]28]. The ground state density matrix D in Eq. (4) is given from the Heaviside step function, θ, of the converged Fockian or Kohn-Sham Hamiltonian, i. e. for lim n→∞ H(P n ), in an orthogonal representation, Z T HZ, with the step formed at the chemical potential, µ 0 , separating the occupied from the unoccupied states. The congruence transformation matrix Z is given from the inverse Cholesky or Löwdin factorization of the overlap matrix, S, determined by ZSZ T = I.

A. Fast quantum mechanical molecular dynamics
As n → 0 in Eq. (4), i.e. in the limit of vanishing self-consistent field optimization, the equations of motion for the extended Lagrangian formulation of Born-Oppenheimer molecular dynamics, Eq. (2), are given by By avoiding the self-consistent-field optimization of P 0 , these equations of motion require only one single diagonalization per time step in the construction of D(P 0 ) and therefore provide a computationally fast method for first principles quantum mechanical molecular dynamics (fast-QMMD) [26]. An alternative derivation of the fast dynamics represented by Eq. (6) that is motivated through a different set of arguments is given in Ref. [26].
To guarantee stability in the integration of the electronic degrees of freedom in Eq. (6), using an integration time step of δt, the dimensionless integration parameter δt 2 ω 2 , typically needs to be rescaled by a factor c ∈ [0, 1] [26] compared to the original integration of extended Lagrangian Born-Oppenheimer molecular dynamics [30]. This stability condition further assumes convexity of the total energy functional between P 0 and D(P 0 ) [26].
The definition of D ≡ lim n→∞ D(P n ) in Eq. (4) and our particular choice of sequence of limits both for µ → 0 and n → 0 are important. For example, if we instead use D ≡ P n and let n → 0 in the Lagrangian (before deriving the Euler-Lagrange equations of motion), we end up with a µ-dependent set of unconstrained Car-Parrinellolike equations [19] and if µ → 0 already in the initial Lagrangian, but with full self-consistency convergence, we recover (trivially) regular Born-Oppenheimer molecular dynamics. For our particular sequence of limits of µ and n, the fast-QMMD defined by Eq. (6) is formally neither an extended Lagrangian nor a Born-Oppenheimer molecular dynamics. However, as will be demonstrated in our examples, the first principles fast-QMMD in Eq. (6) is a very close approximation of "exact", fully converged, extended Lagrangian Born-Oppenheimer molecular dynamics.

A. Implementation
Our fast-QMMD, Eq. (6), has been implemented based on Hartree-Fock theory in the Uppsala Quantum Chemistry (UQuantChem) simulations package [29], which is a freely available suite of programs for parallel ab initio electronic structure calculations using Gaussian basis sets, including Hartree-Fock and Møller-Plesset perturbation theory, configuration interaction, variational and diffusion Monte-Carlo, structural optimization, and first principles molecular dynamics. The nuclear coordinates are integrated using the velocity Verlet scheme and the electronic degrees of freedom with a modified Verlet algorithm, including a weak dissipation term to remove the accumulation of numerical noise [25,30].
Since P 0 appears as a dynamical variable in Eq. (6), a Hellmann-Feynman-like expression for the nuclear forces, under the constraint of P 0 being constant, can be applied. Thus, even if the ground state condition necessary for Hellmann-Feynman forces are not fulfilled, we still have a force expression of similar simplicity. For the basis-set dependent contribution we use the original expression of the Pulay force term [31], which provides a sufficiently accurate approximation [26]. Our first principles dynamics is implemented based on Hartree-Fock theory [27,32]. The Hartree-Fock method is the starting point for correlated wavefunction methods and can be used as the computational prototype for density functional theory [10,11,33,34] and hybrid schemes [35]. Our optimization-free Hartree-Fock molecular dynamics therefore demonstrates applicability for a broad class of first principles methods. Extensions to plane wave schemes should also be straightforward [25].  a STO-3G basis set in (a) and a 6-31G * * in (b). Regular Born-Oppenheimer molecular dynamics, where the density matrix from the previous time step was used as the initial guess to the iterative ground-state optimization, exhibits an unphysical systematic drift in the total energy because of the broken time-reversal symmetry [20,24]. This drift is avoided in the "exact" fully optimized extended Lagrangian Born-Oppenheimer molecular dynamics (XL-BOMD), Eq. (2), which is very close   2), for the nuclear forces (red squares) and for the total energy (green circles) calculated for four different molecules at different time steps. For comparison the local error of the total energy (black filled circles) has been calculated. Simulations were performed with a 6-31G * * basis set, using URHF theory as implemented in the UQuantChem code [29]. , after perturbing the density matrix, P0(t), at t = 500 fs, by resetting P0 at t = 500 fs to the initial density matrix at t = 0.5 fs. Simulations were performed for a single water molecule at room temperature with a time step of 0.5 fs using URHF theory as implemented in the UQuantChem code [29].

B. Molecular dynamics simulations
to the results from the optimization-free fast first principles QMMD (red circles), Eq. (6), as seen in the insets.
In particular, any deviations between the optimizationfree and the fully optimized Born-Oppenheimer molecular dynamics simulations are small compared to the local truncation error, i.e. the amplitude of the total energy fluctuations that are caused by the finite size of the inte-gration time step δt = 10 a. u. As in classical molecular dynamics, the dominating integration error is thus determined by our choice of integration scheme and the size of the time step. Figure 2 demonstrates the longterm stability of our first principles fast-QMMD, which shows no systematic drift in the energy over 120 ps of simulation time. However, for longer integration time steps we have occasionally noticed a small drift that seems to be caused by the dissipation force of the modified Verlet scheme. This sensitivity, which not yet is completely understood, is not found in partially or fully SCF optimized versions of extended Lagrangian Born-Oppenheimer molecular dynamics. The next figure, Fig.  3, shows the corresponding simulation for a water cluster containing 10 water molecules simulated for a shorter simulation time. Because of the chaotic movements of the larger system, a direct comparison with respect to the total energy is harder. The inset shows a comparison of the kinetic energy fluctuations given by the temperature over the first 200 fs of simulation time, shortly before they eventually get out of phase. The total energy fluctuations of the fast-QMMD simulation shows a noisy behavior similar to a random walk compared to the "exact" Born-Oppenheimer simulation (XL-BOMD). Similar random walk-like noise have been seen in linear scaling XL-BOMD simulations [28]. Figure 4 shows the behavior of the total energy, Eq. (5), in simulations of a C 2 H 6 molecule using a 6-31G * * basis set, which represents a slightly larger and more complex system compared to the water molecule. As a comparison two different time steps were used, in panel (a) ∆t = 5 au and in panel (b) ∆t = 10 au.
In Figure 5, simulations of a H 2 O and a CF 4 molecule, showing the behavior of the total energy, Eq. (5), using two times respectively four times as long time steps as the maximum time step used in the previous examples. Figure 6 illustrates the behavior of the interatomic distances, in simulations of a H 2 O and a CF 4 molecule. Figure 7 shows the convergence toward "exact" Born-Oppenheimer molecular dynamics as the length of the integration time step δt is reduced. This scaling demonstrates how the fast-QMMD scheme provides a well defined approximation to exact Born-Oppenheimer molecular dynamics whith an error of order δt 2 , i.e.
The corresponding behavior was recently found in our studies based on self-consistent-charge tight-binding simulations [26].
To illustrate the stability of first principles fast-QMMD we perturb a simulation by resetting the auxiliary density matrix P 0 (t) to its t 0 initial value after 500 fs of simulation time. During the continued simulation, the perturbation slowly disappears as P 0 (t) converges toward the electronic ground state, as seen in Figure 8, where the deviation of the fast-QMMD density matrix P 0 relative to the "exact" density matrix P 5 is plotted as a function of time. This behavior demonstrates a key mechanism of our method. Instead of optimizing to ground state in each iteration as in regular Born-Oppenheimer molecular dynamics, the time evolution of the electronic degrees of freedom makes P 0 (t) converge toward the ground state dynamically. At convergence, the auxiliary density matrix P 0 (t) oscillates around the exact ground state with an amplitude that is of the order δt 2 .

IV. CONCLUSIONS AND SUMMARY
The extended Lagrangian approach to first principles molecular dynamics, as pioneered by Roberto Car and Michele Parrinello [4], in its modern formulation of extended Lagrangian Born-Oppenheimer molecular dynamics [19,24], provides an efficient and versatile framework for first principles molecular dynamics simulations.
Here we have shown how the ground state optimization requirement can be simplified and reduced to a mini-mum without causing any significant loss of accuracy or long-term stability. This has been demonstrated using Hartree-Fock theory and should be applicable to a broad class of first principles methods. The optimization-free first principles molecular dynamics requires only one single diagonalization per time step and yields trajectories that are very close to an "exact", time-reversible, first principles Born-Oppenheimer molecular dynamics simulation.