Density functional theory (DFT) combined with the pseudopotential method has been established as an important theoretical tool for studying a wide range of problems in condensed matter physics [1]. However, the computational cost of performing a total-energy calculation on a system scales asymptotically as the cube of the system size. Consequently, plane-wave pseudopotential DFT can only be used to study systems of up to about one hundred atoms on a single workstation and up to a few hundred atoms on parallel supercomputers. As a result there has been considerable recent effort in the development of methods whose computational cost scales linearly with system size [2].

A common feature of many of the linear-scaling strategies is the
expansion of the single-particle density matrix in terms of a set
of localised functions. We refer to these functions
as `support functions' [3]. A support function is required
to be non-zero only within a spherical region, which we refer to
as a `support region', centred on an atomic position. Here we consider
a representation of the support functions in terms of a regular
real space grid, which constitutes our basis set. If the set
of support functions is
, the single-particle
density matrix is expressed as

(1) |

Real space methods have the advantage that they provide a clear spatial
partitioning of all quantities encountered in a density functional
calculation, a property that is ideal for code parallelisation. As
a result, this approach has gained popularity in recent years
and a number of such density functional calculations
have been reported by different
authors [4,5,6]. These
approaches use finite difference (FD) methods [7]
for the calculation of the kinetic energy. In terms of the
support functions the kinetic energy is

The evaluation of the kinetic energy matrix elements requires the action of the Laplacian operator on the support functions. Here we will show that in the case of localised support functions, fast Fourier transform (FFT) methods can be adapted for the application of the Laplacian, providing an algorithm with essentially the same computational cost as FD but with higher accuracy and also ready applicability to any grid symmetry.

In the following two sections we present the FD method and our new FFT-based method and compare them both in theory and in practice.