2. Origin of the basis functions

In the pseudopotential approximation, the core electrons and strong ionic potential of the atom are replaced by a much weaker potential in which the remaining pseudovalence electrons move. The pseudovalence states no longer have to be orthogonal to lower-lying core states and hence are much smoother than the all-electron valence states in the core region and have less kinetic energy. Thus the pseudovalence states can be accurately represented by a much smaller set of plane-wave basis functions than the all-electron states.

The plane-wave basis state $\mathrm{e}^{{\mathrm{i}} {\bf q} \cdot {\bf r}}$ is a solution of the Helmholtz equation (the time-independent free-electron Schrödinger equation)

$$\left( \nabla^2 + q^2 \right) \psi ({\bf r}) = 0$$ (1)

subject to periodic boundary conditions, with energy $E = {\textstyle{1 \over 2}}q^2$ (we use atomic units throughout.)

If instead we wish to localise the basis functions, say within spherical regions of radius $a$, so that the function vanishes outside these regions, then appropriate conditions would be to require the functions to be finite within the regions and to vanish on the boundary. The solutions to the Helmholtz equation (1) subject to these conditions are then truncated spherical-waves

$$\psi ({\bf r}) = \left\{ \begin{array}{ll} j_{\ell}(q r)~Y_{... ...\varphi), \qquad & r < a \\ 0, & r \geq a \end{array} \right.$$ (2)

where $(r,\vartheta,\varphi)$ are spherical polar coordinates with the origin at the centre of the spherical region, ${\ell}$ is a non-negative integer, $m$ is an integer satisfying $-{\ell} \leq m \leq {\ell}$ and $q$ is chosen to satisfy $j_{\ell}(q a) = 0$. $j_{\ell}(x)$ is a spherical Bessel function and $Y_{\ell m}(\Omega)$ is a spherical harmonic. Solutions involving the spherical von Neumann function $n_{\ell}(x)$ have been rejected because they are not finite at the centre of the sphere.

We note that these functions solve the same equation as the plane-wave basis functions, so that within the pseudopotential approximation the wave functions will be well-described by a truncated set of these basis functions. Moreover, these functions are eigenstates of the kinetic energy operator within the localisation region $r < a$ (i.e. in the region in which they will be used to describe the wave functions) with eigenvalue ${\textstyle{1 \over 2}}q^2$ so that the same kinetic energy cut-off used to truncate the plane-wave basis can be used here to restrict the values of ${\ell}$ and $q$.

Since the Laplacian is a self-adjoint operator under these boundary conditions, application of Sturm-Liouville theory proves that all states within the same spherical region are mutually orthogonal.

In a calculation, the electronic states are described by covering the simulation cell with overlapping spheres, usually chosen to be centred on the ions at positions ${\bf R}_{\alpha}$, and expanding the wave functions ${\phi}_{\alpha}$ within these spheres in this basis:

$${\phi}_{\alpha}({\bf r}) = \sum_{n \ell m} c_{n \ell m}^{\al... ...\right\vert)~Y_{\ell m}(\Omega_{{\bf r} - {\bf R}_{\alpha}}) .$$ (3)

The notation $\Omega_{\bf r}$ is introduced as shorthand for the polar and azimuthal angles of the vector ${\bf r}$ used to represent that vector in spherical polar coordinates. We denote the radius of the sphere by $r_{\alpha}$ so that the $\left\{ q_{n \ell} \right\}$ are defined by $j_{\ell} (q_{n \ell} r_{\alpha}) = 0$.

The expansion (3) is frequently written down formally, but rarely used computationally because of the inconvenience of using spherical Bessel functions in numerical work. However, the analytic results derived in the following sections offset this disadvantage.

$O(N)$ methods are aimed at large systems, and so the Brillouin zone sampling of the electronic states is usually restricted to the states at the $\Gamma$-point only. The wave functions can then be made real without loss of generality, and so in practice we use real linear combinations of the spherical harmonics defined below, which does not alter any of the analysis here.

$$\{Y_{\ell m}\} \rightarrow \{{\bar Y}_{\ell m}\} = \left\{ \... ...ega) - (-1)^m Y_{\ell,m} (\Omega) \right] \end{array} \right\}$$ (4)

These real combinations of spherical harmonics, which we denote $\bar{Y}_{\ell m}$, can be written down as real functions of the variables $\left\{{x \over r},{y \over r},{z \over r}\right\}$ where $(x,y,z)$ are Cartesian coordinates with origin at the centre of the sphere, and are familiar as the angular components of s, p, d etc. orbitals.

We introduce $\chi_{n \ell m}^{\alpha}({\bf r})$ to represent a truncated spherical-wave basis function centred at the origin and confined to a sphere of radius $r_{\alpha}$:

$$\chi_{n \ell m}^{\alpha}({\bf r}) = \left\{ \begin{array}{ll... ... r \leq r_{\alpha}, \\ 0, & r > r_{\alpha}. \end{array}\right.$$ (5)

Equation (3) can then be written:

$${\phi}_{\alpha}({\bf r}) = \sum_{n \ell m} c_{n \ell m}^{\alpha} ~\chi_{n \ell m}^{\alpha}({\bf r} - {\bf R}_{\alpha}) .$$ (6)