Tests and discussion

We have performed tests of the FD and FFT box methods for calculating the kinetic energy of localised functions. Choosing a particular type of support function $\phi$ with spherical symmetry, placing one at $\mathbf{R}_\alpha$ and another at $\mathbf{R}_\beta$, we rewrite the integral of equation (3) as

$$T(\vert\mathbf{R}_\alpha-\mathbf{R}_\beta\vert) = -\frac{1}{... ...^2 \phi(\mathbf{r}-\mathbf{R}_\beta) \mathrm{d} \mathbf{r} \ .$$ (6)

For our first test we calculate the following quantity as a function of the distance $d$ between the centres $\mathbf{R}_\alpha$ and $\mathbf{R}_\beta$

$$\eta_{1}(d)=T_{ap}(d)-T_{ex}(d) \ ,$$ (7)

where $T_{ex}(d)$ is the exact value of the integral in the continuous representation of the support functions and $T_{ap}(d)$ is its approximation on the real space grid, either by FD or the FFT box method. We chose $\phi(\mathbf{r})$ to be a 2s valence pseudo-orbital for a carbon atom, generated using an atomic norm-conserving carbon pseudopotential [13] within the local density approximation. The pseudo-orbital is confined in a spherical region of radius $6.0 \: a_{0}$, and vanishes exactly at the region boundary [14]. It is initially generated as a linear combination of spherical Bessel functions, which are the energy eigenfunctions of a free electron inside a spherical box. Our functions are limited up to an energy of 800 eV, resulting in a combination of fourteen Bessel functions. The formula for calculating kinetic energy integrals between Bessel functions is known [15] and we used it to obtain $T_{ex}(d)$ for our valence pseudo-orbital. We then calculated $\eta _{1}(d)$ with a grid spacing of $0.4 \: a_0$ (corresponding to a plane-wave cut-off of 839 eV) in an orthorhombic simulation cell, as we are restricted to do so by the FD method. With these parameters $N_{box}$ is $60^{3}$, and hence it is trivial to perform the FFT of one support function on a single node. $\eta _{1}(d)$ is plotted for the FFT box method and for various orders of the FD method in the top graph of figure 2.

It can be seen that low order FD methods are inaccurate as compared to the FFT box method, and only when order 28 FD is used does the accuracy approach that of the FFT box method. The A = 12 FD scheme, the highest order that has been used in practice for calculations [9], gives an error of $-3.97 \times 10^{-5} \: \mbox{Hartree}$ at $d=0$ as compared to $1.027 \times 10^{-5} \: \mbox{Hartree}$ for the FFT box method. The feature that occurs in the top graph of figure 2, between $d=5 \: a_0$ and $d=7 \: a_0$, is an artefact of the behaviour of our pseudo-orbitals at the support region boundaries where they vanish exactly, but with a finite first derivative. This causes an enhanced error in all the methods when the edge of one support function falls on the centre of another.

The error in the FFT box method is small, yet non-zero, and we attribute this to the inherent discretisation error associated with representing functions that are not bandwidth limited on a discrete real space grid. Convergence to the exact result is observed as the grid spacing is reduced, as expected.

Figure 2: Top panel: $\eta _{1}(d)$ for a carbon 2s valence pseudo-orbital in a spherical support region with a radius of $6.0 \: a_{0}$. Bottom panel: $\eta _{2}(d)$ for the same pseudo-orbital. The insets show a magnification of the plots for A = 28 FD and the FFT box method.

As our next comparison of the FFT box and FD methods we used the same pseudo-orbitals as before, but considered the quantity

$$\eta_{2}(d)=T_{ap}(d)-T_{PW}(d)$$ (8)

as the measure of the error, where $T_{PW}$ is the result obtained by Fourier transforming the support functions using the periodicity of the entire simulation cell. One may think of $T_{PW}$ as being the result that would be obtained from a plane-wave code: the support functions may be considered as generalised Wannier functions. Calculating the kinetic energy integrals by performing a discrete Fourier transform on the support functions over the entire simulation cell (an $O(N_{s}^{2}\log N_{s})$ process for all support functions) is equivalent to summing the contributions to the kinetic energy from all of the plane-waves up to the cut-off energy determined by the grid spacing. Thus our FFT box method can be viewed as equivalent to a plane-wave method that uses a contracted basis set (i.e. a coarse sampling in reciprocal space). In some ways $\eta _{2}(d)$ is a better measure of the relative accuracy of the FD and FFT box methods as our goal is to converge to the `exact' result as would be obtained using a plane-wave basis set over the entire simulation cell. $\eta _{2}(d)$ is plotted in the bottom graph of figure 2.

$T_{PW}$ was computed using a cell that contained 256 grid points in each dimension. Increasing the cell size further had no effect on $T_{PW}$ up to the eleventh decimal place ( $10^{-11} \: \mbox{Hartree}$). The plots show that the FFT box method performs significantly better than all orders of FD that were tested. For example at $d=0$ the error for A = 28 FD is $-3.49 \times 10^{-6} \: \mbox{Hartree}$ as compared to $-1.09 \times 10^{-9} \: \mbox{Hartree}$ for the FFT box method. The fact that the FFT box error is so small shows that coarse sampling in reciprocal space has little effect on accuracy, as one would expect for functions localised in real space.

Our implementation can produce similar FFT box results to the above in regular grids of arbitrary symmetry (non-orthogonal lattice vectors) as long as we include roughly the same number of grid points in the support region sphere. As we described earlier the application of the FD method to grids without orthorhombic symmetry is not straightforward.

Furthermore, in our implementation the kinetic energy matrix elements $T_{\alpha \beta}$ for both the FFT box method and the FD method (of any order) are Hermitian to machine precision. This is a direct consequence of $\hat{T}$, our representation of the Laplacian operator $\nabla^{2}$ on the grid, being Hermitian. As mentioned earlier, this is an important point. The matrix elements $T_{\alpha \beta}$ may always be made Hermitian by construction without $\hat{T}$ itself being an Hermitian operator. This would ensure real eigenvalues, as is required. However, when it comes to optimisation of the support functions during a total-energy calculation, we require the derivative of the kinetic energy with respect to the support function values [12]:

$$\frac{\partial E_T}{\partial \phi^{*}_{\alpha} (\mathbf{r}_i... ...ta K^{\alpha\beta}\hat{T} \phi^{*}_{\beta} (\mathbf{r}_i ) \ ,$$ (9)

where the $\mathbf{r}_i$ are grid points belonging to the support region of $\phi_\alpha$. These relations both hold only if $\hat{T}$ is an Hermitian operator, and support function optimisation can only be performed in a consistent manner if there is one unique representation of $\hat{T} \phi_{\beta}$ for each support function $\phi_{\beta}$. It is also worth noting that the evaluation of these derivatives is the reason why we prefer to perform the sum of equation (4) for the FFT box method in real space, rather than in its equivalent form in reciprocal space. Applying the FFT box method in reciprocal space would be no more costly as far as integral evaluation is concerned but we would require an extra FFT per support function for the subsequent evaluation of equation (9).

For all the methods we describe in this paper we observe variation in the values of the kinetic energy integrals when we translate the system of the two support functions with respect to the real space grid. This is to be expected as the discrete representation of the support functions changes with the position of the support region with respect to the grid. Such variations may have undesirable consequences when it comes to calculating the forces on the atoms. In FFT terminology, they result from irregular aliasing of the high frequency components of our support functions as they are translated in real space. Ideally, in order to avoid this effect, the reciprocal representation of the support functions should contain frequency components only up to the maximum frequency that corresponds to our grid spacing, in other words it should be strictly localised in reciprocal space. Unfortunately this constraint is not simultaneously compatible with strict real space localisation. It should be possible however to achieve a compromise, thus controlling the translation error by making it smaller than some threshold. Such a compromise should involve an increase in the support region radii of our functions by a small factor. This situation is similar to the calculation of the integrals of the nonlocal projectors of pseudopotentials in real space with the method of King-Smith et al. [16] which requires an increase of the core radii by a factor of 1.5 to 2. For example, if we consider two carbon valence pseudo-orbitals of support radius $6.0 \: a_0$ and with $d=5.0 \: a_0$ and translate them both in a certain lattice vector direction over a full grid spacing, the maximum variation in the value of the integral with the FFT box method is $8.28 \times 10^{-6}$ Hartree. If we then do the same with carbon pseudo-orbitals generated with precisely the same parameters but instead with a support radius of $10.0 \: a_0$, the maximum variation with respect to translation is reduced to $2.05 \times 10^{-8}$ Hartree.