The FFT box technique

Our discussion so far has demonstrated how to localise the NGWFs in real-space in a manner that ensures that they are composed by a number of delta functions that is constant with system-size. On the other hand, each delta function is expanded in the plane-waves that can be supported by the regular grid of the simulation cell. The number of these plane-waves is proportional to the size of the system. As a result, the cost of a calculation still scales cubically with system-size as in the traditional plane-wave case.

In order to reduce the computational cost we must restrict the number of plane-waves that contribute to each delta function so that it is independent of system-size. This is not such a straightforward matter as the restriction of the NGWFs in real-space. Several factors have to be taken into account the most important of which is the Hermiticity of operators in integrals between NGWFs and also a representation of operators that is consistent when they are acting on different NGWFs. We have investigated these matters in detail in the context of the evaluation of kinetic energy integrals with NGWFs in a previous paper [43] where we have proposed the `FFT box' technique as an accurate and efficient solution. It is shown there that all the imposed conditions are satisfied if the plane-waves that are used to expand the delta functions are restricted to belong to a miniature simulation cell which we call the `FFT box'. The FFT box must be large enough to contain any possible orientation of overlapping NGWFs, and must be a parallelepiped with a shape commensurate with the simulation cell. It should contain a regular grid which is a subset of the simulation cell grid and thus the origin of the grid of the FFT box should coincide with a particular grid point of the simulation cell. FFT techniques using smaller boxes have been used in the past to take advantage of localised functions: Pasquarello et al. [44,45] use them to efficiently deal with the augmentation charges that arise when using ultrasoft pseudopotentials; Hutter et. al [46] use them to FFT localised Gaussian orbitals. However, to the best of our knowledge, this is the first time that a small FFT box has been used to define a systematic basis set for the entire total energy calculation.

Figure 1: The FFT box as defined for the pair of functions $\phi _{\alpha }$ and $\phi _{\beta }$.

Figure 1 demonstrates in two dimensions the FFT box inside the simulation cell as defined for two overlapping NGWFs $\phi _{\alpha }$ and $\phi _{\beta }$. In the same Figure the regular grid is also visible and we can observe that a portion of it is included in the FFT box. For the subset of grid points inside the FFT box we can define delta functions as we did for the simulation cell. We represent these FFT box delta functions by $d_{klm}(\mathbf{r})$ and in general we follow the convention of using lowercase letters to denote quantities associated with the FFT box and uppercase letters for quantities associated with the simulation cell.

Based on the knowledge gained through the use of the FFT box to compute kinetic energy integrals[43], we extend here its use to the calculation of all the terms that constitute the total energy (section 2) and all the terms needed for its subsequent optimisation with respect to the density kernel and NGWF coefficients (section 3). For this purpose we need to be able to express or `project' a function from the simulation cell to the FFT box and back. Assuming the function is expressed entirely by delta functions $D_{KLM}(\mathbf{r})$ on grid points common to the simulation cell and the FFT box, this task is straighforward and all we have to do is to rewrite the function as a linear combination of the FFT box delta functions that lie on the same centres as the simulation cell delta functions. This task is expressed formally by the $\hat{P}(\alpha \beta)$ projection operator. When this operator acts on a function in the simulation cell, it maps it onto a function in the FFT box. This is demonstrated for the function of equation (7) by the following expression which can also be considered as a definition of $\hat{P}(\alpha \beta)$ (a more explicit expression for $\hat{P}(\alpha \beta)$ is given in Appendix B):

$$\hat{P}(\alpha \beta) \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})$$

$$= \sum_{k=0}^{n_1-1} \sum_{l=0}^{n_2-1} \sum_{m=0}^{n_3-1} C_{(k+K_{\alpha \beta})(l+L_{\alpha \beta})(m+M_{\alpha \beta})} d_{klm}(\mathbf{r})$$ (25)

while the adjoint operator $\hat{P}^{\dag }(\alpha \beta)$ can act on a function in the FFT box and turn it into a function in the simulation cell in an analogous way.

With this compact notation it is relatively straightforward to devise a way of calculating the total energy by using only the delta functions (and hence the plane-waves) periodic in the FFT box as the basis set for each NWGF. It is also equally straightforward to write formulae that represent this process in a concise way.

We start with the charge density of equation (2), which should of course extend over the whole simulation cell, however its contributions $\rho_{\alpha \beta}(\mathbf{r})$ from pairs of NGWFs need not, and do not, when they are calculated with the FFT box technique. Therefore the charge density is calculated by replacing these pair contributions in (10) by the following expression:

$$\rho_{\alpha \beta}^{\mathrm{box}}(\mathbf{r})= \hat{Q}^{\da... ...(\alpha \beta)\phi_{\beta}(\mathbf{r}) \right) \, \, \right],$$ (26)

which involves multiplying $\phi_{\alpha}(\mathbf{r})$ with $\phi_{\beta}(\mathbf{r})$ after they have been transferred into the FFT box and as a consequence their product is limited to extend only over the volume of the FFT box. Here we have made use of the fine grid delta function projection operator $\hat{Q}(\alpha \beta)$ which defines a mapping between fine grid delta functions $B_{XYZ}(\mathbf{r})$ of the simulation cell and the fine grid delta functions $b_{xyz}(\mathbf{r})$ of the FFT box. It is defined in an analogous manner to $\hat{P}(\alpha \beta)$ of equation (25).

The matrix elements of one-electron operators, such as the kinetic energy and non-local potential can be easily calculated in the FFT box rather than in the simulation cell by simply transferring the NGWFs in the FFT box before evaluating the integrals. The notation for this set of operations simply involves `sandwiching' the one-electon operator between the standard grid projector operator and its adjoint. The kinetic energy expression of equation (12) becomes

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

and a similar expression can be written for the non-local potential energy.

The matrix elements of the Hartree and local potentials can be calculated with the same ease provided however we have a way to transfer these potentials to the FFT box. This is indeed possible as these local potentials are expressed in terms of the fine grid delta functions and therefore we can transfer them if we use the fine grid delta function projector $\hat{Q}(\alpha \beta)$. Therefore the Hartree energy of equation (14) becomes

$$E_{H}[n]=K^{\alpha \beta} \langle \phi_{\beta} \vert \hat{P}... ...t{P}^{\dag }(\alpha \beta) \vert \phi_{\alpha} \rangle \, \, .$$ (28)

Nothing changes in the evaluation of the exchange-correlation energy which is simply an integral over the fine grid of the whole simulation cell except of course for the fact that the charge density is now calculated by summing $\rho_{\alpha \beta}^{\mathrm{box}}(\mathbf{r})$ terms in equation (10) in place of the $\rho_{\alpha \beta}(\mathbf{r})$.

As far as the optimisation of the energy with respect to the density kernel is concerned, the results of section 3 are still valid provided the Hamiltonian matrix elements that constitute the density kernel gradient (23) are calculated with the FFT box method as we have described above.

In the same way, the calculation of the total energy gradient with respect to the NGWF expansion coefficients of equation (24) can be calculated in the FFT box provided we use only the part of the Kohn-Sham Hamiltonian that exists in the FFT box. Therefore all we have to do is to substitute $\hat{H}$ of equation (24) by

$$\hat{H}(\alpha \beta)= \hat{P}^{\dag }(\alpha \beta) \left[... ...) V_{Hlxc} \right) (\mathbf{r}) \right] \hat{P}(\alpha \beta)$$ (29)

where $V_{Hlxc}(\mathbf{r})$ is the sum of the Hartree potential, local pseudopotential and exchange-correlation potential and of course now this Kohn-Sham operator is in general different for each pair of functions since it contains only the parts of the local potentials that fall inside the FFT box and therefore changes with the location of the FFT box relative to the simulation cell.

As far as implementation is concerned, there are many more issues and algorithmic details about the use of the FFT box that are beyond the scope of this paper. We will describe these in another paper [47].