{'fraction_non_alphanumeric': 0.03218500417213017, 'fraction_numerical': 0.011681964477291692, 'mean_word_length': 4.821998612074948, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, '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}}
Wavefunction extended Lagrangian Born-Oppenheimer molecular dynamics

Extended Lagrangian Born-Oppenheimer molecular dynamics [Niklasson, Phys. Rev. Lett. 100 123004 (2008)] has been generalized to the propagation of the electronic wavefunctions. The technique allows highly efficient first principles molecular dynamics simulations using plane wave pseudopotential electronic structure methods that are stable and energy conserving also under incomplete and approximate self-consistency convergence. An implementation of the method within the planewave basis set is presented and the accuracy and efficiency is demonstrated both for semi-conductor and metallic materials.


I. INTRODUCTION
As the available computational capacity for scientific computing is growing, first principles Born-Oppenheimer (BO) molecular dynamics (MD) is becoming an increasingly important tool for studying a wide range of material problems. First principles BOMD delivers a very accurate approach to atomistic simulations without relying on a fitted parameterization of the atomic interactions as in classical molecular dynamics. Unfortunately, applications of BOMD that are based on self-consistent field (SCF) calculations such as density functional theory [1][2][3] are often limited by a very high computational cost or by fundamental shortcomings such as unbalanced phase space trajectories, numerical instabilities and a systematic long-term energy drift [4][5][6] .
Recently an extended Lagrangian BOMD (XL-BOMD) was introduced that avoids some of the most serious problems of regular BOMD and enables computationally efficient and stable simulations of energy conserving (microcanoncial) ensembles [6][7][8][9] . In XL-BOMD, auxiliary electronic degrees of freedom are included, in addition to the nuclear coordinates and velocities. In contrast to the popular extended Lagrangian Car-Parrinello molecular dynamics methods 3,10-12 , the nuclear forces are calculated at the ground state BO potential energy surface and the total BO energy is a constant of motion.
So far XL-BOMD has been limited to density matrix formulations of the extended electronic degrees of freedom. This excludes any practical implementation in widely used plane-wave pseudopotential schemes, since it would lead to unmanageable large density matrices. Because of the arbitrary phase of the electronic wavefunctions [13][14][15] it is difficult to use wavefunctions as the extended electronic degrees of freedom in a stable time-reversible or geometric integration of the equations of motion. Here we show how the electronic wavefunctions can be included in XL-BOMD. Our formulation allows a time-reversible integration of both the nuclear and the electronic degrees of freedom and it provides a highly efficient BOMD for plane-wave pseudopotential methods that is stable and energy conserving also under incomplete and approximate SCF convergence. The wavefunction XL-BOMD method was implemented in the Vienna Ab-initio Simulation Package (VASP) [16][17][18] and its accuracy and efficiency are demonstrated both for semiconductor and metallic materials.

II. FIRST PRINCIPLE MOLECULAR DYNAMICS
A. Born-Oppenheimer molecular dynamics First principle BOMD based on density functional theory (DFT) is given by the Lagrangian, where R = {R i } are the nuclear coordinates and the dot denotes the time derivative. The potential U DFT [R; Ψ sc ] is the ground state energy, including ion-ion repulsions, for the density given by the self consistent (sc) electronic wavefunctions, Ψ sc = {ψ sc nk }. Here n and k denote the band and reciprocal lattice vectors, respectively. The Euler-Lagrange equations, give the equations of motion for the dynamical variables R(t) andṘ(t).
The high cost of finding the ground state SCF solution Ψ sc (t) is significantly reduced by using an initial guess that is extrapolated from previous time steps 4,5,[19][20][21][22] , In the SCF optimization in Eq. (3) above we assume a full optimization, which may include several iterative cycles based on, for example, simple linear mixing, Broyden mixing, or the direct inversion of the iterative subspace (DIIS) method [23][24][25][26] . However, since the SCF optimization in practice never is complete, the extrapolation procedure in Eq. (3) leads to an irreversible evolution of the ground state electronic wavefunctions. The nuclear forces are therefore calculated with an underlying electronic degrees of freedom that behave unphysically. The irreversibility appears most strikingly as a systematic long-term energy drift 4,5 . By using thermostats, e.g. an artificial interaction with an external heat bath, these shortcomings of BOMD may not be noticed. However, a thermostat requires an underlying dynamics that is physically correct, and the problems are therefore never removed. Only by improving the SCF convergence, which is increasing the computational cost, is it possible to suppress the energy drift, though the problem never fully disappears.

B. Wavefunction Extended-Lagrangian Born-Oppenheimer MD
In our wavefunction XL-BOMD, proposed here, the dynamical variables of the BO Lagrangian are extended with a set of auxiliary wavefunctions Φ = {φ nk } evolving in harmonic oscillators centered around the selfconsistent ground state wavefunctions Ψ sc (t), Here µ is a fictitious electron mass parameter and ω is a frequency or curvature parameter for the harmonic potentials. Applying the Euler-Lagrange equations to the extended Lagrangian in Eq. (4) gives In the limit µ → 0 we get Thus, in the limit of vanishing fictitious mass parameter, µ, we recover the regular BO equations of motion in Eq. (7), with the total BO energy as a constant of motion. Equation (8) determines the dynamics of our auxiliary wavefunctions Φ(t). Since µ is set to zero, the only remaining undetermined parameter is the frequency or curvature ω of the extended harmonic potentials. As will be shown below, ω occurs in the integration of Eq. (8) only as a dimensionless factor δt 2 ω 2 and therefore affects the dynamics in the same way as the finite integration time step δt.
Since the auxiliary wavefunctions Φ(t) are dynamical variables, they can be integrated by, for example, the time-reversible Verlet algorithm 27 . Morevover, since the auxiliary wavefunctions evolve in a harmonic well centered around the ground state solution, Φ(t) will stay close Ψ sc (t). By maximizing the curvature ω 2 of the harmonic extensions we can minimize their separation. Using the auxiliary dynamical variables Φ(t) in the initial guess to the SCF optimization, therefore provides an efficient SCF procedure that can be used within a time-reversible framework. The nuclear forces will then be calculated with an underlying electronic degrees of freedom with the correct physical time-reversal symmetry. This is in contrast to conventional BOMD, where the SCF optimization is given from an irreversible propagation of the underlying electronic degrees of freedom as in Eq. (3). Hence the system will be propagated reversibly and should not suffer from any systematic drift in the total energy and phase space.

C. Integration
Both the nuclear and electronic degrees of freedom in Eqs. (7) and (8) can be integrated with the Verlet algorithm, or with other geometric integration schemes that preserve properties of the exact underlying flow of the dynamics 7,9,28 . The Verlet integration of Eq. (8), including a weak external dissipative electronic force that removes accumulation of numerical noise 8 , has the following form where α determines the magnitude of the dissipative force term with the c m coefficients given in Ref. 8 . The additional electronic force introduces dissipation of numerical noise that would accumulate in a perfectly reversible and lossless propagation. The dissipation breaks timereversibility, but only to a high order in δt 8,29 . In this way numerical errors can be removed without causing any significant drift in the total energy. The electronic ground state wavefunctions are unique except with respect to their phase. This presents a problem for the accuracy and stability of the Verlet integration above. The wavefunctions need to be aligned to a common orientation to allow an accurate and stable integration. Aligning the wavefunctions backwards in time as in some previous integration schemes for regular BOMD [13][14][15]30 is not possible, since it would break the time-reversal symmetry. We solve this problem by including a unitary rotation transform U in the SCF optimization, which rotates Ψ sc (t) such that the deviation from Φ(t) is minimized in the Frobenius norm, i.e.
U can be calculated from U = (OO † ) −1/2 O where O = Ψ sc |Φ is the overlap matrix between Ψ sc (t) and Φ(t) 13,14,30 . Since the rotation is only applied to Ψ sc (t) and not to previous auxiliary wavefunctions, the reversibility is not affected. The redefined Verlet integration is Note that good initial values for the axillary variables are important. If a poor initial guess are used the weak dissipation will eventually relax the auxiliary dynamics to a similar dynamics, but it would take time, and meanwhile we would have bad initial guesses for the SCF optimization. In our implementation the initial values of the auxiliary variables, are set to the SCF optimized ground states, i.e. as intial conditions for Φ(t) we chose to set Φ(t) ≡ Ψ sc (t) for the first K + 1 time steps, where we perform phase alignements to the first optimized wavefunctions Ψ sc (t = t 0 ). It may be preferable to run with a stronger convergence criteria during the initial steps to get a good starting guess. In the first K + 1 initial steps we therefore chose to have a higher degree of SCF convergence than in later time steps.

III. STABILITY AND NOISE DISSIPATION
By aligning the phase in the SCF optimization, the stability of the Verlet integration in Eq. (12), under the condition of an approximate and incomplete SCF convergence, can be analyzed from the roots λ of the characteristic equation of the homogeneous (steady state) part of the Verlet scheme, in the same way as for the density matrix 7,8 . Assume a linearization of an approximate SCF optimization, Eq. (9), around the hypothetical exact solution Ψ * , where Let γ be the largest eigenvalue of the SCF response kernel Γ SCF . Inserting Eq. (13) in the Verlet scheme, Eq. (12), with Γ SCF replaced by γ, for the homogeneous steady state solution for which Ψ * ≡ 0, gives the characteristic equation Here the dimensionless constant κ = δt 2 ω 2 , and γ ∈ [−1, 1] is proportional to the amount of convergence in the SCF optimization. As long as the initial guess Φ(t) is brought closer to the ground state solution by the SCF procedure, |γ| will be smaller than 1. If the characteristic roots have a magnitude |λ| max > 1, the integration is unstable (even if the accuracy is good), whereas it is stable if |λ| max ≤ 1 (even if the optimization is approximate). For |λ| max < 1 the accumulation of numerical noise will be suppressed through dissipation. By optimizing κ = δt 2 ω 2 under the condition of stability under incomplete SCF convergence with γ ∈ [−1, 1], the curvature ω 2 of the extended harmonic wells will be maximized, which keeps the auxiliary wavefunctions Φ(t) as close as possible to the ground state solutions Ψ sc (t) 9 . This optimization is performed under the additional condition of maximum dissipation. Our optimized values of α and κ and the c m coefficients can be found in Ref. 8

IV. PLANE WAVE PSEUDO-POTENTIAL IMPLEMENTATION (VASP)
Our wavefunction XL-BOMD method has been implemented in the Vienna Ab-initio Simulation Package (VASP) [16][17][18] and the projector augmented wave method 31,32 . These particular methods not only require the integration of the wavefunctions, but also the electron density and the Kohn-Sham eigenvalues that are used in the SCF optimization. In this work these additional quantities has been be added to the Lagrangian as extended dynamical variables evolving in harmonic oscillators centered around their own optimized values in the same way as the auxiliary wavefunctions. For example, an auxiliary density, ρ(r), can be included as a dynamical variable through the extended Lagrangian, Here ρ(r) follows the SCF optimized ground state density n sc (r). Also for the auxiliary density and other extended variables the initial values are set equal to the optimized ground state values for the first K + 1 steps.

V. APPLICATIONS
To demonstrate the accuracy and efficiency of the wavefunction XL-BOMD scheme we simulate two different model systems with qualitatively different bonding, metallic sodium and semiconducting silicon.

A. Sodium
A unit cell of 16 bcc Na atoms was simulated for a total of 10,000 steps with an ionic temperature fluctuating around 500 K with a time step of 4 fs. The regular BOMD integration scheme was based on a 2nd-order extrapolation of the wavefunctions from three previous time steps 13,14 and the XL-BOMD scheme used K = 5 for the dissipation 8 . Both methods used the velocity Verlet integration for the nuclear degrees of freedom and were run with the same SCF energy convergence criterion, δE = 5 µeV, resulting in about 2 SCF iterations per time step for both methods. Each SCF cycle includes one single construction and solution of the Hamiltonian eigenvalue problem. As SCF convergence accelerating algorithm we used the DIIS scheme 26 . A plane-wave energy cutoff of 102 eV and a grid of 64 k-points was used and the exchange-correlation energy was given by the local density approximation (LDA) 33 .
The fluctuations in the total energy can be seen in Fig.  2. For regular BOMD we see a small but systematic drift in the total energy of the order of 0.25 meV/ps. In comparison, XL-BOMD shows no drift and the magnitude of the energy fluctuations due to the local truncation errors, occurring because of the finite time steps and the approximate SCF convergence, is the same. In fact, we have found that XL-BOMD is stable even when only 1 SCF cycle per time step is used. This would be a general statement if the SCF procedure systematically improves the convergence in a single step 7 . Unfortunately, this is not always the case.

B. Silicon
Next a Si system with 8 atoms per unit cell in a diamond structure was simulated for a total of 10,000 steps at an ionic temperature fluctuating around 500 K with a time step of 1 fs. The SCF convergence threshold δE was set to 5 µeV with a plane wave cutoff of 246 eV and a grid of 64 k-points was used. Otherwise the same settings as for sodium were applied. In Fig. 3 the fluctuations in the total energy ∆E is plotted. Also in this case, we find that XL-BOMD restores balance to the unphysical trajectories of regular BOMD that shows a significant systematic drift in the total energy.

VI. DISCUSSION AND SUMMARY
In a direct comparison using the same time step and convergence criteria, we find that XL-BOMD and regular BOMD have the same local truncation error, as measured by the local amplitude of the oscillations in the total energy, both for the metallic and the non-metallic system. However, the unphysical behavior of regular BOMD, which has a systematic long-term energy drift, is removed in XL-BOMD. Only by significantly increasing the computational cost of regular BOMD with a higher degree of SCF convergence, or shorter time steps, is it possible to reduce the long-term energy drift. Table  II summarizes the results of a comparison between XL-BOMD and regular BOMD for the Na simulation. The results clearly show that even though the energy drift can be substantially reduced in conventional BOMD, a large performance penalty has to be paid. XL-BOMD requires in general more memory, 1.5 to 2.5 times the temporary storage used in regular BOMD depending on the dissipation scheme used. However, for most practical situations wall time is the limiting factor when running first principal BOMD, not memory usage. XL-BOMD therefore combines a more correct physical description with a lower computational cost.
In many ways, XL-BOMD integrates some of the best features of regular BOMD and Car-Parrinello molecular dynamics, i.e the parameter-free rigor of BOMD and an efficient extended Lagrangian framework as in Car-Parrinello molecular dynamics, where both nuclear and electronic degrees of freedom are included as dynamical variables.
In summary, we have proposed and demonstrated a wavefunction XL-BOMD scheme that allows highly efficient first principles molecular dynamics simulations using plane wave pseudopotential electronic structure methods that are stable and energy conserving also under incomplete and approximate self-consistency convergence. This extends the capability and accuracy of mod-