Comparing BOMD and CPMD

The comparison of the overall performance of Car-Parrinello and Born-Oppenheimer molecular dynamics in terms of computer time is a delicate issue. For instance it depends crucially on the choice made concerning the accuracy of the conservation of the energy $E_{\rm cons}$ . Thus, this is to some extend subject of ``personal taste'' as to what is considered to be a ``sufficiently accurate'' energy conservation. In addition, this comparison might lead to different conclusions as a function of type and size of the system investigated, and how much the specific idiosyncrasies of each method affects the investigated properties. One major advantage of CP-dynamics is the fact that there is determistically only one ``wavefunction optimization'' step per time step and a significant disadvantage is the red-shift in the dynamics of light atoms (``drag'') due to the fictitious electron dynamics. BO-dynamics does not suffer from the latter, but for systems where convergence of the wavefunction is difficult the computational cost can be much higher. Recent developments in using wavefunction extrapolation [41,10] to improve the quality of the initial guess wavefunction in terms of forces and energy conservation have made BO-dynamics much more attractive. Nevertheless, in order to shed light on this point, microcanonical simulations of 8 silicon atoms were performed with various parameters using Car-Parrinello and Born-Oppenheimer molecular dynamics. This large-gap system was initially extremely well equilibrated and the runs were extended to 8 ps (and a few to 12 ps with no noticeable difference) at a temperature of about 360-370 K (with $\pm 80$  K root-mean-square fluctuations). The wavefunction was expanded up to $E_{\rm cut} = 10$  Ry at the $\Gamma$ -point of a simple cubic supercell and LDA was used to describe the interactions. In both cases the velocity Verlet scheme was used to integrate the equations of motion. In Car-Parrinello molecular dynamics two different time steps were used, 5 a.u. and 10 a.u. (corresponding to about 0.24 fs), in conjunction with a fictitious electron mass of $\mu = 400$  a.u.. Within Born-Oppenheimer molecular dynamics the minimization of the energy functional was done using the highly efficient DIIS (direct inversion in the iterative subspace) scheme [42]. In this case, the time step was either 10 a.u. or 100 a.u. and three convergence criteria were used; note that the large time step corresponding to 2.4 fs is already at the limit to be used to investigate typical molecular systems (with frequencies up to 4000 cm$^{-1}$ ). The convergence criterion is based on the largest element of the wavefunction gradient which was required to be smaller than 10$^{-6}$ , 10$^{-5}$ or 10$^{-4}$  a.u..

Figure 2: Conserved energy $E_{\rm cons}$ from Car-Parrinello (CP) and Born-Oppenheimer (BO) molecular dynamics simulations of a model system for various time steps and convergence criteria; see text for further details and Table 1 for the corresponding timings. Top: solid line: CP, 5 a.u.; open circles: CP, 10 a.u.; filled squares: BO, 10 a.u., 10$^{-6}$ . Middle: open circles: CP, 10 a.u.; filled squares: BO, 10 a.u., 10$^{-6}$ ; filled triangles: BO, 100 a.u., 10$^{-6}$ ; open diamonds: BO, 100 a.u., 10$^{-5}$ . Bottom: open circles: CP, 10 a.u.; open diamonds: BO, 100 a.u., 10$^{-5}$ ; dashed line: BO, 100 a.u., 10$^{-4}$ .

The outcome of this comparison is shown in Fig. 2 in terms of the time evolution of the conserved energy $E_{\rm cons}$ on scales that cover more than three orders of magnitude in absolute accuracy. Within the present comparison ultimate energy stability was obtained using Car-Parrinello molecular dynamics with the shortest time step of 5 a.u., which conserves the energy of the total system to about 6$\times$ 10$^{-8}$  a.u. per picosecond, see solid line in Fig. 2(top). Increasing the time step to 10 a.u. leads to an energy conservation of about 3$\times$ 10$^{-7}$  a.u./ps and much larger energy fluctuations, see open circles in Fig. 2(top). The computer time needed in order to generate one picosecond of Car-Parrinello trajectory increases linearly with the time step, see Table 1. The most stable Born-Oppenheimer run was performed with a time step of 10 a.u. and a convergence of 10$^{-6}$ . This leads to an energy conservation of about 1$\times$ 10$^{-6}$  a.u./ps, see filled squares in Fig. 2(top).

Table: Timings in CPU seconds and energy conservation in a.u. / ps for Car-Parrinello (CP) and Born-Oppenheimer (BO) molecular dynamics simulations of a model system for 1 ps of trajectory on an IBM RS6000 / model 390 (Power2) workstation using the CPMD package; see Fig. 2 for corresponding energy plots.

Method

Time step (a.u.)

Convergence (a.u.)

Conservation (a.u./ps)

Time (s) CP

As the maximum time step in Born-Oppenheimer dynamics is only related to the time scale associated to nuclear motion it could be increased from 10 to 100 a.u. while keeping the convergence at the same tight limit of 10$^{-6}$ . This worsens the energy conservation slightly (to about 6$\times$ 10$^{-6}$  a.u./ps), whereas the energy fluctuations increase dramatically, see filled triangles in Fig. 2(middle) and note the change of scale compared to Fig. 2(top). The overall gain is an acceleration of the Born-Oppenheimer simulation by a factor of about seven to eight, see Table 1. In the Born-Oppenheimer scheme, the computer time needed for a fixed amount of simulated physical time decreases only sub linearly with increasing time step since the initial guess for the iterative minimization degrades in quality as the time step is made larger. Further savings of computer time can be easily achieved by decreasing the quality of the wavefunction convergence from 10$^{-6}$ to 10$^{-5}$ and finally to 10$^{-4}$ , see Table 1. This is unfortunately tied to a significant decrease of the energy conservation from 6$\times$ 10$^{-6}$  a.u./ps at 10$^{-6}$ (filled triangles) to about 1$\times$ 10$^{-3}$  a.u./ps at 10$^{-4}$ (dashed line) using the same 100 a.u. time step, see Fig. 2(bottom) but note the change of scale compared to Fig. 2(middle).

In conclusion, Born-Oppenheimer molecular dynamics can be made as fast as (or even faster than) Car-Parrinello molecular dynamics (as measured by the amount of CPU time spent per picosecond) at the expense of sacrificing accuracy in terms of energy conservation.