7.7 Practical details

In this section we outline a number of details concerning the implementation of the penalty method. These concern the derivatives of the functional with respect to the expansion coefficients for the support functions $\{ c^{n \ell m}_{(\alpha)} \}$, the imposition of the normalisation constraint and the general outline of the scheme.

7.7.1 Expansion coefficient derivatives

The support functions are expanded in spherical-wave basis functions:

$${\phi}_{\alpha}({\bf r}) = \sum_{n \ell m} c^{n \ell m}_{(\alpha)} ~\chi_{\alpha , n \ell m}({\bf r}) .$$ (7.71)

For a functional of the support functions $f[\{ \phi_{\alpha} \}]$, the derivative with respect to the expansion coefficients is

$$\frac{\partial f[\{ \phi_{\alpha} \}]}{\partial c^{n \ell m}_{(\beta)}} = \sum_{\gamma} \int {\mathrm d}{\bf r}~ \frac{\delta f[ \{ \phi_{\... ...bf r})} \frac{\partial \phi_{\gamma}({\bf r})}{\partial c^{n \ell m}_{(\beta)}}$$

$$= \sum_{\gamma} \int {\mathrm d}{\bf r}~ \frac{\delta f[ \{ \phi_{\... ...hi_{\gamma}({\bf r})} \delta_{\gamma}^{\beta} \chi_{\gamma , n \ell m}({\bf r})$$

$$= \int {\mathrm d}{\bf r}~\frac{\delta f[ \{ \phi_{\alpha} \} ]} {\delta \phi_{\beta}({\bf r})} \chi_{\beta , n \ell m}({\bf r}) .$$ (7.72)

For example, for the derivative of the total energy,

$$\frac{\delta E[\rho]}{\delta \phi_{\alpha}({\bf r})} = 4 K^{\alpha \beta} {\hat H} \phi_{\beta}({\bf r}) ,$$ (7.73)

we obtain

$$\frac{\partial E[\rho]}{\partial c^{n \ell m}_{(\beta)}} = 4 K^{\beta \gamma} \int {\mathrm d}{\bf r}~ \chi_{\beta , n \ell m}({\bf r}) {\hat H} \phi_{\gamma}({\bf r})$$

$$= 4 K^{\beta \gamma} \sum_{n' \ell' m'} \langle \chi_{\beta , n \el... ...t {\hat H} \vert \chi_{\gamma , n' \ell' m'} \rangle c^{n' \ell' m'}_{(\gamma)}$$ (7.74)

in which $\langle \chi_{\beta , n \ell m} \vert {\hat H} \vert \chi_{\gamma , n' \ell' m'} \rangle = {\cal H}_{\beta , n \ell m; \gamma , n' \ell' m}$ denotes the matrix element of the Kohn-Sham Hamiltonian with respect to the spherical-wave basis functions.

7.7.2 Normalisation constraint

We choose to minimise the total functional whilst constraining the normalisation of the density-matrix to the correct value. This is achieved firstly by projecting all gradients to be perpendicular to the gradient of the electron number, thus maintaining the electron number to first order, and secondly by re-converging the electron number to its correct value before each evaluation of the total functional.

7.7.2.1 Density-kernel variation

In section 7.3 the electron number gradient with respect to the density-kernel, which is here denoted ${\mit\Delta}$, is given as

$$\frac{\partial N}{\partial K^{\alpha \beta}}= 2 S_{\beta \alpha}={\mit\Delta}_{\alpha \beta} .$$ (7.75)

The density-kernel along this search direction ${\mit\Delta}$ is parameterised by $\lambda$ as

$$K(\lambda) = K(0) + \lambda {\mit\Delta}$$ (7.76)

where $K(0)$ denotes the initial density-matrix. The electron number, given by $N = 2 {\rm Tr}(KS)$ thus behaves linearly:

$$N(\lambda) = N(0) + 2 \lambda {\rm Tr}({\mit\Delta}S)$$

$$= N(0) + 4 \lambda {\rm Tr}(S^2)$$ (7.77)

and it is a trivial matter to calculate the required value of $\lambda$ to return the electron number to its correct value.

In general during the minimisation, the search direction is ${\mit\Lambda}$, and again this search can be parameterised by a single parameter $\lambda$:

$$K(\lambda) = K(0) + \lambda {\mit\Lambda} .$$ (7.78)

We wish to project out from ${\mit\Lambda}$ that component which is parallel to ${\mit\Delta}$. The modified search direction ${\tilde {\mit\Lambda}}$ can be written as

$${\tilde {\mit\Lambda}} = {\mit\Lambda} - \omega {\mit\Delta} .$$ (7.79)

The variation of the electron number along this modified direction is

$$N(\lambda) = N(0) + 2 \lambda {\rm Tr}({\mit\Lambda}S) - 2 \omega \lambda {\rm Tr}({\mit\Delta}S)$$ (7.80)

and we wish the coefficient of the linear term in $\lambda$ to vanish, which defines the required value of $\omega$ to be

$$\omega = \frac{{\rm Tr}({\mit\Lambda}S)}{{\rm Tr}({\mit\Delta}S)} = \frac{{\rm Tr}({\mit\Lambda}S)}{2{\rm Tr}(S^2)} .$$ (7.81)

Since the electron number depends linearly upon the density-kernel, after this projection, the electron number is constant along the modified search direction, and the electron number need not be corrected after a trial step is taken.

7.7.2.2 Support function variation

The electron number gradient with respect to the support functions is given in section 7.3 and denoted by $\{ \zeta^{\alpha}({\bf r})\}$:

$$\frac{\delta N}{\delta \phi_{\alpha}({\bf r})}= 2 K^{\alpha \beta} \phi_{\beta}({\bf r}) = \zeta^{\alpha}({\bf r}) .$$ (7.82)

The support function variation is again parameterised by the parameter $\lambda$:

$$\phi_{\alpha}({\bf r};\lambda) = \phi_{\alpha}({\bf r};\lambda=0) + \lambda \zeta^{\alpha}({\bf r})$$ (7.83)

and this results in the following quadratic variation of the overlap matrix

$$S_{\alpha \beta}(\lambda) = S_{\alpha \beta}(0) + \lambda \langle \phi_{\alpha}(\lambda=0) \v... ...ambda=0) \rangle + \lambda^2 \langle \zeta^{\alpha} \vert \zeta^{\beta} \rangle$$

$$= S_{\alpha \beta}(0) + \lambda S'_{\alpha \beta} +\lambda^2 S''_{\alpha \beta}$$ (7.84)

which defines the matrices $S'$ and $S''$. The variation of the electron number is therefore also quadratic

$$N(\lambda) = N(0) + 2 \lambda {\rm Tr}(KS') + 2 \lambda^2 {\rm Tr}(KS'')$$ (7.85)

and the roots of this expression can be found to correct the electron number.

We consider a general search direction for the localised functions denoted $\{ \xi^{\alpha}({\bf r}) \}$ and modify this search direction to obtain the direction which maintains the electron number to first order:

$${\tilde \xi}^{\alpha}({\bf r}) = \xi^{\alpha}({\bf r}) - \omega \zeta^{\alpha}({\bf r}) .$$ (7.86)

Now varying the localised functions according to

$$\phi_{\alpha}({\bf r};\lambda) = \phi_{\alpha}({\bf r};\lambda=0) + \lambda {\tilde \xi}^{\alpha}({\bf r})$$ (7.87)

results in the following variation of the electron number:

$$N(\lambda) = N(0) + 2 \lambda {\rm Tr}(K{\tilde S}') - 2 \lambda \omega {\rm Tr}(KS') + {\cal O}(\lambda^2)$$ (7.88)

where ${\tilde S}'_{\alpha \beta} = \langle \phi_{\alpha}(\lambda=0) \vert \xi^{\beta} \rangle + \lambda \langle \xi^{\alpha} \vert \phi_{\beta}(\lambda=0) \rangle$. Thus to maintain the electron number to first order, we choose

$$\omega = \frac{{\rm Tr}(K{\tilde S}')}{{\rm Tr}(KS')} .$$ (7.89)

In this case, the electron number is not constant along the search direction, but will still vary quadratically, so that it is necessary to correct the density-matrix before evaluating the total functional.

7.7.3 General outline of the scheme

The initialisation of the density-matrix is described in chapter 8. The functional minimisation consists of two nested loops. The inner loop consists of the minimisation with respect to the density-kernel, while keeping the support functions constant. As shown in section 7.4, this corresponds to making the density-matrix commute with the Hamiltonian in the representation of the current support functions. The outer loop consists of the support function minimisations, which in section 7.4 were shown to correspond to solving the Kohn-Sham equations. In general, we find that two or three cycles of the inner loop for each cycle of the outer loop suffice, and this is demonstrated in figure 7.1 in which three cycles of the inner loop appears to give the best performance to computational cost ratio.

Figure 7.1: Rate of convergence for different numbers of inner cycles. In the legend, ${\rm DM}_n$ corresponds to $n$ cycles of the inner loop for each iteration of the outer loop (horizontal axis).

The conjugate gradients scheme is used to determine the search directions from the gradients, and these gradients are then projected perpendicular to the electron number gradient as described in section 7.7.2. We approximate the total functional by a parabola along the search direction, using the initial value of the functional, the first derivative of the functional at the initial position (which is simply the scalar product of the steepest descent and search directions) and the value of the functional at some trial position. The value of the functional at the predicted minimum is then evaluated. If this value deviates significantly from the value predicted by the quadratic fit, a cubic fit is constructed using this new value of the functional. In general, this is only necessary for the first few steps, and the parabolic fit is very good. In the case of the support function variation, the support functions are altered to give the correct number of electrons before each evaluation of the functional. The conjugate gradients procedure is reset after a finite number of steps, and this is illustrated in figure 7.2 for the inner and outer loops. In both cases, we see that there is little advantage in conjugating more than eight gradients before resetting.

Figure 7.2: Performance of the conjugate gradients algorithm in the density-kernel variation (left) and the support function variation (right).

Once the functional has been minimised, the correction to the total energy is calculated. The whole calculation is generally repeated for a few different values of $\alpha$ to ensure that the corrected energy has indeed converged.