8.2 Density-matrices from Kohn-Sham orbitals

This method is based upon work on the projection of plane-wave calculations onto atomic orbitals [166], which has been used to analyse atomic basis-sets [167] and obtain local atomic properties from the extended Kohn-Sham orbitals [168]. Here we review the method, for the special case of a $\Gamma$-point Brillouin zone sampling.

8.2.1 Projecting plane-wave eigenstates onto support functions

The plane-wave eigenstates are denoted $\vert \psi_i \rangle$ and the support functions are denoted $\vert\phi_{\alpha}\rangle$. The states obtained by projecting the plane-wave eigenstates onto the space spanned by the support functions are denoted $\vert \xi_{\alpha} \rangle$. As in section 4.6 we also introduce the dual states $\vert \phi^{\alpha} \rangle$ and $\vert \xi^{\alpha} \rangle$ with the properties outlined below.

$$S_{\alpha \beta} = \langle \phi_{\alpha} \vert \phi_{\beta} \rangle \qquad \qquad \Sigma_{\alpha \beta} = \langle \xi_{\alpha} \vert \xi_{\beta} \rangle$$ (8.5)

$$\vert \phi^{\alpha} \rangle = \vert \phi_{\beta} \rangle S_{\beta \alpha}^{-1} \qquad \qquad \vert \xi^{\alpha} \rangle = \vert \xi_{\beta} \rangle \Sigma_{\beta \alpha}^{-1}$$ (8.6)

$$\langle \phi^{\alpha} \vert \phi_{\beta} \rangle = \langle \phi_{\alpha} \vert \phi^{\beta} \rangle = \delta_{\alpha}^{\beta} \qquad \qquad \langle \xi^{\alpha} \vert \xi_{\beta} \rangle = \langle \xi_{\alpha} \vert \xi^{\beta} \rangle = \delta_{\alpha}^{\beta}$$ (8.7)

The projection operator onto the subspace spanned by the support functions is defined by

$${\hat P} = \vert \phi_{\alpha} \rangle \langle \phi^{\alpha}... ...ha} \rangle S_{\alpha \beta}^{-1} \langle \phi_{\beta} \vert .$$ (8.8)

A spilling parameter ${\cal S}$ can be defined to measure how much the subspace spanned by the plane-wave eigenstates falls outside the subspace spanned by the support functions. Minimising this quantity is one method of optimising the choice of support functions, and is described for the case of the spherical-wave basis (chapter 5) in section 8.2.3.

$${\cal S} = \frac{1}{N_{\mathrm b}} \langle \psi_i \vert \bigl( 1 - \hat P \bigr) \vert \psi_i \rangle$$ (8.9)

where $N_{\mathrm b}$ is the number of bands (labelled $i$) considered. The density-operator is then defined by

$$\hat \rho = \sum_{\alpha}^{\mathrm{occ}} \vert \xi_{\alpha} \rangle \langle \xi^{\alpha} \vert$$ (8.10)

where the sum is taken over occupied bands only. Substitution of the results given above then yields the following expression for the density-kernel:

$$K^{\alpha \beta} = \langle \phi^{\alpha} \vert \hat \rho \ve... ...} \langle \psi_j \vert \phi_{\mu} \rangle S_{\mu \beta}^{-1} .$$ (8.11)

Defining the rectangular matrix $L$ as

$$L_{\lambda i} = \langle \phi_{\lambda} \vert \psi_i \rangle$$ (8.12)

we give an expression for the matrix $\Sigma$ in terms of $L$ and $S$:

$$\Sigma_{\alpha \beta} = \langle \xi_{\alpha} \vert \xi_{\bet... ...nu}^{-1} S_{\nu \beta} = (L^{\dag } S^{-1} L)_{\alpha \beta} .$$ (8.13)

We can thus minimise the spilling parameter ${\cal S}$ to optimise our choice of support functions, and then calculate $K$ to obtain all of the information required to start a linear-scaling calculation.

8.2.2 Obtaining auxiliary matrices

In the case when the density-kernel $K$ is expanded in terms of an auxiliary matrix $T$ e.g. in order to construct a positive semi-definite density-matrix, it is necessary to be able to calculate the auxiliary matrix $T$ which corresponds to a given density-kernel $K$ by

$$K = T T^{\dag } .$$ (8.14)

This can be achieved by minimising the function ${\cal I}(T)$ given by

$${\cal I}(T) = {\rm Tr} \left[ (K - T T^{\dag })^2 \right]$$ (8.15)

whose derivative with respect to $T$ is

$$\frac{\partial {\cal I}(T)}{\partial T^{\alpha \beta}} = -4 \left[ T^{\dag } (K - T T^{\dag }) \right]_{\beta \alpha} .$$ (8.16)

This derivative vanishes at the minimum, and so we find that the matrix $T$ which minimises ${\cal I}(T)$ is the desired auxiliary matrix (the solution $T=0$ corresponds to a local maximum). We therefore choose to minimise ${\cal I}(T)$ by the conjugate gradients method to obtain the auxiliary matrix.

8.2.3 Optimising the support functions

As mentioned above, we can optimise our choice of support functions by minimising the spilling parameter ${\cal S}$. We describe this process here when the support functions are themselves described in terms of a localised basis:

$$\vert {\phi}_{\alpha} \rangle = \sum_{n \ell m} c^{n \ell m}_{(\alpha)} \vert \chi_{\alpha , n \ell m} \rangle .$$ (8.17)

The spilling parameter can be written in terms of the matrices $L$ and $S$ by:

$${\cal S} = \frac{1}{N_{\mathrm b}} \langle \psi_i \vert \bigl( 1 - \hat P \b... ...e \psi_i \vert \phi_{\alpha} \rangle \langle \phi^{\alpha} \vert \psi_i \rangle$$

$$= 1 - \frac{1}{N_{\mathrm b}} L_{i \alpha}^{\dag } S^{-1}_{\alpha \beta} L_{\beta i} = 1 - \frac{1}{N_{\mathrm b}} {\rm Tr}[L^{\dag } S^{-1} L]$$ (8.18)

and we wish to obtain the gradients of ${\cal S}$ with respect to the expansion coefficients $\{ c^{n \ell m}_{(\alpha)} \}$.

$$\frac{\partial {\cal S}}{\partial c^{n \ell m}_{(\alpha)}} =... ...c{\partial L_{ki}}{\partial c^{n \ell m}_{(\alpha)}} \right] .$$ (8.19)

We obtain the derivative of the inverse matrix by differentiating $S^{-1} S = 1$ i.e.

$$\frac{\partial (S^{-1}S)_{\alpha \beta}}{\partial x} = \frac... ... \gamma}^{-1} \frac{\partial S_{\gamma \beta}}{\partial x} = 0$$ (8.20)

which can be rearranged to give

$$\frac{\partial S_{\alpha \beta}^{-1}}{\partial x} = - S_{\al... ...partial S_{\gamma \delta}}{\partial x} S_{\delta \beta}^{-1} .$$ (8.21)

Therefore (no summation over $\alpha$)

$$\frac{\partial {\cal S}}{\partial c^{n \ell m}_{(\alpha)}} = - \frac{1}{N_{\mathrm b}} \left[ \langle \psi_i \vert \chi_{\alph... ... \ell m} \rangle \delta_{\beta \gamma} \right) S_{\gamma k}^{-1} L_{ki} \right.$$

$$\qquad \left. + L_{ij}^{\dag } S_{j \alpha}^{-1} \langle \chi_{\alpha , n \ell m} \vert \psi_i \rangle \right]$$

$$= - \frac{2}{N_{\mathrm b}} {\rm Re} \left[ L_{\alpha \beta}^{\dag ... ...hi_{\beta , n \ell m} \rangle S_{\beta \gamma}^{-1} L_{\gamma \alpha} \right] .$$ (8.22)

In the case of the set of basis functions introduced in chapter 5, the overlap between plane-wave eigenstates and localised basis functions, e.g. $\langle \psi_{\alpha} \vert \chi_{\beta , n \ell m} \rangle$, can be calculated using the expression for the basis function Fourier transform (5.9).

We can use these gradients to minimise the spilling parameter (by the conjugate gradients method) to obtain the set of optimal coefficients $\{ c^{n \ell m}_{(\alpha)} \}$ which define the set of support functions which best span the space of the occupied plane-wave orbitals. The final minimum spilling parameter value also gives an estimate of the quality of the basis-set being used.