Charge density and total electronic energy with Non-orthogonal Generalised Wannier functions

Linear-scaling DFT calculations are aimed at large systems, and in particular, large unit cells. Thus in this work we will be concerned with calculations only at the $\Gamma$-point, i.e. $\mathbf{k}=0$. This means that the Bloch bands and therefore the NGWFs can be chosen to be real. We can also drop the dependence of the NGWFs on $\mathbf{R}$, so that $\phi_{\alpha\mathbf{R}}(\mathbf{r})=\phi_{\alpha}(\mathbf{r})$.

Our basis set is the set of periodic bandwidth limited delta functions that are centred on the points $\mathbf{r}_{KLM}$ of a regular real-space grid:

$$D_{KLM}(\mathbf{r})= \frac{1}{N_1 N_2 N_3} \sum_{P= -J_1}^{... ...2 + R\mathbf{B}_3 ) \cdot(\mathbf{r}-\mathbf{r}_{KLM}) } \,\,,$$ (5)

where $\mathbf{B}_1$ is one of the reciprocal lattice vectors of the simulation cell. $N_1$ is the number of grid points in the direction of direct lattice vector $\mathbf{A}_1$, and $N_1=2J_1+1$. The delta function basis is equivalent to the plane-waves that can be represented by the real-space grid since it is related to them via a unitary transformation. An important property of the basis set is that the projection of a function $f(\mathbf{r})$ on $D_{KLM}(\mathbf{r})$ is

$$\int_V D_{KLM}(\mathbf{r})f(\mathbf{r}) d\mathbf{r}= W f_{D}(\mathbf{r}_{KLM})$$ (6)

where $W$ is the volume per grid point and $f_{D}(\mathbf{r})$ is the result of bandwidth limiting the function $f(\mathbf{r})$ to the same plane-wave components as in (5).

We represent the NGWFs in the delta function basis by

$$\phi_{\alpha}(\mathbf{r})= \sum_{K=0}^{N_1-1} \sum_{L=0}^{N_2-1} \sum_{M=0}^{N_3-1} C_{KLM,\alpha} D_{KLM}(\mathbf{r}) \,\,,$$ (7)

and in the plane-wave basis by

$$\phi_{\alpha}(\mathbf{r})=\frac{1}{V} \sum_{P=-J_1}^{J_1} \... ...B}_1+ Q \mathbf{B}_2 + R \mathbf{B}_3) \cdot\mathbf{r}} \,\, .$$ (8)

where it is straightforward to show that the amplitudes $\tilde{\phi}_\alpha(P \mathbf{B}_1+ Q \mathbf{B}_2 + R\mathbf{B}_3 )$ are the result of a discrete Fourier transform on the delta function expansion coefficients $C_{KLM,\alpha}$.

In (7) the sum over the $K$, $L$ and $M$ indices formally goes over the grid points of a regular grid that extends over the whole simulation cell. From now on however, we will restrict all NGWFs to have contributions only from delta functions centred inside a predefined spherical region. This spherical region is in general different for each NGWF. Thus we impose on (7) the condition:

$$C_{KLM,\alpha}=0 \,\,\,\, \mbox{if }\mathbf{r}_{KLM} \,\,\, \mbox{does not belong to the sphere of }\phi_{\alpha} \, .$$ (9)

This of course does not affect the form or the applicability of equation (8).

The charge density of equation (2) with our NGWFs becomes (from now on we will use the summation convention for repeated Greek indices)

$$n(\mathbf{r})= 2\rho(\mathbf{r},\mathbf{r}) = 2\phi_{\alpha}(\mathbf{r}) K^{\alpha \beta} \phi_{\beta}(\mathbf{r}) = 2 K^{\alpha \beta} \rho_{\alpha \beta}(\mathbf{r})$$

$$= 2 \sum_{X=0}^{(2N_1-1)} \sum_{Y=0}^{(2N_2-1)} \sum_{Z=0}^{(2N_3-1... ...^{\alpha \beta} \rho_{\alpha \beta}(\mathbf{r}_{XYZ}) B_{XYZ}(\mathbf{r}) \,\,,$$ (10)

which involves the fine grid delta functions $B_{XYZ}(\mathbf{r})$ that are defined in a similar way to the $D_{KLM}(\mathbf{r})$ of equation (5) but include up to twice the maximum wavevector of $D_{KLM}(\mathbf{r})$ in every reciprocal lattice vector direction (see also Appendix A). This is necessary because a product of two $D_{KLM}(\mathbf{r})$ delta functions is a linear combination of fine grid delta functions $B_{XYZ}(\mathbf{r})$, a result reminiscent of the Gaussian function product rule [30].

The expressions for the various contributions to the total electronic energy with the NGWFs are simple to derive from (10). The total energy is the sum of the kinetic energy $E_K$, the Hartree energy $E_H$, the local pseudopotential energy $E_{loc}$, the non-local pseudopotential energy $E_{nl}$ and the exchange and correlation energy $E_{xc}$

$$E[n]= E_{K}[n] + E_{H}[n]+ E_{loc}[n] + E_{nl}[n] + E_{xc}[n] \,\, .$$ (11)

The kinetic energy is written as a trace of the product of the density kernel and of the matrix elements of the kinetic energy operator $\hat{T}=-(1/2)\nabla^2$

$$E_{K}[n] = 2 K^{\alpha \beta} \langle \phi_{\beta} \vert \hat{T} \vert \phi_{\alpha} \rangle .$$ (12)

To compute these matrix elements we can apply $\hat{T}$ to the plane-wave representation (8) of $\phi_{\alpha}(\mathbf{r})$ and then evaluate the integral in real-space where it is equal to a discrete sum over grid points where $\hat{T}\phi_{\alpha}(\mathbf{r})$ obviously plays the role of $f_{D}(\mathbf{r})$ of equation (6).

Calculation of the Hartree energy requires first the Hartree potential. From equation (10) we see that the charge density is a fine grid delta function expansion, thus the same should be true for the Hartree potential, which is a convolution of the charge density with the Coulomb potential. Therefore, $V_{H}(\mathbf{r})$ can be written as a linear combination of fine grid delta functions and extends over the whole simulation cell:

$$V_{H}(\mathbf{r})= \sum_{X=0}^{(2N_1-1)} \sum_{Y=0}^{(2N_2-... ...^{(2N_3-1)} V_{H}(\mathbf{r}_{XYZ} ) B_{XYZ}(\mathbf{r}) \,\,.$$ (13)

The Hartree energy is

$$E_{H}[n]= \frac{1}{2}\int V_{H}(\mathbf{r}) n(\mathbf{r}) d\... ...langle \phi_{\beta} \vert V_{H} \vert \phi_{\alpha} \rangle .$$ (14)

This quantity can be calculated as a discrete summation on the fine grid of the product of $V_{H}(\mathbf{r})$ with $n(\mathbf{r})$ or equivalently as a trace of the product of the density kernel and the potential matrix elements. The local potential matrix elements are integrals that are identically equal to discrete sums on the regular grid provided of course that $(V_{H}\phi_{\alpha})_{D}(\mathbf{r})$ is first put on the regular grid.

The local pseudopotential energy is calculated in an entirely analogous manner to the Hartree energy and can be represented by equation (14) if we put $V_{loc}$ in place of $V_H$ and multiply it by a factor of 2 to take into account the lack of self-interaction in this case.

The non-local pseudopotential energy $E_{nl}[n]$ is the expectation value of the non-local potential operator $\hat{V}_{nl}$ in the Kleinman-Bylander form [31]:

$$\hat{V}_{nl}= \sum_A \sum_{lm(A)} \frac{\vert\delta \hat{V}^... ...\delta \hat{V}^{(A)}_{l} \vert \Psi^{(A)}_{lm} \rangle } \,\,,$$ (15)

where the $A$-summation runs over the atoms in the system and the $lm$-summation runs over the pseudo-atomic orbitals of a particular atom. The $\delta \hat{V}^{(A)}_{l}$ is an angular momentum dependent component of the non-local potential of a pseudo-Hamiltonian for a particular atom and the $\Psi^{(A)}_{lm}$ are the atomic pseudo-orbitals associated with it. In the NGWF representation the non-local potential energy is again expressed as a matrix trace:

$$E_{nl}[n]=2 K^{\alpha \beta} \langle \phi_{\beta} \vert \hat{V}_{nl} \vert \phi_{\alpha} \rangle \, .$$ (16)

The $\langle \phi_{\beta} \vert\hat{V}_{nl} \vert \phi_{\alpha} \rangle$ matrix elements require the calculation of overlap integrals $\langle \phi_{\alpha} \vert\delta \hat{V}^{(A)}_{l} \Psi^{(A)}_{lm} \rangle$ between the NGWFs and the non-local projectors $\vert\delta \hat{V}^{(A)}_{l} \Psi^{(A)}_{lm} \rangle$. These are simple to compute as discrete summations on the regular grid, starting from the plane-wave representation of the non-local projectors which is analogous to the plane-wave representation of the NGWFs in equation (8). These integrals need only be calculated when the sphere of function $\phi_{\alpha}(\mathbf{r})$ overlaps the core of atom $A$.

The exchange-correlation energy is obtained by approximating the exchange-correlation functional expression as a direct summation on the fine grid, which first involves the evaluation of a function $F(n(\mathbf{r}))$ whose particular form depends on our choice of exchange-correlation functional[2]:

$$E_{xc}[n]=\int_{V} F(n(\mathbf{r})) d \mathbf{r} \simeq \frac{V}{8N_1 N_2 N_3} \sum_{XYZ} F(n(\mathbf{r}_{XYZ}) ) \,\, .$$ (17)

This is the only approximation in integral evaluation in our method as all direct summations described up to now were exactly equal to the analytic integrals. However, in the case of the exchange-correlation energy, the exchange-correlation functionals usually contain highly non-linear expressions that can not be represented without any aliasing even when we use the delta functions of the fine grid. The resulting errors however will be of the same nature as in conventional plane-wave codes and therefore negligible [1].