Introduction

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 [1]. Kohn-Sham DFT maps the interacting system of electrons to a fictitious system of non-interacting particles [2], which can be fully described by the single-particle density matrix $\rho(\mathbf{r},\mathbf{r}')$, expressed as a sum of contributions from single-particle Bloch states $\psi_{n \mathbf{k}}(\mathbf{r})$:

$$\rho(\mathbf{r}, \mathbf{r}')= \sum_n f_n \frac{V}{(2\pi)^3... ...thbf{r}) \psi_{n \mathbf{k}}^*(\mathbf{r}') d\mathbf{k} \,\,.$$ (1)

We have assumed that we are dealing with an insulator with completely filled (occupation number $f_n=1$) or empty ($f_n=0$) states and $V$ is the volume of the simulation cell. The $\mathbf{k}$-point integration is carried out in the first Brillouin zone (1BZ). The single-particle states are the eigenfunctions of the Kohn-Sham Hamiltonian at each $\mathbf{k}$-point. They are required to be orthonormal and, in general, extend over the whole simulation cell. The consequence of this orthogonality requirement is that the cost of a DFT calculation involving the $\{ \psi_{n \mathbf{k}} \}$ grows cubically with the system-size. The electronic charge density is equal to the diagonal part of the density matrix multipied by a factor of 2 to take into account the spin degeneracy and is commonly abbreviated as $n(\mathbf{r})= 2\rho(\mathbf{r},\mathbf{r})$.

The most general representation of the density matrix, equivalent to (1), and first applied to linear-scaling DFT calculations by Hernández and Gillan [3], is in terms of a set of localised non-orthogonal functions $\{ \phi_{\alpha \mathbf{R} } \}$:

$$\rho(\mathbf{r}, \mathbf{r}')= \sum_{\alpha \beta} \sum_{\ma... ...r}) K^{\alpha \beta} \phi_{\beta \mathbf{R}}^* (\mathbf{r}'),$$ (2)

where the sum over $\mathbf{R}$ runs over the lattice vectors of the crystal and the matrix $K^{\alpha\beta}$ is called the density kernel, a generalisation of the occupation numbers $\{ f_n \}$. We will call the $\{ \phi_{\alpha \mathbf{R} } \}$ `Non-orthogonal Generalised Wannier Functions' (NGWFs) as they can be derived from a subspace rotation $\mathbf{M}$ between a set of Bloch orbitals at each $\mathbf{k}$-point and a unitary transformation of the results in $\mathbf{k}$-space (Wannier transformation):

$$\phi_{\alpha \mathbf{R} } (\mathbf{r}) = \frac{V}{(2\pi)^3}... ...\mathbf{k}} (\mathbf{r}) M_{n \alpha} \right ] d \mathbf{k} ,$$ (3)

where the number of NGWFs at each lattice cell $\mathbf{R}$ can be greater than or equal to the number of occupied Bloch bands at each $\mathbf{k}$-point. The density kernel is the result of applying the inverse of these transformations on each side on the diagonal occupation number matrix diag$( \{ f_n \} )$:

$$K^{\alpha \beta}= \sum_n N^{\alpha}_{\hspace{1ex} n} \, f_n \, (N^{\dag })_n^{\hspace{1ex} \beta},$$ (4)

where $\mathbf{N}=\mathbf{M}^{-1}$.

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 [4] 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 $\phi_{\alpha \mathbf{R}}(\mathbf{r})$ and $\phi_{\beta \mathbf{R}}(\mathbf{r})$ 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 [11]. Our approach is closely related to that of Hernández et al. [12] who developed a method to do calculations by optimising both the density kernel and the functions $\{ \phi_{\alpha \mathbf{R} } \}$, 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. [27]; 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.