Analysis of the Time Reversible Born-Oppenheimer Molecular Dynamics

We analyze the time reversible Born-Oppenheimer molecular dynamics (TRBOMD) scheme, which preserves the time reversibility of the Born-Oppenheimer molecular dynamics even with non-convergent self-consistent field iteration. In the linear response regime, we derive the stability condition as well as the accuracy of TRBOMD for computing physical properties such as the phonon frequency obtained from the molecular dynamic simulation. We connect and compare TRBOMD with the Car-Parrinello molecular dynamics in terms of accuracy and stability. We further discuss the accuracy of TRBOMD beyond the linear response regime for non-equilibrium dynamics of nuclei. Our results are demonstrated through numerical experiments using a simplified one dimensional model for Kohn-Sham density functional theory.

2. An illustrative model 70 To start, let us illustrate the main idea for a simple model problem, which provides the essence of TRBOMD in a much simplified setting. Consider the following nonlinear ODË where we assume that the right hand side f (x) is difficult to compute, and it can be approximated by an iterative procedure. Starting from an initial guess s ≈ f (x), the final approximation via the iterative procedure is denoted by g(x, s). We assume the approximation g(x, s) is consistent, i.e.
g(x, f (x)) = f (x). ( To numerically solve the ODE (1), we discretize it by some numerical scheme, then it remains to decide the initial guess s at each time step. A natural choice of s would be g(x, s) from the previous step, as x does not change much in successive steps. For instance, if the Verlet algorithm is used and t k = k∆t with ∆t being the time step, the discretized ODE becomes x k+1 = 2x k − x k−1 + (∆t) 2 g(x k , s k ), We immediately observe that the discretization scheme (3) breaks the time reversibility of the original ODE (1). In other words, for the original ODE (1), we propagate the system forward in time from (x(t 0 ),ẋ(t 0 )) to (x(t 1 ),ẋ(t 1 )). Then if we use (x(t 1 ),ẋ(t 1 )) as the initial data at t = t 1 and propagate the system backward in time to time t = t 0 , we will be at the state (x(t 0 ),ẋ(t 0 )). The loss of the time reversible structure can introduce large error in long time numerical simulation [19]. This is the main reason why BOMD with non-convergent SCF iteration fails for long time simulations [13]. To overcome this obstacle, the idea of TRBOMD is to introduce a fictitious dynamics for the initial guess s. Namely, we consider the time reversible coupled system x(t) = g(x(t), s(t)), s(t) = ω 2 (g(x(t), s(t)) − s(t)), where ω is an artificial frequency. We analyze now the accuracy and stability of Eq. (4) in the linear response regime by assuming that the trajectory x(t) oscillates around a equilibrium position x * . We denote by x(t) = x(t) − x * the deviation from the equilibrium position and s(t) = s(t) − f (x(t)) the deviation of the initial guess from the exact force term. Consequently, the equation of motion (4) can be rewritten as (for simplicity we suppress the t-dependence in the notation for the rest of the section) x = g(x, s), where the term −f ′′ (x)(ẋ) 2 − f ′ (x)ẍ comes from the term f (x) in s by the chain rule.

71
In the linear response regime, we assume the linear approximation of force for x around x * : where Ω is the oscillation frequency of x in the linear response regime. We also linearize g with respect to s and x and dropping all higher order terms as g(x, s) = g(x, f (x) + s) ≈ g(x, f (x)) + g s (x, f (x)) s ≈ −Ω 2 x + g s (x * , f (x * )) s, where g s denotes the partial derivative of g with respect to s and the consistency condition (2) is applied. We then have In accord with notations used in later discussions, let us denote with which the linearized system of Eq. (5) becomes Note that when the force is computed accurately, i.e.
g(x, s) = f (x), ∀s, we have L = 0, K = 1, meaning that the motion of x is decoupled from that of s, and x follows the exact harmonic motion in the linear response regime with the accurate frequency Ω. When the force is computed inaccurately, x is coupled with s in Eq. (10). Actually, we can solve (10) analytically and the eigenvalues of A are Then the frequencies of the normal modes of the ODE are Ω = −λ Ω and ω = √ −λ ω respectively.
Assume ω 2 ≫ Ω 2 and expand the solution to the order of O(1/ω 2 ), we have Similarly the frequency for the other normal mode which is dominated by the motion of s is It is found that one of the normal mode of Eq. (10) has frequency Ω ≈ Ω. We can therefore measure the 72 accuracy of Eq. (4) using the relative error between Ω and Ω. Furthermore, if the dynamics (4) then K ≈ 1, and the accuracy of Ω is determined by L or g s (x * , f (x * )), which indicates the sensitivity 80 of the computed force with respect to the initial guess, or the accuracy of the iterative procedure for 81 computing the force. If a "good" iterative procedure is used, g s (x * , f (x * )) will be small. Therefore the 82 presence of the term L allows one to obtain relatively accurate approximation to the frequency Ω without 83 using a large ω. The same behavior can be observed when using TRBOMD to approximate BOMD (vide 84 post). 85 Finally, we remark that even though Eq. (1) is a much simplified system, it will be seen below that 86 for BOMD with M atoms and N interacting electrons, the analysis in the linear response regime follows 87 the same line, and the result for the frequency is similar to Eq. (14).

89
Consider a system with M atoms and N electrons. The position of the atoms at time t is denoted by R(t) = (R 1 (t), . . . , R M (t)) T . In BOMD, the motion of atoms follows Newton's law where E(R(t)) is the total energy of the system at the atomic configuration R(t). In KSDFT, the total energy is expressed as a functional of a set of Kohn-Sham orbitals To illustrate the idea with minimal technicality, let us consider for the moment a system of N electrons at zero temperature. The energy functional in KSDFT takes the form The first term in the energy functional is the kinetic energy of the electrons. The second term contains the electron-ion interaction energy. The ion-ion interaction energy usually takes the form I<J where Z I is the charge for the nucleus I. The ion-ion interaction energy does not depend on the electron density ρ. To simplify the notation, we include the ion-ion interaction energy in the V ion term as a constant shift that is independent of the x variable. The third term does not explicitly depend on the atomic configuration R, and is a nonlinear functional of the electron density ρ. It represents the Hartree part of electron-electron interaction energy (h), and the exchange-correlation energy (xc) characterizing many body effects. The energy E(R) as a function of atomic positions is given by the following minimization problem We denote by {ψ i (x; R)} N i=1 the (local) minimizer, and ρ * (x; R) = N i=1 |ψ i (x; R)| 2 the converged electron density corresponding to the minimizer (here we assume that the minimizing electron density is unique). Then the force acting on the atom I is In physics literature the force formula in Eq. (19) is referred to as the Hellmann-Feynman force. The 90 validity of the Hellmann-Feynman formula relies on that the electron density ρ * (x; R) corresponds to 91 the minimizers of the Kohn-Sham energy functional. Since E hxc [ρ] is a nonlinear functional of ρ, the 92 electron density ρ is usually determined through the self-consistent field (SCF) iteration as follows.
Starting from an inaccurate input electron density ρ in , one first computes the output electron density by solving the lowest N eigenfunctions of the problem and the output electron density ρ out is defined by Here the operator F is called the Kohn-Sham map. ρ out can be used directly as the input electron density ρ in in the next iteration. This is called the fixed point iteration. Unfortunately, in most electronic structure calculations, the fixed point iteration does not converge even when ρ in is very close to the true electron density ρ * . The fixed point iteration can be improved by the simple mixing method, which takes the linear combination of the electron density as the input density for the next iteration with 0 < α ≤ 1. Simple mixing can greatly improve the convergence properties of the SCF iteration over the fixed point iteration, but the convergence rate can still be slow in practice. . More detailed discussion on convergence properties of these SCF schemes can be found in [24]. In the following discussions, we denote by ρ SCF (x; R, ρ) the final electron density after the SCF iteration starting from an initial guess ρ. We assume that ρ SCF satisfies the consistency condition If a non-convergent SCF iteration procedure is used, ρ SCF (x; R, ρ) might deviate from ρ * (x; R). Such  The map ρ SCF is usually highly nonlinear, which makes it difficult to correct the error in the force. The TRBOMD scheme avoids the direct correction for the inaccurate ρ SCF , but allows the initial guess to dynamically evolve together with the motion of the atoms. We denote by ρ(x, t) the initial guess for the SCF iteration at time t. When ρ(·, t) is used as an argument, we also write ρ SCF (x; R(t), ρ(t)) := ρ SCF (x; R(t), ρ(·, t)). The Hellmann-Feynman formula (19) is used to compute the force at the electron density ρ SCF (x; R(t), ρ(t)) even though ρ * (x; R(t)) is not available. Thus, the equation of motion in TRBOMD reads It is clear that TRBOMD is time reversible. The discretized TRBOMD is still time reversible if the numerical scheme is time reversible. For instance, if the Verlet scheme is used, the discretized equation of motion becomes which is evidently time reversible. The artificial frequency ω controls the frequency of the fictitious 98 dynamics of ρ(x, t) and is generally chosen to be larger than the frequency of motion of the atoms. The

102
Let us mention that TRBOMD is closely related to CPMD. In CPMD, the equation of motion is given by where µ is the fictitious electron mass for the fake electron dynamics in CPMD, and Λ's are the Lagrange multipliers determined so that {ψ i (t)} is an orthonormal set of functions for any time. The CPMD scheme (27) can be viewed as the equation of motion with an extended Lagrangian which contains both ionic and electronic degrees of freedom. Therefore, CPMD is a Hamiltonian 103 dynamics and thus time reversible.

104
Note that the frequency of the evolution equation for {ψ i } in CPMD is adjusted by the fictitious mass 105 parameter µ. Comparing with TRBOMD, the parameter µ plays a similar role as ω −2 which controls the 106 frequency of the fictitious dynamics of the initial density guess in SCF iteration. This connection will be 107 made more explicit in the sequel. 108 We remark that the papers [15,16] made a further step in viewing TRBOMD by an extended 109 Lagrangian approach in a vanishing mass limit. However, unless very specific and restrictive form of the 110 error due to non-convergent SCF iterations is assumed, the equation of motion in TRBOMD does not 111 have an associated Lagrangian in general. The connection remains formal, and hence we will not further 112 explore here.

Analysis of TRBOMD in the linear response regime 114
In this section we consider Eq. (25) in the linear response regime, in which each atom I oscillates around its equilibrium position R * I . The displacement of the atomic configuration R from the equilibrium position is denoted by R(t) := R(t) − R * , and the deviation of the electron density from the converged density is denoted by ρ(x, t) := ρ(x, t) − ρ * (x; R(t)). Both R(t) and ρ(x, t) are small quantities in the linear response regime, and contain the same information as R(t) and ρ(x, t). Using R(t) and ρ(x, t) as the new variables and noting the chain rule due to the R-dependence in ρ * (x; R(t)), the equation of motion in TRBOMD becomes To simplify notation from now on we suppress the t-dependence in all variables, and Eq. (29) becomes In the linear response regime, we expand Eq. (30) and only keep terms that are linear with respect to R and ρ. All the higher order terms, including all the cross products of R I ,˙ R I , and ρ will be dropped. First we linearize the force on atom I with respect to ρ as Next we linearize with respect to R, we have Here the matrix {D IJ } is the dynamical matrix for the atoms. For the last term in Eq. (31) we have The last equation in Eq. (33) defines a linear functional L I , with δρ SCF δρ (x, y; R * ) and ∂V ion (x;R * ) at the fixed equilibrium point R * .

116
In the linear response regime, the operator δρ SCF δρ (x, y; R * ) carries all the information of the SCF iteration scheme. Let us now derive the explicit form of δρ SCF δρ (x, y; R * ) for the k-step simple mixing scheme with mixing parameter (step length) α (0 < α ≤ 1). If k = 1, the simple mixing scheme reads Here δ(x) is the Dirac δ-function, and the operator δ(x − y) − δF δρ (x, y) := ε(x, y) is usually refereed to as the dielectric operator [26,27]. To simplify the notation we would not distinguish the kernel of an integral operator from the integral operator itself. For example ε(x, y) is denoted by ε. Neither will we distinguish integral operators defined on continuous space from the corresponding finite dimensional matrices obtained from certain numerical discretization. This slight abuse of notation allows us to simply denote f (x) = A(x, y)g(y) dy by f = Ag as a matrix-vector multiplication, and to denote the composition of kernels of integral operators C(x, y) = dzA(x, z)B(z, y) by C = AB as a matrix-matrix multiplication. Using such notations, Eq. (35) can be written in a more compact form Similarly for the k-step simple mixing method, we have In general the dielectric operator is diagonalizable and all eigenvalues of ε are real. Therefore the 117 linear response operator δρ SCF δρ for the k-th step simple mixing method is also diagonalizable with real 118 eigenvalues.
From Eq. (30b) we have Here we have used the consistency condition (24). The last line of Eq. (38) defines a kernel which is an important quantity for the stability of TRBOMD as will be seen later. Using Eqs. (33) and (38), the equation of motion (30) can be written in the linear response regime as Define then Eq. (40) can be rewritten in a more compact form as Now if the self-consistent iteration is performed accurately regardless of the initial guess, i.e.
The linearized equation of motion (42) becomes Therefore in the case of accurate SCF iteration, according to Eq. (45a), the equation of motion of atoms follows the accurate linearized equation, and is decoupled from the fictitious dynamics of ρ. The normal modes of the equation of motion of atoms can be obtained by diagonalizing the dynamical matrix D as The frequencies {Ω l } (Ω l > 0) are known as phonon frequencies. When the SCF iterations are performed inaccurately, it is meaningless to assess the accuracy of the approximate dynamics (42) by direct investigation of the trajectories R(t), since small difference in the phonon frequency can cause large error in the phase of the periodic motion R(t) over long time. However, it is possible to compute the approximate phonon frequencies { Ω l } from Eq. (42), and measure the accuracy of TRBOMD in the linearized regime from the relative error The operator K(x, y) in Eq. (39) is directly related to the stability of the dynamics. Eq. (42b) also 120 suggests that in the linear response regime, the spectrum of K(x, y) must be on the real line, which 121 requires that the matrix δρ SCF δρ (x, y; R * ) be diagonalizable with real eigenvalues. This has been shown for 122 the simple mixing scheme. However, we remark that the condition that all eigenvalues of K(x, y) are real 123 may not hold for general preconditioners or for more complicated SCF iterations (for instance, Anderson 124 mixing). This is one important restriction of the linear response analysis. Of course, this may not be a 125 restriction for practical TRBOMD simulation for real systems. We will leave further understanding of 126 this to future works.

127
Let us now assume that all eigenvalues of K are real. The lower bound of the spectrum of K, denoted by λ min (K), should satisfy λ min (K) > 0.
Eq. (48) is a necessary condition for TRBOMD to be stable, which will be referred to as the stability condition in the following. Furthermore, ω should be chosen large enough in order to avoid resonance between the motion of R and ρ. Therefore the adiabatic condition should also be satisfied. Due to Eq. (49), we may assume ǫ = 1/ω 2 is a small number, and expand Ω l in the perturbation series of ǫ to quantify the error in the linear response regime. Following the derivation in the appendix, we have where K −1 is the inverse operator of K (K is invertible due to the stability condition). Since ω = √ κ/∆t, mainly determined by L, i.e. the accuracy of the SCF iteration.
Let us compare TRBOMD with CPMD. It is well known that CPMD accurately approximates the results of BOMD, provided that the electronic and ionic degrees of freedom remain adiabatically separated as well as the electrons stay close to the Born-Oppenheimer surface [11,12]. More specifically, the fictitious electron mass should be chosen so that the lowest electronic frequency is well above ionic frequencies where E gap is the spectral gap (between highest occupied and lowest unoccupied states) of the system, and recall that Ω l is the vibration frequency of the lattice phonon. For CPMD, a similar analysis in the linear response regime as above (which we omit the derivation here) shows that under the assumption (51).

131
Note that the condition (51) implies that CPMD no longer works if the system has a small gap or 132 is even metallic. The usual work-around for this is to add a heat bath for the electronic degrees of 133 freedom in CPMD [28], so that it maintains a fictitious temperature for the electronic degree of freedom.

134
Nonetheless the adiabaticity is lost for metallic systems and CPMD is no longer accurate over long time 135 simulation. In contrast, as we have discussed previously, TRBOMD may work for both insulating and 136 metallic systems without any modification, provided that the SCF iteration is accurate and no resonance 137 occurs. This is an important advantage of TRBOMD, which we will illustrate using numerical examples 138 in the next section.

139
When the system has a gap we can take µ sufficiently small to satisfy the adiabatic separation The total energy functional in our 1D density functional theory (DFT) model is given by The associated Hamiltonian is given by where Z I is an integer representing the charge of the i-th nucleus. This can be understood as a local pseudopotential approximation to represent the electron-ion interaction. The second term on the right hand side of Eq. (53) represents the electron-ion, electron-electron and ion-ion interaction energy. The parameter σ I represents the width of the nuclei in the pseudopotential theory. Clearly as σ I → 0, m I (x) → −Z I δ(x) which is the charge density for an ideal nucleus. In our numerical simulation, we set σ I to a finite value. The corresponding m I (x) is called a pseudo charge density for the I-th nucleus. We refer to the function m(x) as the total pseudo-charge density of the nuclei. The system satisfies charge neutrality condition, i.e.
Since m I (x) dx = −Z I , the charge neutrality condition (56) implies where N is the total number of electrons in the system. To simplify discussion, we omit the spin degeneracy here. The Hellmann-Feynman force is given by Instead of using a bare Coulomb interaction, which diverges in 1D, we adopt a Yukawa kernel which satisfies the equation As κ → 0, the Yukawa kernel approaches the bare Coulomb interaction given by the Poisson equation.

156
The parameter ǫ 0 is used to make the magnitude of the electron static contribution comparable to that of   In the linear response regime, we measure the error of the phonon frequency calculated from TRBOMD. This can be done in two ways. The first is given by Eq. (50), namely, all quantities in the big parentheses in Eq. (50) can be directly obtained by using the finite difference method at the equilibrium position R * . The second is to explore the fact that in the linear response regime, there is linear relation between the force and the atomic position as in Eq. (32), i.e. Hooke's law holds approximately at each time step. Here {f I (t l )} and { R I (t l )} are obtained from the trajectory of the TRBOMD simulation directly. To numerically compute D IJ , we solve the least square problem where The frequencies { Ω l } can be obtained by diagonalizing the matrix D. Similarly one can perform the 175 calculation for the accurate BOMD simulation and obtain the exact value of the frequencies {Ω l }.

176
In order to compare the performance among BOMD, TRBOMD and CPMD, we define the following relative errors err Hooke where the results from BOMD with convergent SCF iteration are taken to be corresponding reference   with converged SCF iteration, and BOMD(n) (resp. TRBOMD(n)) represents the BOMD (resp. TRBOMD) simulation with n SCF iterations per time step. It shows clearly that BOMD(5) produces large drift for both insulator (a) and metal (b), but TRBOMD(5) does not.  Table 1. However, for BOMD(5), the atom will cease oscillation after a while. A similar 193 phenomena occurs for other atoms. In Table 1, we present more results for TRBOMD(n) with n = 194 3, 5, 7. We observe there that TRBOMD(n) gives more accurate results with larger n, and err Hooke , err E , err L 2 R for TRBOMD(3) as a function of 1/ω 2 in logarithmic scales. When Fig. 4 shows clearly that all of |err Hooke Ω |, |err E |, |err L 2 R | depend linearly on 206 1/ω 2 . The error err L ∞ R has a similar behavior to err L 2 R and is skipped here for saving space.   Fig. 2. It shows clearly that the trajectory from TRBOMD(5) almost coincides with that from BOMD(c). However, for BOMD(5), the atom will cease oscillation after a while.

217
We now present some numerical examples for CPMD illustrating the difference between CPMD and 218 TRBOMD. As we have discussed, TRBOMD is applicable to both metallic and insulting systems, while 219 CPMD becomes inaccurate when the gap vanishes. To make this statement more concrete, we apply 220 CPMD to the same atom chain system. We implement CPMD using standard velocity Verlet scheme 221 combined with RATTLE for the orthonormality constraints [30][31][32]. 222 We present in Fig. 6 the error of CPMD simulation for different choices of fictitious electron mass µ. 223 We study the relative error of the phonon frequency err Hooke ρ SCF , and hence TRBOMD also works for metallic system as Fig. 4 indicates.    The different behavior of CPMD for insulating and metallic systems is further illustrated by Fig. 7 232 which shows the trajectory of the position of the left-most atom during the simulation. The phase error is 233 apparent from the two subfigures. While the phase error decreases so that the trajectory approaches that 234 of BOMD for insulator in Fig. 7(a), the result in Fig. 7(b) shows a systematic error for metallic system.  6. Beyond the linear response regime: Non-equilibrium dynamics 236 The discussion so far has been limited to the linear response regime so that we can make linear 237 approximations for the degrees of freedom of both nuclei and electrons. In this case, as the system 238 becomes linear, explicit error analysis has been given. For practical applications, we will be also 239 interested in non-equilibrium nuclei dynamics so that the deviation of atom positions is no longer small.

240
In this section, we will investigate the non-equilibrium case using averaging principle (see e.g. [33,34] 241 for general introduction on averaging principle).

242
Let us first show numerically a non-equilibrium situation for the atom chain example discussed before.

243
Initially, the 32 atoms stay at their equilibrium position. We set the initial velocity so that the left-most is "soft" without hard-core repulsion. In Fig. 9, we plot the difference between ρ SCF and the converged electron density of the SCF iteration (denoted by ρ KS ) along the TRBOMD simulation. We see that the 252 electron density used in TRBOMD stays close to the ground state electron density corresponds to the 253 atom configuration.
254 Figure 8. Comparison of trajectories of the first three atoms from the left for a non-equilibrium system. Different atoms are distinguished by color (blue for the initially left-most atom; green for the initially second left-most atom; red for the initially third left-most atom). Solid lines are results from BOMD(c); circled lines are results from TRBOMD(7); dashed lines are results from BOMD(7). It is evident that while results from BOMD with a non-convergent SCF iteration have a huge deviation, the results from TRBOMD are hardly distinguishable from the "true" results from BOMD. To understand the performance of TRBOMD, recall that the equations of motion are given by To satisfy the adiabatic condition (49) from the linear analysis, ω here is a large parameter. As a result,

257
Let us consider the limit ω → ∞. In this case, we may freeze the R degree of freedom in the equation of motion for ρ, as ρ changes on a much faster time scale. To capture the two time scale behavior, we introduce a heuristic two-scale asymptotic expansion with faster time variable given by τ = ωt (with some abuse of notations): and henceρ (x, t) = ω 2 ∂ 2 τ ρ(x, t, τ ) + 2ω∂ τ ∂ t ρ(x, t, τ ) + ∂ 2 t ρ(x, t, τ ).
(71) Figure 9. The difference of ρ SCF with the converged electron density of SCF iteration (denoted by ρ KS ) measured in L 1 norm along the TRBOMD simulation for a non-equilibrium system. Therefore, to the leading order, after neglecting terms of O(ω −1 ), we obtain For the equation of motion for ρ, note that as R only depends on t, the nuclear positions are fixed 258 parameters in Eq. (73).

259
To proceed, we consider the scenario that ρ(t, τ ) is close to the ground state electron density corresponding to the current atom configuration ρ * (R(t)). We have seen from numerical examples ( Fig. 9) that this is indeed the case for a good choice of SCF iteration, while we do not have a proof of this in the general case. Hence, we linearize the map ρ SCF .
and Eq. (73) becomes where K(R) is the same as in (39) except it is now defined for each atom configuration R. Let us 260 emphasize that here we have only taken the linear approximation for the electronic degrees of freedom, 261 while keeping the possibly nonlinear dynamics of R. This is different from the linear response regime 262 considered before, where the nuclei motion is also linearized.
Under the stability condition (48), it is easy to see that for ρ(t, τ ) satisfying Eq. (75), the limit of time average (76) Take the average of Eq. (72) in τ , we have Because of Eq. (76), the above dynamics is given by Remark. If we do not make the linear approximation for the electronic degree of freedom, as the map ρ SCF is quite nonlinear and complicated, the analysis of the long time (in τ ) behavior of Eq. (73) is not as straightforward. In particular, it is not clear to us whether the limit exists or how close the limit is to ρ * (x; R(t)) in a fully nonlinear regime. One particular difficulty lies in the fact that unlike BOMD or CPMD, we do not have a conserved Lagrangian for the TRBOMD. Actually, it is easy to construct much simplified analog of Eq. (73) that the average is different from ρ * . For example, if we consider the following analog which only has one degree of freedom ξ where (ξ/2 + aξ 2 ) is the analog of ρ SCF here and a > 0 is a small parameter which characterizes the nonlinearity of the map. Note thaẗ The motion of ξ is equivalent to a motion of a particle in an anharmonic potential. It is clear that if 267 initially ξ(0) = 0, the long time average of ξ will not be 0. Furthermore, if initially, ξ(0) is too large, the 268 orbit is not closed (ξ escapes the well around ξ = 0). If phenomena similar to this occur for a general 269 ρ SCF , then even in the limit ω → ∞, there will be a systematic uncontrolled bias between BOMD and TRBOMD. This is in contrast with Car-Parrinello molecular dynamics, which agrees with BOMD in the 271 limit fictitious mass goes to zero (µ → 0) if the adiabatic condition holds.

272
As a result of this discussion, in practice, when we apply TRBOMD to a particular system, we need to 273 be cautious whether the electronic degree of freedom remains around the converged Kohn-Sham electron 274 density, which is not necessarily guaranteed (in contrast to CPMD for systems with gaps). stability condition in TRBOMD is directly associated with the quality of the SCF iteration procedure.

285
In particular, we demonstrate in the case when the SCF iteration procedure is not very accurate, the 286 stability condition can be violated and TRBOMD becomes unstable. We also compare TRBOMD with 287 the Car-Parrinello molecular dynamics (CPMD) scheme. CPMD relies on the adiabatic evolution of 288 the occupied electron states and therefore CPMD works better for insulators than for metals. However,

289
TRBOMD may be effective for both insulating and metallic systems. The present study is restricted to 290 NVE system and to simplified DFT models. The performance of TRBOMD for NVT system and for 291 realistic DFT systems will be our future work.  Here we derive the perturbation analysis result in Eq. (50). When deriving the perturbation analysis below, we use linear algebra notation and do not distinguish matrices from operators. We use the linear algebra notation, replace all the integrals by matrix-vector multiplication, and drop all the dependencies of the electron degrees of freedom x and y. For instance, K ρ should be understood as K(x, y) ρ(y) dy. We also denote ∂ρ * ∂R (x; R * ) simply by ∂ρ * ∂R , then Eq. (42) can be rewritten as Here is a block diagonal matrix, and is a rank-M matrix. I is a M × M identity matrix. Now assume the eigenvalues and eigenvectors of A follows the expansion λ = λ 0 + ǫλ 1 + · · · , v = v 0 + ǫv 1 + · · · .
Match the equation up to O(ǫ), and Eq. (86a) implies that v 0 ∈ KerA 1 . Apply the projection operator P KerA 1 to both sides of Eq. (86b), and use that v 0 = P KerA 1 v 0 , we have Finally we apply v 0 to both sides of Eq. (86c) we have Therefore In other words, the phonon frequency Ω l = √ −λ up to the leading order is which is Eq. (50).