# Practical Methods for Ab Initio Calculations on Thousands of Atoms

## I Introduction

Over the last thirty years, the use of ab initio electronic structure techniques has become widespread in chemistry and physics. However, all traditional techniques are limited in the system size which they can treat (often defined by the number of atoms, $`N`$) by poor scaling of the computational effort required, which is generally at least $`𝒪(N^2)`$, if not worse. In the context of standard self-consistent field (SCF) quantum chemical methods such as Hartree-Fock theory or density functional theory (DFT), there are two demanding parts of a calculation: first, the build of the Fock matrix (or, in DFT, the Hamiltonian matrix), which can scale as $`𝒪(N^2)`$; second, the solution for eigenvectors of the Fock matrix, which scales as $`𝒪(N^3)`$ if performed as a diagonalisation; this only dominates for large values of $`N`$. An alternative to diagonalisation is iterative minimisation, which has $`𝒪(N^2)`$ scaling (one $`N`$ dependence comes from the eigenvectors spreading over all space, and the other from the number of eigenvectors, which depends on $`N`$); if, as is often the case, the eigenvectors must be orthogonalised to each other, this leads asymptotically to $`𝒪(N^3)`$ scaling. Whatever the cause of poor scaling, it results in a practical limit of a few hundred atoms in conventional ab initio techniques. The desire to model large systems, e.g. biomolecules or nanostructures, has seen a strong push in recent years to achieve linear scaling with system size.

The building of the Fock matrix (and specifically the Coulomb and exchange terms which are the most expensive) with linear scaling has been addressed recently by several groups; however, this aspect will not concern us here. Rather, we will consider linear scaling techniques (which also rely on iterative minimisation) for finding the self-consistent ground state of the system.

Recently, many linear scaling techniques have been proposed, which are all based on the search for the density matrix. They start from the observation that the density matrix between two points in space decays in some manner as those two points increase in separation. This is intuitively clear, as it is well-known that bonding is local. The result is that the electronic structure of an atom depends only on its local environment, so that if the overall size of the system changes, there is no effect on the local electronic structure; thus the effort required to solve for the whole system should be proportional to $`N`$. This is the foundation of all linear scaling electronic structure methods (with a few exceptions which will not concern us here).

The paper is divided up as follows: Section II describes the basic theory behind our $`𝒪(N)`$ DFT method, while Section III presents recent advances we have made in searching for the electronic ground state. The future directions the work will take are shown in Section IV, and the paper is concluded in Section V.

## II $`𝒪(N)`$ density functional theory

The density matrix within density functional theory (DFT) can be written as:

$$\rho (𝐫,𝐫^{})=\underset{i}{}f_i\psi _i(𝐫)\psi _i^{}(𝐫^{}),$$
(1)   
where $`\psi _i`$ is a Kohn-Sham eigenfunction, and $`f_i`$ is the occupancy of that eigenfunction. The key observation which underpins $`𝒪(N)`$ DFT is that DFT can be formulated in terms of $`\rho (𝐫,𝐫^{})`$, and that the ground state can be found by minimising the total energy, $`E_{\mathrm{tot}}`$, with respect to $`\rho (𝐫,𝐫^{})`$ subject to the condition that $`\rho (𝐫,𝐫^{})`$ is idempotent (these statements are proved elsewhere). Idempotency means that $`\rho \rho =\rho `$, which is equivalent to the eigenvalues of $`\rho `$ (which are the occupation numbers $`f_i`$) being either zero or one. This is a crucial property for the density matrix, as it is ensures that the density matrix is a projector – it is the operator which projects onto the space of occupied states.

Another important property of the density matrix is that it decays as the separation between points increases:

$$\rho (𝐫,𝐫^{})0\mathrm{as}𝐫𝐫^{}\mathrm{}.$$
(2)   
The fundamental reason for this decay is the loss of quantum phase coherence between distant points. The net result is that the local environment is all that is important in determining $`\rho (𝐫,𝐫^{})`$, and thus that the amount of information contained in $`\rho `$ is proportional to $`N`$. This decay property of the density matrix can be used to make the amount of information in the system strictly linear with the system size by imposing a constraint, and setting the density matrix to zero beyond a cutoff:

$$\rho (𝐫,𝐫^{})=0,𝐫𝐫^{}>R_c,$$
(3)   
where $`R_c`$ is some cut-off radius. Solving for the energy with this constraint imposed, and that of idempotency, will lead to an upper bound on the ground state energy<sup>*</sup><sup>*</sup>*The idempotency constraint on $`\rho `$ will make it the projection operator onto the occupied states. The energy of the ground state given by these states will be the energy for the system under the constraint of a localised $`\rho `$. Since DFT is variational, this extra constraint will raise the energy above the true ground state energy (i.e. without localisation), and so will give an upper bound to the true ground state energy.; as $`R_c`$ is increased, it will converge to the true ground state. Clearly there is a balance to be struck between accuracy, which increases as $`R_c`$ is increased, and the complexity of the computational problem (e.g. number of variational degrees of freedom associated with $`\rho `$), which also increases with $`R_c`$.

This exact formulation cannot be followed for practical methods, as $`\rho `$ is dependent on two vector positions. Instead, a further approximation is introduced, namely that $`\rho `$ be separable, so that it can be written in the form:

$$\rho (𝐫,𝐫^{})=\underset{i\alpha ,j\beta }{}\varphi _{i\alpha }(𝐫)K_{i\alpha ,j\beta }\varphi _{j\beta }(𝐫^{}),$$
(4)   
which is equivalent to requiring that $`\rho `$ only have a finite number of non-zero eigenvalues. The functions $`\varphi _{i\alpha }(𝐫)`$ are known as ‘localised orbitals’, where $`i`$ runs over atoms and $`\alpha `$ over localised orbitals on each atom. The matrix $`K_{i\alpha ,j\beta }`$ is the density matrix in the basis set of $`\{\varphi _{i\alpha }(𝐫)\}`$ (and is identical to the density matrix commonly seen in $`𝒪(N)`$ tight binding schemes; indeed many of these schemes can be used to solve for this matrix within $`𝒪(N)`$ DFT).

The localisation of $`\rho `$ is accomplished by setting the localised orbitals to be non-zero only inside a certain radius, $`R_{\mathrm{reg}}`$, centred on the atoms $`i`$. A spatial cutoff (not generally of the same value) must also be imposed on the matrix $`K`$, so that $`K_{i\alpha ,j\beta }=0`$ once atoms $`i`$ and $`j`$ are more than a certain distance apart.

### A A Specific Implementation: CONQUEST

The framework described above is completely general; we will now concentrate on our specific implementation, Conquest. This is based on the pseudopotential approach to DFT (described briefly in Appendix A), and has been constructed so as to be as accurate as conventional ab initio pseudopotential calculations, which use plane waves as a basis set.

In practice, there are various issues which must be addressed: minimising the total energy with respect to $`K_{i\alpha ,j\beta }`$ while maintaining idempotency and spatial localisation; representing the localised orbitals; the cutoffs required to achieve good convergence to the true ground state; and practical questions, including integration and implementation on parallel computers. These have been addressed in detail elsewhere; we present a brief summary here.

The imposition of idempotency on $`K`$ during the minimisation of $`E_{\mathrm{tot}}`$ with respect to $`K_{i\alpha ,j\beta }`$ is the hardest constraint to maintain. There are several proposed means of accomplishing this; the method described here is based on the purification technique of McWeeny, recently used in tight binding calculations by Li, Nunes and Vanderbilt, and described in detail in Section III A. It requires $`K`$ to be written in terms of an ‘auxiliary’ density matrix, $`L`$:

$`K=3LSL2LSLSL,`$ (5)   
where $`S`$ is the overlap matrix:

$$S_{i\alpha ,j\beta }=d𝐫\varphi _{i\alpha }(𝐫)\varphi _{j\beta }(𝐫).$$
(6)   
The localisation of $`K`$ is then imposed as a spatial cutoff on $`L`$:

$$L_{i\alpha ,j\beta }=0,𝐑_i𝐑_j>R_\mathrm{L},$$
(7)   
where $`𝐑_i`$ is the position of atom $`i`$ and $`R_\mathrm{L}`$ is a cutoff radius. The energy is then minimised with respect to the matrix elements $`L_{i\alpha ,j\beta }`$ using the standard conjugate gradients technique. The effect of the purification is to make $`K`$ more nearly idempotent given an $`L`$ which is nearly idempotent.

The next issue to consider is the representation of the localised orbitals. The lessons learnt from the use of plane waves in conventional pseudopotential calculations are helpful to remember here. First, they offer a systematic convergence of the energy with respect to the basis set completeness, and this is achieved with a single parameter (the cutoff energy). Second, they are bias free - that is, they are completely flexible, and no knowledge of the kind of bonding in the system is required. If possible, the choice of basis set would reflect these qualities.

Conquest represents the localised orbitals in a real-space basis. There are various possibilities for an efficient, real-space basis. One is to use the spherical equivalent of plane waves, that is spherical Bessel functions $`j_l(𝐫)`$ combined with spherical harmonics $`Y_m^l(𝐫)`$, within each localisation region. This representation has been discussed by Haynes and Payne, but practical results have not yet been reported. Another $`𝒪(N)`$ DFT scheme uses pseudo-atomic orbitals as the basis functions, with considerable success. An alternative is to represent the $`\varphi _{i\alpha }`$ by their values on a grid, and to calculate matrix elements of the kinetic energy by taking finite differences. This technique is well established in conventional first principles calculations, and has been investigated in the context of $`𝒪(N)`$ techniques by us and recently by Hoshi and Fujiwara. At present we use a basis of B-splines, $`\mathrm{\Theta }(𝐫)`$, (also called blip functions) in the expansion:

$$\varphi _{i\alpha }(𝐫)=\underset{s}{}b_{i\alpha s}\mathrm{\Theta }(𝐫𝐑_{is}),$$
(8)   
where the B-splines are piecewise polynomial functions (continuous up to the third derivative) which are strictly localised on the points of a grid (notated as $`𝐑_{is}`$ above) which is rigidly attached to each atom (known as the blip grid). The energy is then minimised with respect to the coefficients of the B-splines, $`b_{i\alpha s}`$.

Having described the minimisation of the total energy with respect to both the $`K`$ matrix and the localised orbitals, it is now appropriate to consider the practical performance of the method: are the cutoffs required to achieve convergence to the ground state small enough to be practical ? Tests on a model, local pseudopotential and standard non-local pseudopotentials have already been reported, and show that for reasonable cutoffs, good convergence is obtained. We reproduce some of these results in Figure 1. Fig. 1(a) shows the calculated total energy as a function of $`R_{\mathrm{reg}}`$ for Si. The results show that $`E_{\mathrm{tot}}`$ converges to the correct value extremely rapidly once $`R_{\mathrm{reg}}`$ is greater than 4 Å. For this radius, each localisation region contains 17 neighbouring atoms, and the calculations are perfectly manageable. Fig. 1(b) shows the total energy for $`R_{\mathrm{reg}}`$ = 2.715 Å, as a function of $`R_L`$. Rather accurate convergence to the $`R_\mathrm{L}=\mathrm{}`$ value is obtained for $`R_\mathrm{L}`$ 8 Å, which again is acceptable. No value is shown for exact diagonalisation because of technical difficulties in performing comparisons.

To perform integrations such as $`S_{i\alpha j\beta }=d𝐫\varphi _{i\alpha }(𝐫)\varphi _{j\beta }(𝐫)`$ numerical integration on a grid is used. This integration grid is generally of different spacing to the blip grid (and normally about twice as fine). Most matrix elements are found by integration on this grid (with the exception of fast-varying quantities which are calculated analytically), and the localised orbitals are projected from the blip grid to the integration grid in a manner similar to a fast Fourier transform (FFT), called a blip-to-grid transform. Once the charge density is calculated on the grid (as $`n(𝐫)=\rho (𝐫,𝐫)`$), the Hartree potential and energy are found using FFTs.

Conquest has been designed with parallel computers in mind; here we summarize the strategy. Each processor has three responsibilities. First, a group of atoms: it holds the blip coefficients, $`b_{i\alpha s}`$ and their derivatives of the energy, $`E_{\mathrm{tot}}/b_{i\alpha s}`$ and is responsible for performing the blip-to-grid transforms for these atoms. It also stores the rows of matrices corresponding to these atoms and performs the matrix multipications for these rows. Second, a domain of integration grid points: it has responsibility for calculating contributions to matrix elements arising from sums over these points, and for holding the electron density and the Kohn-Sham potential on these points. Third, part of the spatial FFT for the Hartree potential: it deals with a set of columns in the $`x,y`$ or $`z`$ directions. All processors switch between tasks in a concerted manner.

To test the efficiency of the scheme, we have extensively tested its scaling properties. There are two completely different kinds of scaling: parallel scaling (i.e. scaling of CPU time for a given system with varying numbers of processors); and inherent scaling (i.e. scaling of CPU time for a fixed number of processors as the system size varies. In the present implementation of Conquest, both types of scaling are excellent.

The overall Conquest scheme can be summarised as follows: the ground state energy and density matrix of the system are found by minimising the energy $`E_{\mathrm{tot}}`$ with respect to the elements of the auxiliary density matrix, $`L_{i\alpha ,j\beta }`$, and the localised orbitals, $`\varphi _{i\alpha }(𝐫)`$, subject to the spatial cutoffs $`R_L`$ and $`R_{\mathrm{reg}}`$. This yields an upper bound to the true ground state, which improves as the cutoffs are increased.

## III Strategies for reaching the ground state

Having described the specific manner in which Conquest is implemented, we now consider ways of reaching the ground state efficiently and robustly. At present, the minimisation is carried out in three separate stages: first, the ground state density matrix is found for a given set of localised orbitals; second, self-consistency is achieved between the charge density and the potential (which includes further density matrix minimisation in response to the new potential); finally, the form of the localised orbitals is changed in accordance with the gradient of the energy. The inner loops (density matrix minimisation and self-consistency) are then repeated. Schemes for efficiency and robustness for these three stages are now discussed.

### A Density matrix minimisation

As has already been emphasised, the density matrix of the true electronic ground state is idempotent. This important property is hard to impose during a minimisation; however, a number of schemes which drive the matrix towards idempotency have been suggested. For simplicity, we will consider these schemes in the framework of orthogonal tight binding theory; the extension to the non-orthogonal case and DFT is simple enough. The first scheme was proposed by McWeeny, who noted that if a matrix $`\rho `$ is close to idempotency, then the matrix $`\stackrel{~}{\rho }`$ given by:

$$\stackrel{~}{\rho }=3\rho ^22\rho ^3$$
(9)   
will be more nearly idempotent. It has the effect of driving the eigenvalues of $`\rho `$ towards zero and one (this can be seen by considering the function $`3\lambda ^22\lambda ^3`$, which is shown in Figure 2). If this transformation (often called the McWeeny transformation or purification transformation) is applied iteratively (writing $`\rho _{n+1}=3\rho _n^22\rho _n^3`$ for iteration $`n+1`$), then the sequence of matrices generated will converge on an idempotent matrix. In fact, this transformation is quadratically convergent (i.e. if the idempotency error in $`\rho `$ is $`\delta \rho `$, then the idempotency error in $`\stackrel{~}{\rho }`$ is $`\delta \rho ^2`$).

Palser and Manolopoulos have recently suggested using this iterative scheme in an $`𝒪(N)`$ manner. They point out that if the initial density matrix is a linear function of the Hamiltonian, with eigenvalues between zero and one, then the iteration will converge to the correct ground state density matrix (given by $`\theta (\mu H)`$, where $`\theta (x)`$ is the Heaviside step function ($`\theta =1`$ for $`x>0`$ and $`\theta =0`$ for $`x<0`$) and $`\mu `$ is the chemical potential for electrons, or the Fermi energy), and that the energy, $`E=2\mathrm{T}\mathrm{r}[\rho _nH]`$, will decrease monotonically at each step. This procedure has the advantage of being fast (it only requires two matrix multiplies per iteration) and efficient (it converges quadratically). Unfortunately, when a localisation criterion (also called a truncation) is applied to the density matrix to achieve linear scaling, the monotonic decrease of energy will fail at some point in the iterative search. This can be taken as an indication that truncation errors are dominating the calculation, and that the search should be stopped; indeed, if it is not stopped, there is no guarantee that it will continue to converge towards an idempotent matrix. This is a heuristic criterion for stopping the iteration, and has the drawback that the method will not be variational, so that analytic forces will not be in agreement with the numerical gradient of the energy.

The Li, Nunes and Vanderbilt (LNV) scheme for achieving the ground state density matrix also uses the McWeeny transformation, though in a rather different manner. Here, the energy is written as $`E=2\mathrm{T}\mathrm{r}[\stackrel{~}{\rho }H]`$, with $`\stackrel{~}{\rho }`$ given by equation 9. Then the energy is minimised with respect to the elements of $`\rho `$, typically using a scheme such as conjugate gradients to generate a sequence of search directions. The localisation of the density matrix is achieved by applying a spatial cutoff to the elements of $`\rho `$. This scheme has at least two advantages: first, each line minimisation can be performed analytically, as the energy is cubic in $`\rho `$; second, it is variational, so that the energy found is always an upper bound to the ground state, and forces obey the Hellmann-Feynman theorem and are in exact agreement with the numerical derivative of the energy.

However, there are drawbacks to the LNV technique. First, it is unclear what the best initial value should be for the density matrix; typically, it is taken to be $`\frac{1}{2}`$I, or $`\frac{1}{2}`$S<sup>-1</sup> in a non-orthogonal basis set. Second, as the McWeeny transformation is a cubic, it is unbounded from below, and a poor starting choice for the minimisation can lead to runaway solutions; a sign of this is typically that the cubic has complex extrema. Third, the scheme can be poorly convergent, and is not guaranteed the quadratic convergence of the McWeeny method.

We have recently proposed a hybrid between the McWeeny and LNV schemes which builds on the complementarity between these two methods. It is based on the observation that the sequence of matrices generated during a McWeeny iterative search converges quadratically towards idempotency, and that the LNV search direction maintains idempotency in the density matrix to first order. Thus the McWeeny scheme is used as an initialisation to find an idempotent density matrix (but one which is not the ground state matrix because of truncation errors); this matrix is then used as the input to the LNV scheme, which maintains the idempotency while searching efficiently for the ground state density matrix. The combination of the two methods is both variational and robust - two highly desirable attributes for the inner loop of an $`𝒪(N)`$ DFT method.

As an example of the improved speed of convergence given by the hybrid scheme, Figure 3 shows the convergence to the ground state energy in diamond carbon for the LNV stage of the hybrid scheme and pure LNV (initialised from $`\rho =\frac{1}{2}𝐈`$). These results show that the McWeeny stage of the hybrid scheme gets closer to the ground state as the radius is increased, as expected, and that it acts as an excellent initial density matrix for the LNV scheme. The method has also been tested on a vacancy in diamond C, the Si(001) surface and liquid Si, as reported elsewhere.

### B Non-orthogonality

The previous section focuses on the traditional tight binding scheme where the localised orbitals are taken to be orthogonal. However, in our $`𝒪(N)`$ method, the orbitals are non-orthogonal, and this introduces a significant degree of complication. In strict mathematical terms, the metric for the space spanned by the localised orbitals must be chosen with care; this means defining a scalar product in a specific way, as explained in detail in Appendix B, along with other implications of using non-orthogonal orbitals. We will consider a few simple implications of this theory here.

For any given set of non-orthogonal orbitals, {$`\overline{\varphi }_i`$}, an orthogonal set can be defined by using the overlap matrix, $`S`$:

$$\varphi _i=\underset{j}{}(S^{1/2})_{ij}\overline{\varphi }_j,$$
(10)   
where we use an over-bar to indicate the quantities in the non-orthogonal case. If the metric is chosen suitably, then similar transformations between the Hamiltonian and density matrices in the two spaces can be defined:

$`\overline{H}`$ $`=`$ $`\overline{S}^{1/2}H\overline{S}^{1/2}`$ (11)   
$`\overline{\rho }`$ $`=`$ $`\overline{S}^{1/2}\rho \overline{S}^{1/2}.`$ (12)   
There are various points which can be drawn from the above equations. First, the McWeeny transformation (and the other quantities associated with a minimisation such as the gradient of the energy) will change in the new basis set; in fact the McWeeny transformation becomes $`\stackrel{~}{\overline{\rho }}=3\overline{\rho }S\overline{\rho }2\overline{\rho }S\overline{\rho }S\overline{\rho }`$. Second, there are different types of matrix in the non-orthogonal case, one of which transforms with $`\overline{S}^{1/2}`$ and the other with $`\overline{S}^{1/2}`$ (in fact there are two types of orbital and hence four types of matrix); great care must be taken to combine these matrices in the correct fashion, as pointed out by White et al. for the case of the gradient of the energy with the non-orthogonal density matrix. Third, there may be a considerable advantage in choosing the metric so that the transformations shown in equation 12 apply, as this will enable direct comparison with the orthogonal case. The interested reader is referred to Appendix B and references therein for more details.

### C Charge mixing and self-consistency

The question of achieving self-consistency between the charge density and the potential has been examined in great detail over many years, and much is known about efficient implementation. Within Conquest, the direct inversion of the iterative subspace (DIIS) method of Pulay is used. At each iteration, a residual can be defined as:

$$R[\rho ^{\mathrm{in}}]=\rho [\rho ^{\mathrm{in}}]\rho ^{\mathrm{in}},$$
(13)   
where $`\rho [\rho ^{\mathrm{in}}]`$ is the output charge density: that is, the potential arising from $`\rho _{\mathrm{in}}`$ is found (consisting of Hartree and exchange-correlation parts), the Schrödinger equation is solved, and the output charge density generated from the resultant wavefunctions. Clearly, the aim is to reduce $`R[\rho ^{\mathrm{in}}]`$ to as close to zero as possible in the least number of iterations. The simplest possible technique is to use the output charge density from one cycle ($`\rho _n^{\mathrm{out}}=\rho [\rho _n^{\mathrm{in}}]`$) as the input for the next cycle: $`\rho _{n+1}^{\mathrm{in}}=\rho _n^{\mathrm{out}}`$; however, this is potentially rather slow, and prone to the phenomenon known as ‘charge sloshing’, where long wavelength variations of the charge in the unit cell dominate the self-consistency procedure. There are in fact cases where this simple method fails to work at all, and self-consistency is never reached. This is clearly unacceptable.

Better than this is to perform a linear mix of the two previous charge densities, so that:

$$\rho _{n+1}^{\mathrm{in}}=(1\lambda )\rho _{n1}^{\mathrm{out}}+\lambda \rho _n^{\mathrm{out}}.$$
(14)   
The value of lambda can be found easily. If the residuals (defined above in equation 13) are treated as vectors (considering the value at each spatial position as an entry in the vector), then scalar products can be formed between residuals, and the norm of a residual can be found as $`\left(R[\rho _{n+1}^{\mathrm{in}}]R[\rho _{n+1}^{\mathrm{in}}]\right)^{1/2}`$. The optimum value of $`\lambda `$ is found by minimising this norm with respect to $`\lambda `$. This gives:

$$\lambda =\frac{R[\rho _n^{\mathrm{in}}]R[\rho _n^{\mathrm{in}}]R[\rho _{n1}^{\mathrm{in}}]}{R[\rho _n^{\mathrm{in}}]R[\rho _{n1}^{\mathrm{in}}]R[\rho _n^{\mathrm{in}}]R[\rho _{n1}^{\mathrm{in}}]}.$$
(15)   
This procedure can be generalised to more than two previous densities, which can give significant benefits, as described elsewhere. It is often important to mix in a small amount of the input charge densities, as well as performing the mixing given above.

As well as combining charge densities in the optimum manner, it is important to suppress the phenomenon of ‘charge sloshing’; an ideal way to do this is to use Kerker preconditioning. Here a scaling is applied in reciprocal space to the residual:

$$R[\rho _j^{\mathrm{in}}]=R[\rho _j^{\mathrm{in}}]\times \frac{q^2}{q^2+q_0^2},$$
(16)   
where $`q`$ is a reciprocal space vector, and $`q_0`$ is chosen suitably (a value close to $`2\pi /a_0`$, where $`a_0`$ is a lattice vector, is appropriate). This scaling is an approximation to the inverse dielectric function, and enables fast and robust iteration to a self-consistent charge density and potential.

### D Pre-conditioning localised orbital variation

Now that we have described the robust and efficient search for the ground state density matrix (for a given set of localised orbitals) and the fast iteration to a self-consistent charge density and potential, we must consider the variation of the localised orbitals.

As is the case for minimisation problems in many areas of science, Conquest suffers from ill conditioning in the search for the ground state when varying the localised orbitals. Ill conditioning occurs if the function being minimised has a wide range of curvatures. For a general function $`f(x_1,x_2,\mathrm{},x_N)`$, dependent upon variables $`\{x_i\}`$, the curvature matrix (or Hessian) can be defined as $`C_{ij}=^2f/x_ix_j`$. If the eigenvalues $`\lambda _n`$ of $`C_{ij}`$ span a wide range, then the surfaces of constant $`f`$ are elongated (illustrated in Figure 4), and conventional techniques such as conjugate gradients will become very inefficient. It is known that the number of iterations required by conjugate gradients is proportional to $`(\lambda _{\mathrm{max}}/\lambda _{\mathrm{min}})^{1/2}`$, where $`\lambda _{\mathrm{max}}`$ and $`\lambda _{\mathrm{min}}`$ are the maximum and minimum eigenvalues of $`C_{ij}`$.

While ill conditioning is a widespread problem, the solution depends on the specific situation. Conventional first principles calculations which use plane waves as a basis set have been recognised for many years to suffer from ill conditioning, and it turns out that ill conditioning found in $`𝒪(N)`$ calculations is closely related to this; we will review the plane wave ill conditioning before describing the $`𝒪(N)`$ ill conditioning.

The plane wave energy functional (eq. A4) has large curvatures associated with high wavevector G, because of the form of the kinetic energy:

$$E_{\mathrm{kin}}=2\underset{i}{}f_i\underset{𝐆}{}\frac{\mathrm{}^2G^2}{2m}c_{i𝐆}^2,$$
(17)   
so that the energy is proportional to $`G^2`$. This first type of ill conditioning is easily cured (essentially by scaling $`c_{i𝐆}`$ by a factor of $`(1+G^2/G_0^2)^{1/2}`$), and is referred to as ‘length scale ill conditioning’, as it comes from the variation of curvature with length scale.

Another type of ill conditioning seen in conventional techniques is associated with the invariance of $`E_{\mathrm{tot}}`$ under a unitary transformation of the orbitals. If the occupation numbers $`f_i`$ are all either zero or one, then $`E_{\mathrm{tot}}`$ is exactly invariant under transformations such as:

$$\psi _i^{}=\underset{j}{}U_{ij}\psi _j,$$
(18)   
where $`U_{ij}`$ is a unitary matrix. If the occupancies deviate slightly from zero or one, however, the exactness of the invariance is broken, and the energy changes slightly. Some of the eigenvalues of the Hessian will go from exactly zero (under the exact transformation) to very small, leading to poor conditioning. We shall refer to this as ‘superposition ill conditioning’. In conventional techniques this is cured by performing a rotation of the wave functions so that the Hamiltonian becomes diagonal in the subspace spanned by the occupied states.

A final type of ill conditioning found in conventional methods arises with variable occupation numbers, and is associated with eigenvalues whose energies are well above the Fermi energy, which will have very small occupation numbers. Variations of the $`\psi _i`$ associated with these small occupation numbers will have little effect on the value of $`E_{\mathrm{tot}}`$, and lead to small values of the curvature. Since the variations of these eigenvalues are almost redundant in the minimisation, we refer to this as ‘redundancy ill conditioning’.

All three of these forms can cause difficulties within $`𝒪(N)`$ techniques, though typically the specific solution will vary. For instance, it is clear that variations of the localised orbitals, $`\{\varphi _{i\alpha }\}`$, will have different length scales, and will suffer from length scale ill conditioning. This is easily cured in the same way as for plane waves, as has been recently demonstated, though the methodology is somewhat different. As a demonstration of the efficacy of this preconditioning, Figure 5 shows the convergence of Conquest with and without length scale preconditioning for three different region radii for the localised orbitals. Clearly this problem becomes significantly worse for larger regions.

Superposition ill conditioning is associated with the linear mixing of localised orbitals. It is easily shown that linear mixing of the functions on the same atom leaves $`E_{\mathrm{tot}}`$ unchanged, and so will not cause ill conditioning. Variations of the localised orbital $`\varphi _{i\alpha }`$ such as:

$$\varphi _{i\alpha }^{}=\varphi _{i\alpha }+\underset{j\beta ,ji}{}c_{j\beta }\varphi _{j\beta }$$
(19)   
are rather more interesting. Strictly, these are not possible, as the localised orbitals are constrained to be zero outside their localisation regions. However, once the region radii become large, there will be variations which almost fulfil this criterion. It is the small eigenvalues of the Hessian of $`E_{\mathrm{tot}}`$ associated with this almost perfect mixing which will cause superposition ill conditioning. It is perfectly possible to cure this, however, as the form of the variations can be written down. We have developed a method to precondition these variations, and are testing it. It will be described in a future publication.

Finally, we come to redundancy ill conditioning. Just as in conventional calculations this occurs when the occupation numbers $`f_i`$ are very small, this may occur in $`𝒪(N)`$ techniques when the number of localised orbitals $`\varphi _{i\alpha }`$ is more than half the number of electrons. It is desirable, if not essential, to be able to work with an extended number of orbitals; for instance, in group IV elements, the natural basis will consist of four orbitals, roughly corresponding to one $`s`$ and three $`p`$ orbitals. (It is worth noting that Kim, Mauri and Galli have found that it is essential to have more orbitals than filled bands to avoid local minima in the energy functional in a related $`𝒪(N)`$ scheme.) As before, we believe that this form of ill conditioning can be removed by suitable preconditioning, but detailed techniques have yet to be formulated.

## IV Future directions

Having summarised the techniques involved in Conquest, it is appropriate to look to the near future, and consider the directions in which the project is going.

### A Forces

The issue of forces is a key one for any electronic structure technique; if relaxation of ions or molecular dynamics are to be performed then the analytical forces must agree with the gradient of the energy. Conquest has been specifically constructed so that, provided small Pulay-type corrections are included, the forces are guaranteed to be consistent with the energy. Pulay corrections are required because the B-spline basis functions move with the atoms, and are easily found, as will be described elsewhere. This means that the relaxation of the system to mechanical equilibrium and the generation of time-dependent ionic trajectories will be feasible in $`𝒪(N)`$ DFT calculations.

### B Efficient choices for representation of localised orbitals

The present choice of basis for representing the localised orbitals has been described in Section II A. However, there are two good reasons for changing this basis somewhat. First, there is the problem of ripples in the energy caused by the numerical integration; this is due to a lack of translational invariance with respect to the integration grid. If the rapidly varying parts of the localised orbitals (i.e. the core regions, which do not alter greatly during a calculation) could be represented in a more efficient form, then the integrals could be performed analytically, significantly reducing ripples. Second, there would be great value in being able to perform ab initio tight binding calculations with the code (or even to have certain parts of a supercell treated with full DFT, while others were treated with ab initio tight binding). For these reasons, we are planning to move to a mixed basis (seen recently in conventional pseudopotential calculations) which combines pseudo-atomic orbitals (possibly of the form of Sankey and Niklewski) with a coarser blip grid. This will suppress the ripples with respect to the integration grid, and give the flexibility to model different parts of the system with appropriate accuracy.

### C Finding the density matrix for metals

The question of modelling metals is much harder than insulators or semiconductors for the $`𝒪(N)`$ methods described above, simply because the density matrix is more delocalised in metals, meaning that the cutoff applied to $`L`$ (and hence to $`K`$) must be much larger to obtain accurate results. If the metal is close packed, as is frequently the case, this entails a rapid increase in the number of elements in the density matrix, and a slowing down of variational techniques; this is discussed elsewhere.

One approach is to reduce the range of the density matrix in metals by introducing an artificial electronic temperature, which broadens the Fermi occupation function, and localises the density matrix. The drawback is a potentially large electronic entropy contribution to the energy; however, there is a scheme for extrapolating the results back to zero electronic temperature. It has been shown that even with such a scheme the variational density matrix method described above is inefficient; the hybrid method described in Section III A improves the efficiency. More efficient are recursion-based methods, such as the Fermi Operator Expansion or the Bond Order Potential. However, these have the disadvantage of not being strictly variational.

A further possibility which has emerged recently is to use a series of nested Hilbert spaces to find the exact zero temperature density matrix. This method has the advantage that no approximation is being made to remove the electronic entropy contributions, and allows high precision calculations on metals. It is, however, still in development, and practical demonstrations have yet to be published.

## V Conclusions

The recent developments which have been outlined above show that the future of $`𝒪(N)`$ ab initio techniques is extremely bright. We have shown that the localisation of the density matrix gives the framework within which these methods can be constructed, and have given details of the implementation of one such code, Conquest. The examples presented show that this method is practical, and that the spatial cutoffs required for accuracy are small enough to make the calculations perfectly feasible. The search for the ground state has been addressed, and methods for making this search more robust and efficient have been discussed. The remaining tasks for the Conquest code have been described, and the way forward for all of them is clear. The most important conclusion to draw from this body of work is that $`𝒪(N)`$ DFT methods actually work. Indeed, these methods are being demonstrated in practical calculations. Our group is working towards practical application of the Conquest code to large-scale problems.

## Acknowledgements

We are happy to acknowledge useful discussions with David Manolopoulos, Peter Haynes, Chris Goringe and Ed Hernández. Conquest has been developed within the framework of the UK Car-Parrinello consortium, which is supported by the EPSRC grant GR/M01753. The work of MJG is financially supported by CCLRC and GEC.

## A The pseudopotential method

As Conquest uses the standard pseudopotential method, it is important to recall the salient facts within this formalism; there are many excellent reviews elsewhere which give more detail.

In the pseudopotential method, only the valence electrons are considered, and their interaction with the ionic cores is replaced by a pseudopotential, $`v(r)`$. This means that in fact the solutions of the Schrödinger equation are pseudo-wavefunctions, and the charge density, $`n(𝐫)`$, is the pseudo-density of the valence electrons. The energy arising from the interaction between the cores and the valence electrons is given by:

$$E_{ei}=d𝐫V(𝐫)\rho (𝐫),$$
(A1)   
where $`V(𝐫)`$ is found as the sum over the ionic pseudopotentials:

$$V(𝐫)=\underset{i}{}v(𝐫𝐑_i),$$
(A2)   
with $`𝐑_i`$ the core positions. In general, the pseudopotential is non-local, that is $`v(𝐫,𝐫^{})`$.

In conventional pseudopotential ab initio techniques, the wavefunctions are often expanded in terms of plane waves:

$$\psi _i=\underset{𝐆}{}c_{i𝐆}\mathrm{exp}(i𝐆𝐫),$$
(A3)   
where G is a reciprocal lattice vector. The total energy is then minimised with respect to the coefficients $`c_{i𝐆}`$. Frequently, particularly in metals, variable occupation numbers are allowed, so that the wave function $`\psi _i`$ has an occupation $`f_i`$. This gives a total energy function:

$$E_{\mathrm{tot}}=E_{\mathrm{tot}}(\{c_{i𝐆}\},\{f_i\}).$$
(A4)   
## B Metrics and minimisation in non-orthogonal basis sets

When working with a non-orthogonal basis, care must be taken with notation. It is common to use raised and lowered indices to distinguish between vectors and matrices which transform differently; this has been introduced to electronic structure calculations by Ballentine and Kolář, who also describe the general formalism.

The eigenfunctions of the Hamiltonian are expanded in a set of non-orthogonal, localised, atom-centred orbitals $`\{\varphi _\alpha (𝐫)\}`$, where $`\alpha `$ runs over all orbitals on all atoms. These orbitals span a Hilbert space $`𝒱`$, and define an overlap matrix which is given by $`S_{\alpha \beta }=\varphi _\alpha \varphi _\beta `$. The inverse overlap matrix, $`\left(S^1\right)^{\alpha \beta }`$, is defined by the relation:

$`{\displaystyle \underset{\beta }{}}\left(S^1\right)^{\alpha \beta }S_{\beta \gamma }`$ $`=`$ $`\delta _\gamma ^\alpha `$ (B1)   
$`=`$ $`1\mathrm{i}\mathrm{f}\alpha =\gamma `$ (B2)   
$`=`$ $`0\mathrm{i}\mathrm{f}\alpha \gamma .`$ (B3)   
A dual space, $`𝒱^{}`$, exists which is spanned by the orbitals $`\varphi ^\alpha =_\beta \left(S^1\right)^{\alpha \beta }\varphi _\beta `$. These two sets of vectors are bi-orthogonal, that is $`\varphi ^\alpha \varphi _\beta =\delta _\beta ^\alpha `$. It is important to note that contraction can only be carried out over indices which are opposed, while addition can only be carried out between tensors for which indices agree. The vectors in the original space are called covariant, while the vectors in the dual space are contravariant.

A convenient choice for the metric of $`𝒱`$ is $`S^1`$, so that $`_\beta \varphi _\alpha \left(S^1\right)^{\gamma \beta }\varphi _\beta =\delta _\alpha ^\gamma `$. An equivalent choice of metric for $`𝒱^{}`$ is $`S_{\alpha \beta }`$. These operate to change a vector in one space to the vector in another; thus a proper scalar product within a space can be formed by incorporation of the metric. A covariant operator can be represented as an outer product:

$$\widehat{A}=\underset{\alpha ,\beta }{}\varphi _\alpha A^{\alpha \beta }\varphi _\beta .$$
(B4)   
Then the scalar product of two covariant operators is written:

$`(𝐀,𝐁)`$ $`=`$ $`{\displaystyle \underset{\alpha ,\beta ,\gamma ,\delta }{}}\varphi _\delta \varphi _\beta A^{\beta \alpha }\varphi _\alpha \varphi _\gamma B^{\gamma \delta },`$ (B5)   
$`=`$ $`{\displaystyle \underset{\alpha ,\beta ,\gamma ,\delta }{}}S_{\delta \beta }A^{\beta \alpha }S_{\alpha \gamma }B^{\gamma \delta }`$ (B6)   
$`=`$ $`Tr[A^{}SBS].`$ (B7)   
It is important to note that the product of a covariant and a contravariant pair (such as $`H`$ and $`\rho `$) is invariant with basis set. Traditionally, the Hamiltonian is taken as covariant and the density matrix as contravariant.

### 1 Variation of L

The innermost part of Conquest consists of the minimisation of the energy with respect to the elements of the matrix $`L^{i\alpha j\beta }`$ with fixed support functions. This is achieved by performing line minimisations along directions supplied by the conjugate gradients algorithm (with, on occasion, a correction for maintaining the electron number constant). As pointed out recently by White et al., the gradient of the energy with respect to $`L^{i\alpha j\beta }`$, $`\mathrm{\Omega }`$, (which is used as the search direction in the minimisation) is actually covariant, while $`L^{i\alpha j\beta }`$ is contravariant; this means that the gradient must be transformed to a contravariant matrix before being combined with the density matrix: $`(S^1)\mathrm{\Omega }(S^1)`$. But there is more to the problem of minimisation than just this; conjugate gradients assumes the following relations:

$`𝐠_{i+1}`$ $`=`$ $`f(𝐏_{i+1}),`$ (B8)   
$`𝐡_{i+1}`$ $`=`$ $`𝐠_{i+1}+\gamma _i𝐡_i,`$ (B9)   
$`\gamma _i`$ $`=`$ $`{\displaystyle \frac{𝐠_{i+1}𝐠_{i+1}}{𝐠_i𝐠_i}},`$ (B10)   
where $`𝐠_i`$ is the gradient of the function at step $`i`$ and $`𝐡_i`$ is the search direction (which is conjugate to the previous search directions) at step $`i`$. Clearly in the formation of $`\gamma _i`$, care must be taken to ensure that the product is tensorially correct, otherwise the choice of the new search direction will be wrong, so the correct formula for $`\gamma _i`$ is:

$$\gamma _i=\frac{𝐠_{i+1}(S^1)𝐠_{i+1}(S^1)}{𝐠_i(S^1)𝐠_i(S^1)},$$
(B11)   
(This is discussed at more length in Section 2.7 of Ref., where the bi-conjugate gradient method is described. This degree of complexity is not needed here.) Similar care both with the correct nature of gradients and with the search directions must be taken when varying the localised orbitals.