{'fraction_non_alphanumeric': 0.03485140872250096, 'fraction_numerical': 0.009340023157082208, 'mean_word_length': 4.53772173541355, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the existence of the optimal order for wavefunction extrapolation in Born-Oppenheimer molecular dynamics

Wavefunction extrapolation greatly reduces the number of self-consistent field (SCF) iterations and thus the overall computational cost of Born-Oppenheimer molecular dynamics (BOMD) that is based on the Kohn-Sham density functional theory. Going against the intuition that the higher order of extrapolation possesses a better accuracy, we demonstrate, from both theoretical and numerical perspectives, that the extrapolation accuracy firstly increases and then decreases with respect to the order, and an optimal extrapolation order in terms of minimal number of SCF iterations always exists. We also prove that the optimal order tends to be larger when using larger MD time steps or more strict SCF convergence criteria. By example BOMD simulations of a solid copper system, we show that the optimal extrapolation order covers a broad range when varying the MD time step or the SCF convergence criterion. Therefore, we suggest the necessity for BOMD simulation packages to open the user interface and to provide more choices on the extrapolation order. Another factor that may influence the extrapolation accuracy is the alignment scheme that eliminates the discontinuity in the wavefunctions with respect to the atomic or cell variables. We prove the equivalence between the two existing schemes, thus the implementation of either of them does not lead to essential difference in the extrapolation accuracy.


I. INTRODUCTION
Among ab-initio molecular dynamics (AIMD) simulation methods, Born-Oppenheimer molecular dynamics (BOMD) is considered to have good reliability and predictive power, especially for material systems with a small or vanishing band gap [1][2][3], and thus has been accepted as the default choice in several well-known AIMD simulation packages [4][5][6][7][8]. BOMD involves solving the quantum equation (for us, the Kohn-Sham equation [9][10][11]) for each instantaneous atomic configuration to obtain the atomic and cell forces, and propagating the atomic and cell variables according to classical Newtonian dynamics [1].
The computational cost of this very time-consuming procedure can be strongly reduced by the extrapolation technique [1,12], which combines results from previous MD time steps to predict a good initial guess for the self-consistent field (SCF) solution of the Kohn-Sham equation at the current MD time step. In the plane-wave method, the Kohn-Sham Hamiltonian matrix is never assembled explicitly and generally diagonalized by iterative schemes [5,11,13], so both the extrapolation of electron density and that of the Kohn-Sham eigenvalues and wavefunctions/orbitals are necessary. In this work, we discuss only the wavefunction extrapolation because it covers all difficulties in the extrapolation technique.
The existing extrapolation schemes could be classified into two categories [14]. The first is the time extrapolation (TX), which directly extrapolates the wavefunctions as functions of time. It can be implemented either by the Lagrange's polynomial extrapolation formula [15][16][17], or by the polynomial least-squares fitting proposed by literature on the Fock matrix dynamics [18,19]. In the following text, we refer to the first implementation when we use the terminology "TX". The second category is the geometric extrapolation (GX), which determines the wavefunction extrapolation coefficients by minimizing the distance between the current atomic coordinates and a linear combination of atomic coordinates from the previous time steps [20][21][22]. These schemes have been implemented in several popular AIMD simulation packages using the Kohn-Sham density functional theory and a plane-wave basis set, see Tab. I for an incomplete summary. It is worth noting that the order of the TX scheme usually refers to the order of the polynomial used to fit the wavefunctions, but the meaning of "order" in the context of GX scheme is not very clear. However, both TX and GX schemes can be uniformly characterized by the number of wavefunctions from previous MD time steps used in the extrapolation. Therefore, we accept this convention and denote The effectiveness of the extrapolation is calibrated by the number of iterations used to achieve a converged SCF solution under some user-defined convergence criterion. At first glance, the scheme using more previous-step information is more effective in predicting the wavefunctions. However, in practice, the number of SCF iterations was observed to firstly decrease and then to increase when extrapolating from an increasing number of previous wavefunctions [18,19]. That is to say, an optimal value of M that minimizes the number of SCF iterations was observed in these numerical tests. Moreover, it was reported, the optimal value of M gets larger when tightening the SCF convergence criterion [19]. These phenomena, as far as we see, lack proper theoretical explanation [23]. Most of the wellestablished AIMD packages provide very limited choices on the value of M (see Tab. I), and it is still questionable whether these configurations are optimal in terms of the number of SCF iterations.
In this work, we present an error analysis for the TX scheme in terms of the MD time step size and the SCF convergence criterion. Based on the error analysis, we claim that the optimal value of M always exists, and explain why the optimal M increases when using larger time step sizes or more strict SCF convergence criteria. Moreover, we theoretically elucidate the similarity between TX and GX schemes, so the theoretical findings on the TX scheme can be safely extended to the GX scheme.
As numerical examples, we select a solid copper system with 64 atoms, and carry out microcanonical simulations under different time step sizes and SCF convergence criteria. The appearance of optimal values of M well validate the error analysis, and a comparison between TX and GX confirms their similarity.
We demonstrate that the optimal value of M takes 2 in some cases, and takes as large as 5 in other cases. Therefore, we may suggest the necessity for AIMD simulation packages to open the user interface and provide more choices on the parameter M.
It should be noted that the extrapolations are linear combinations of wavefunctions from several adjacent time steps. An important underlying assumption, which is usually not satisfied in practice, is the continuity of these wavefunctions with respect to the atomic or cell variables [20,[24][25][26]. Therefore, a prerequisite step for the extrapolation is removing this discontinuity by using the wavefunction alignment scheme proposed either by Arias, Payne and Joannopoulos (APJ) [20], or by Mead [24]. Table I summarizes their implementations in the AIMD packages. We show that the corresponding extrapolated wavefunctions by using these two alignment schemes differ only up to a unitary transformation, so a particular choice between them does not lead to substantial difference in the extrapolation accuracy.
The remained part of this paper is organized as follows: The TX and GX schemes with an arbitrary value of M are given in Sec. II. Then Sec. III presents an error analysis for the TX scheme. As a prerequisite for the extrapolation, the wavefunction alignment is discussed in Sec. IV. In Sec. V, we validate our theoretical analysis by a set of microcanonical simulation tests on a solid copper system. Finally, the conclusion is made in Sec. VI.

II. WAVEFUNCTION EXTRAPOLATION SCHEMES
We consider the following extrapolation scheme, which is a linear combination of wavefunctions from M previous MD time steps, i.e. t n−M +1 , · · · , t n−1 , t n , and provides an initial guess for the SCF iteration at the (n + 1)-th time step: where X represents any of the N band wavefunctions, the superscript "p" indicates the predicted function, {c k } are extrapolation coefficients, and atomic positions R (t) ≡ (R 1 (t), . . . , R natom (t)) with row vector R I (t) denoting coordinate of the I-th atom at time t.

A. The time extrapolation scheme
The TX scheme fits the wavefunctions from M previous time steps, by an (M − 1)-th order Lagrange polynomial [27], and evaluates the polynomial at the (n + 1)-th time step to provide the initial guess for the SCF iteration. If a fixed time step is used, then the extrapolation coefficients {c k } are computed by The superscript "T" means that the coefficients are derived from TX scheme. Table II lists extrapolation coefficients of schemes with M from 2 to 5.
It is noted that the extrapolation coefficients can be easily computed for simulations with varying time steps.

B. The geometric extrapolation scheme
The extrapolation coefficients {c k } of the GX scheme are determined by minimizing the distance between atomic coordinates of the (n + 1)-th time step and the extrapolated atomic coordinates from M previous time steps: where the extrapolated atomic coordinates R ′ are given by The solution of the minimization problem (3) with respect to the extrapolation coefficients {c k } is given by the solution of the following linear system: where A is a (M − 1) × (M − 1) matrix, the elements of which are with ∆R j ≡ R(t j ) − R(t j−1 ). The solution of the optimization problem (3) is denoted by c G k .
Remark: When R(t) is exactly an (M −1)-th order polynomial function with respect to time variable t, then the extrapolation coefficients of the GX scheme are identical to those of the TX scheme, i.e. c G k = c T k , k = 1, · · · , M −1, because it is obvious that g(c T k ) = 0. In practice, the R(t) is not likely to be a polynomial, however, it can be very close to a polynomial in the time interval [t n−M +1 , t n+1 ] for sufficiently small time steps. For practically used time steps in AIMD simulations, we observe that the extrapolation coefficients of the TX and GX schemes are very similar. Therefore, the theoretical results derived for the TX scheme can be safely used to analyze the numerical phenomena of the GX scheme. When the R(t) is close to a polynomial with an order strictly lower than M − 1, the linear system (4) will be underdetermined, and has infinitely many solutions. What is observed in simulations is that the determinant of matrix A becomes even smaller than the machine precision. In this case we do not solve the system (4), but simply let c G k to be the same as the TX coefficients.

III. ERROR ANALYSIS OF THE WAVEFUNCTION EXTRAPOLATION
In practice the SCF procedure is always stopped at a finite number of iterations by some precision criterion (denoted by ǫ scf ), so the not-fully-converged SCF procedure inevitably introduces error in the wavefunctions. We denote this error by E(R (t)), and the wavefunction used in the simulations is the composition of two parts: whereX(R (t)) denotes the fully converged wavefunctions derived from the SCF procedure stopped at an arbitrarily high precision. Now we compute the difference between the predicted and fully converged wavefunctions at step t n+1 by in which the error of the wavefunction prediction is also composed by two parts, which are The first term on the R.H.S. of Eq. (8) is the error introduced by approximating the fully converged wavefunction at step t n+1 by extrapolating the fully converged wavefunctions at steps {t n , · · · , t n−M +1 }. We assume that the fully converged wavefunctions are smooth with respect to t, in the sense that they are infinitely differentiable in the time domain. For the Lagrange extrapolation, this error is expressed by tive of the wavefunction with respect to t and evaluated at t = ξ. Now if it is further assumed that the derivatives are uniformly bounded by some constant |d MX (R (t))/d t M | ≤ C, then the extrapolation error is upper bounded by The second term E scf on the R.H.S. of Eq. (8) is the extrapolation of the SCF error. It should be noted that the SCF error is, in general, not a continuous function with respect to time. Therefore an estimate of E scf like Eq. (11) would NOT hold for any M ≥ 1. This means that the Lagrange extrapolation of the SCF error would not be able to predict its value at future times. Actually, if the SCF error at different time steps are assumed to be independent random variables with identical distribution, we claim that (details of the proof are provided in Appendix B) the root mean-square magnitude of E scf grows exponentially with respect to the value of M, and has an estimate of where σ 2 scf = Var[E(R(t))] is the variance of the SCF error. In practice, the variance σ 2 scf decays monotonically with more strict SCF convergence criterion ǫ scf .
The error estimates Eqs. (11) and (13)  1. There always exists an optimal choice of M that minimizes the difference between the predicted and converged wavefunctions, and thus minimizes the number of SCF iterations. This is because as M increases, the magnitude of E xpl decreases exponentially, while the magnitude of E scf increases exponentially. Thus the optimal M is achieved at the cross-over of these two errors. In general, it is difficult to predict, a priori, the optimal value of M. The knowledge, at least for an equilibrated system, may be numerically obtained from short testing simulations, because the number of SCF iterations is also in equilibrium and does not significantly change with respect to time.
2. The optimal value of M appears earlier when using smaller time steps. The convergence of E xpl becomes fast when using smaller time steps, and the magnitude of E scf does not depend on the time step size, so the cross-over of the two errors occurs earlier.
That is to say a smaller M is preferred in simulations using smaller time steps.

3.
A larger M is preferred in simulations that use more strict SCF convergence criteria (i.e. smaller ǫ scf ). It is noticed that the magnitude of E scf is proportional to the standard deviation of the SCF error, and E xpl is independent with the SCF convergence criterion.
Therefore, smaller ǫ scf postpones the cross-over of the two errors, and the optimal M appears at a larger value.
The correctness of these qualitative guides will be checked by the numerical examples presented in Sec. V.
The stability of BOMD in presence of the incomplete SCF iterations can be analyzed by calculating the roots of the characteristic equations of the extrapolation schemes [15]. In this type of analysis, the SCF error was linearized and characterized by the largest eigenvalue of the SCF response kernel [28]. It was concluded that by using higher orders of extrapolation schemes, the region of the stable SCF response kernel decreases, which indicates that a more strict SCF convergence criterion should be used to guarantee the stability. What is predicted from our error estimate is consistent with the finding obtained from the stability analysis.

IV. WAVEFUNCTION ALIGNMENT
In this section, we denote two sets of wavefunctions from two successive MD time steps N band }, respectively, with overlap matrix ele- and suppose that S (0) is non-singular. The wavefunctions can not be combined directly in general, since the change from |Φ As discussed in [20,[24][25][26], it may be caused by various reasons including phase indeterminacies, band crossing, energy level splitting, etc. From a mathematical perspective, we understand the essence of these reasons as the change in the choice of orthonormal basis in the wavefunction subspace, which is defined as the linear subspace spanned by N band wavefunctions. In order to deal with this discontinuity problem, the orthonormal bases needs to be adjusted by unitary transformations on the original ones: where U Φ and U Ψ are both unitary matrices. It is easy to derive that the overlap matrix after the transformations becomes S = U † Ψ S (0) U Φ . The APJ scheme [20] minimizes the distance between the transformed bases, and computes the unitary transformations by where the superscript "A" indicates APJ scheme, and D is a diagonal matrix. Detailed derivation of the transformations is provided in Appendix A.
Mead, in a review paper on Berry phase [24], defined two sets of kets as in phase or parallel if their overlap matrix was both Hermitian and positive-definite, and any two sets of kets with a non-singular overlap matrix can be made parallel by a unitary transformation on a certain one of them. If letting U M Ψ = I, then is uniquely determined [24], where the superscript "M" denotes Mead's scheme . As discussed in Appendix A, (S (0) † S (0) ) −1/2 can be implemented through a singular value decomposition of S (0) † or a matrix diagonalization of S (0) † S (0) , and the latter is computationally more effective.
In the APJ scheme, the overlap matrix after transformation U A Ψ † S (0) U A Φ is a diagonal matrix with positive elements (see Appendix A). Consider that we keep { |Ψ is obviously both Hermitian and positive-definite. According to the uniqueness of U M Φ , we conclude that By using (16), it is straightforward to show that the extrapolated wavefunctions using these two alignment schemes differ only by a unitary transformation: Noticing that the extrapolated wavefunctions should be orthonormalized before using as the initial guess for the SCF iteration, the theoretical comparison of extrapolation accuracies by using the two alignment schemes is indeed a subtle problem. The numerical performance of these two alignment schemes are similar due to the presented equivalency (17), and will be discussed in Sec. V.

V. NUMERICAL RESULTS
We have implemented the TX and GX schemes mentioned in Sec. II and also the APJ and Mead's wavefunction alignment schemes in Sec. IV, in an in-house ab-initio simulation package CESSP developed on infrastructure JASMIN [29], and carried out systematic tests on a solid copper system with 64 atoms to validate our theoretical analysis. At each MD time step, the Kohn-Sham equation under the generalized gradient approximation (GGA) [30] was solved by the plane-wave pseudopotential method [5,11,13]. The kinetic energy cutoff for the plane-wave expansion was chosen as 350 eV, and a 2 × 2 × 2 k-points mesh was used. In the SCF iteration, the Pulay's mixing scheme [31] was employed to mix the electronic charge densities, and the SCF convergence was considered to be reached when differences in both the potential energy and band structure energy between successive SCF steps are smaller than the criterion ǫ scf . The iterative matrix diagonalization of the Kohn-Sham Hamiltonian was realized by the block Davidson scheme [32].
The initial coordinates and velocities were obtained by a warm-up microcanonical simulation, which started from a perfect FCC configuration with atomic velocities sampled from a Maxwell-Boltzmann distribution at temperature 1400 K. The system was equilibrated after  ǫ scf = 10 −3 eV ǫ scf = 10 −4 eV ǫ scf = 10 −5 eV ǫ scf = 10 −6 eV ǫ scf = 10 −7 eV ǫ scf = 10 −8 eV  As a reference, we also provide the energy drift results in Tabs. III and IV. They were evaluated by fitting the total energies to a linear function of time. This is a reasonable ansatz because the time evolution of the total energy drift presents a linear dependence on the simulation time, which is consistent with earlier results on the energy conservation of BOMD reported in e.g. Refs. [1,18,19]. It should be noted that neither the TX nor GX scheme is designed to conserve the total energy of the system, therefore, not surprisingly, a systematic  which cannot be entirely removed even when using a very strict SCF criterion (e.g. ǫ scf = 10 −8 eV), is considered to stem from a broken time reversal symmetry along with the implicit electronic propagation composed of the extrapolation step and the SCF procedure [12,34].
It has been remarked in Sec. II that extrapolation coefficients of the GX scheme are in fact close to those of the TX scheme, especially for small time step sizes. To test this point,   a steady state only after a few tens of time steps. As an example, the case of ∆t = 1.0 fs and ǫ scf = 10 −5 eV is shown in Fig. 1.

VI. CONCLUSIONS AND DISCUSSIONS
In this work, we investigated both theoretically and numerically the optimal choice of the critical parameter M in the wavefunction extrapolation, which denotes the number of previous MD time steps used to construct the initial guess for the SCF iteration. We focused on two types of extrapolations schemes, i.e. time extrapolation (TX) and geometric extrapolation (GX) that are most frequently implemented in the popular AIMD simulation packages using the Kohn-Sham density functional theory and a plane-wave basis set.
We showed that the extrapolation error is divided into two parts: The first stems from extrapolating the fully converged wavefunctions and decays exponentially with respect to the parameter M. The second is a linear combination of the SCF errors and explodes exponentially with respect to the parameter M. Therefore, an optimal value of M always exists. Since the convergence speed of the first term is governed by the time step and the magnitude of the second term is proportional to the SCF error, we qualitatively predicted that the optimal M tends to get larger when increasing the MD time step size or tightening the SCF convergence criterion. Moreover, we theoretically explained the reason why the TX and GX schemes have much similarity, therefore the analysis for the TX scheme can be extended to the GX scheme without substantial difficulty. A For the prerequisite alignment step, we theoretically proved that the extrapolated wavefunctions by using APJ and Mead's schemes differ only by a unitary transformation, and illustrated numerically the essential equivalence between them. The only concern is for the numerical implementation, that a matrix diagonalization of S (0) † S (0) other than a singular value decomposition of S (0) † is recommended due to the computational efficiency.
Although the TX and GX schemes could significantly accelerate SCF convergence, it should be noted that neither of them is designed for total energy conservation. Correspondingly, a systematic total energy drift in microcanonical trajectories exists, and increasing M or decreasing the time step does not guarantee a smaller drift. Therefore, in the case when the energy conservation is critical, we suggest either a very strict SCF convergence criterion, or new techniques like the ASPC or XL-BOMD method.
In the end, we would like to point out that our work only prove the existence of the optimal value M, and do not tell, a priori, how large M is for a specific system and how much the advantage is by using the optimal M instead of the widely accepted value of 3.
We notice that for the tested copper system, M = 3 is already optimal under the normally used simulation setups: A time step of 1 ∼ 2 fs and a SCF convergence criterion of 10 −4 ∼ 10 −5 eV. However, this experience cannot be directly extended to other systems. Therefore the value of this work is to point out the existence of the optimal M, and the possibility of improving the extrapolation quality, marginally or dramatically, by tuning the value of M. In the simulations of equilibrium systems, the optimal M can be estimated by short testing simulations of a few dozen time steps. Moreover, the direction of changing M is also predicted by our theory when using a different time step size or SCF convergence criterion.
In the simulations of non-equilibrium systems, the optimal choice of M is still an open question, and will be investigated in the future.

VII. ACKNOWLEDGMENT
The authors acknowledge the valuable discussions with Dr. Deye Lin, and gratefully ac- We need to solve the following constrained minimization problem: ν (U A Ψ ) νj are unitary transformed functions from two sets of wavefunctions given by successive MD time steps, w j > 0 the weight of the j-th function, and Ω the computational domain. It is convenient to rewrite objective function f in matrix form: where tr represents the trace of a matrix, W ≡ diag{w 1 , . . . , w N band }, I is the identity matrix, and S (0) and S are overlap matrices before and after the unitary transformations, respectively.
where Lagrange multiplier matrices Λ Φ and Λ Ψ are both Hermitian. Taking partial derivatives of L with respect to real and imaginary parts of matrix elements of U A Φ and U A Ψ , we get It is easy to see from (A2) that which means Λ Φ and Λ Ψ are both commutative with W 2 . Since this is true for any W = diag{w 1 , . . . , w N band } (w j > 0), we assume that any two matrix elements of W are not equal, and conclude that Λ Φ and Λ Ψ are diagonal matrices. Thus, the second equation in (A2) implies that i.e., U A Φ is determined as the unitary matrix that diagonalizes S (0) † S (0) . Further, according to (A2) we have Because |s j | = | Ψ j |Φ j | ≤ 1, f reaches its minimum when s j ≥ 0. Thus, we get

Φ
One choice is to carry out singular value decomposition of S (0) † : where U 1 and U 2 are both unitary matrices, and Σ is a diagonal matrix with non-negative elements. Since S (0) is assumed to be non-singular, matrix elements of Σ is in fact positive. Then It is easy to see that (U 1 Σ 2 U † 1 ) −1/2 = U 1 Σ −1 U † 1 , so we have The other choice is matrix diagonalization: where U is a unitary matrix and D is a diagonal matrix with positive elements. It is easy to see that Since singular value decomposition is more computationally consuming than diagonalization, we suggest the latter choice in numerical implementation.

Appendix B: Estimate the error E scf
In this section we prove Eq. (13). We assume that the SCF errors E(R(t n )) are independent random variables that have vanishing mean and identical variance denoted by σ 2 scf . As a short-hand notation, we write E n = E(R(t n )). Noticing Eq. (10), it is clear that the mean of E scf vanishes, and the variance of E scf is estimated by