An Efficient and Accurate Car-Parrinello-like Approach to Born-Oppenheimer Molecular Dynamics

We present a new method which combines Car-Parrinello and Born-Oppenheimer molecular dynamics in order to accelerate density functional theory based ab-initio simulations. Depending on the system a gain in efficiency of one to two orders of magnitude has been observed, which allows ab-initio molecular dynamics of much larger time and length scales than previously thought feasible. It will be demonstrated that the dynamics is correctly reproduced and that high accuracy can be maintained throughout for systems ranging from insulators to semiconductors and even to metals in condensed phases. This development considerably extends the scope of ab-initio simulations.

Density functional theory (DFT) [1] based ab-initio molecular dynamics [2], in which interactions are computed on the fly from electronic structure calculations, has been very successful in describing a large variety of physical phenomena and has proven its relevance in many fields.However, its computational cost has limited the attainable length and time scales in spite of substantial progress.For a while it was believed that the development of linear scaling methods [3,4,5] could have offered a solution.Unfortunately, the crossover point at which linear scaling methods become advantageous has remained fairly large, especially if high accuracy is needed.
On the other hand, it would be very desirable to accelerate ab-initio simulations with up to thousands of atoms, such that simulations as long as tens or even hundreds of picoseconds can be routinely performed, thus making completely new phenomena accessible to ab-initio simulations.Born-Oppenheimer molecular dynamics (BOMD) simulations, in which at each time step the DFT functional is fully minimized, do not seem to offer much room for further improvement.For this reason another direction has recently been followed to improve the efficiency at current system sizes.In the spirit of Car-Parrinello molecular dynamics (CPMD) [2], some form of dynamics for the electronic degrees of freedom is implemented, which automatically maintains the system close to the BO surface, but at variance with the original proposal in a localized orbital representation [6,7,8,9].The acceleration stems from the ability to reduce or fully bypass the self-consistency cycle.However, just like in CPMD, these methods suffer from rather short integration time steps.In this Letter we present a new method which overcomes this limitation and combines the accuracy and long time steps of BOMD with the efficiency of CPMD.
We consider the general case in which the DFT Kohn-Sham orbitals are expanded in a non-orthogonal basis.Let M be the dimension of the Hilbert space and S the M × M overlap matrix.As is customary in quantum chemistry we arrange the expansion coefficients of the N lowest occupied orbitals in a rectangular M × N matrix C. The density matrix is then written as P = CC T and obeys the idempotency condition P = PSP.Maintaining the idempotency of P is crucial in any electronic structure calculation and is one of its algorithmic challenges.The potential energy surface on which the ions move is defined by the minimum of an appropriately chosen energy functional E DFT [C, R I ], which we express as a functional of C and a function of the ionic coordinates R I .In this notation the Born-Oppenheimer equation of motions reads as follows: where the search for the minimum is restricted to the C's that satisfy the orthonormality condition C T SC = I, which is equivalent to imposing the idempotency condition on P. The forces of Eq. ( 1) can be divided into three contributions, the Hellmann-Feynman forces [10,11], the Pulay forces [12], which are present whenever the basis set depends on the ionic positions, and a residual term [13] which is non-zero except when full self-consistency is reached.This term leads to poor energy conservation in BOMD unless a very tight convergence criterion is imposed.In Car-Parrinello-like approaches this is circumvented by the design of a coupled electron-ion dynamics which maintains the system close to the BO surface, but at the cost of small time steps.Based on ideas of the original Car-Parrinello approach we design an improved dynamics for the coupled system of electrons and ions.Contrary to the Car-Parrinello scheme we will not write an explicit equation of motion for the C's, but rather an integration scheme for the electronic degrees of freedom.The knowledge of the previous K values of C (t n−l ), where l ∈ [1, K], determines the value of C (t n ), such that at any instant of time the C's are as close as possible to the instantaneous ground state.As for the short-term integration of the electronic degrees of freedom, accuracy is crucial, a highly accurate and efficient algorithm is required.Here we have chosen the always stable predictor-corrector (ASPC) method [14] of Kolafa.This method was originally designed to deal with classical polarization, so that additional care must be taken that during the evolution the idempotency condition is satisfied.The predictor uses the extrapolated contra-covariant density matrix PS [15] as an approximate projector on to the occupied subspace C (t n−1 ).In this way, we take advantage of the fact that the physically relevant contra-covariant density matrix PS evolves more smoothly and is therefore easier to predict than C.This is followed by a corrector step to minimize the error and to further reduce the deviation from the ground state.The corrector consists of a single minimization step MIN [C p (t n )] of a properly selected minimization procedure.Alternatively, the predictor can also be repeatedly applied.In this case the ground state is even more closely approached, but in general this is not necessary.The numerical coefficients of Eq. ( 2) were selected by Kolafa in order to ensure time-reversibility up to O(h K+2 ), while ω was chosen to guarantee a stable relaxation towards the minimum.Because the energy is invariant under a unitary transformation within the subspace of occupied orbitals C, it must be ensured that this gauge transformation is not strongly changed by MIN [C p (t n )], as in this case continuity between the C's may be lost.Moreover the minimization scheme must be very efficient in bringing the system close to the ground state and preserving the idempotency of the density matrix.For these reasons we have chosen the orbital transformation (OT) method of VandeVondele and Hutter [16].Inspired by the form of the exponential transformation [17] an auxiliary variable X is introduced, to parameterize the occupied orbitals where U = X T SX 1/2 and the variable X has to obey the linear constraint: Under this condition C (X) leads to an idempotent density matrix for any choice of X, provided that the reference orbitals C p (t n ) are orthonormal.Thus a finite minimization step along the gradient direction will still fulfill the idempotency constraint exactly.Due to the linear constraint the minimization with respect to X is performed in an auxiliary tangent space.Because this is a linear one no curved geodesics must be followed, as is the case for variables such as C, which are nonlinearly constrained.In this way large minimization steps can be taken, especially if a good preconditioner is used.In fact using a very efficient, idempotency conserving direct minimizer like OT has been decisive for the success of this approach.Since the ASPC integrator preserves the orthonormality constraint only approximately, it occasionally has to be explicitly enforced, either by a Cholesky decomposition or by a few purification iterations [18].
Having obtained the new wavefunction we can now evaluate the energy and the nuclear forces, which are derived from the following approximate energy functional: where ρ p (r) is the density associated with C p (t n ).E PC [ρ] can be thought of as an approximation to the Harris-Foulkes functional [19] and maintains the predictorcorrector flavor of the present work.The validity of E PC [ρ] depends only on the efficiency of the minimizer and on the quality of the propagation scheme.The ionic forces are calculated by evaluating the analytic gradient of E PC [ρ], with respect to the nuclear coordinates.However, as ∆ρ = ρ − ρ p = 0, besides the usual Hellmann-Feynman and Pulay forces an extra term appears: where ρ is the corrected density evaluated using C(t n ) and ρ p is the predicted density calculated from C p (t n ).However, as a single minimization step is performed, C(t n ) is only an approximate eigenfunction of H[ρ p ] within the subspace spanned by the finite basis set used [20].This leads to an insignificant error in the forces, provided that C(t n ) is very close to the ground-state.The ability of this dynamics to maintain the system on the BO surface may vary considerably.It is almost perfect in systems like water, but somewhat less satisfactory in liquid Si at high temperature, where swift bonding and rebonding processes take place.However, in all cases the dynamics is dissipative, most likely due to the fact that the propagation scheme employed is not symplectic.It is possible to remedy this downward drift if we assume that the forces arising from our dynamics F PC can be modeled as F PC = F BO − γ D ṘI which, as we shall see later, is an excellent assumption.The value of the intrinsic friction coefficient γ D does not need to be known but it can be bootstrapped by taking a cue from the work of Krajewski and Parrinello [21].We sample the canonical distribution by using the following Langevin-type equation where M I is the ionic mass, γ L is a Langevin friction coefficient and Ξ I = Ξ D I + Ξ L I a random noise.Using the assumption above this can be equally written as: In order to guarantee an accurate sampling of the Boltzmann distribution, the noise has to obey the fluctuation dissipation theorem: The choice of γ L is arbitrary, while the unknown γ D has to be determined by requiring that the aggregate noise term generate the correct average temperature, as measured by the equipartition theorem 1 2 M I Ṙ2 I = 3 2 k B T .We will show that this leads to correct sampling.Since our initial dynamics is quite accurate γ D is small, so that dynamical properties are also well reproduced.
For the purpose of demonstrating our new method, we have implemented it in the mixed Gaussian Plane Wave (GPW) [22] code Quickstep [23] which is part of the publicly available suite of programs CP2K [24].We present here calculations on metallic liquid silicon and liquid silica, to illustrate that our method works well irrespective of band-gap, system size and system type.Both systems are known to be very difficult, and are examples of liquid metals (Si) and of complex, highly polarizable, ionic liquids (SiO 2 ).In addition the simulations have been performed at 3000 K and 3500 K respectively, which leads to rapidly varying density matrix elements, thus making the propagation of the electronic degrees of freedom quite challenging.Hence, the selected tests can be considered as worst-case scenarios for our method.
All simulations have been performed at their experimental liquid densities using double-zeta valence polarization (DZVP) basis sets, adequate density cutoffs, the Goedecker-Teter-Hutter pseudopotentials [25] and the local-density approximation to the exact exchange and correlation functional.For simplicity the Brillouin zone is sampled at the Γ-point only.Eq. ( 8) is integrated using the algorithm of Ricci and Ciccotti [26], with a time step of h = 1.0 fs.The friction coefficient γ L was set equal to zero, while the values for γ D turned out to be in the range of 10 −4 fs −1 .We predict the new C's using K = 4 in Eq. ( 2), which ensures time-reversibility up to O(h 6 ).The OT minimizer is preconditioned, as proposed in [16].
We first consider the accuracy in terms of deviation from the BO surface.As can be seen in FIG. 1 our energies are an upper bound to the ground state and are displaced by a small and approximately constant amount.It is also shown that the deviation from the BO surface can be even further reduced by increasing the number of corrector steps.In fact the deviation can be controlled by varying the number of corrector steps in order to achieve a preassigned accuracy level.In the following only simulations based on a single corrector step will be reported.
We turn now to more realistic problems such as those shown in FIG. 2. Although these simulations have been performed with only a single corrector step, they are still amazingly close to the BOMD results.It should be emphasized that even in liquid Si, which poses problems when using an ordinary Car-Parrinello scheme due to its metallicity, a single corrector step is sufficient.This establishes the efficiency of our method, since only a single preconditioned gradient calculation and no selfconsistent iteration has to be performed.The possible acceleration, in comparison with regular BOMD calculations, depends crucially on the system studied.In these difficult cases a speed-up of two orders of magnitude compared to using the pure extrapolation scheme have been observed.For simpler problems an increase in efficiency of at least one order of magnitude can be expected.
In FIG. 3 we present results which prove that also dynamical properties can be evaluated with accuracy.This time we look at the velocity autocorrelation function and its Fourier transform at 325 K.The results are in good agreement with our reference calculations and are consistent with experiment and ab-initio all-electron calculations [27], showing that in spite of the stochastic nature of Eq. ( 8) dynamical properties can also be simulated.This implies, that also non-equilibrium processes and chemical reactions can be handled.In the same picture we verify that our assumptions are justified, and we are indeed performing a canonical sampling, by showing that the kinetic energy distribution is Maxwellian.To this effect, we have carried out a 64 atom liquid Si simulation for as long as 1 ns, to reduce the noise and to ensure a proper sampling of the kinetic energy distribution tails.
Because of space considerations we have reported here only a fraction of the systems studied.In all cases our method has proven to be accurate and the gain in speed has always been remarkable.Also structure relaxation via annealing and geometry optimization have been successfully performed.Contrary to CPMD and related methods at least BOMD integration time steps can be used.Thanks to this development it is now possible to perform ab-initio molecular dynamics on medium-sized systems up to a few nanoseconds, thus making a new class of problems accessible to ab-initio simulations.The computational scaling of the algorithm is in principle linear.However, an efficient parallel sparse matrix algebra has not been implemented yet, so that the sustained scaling is still O(M N 2 ), albeit with a much reduced prefactor.
In conclusion we wish to highlight that our ideas are of general interest, as our method is independent of the particular choice of the minimizer and propagation scheme, and can be further improved.In fact, during the submission process we became aware of a time-reversible extrapolation scheme [28], which could be used as an alternative to ASPC.The possibility to apply the presented ideas with benefit to plane wave based CPMD is also to be underlined and will be presented elsewhere.
We would like to thank the whole CP2K team, in particular J. Hutter and J. VandeVondele for sharing their deep insights, and F. R. Krajewski and J. Kolafa for fruitful discussions.The generous allocation of computer time and support from CSCS Manno, the ICT Services of ETH Zurich and the DEISA consortium is acknowledged.

FIG. 1 :FIG. 2 :
FIG.1: Deviations from the BO surface of liquid SiO2 with respect to total energies (upper panel) and mean force deviations (lower panel).The deviation in the energies corresponds to a constant shift of 4.16×10−4 Hartree per atom for one corrector step and 3.5 × 10 −5 Hartree per atom for two corrector steps.The average mean force deviation is unbiased.

FIG. 3 :
FIG.3:The kinetic energy distribution calculated from a 1 ns trajectory with a time step of 3.25 fs of metallic liquid Si64 using a DZVP Gaussian basis set and a density cutoff of 100 Ry (left panel).Velocity autocorrelation function (upper right) and its Fourier transform (lower right) of 32 Water at 325 K using a TZV2P Gaussian basis set, a density cutoff of 280 Ry and the BLYP exchange-correlation functional.The Langevin friction coefficients are γL = 0 and γD ∼ 10 −8 fs −1 .