3.3 The pseudopotential approximation

In this section we outline a further approximation which is based upon the observation that the core electrons are relatively unaffected by the chemical environment of an atom. Thus we assume that their (large) contribution to the total binding energy does not change when isolated atoms are brought together to form a molecule or crystal. The actual energy differences of interest are the changes in valence electron energies, and so if the binding energy of the core electrons can be subtracted out, the valence electron energy change will be a much larger fraction of the total binding energy, and hence much easier to calculate accurately. We also note that the strong nuclear Coulomb potential and highly localised core electron wave-functions are difficult to represent computationally.

Since the atomic wave-functions are eigenstates of the atomic Hamiltonian, they must all be mutually orthogonal. Since the core states are localised in the vicinity of the nucleus, the valence states must oscillate rapidly in this core region in order to maintain this orthogonality with the core electrons. This rapid oscillation results in a large kinetic energy for the valence electrons in the core region, which roughly cancels the large potential energy due to the strong Coulomb potential. Thus the valence electrons are much more weakly bound than the core electrons.

It is therefore convenient to attempt to replace the strong Coulomb potential and core electrons by an effective pseudopotential which is much weaker, and replace the valence electron wave-functions, which oscillate rapidly in the core region, by pseudo-wave-functions, which vary smoothly in the core region [56,57]. We outline two justifications for this approximation below; for further details see [58] and also [59,60] for recent reviews.

Figure 3.2: Schematic diagram of the relationship between all-electron and pseudo- potentials and wave-functions.

3.3.1 Operator approach

Following the orthogonalised plane-waves approach [61], we consider an atom with Hamiltonian ${\hat H}$, core states $\{ \vert \chi_n \rangle \}$ and core energy eigenvalues $\{ E_n \}$ and focus on one valence state $\vert \psi \rangle$ with energy eigenvalue $E$. From these states we attempt to construct a smoother pseudo-state $\vert \varphi \rangle$ defined by

$$\vert \psi \rangle = \vert \varphi \rangle + \sum_n^{\mathrm{core}} a_n \vert \chi_n \rangle .$$ (3.55)

The valence state must be orthogonal to all of the core states (which are of course mutually orthogonal) so that

$$\langle \chi_m \vert \psi \rangle = 0 = \langle \chi_m \vert \varphi \rangle + a_m$$ (3.56)

which fixes the expansion coefficients $\{ a_n \}$. Thus

$$\vert \psi \rangle = \vert \varphi \rangle - \sum_n^{\mathrm... ...}} \vert \chi_n \rangle \langle \chi_n \vert \varphi \rangle .$$ (3.57)

Substituting this expression in the Schrödinger equation ${\hat H} \vert \psi \rangle = E \vert \psi \rangle$ gives

$${\hat H} \vert \varphi \rangle - \sum_n^{\mathrm{core}} E_n ... ...um_n \vert \chi_n \rangle \langle \chi_n \vert \varphi \rangle$$ (3.58)

which can be rearranged in the form

$${\hat H} \vert \varphi \rangle + \sum_n^{\mathrm{core}} (E -... ...\langle \chi_n \vert \varphi \rangle = E \vert \varphi \rangle$$ (3.59)

so that the smooth pseudo-state obeys a Schrödinger equation with an extra energy-dependent non-local potential ${\hat V}_{\mathrm{nl}}$;

$$\left[ {\hat H} + {\hat V}_{\mathrm{nl}} \right] \vert \varphi \rangle = E \vert \varphi \rangle$$ (3.60)

$${\hat V}_{\mathrm{nl}} = \sum_n^{\mathrm{core}} (E - E_n) \vert \chi_n \rangle \langle \chi_n \vert .$$ (3.61)

The energy of the smooth state described by the pseudo-wave-function is the same as that of the original valence state. The additional potential $V_{\mathrm{nl}}$, whose effect is localised in the core, is repulsive and will cancel part of the strong Coulomb potential so that the resulting sum is a weaker pseudopotential. Of course, once the atom interacts with others, the energies of the eigenstates will change, but if the core states are reasonably far from the valence states in energy (i.e. $\delta E \ll E - E_n$) then fixing $E$ in $V_{\mathrm{nl}}$ to be the atomic valence eigenvalue is a reasonable approximation. In fact we would like to make the behaviour of the pseudopotential follow that of the real potential to first order in $E$, and this can be achieved by constructing a norm-conserving pseudopotential (see section 3.3.3).

3.3.2 Scattering approach

For a fuller discussion of the theory of scattering see [62]. Consider a plane-wave with wave-vector ${\bf k}$ scattering from some spherically-symmetric potential localised within a radius $r_{\mathrm c}$ and centred at the origin. The incoming plane-wave can be decomposed into spherical-waves by the identity

$$\exp({\mathrm i} {\bf k} \cdot {\bf r}) = 4 \pi \sum_{\ell=0... ...) Y_{\ell m}^{\ast}({\hat {\bf k}}) Y_{\ell m}({\hat {\bf r}})$$ (3.62)

where ${\hat {\bf k}}$ denotes a unit vector in the direction of ${\bf k}$. These spherical- or partial-waves are then elastically scattered by the potential which introduces a phase-shift $\delta_{\ell}$, which is related to the logarithmic derivative of the exact radial solution for given $\ell$ and energy $E={\textstyle{1 \over 2}}k^2$ within the core, evaluated on the surface of the core region:

$$L_{\ell}(E) = \left[ \frac{\mathrm d}{{\mathrm d}r} \log[R_... ...athrm c}) - \tan({\delta_{\ell}}) n_{\ell}(k r_{\mathrm c})} .$$ (3.63)

$j_{\ell}$ and $n_{\ell}$ denote the spherical Bessel and von Neumann functions respectively, and the radial wave-function $R_{\ell}(r,E)$ is related to the solution of the Schrödinger equation with angular momentum state determined by the good quantum numbers $\ell$ and $m$, and energy $E$, within the core region, $\psi_{\ell m}({\bf r},E)$ by

$$\psi_{\ell m}({\bf r},E) = R_{\ell}(r,E) Y_{\ell m}({\hat {\bf r}}) .$$ (3.64)

The phase-shifted spherical-waves can then be recombined to form the total scattered wave. We can define a reduced phase-shift $\eta_{\ell}$ by

$$\delta_{\ell} = n_{\ell} \pi + \eta_{\ell}$$ (3.65)

which has the same effect (the scattering amplitude depends on $\exp (2 {\mathrm i} \delta_{\ell})$ so that factors of $\pi$ in $\delta_{\ell}$ have no effect) and fix $n_{\ell}$ by requiring $\eta_{\ell}$ to lie in the interval $0 \leq \eta_{\ell} \leq \pi$. The integer $n_{\ell}$ counts the number of radial nodes in $R_{\ell}(r,E)$, two in the case of figure 3.2, and is thus equal to the number of core states with angular momentum ${\ell}$.

The pseudopotential is then defined as the potential whose complete phase-shifts are the reduced shifts $\eta_{\ell}$ so that the radial pseudo-wave-function has no nodes and thus the potential has no core states. The scattering effect of this potential is the same as the original potential. We note again the energy-dependence of the phase-shifts so that for a good approximation it will be necessary to match these phase-shifts to first order in the energy so that it is accurate over a reasonable range of energies, a property which results in good transferability of the pseudopotential i.e. it is accurate in a variety of different chemical environments. The non-local nature is also exhibited since different angular momentum states are scattered differently.

3.3.3 Norm conservation

The conditions of a good pseudopotential are that it reproduces the logarithmic derivative of the wave-function (and thus the phase-shifts) correctly for the isolated atom, and also that the variation of this quantity with respect to energy is the same to first order for pseudopotential and full potential^3.2. Having replaced the full potential by a pseudopotential, we can once again solve the Schrödinger equation in the core region to obtain the pseudo-wave-function, with radial part $R_{\mathrm{ps},\ell}(r,E)$.

Pseudopotential generation has itself been the subject of a great deal of study in the past (see [65,66,67,68,69,70,71,72,73]) and in this work we have chosen to use those pseudopotentials generated by the method of Troullier and Martins [74]. With one notable exception [75], all of the recent methods have used norm-conservation to guarantee that the phase-shifts are correct to first order in the energy (correction to higher orders is also possible [76]).

Consider the following second-order ordinary differential equations which are eigenvalue equations for the same differential operator but with different eigenvalues:

$${y_1}''(x) + p(x) {y_1}'(x) + q(x) y_1(x) = \lambda_1 y_1(x)$$

$${y_2}''(x) + p(x) {y_2}'(x) + q(x) y_2(x) = \lambda_2 y_2(x) .$$ (3.66)

In the context of homogeneous differential equations, the quantity known as the Wronskian is defined by

$$W(x) = y_1(x) {y_2}'(x) - y_2(x) {y_1}'(x)$$ (3.67)

and can be calculated according to

$$W(x) = W_0 \exp \left[ - \int^x {\mathrm d}x' p(x') \right]$$ (3.68)

in which the constant $W_0$ is arbitrary and of no consequence.

Following a similar analysis which leads to equation 3.68 for the quantity defined in equation 3.67 but for the functions which solve equations 3.66 we obtain

$$W(x) = \left[ (\lambda_2 - \lambda_1) \int^x {\mathrm d}x' y... ...+ W_0 \right] \exp \left[ - \int^x {\mathrm d}x' p(x') \right]$$ (3.69)

and note that the Wronskian can also be rewritten in terms of logarithmic derivatives:

$$W(x) = y_1(x) y_2(x) \frac{\mathrm d}{{\mathrm d}x} \Bigl\{ \log[y_2(x)] - \log[y_1(x)] \Bigr\} .$$ (3.70)

Using equations 3.69 and 3.70 in the case of the Schrödinger equation for the radial wave-function $R_{\ell}(r,E)$, by making the replacements

$$x \rightarrow r \quad ; \quad p(x) \rightarrow \frac{2}{r} \... ...r^2} \right] \quad ; \quad \lambda \rightarrow -2E , \nonumber$$

and using limits $r=0,r_{\mathrm c}$ we obtain

$$\left[ r^2 R_{\ell,1}(r) R_{\ell,2}(r) \frac{\mathrm d}{{\ma... ...r_{\mathrm c}} {\mathrm d}r ~r^2 R_{\ell,1}(r) R_{\ell,2}(r) .$$ (3.71)

Rearranging, multiplying by $-2 \pi$ and noting that the lower limit on the left-hand side contributes nothing because of the $r^2$ factor:

$$-2 \pi r_{\mathrm c}^2 \frac{R_{\ell,1}(r_{\mathrm c}) R_{\... ...r_{\mathrm c}} {\mathrm d}r ~r^2 R_{\ell,1}(r) R_{\ell,2}(r) .$$ (3.72)

Finally, taking the limit $E_2 \rightarrow E_1$ so that $R_{\ell,2}(r) \rightarrow R_{\ell,1}(r)$ and interpreting the left-hand side as a derivative with respect to energy we obtain the desired result:

$$-2 \pi r_{\mathrm c}^2 R_{\ell}^2(r_{\mathrm c}) \frac{\math... ...= 4 \pi \int_0^{r_{\mathrm c}} {\mathrm d}r ~r^2 R_{\ell}^2(r)$$ (3.73)

i.e. the first energy-derivative of the logarithmic derivative evaluated at the core radius (and hence the phase-shift) is related directly to the norm of the radial wave-function within the core region. Thus if the pseudo-wave-function is norm-conserving such that

$$4 \pi \int_0^{r_{\mathrm c}} {\mathrm d}r ~r^2 R_{\ell}^2(r)... ...^{r_{\mathrm c}} {\mathrm d}r ~r^2 R_{{\mathrm{ps}},\ell}^2(r)$$ (3.74)

then the phase-shifts of the pseudopotential will be the same as those of the real potential to first order in energy, and this can be achieved by making the pseudo-wave-function identical to the original all-electron wave-function outside the core region.

3.3.4 Kleinman-Bylander representation

We have seen that it is necessary to use a non-local pseudopotential to accurately represent the combined effect of nucleus and core electrons, since different angular momentum states (partial waves) are scattered differently. In general we can express the non-local pseudopotential in semi-local form

$${\hat V}_{\mathrm{ps}} = {\hat V}_{\mathrm{loc}} + \sum_{\el... ...ert \ell m \rangle \delta {\hat V}_{\ell} \langle \ell m \vert$$ (3.75)

in which $\vert\ell m \rangle$ denotes the spherical harmonic $Y_{\ell m}$. The choice of local potential ${\hat V}_{\mathrm{loc}}$ is arbitrary, but in general the sum over $\ell$ is truncated at a small value (e.g. $\ell=2$) so that the local part is required to represent the potential which acts on higher angular momentum components.

This semi-local form suffers from the disadvantage that it is computationally very expensive to use, since the number of matrix elements which need to be calculated scales as the square of the number of basis states, and this is generally too costly. In section 5.6.1 we will describe how this problem can be overcome analytically with a certain set of localised basis functions, but the most common solution, and one which we have also implemented for consistency, is to use the Kleinman-Bylander separable form [77]

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

where $\vert \phi_{\ell m} \rangle$ is an eigenstate of the atomic pseudo-Hamiltonian. This operator acts on this reference state in an identical manner to the original semi-local operator ${\hat V}_{\mathrm{ps}}$ so that it is conceptually well-justified, but now the number of projections which need to be performed scales only linearly with the number of basis states. This separable form can in fact be viewed as the first term of a complete series [78].