1. Introduction

Traditional total-energy calculations using density-functional theory (DFT) require a computational effort and quantity of memory which scale as the cube and square of the system size $N$ (i.e. the number of atoms or the volume of the system) respectively. Therefore as systems of increasing size are considered, the computational resources are rapidly exhausted, and a ten-fold increase in computing power will roughly only double the size of system which can be studied. However, the complexity of the problem within DFT scales only linearly with $N$, and there has therefore been considerable interest in developing new schemes for performing these calculations in which the computational effort and memory required also scale linearly: so-called $O(N)$ methods.

One elegant and popular choice of basis in $O(N^3)$ calculations has been the plane-wave basis. However, because of the extended nature of these basis functions they cannot be used in $O(N)$ calculations, and a different choice has to be made, in which the basis functions are localised in real space. Examples include truncated Gaussian orbitals, orbitals based on pseudoatomic wave functions and representing the functions on a real-space grid. An $O(N)$ method results from the combination of a localised basis set and exploitation of the fact that local properties of a system (e.g. the density $n({\bf r})$) depend only upon the electronic states in the vicinity of the point of interest [1].

In this paper we present a set of localised functions which are related to the plane-wave basis set and share some of its attractive features. A significant problem associated with localised basis functions is that they are not in general orthogonal, so that as the size of the basis is increased, the overlap matrix becomes singular. We demonstrate that the basis functions introduced here are orthogonal, by construction, to others centred on the same site, and that the overlap matrix elements for functions centred on different sites can be calculated analytically, and hence evaluated efficiently and accurately when implemented computationally.

Another disadvantage of using basis functions localised in real-space arises in the calculation of the action of the kinetic energy operator. To take advantage of the localisation it is necessary to focus on real-space and calculate all quantities in that representation. However, since the kinetic energy operator is diagonal in reciprocal-space, the kinetic energy matrix elements are most naturally calculated in reciprocal-space. Methods to evaluate the kinetic energy using finite-difference schemes can be inaccurate. With this new choice of basis, the matrix-elements of the kinetic energy operator between any two functions can also be calculated analytically, thus overcoming this problem.

One final advantage arises in the inclusion of non-local pseudopotentials which traditionally required significant computational effort. We present a method of obtaining the matrix-elements of the non-local pseudopotential operator by performing the projection of the basis function onto a core angular-momentum state analytically.