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 with spherical symmetry, placing
one at
and another at
, we
rewrite the integral of equation (3) as
(6) 
(7) 
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 at as compared to for the FFT box method. The feature that occurs in the top graph of figure 2, between and , is an artefact of the behaviour of our pseudoorbitals 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 nonzero, 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.

As our next comparison of the FFT box and FD methods we
used the same pseudoorbitals as before, but considered the quantity
(8) 
was computed using a cell that contained 256 grid points in each dimension. Increasing the cell size further had no effect on up to the eleventh decimal place ( ). The plots show that the FFT box method performs significantly better than all orders of FD that were tested. For example at the error for A = 28 FD is as compared to 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 (nonorthogonal 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
for both the FFT box
method and the FD method (of any order) are Hermitian
to machine precision. This is a direct consequence
of , our representation of the Laplacian
operator on the grid, being Hermitian. As
mentioned earlier, this is an important point. The
matrix elements
may always
be made Hermitian by construction without 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 totalenergy calculation, we require
the derivative of the kinetic energy with respect
to the support function values [12]:
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 KingSmith 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 pseudoorbitals of support radius and with 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 Hartree. If we then do the same with carbon pseudoorbitals generated with precisely the same parameters but instead with a support radius of , the maximum variation with respect to translation is reduced to Hartree.