Normconserving Pseudopotentials

Normconserving Pseudopotentials

Basic idea

A substantial fraction of present-day electronic structure calculations, covering such a variety of topics as magnetic materials, vacancies and impurities, quasi-crystals, surfaces, clusters and biomolecules (see e.g. [1-6]), relies on the pseudopotential (PP) framework. The PP approach is based on the observation that most physical and chemical properties of atoms are determined by the structure and dynamics of the atomic valence states. This is true in particular for the formation of chemical bonds, but also for the magnetic behavior and for low-energy excitations. On the other hand, these properties are affected only indirectly by the electrons filling the deeper-lying levels, the so-called core electrons: The presence of the core electrons essentially restricts the Hilbert space available to the valence electrons in some binding or excitation process via the orthogonality requirement.

The main reason for the limited role of the core electrons in many processes is the spatial separation of the core and valence shells which originates from the comparatively strong binding of the core electrons to the nucleus. This effect is illustrated in Figure 1 where the radial orbitals Pnl of the oxygen atom are plotted.

oxygen: radial all-electron orbitals
Figure 1: Radial orbitals of oxygen: Spatial separation of K-, L- and M-shell.

One immediately recognizes the separation of the 1s-state with its r-expectation value of 0.2 Bohr from the L-shell orbitals with r-expectation values of roughly 1.2 Bohr. The corresponding eigenvalues differ by an order of magnitude. Similarly, the L-shell is well separated from the unoccupied M-shell orbitals. The spatial separation allows a clear decomposition of the complete spectrum into core, valence and higher-lying excited states.

The interatomic distances in molecules or solids, however, are controlled by the spatial extension of the most weakly bound occupied orbitals. As a consequence, the much more localized core electrons have little overlap with the wavefunctions of electrons from neighboring atoms and thus essentially do not interact with them. They experience the presence of neighboring atoms only indirectly via the deformation of the valence orbitals. This repercussion, however, is rather weak due to the strong binding of the core states to the nucleus, so that the form of the core orbitals remains almost unchanged. There are only few situations in which the core states are directly affected by neighboring atoms, so that the PP concept breaks down. This is typically the case if charge centers are present and the core states are not sufficiently attracted by their own nucleus. The deformation of the 1s-orbital of lithium in LiF may serve as an example.

The precise decomposition of the bound state spectrum into core, valence and excited states is sometimes more difficult than for oxygen. From the arguments given above it is clear that problems will come up as soon as the spatial and/or energetic separation between two shells breaks down. This is, for instance, the case for the 3d transition metal elements, for which the most weakly bound 4s-state is not well-separated from the 3d-states. As a result the overlap of the 3d-orbitals with the orbitals from neighboring atoms is too large to be ignored. In this situation it is usually necessary to consider the complete shell, to which the critical d- or f-state belongs, as valence space, rather than just the critical state itself. These more strongly bound, but nevertheless relevant states are called semi-core states.

It is obvious that the quantum mechanical description of poly-atomic systems is simplified considerably if part of the spectrum is kept out of the calculation: As such calculations inevitably utilize some form of basis set for the representation of the wavefunctions, the computational effort is reduced substantially if only the valence states have to be expanded in terms of basis functions. The elimination of the core electrons and, at the same time, the small scale structure of the valence states allows the use of much smaller and much simpler basis sets, as the plane-wave basis usually employed in the context of quantum molecular dynamics. This is the basic idea behind the PP concept: Only the valence electrons are treated dynamically, i.e. represented by wavefunctions for which some Schrödinger-type equation is solved, while the core states are completely eliminated from the poly-atomic calculation. The presence of the core electrons, which essentially restrict the Hilbert space available to the valence electrons, must therefore be simulated by some other means, the pseudopotential: The orthogonality requirement is recast in the form of an effective potential. For each atom in the poly-atomic system there is a corresponding effective potential, responsible for representing the core states and the nucleus of this particular atom. The molecular or crystalline PP is then obtained by a linear superposition of the atomic PPs.

Ideally the atomic PP simply projects out all core states from the complete Fock space, while leaving the remaining spectrum unchanged. The exact formulation of such a projection potential is possible in principle, but the projection operators involved rely on an explicit representation of the core states. In the standard form of the PP approach the ambitions are thus cut back a little bit: Rather than requiring the PP to generate exactly the same atomic spectrum as the all-electron (AE) Hamiltonian (apart from the core states), one is satisfied with an accurate reproduction of just the valence states.

The term 'accurate reproduction' requires some explanation: Clearly, in the region of space in which most of the electronic norm is concentrated both orbitals must be very close, if not identical. On the other hand, the form of the valence orbitals in the core region, where the core electrons are moving, is less relevant: Otherwise the core states themselves would play a more important role. In the core region one can thus allow the pseudoorbital (PO) to differ from the all-electron orbital (AO) without loosing too much accuracy. In the case of our example, the oxygen atom, the core region roughly extends to 0.5 Bohr (see Figure 1).

Of course, one can not expect a single multiplicative PP to reproduce several valence AOs, i.e. to simulate several projection operators simultaneously. In order to give this statement a more precise form one considers some atomic valence shell n. The standard form of the AOs is

\begin{displaymath} \phi_{nlm}(\bm{r}) = {P_{nl}(r) \over r} \, Y_{lm}(\Omega) \end{displaymath}

This form implies a spherically symmetric potential, which is ensured by a suitable spherical average in the case of open sub-shells. The associated POs must have the same form - any deviation from the angular characteristics of the original all-electron state would lead to a fundamentally different behavior. As a consequence a multiplicative PP necessarily has to be spherically symmetric, so that the single-particle Schrödinger equation for the POs has the usual form,

\begin{displaymath} \Bigg\lbrace - {1\over2m} {\partial^2\over\partial r^2} + {l... ...2mr^2} + v(r) \Bigg\rbrace P_{nl}(r) = \epsilon_{nl} P_{nl}(r) \end{displaymath}

However, for given l the lowest energy solution of this radial differential equation is always nodeless (see e.g.[7]). In a PP approach thus all radial orbitals of the valence shell must be nodeless, as for each l all lower-lying states have to be projected out by the PP. In the case of oxygen, for instance, the 2s-PO must be nodeless. An example for such a nodeless Pnlps which nevertheless agrees with the corresponding AO Pnlae in the relevant part of space, is shown in Figure 2.
oxygen: pseudo versus all-electron 2s orbital
Figure 2: Radial 2s orbital of oxygen: Pseudo- versus all-electron orbital.


At the same time the 2p-PO is nodeless (see Figure 3).
oxygen: pseudo versus all-electron 2p orbital
Figure 3: Radial 2p orbital of oxygen: Pseudo- versus all-electron orbital.


It is, however, almost impossible to produce two (or even more) nodeless orbitals in the same energy range with only a single spherical potential v(r), as for fixed v(r) only the angular momentum term can generate differences. Consequently, the total atomic PP usually consists of several components, one for each angular momentum present in the valence space. The standard form utilized in the context of density functional theory (DFT) is

\begin{displaymath} \langle\bm{r}\vert \hat v_{\rm ps} \vert\bm{r}'\rangle = v_{... ...(r) \Big\rbrack \sum_{m=-l}^l Y_{lm}(\Omega) Y_{lm}^*(\Omega') \end{displaymath}

where lmax represents the highest angular momentum present among the valence states and a suitably chosen potential vloc has been extracted from the individual l-components. It is the PP experienced by all states with angular momentum beyond lmax. Usually one of the PP components is used for vloc.

Construction of normconserving pseudopotentials

There exists a variety of schemes for the explicit construction of PPs. In the context of DFT one usually relies on the concept of norm-conserving PPs [8]. On the basis of the results of an all-electron DFT calculation for the atom of interest, the norm-conserving PP for angular momentum l is chosen so that
  1. the resulting atomic valence PO is identical with the corresponding AO for all r larger than some l-dependent cut-off radius rc,l,

    \begin{displaymath} P_{nl}^{\rm ps}(r) = P_{nl}^{\rm ae}(r) \hskip20pt\forall~r\ge r_{{\rm c},l} \end{displaymath}

  2. the norm of the orbital is conserved,

    \begin{displaymath} \int_0^{r_{{\rm c},l}}dr\, P_{nl}^{\rm ps}(r)^2 = \int_0^{r_{{\rm c},l}}dr\, P_{nl}^{\rm ae}(r)^2 \end{displaymath}

Condition 1 automatically implies that the eigenvalue of the PO is identical with that of AO, as the eigenvalue determines the asymptotic decay of the orbitals.

Given the AE information, the actual construction of the PPs proceeds in two steps. In the first step a screened radial PP vps,lsc is generated: vps,lsc is the total single-particle potential which has the desired pseudo valence state as lowest eigenstate. The identity of the pseudo with the AE orbital for r>rc,l requires the screened PP to coincide with the total AE Kohn-Sham (KS) potential in this valence region (as both orbitals satisfy the same radial differential equation). Consequently only the form of the PP in the core region remains to be determined. For the extension of the PP into the core region different procedures have been suggested [8-12] among which the Troullier-Martins (TM) approach [12] seems to be the most widely used. Compared with the orginal concept of Hamann, Schlüter and collaborators [8] the TM scheme produces somewhat smoother PPs, so that smaller plane-wave basis sets can be utilized in applications.

In the TM scheme an ansatz for the PO in the core region is made, which

  1. has the correct rl+1 behavior at the origin,
  2. is smooth otherwise,
  3. depends on a number of adjustable parameters.
The parameters are determined by requiring continuity of the PO and its first four derivatives at the cut-off radius as well as norm-conservation. For oxygen the scheme described in this Section leads to the results plotted in Figures 2 and 3. The screened PPs associated with this PO are then obtained by inversion of the radial Schrödinger equation. An example for the screened PP and its relation to the original AE potential is given in Figure 4.
oxygen: screened pseudopotentials
Figure 4: Pseudopotential of oxygen: Screened pseudopotentials versus all-electron potential. As cut-off radii 1.0 Bohr and 0.9 Bohr have been used for the 2s and the 2p state, respectively.


In many cases the valence shell of an atom does not only consist of occupied states, but also contains unoccupied states. The most obvious example are the second row elements, for which the 3s and 3p orbitals are at least partially occupied, while the 3d orbital is empty. Unfortunately, the standard exchange-correlation (xc) functionals of DFT, the LDA and the GGA, do not yield bound 3d-levels for the ground state configurations of these atoms, as the corresponding xc-potentials decay exponentially for large r. One way to solve this problem is the use of a fractionally ionized configuration for the atom [8]: If the net positive charge on the atom is sufficiently large the critical orbital is bound. An alternative approach has been suggested by Hamann [11]. It is based on the observation that the actual form of the AE-orbital is only required for r<rc,l - for r>rc,l the PP is identical with the AE potential anyway. Moreover, the normalization of the orbital does not play any role for the final PP. The complete TM procedure can thus be carried through with any function obtained by outward-integration of the radial Kohn-Sham equations with the AE-potential, starting with appropriate boundary conditions at the origin. In the Hamann prescription one thus chooses some energy value in the appropriate range between the highest occupied KS state and zero and then performs an outward integration from 0 to some r-value beyond rc,l, using an arbitrary normalization. With the resulting function the standard procedure for the construction of the associated PP is gone through.

In a second step the interaction among the valence electrons (and their self-interaction) is eliminated from the screened PP by the so-called unscreening procedure. This step is required as PPs are usually applied to molecules or solids within a selfconsistent framework, typically nonrelativistic spin-density functional theory. In any non-perturbative framework, however, the differences between the atomic and the poly-atomic valence states are of crucial importance. The same is true for the interaction potentials between the valence states, most notably the direct electrostatic interaction. The screened PP contains the interaction among the atomic valence states, which is quite different from the interaction potential obtained with the molecular or bulk valence states. The atomic interaction potential must thus be eliminated from vps,lsc. The result is the unscreened PP which plays the role of the external potential in the poly-atomic calculation, analogous to the nuclear potential in the AE situation.

The simplest approach for this elimination is linear unscreening,

\begin{displaymath} v_{{\rm ps},lj}(r) = v_{{\rm ps},lj}^{\rm sc}(r) -v_{\rm H}([n_{\rm v,ps}];r) -v_{\rm xc}([n_{\rm v,ps}];r) \end{displaymath}

where the Hartree potential vH and the xc-potential vxc,

\begin{eqnarray*} v_{\rm H}([n];r) &=& 4\pi e^2 \Bigg\lbrace {1\over r} \int_0^r... ..._{\rm xc}([n];r) &=& {\delta E_{\rm xc}[n]\over\delta n(\bm{r})} \end{eqnarray*}

are evaluated with the pseudo valence density nv,ps constructed from the occupied valence POs. Of course, the xc-functional utilized for unscreening must be the same as that used in the underlying AE calculation. For oxygen linear unscreening produces the PPs of Figure 5.
oxygen: unscreened pseudopotentials
Figure 5: Pseudopotential of oxygen: Screened versus unscreened pseudopotential.

As is obvious from this plot, the unscreening potential sometimes introduces additional small-scale structure into the final PP. In fact, this is quite often the case if GGAs are used for the xc-energy functional. In the present example the peak structure at the origin is a consequence of the exact exchange functional which has been utilized both in the AE calculation and the unscreening procedure. This small-scale structure, however, is irrelevant for applications, as it can not be resolved by standard (plane-wave) basis sets.

Linear unscreening provides a complete elimination of the direct Coulomb interaction among the valence electrons, as the Hartree potential is linear in the density. The xc-energy functional Exc of DFT, on the other hand, is a nonlinear functional of the density. Linear unscreening thus implies a linearization of Exc with respect to the core-valence interaction. The quality of the linearization obviously depends on the spatial overlap between the core density and the valence xc-potential. In the case of the (semi-)local functionals LDA and GGA this overlap is directly determined by the overlap of core and valence density, which is nonzero only in the transition region between valence and core shells. The overlap does not vanish completely even for first row elements (compare Figure 1), so that the nonlinearity of Exc shows up even for these ideal PP candidates. Nevertheless, linear unscreening is often sufficient, at least as long as the accuracy requirements are not particularly high.

The nonlinearity of Exc can be taken into account via the so-called nonlinear core corrections (nlccs) [13]. In this approach the atomic core density nc is included in the xc-component of the unscreening potential, so that the complete atomic core-valence interaction is eliminated. As a consequence nc must then be added to the valence density of the poly-atomic system when the xc-potential of this system is calculated. Fortunately, one needs only that part of nc which actually has overlap with the valence density. This allows a smooth truncation of nc below some core-cut-off radius rnlcc. The unscreened PP including nlccs is then given by

\begin{displaymath} v_{{\rm ps},lj}^{\rm nlcc}(r) = v_{{\rm ps},lj}^{\rm sc}(r)... ...}([n_{\rm v,ps}];r) -v_{\rm xc}([n_{\rm v,ps}+n_{\rm c,ps}];r) \end{displaymath}

rnlcc is chosen so that most of the overlap between core and valence density is found for r>rnlcc.

The example of oxygen is given in Figure 5. One immediately observes the different asymptotics, reflecting the fact that vps,l asymptotically approaches -Zion/r, while the screened PP has the same behavior as the AE potential (which falls off like -1/r on the exact level and decays exponentially within the standard approximations of DFT).

So far, the formulation of PPs has been strictly nonrelativistic, so that the issue of heavy elements still remains to be addressed. In fact, the PP treatment is particularly attractive for heavy atoms for which many shells belong to the core states. In order to cover the complete periodic table the PP construction is thus usually based on relativistic AE results [14]. This leads to j-dependent PPs vps,lj, in accordance with the j-dependence of the individual valence states. However, while the core states directly experience the relativistic effects, the valence states are in most cases only indirectly affected by relativity via the orthogonality requirement (which is simulated by the PP). Consequently, poly-atomic PP calculations are usually performed in a nonrelativistic framework, even if heavy ions are present. The j-independent PPs required for these nonrelativistic calculations are then obtained from a j-average,

\begin{displaymath} v_{{\rm ps},l}(r) = \sum_{j=l\pm1/2} {2j+1\over4l+2} \, v_{{\rm ps},lj}(r) \end{displaymath}

A corresponding relativistic TM PP construction scheme has been developed in our group [15].

Transferability

The PP concept has been motivated by the inertness of the atomic core states in binding and low-energy excitation processes. Given this inertness, the ionic core of the atom provides a fixed potential in which the valence electrons are moving, independently of the system (atom, molecule or solid) which is considered. However, in poly-atomic systems the valence states undergo obvious modifications compared to the atomic valence orbitals, even if the poly-atomic core potential is given by a simple linear superposition of atomic core potentials. Most notably, the eigenenergies change when packing atoms together, which leads to bonding and anti-bonding states in molecules and to energy bands in solids. Thus, while PPs are designed to reproduce the valence AOs of some chosen atomic reference configuration (usually the ground state), it is not clear a priori that they will have the same property for all kinds of poly-atomic systems and for other atomic configurations. Consequently, one has to make sure that the PP is transferable from its atomic reference state to the actually interesting environment.

In order to address this question one has to analyze the sensitivity of the agreement between atomic POs and AOs to the specific eigenenergy in the single-particle equation. One finds that the variation of the logarithmic derivative P'nl(r)/Pnl(r) with the single-particle energy is determined by the norm contained in the sphere between the origin and r. This is true in particular in the neighborhood of one of the actual atomic eigenvalues εnl, i.e. for a bound atomic eigenstate. Thus, as soon as norm-conservation is ensured, the POs exactly reproduce the energy dependence of the logarithmic derivative of the AOs for r>rc,l. Consequently, one expects the POs to react as the AOs when the valence states experience some energy shift in a poly-atomic environment -- provided the underlying PPs are norm-conserving. This argument supporting the transferability of PPs emphasizes the importance of norm-conservation in a very explicit way.

In practice, it is nevertheless always recommendable to check the transferability explicitly by examination of some suitable atomic excitation process and of the binding properties of simple molecular or crystalline systems.

Convergence issues

The continuity of the PO and its first four derivatives at the cut-off radius ensures that the pseudo and AE orbitals agree well below the cut-off radius, which allows a more relaxed choice of rc,l for the TM-PPs than for the original BHS-scheme. An exemplary list of rather conservative, i.e. small, cut-off radii is given in Table 1.

Table 1: Conservative cut-off radii for first and second row atoms within the Troullier-Martins scheme.
  rc,l    [Bohr] rnlcc
         s               p               d           [Bohr]   
Li 2.8 3.0 -- 1.2
Be 1.7 1.5 -- 0.7
B 1.3 1.3 -- 0.6
C 1.2 1.1 -- 0.4
N 1.1 1.0 -- 0.3
O 1.0 0.9 -- 0.25
F 0.6 0.5 -- 0.2
Na 2.7 2.7 2.7 2.0
Mg 1.9 2.3 2.3 1.2
Al 1.9 2.3 2.3 1.0
Si 1.8 2.0 2.0 1.0
P 1.7 1.7 1.7 0.8
S 1.4 1.5 1.5 0.7
Cl 1.2 1.2 1.2 0.6

In computationally demanding applications often slightly larger cut-off radii are utilized. This leads to smoother PPs and thus reduces the required basis set size. Table 2 demonstrates the limited sensitivity of PP results to the choice of the cut-off radii.

Table 2: Sensitivity of PP results to choice of cut-off radius: Bond length Re, dissociation energy De (including zero-point energy) and harmonic frequency ωe of the nitrogen dimer for different sets of rc,l. The calculations are based on the LDA for the xc-functional [16] and linear unscreening.
rc,l     [Bohr] Re De ωe
       s               p           [Bohr]          [eV]           [cm-1]    
1.1 1.0 2.065 10.967 2381
1.3 1.2 2.068 10.935 2371
1.5 1.4 2.073 10.827 2342

On the other hand, the corresponding gain in computational efficiency is often not dramatic.

The second quality parameter of interest is the truncation radius of the core density, required for nlccs. The truncated core density nps,c becomes more and more strongly peaked, the smaller rnlcc is chosen. The spacing of the spatial grid in poly-atomic calculations, however, is often too low to resolve this small-scale structure of nps,c. A smaller rnlcc thus requires a finer grid, which means a higher energy cut-off in the case of plane-wave basis sets. This effect is demonstrated in Table 3.

Table 3: Convergence of pseudopotential results with cut-off energy Emax of plane wave basis set: Total energy of atomic nitrogen as well as bond length Re, dissociation energy De (including zero-point energy) and harmonic frequency ωe of the nitrogen dimer with and without nlccs as a function of cut-off energy and spatial resolution. All calculations are based on the plane-wave pseudopotential framework, utilizing a cubic supercell with a = 10 Bohr and the LDA for the xc-functional [16].
cut-off results
   Emax = Gmax2 / (2m)       π/ (2Gmax)       -Etot(N)    Re De ωe
  [Rydberg]     [Bohr]     [Hartree]     [Bohr]         [eV]           [cm-1]    
  without nlccs
50 0.22 9.669 2.085 10.508 2348
100 0.16 9.749 2.065 10.830 2354
200 0.11 9.751 2.065 10.830 2354
400 0.08 9.751 2.065 10.830 2354
  with nlccs
50 0.22 12.662 2.009 1.805 2154
100 0.16 12.197 2.060 10.452 2386
200 0.11 12.147 2.067 11.262 2368
400 0.08 12.142 2.067 11.346 2368

It is obvious that for nitrogen one needs a cut-off of at least 100Rydberg to obtain binding energies with an error below 1eV. On the other hand, the structural data like bond lengths and force constants are converged much earlier. The situation is more favorable if PPs without nlccs are used. In this case the cut-off energy of 50Rydberg is completely sufficient to obtain converged results. For pseudopotentials including nlccs it is thus not the representation of the POs which limits the accuracy of plane-wave calculations, but rather that of the core density.

The convergence properties can also be analyzed from a different perspective. In the standard plane-wave pseudopotential scheme the Fourier representation of the POs is truncated at some cut-off momentum Gmax. This truncated Fourier transform deviates from the original PO. The difference between the truncated Fourier and the original PO is illustrated in Figure 6 for the case of the 2p-state of oxygen.

oxygen: exact 2p pseudoorbital versus Fourier representation
Figure 6: Radial 2p orbital of oxygen: Original pseudoorbital versus Fourier representations obtained for different cut-off momenta Gmax.

It is obvious that for an accurate Fourier representation of this orbital a cut-off energy Emax of roughly Gmax2 / (2m) = 100 Rydberg is required. Particularly critical is the Fourier representation of the rather localized d- and f-states of the transition metal elements. This is demonstrated in Figure 7 in which the 3d state of iron is plotted. One needs an Emax of more than 200 Rydberg to achieve good agreement between the original PO and its Fourier transform.

iron: exact 3d pseudoorbital versus Fourier representation
Figure 7: Radial 3d orbital of iron: Original pseudoorbital versus Fourier representations obtained for different cut-off momenta Gmax (rc,3d = 0.8 Bohr).



Illustrative Results

When coming across the PP concept for the first time, often scepticism prevails. In particular the fact that the nodal structure of the valence states in the core region is irrelevant is difficult to accept. In order to overcome this scepticism the high level of accuracy which can be achieved with PPs is demonstrated in this Section by providing a number of prototype results. In order to focus on the properties of the PPs mainly atoms and diatomic molecules are considered, for which the technical details of the calculations can be chosen so that they are beyond doubt. At the end of the Section also some results for solids are given, which, however, suffer a litte bit from the different technical details of the PP and AE calculations.

All calculations are based on the LDA for the xc-functional, using the parametrization of Vosko, Wilk and Nusair [16]. In order to be consistent with the nonrelativistic AE data, both the PP construction and all calculations are strictly nonrelativistic. Whereever present, spin-polarization is taken into account.

The first class of results which is of interest in this context are excitation energies. In Table 4 we list the s-p transfer energies of first and second row atoms. All calculations utilize spherically averaged densities.

Table 4: s-p transfer energies of selected first and second row atoms: PP versus AE LDA results [17] (all energies in eV).
atom ns↓ → np↑ ns↓ → np↓
      AE         PP     PP+nlcc     AE         PP     PP+nlcc
Be 2.47 2.37 2.46 3.49 3.49 3.49
B 3.23 3.00 3.21 5.57 5.57 5.56
C 4.05 3.68 4.01 8.02 8.02 8.00
N       10.86 10.87 10.84
O       14.38 14.39 14.37
F       18.22 18.24 18.22
Mg 2.80 2.69 2.80 3.44 3.40 3.44
Al 3.61 3.41 3.61 5.01 4.96 5.01
Si 4.31 3.98 4.29 6.57 6.48 6.57
P       8.14 8.02 8.14
S       9.95 9.89 9.95
Cl       11.79 11.76 11.78
           

Whereever possible, two types of excitations are shown: In the first case, the minority spin (↓) s-electron of the valence shell is transferred to the majority spin (↑) p-orbital, so that the excited state has a larger spin-polarization (magnetic moment) than the ground state. In the second case, the s↓ electron goes into the p↓ state, leaving the magnetic moment unchanged. For the latter transition even the PPs not including nlccs yield excitation energies which agree excellently with the AE data. The PPs perform particularly well if the net spin-polarization is small. On the other hand, for the transitions involving a spin-flip the AE data are somewhat underestimated by the PPs without nlccs. The maximum error is found for C and Si, for which the excited state exhibits the highest possible spin-polarization -- the stability of the high-spin excited state is slightly overestimated by the PP calculation. After inclusion of nlccs the PP excitation energies are more or less identical with the AE numbers. Taking into account the (unpolarized) core density in the xc-functional reduces the energy gain obtained by an increase of spin-polarization, which leads to improved ground state and excited state energies. Nevertheless, even the level of accuracy reached without nlccs is completely sufficient for many purposes.

As a third example we consider molecular ground state properties. The spectroscopic constants obtained for first row diatomic molecules are listed in Table 5.

Table 5: Bond length Re, dissociation energy De (including zero-point energy) and harmonic frequency ωe of first row diatomic molecules: PP versus AE LDA results [17]. The p-PP component has been used for vloc.
  mode Re De ωe
    [Bohr] [eV] [cm-1]
B2 AE 3.032 3.855 1042
3Σ PP 3.023 3.790 1037
PP+nlcc 3.034 3.852 1040
C2 AE 2.353 7.251 1890
1Σ PP 2.350 6.923 1877
PP+nlcc 2.355 7.223 1885
N2 AE 2.068 11.593 2390
1Σ PP 2.065 10.967 2381
PP+nlcc 2.068 11.494 2387
O2 AE 2.273 7.591 1626
3Σ PP 2.279 7.240 1613
PP+nlcc 2.279 7.506 1612
F2 AE 2.614 3.397 1062
1Σ PP 2.619 3.234 1058
PP+nlcc 2.618 3.375 1059
LiH AE 3.032 2.643 1375
1Σ PP 2.917 2.720 1320
PP+nlcc 3.037 2.597 1322
CO AE 2.128 12.967 2183
1Σ PP 2.125 12.589 2164
PP+nlcc 2.128 12.900 2181
average PP 0.026 0.235 15
deviation PP+nlcc 0.003 0.044 12

For all molecules the PP results obtained with and without nlccs are listed together with the AE data. For both sets of molecules the average absolute deviations of the PP results from the AE reference values are given as well.

Table 5 documents the well-known success of norm-conserving PPs. Even without nlccs the average error in the bond length is no larger than 0.03 Bohr. This deviation is reduced to less than 0.01 Bohr by inclusion of nlccs. The nlccs are somewhat more important for the dissociation energy. Although the average error of 0.24 eV found for the first row molecules without nlccs is sufficient for many purposes, it is not completely satisfactory. This error mainly originates from molecules whose constituents are ``strongly'' spin-polarized, with Nitrogen being the prime example. For these systems the dissociation energy is generally underestimated without nlccs, as a consequence of an overestimation of the atomic ground state energy. This is consistent with the observation that quite generally the energy of high-spin states is overestimated by PPs without nlccs. These difficulties are resolved by inclusion of nlccs: The resulting average error of 0.05 eV is much smaller than the level of accuracy which can be achieved in calculations for more complex molecules.

As pointed out earlier, the PP-treatment of transition metal elements is complicated by the large overlap of the (n-1)d orbital with the ns orbital. In this case it is necessary to include the semi-core states in the valence space, if fully reliable results are desired. The success of this procedure is demonstrated in Table 6 in which the spectroscopic constants for some diatomic molecules containing transition metal elements are listed.

Table 6: Bond length Re, dissociation energy De (including zero-point energy) and harmonic frequency ωe of diatomic molecules containing transition metal elements: PP versus AE LDA results [15]. The s-PP component has been used for vloc in the case of the transition metal elements, while the p-PP has been utilized for oxygen.
  mode Re De ωe
    [Bohr] [eV] [cm-1]
Fe2 PP 3.68 3.91 440
7Δ PP+nlcc 3.68 4.31 440
AE 3.68 4.38 497
Ni2 PP 3.85 3.56 363
3Σ PP+nlcc 3.85 3.49 362
AE 3.87 3.64 354
FeO PP 2.99 6.65 962
5Δ PP+nlcc 2.99 7.00 968
AE 3.01 7.06 957
CrO PP 3.01 6.06 952
5Π PP+nlcc 2.99 6.38 978
AE 3.02 6.15 976

For the d-elements always the complete M-shell is included in the valence space, i.e. for iron the valence space is given by 3s2 3p6 3d6 4s2. Again good agreement between AE and PP results is observed, if one takes into account that the AE data have been generated with completely different basis sets as the PP results. As to be expected the inclusion of nlccs (for the L-shell) is particularly important for these elements.

In order to illustrate the accuracy of the charge densities resulting from PP calculations, the dipole moment of diatomic FeO is considered in Table 7.

Table 7: Dipole moment for the 5Δ state of FeO: AE versus PP (without nlccs) LDA results for various internuclear distances R.
R AE PP
[Bohr] [Debye]
2.90 3.92 3.93
3.00 4.27 4.29
3.10 4.59 4.61

Even without nlccs the dipole moments of the PP and the AE calculation are essentially identical.

Finally, in Table 8 some exemplary results for the cohesive properties of metals and semiconductors are given.

Table 8: Equilibrium lattice constant a, cohesive energy Ecoh, bulk modulus B and magnetic moment M (at the equilibrium lattice constant) of selected solids: PP versus AE linearized-augmented-plane-wave LDA results [18,19].
system mode a Ecoh B M
     [Bohr]   [eV/atom]   [GPa]   [μB/atom] 
Na AE [18] 7.65   9.2  
  PP [19] 7.52 1.28 8.7  
  PP+nlcc [19] 7.65 1.22 9.1  
Al AE [18] 7.52   84  
  PP 7.47 4.08 86  
PP+nlcc 7.48 4.08 88  
Fe AE (rel.) 5.21   266 2.08
(bcc) PP (rel.) 5.21 6.36 276 2.15
PP+nlcc (rel.) 400/26 5.18 6.63 192 1.97
Si AE [18] 10.22 5.28 96  
  PP [19] 10.17 5.34 94  
  PP+nlcc [19] 10.19 5.32 94  
Ge AE [18] 10.63 4.54 78  
  PP [19] 10.51 4.75 73  
  PP+nlcc [19] 10.58 4.58 71  
GaAs AE [18] 10.62 7.99 74  
  PP [19] 10.39 8.68 79  
  PP+nlcc [19] 10.53 8.15 75  

It seems worthwhile to emphasize that the comparison between AE and PP results is much more difficult in the case of solids than for diatomic molecules: While the same code and high quality basis sets have been employed for the AE and PP entries in Tables 4-7 the computational details of the AE and PP calculations of Table 8 are quite different (in particular, the basis sets). Nevertheless the overall agreement between PP and AE data is rather good. The results are consistently improved by the inclusion of nlccs.

Exact exchange

The concept of linear unscreening relies on the locality of the density-dependence of the xc-functional: Otherwise the local subtraction of the xc-potential obtained from the pseudo valence density does not eliminate all valence xc-interaction effects, as vxc[nv,ps] is not identical with the corresponding xc-potential resulting from the valence AOs even for r>rc,l. While the standard xc-functionals LDA and GGA have such a local form, this locality is lost as soon as the exact exchange of DFT is used. The fully nonlocal structure of its kernel lets the valence density in the core region affect the exchange potential in the valence region. As a consequence, the differences between the POs and the corresponding AOs for r<rc,l show up in the unscreened PP. Unfortunately, the net effect is a spurious long-range structure in the unscreened PP, which decays so slowly that, in the typical range of interatomic distances, it mimics a fractional ionic charge [21]. This artificial charge then spoils any structural optimization. This is true in particular for larger molecules and solids in which each constituent atom has several neighbors.

The long-range structure thus has to be damped which is most easily done by multiplication of the unscreened PP with some damping function [21]. However, any a posteriori damping of the unscreened PP introduces additional parameters for which no physically motivated values are available. In fact, suitable choice of the damping parameters allows to generate essentially any desired result for poly-atomic systems [17]. A systematic scheme for the elimination of the long-range structure has been suggested by our group [17]. In this scheme a parameter-free, self-consistent procedure based on the standard norm-conserving PPs is utilized for the PP construction. With these asymptotically correct PPs the quality of PP results obtained with the exact exchange is as high as that of the corresponding LDA PP data [17]. As an illustration Table 9 lists the spectroscopic constants of the same first row diatomic molecules as shown in Table 5.

Table 9: Bond length Re, dissociation energy De (including zero-point energy) and harmonic frequency ωe of first row diatomic molecules: Standard and selfconsistent (SC-) PP versus AE results for the exact exchange. In the case of the PP data only the deviation from the corresponding AE value is given.
  mode Re De ωe
    [Bohr] [eV] [cm-1]
Li2 AE 5.266 0.168 338
  PP -0.004 0.001 -1
  SC-PP -0.009 0.003 -1
B2 AE 3.068 0.608 972
  PP 0.006 -0.028 -5
  SC-PP -0.005 0.046 0
C2 AE 2.332 0.281 1933
  PP -0.001 0.173 7
  SC-PP -0.009 0.098 6
N2 AE 2.011 4.972 2736
  PP -0.003 0.381 24
  SC-PP -0.008 0.170 13
O2 AE 2.184 1.441 1981
  PP -0.017 0.365 66
  SC-PP -0.002 0.025 -10
F2 AE 2.496 -1.607 1283
  PP -0.034 -0.097 63
  SC-PP 0.002 -0.034 -2
LiH AE 3.037 1.483 1427
  PP 0.028 -0.036 -37
  SC-PP 0.000 -0.022 -34
FH AE 1.694 4.203 4501
  PP 0.001 -0.014 -16
  SC-PP 0.000 -0.040 -22
CO AE 2.080 7.530 2444
  PP -0.002 0.262 16
  SC-PP -0.009 0.162 12
average PP 0.011 0.150 26
deviation SC-PP 0.005 0.067 11

Even for diatomic molecules, for which the effect of the spurious long-range component in the unscreened PP is minimal, its impact is not negligible. It is, however, much more pronounced for solids, as Table 10 shows.

Table 10: Equilibrium lattice constant a, cohesive energy Ecoh and bulk modulus B of fcc aluminum obtained from exact exchange (without correlation functional): Standard versus selfconsistent (SC) PP and AE-LDA [18] results (experimental values taken from [19]).
mode a Ecoh B
   [Bohr]   [eV/atom]   [GPa] 
PP 7.10 3.98 135
SC-PP 7.79 1.37 71
AE-LDA 7.52   84
Expt. 7.65 3.39 77

Lacking AE results with the exact exchange, AE-LDA and experimental values are given for comparison. It is obvious that the self-consistent damping of the PPs leads to much more realistic data for the lattice constant and the bulk modulus. On the other hand, the PP calculations only utilize the exact exchange, while correlation is neglected. Consequently, an important contribution to the cohesive energy is missing. The remaining discrepancy between the self-consistent PP results and the experimental data should thus not be misinterpreted as a failure of the SC-PP.

References

  1. E. G. Moroni, G. Kresse, J. Hafner, and J. Furthmüller, Phys. Rev. B 56, 15629 (1997).
  2. N. Chetty, M. Weinert, T. S. Rahman, and J. W. Davenport, Phys. Rev. B 52, 6313 (1995).
  3. M. Krajci, J. Hafner, and M. Mihalkovic, Phys. Rev. B 56, 3072 (1997); M. Krajci and J. Hafner, Phys. Rev. B 58, 5378 (1998).
  4. C. G. Morgan, P. Kratzer, M. Scheffler, Phys. Rev. Lett. 82, 4886 (1999).
  5. R. N. Barnett and U. Landman, Phys. Rev. B 48, 2081 (1993).
  6. M. Eichinger, P. Tavan, J. Hutter, and M. Parinello, J. Chem. Phys. 110, 10452 (1999); D. Marx, M. E. Tuckerman, J. Hutter, and M. Parinello, Nature 397, 601 (1999).
  7. P. M. Morse and H. Feshbach, Methods of Theoretical Physics (McGraw-Hill, New York, 1953).
  8. D. R. Hamann, M. Schlüter, and C. Chiang, Phys. Rev. Lett. 43, 1494 (1979); G.B. Bachelet, D.R. Hamann, and M. Schlüter, Phys. Rev. B 26, 4199 (1982).
  9. G. P. Kerker, J. Phys. C 13, L189 (1980).
  10. D. Vanderbilt, Phys. Rev. B 32, 8412 (1985).
  11. D. R. Hamann, Phys. Rev. B 40, 2980 (1989).
  12. N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
  13. S. G. Louie, S. Froyen, and M. L. Cohen, Phys. Rev. B 26, 1738 (1982).
  14. E. Engel, in: Relativistic Electronic Structure Theory, Part 1. Fundamentals, edited by P. Schwerdtfeger (Elsevier, Amsterdam, 2002), p.524-624.
  15. E. Engel, A. Höck, and S. Varga, Phys. Rev. B 63, 125121 (2001).
  16. S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980).
  17. E. Engel, A. Höck, R. N. Schmid, R. M. Dreizler, and N. Chetty, Phys. Rev. B 64, 125111 (2001).
  18. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Phys. Rev. B 46, 6671 (1992); C. S. Wang, B. M. Klein, and H. Krakauer, Phys. Rev. Lett. 54, 1852 (1985); C. Filippi, D. J. Singh, and C. J. Umrigar, Phys. Rev. B 50, 14947 (1994).
  19. M. Fuchs, M. Bockstedte, E. Pehlke, and M. Scheffler, Phys. Rev. B 57, 2134 (1998).
  20. E. Engel, A. Facco Bonetti, S. Keller, I. Andrejkovics, and R. M. Dreizler, Phys. Rev. A 58, 964 (1998).
  21. D. M. Bylander and L. Kleinman, Phys. Rev. Lett. 74, 3660 (1995); Phys. Rev. B 52, 14566 (1995); Phys. Rev. B 54, 7891 (1996); Phys. Rev. B 55, 9432 (1997).


Home Top

Last modified: April 19, 2004

Disclaimer

ee