Next: NGWF optimisation Up: Overview Previous: Density-matrix formulation

## Density-kernel optimisation

Most linear-scaling methods fall into two categories, typically along the same lines as the basis set division above. The first is equivalent to optimising the density-kernel only [22,37,38,39] for a fixed, but potentially large, set of local orbitals and the second involves both density-kernel and NGWF optimisation [26,34,40]. These two approaches are further compared in the next section, while this section concerns the common point of density-kernel optimisation.

Functionally at least, this stage of the calculation proceeds along the lines of an *ab initio* tight-binding [41] calculation and linear-scaling methods for tight-binding [42] may be applicable. While self-consistency may be treated independently of density-kernel and NGWF optimisation [43], the approach in ONETEP is closely related to the ensemble DFT method for metallic systems [44] and performs both density-kernel and NGWF optimisation self-consistently.

To impose the idempotency constraint a combination of two methods is used: a penalty-functional method [45] and the method of Li, Nunes and Vanderbilt [47], and independently Daw [48], based on McWeeny's purification transformation [35]. The latter has been implemented in both its original non-orthogonal formulation [49] and the Millam and Scuseria variant [50].

Figure 3: The purification transformation illustrated in terms of its effect on the density-matrix eigenvalues (occupation numbers).

The purification transformation may be defined in terms of an auxiliary matrix σ(**r**,**r**') by:

ρ(r,r') = 3σ^{2}(r,r') - 2σ^{3}(r,r'). | (12) |

Clearly ρ(**r**,**r**') and σ(**r**,**r**') commute and so the transformation is best illustrated as in Fig. 3 in terms of their eigenvalues, which for ρ(**r**,**r**') are the occupation numbers {f_{n}} of Eq. 4. Under repeated iteration of this transformation, eigenvalues around zero vanish quadratically and those close to unity converge to the maximum there. Expressing the density-matrix in this way and optimising the auxiliary matrix thus applies the constraint of idempotency, so long as the occupation numbers remain in the interval {$\left[ \frac{1-\sqrt{5}}{2} , \frac{1+\sqrt{5}}{2} \right]$} where the transformation is stable. Should the calculation become unstable (generally at the start of a calculation) then the globally convergent penalty-functional method is used instead.

The purification transformation may also be derived as the result of a steepest descents minimisation of the positive definite functional P[ρ] given by

(13) |

which gives a quantitative measure of the departure of ρ(**r**,**r**') from idempotency (for nearly-idempotent density-matrices it is proportional to the mean square deviation of the occupation numbers from zero or unity). P[ρ] may thus be used as a penalty functional, some proportion of which is added to the energy functional to enforce idempotency [45].

Figure 4: Total energy convergence with respect to density-kernel cut-off r_{K} for a Ti_{38}O_{76} cluster [46].

Figure 4 shows a particular example of how the total energy varies as the density-kernel cut-off r_{K} is increased, and this behaviour is qualitatively similar to that observed in all insulating systems. Physical properties based on energy differences converge even more rapidly with r_{K}. For example, in crystalline silicon the errors from truncating the density-kernel to r_{K} = 10 Å are 0.3% and 0.1% in the lattice parameter and bulk modulus respectively [13].

Next: NGWF optimisation Up: Title page Previous: Density-matrix formulation