Derivation of the correction

In this section we present a new method to perform total energy calculations using a penalty functional to enforce idempotency approximately. We define a generalised energy functional $Q[{\rho};\alpha]$ for trial density-matrices ${\rho}$ by

$$Q[{\rho};\alpha] = E[\rho] + \alpha P[\rho] .$$ (6.15)

$P[\rho]$ represents the penalty functional, which in this new method is required to be analytic at all points so that efficient minimisation schemes such as conjugate gradients may be applied. The simplest example is to take the square of Kohn's penalty functional i.e.

$$P[{\rho}] = \int {\mathrm d}{\bf r}~\left[ {\rho}^2 \left( 1... ...ht]({\bf r},{\bf r}) = \sum_i f_i^2 \left( 1 - f_i \right)^2 ,$$ (6.16)

but other choices are also possible (section 6.2.2). This penalty functional is sketched in figure 6.5.

Figure 6.5: One possible choice of analytic penalty functional.

Since the penalty functional becomes large as the density-matrix becomes less idempotent, minimisation of the total functional $Q[{\rho};\alpha]$ is stable against run-away solutions in which $f_i \rightarrow \infty$ for occupied bands and $f_i \rightarrow -\infty$ for unoccupied bands. However, as illustrated in the sketch in figure 6.6, the minimising density-matrix is only approximately idempotent. In this particular case of an unoccupied band, the total functional is minimised when this band is negatively occupied i.e. for $f_i < 0$. In general the minimising density-matrix ${\bar \rho}$ will have eigenvalues lying outside the interval $[0,1]$, so that the energy calculated from such a density-matrix will be below the true ground-state energy.

Figure 6.6: Schematic illustration of the analytic penalty functional: behaviour of the total energy (black) and total functional (red) with respect to a single occupation number for an unoccupied band.

We denote the set of occupation numbers which minimise the total functional $Q[{\rho};\alpha]$ by $\{ {\bar f}_i \}$ and the errors in these occupation numbers with respect to idempotency by ${\bar f}_i = f_i^{(0)} + \delta f_i$ where $f^{(0)}_i = 1$ for occupied bands and $f^{(0)}_i = 0$ for unoccupied bands. Since $Q[{\rho};\alpha]$ is minimised by ${\bar \rho}$, it is a minimum with respect to all changes in the occupation numbers which maintain the normalisation constraint

$$2 \int {\mathrm d}{\bf r}~{\rho}({\bf r},{\bf r}) = 2 \sum_i f_i = N$$ (6.17)

which is imposed by introducing a Lagrange multiplier $\lambda$:

$$\frac{\partial}{\partial f_i} \left[ Q[{\rho};\alpha] - \lam... ...eft( 2 \sum_j f_j - N \right) \right]_{f_i = {\bar f}_i} = 0 .$$ (6.18)

Using Janak's theorem for the derivative of the energy functional we obtain

$${\bar \varepsilon}_i + \alpha {\bar f}_i \left( 1 - {\bar f}_i \right) \left( 1 - 2 {\bar f}_i \right) - \lambda = 0$$ (6.19)

in which $\{ {\bar \varepsilon}_i \}$ are the eigenvalues of the Hamiltonian obtained from the electronic density ${\bar n}({\bf r}) = 2 {\bar \rho} ({\bf r},{\bf r})$ and are therefore different from the true ground-state eigenvalues. Assuming $\alpha$ to be sufficiently large so that the errors $\{ \delta f_i \}$ are small,

$$\delta f_i \approx -\frac{{\bar \varepsilon}_i - \lambda}{\alpha} .$$ (6.20)

Application of the normalisation constraint 6.17 requires $${\sum_i \delta f_i = 0}$$ from which we obtain the value of the Lagrange multiplier $\lambda$:

$$\lambda = \frac{2}{N} \sum_i {\bar \varepsilon}_i$$ (6.21)

which is the mean energy eigenvalue. The variance of the errors in the occupation numbers is thus related to the variance of the energy eigenvalues, scaled by the parameter $\alpha$. Therefore, as is intuitively expected, the errors in the occupation numbers decrease as $\alpha$ is increased. More precisely, $\delta f_i \propto \alpha^{-1}$ and this behaviour is confirmed numerically in figure 6.7.

Figure 6.7: Variation of the occupation number errors with $\alpha$. Lines show the best fit to $\alpha^{-1}$ behaviour.

For small deviations from idempotency, the penalty functional $${P[{\bar \rho}] \approx \sum_i \delta f_i^2}$$ so that the penalty contribution to the minimised total functional, $\alpha P[{\bar \rho}]$ also decreases in proportion to $\alpha^{-1}$. Hence the energy approaches the true ground-state energy as $\alpha \rightarrow \infty$, again with an error which decreases as $\alpha^{-1}$. This $\alpha^{-1}$ convergence is unsatisfactory for practical applications, since it requires large values of $\alpha$ to obtain accurate estimates of the ground-state energy, and for large values of $\alpha$ the penalty term dominates the total functional and hinders efficient minimisation of the energy term. We now proceed to derive a correction to the estimated energy which allows accurate values for the ground-state energy to be obtained from the approximately idempotent density-matrices which minimise the total functional.

At the minimum of the total functional,

$$\left. \frac{\partial Q[{\rho};\alpha]}{\partial f_i} \right... ... \left( 1 - {\bar f}_i \right) \left( 1 - 2 {\bar f}_i \right)$$ (6.22)

and this expression can be used to construct a first order Taylor expansion for the total energy with respect to the occupation numbers. We thus estimate the true ground-state energy $E_0$ to be

$$E_0 = E[\rho_0] \approx E[{\bar \rho}] + 2 \alpha \sum_i {\b... ...{\bar f}_i \right) \left( 1 - 2 {\bar f}_i \right) \delta f_i.$$ (6.23)

For occupied bands, $\delta f_i = {\bar f}_i - 1$ whereas for unoccupied bands $\delta f_i = {\bar f}_i$ so that

$$E_0 \approx E[{\bar \rho}] - 2 \alpha \sum_i^{\rm all} {\bar... ...left( 1 - {\bar f}_i \right) \left( 1 - 2 {\bar f}_i \right) .$$ (6.24)

The first term of the correction has been written as a sum over all bands so that it can be expressed in terms of the trace $2 \alpha {\rm Tr} \left[ {\bar \rho} \left( 1 - {\bar \rho} \right)^2 \left( 1 - 2 {\bar \rho} \right) \right]$ which can always be evaluated in ${\cal O}(N)$ operations. The second term only contributes when unoccupied bands are included in the calculation, which is not necessary for insulators. Since a single eigenvalue of the (sparse) density-matrix can always be obtained in ${\cal O}(N)$ operations, it is possible to evaluate the correction for a small number ($\ll N$) of unoccupied bands and retain the linear-scaling. However, we note that the correction need only be calculated once the minimum of the total functional has been found, so that a single ${\cal O}(N^2)$ step to obtain all of the occupation numbers will still be a tiny fraction of the total computational effort.

The error in a Taylor expansion is generally estimated by considering the lowest order term omitted, which in this case is

$$\frac{1}{2} \sum_{ij} \delta f_i \left. \frac{\partial^2 E[{... ...elta f_j \approx \sum_{ij} \delta f_i {\cal H}_{ij} \delta f_j$$ (6.25)

where ${\cal H}_{ij}$ is the chemical hardness matrix. Unfortunately this matrix is not guaranteed to be either positive or negative definite, and so the estimate of the ground-state energy 6.24 is not a strict upper or lower bound to the true energy. As will be seen shortly, this error is generally much smaller than other sources of error (such as the finite support region size, with respect to which the energy does behave variationally) so that this is not an issue in practice.

In figure 6.8 the energy, total functional and corrected energy are plotted for different values of the parameter $\alpha$. These results confirm the $\alpha^{-1}$ behaviour of the energy and penalty functionals, and the error in the corrected energy even for $\alpha = 50$ eV is smaller than $10^{-4}$ eV per atom. Thus we have achieved our aim of being able to obtain accurate estimates of the ground-state energy from approximate density-matrices which minimise the total functional for values of $\alpha$ for which efficient minimisation is possible.

Figure 6.8: Total energy, total functional and corrected energy versus $\alpha$. Lines are best fits to $\alpha^{-1}$ behaviour.