The pseudopotential plane-wave method for density-functional
theory (DFT) calculations has been developed and perfected over many years
into a reliable tool for predicting static and dynamic
properties of molecules and solids .
Kohn-Sham DFT maps the interacting system of electrons to
a fictitious system of non-interacting particles ,
which can be fully described by the single-particle density matrix
, expressed as a sum of contributions
The most general representation of the density matrix, equivalent to
(1), and first applied to linear-scaling DFT
calculations by Hernández and Gillan , is in terms
of a set of localised non-orthogonal functions
In our presentation so far, we have considered the Bloch states as the natural representation and starting point from which to construct the charge density and NGWFs. This is also the usual order which has been followed in discussions and derivations of Wannier functions  in the literature. Generalised Wannier functions, orthonormal or not, are most commonly constructed in a post-processing fashion [5,6,7,8] after the end of a plane-wave band structure calculation.
Since the energy is variational with respect to the charge density, directly varying the NGWFs and the density kernel of equation (2) is equivalent to varying the Bloch states of equation (1). The localisation properties of the NGWFs are set a priori, i.e. each NGWF is non-zero only within a pre-defined localisation region. As a result, the computational cost of a density-functional calculation scales only quadratically with system-size rather than cubically. Furthermore, it can be made to scale linearly with one extra variational approximation: the truncation of the density kernel in equation (2) when the centres of the functions and lie beyond some cutoff distance [9,10].
The NGWFs are expanded in terms of a basis of periodic, bandwidth limited delta functions (Appendix A). These are centered on the points of a regular real-space grid and are related to an equivalent plane-wave basis through a unitary transformation (the Fourier transform). Hence, the method we present is directly equivalent to the Bloch state pseudopotential plane-wave density-functional approach. The delta function basis allows us to restrict our NGWF expansions to contain only the delta functions that are included in some spherical localisation region. This should be an accurate assumption for insulators in which the NGWFs decay exponentially . Our approach is closely related to that of Hernández et al.  who developed a method to do calculations by optimising both the density kernel and the functions , which they call `support functions'.
Real-space methods, as an alternative to pure reciprocal-space plane-wave methods, have been used by many other authors for DFT calculations [13,14,15,16,17,18,19,20] in the past in order to take advantage of the benefits of localisation in real-space. In particular, approaches have been developed that use functions strictly localised in spherical regions on real-space grids [21,22]. These functions play the same role as the NGWFs we present here. A different approach is taken by the SIESTA program[23,24], which uses a basis set of numerical atomic orbitals. These are generated as described by Artacho et al. [25,26] and are not optimised during the calculations.
Even though these real-space methods have led to some important methodological developments, it would be very desirable if they could be directly comparable with plane-wave pseudopotential DFT. In other words, we wish to have a method that rigorously adheres to a basis set we can improve systematically, such as the plane-wave basis where its convergence towards completeness is controlled by the kinetic energy cutoff parameter. Our method achieves this by working both in real- and reciprocal-space. We demonstrate that our approach can actually be viewed as an alternative way of performing plane-wave DFT calculations which is easy to turn into a linear-scaling method in the future with only trivial modifications. It is thus also directly applicable to any Bravais lattice symmetry, in contrast to common finite difference methods that are usually restricted to orthorhombic lattices.
We should note at this point that there exist other basis sets, apart from plane-waves, that can be improved systematically: the B-splines of Hernández et al. ; and the polynomial basis in the finite-element approach of Pask et al. [28,29], who have done some pioneering work using a technique for solving differential equations common in engineering applications and adapting it for electronic structure calculations.
In what follows, we begin by describing the calculation of the total energy in our scheme, directly with NGWFs, in section 2. In section 3 we describe our strategy for total energy optimisation, i.e. minimisation of the energy with respect to both the density kernel and the NGWFs. In section 4 we present the FFT box technique, an essential ingredient for lowering the cost of the calculations and for eventually achieving linear-scaling behaviour. In section 5 we present tests on a variety of systems showing the accuracy and efficiency of this method and finally we conclude and mention what we see as future developments.