6. Non-local pseudopotential

The general form for a semi-local pseudopotential operator (i.e. one which is non-local in the angular but not radial coordinates) for an ion is

$${\hat V}_{\mathrm{NL}} = \sum_{\ell m} \vert \ell m \rangle {\hat V}_{\ell} \langle \ell m \vert$$ (30)

where $\langle {\bf r} \vert \ell m \rangle = {\bar Y}_{\ell m}(\Omega)$ and ${\bar Y}_{\ell m}$ is centred on the ion.

The pseudopotential components $V_{\ell}$ are themselves short-ranged in real-space, and vanish beyond the core radius $r_{\mathrm{c}}$. Therefore the action of the non-local pseudopotential depends only upon the form of the wave functions within this core region. We require the matrix elements of the non-local pseudopotential between localised basis functions which are not necessarily centred on the ion.

We therefore need to find an expansion of the basis functions in terms of functions localised within the pseudopotential core. Since the basis functions are all solutions of the Helmholtz equation, we invoke the uniqueness theorem which states that the expansion we seek is uniquely determined by the boundary conditions on the surface of the core region and solve the Helmholtz equation subject to these inhomogeneous boundary conditions by the standard method using the formal expansion of the Green's function. The result is

$$\chi_{n \ell m}^{\alpha}({\bf r}) = \sum_{\ell' m'} f_{\ell'... ...(q_{n' \ell'} r') \right] {\bar Y}_{\ell' m'}(\Omega_{\bf r'})$$ (31)

and is valid for points ${\bf r'} = {\bf r} - {\bf R}_{\mathrm{ion}} + {\bf R}_{\alpha}$ within the core region (i.e. for $r' \leq r_{\mathrm{c}}$.)

The coefficients $f_{\ell' m'}^{n \ell m}$ and $a_{n' \ell'}^{n \ell m}$ are defined by:

$$f_{\ell' m'}^{n \ell m} = \int_{r' = r_{\mathrm{c}}} \mathrm{d}\Omega_{\bf r'} {\bar Y}_{\e... ...chi_{n \ell m}^{\alpha}({\bf r'} + {\bf R}_{\mathrm{ion}} - {\bf R}_{\alpha}) ,$$ (32)

$$a_{n' \ell'}^{n \ell m} = \frac{2}{r_{\mathrm{c}}^3~j_{\ell' - 1}^2 \left( \lambda_{n' \ell... ...rm{c}}} \mathrm{d}r~r^2 j_{\ell'} \left( \lambda_{n' \ell'} r \right) \right] .$$

The $\{ \lambda_{n \ell} \}$ are chosen by $j_{\ell} ( \lambda_{n \ell} r_{\mathrm{c}} ) = 0$ and play the same role as the $\{ q_{n \ell} \}$ in the expansion of the wave functions. The integral in equation (33) is straightforward to evaluate for given ${\ell'}$.

The surface integral in equation (32) is evaluated by first rotating the coordinate system so that the new $z$-axis is parallel to ${\bf R}_{\alpha} - {\bf R}_{\mathrm{ion}}$, thus mixing the spherical harmonics [5]. The elements of the orthogonal spherical harmonic mixing matrices $C_{m m'}^{\ell}$ are defined by the elements of the rotation matrix for the coordinate system. In the new coordinate system, the surface integral is written in terms of a one-dimensional integral

$$K_{\ell' n \ell m}(u,q_{n \ell}) = \frac{1}{2 u} \left[ (2 \ell' ... ...t)!}{(\ell' + \vert m\vert)! (\ell + \vert m\vert)!} \right]^{1 \over 2} \times$$

$$\quad \int_{\vert u-1\vert}^{{\mathrm{min}}(u+1,r_{\alpha}/r_{\ma... ...athrm{c}} z) P_{\ell}^{\vert m\vert} \left( \frac{1 - u^2 - z^2}{2 u z} \right)$$

in which the dimensionless variable $u = {\vert{\bf R}_{\alpha} - {\bf R}_{\mathrm{ion}}\vert \over r_{\mathrm{c}}}$ is introduced. $P_{\ell}^{\vert m\vert} (x)$ denotes an associated Legendre polynomial, and these integrals can all be calculated indefinitely using elementary methods once the integrand is expanded into trigonometric functions.

The final result for $f_{\ell' m'}^{n \ell m}$ is then

$$f_{\ell' m'}^{n \ell m} = \sum_{m'' = -{\mathrm{min}(\ell, \... ...l')} C_{m'' m'}^{\ell'} K_{\ell' n \ell m} C_{m'' m} ^{\ell} .$$ (33)

Defining the core matrix elements

$$V_{n n'}^{\ell} = \left\{ \begin{array}{ll} \displaystyle{ \... ...athrm{d}r~r^2~V_{\ell}(r) }& n = n' = 0 \\ \end{array}\right.$$ (34)

the matrix element of the non-local pseudopotential operator between any two basis functions overlapping the core ( $\chi_{n \ell m}^{\alpha}$ and $\chi_{n' \ell' m'}^{\beta}$) can be written as the sum:

$$V_{\mathrm{NL}, \alpha \beta} = \sum_{\ell'' m''} f_{\ell'' m''}^{n \ell m} f_{\ell'' m''}^{n' \e... ...'}^{n \ell m} + a_{n'' \ell''}^{n' \ell' m'} \right) V_{n'' 0}^{\ell''} \right.$$

$$+ \left. \sum_{n'' n'''} a_{n'' \ell''}^{n \ell m} a_{n''' \ell''}^{n' \ell' m'} V_{n'' n'''}^{\ell''} \right] .$$ (35)

The non-local pseudopotential data is therefore stored in terms of the core matrix elements defined in equation (37). In figure 2 we plot the non-local pseudopotential energy against the number of core Bessel functions for an s-local silicon pseudopotential generated according to the scheme of Troullier and Martins [3]. We see that the energy converges rapidly with the number of core Bessel functions used (the dashed line is the energy calculated with fifty core functions.) Increasing the number of core functions only increases the number of $a_{n' \ell'}^{n \ell m}$ coefficients required, and the separable nature of the calculation means that even using fifty core functions requires very little computational effort.

Figure 2: Non-local pseudopotential energy against number of core Bessel functions.