4.5 Requirements for linear-scaling methods

We have shown how it is possible to reformulate density-functional theory in terms of the single-particle density-matrix, and the constraints which must be obeyed by ground-state density-matrices. However, in the coordinate representation, we note that the density-matrix is a function of two position variables $\rho({\bf r},{\bf r'})$ and thus contains an amount of information which scales as the square of the system-size (as of course it must since it contains all of the information in the Kohn-Sham orbitals, which are functions of one position but with the number of occupied orbitals also scaling linearly with system-size). To obtain a linear-scaling method, it is necessary to impose some further restrictions on the density-matrix.

4.5.1 Separability

In practice we do not wish to deal with a function of six variables (i.e. two three-dimensional positions). From the factorisation property of idempotent density-matrices, or the definition of the ground-state density-matrix in terms of the Kohn-Sham orbitals;

$$\rho({\bf r},{\bf r'}) = \sum_i \psi_i({\bf r}) f_i~\psi_i^{\ast}({\bf r'}) ,$$ (4.58)

we see that it is possible to consider separable density-matrices described in terms of some auxiliary orbitals. The general form of a separable density-matrix in terms of orbitals $\{ \varphi_i \}$ is

$$\rho({\bf r},{\bf r'}) = \sum_{ij} \varphi_i({\bf r}) R_{ij} \varphi_j^{\ast}({\bf r'}) .$$ (4.59)

Although it is not necessary for the auxiliary orbitals to be orthonormal, in the case when they are, we can consider this general form as simply a unitary transformation of the Kohn-Sham expression (4.58):

$$\psi_i({\bf r}) = \sum_j \varphi_j({\bf r}) U_{ji}$$ (4.60)

where $U$ is the unitary matrix which diagonalises $R$:

$$f_i = (U^{\dag }RU)_{ii} \qquad \rm {(no~summation).}$$ (4.61)

When the auxiliary orbitals are not orthonormal, then they can be viewed as a more general linear combination of the Kohn-Sham orbitals (involving both a unitary and Löwdin transformation) which is described in section 4.6. Whichever case applies, there is no loss of generality here as all idempotent matrices can always be expressed in this way, and these are the density-matrices of interest to us.

4.5.2 Spatial localisation

Kohn [136] has proved that in one-dimensional systems with a gap, a set of exponentially decaying Wannier functions can be found in the tight-binding limit, and that this localisation is related to the square-root of the gap. His method is not easily generalised to higher numbers of dimensions, and so until recently the exact nature of the Wannier functions in general three-dimensional systems was unknown, although it was anticipated that they would decay exponentially [137,138,139]. More recent numerical and analytical studies of the localisation of the density-matrix showed the decay to be exponential and again related to the square-root of the gap [140,141], thus supporting the general validity of Kohn's result. Very recently, however, Ismail-Beigi and Arias [142] have argued that in the weak-binding limit the exponential decay varies linearly with the gap. What is now certain is that the Wannier functions and density-matrix decay exponentially in systems with a gap, and that this decay is more rapid in systems with larger gaps.

Wannier functions are simply a unitary transformation of Bloch wave-functions with respect to the complementary variables of Bloch wave-vector and lattice vector. Let $\psi_{n{\bf k}}({\bf r})$ be the normalised Bloch wave-function for the $n$th band with wave-vector ${\bf k}$. Then the corresponding Wannier function for that band $w_{n{\bf R}}({\bf r})$ is defined by [143]

$$w_{n{\bf R}}({\bf r}) = \left( \frac{\Omega_{\mathrm{cell}}}... ...i_{n{\bf k}}({\bf r}) \exp(-{\mathrm i} {\bf k} \cdot {\bf R})$$ (4.62)

and naturally the inverse relation holds:

$$\psi_{n{\bf k}}({\bf r}) = \left( \frac{\Omega_{\mathrm{cell... ..._{n{\bf R}}({\bf r}) \exp({\mathrm i} {\bf k} \cdot {\bf R}) .$$ (4.63)

The properties of Wannier functions are that they are localised in different cells (labelled by lattice vector ${\bf R}$) and are orthonormal:

$$\int {\mathrm d}{\bf r}~w_{n{\bf R}}^{\ast}({\bf r})~w_{n{\bf R'}}({\bf r}) = \delta_{{\bf R}{\bf R'}} .$$ (4.64)

The single-particle density-matrix in the case of full ${\bf k}$-point sampling $\varrho({\bf r},{\bf r'})$ is given by

$$\varrho({\bf r},{\bf r'}) = \sum_i f_i ~\frac{1}{\Omega_{\ma... ...k} ~\psi_{i{\bf k}}({\bf r}) \psi_{i{\bf k}}^{\ast} ({\bf r'})$$ (4.65)

(in which we have assumed that we are dealing with an insulator with completely full or empty bands) which is a trace over wave-vector ${\bf k}$ and thus invariant under unitary transformation so that it can also be written

$$\varrho({\bf r},{\bf r'}) = \sum_i f_i \sum_{\bf R} w_{i{\bf R}}({\bf r}) w_{i{\bf R}}^{\ast}({\bf r'}) .$$ (4.66)

Thus if the Wannier function $w_{i{\bf R}}({\bf r})$ is vanishing when $\vert{\bf r}-{\bf R}\vert$ is large, then when $\vert{\bf r}-{\bf r'}\vert$ is large the density-matrix must also vanish in the same way since it is impossible in that case for both $\vert{\bf r}-{\bf R}\vert$ and $\vert{\bf r'}-{\bf R}\vert$ to be small. Thus we expect that

$$\rho({\bf r},{\bf r'}) \rightarrow 0 \qquad {\rm as}~\vert{\bf r}-{\bf r'}\vert \rightarrow \infty$$ (4.67)

where at zero temperature the decay is exponential in insulators and algebraic in metals.

We can exploit this long-range behaviour to obtain a linear-scaling method: we introduce a spatial cut-off $r_{\mathrm{cut}}$ and require that the density-matrix be strictly zero when the separation of its arguments exceeds this cut-off;

$$\rho({\bf r},{\bf r'}) = 0 , \qquad \vert{\bf r}-{\bf r'}\vert > r_{\mathrm{cut}} ,$$ (4.68)

so that the density-matrix now only contains an amount of information which scales linearly with system-size. Imposing this cut-off naturally restricts the variational freedom of the density-matrix, so it will be necessary to converge the ground-state energy with respect to this parameter in real calculations. Using the separable form above, we can impose this restriction by requiring the auxiliary orbitals to be localised in space (i.e. vanishing outside a certain region of space) and by making the matrix $R$ sparse, so that elements of $R$ corresponding to orbitals localised in regions separated by more than the spatial cut-off $r_{\mathrm{cut}}$ are automatically set to zero. The localised nature of the auxiliary orbitals requires a localised basis-set to describe them, and this is the subject of chapter 5.