Velocity Verlet Equations for CPMD

The appearance of the constraint terms complicates the velocity Verlet method slightly for CPMD. The Lagrange parameters $\Lambda_{ij}$ have to be calculated to be consistent with the discretization method employed. How to include constraints in the velocity Verlet algorithm has been explained by Andersen [25]. In the following we will assume that the overlap is not position dependent, as it is the case for plane wave basis sets. The more general case will be explained when needed in a later section. These basic equations and generalizations thereof can be found in a series of papers by Tuckerman et al. [33].

For the case of overlap matrices that are not position dependent the constraint term only appears in the equations for the orbitals.

$$\mu \mid {\ddot \Phi}_i(t) \rangle = \mid \varphi_i(t) \rangle + \sum_{j} \Lambda_{ij} \mid \Phi_j(t) \rangle$$ (50)

where the definition

$$\mid \varphi_i(t) \rangle = - f_i H^{\rm KS} \mid \Phi_i(t) \rangle$$ (51)

has been used. The velocity Verlet scheme for the wavefunctions has to incorporate the constraints by using the RATTLE algorithm. The explicit formulas will be derived in the framework of plane waves in a later section. The structure of the algorithm for the wavefunctions is given below

       CV(:) := CV(:) + dt/(2*m)*CF(:)
       CRP(:) := CR(:) + dt*CV(:)
       Calculte Lagrange multiplier X
       CR(:) := CRP(:) + X*CR(:)
       Calculate forces CF(:) = HKS*CR(:)
       CV(:) := CV(:) + dt/(2*m)*CF(:)
       Calculte Lagrange multiplier Y
       CV(:) := CV(:) + Y*CR(:)