Car-Parrinello Equations

The Car-Parrinello Lagrangian and its derived equations of motions were introduced before. Here the equations are specialized to the case of a plane wave basis within Kohn-Sham density functional theory using norm-conserving pseudopotentials. Specifically, the functions $\Phi_i$ are replaced by the expansion coefficients $c_i({\bf G})$ and the orthonormality constraint only depends on the wavefunctions, not the nuclear positions. The equations of motion for the Car-Parrinello method are derived from this specific extended Lagrangian

$${\cal L}_{\rm CP} = \mu \sum_{i} \sum_{\bf G} \vert \dot c_{i}({\... ...M_{I} \dot {\bf R}_{I}^2 - E_{\rm KS} \left[ \{\bf G\}, \{{\bf R}_{I}\} \right]$$

$$+ \sum_{ij} \Lambda_{ij} \left( \sum_{\bf G} c^\star_{i}({\bf G}) c_{j}({\bf G}) -\delta_{ij} \right) \enspace ,$$ (96)

where $\mu$ is the fictitious electron mass, and $M_{I}$ are the masses of the nuclei. Because of the expansion of the Kohn-Sham orbitals in plane waves, the orthonormality constraint does not depend on the nuclear positions. The Euler-Lagrange equations derived from Eq.( 96) are

$$\mu \ddot c_{i}({\bf G}) = -{\partial E_{\rm KS} \over \partial c_{i}^\star({\bf G}) } + \sum_{j} \Lambda_{ij} c_{j}({\bf G})$$ (97)

$$M_{I} \ddot {\bf R}_{I} = -{\partial E_{\rm KS} \over \partial {\bf R}_{I}}.$$ (98)

The two sets of equations are coupled through the Kohn-Sham energy functional and special care has to be taken for the integration because of the orthonormality constraint. The velocity Verlet integration algorithm for the Car-Parrinello equations is defined as follows

$${\dot{\bf\tilde R}_{I}}(t + \delta t) = \dot{\bf R}_{I}(t) + {\delta t \over 2 M_{I}} {\bf F}_{I}(t)$$

$${\bf R}_{I}(t + \delta t) = {\bf R}_{I}(t) + \delta t \: {\dot{\bf\tilde R}_{I}}(t + \delta t)$$

$${\dot{\bf\tilde c}_{I}}(t + \delta t) = \dot{\bf c}_{I}(t) + {\delta t \over 2 \mu} {\bf f}_{i}(t)$$

$${\bf\tilde c}_{i}(t + \delta t) = {\bf c}_{i}(t) + \delta t \: {\dot{\bf\tilde c}_{i}}(t + \delta t)$$

$${\bf c}_{i}(t + \delta t) = {\bf\tilde c}_{i}(t + \delta t) + \sum_{j} {\bf X}_{ij} \: {\bf c}_{j}(t)$$

$${\bf F}_{I}(t + \delta t)$$ calculate

$${\bf f}_{i}(t + \delta t)$$ calculate

$${\dot{\bf R}_{I}}(t + \delta t) = {\dot{\bf\tilde R}_{I}}(t + \delta t) + {\delta t \over 2 M_{I}} {\bf F}_{I}(t + \delta t)$$

$${\dot{\bf c'}_{i}}(t + \delta t) = {\dot{\bf\tilde c}_{i}}(t + \delta t) + {\delta t \over 2 \mu} {\bf f}_{i}(t + \delta t)$$

$${\dot{\bf c}_{i}}(t + \delta t) = {\dot{\bf c'}_{i}}(t + \delta t) + \sum_{j} {\bf Y}_{ij} \: {\bf c}_{j}(t + \delta t) \enspace ,$$

where ${\bf R}_{I}(t)$ and ${\bf c}_{i}(t)$ are the atomic positions of particle $I$ and the Kohn-Sham orbital $i$ at time $t$ respectively. Here, ${\bf F}_{I}$ are the forces on atom $I$ , and ${\bf f}_{i}$ are the forces on Kohn-Sham orbital $i$ . The matrices ${\bf X}$ and ${\bf Y}$ are directly related to the Lagrange multipliers by

$${\bf X}_{ij} = {\delta t^2 \over 2 \mu} \Lambda^{\rm p}_{ij}$$ (99)

$${\bf Y}_{ij} = {\delta t \over 2 \mu} \Lambda^{\rm v}_{ij} \enspace .$$ (100)

Notice that in the RATTLE algorithm the Lagrange multipliers to enforce the orthonormality for the positions $\Lambda^{\rm p}$ and velocities $\Lambda^{\rm v}$ are treated as independent variables. Denoting with ${\bf C}$ the matrix of wavefunction coefficients $c_{i}({\bf G})$ , the orthonormality constraint can be written as

$${\bf C}^{\dagger}(t + \delta t) {\bf C}(t + \delta t) - {\bf I} =$$ 0 (101)

$$\left[ {\bf\tilde C} + {\bf X} {\bf C} \right]^{\dagger} \left[ {\bf\tilde C} + {\bf X} {\bf C} \right] - {\bf I} =$$ 0 (102)

$${\bf\tilde C}^{\dagger} {\bf\tilde C} + {\bf X} {\bf\tilde C}^{\d... ...}^{\dagger}{\bf\tilde C} {\bf X}^{\dagger} + {\bf X}{\bf X}^{\dagger} - {\bf I} =$$ 0 (103)

$${\bf X}{\bf X}^{\dagger} + {\bf X}{\bf B} + {\bf B}^{\dagger}{\bf X}^{\dagger} = {\bf I} - {\bf A} \enspace ,$$ (104)

where the new matrices ${\bf A}_{ij} = {\bf\tilde c}^\dagger_{i} (t + \delta t) {\bf\tilde c}_{j} (t + \delta t)$ and ${\bf B}_{ij} = {\bf c}^\dagger_{i} (t) {\bf\tilde c}_{j} (t + \delta t)$ have been introduced in Eq. (104). The unit matrix is denoted by the symbol ${\bf I}$ . By noting that ${\bf A} = {\bf I} + {\cal O}(\delta t^2)$ and ${\bf B} = {\bf I} + {\cal O}(\delta t)$ , Eq. (104) can be solved iteratively using

$${\bf X}^{(n+1)} = {1 \over 2} \left[ {\bf I} - {\bf A} + {\bf X}^{(n)} \left( {\bf ... ...{\bf I} - {\bf B} \right) {\bf X}^{(n)} - \left({\bf X}^{(n)} \right)^2 \right]$$ (105)

and starting from the initial guess

$${\bf X}^{(0)} = {1 \over 2} ({\bf I} - {\bf A}) \enspace .$$ (106)

In Eq. (105) it has been made use of the fact that the matrices ${\bf X}$ and ${\bf B}$ are real and symmetric, which follows directly from their definitions. Eq. (105) can usually be iterated to a tolerance of $10^{-6}$ within a few iterations.

The rotation matrix ${\bf Y}$ is calculated from the orthogonality condition on the orbital velocities

$$\dot {\bf c}^\dagger_{i}(t + \delta t) {\bf c}_{j}(t + \delta t) + {\bf c}^\dagger_{i}(t + \delta t) \dot {\bf c}_{j}(t + \delta t) = 0.$$ (107)

Applying Eq. (107) to the trial states $\dot {\bf C}' + {\bf Y} {\bf C}$ yields a simple equation for ${\bf Y}$

$${\bf Y} = - {1 \over 2} ({\bf Q} + {\bf Q}^\dagger),$$ (108)

where ${\bf Q}_{ij} = {\bf c}^\dagger_{i}(t + \delta t) \dot {\bf c'}^\dagger_{i}(t + \delta t)$ . The fact that ${\bf Y}$ can be obtained without iteration means that the velocity constraint condition Eq. (107) is satisfied exactly at each time step.