Introduction

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 $\{\phi_\alpha\}$, the single-particle density matrix is expressed as

$$\rho(\mathbf{r}',\mathbf{r})=\sum_{\alpha,\beta}\phi_\alpha(\mathbf{r}') K^{\alpha \beta} \phi^{*}_\beta(\mathbf{r}) \ ,$$ (1)

where $K^{\alpha\beta}$ is the matrix representation of the density matrix in terms of the duals of the support functions. In general the support function set is not orthonormal.

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

$$E_T[\rho]=\sum_{\alpha, \beta} T_{\alpha \beta } K^{\beta \alpha} \ ,$$ (2)

where $T_{\alpha \beta}$ denotes kinetic energy matrix elements between support functions, given in Hartree atomic units by

$$T_{\alpha \beta}= -\frac{1}{2} \int \phi^*_{\alpha}(\mathbf{r}) \nabla^2 \phi_{\beta}(\mathbf{r}) \mathrm{d} \mathbf{r} \ .$$ (3)

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.