Chapter 9. From Atoms to Solids: The Pseudopotential Method.

[Note]Development note

Development of this chapter has been postponed.

The inner core of electrons is tightly bound to the nucleus and have a relatively high kinetic energy. Indeed, for atoms with high atomic number relativistic effects can be substantial. Tightly bound electrons with large velocities translates into wavefunctions with a high degree of curvature and a broad spectrum in momentum space. To represent these wavefunctions would require a large number of basis functions (especially in the case of plane-waves), but at the same time these states have only an indirect role in chemical processes and are insensitive to the chemical environment. The valence state wavefunctions also have a high degree of curvature in the core region, which manifests itself as oscillatory behavior. Unlike the core states, however, the amplitude of the valence wavefunctions decrease rapidly near the core. A key assumption of the pseudopotential method is that the core states in a chemical environment are the same as in a free atom. The main objective is to replace the need for calculating the core electron states with a pseudopotential that reproduces the effects of the core electrons.

Before we start a general discussion we would like to mention a few key terms and points to keep in mind, which we hope will become clearer as one progresses. These are:

  • Transferability: the idea that a pseudopotential should reproduce all-electron results in a variety of chemical environments.

  • Norm-conserving [132]: the notion that outside the core region the pseudo valence wavefunctions and the all-electron wavefunctions should be identical; and the integrated charge within the "core" region calculated with the pseudo wavefunctions should equal the real charge within that region. An idea introduced by Topp and Hopfield [131] to improve transferability.

  • Consistency: pseudopotentials constructed in the framework of the local-density approximation (LDA) should not be used in calculations using the generalized gradient approximation (GGA) and vice versa.

  • Local, semilocal and nonlocal potentials: A local potential is only a function of the distance from the nucleus, r. Both the nonlocal and the semilocal pseudopotentials are angular momentum dependent, but a semilocal pseudopotential is only nonlocal in angular coordinates, not in the radial coordinate [134].

In most formulations both the pseudopotential and the pseudo wavefunctions are equal to their all-electron counterparts beyond some cut-off radius   r c . At radial distances less the cut-off radius   r c , the oscillatory behavior of the all-electron (AE) valence wavefunction is replaced by a smoothly varying function. Thus,   r c   must be large enough that all the nodes found in the AE valence wavefunction lie within this radius. At the same time, however, with a smaller   r c   it is easier to replicate the behavior of the AE valence wavefunction and therefore improve transferability.

For hydrogenlike atoms, the root-mean-square velocity is given by [Gasiorowicz,p.271;Levine,p.525]

where   β   is related to the Lorentz factor   γ , a measure of the importance of relativistic effects, as follows

Thus, for atoms with a large nuclear charge Z relativistic effects cannot be ignored if accurate results are required. In 1926, Schrödinger presented a simple relativistic extension to his wave equation valid for spin zero particles. In 1928, Dirac presented a new relativistic formulation of quantum mechanics that incorporated four component wavefunctions (spinors), which accurately predicted the existence of positrons and provided a stepping stone to the development of quantum field theory. The Dirac equation with a central field can be written



both   β ^   (no relation to the   β v / c   used previously) and   α ^   are 4×4 dimensional matrices, and   σ ^   are the 2×2 dimensional Pauli spin matrices [Eq. (2.21)]. A solution to the Dirac equation can be written as

where   ψ   is a four-component spinor and   ϕ jm l   are r-independent two-component spinors. By writing   ϕ jm l   in terms of spherical harmonics [Sakurai,p.124] and applying the Dirac equation to   ψ   one can obtain the following radial wave equations [Schiff,p.484; Sakurai,p.125; Itzykson,p.78], [137]

where, following Bachelet, Hamann, and Schüter (BHS) [137], the rest mass m is not included in the energy   ε i , i.e. ε i = E - α - 2 . We have also used atomic units with   = m = e = 1   and   c = α - 1 = 137.04 .

In the case of hydrogen-like atoms, an exact solution to the Dirac equations exist [Itzykson,p.79]. Here however, we will consider the nonrelativistic limit in which case we can write Eq. (9.8) as

valid for valence electrons (small   ε i ) outside the core region [small V(r)]. By substituting Eq. (9.8) into Eq. (9.7) we obtain the Schrödinger-like approximation to the Dirac equation for a central field [135] [136] [137]

which is valid up to, but not including, terms of order   α 2 .   - ℏκ   are eigenvalues of the operator   Κ ^ ,  which is given by [Schiff,p.483; Sakurai,p.123; Baym,p.565]


The operator   Κ ^   commutes with the relativistic Hamiltonian and the total momentum operator   J ^   and   κ   is a good quantum number.   Κ ^   describes the extent that the spin and orbital momentum are aligned. The values of   κ   can be related to the total angular momentum quantum number j by writing

where we have used the identities


If we recall that

then we can use Eq. (9.13) to find the following relationship between the eigenvalues of   Κ ^ 2  and  J ^ 2

which has the two solutions

Although relativistically l is no longer a good quantum number it can be shown that in the nonrelativistic limit the quantum number,   κ ,  is related to the orbital angular momentum quantum number, l. The two components of   ψ   in Eq. (9.6),   ψ A  and  ψ B , are separately eigenfunctions of the operators   L ^ 2 , L ^ z , and  S ^ z .  Applying the two block diagonal elements of Eq. (9.11) to the two components of   ψ

together with

leads to

The minus sign in the last equation corresponds to the upper component of the wave function (i.e.   ψ A ) and the plus sign corresponds to the lower part (i.e.   ψ B ). For G(r) we are interested in the solution with the minus sign and thus

Next, we will need to find solutions to Eq. (9.10) with Eq. (9.22) using pseudopotential methods.

There are several nice review articles of which Fuchs and Scheffler [128] is one of the most recent and up-to-date. In that review, Fuchs and Scheffler [128] also introduce their ab initio pseudopotential package fhi98PP, which can be used to produce pseudopotentials. Other nice review articles include ones by Pickett [127], Heine [125], Cohen [126], which are not as up-to-date but offer a good background discussion. There are also several good reviews that are directed a bit more toward those with a chemistry background. Among these are ones by Krauss and Stevens [147] and Frenking et al. [149].

Articles specific to this implementation:

[Back to Table of Contents] [Top of Chapter 9] []

Troullier and Martins introduced additional constraints to produce smoother pseudopotentials, which require a smaller set of basis functions resulting in greatly reduced computation times. Following Kerker [133] they wrote the radial pseudo wavefunction as

where the polynomial p(r) is given by

The order n of the p(r) was taken to be n = 10 by Troullier and Martins [141] compared to n = 4 in Kerker's work [133].

[Back to Table of Contents] [Top of Chapter 9] []

The total HGH pseudopotential is given by

where the local part of the pseudopotential is

and the nonlocal part is

The contribution from the spin-orbit coupling is

where   Y l , m   are the spherical harmonics, l is the angular momentum quantum number,   r loc   gives the range of the Gaussian ionic charge distribution,   Z ion   is the ionic charge (the charge of the nucleus minus the charge of the valence electrons), and the projectors are given by

The parameters for the various chemical elements are given in tables with the following format

For example, the LDA pseudopotential parameters for hydrogen in atomic units are

and the LDA pseudopotential parameters for carbon are

In the octopus program these tables are stored in files labeled element.hgh, or element_sc.hgh for semicore pseudopotentials [145][146]. ABINIT has six options available for pseudopotentials including Hartwigsen-Goedecker-Hutter LDA pseudopotentials, and Hartwigsen-Goedecker-Hutter GGA (PBE) pseudopotentials for a limited number of chemical elements. The parameters are stored in files labeled element.psphgh and element.pbe_hgh, respectively, and a couple examples are included with the installation and can be found in the ~ABINIT/Psps_for_tests directory.

[Back to Table of Contents] [Top of Chapter 9] []