A number of first-principles molecular dynamics simulations of biologically important systems, e.g. molecular liquids [1][2][3] and biomolecules in aqueous solution, 4) are emerging in the past few years. Since these calculations are computationally very demanding, clever use of CPU time is crucial. One way to accelerate the simulations is to use the largest possible timestep, which is mainly limited by the highest vibrational frequency of the system. In the systems mentioned above, the highest frequencies arise from stretching motion of the covalent bonds, which are only weakly coupled to other modes in most cases. Therefore, by freezing these fast modes, we can use much larger timesteps without changing the long-time behavior of the system. 5, 6) Bond-angle constraints, on the other hand, will not be considered in this article, since they can significantly alter the properties of the system, with little gain in efficiency.
5, 6)Although the constrained molecular dynamics (CMD) is a well-established technique for classical simulations, 5, 6) its applications within the first-principles calculations are limited to the estimation of free-energy differences associated with chemical reactions.7-9) In these simulations, constraints were imposed only on the reaction coordinates, and the time evolution of the system was followed by the fictitious dynamics.10) In contrast, CMD has rarely been used for the Born-Oppenheimer dynamics (BOD), in which the electronic ground state is reached at every ionic step. In general, the number of iterations for electronic minimization increases for larger timesteps in BOD, since the quality of the extrapolated wavefunctions deteriorates.11) Therefore, it is not obvious if CMD provides any improvement of the performance. In this article, we investigate the efficiency of CMD within BOD as well as the effects of using deuterium atoms in place of hydrogen atoms.We have implemented CMD within our first-principles molecular dynamics code, FEMTECK,[12][13][14] along the lines of RATTLE algorithm.15) Note that the computational efforts to calculate the constraint forces are negligible compared with those for the (unconstrained) ionic forces in the first-principles calculations. Therefore, we solve the nonlinear equations for the bond lengths (eqs. (A1) and (A3) gives constraint forces converged to machine precision in three to four iterations. Linear equations for the velocities (eqs. (A2) and (A4) of ref. 15) are also solved directly by matrix inversion.The performance of CMD has been examined by constant-NVE molecular dynamics simulations of a cytosine molecule in the amino-oxo form (Fig. 1) in a cubic supercell of (16 a.u.) 3 . The local density approximation 17) of the density functional theory 18) was used with an average cutoff energy of 39.5 Ry. The equations of motion for the ions were integrated with the velocity-Verlet algorithm.19) The wavefunctions were extrapolated from previous ionic steps 11) during the simulations. Bond lengths were fixed at their equilibrium values when con...