Next: Density-kernel optimisation Up: Overview Previous: Overview

## Density-matrix formulation

The set of Kohn-Sham orbitals {ψn(r)} provides a complete description of the fictitious system of non-interacting particles in DFT. An equivalent description may also be given by the single-particle density-matrix ρ(r,r'),

 (4)

which possesses the property of idempotency, namely

 (5)

which requires the orthonormality of the orbitals {ψn(r)} according to Eq. 3 and the Aufbau principle of singly occupying all states up to the chemical potential, which itself follows straightforwardly from the Pauli exclusion principle. The density-matrix is thus the position representation of the projection operator onto the space of occupied states {$\hat{\rho}$}.

The density n(r) may be obtained from the diagonal elements of the density-matrix,

 (6)

where the factor of two again accounts for spin degeneracy (assuming no spin polarisation). The total energy of the non-interacting system Es may be obtained via

 (7)

and thus the energy of the real interacting system calculated by applying the usual double-counting corrections to the Hartree and exchange-correlation terms. The solution to Eq. 1 may therefore be found by minimising the energy with respect to the density-matrix, subject to the constraints of idempotency (5) and normalisation,

 (8)

i.e. the density-matrix corresponds to a system of Ne electrons.

Since the number of occupied states is directly proportional to N and each state extends over the whole system (Fig. 1), the amount of information in the density-matrix defined by Eq. 4 therefore scales as N2. Any calculation involving manipulation of this density-matrix will therefore scale quadratically with system-size at best. In order to obtain a linear-scaling method, it is necessary to exploit the nearsightedness of many-body quantum mechanics [15] mentioned above.

Both analytical [30,31] and numerical [32] studies have demonstrated that, for an insulating system, both the Wannier functions and density-matrix decay exponentially, so that

 (9)

This means that for any given position r the density-matrix ρ(r,r') differs significantly from zero only for points r' within a finite volume around r. Since the decay rate γ depends only on the energy gap between the highest occupied and lowest unoccupied states, this volume is independent of system-size. Thus the total amount of significant information in the density-matrix only scales linearly with N.

Figure 2: Three localised orbitals (non-orthogonal generalised Wannier functions) generated by ONETEP for the same oligopeptide molecule as Fig. 1.

In practice this is exploited by writing the density-matrix in separable form:

 (10)

where the {φα(r)}, illustrated in Fig. 2, are a set of spatially localised non-orthogonal functions which span a superspace of the Hilbert space of the set of occupied Kohn-Sham orbitals {ψn(r); εn < μ} and have been called support functions [33] or non-orthogonal generalised Wannier functions (NGWFs) [34] in this work. The matrix Kαβ, known as the density-kernel [35] is in fact the representation of the density-matrix in the set of duals of the NGWFs {φα(r)} defined by

 (11)

The total amount of information contained in the density-matrix defined by Eq. 10 must of course still scale as N2. However the advantage of this form is that it allows the nearsightedness to be exploited with the use of spatial cut-offs.

First, the NGWFs which are exponentially localised are truncated, by allowing them to be non-vanishing only in spherical regions of fixed radii {rα} and centred at positions {Rα. In ONETEP, a number of NGWFs are associated with each atom in the system, so that the regions are centred on atoms and their radii essentially depend only on the atomic species.

Second, the density-kernel is required to be a sparse matrix, by discarding elements Kαβ corresponding to NGWFs centred further apart than some cut-off rK. Note that since the density-kernel is related directly to the duals of the NGWFs, rather than the NGWFs themselves, the density-kernel cut-off rK is not simply the sum of the radii rα+rβ. This is a consequence of the non-orthogonality of the NGWFs.

The approach adopted in ONETEP of truncating the NGWFs and predetermining the sparsity pattern of the density-kernel contrasts with approaches based on thresholding [36] which are common in methods based on Gaussian basis sets. The former approach has the advantage that the overlap matrix for NGWFs remains strictly positive definite, whereas the latter approach is more flexible and hence potentially more efficient.

Imposing spatial cut-offs on the NGWFs and the density-kernel results in a density-matrix whose information content scales linearly with system-size, an approximation which is controlled by adjusting the {rα} and rK. In practice, these cut-offs are increased until the desired physical properties of the system converge.

Next: Density-kernel optimisation Up: Title page Previous: Overview