next up previous contents index
Next: Metals; Free Energy Functional Up: Calculating the Electronic Structure Previous: Exchange and Correlation Energy   Contents   Index


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

$\displaystyle {\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]$      
$\displaystyle + \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
$\displaystyle \mu \ddot c_{i}({\bf G})$ $\displaystyle =$ $\displaystyle -{\partial E_{\rm KS} \over \partial c_{i}^\star({\bf G}) }
+ \sum_{j} \Lambda_{ij} c_{j}({\bf G})$ (97)
$\displaystyle M_{I} \ddot {\bf R}_{I}$ $\displaystyle =$ $\displaystyle -{\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
$\displaystyle {\dot{\bf\tilde R}_{I}}(t + \delta t)$ $\displaystyle =$ $\displaystyle \dot{\bf R}_{I}(t) + {\delta t \over 2 M_{I}}
{\bf F}_{I}(t)$  
$\displaystyle {\bf R}_{I}(t + \delta t)$ $\displaystyle =$ $\displaystyle {\bf R}_{I}(t) + \delta t \:
{\dot{\bf\tilde R}_{I}}(t + \delta t)$  
$\displaystyle {\dot{\bf\tilde c}_{I}}(t + \delta t)$ $\displaystyle =$ $\displaystyle \dot{\bf c}_{I}(t) + {\delta t \over 2 \mu}
{\bf f}_{i}(t)$  
$\displaystyle {\bf\tilde c}_{i}(t + \delta t)$ $\displaystyle =$ $\displaystyle {\bf c}_{i}(t) + \delta t \:
{\dot{\bf\tilde c}_{i}}(t + \delta t)$  
$\displaystyle {\bf c}_{i}(t + \delta t)$ $\displaystyle =$ $\displaystyle {\bf\tilde c}_{i}(t + \delta t) + \sum_{j} {\bf X}_{ij} \:
{\bf c}_{j}(t)$  
calculate   $\displaystyle {\bf F}_{I}(t + \delta t)$  
calculate   $\displaystyle {\bf f}_{i}(t + \delta t)$  
$\displaystyle {\dot{\bf R}_{I}}(t + \delta t)$ $\displaystyle =$ $\displaystyle {\dot{\bf\tilde R}_{I}}(t + \delta t) + {\delta t \over 2 M_{I}}
{\bf F}_{I}(t + \delta t)$  
$\displaystyle {\dot{\bf c'}_{i}}(t + \delta t)$ $\displaystyle =$ $\displaystyle {\dot{\bf\tilde c}_{i}}(t + \delta t) + {\delta t \over 2 \mu}
{\bf f}_{i}(t + \delta t)$  
$\displaystyle {\dot{\bf c}_{i}}(t + \delta t)$ $\displaystyle =$ $\displaystyle {\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
$\displaystyle {\bf X}_{ij}$ $\displaystyle =$ $\displaystyle {\delta t^2 \over 2 \mu} \Lambda^{\rm p}_{ij}$ (99)
$\displaystyle {\bf Y}_{ij}$ $\displaystyle =$ $\displaystyle {\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
$\displaystyle {\bf C}^{\dagger}(t + \delta t) {\bf C}(t + \delta t) - {\bf I}$ $\displaystyle =$ 0 (101)
$\displaystyle \left[ {\bf\tilde C} + {\bf X} {\bf C} \right]^{\dagger}
\left[ {\bf\tilde C} + {\bf X} {\bf C} \right] - {\bf I}$ $\displaystyle =$ 0 (102)
$\displaystyle {\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}$ $\displaystyle =$ 0 (103)
$\displaystyle {\bf X}{\bf X}^{\dagger} + {\bf X}{\bf B} + {\bf B}^{\dagger}{\bf X}^{\dagger}$ $\displaystyle =$ $\displaystyle {\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
$\displaystyle {\bf X}^{(n+1)}$ $\displaystyle =$ $\displaystyle {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

$\displaystyle {\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

$\displaystyle \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}$

$\displaystyle {\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.


next up previous contents index
Next: Metals; Free Energy Functional Up: Calculating the Electronic Structure Previous: Exchange and Correlation Energy   Contents   Index
Costas Bekas 2008-07-04