Introduction

Density functional theory (DFT) used in conjunction with the plane-wave pseudopotential method has been widely used as a theoretical tool for studying a variety of condensed matter systems [1]. The original formulation of Kohn and Sham [2] in terms of extended orthonormal wavefunctions has a computational cost that scales asymptotically as the cube of the system size. This bottleneck restricts the approach to the study of no more than a few hundred atoms, even on parallel supercomputers. As a consequence, over the last decade there has been much effort devoted to the development of methods whose computational cost scales linearly with system size [3]. Such schemes often rely on a reformulation of the problem in terms of localised functions in real space [4]. One such set of functions are the orthogonal Wannier functions: a unitary transformation of the extended Bloch wavefunctions [5]. These (and as a result the density matrix) are known to be exponentially localised in insulators, with the degree of localisation determined by the band gap [6,7,8,9]. Furthermore, it is known that an essentially equivalent representation that is non-orthogonal can be better localised [10] and hereafter we shall refer to such localised functions as non-orthogonal generalised Wannier functions (NGWFs).

In our approach we formulate the total energy of the system in terms of NGWFs that are represented on a real space grid subject to periodic boundary conditions. The localisation properties of the NGWFs are set a priori, i.e., each NGWF (usually atom centered) is confined within a spherical localisation region (LR) whose radius is pre-defined. In general, each NGWF may have a different localisation radius. Our basis set consists of a mesh of cell periodic functions which are a unitary transformation of an equivalent plane-wave basis. As a result, we can use fast Fourier transform (FFT) methods familiar from plane-wave DFT calculations to switch easily between real and reciprocal space representations.

Real space techniques have been used by many authors for performing DFT calculations. In particular, there are methods which use functions that are strictly localised within spherical regions on real space grids [11,12]. These approaches have given rise to some important methodological developments. However, it would be desirable to have a method with a basis set that can be improved systematically [13,14], just as in plane-wave DFT. Our approach uses a basis that is directly comparable to a plane-wave basis.

We begin by describing our basis set and its properties in section 2. In section 3 we describe how to calculate the total energy of a system with a set of localised functions and in section 4 we introduce our novel FFT box technique and demonstrate how it can be used to calculate the total energy with a computational cost that scales linearly with system size. Section 5 describes how we minimise the total energy, and we present our results and conclusions in sections 6 and 7 respectively.