5.5 Kinetic energy matrix elements

The kinetic energy matrix elements for any two basis functions $\chi_{\alpha , n \ell m}$ and $\chi_{\beta , n' \ell' m'}$ centred at ${\bf R}_{\alpha}$ and ${\bf R}_{\beta}$ respectively are defined by

$${\cal T}_{\alpha, n \ell m ; \beta, n' \ell' m'} = -{\textstyle{1 \over 2}} \int {\mathrm d}{\bf r}~\chi_{\alpha , n \ell m}({\bf r}) \nabla^2 \chi_{\beta , n' \ell' m'}({\bf r})$$

$$= {1 \over 2(2 \pi)^3} \int {\mathrm d}{\bf k}~k^2 \exp[-{\mathrm{i... ...i}_{\alpha , n \ell m}({\bf k}) {\tilde \chi}_{\beta , n' \ell' m'}(-{\bf k}) .$$ (5.26)

Because of the discontinuity in the first derivatives of the basis functions at the sphere boundaries, a delta-function arises when the Laplacian operates on a basis function. This is integrated out when the matrix element is calculated and this contribution is included when transforming the real-space integral to reciprocal-space in equation 5.26.

The second line of equation 5.26 is identical to equation 5.11 apart from a factor of ${\textstyle{1 \over 2}} k^2$. The same separation into individually regular terms can be applied here, and the result is that we need to calculate the contour integral (5.17) as before, except that the integer $p$ must be replaced by $(p-2)$ and a numerical factor of ${1 \over 2}$ is introduced. The calculation of the residues is identical to that presented in the previous section, except that the integrand no longer always has a pole at $z = 0$ in every term.

The results for ${\cal T}_{\alpha, n \ell m ; \beta, n' \ell' m'}$ when ${\bf R}_{\alpha \beta} = 0$ are

$$\frac{{\textstyle{1 \over 2}}\delta_{\ell \ell'}\delta_{m m'}}{q_... ... r_{\alpha} \geq r_{\beta} \end{array} \right\} q_{n \ell} \not= q_{n' \ell'} ,$$ (5.27)

$$\textstyle{1 \over 4} \delta_{\ell \ell'} \delta_{m m'} q_{n \ell... ...quad r_{\alpha} \geq r_{\beta} \end{array} \right\} q_{n \ell} = q_{n' \ell'} .$$

The calculation of the kinetic energy has been checked by projecting a set of wave-functions expanded in the spherical-wave basis onto the plane-wave basis using equation 5.9a. As the kinetic energy cut-off for the plane-wave basis is increased, so the description of the wave-functions becomes more accurate. The kinetic energy calculated using the results above can then be compared against the kinetic energy calculated by the plane-wave ${\cal O}(N^3)$ CASTEP code [82].

From the asymptotic behaviour of the spherical Bessel functions, the Fourier transform (5.9a) for large $k$ is

$${\tilde \chi}_{\alpha , n \ell m}({\bf k}) \propto \frac{\si... ...- {\ell \pi \over 2})}{k^3}~{\bar Y}_{\ell m} (\Omega_{\bf k})$$ (5.28)

and so the error in the kinetic energy due to truncating the plane-wave basis with cut-off $E_{\mathrm{cut}} = {\textstyle{1 \over 2}} k_{\mathrm{cut}}^2$ is

$$\Delta T \propto \int_{k_{\mathrm{cut}}}^{\infty} {\mathrm d... ...{k_{\mathrm{cut}}} \propto \frac{1}{\sqrt{E_{\mathrm{cut}}}} .$$ (5.29)

In figure 5.1 the kinetic energy as calculated by the plane-wave code has been plotted against $1 / \sqrt{E_{\mathrm{cut}}}$ and yields a straight line as expected, which can then be extrapolated to obtain an estimate of the kinetic energy calculated for infinite cut-off: $60.66 \pm 0.01$ eV. This is in agreement with the value calculated analytically of $60.65$ eV.

Figure 5.1: Plot of asymptotic fit to kinetic energy data.