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 {fn} 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 rK for a Ti38O76 cluster [46].

Figure 4 shows a particular example of how the total energy varies as the density-kernel cut-off rK 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 rK. For example, in crystalline silicon the errors from truncating the density-kernel to rK = 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