next up previous contents
Next: 7.2 Energy gradients Up: 7. Computational implementation Previous: 7. Computational implementation   Contents

Subsections

7.1 Total energy and Hamiltonian

A quantity which is frequently required is the overlap matrix $S_{\alpha \beta}$ defined by

\begin{displaymath}
S_{\alpha \beta} = \int {\mathrm d}{\bf r}~
\phi_{\beta}({\bf r}) \phi_{\alpha}({\bf r}) .
\end{displaymath} (7.3)

The overlap matrix elements between the spherical-wave basis functions can be calculated analytically (section 5.4), and are denoted ${\cal S}_{\alpha, n \ell m ; \beta, n' \ell' m'}$ where
\begin{displaymath}
{\cal S}_{\alpha , n \ell m ;\beta , n' \ell' m'} =
\langle ...
..._{\alpha , n \ell m} \vert \chi_{\beta , n' \ell' m'} \rangle
\end{displaymath} (7.4)

so that the overlap matrix elements are given by
\begin{displaymath}
S_{\alpha \beta} = \sum_{n \ell m, n' \ell' m'} c^{n \ell m}...
...pha , n \ell m ;\beta , n' \ell' m'} c^{n' \ell' m'}_{(\beta)}
\end{displaymath} (7.5)

recalling that the support functions may be assumed real in the case of $\Gamma$-point Brillouin zone sampling.

7.1.1 Kinetic energy

The kinetic energy of the non-interacting Kohn-Sham system is given by

\begin{displaymath}
T_{\mathrm s}^{\mathrm J} = - \int {\mathrm d}{\bf r'}~
\lef...
...ight]_{{\bf r}={\bf r'}}
= 2 K^{\alpha \beta} T_{\beta \alpha}
\end{displaymath} (7.6)

in which
\begin{displaymath}
T_{\alpha \beta} = -\frac{1}{2} \int {\mathrm d}{\bf r}~
\phi_{\alpha}({\bf r}) \nabla_{\bf r}^2 \phi_{\beta}({\bf r})
\end{displaymath} (7.7)

are the matrix elements of the kinetic energy operator in the representation of the support functions. Since all of the matrix elements between the spherical-wave basis functions can be calculated analytically,
\begin{displaymath}
T_{\alpha \beta} = \sum_{n \ell m, n' \ell' m'} c^{n \ell m}...
...pha , n \ell m ;\beta , n' \ell' m'} c^{n' \ell' m'}_{(\beta)}
\end{displaymath} (7.8)

where ${\cal T}$ denotes matrix elements of the kinetic energy operator between spherical-wave basis functions.

7.1.2 Hartree energy and potential

The Hartree and exchange-correlation terms are calculated by determining the electronic density on a real-space grid $n({\bf r})$, and Fast Fourier Transforms (FFTs) are used to transform between real- and reciprocal-space7.1to obtain ${\tilde n}({\bf G})$. The Hartree energy is then given by

\begin{displaymath}
E_{\mathrm H} = \frac{1}{2} \int {\mathrm d}{\bf r}~{\mathrm...
...t= 0}
\frac{\left\vert {\tilde n}({\bf G}) \right\vert^2}{G^2}
\end{displaymath} (7.9)

where $\Omega_{\mathrm{cell}}$ is the volume of the supercell and the (infinite) ${\bf G} = 0$ term is omitted because the system is charge neutral overall. This term is therefore cancelled by similar terms in the ion-ion and electron-ion interaction energies. The Hartree potential in real-space is given by
\begin{displaymath}
V_{\mathrm H}({\bf r}) = \int {\mathrm d}{\bf r'}~ \frac{n({\bf r'})}{\left\vert
{\bf r} - {\bf r'} \right\vert}
\end{displaymath} (7.10)

but is calculated in reciprocal-space as
\begin{displaymath}
{\tilde V}_{\mathrm H}({\bf G}) = \frac{4 \pi {\tilde n}({\bf G})}{\Omega_{\mathrm {cell}} G^2}
\end{displaymath} (7.11)

and then transformed back into real-space by a FFT.

7.1.3 Exchange-correlation energy and potential

Having calculated the electron density on the grid points, the exchange-correlation energy is obtained by summing over those grid points

\begin{displaymath}
E_{\mathrm {xc}} = \delta \omega \sum_{\bf r} n({\bf r}) \epsilon_{\mathrm {xc}}(n({\bf r}))
\end{displaymath} (7.12)

in the local density approximation. $\delta \omega$ is the volume of the supercell divided by the number of grid points. The exchange-correlation potential is similarly calculated at each grid point as
\begin{displaymath}
V_{\mathrm {xc}}({\bf r}) = \left[ \frac{\mathrm d}{{\mathrm...
...\epsilon_{\mathrm {xc}}(n) \right\} \right]_{n = n({\bf r})} .
\end{displaymath} (7.13)

In practice, the values of $\epsilon_{\mathrm {xc}}(n)$ and $\frac{\mathrm d}{{\mathrm d}n} \left[
n \epsilon_{\mathrm {xc}}(n) \right]$ are tabulated for various values of the electronic density $n$ and then interpolated during the calculation.

7.1.4 Local pseudopotential

Like the Hartree potential, the local pseudopotential is also calculated in reciprocal-space as

\begin{displaymath}
{\tilde V}_{\mathrm {ps,loc}}({\bf G}) = \sum_s {\tilde v}_{\mathrm {ps,loc}}^s
({\bf G}) {\mbox{\Bbb S}}^s({\bf G})
\end{displaymath} (7.14)

where the summation is over ionic species $s$, ${\tilde v}_{\mathrm {ps,loc}}^s
({\bf G})$ is the local pseudopotential for an isolated ion of species $s$ in reciprocal-space and ${\mbox{\Bbb S}}^s({\bf G})$ is the structure factor for species $s$ defined by
\begin{displaymath}
{\mbox{\Bbb S}}^s({\bf G}) = \sum_{\alpha} \exp[-{\mathrm i}{\bf G} \cdot {\bf r}_{\alpha}^s]
\end{displaymath} (7.15)

where the sum is over all ions $\alpha $ of species $s$ with positions ${\bf r}_{\alpha}^s$. We note that in general the calculation of the structure factor is an ${\cal O}(N^2)$ operation, but since it only has to be calculated once for each atomic configuration, it is not a limiting factor of the overall calculation at this stage. Within the quantum chemistry community, work on generalised multipole expansions and new algorithms [156,157,158,159,160,161,162] has led to the development of methods to calculate Coulomb interaction matrix elements which scale linearly with system-size. The local pseudopotential energy can be calculated in reciprocal-space as
$\displaystyle E_{\mathrm{ps,loc}}$ $\textstyle =$ $\displaystyle \int {\mathrm d}{\bf r}~ V_{\mathrm{ps,loc}}({\bf r}) n({\bf r})$  
  $\textstyle =$ $\displaystyle \sum_{{\bf G} \not= 0} \sum_s
{\tilde v}_{\mathrm{ps,loc}}^s({\bf...
...b S}}^s({\bf G}) {\tilde n}^{\ast}
({\bf G})+\sum_s N_s E_{\mathrm{ps, core}}^s$  
  $\textstyle =$ $\displaystyle 2 K^{\alpha \beta} V_{\mathrm{loc},\beta \alpha}$ (7.16)

where $E_{\mathrm{ps, core}}^s$ is the pseudopotential core energy, and $N_s =
{\mbox{\Bbb S}}^s({\bf G}=0)$ the number of ions of species $s$. The matrix elements $V_{\mathrm{loc},\alpha \beta}$ are defined by
\begin{displaymath}
V_{\mathrm{loc},\alpha \beta} = \int {\mathrm d}{\bf r}~ \ph...
...({\bf r})
V_{\mathrm{ps,loc}}({\bf r}) \phi_{\beta}({\bf r}) .
\end{displaymath} (7.17)

The Hartree potential and local pseudopotential can be summed and then transformed back together into real-space and added to the exchange-correlation potential to obtain the local part of the Kohn-Sham potential in real-space.

We note that the FFT is not strictly an ${\cal O}(N)$ operation but an ${\cal O}(N \log_m N)$ operation (where $m$ is some small number which depends upon the prime factors of the number of grid points), but in practice (section 9.2) this scaling is not observed.

7.1.5 Non-local pseudopotential

The non-local pseudopotential energy is given by

\begin{displaymath}
E_{\mathrm{ps,NL}} = \int {\mathrm d}{\bf r}~{\mathrm d}{\bf...
...\bf r}) =
2 K^{\alpha \beta} V_{{\mathrm {NL}},\beta \alpha} .
\end{displaymath} (7.18)

The matrix elements of the non-local pseudopotential in the representation of the support functions $V_{{\mathrm {NL}},\alpha \beta}$ are calculated by summing over all ions whose cores overlap the support regions of $\phi_{\alpha}$ and $\phi_{\beta}$, and using the method described in section 5.6.2 to calculate the spherical-wave basis function matrix elements ${\cal V}_{\alpha , n \ell m ;\beta , n' \ell' m'}$ analytically. The result is therefore
\begin{displaymath}
V_{{\mathrm {NL}},\alpha \beta} = \sum_{n \ell m, n' \ell' m...
...pha , n \ell m ;\beta , n' \ell' m'} c^{n' \ell' m'}_{(\beta)}
\end{displaymath} (7.19)

which is of exactly the same form as the kinetic energy, so that in practice the basis function matrix elements for the kinetic energy and non-local pseudopotential are summed and the two contributions to the energy combined.
next up previous contents
Next: 7.2 Energy gradients Up: 7. Computational implementation Previous: 7. Computational implementation   Contents
Peter Haynes