Chris-Kriton Skylaris, Oswaldo Diéguez, Peter D. Haynes and Mike C. Payne
Theory of Condensed Matter, Cavendish Laboratory, Madingley Road, Cambridge CB3 0HE, UK
PACS: 31.15.Fx, 71.15.Nc, 71.15.Ap
Electronic structure methods usually require the
solution of Schrödinger's equation
In recent years, the ability to perform calculations using only quantities that are localised in real space, such as the density matrix or Wannier-type orbitals, has led to linear-scaling electronic structure methods that can deal with thousands of atoms . As a consequence, the importance of localised basis sets has grown because they are required for the efficient representation of localised functions. Delocalised basis sets such as plane-waves are not suitable as they make the cost of the electronic structure calculations scale as the cube of the system size.
Basis set expansion is by no means the only
practical approximation in electronic structure theory.
A set of important alternatives are the ``real-space''
methods [12,13] which
use a real-space grid to express their solutions.
The representation of functions directly on a real-space
grid, either regular  or adapted to
the positions of the atoms , simplifies the
application of localisation constraints
and for this reason the importance of
real-space methods in recent years has grown in
parallel with that of the localised basis set methods.
Real-space methods based on a regular grid bear some
similarity to the plane-wave pseudopotential approach as the spacing
of the real-space grid defines
a plane-wave kinetic energy cutoff
and is the parameter with respect to which the
solutions are converged.
Most real-space methods are based on the finite
difference (FD) approach, so
they are not variational, contrary
to the plane-wave basis set method.
This is due to the fact that
the Laplacian operator for the kinetic energy in the
Hamiltonian is represented, or rather approximated,
as a FD expansion
Maragakis et al.  correctly recognised that the FD coefficients of equation (2) systematically underestimate the kinetic energy. They have redetermined them so that they always slightly overestimate it, provided one assumes that the solution on the grid is the real-space representation of a plane-wave expansion. In practice, this means that with their new FD coefficients, energies always converge from above as the grid is made finer. So at least one of the desirable features of the variational principle is restored with this amended FD approach.
A completely different approach which attempts to combine the merits of the plane-wave basis set method with the localisation properties of the real-space methods is the FFT box technique [17,18,19]. Here a plane-wave basis is assumed but the functions are kept localised as they are represented on a real-space grid from the outset. A set of plane-waves which is proportional to the number of grid points in the regions of the localised functions and independent of the size of the system is used whenever the functions are represented in reciprocal-space. Functions localised in real-space are smooth and delocalised in reciprocal-space, therefore coarse sampling in reciprocal-space is a very good approximation . There is no modification of the Hamiltonian operator but an intermediate truncation of the basis set is performed, so the variational behaviour cannot be guaranteed, but can be expected as a consequence of the close resemblance in practice of the FFT box method to the full plane-wave approach.
To assess the variational properties of the new methods, we will first use the exactly solvable model of the harmonic oscillator that Maragakis et al. used for their tests.
The fact that none of the FD methods exhibit the quadratic convergence of variational behaviour may not be a serious drawback for practical calculations since computational limitations always preclude high convergence with respect to grid spacing. However, the speed of convergence is important as it will determine the speed at which relative energy differences converge. To examine this aspect we need to consider an example on a real system. For this purpose, we have performed LDA density-functional calculations on a silane molecule inside a large cubic simulation cell of (40 a). The atomic cores are represented by norm-conserving pseudopotentials  and the charge density is expressed in terms of pseudo-atomic orbitals  centred on the atoms and optimised for each pseudopotential in spherical regions of 8.0 a radius. We compare the FD methods (both the conventional one and Maragakis et al. variant) of order and 28 with the FFT box technique, and with the conventional plane-wave basis set approach. In their finite difference method, Chelikowsky et al.  found that order 12 sufficed for their calculations. We have chosen as an example of very high order of discretization. In all FD calculations we solve the Poisson equation for the Hartree potential using the conventional FD coefficients, as Maragakis et al. do, so the only difference between the two forms of FD is in the calculation of the kinetic energy. Figure 2 shows the convergence of the total energy as a function of the number of grid points in every lattice vector direction.
The computational cost of both FD methods is the same: for an order FD applied to a function localised within a sphere grid points in diameter, the number of operations to calculate the Laplacian at all grid points where it is non-zero  is roughly . For the same function, the FFT box has volume and this method therefore requires about operations (assuming radix-2 FFTs for real functions are used). While the scaling of both methods is essentially the same, the FD method has a smaller prefactor if the discretization order is not very high. As an example, consider spheres of radius and a kinetic energy cutoff of corresponding to a grid spacing of , so that . For , whereas , a factor of 4 larger. The FFT box method is more competitive for smaller spheres, but it should of course be remembered that a higher order FD would need to be used in practice to achieve the accuracy of the FFT box method. For example, if we take in the mentioned test case , then , which is a 50% more than the number of operations required using the FFT box method. Both the FD and FFT box methods are suitable for implementation on parallel computers. For the FFT box method, the spherical regions can be partitioned among the nodes so that the FFTs involved in applying the kinetic energy operator can all be performed locally. Some communication would then be required in order to calculate the inner product of this result with functions in other regions, but this would also be necessary in the case of the FD method.
In summary, we have carried out a comparison of methods for electronic structure calculations directly in real-space in order to examine if and to what extent they can display the variational behaviour of basis sets. The FD representation of the Laplacian for the kinetic energy operator systematically underestimates the kinetic energy and is the source of the problems as correctly identified by Maragakis et al. Their proposed modification of FD coefficients systematically overestimates the kinetic energy and causes the energy to converge from above albeit with errors up to an order of magnitude larger than conventional FD. In practice this means that finer grids may be necessary and this could be a high price to pay in terms of computational cost. None of the FD methods display quadratic convergence of the energy with respect to the wavefunction which is the second characteristic of the variational principle. For cases where the functions are strictly localised in regions of the real-space grid, as with most real-space methods intended for linear-scaling calculations [21,7], the FFT box technique could be the method of choice as it is designed to follow closely the behaviour of the fully variational plane-wave method.
C.K.S. would like to thank the EPSRC (grant number GR/M75525) for postdoctoral research funding. O.D. would like to thank the European Commission for a Research Training Network fellowship (grant number HPRN-CT-2000-00154). P.D.H. would like to thank Magdalene College, Cambridge, for a research fellowship.