Total energy

In terms of the set of NGWFs $\{ \phi_{\alpha} \}$, the single-particle density matrix in the co-ordinate representation is expressed as

$$\rho(\mathbf{r},\mathbf{r}') = \phi_{\alpha}(\mathbf{r})K^{\alpha\beta}\phi_{\beta}(\mathbf{r}'),$$ (5)

where the density kernel, $K^{\alpha\beta}$, is the matrix representation of the density matrix in terms of the contravariant duals of the NGWFs [15,16]. The NGWFs are real, as we are concerned only with calculations at the $\Gamma$-point. Summation over repeated Greek indices is assumed throughout.

The charge density is given by,

$$n(\mathbf{r}) = 2\rho(\mathbf{r},\mathbf{r}) = 2K^{\alpha\be... ...hbf{r}) \equiv 2K^{\alpha\beta}\rho_{\alpha\beta}(\mathbf{r}),$$ (6)

where $\rho_{\alpha\beta}(\mathbf{r})$ is defined as,

$$\rho_{\alpha\beta}(\mathbf{r}) \equiv \phi_{\alpha}(\mathbf{... ...3}-1} \rho_{\alpha\beta}(\mathbf{r}_{XYZ})B_{XYZ}(\mathbf{r}).$$ (7)

This quantity is expanded in terms of the basis functions of the fine grid, $B_{XYZ}(\mathbf{r})$ (see Appendix A) because the density needs to be represented in terms of a basis with twice the cut-off frequency of the NGWFs.

Neglecting the contribution due to the fixed background of ionic charges, the total energy, $E[n]$, of the system is the sum of the kinetic energy ( $E_{\mathrm{K}}$), the Hartree energy ( $E_{\mathrm{H}}$), the local pseudopotential energy ( $E_{\mathrm{loc}}$), the non-local pseudopotential energy ( $E_{\mathrm{nl}}$), and the exchange-correlation energy ( $E_{\mathrm{xc}}$).

The kinetic energy is given by

$$E_{\mathrm{K}}[n] = - \int_{V} d\mathbf{r}' \: \left[ \nabla^{2}_{\mathbf{r}} \rho(\mathbf{r},\mathbf{r}')\right]_{\mathbf{r}=\mathbf{r}'}$$

$$= 2K^{\alpha\beta}\langle\phi_{\beta}\vert\hat{T}\vert\phi_{\alpha}\rangle,$$ (8)

where $\hat{T}$ is related to the Laplacian operator by $\hat{T} = -\frac{1}{2}\nabla^{2}$. This is most conveniently applied in reciprocal space, where it is diagonal. We see that the kinetic energy involves matrix elements of the form $\langle\phi_{\beta}\vert\hat{T}\vert\phi_{\alpha}\rangle$. Thus, one way of proceeding would be to FFT $\phi_{\alpha}$ in the whole simulation cell, apply the Laplacian, FFT back and perform the integral over real space as a summation over grid points as suggested by equation (52). If $N$ is the number of grid points in the simulation cell (which is proportional to system size), this procedure has a computational cost that scales as $O(N\log N)$ for each NGWF and hence as $O(N^{2}\log N)$ overall. We will see later how linear-scaling can be achieved for all the terms that compose $E[n]$.

The Hartree energy is given by,

$$E_{\mathrm{H}}[n] = \frac{1}{2}\int\int d\mathbf{r} \: d\mathbf{r}' \frac{n(\mathbf{r... ...}'\vert} = \frac{1}{2}\int d\mathbf{r} V_{\mathrm{H}}(\mathbf{r}) n(\mathbf{r})$$

$$= K^{\alpha\beta}\langle\phi_{\beta}\vert V_{\mathrm{H}}\vert\phi_{\alpha}\rangle,$$ (9)

where the Hartree potential $V_{\mathrm{H}}(\mathbf{r})$ is,

$$V_{\mathrm{H}}(\mathbf{r}) = \int d\mathbf{r}' \: \frac{n(\mathbf{r}')}{\vert\mathbf{r}-\mathbf{r}'\vert}$$

$$= \sum_{X=0}^{2N_{1}-1}\sum_{Y=0}^{2N_{2}-1}\sum_{Z=0}^{2N_{3}-1} V_{\mathrm{H}}(\mathbf{r}_{XYZ})B_{XYZ}(\mathbf{r}).$$ (10)

We see that $V_{\mathrm{H}}$ is the result of a convolution between the charge density and the Coulomb potential, and is thus represented by the fine grid basis functions.

$V_{\mathrm{H}}$ is calculated on the fine grid in reciprocal space, where the real space convolution becomes a simple product, and fast Fourier transformed back to real space. $\phi_{\alpha}$ is Fourier interpolated onto the fine grid and its product with $V_{\mathrm{H}}$ taken. The result is Fourier filtered to the standard grid, and the matrix elements $\langle\phi_{\beta}\vert\left[ V_{\mathrm{H}}\phi_{\alpha} \right]_{D} \rangle$ calculated by summation over the grid points of the standard grid, where the subscript $D$ shows that we only need to consider frequency components represented by the basis functions $\{ D(\mathbf{r}) \}$ as shown in equation (52).

The local pseudopotential energy is given by the integral

$$E_{\mathrm{loc}}[n] = \int d\mathbf{r} \: V_{\mathrm{loc}}(\... ...\left[ V_{\mathrm{loc}}(\mathbf{r}) \right]_{B} n(\mathbf{r}),$$ (11)

where the second equality makes use of equation (56), and the subscript $B$ shows that only frequency components represented by the fine grid basis functions, $B_{XYZ}(\mathbf{r})$, need to be considered. Thus, it may be calculated in the same way as the Hartree energy:

$$E_{\mathrm{loc}}[n] = 2K^{\alpha\beta}\langle\phi_{\beta}\vert \left[ V_{\mathrm{loc}} \right]_{B}\vert\phi_{\alpha}\rangle$$

$$= 2K^{\alpha\beta}\langle\phi_{\beta}\vert \left( \left[ V_{\mathrm{loc}} \right]_{B} \phi_{\alpha} \right)_{D} \rangle.$$ (12)

The non-local pseudopotential energy is given by

$$E_{\mathrm{nl}} = 2K^{\alpha\beta}\langle\phi_{\beta}\vert\hat{V}_{\mathrm{nl}}\vert\phi_{\alpha}\rangle,$$ (13)

where the non-local pseudopotential operator $\hat{V}_{\mathrm{nl}}$ that we use is in Kleinman-Bylander form [17],

$$\hat{V}_{\mathrm{nl}} = \sum_{I}\sum_{lm(I)}\frac{\vert\delt... ...(I)}\vert\delta \hat{V}_{l}^{(I)}\vert\Psi_{lm}^{(I)}\rangle},$$ (14)

where the first summation is over the atoms $I$, and the second is over the angular momentum components for each atom. $\delta \hat{V}_{l}^{(I)}$ is the angular momentum dependent part of the atomic pseudopotential and the $\Psi_{lm}^{(I)}$ are its associated pseudo-orbitals. We see that the matrix elements of the non-local potential involve the calculation of integrals of the form $\langle\phi_{\beta}\vert\delta \hat{V}_{l}^{(I)}\Psi_{lm}^{(I)}\rangle = \langl... ..._{\beta}\vert \left[ \delta \hat{V}_{l}^{(I)}\Psi_{lm}^{(I)} \right]_{D}\rangle$, which may be performed by summation over the standard real space grid.

The exchange-correlation energy within the LDA is given by approximating the integral over the exchange-correlation energy density, $\epsilon_{\mathrm{xc}}(n(\mathbf{r}))$, by a summation over the points of the fine grid:

$$E_{\mathrm{xc}}[n] = \int_{V} d\mathbf{r} \: \epsilon_{\mathrm{xc}}(n(\mathbf{r})) n(\mathbf{r})$$

$$\simeq \frac{V}{8N_{1}N_{2}N_{3}}\sum_{X=0}^{2N_{1}-1}\sum_{Y=0}^{2N_{2}... ...=0}^{2N_{3}-1} \epsilon_{\mathrm{xc}}(n(\mathbf{r}_{XYZ})) n(\mathbf{r}_{XYZ}).$$ (15)

This is approximate as the fine grid cannot represent the exchange-correlation energy density without some aliasing. The errors associated with this approximation are expected to be similar to those encountered in traditional plane-wave codes [1].

We may write the total energy as,

$$E[n] = 2K^{\alpha\beta}H_{\beta\alpha} - E_{\mathrm{DC}}[n],$$ (16)

where the matrix elements of the Kohn-Sham Hamiltonian are given by,

$$H_{\beta\alpha} = \langle\phi_{\beta}\vert\left[ \hat{T} + V... ... V_{\mathrm{xc}}(\mathbf{r})\right] \vert\phi_{\alpha}\rangle,$$ (17)

$V_{xc}(\mathbf{r}) = \frac{\delta E_{\mathrm{xc}}[n]}{\delta n(\mathbf{r})}$ is the exchange-correlation potential, and $E_{\mathrm{DC}}[n]$ is the double-counting correction,

$$E_{\mathrm{DC}}[n] = \frac{1}{2}\int d\mathbf{r} \: n(\mathb... ...(\mathbf{r}) V_{\mathrm{xc}}(\mathbf{r}) - E_{\mathrm{xc}}[n].$$ (18)

As we are dealing with a sparse system, the number of non-zero matrix elements will be proportional to the system size in the limit of large systems. The procedure outlined above for calculating each of these matrix elements has a computational cost that scales as $O(N\log N)$ for each element. This is because FFTs of NGWFs are performed over the entire simulation cell. Thus the cost of computing all non-zero matrix elements scales as $O(N^{2}\log N)$.