5.6 Non-local pseudopotential

In this section we present two methods for obtaining the non-local pseudopotential matrix elements. The first constructs a Green's function in order to expand the basis functions on one site in terms of the same functions centred elsewhere. This method is particularly efficient, although it requires some restrictions to be made in the radial part of the support functions to give a variational method. The second method adopts the Kleinman-Bylander form and uses the results for the overlap matrix elements. This method is slower, but more robust numerically and allows direct comparison with existing pseudopotential codes.

5.6.1 Green's function method

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 \delta {\hat V}_{\ell} \langle \ell m \vert$$ (5.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 $\delta 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. We can write the basis function over all space using the Heaviside step function $H(x)$

$$H(x) = \left\{ \begin{array}{ll} 0, & \qquad x < 0 ,\\ 1, & \qquad x \geq 0 , \end{array}\right.$$ (5.31)

so that a basis function centred in a sphere of radius $r_{\alpha}$ at the origin (i.e. for ${\bf R}_{\alpha}=0$) is

$$\chi_{\alpha , n \ell m}({\bf r}) = j_{\ell}(q_{n \ell}r) {\bar Y}_{\ell m}( \Omega_{\bf r}) H(r_{\alpha} - r)$$ (5.32)

and therefore

$$(\nabla^2 + q_{n \ell}^2) \chi_{\alpha , n \ell m}({\bf r}) ... ...elta(r_{\alpha}-r) \right] {\bar Y}_{\ell m}(\Omega_{\bf r}) .$$ (5.33)

The terms on the right-hand side only contribute on the sphere boundary $r = r_{\alpha}$; everywhere else in space the basis function obeys the homogeneous Helmholtz equation (5.1). These terms will give rise to an extra term in the Green's function solution due to the discontinuity of the first radial derivative of the basis function at the sphere boundary. However, if the radial function for each angular momentum component has a continuous first derivative, then these contributions will cancel. In this case, we can proceed assuming that the basis function obeys the homogeneous equation everywhere. This condition is naturally obeyed by the support functions when the support regions are large enough, since there is a kinetic energy penalty associated with the same discontinuity in the radial wave-function, but in this case we would lose any variational principle if there was a discontinuity at any stage during the minimisation. A better solution is to impose a sum rule on the expansion coefficients in equation 5.3 to maintain a continuous radial derivative i.e.

$$\sum_n c^{n \ell m}_{(\alpha)} q_{n \ell} j_{\ell}'(q_{n \ell} r_{\alpha}) = 0$$ (5.34)

which can be used to fix one coefficient in terms of the rest for each angular momentum component. This is the restriction mentioned above which must be imposed if the Green's function method is to be used. The basis function is assumed to obey the homogeneous Helmholtz equation throughout all space. In the original support region, homogeneous boundary conditions were applied. Now, in an overlapping region of radius $r_{\mathrm c}$ (the core region for the non-local pseudopotential), the basis function must still obey the same homogeneous Helmholtz equation, but it is now subject to inhomogeneous boundary conditions. Standard methods can be used to transform this problem into the solution of an inhomogeneous Helmholtz equation subject to homogeneous boundary conditions, and this new problem has a standard solution in terms of the Green's function, which can be formally expanded in terms of the eigenfunctions of an appropriate self-adjoint operator (here the operator on the left of the Helmholtz equation). The result is

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

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}_{... ...chi_{\alpha , n \ell m}({\bf r'} + {\bf R}_{\mathrm{ion}} - {\bf R}_{\alpha}) ,$$ (5.36)

$$a_{n' \ell'}^{n \ell m} = \frac{2}{r_{\mathrm{c}}^3~j_{\ell' - 1}^2 \left( \lambda_{n' \ell... ...m{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 rôle as the $\{ q_{n \ell} \}$ in the expansion of the wave-functions. The integral in equation 5.37 is straightforward to evaluate for given ${\ell'}$.

The surface integral in equation 5.36 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 [151]. 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 terms of the components $(x,y,z)$ of the vector ${\bf R}_{\alpha} - {\bf R}_{\mathrm{ion}}$ of length $r = \vert{\bf R}_{\alpha} - {\bf R}_{\mathrm{ion}}\vert$, and $v^2 = r(z+r)$, the rotation matrix $O$ is given by:

$$O = \frac{1}{v^2} \left( \begin{array}{ccc} v^2 - x^2 & -x y... ...& -y (z+r) \\ x (z+r) & y (z+r) & z (z+r) \end{array} \right)$$ (5.37)

The singularity which occurs when $z = -r$ i.e. when the ion is ``vertically'' above the support region, is avoided by inverting the calculation through the origin when $z < 0$. The spherical harmonic matrices $C_{m m'}^{(\ell)}$ can be written in terms of the elements of this matrix. Here we write them out explicitly for $\ell = 0,1,2$:

$$C^{(0)} = 1$$ (5.38)

$$C^{(1)} = \left( \begin{array}{ccc} O_{22} & O_{32} & O_{12} \\ O_{23} & O_{33} & O_{13} \\ O_{21} & O_{31} & O_{11} \end{array} \right)$$ (5.39)

$$C^{(2)} = \left( \begin{array}{cccc} O_{12} O_{21} + O_{11} O_{22} & O_{22}... ...\textstyle{\sqrt{3} \over 2} (O_{31}^2 - O_{32}^2) & \ldots \end{array} \right.$$

$$\qquad \left. \begin{array}{ccc} \ldots & O_{12} O_{31} + O_{11} ... ...tyle{1 \over 2} (O_{11}^2 + O_{22}^2 - O_{12}^2 - O_{21}^2) \end{array} \right)$$ (5.40)

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.

Figure 5.2: Non-local pseudopotential energy against number of Bessel functions used in Green's function expansion

The numerical evaluation of the analytic results for $K_{\ell' n \ell m}(u,q_{n \ell})$ is inaccurate when $u \ll 1$ and in this case it is necessary to employ Taylor expansions of the results. The two different cases for the upper limit of the integral also need to be treated separately, and this is one reason why the Kleinman-Bylander method described next is preferred.

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

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

Defining the core matrix elements

$$\delta V_{n n'}^{\ell} = \left\{ \begin{array}{ll} \displays... ...rm d}r~r^2~\delta V_{\ell}(r) }& n = n' = 0 \end{array}\right.$$ (5.42)

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

$${\cal V}_{\alpha , n \ell m ; \beta , n' \ell' m'} = \sum_{L M} f_{L M}^{n \ell m} f_{L M}^{n' \ell' m'} \left[ \sum_{i,j=0} a_{i L}^{n \ell m} \delta V_{i j}^L a_{j L}^{n' \ell' m'} \right] .$$ (5.43)

The non-local pseudopotential data is therefore stored in terms of the core matrix elements defined in equation 5.45. In figure 5.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 [74]. 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.

5.6.2 Kleinman-Bylander form

Figure 5.3: Non-local pseudopotential energy against number of Bessel functions used to describe Kleinman-Bylander projectors.

We reproduce equation 3.76 which demonstrates the Kleinman-Bylander form of the pseudopotential in terms of pseudo-atomic eigenstates $\{ \vert \phi_{\ell m} \rangle \}$.

$${\hat V}_{\mathrm{KB}} = {\hat V}_{\mathrm{loc}} + \sum_{\el... ...l m} \vert \delta {\hat V}_{\ell} \vert \phi_{\ell m} \rangle}$$ (5.44)

For pseudopotentials whose non-local components $\{ \delta V_{\ell} \}$ vanish at the core radius $r_{\mathrm c}$, the projector states $\{ \vert \delta {\hat V}_{\ell} \phi_{\ell m} \rangle \}$ can be expanded in the basis-states $\{ \vert \chi_{n \ell m} \rangle \}$, and then the matrix-elements of the pseudopotential operator are straightforwardly obtained by applying the result for the overlap matrix elements, without resorting to the scheme by King-Smith et al. [152].

In figure 5.3 we show how the non-local pseudopotential energy rapidly converges as the number of Bessel functions used to expand the Kleinman-Bylander projectors is increased. This method is numerically more stable than the Green's function method, although slower, and has the added advantage that it allows a direct comparison to be made between the results calculated using this basis-set and traditional plane-wave codes. For these reasons the Kleinman-Bylander method has been used in computational implementations.