7.1 Total energy and Hamiltonian

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

$$S_{\alpha \beta} = \int {\mathrm d}{\bf r}~ \phi_{\beta}({\bf r}) \phi_{\alpha}({\bf r}) .$$ (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

$${\cal S}_{\alpha , n \ell m ;\beta , n' \ell' m'} = \langle ... ..._{\alpha , n \ell m} \vert \chi_{\beta , n' \ell' m'} \rangle$$ (7.4)

so that the overlap matrix elements are given by

$$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)}$$ (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

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

in which

$$T_{\alpha \beta} = -\frac{1}{2} \int {\mathrm d}{\bf r}~ \phi_{\alpha}({\bf r}) \nabla_{\bf r}^2 \phi_{\beta}({\bf r})$$ (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,

$$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)}$$ (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-space^7.1to obtain ${\tilde n}({\bf G})$. The Hartree energy is then given by

$$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}$$ (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

$$V_{\mathrm H}({\bf r}) = \int {\mathrm d}{\bf r'}~ \frac{n({\bf r'})}{\left\vert {\bf r} - {\bf r'} \right\vert}$$ (7.10)

but is calculated in reciprocal-space as

$${\tilde V}_{\mathrm H}({\bf G}) = \frac{4 \pi {\tilde n}({\bf G})}{\Omega_{\mathrm {cell}} G^2}$$ (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

$$E_{\mathrm {xc}} = \delta \omega \sum_{\bf r} n({\bf r}) \epsilon_{\mathrm {xc}}(n({\bf r}))$$ (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

$$V_{\mathrm {xc}}({\bf r}) = \left[ \frac{\mathrm d}{{\mathrm... ...\epsilon_{\mathrm {xc}}(n) \right\} \right]_{n = n({\bf r})} .$$ (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

$${\tilde V}_{\mathrm {ps,loc}}({\bf G}) = \sum_s {\tilde v}_{\mathrm {ps,loc}}^s ({\bf G}) {\mbox{\Bbb S}}^s({\bf G})$$ (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

$${\mbox{\Bbb S}}^s({\bf G}) = \sum_{\alpha} \exp[-{\mathrm i}{\bf G} \cdot {\bf r}_{\alpha}^s]$$ (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

$$E_{\mathrm{ps,loc}} = \int {\mathrm d}{\bf r}~ V_{\mathrm{ps,loc}}({\bf r}) n({\bf r})$$

$$= \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$$

$$= 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

$$V_{\mathrm{loc},\alpha \beta} = \int {\mathrm d}{\bf r}~ \ph... ...({\bf r}) V_{\mathrm{ps,loc}}({\bf r}) \phi_{\beta}({\bf r}) .$$ (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

$$E_{\mathrm{ps,NL}} = \int {\mathrm d}{\bf r}~{\mathrm d}{\bf... ...\bf r}) = 2 K^{\alpha \beta} V_{{\mathrm {NL}},\beta \alpha} .$$ (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

$$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)}$$ (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.