Quantum Chemistry of Solids / 21-Basis Sets and Pseudopotentials in Periodic LCAO Calculations

8.1 Basis Sets in the Electron-structure Calculations of Crystals

291

consists of 3s2p1d, cc-pVTZ would be 4s3p2d1f , cc-pVQZ would be 5s4p3d2f 1g (see Table 8.3). The advantage of these is that increasing basis-set size is systematic. In fact, if energies are calculated with the cc-pVDZ, cc-pVTZ and cc-pVQZ basis sets, then there is an analytical function that can be used to extrapolate to cc-pVmZ (the complete basis-set limit). Unfortunately, cc-pVQZ calculations are not commonly doable. Moreover, cc-VXZ functions are more di cult to integrate than the Gaussian functions and these calculations take longer even for comparably sized basis sets. Dunning basis sets already contain polarization functions, but they can be augmented with di use functions.

The EMSL (the environmental molecular sciences laboratory) [463] library supplies a wide selection of atomic basis sets optimized for molecules. This library allows extraction of Gaussian basis sets, and any related e ective core potentials (the e ective core potentials are considered in the next sections), from the Molecular Science Research Center’s Basis Set Library. A user may request the basis set be formatted appropriately for a wide variety of popular molecular electronic-structure packages. In addition to the exponents and contraction coe cients that define the basis set, the user can obtain descriptive data that include the overall philosophy behind the basis literature citations, the angular momentum composition of the basis and many other pieces of information.

Molecular AO basis sets can formally be used in periodic LCAO calculations but their adequacy must be carefully checked [23]. We discuss in the next section the modification of molecular basis sets to the study of crystals.

8.1.3 Molecular Basis Sets Adaptation for Periodic Systems

In the periodic systems the basis sets are chosen in such a way that they satisfy the Bloch theorem. Let a finite number mA of contracted GTFs be attributed to the atom A with coordinate rA in the reference unit cell. The same GTFs are then formally associated with all translationally equivalent atoms in the crystal occupying positions

rA + an (an is the direct lattice translation vector). For the crystal main region

of N primitive unit cells there are N mA Gaussian-type Bloch functions (GTBF)

A

constructed according to
1
χki(r) = √N an		exp(ikan)χi(r − rA − an)	(8.9)

The total number M of GTBFs equals M = mA where the summation is made

A

over atoms of the reference cell. Thus, in solids the basis-set functions are modulated over the infinite lattice: any attempt to use large uncontracted molecular or atomic basis sets, with very di use functions can result in the wasting of computational resources [23]. Therefore exponents and contraction coe cients in the molecular and periodic systems are generally rather di erent, and with some exceptions, such as molecular crystals and certain covalent solids, molecular basis sets are not directly transferable to the study of crystalline solids. As in the case of molecules, basis functions in solids are grouped into shells. In general, a shell contains all functions characterized by the same n and l quantum numbers (e.g. all the di erent d-functions in a 3d shell); this allows the partitioning of the total charge density into “shell charge

292 8 Basis Sets and Pseudopotentials in Periodic LCAO Calculations

distributions” and is useful in the selection of bielectronic integrals and in the evaluation of long-range interactions. A feature of the contraction schemes originally used in the basis sets of the Pople type (and often useful in calculations with the CRYSTAL code [23]) is the additional grouping of basis functions with only the same principle quantum number into shells; e.g. a 2sp shell, in which both 2s and 2p functions have the same set of exponents αj but di erent contraction coe cients dj . For STF this is illustrated by Table 8.1. This procedure reduces the number of auxiliary functions to be calculated in the evaluation of electron integrals. Note that if the basis set is restricted to s, p and d basis functions only sp shells may be formed in this way. In certain circumstances it may actually represent an important constraint on the form of the basis functions. For relatively small calculations where the time and storage limitations are not an important factor, some consideration should be given to describing the s and p functions with separate sets of exponents. Most standard molecular codes use what is known as a segmented contraction scheme, in which the transformation from the larger primitive set to the smaller contracted set is restricted in such a way that each Gaussian primitive gj (r) contributes to exactly one contracted GTF. In contrast, the general contraction scheme makes no such assumptions, and allows each Gaussian primitive to contribute to several contracted GTFs. A considerable advantage of the general scheme is that the contracted GTFs reproduce exactly the desired combinations of primitive functions. For example, if an atomic SCF calculation is used to define the contraction coe cients in a general contraction, the resulting minimal basis will reproduce the SCF energy obtained in the primitive basis. This is not the case with segmented contractions. There are other advantages with a general contraction: for example, it is possible to contract inner-shell orbitals to single functions with no error in the atomic energy, making calculations on heavy elements much easier. Another advantage is a conceptual one: using a general contraction, it is possible to perform calculations in which the one-particle space is a set of atomic orbitals, a true LCAO scheme, rather than being a segmented grouping of a somewhat arbitrary expansion basis. The MOs or COs can then be analyzed very simply, just as for the original qualitative MO LCAO or CO LCAO approach, but in terms of “exact AOs” rather than relatively crude approximations to them.

We describe here in more detail some computational aspects of GTFs use in the CRYSTAL code [23], following [458]. The important reason for the usefulness of a Gaussian basis set is embodied in the Gaussian product theorem (GPT), which in its simplest form states that the product of two primitive Gaussian functions with exponents α and β, located at centers A and B, is itself a primitive Gaussian with exponent γ = α + β, multiplied by a constant factor F , located at a point C along

the line segment A − B, where C =	αA+βB	and F = exp −	αβ	(A − B)2 .
	γ		γ

The product of two polynomial GTF, of degree µ and ν and located at points A and B is therefore another polynomial GTF located at C of degree µ + ν in xC , yC and zC , which can be expressed as a short expansion of one-center Gaussians:

			µ+ν
			i
	χAx (x)χBx (x) =		Ciµ+ν ϕCi(x − xC )			(8.10)
			=0
i	2	) and xc =		αxA +βxB	.
where ϕCi(x) = x exp(−αp(x − xc)				α+β

8.1 Basis Sets in the Electron-structure Calculations of Crystals

293

The product of two Gaussians that are functions of the coordinates of the same electron is referred to as an overlap distribution, and all the integrals that must be calculated involve at least one such overlap distribution. The most important consequence of the GPT is that all four-center two-electron integrals can be expressed in terms of two-center quantities. CRYSTAL code actually uses a common and more e cient approach for the evaluation of integrals over Gaussian basis functions, in which Hermite Gaussian functions (HGFs) are used instead of the usual Cartesian Gaussians in the re-expansion (8.10). Hermite Gaussians are defined as derivatives of an s-type Gaussian:

Λ(ξ) = Hi(ξ) exp −αp(x − xp)2		di
	= (−1)i		exp(−ξ2)	(8.11)
		dξi

where Hi(ξ) is a polynomial of order i, and ξ = αp1/2(x − xp). The set of HGFs spans the same space as the expansion functions in (8.10) and as a consequence they can be used for expanding the basis function products:

	lA +lB
	i
χAx(x)χBx(x) =	CilA +lB Λi(ξ)	(8.12)
	=0

where the expansion coe cients must now be redefined. Because of the natural relations between Hermite polynomials and Gaussians, the necessary two-center integrals can be evaluated with very high e ciency. Basis functions with higher quantum numbers can be generated through repeated di erentiation of an s-type Gaussian. Even though the four-center bielectronic integrals can be written in terms of two-center quantities, the cost of evaluating them still scales nominally as N 4, where N is the number of functions in the expansion. This scaling must be reduced in order to treat large systems. One way of doing this that is used in CRYSTAL is the method of prescreening where, rather than attempting to calculate the integrals more e ciently, one seeks, where possible, to avoid their evaluation altogether. The expression for an integral over primitive Gaussians can be formally written as

(ab|cd) = SabScdTabcd

(8.13)

where Sab is a radial overlap between functions χa and χb, and Tabcd is a slowly varying angular factor. In many situations the product SabScd thus constitutes a good estimate of the magnitude of the integral whose product is used as an estimate in screening out small integrals. In order to estimate these overlaps quickly, a single, normalized s-type Gaussian called an adjoined Gaussian is associated with each shell, whose exponent α is the smallest of the exponents in the shell contraction. This function thus reproduces approximately the absolute value of the corresponding AOs at intermediate and long range. The adjoined Gaussian is used in fast algorithms for estimating overlaps on the basis of which integrals are either evaluated exactly, approximately, or not at all. The level of approximation is user-definable through a set of tolerances given in the input. Such algorithms, and a consideration of the crystalline symmetry, mean that the integrals part of CRYSTAL code scales at between N and N 2, depending on the size of the system. The most unpleasant scaling in this code is thus the SCF part that, since it involves diagonalization of the Fock matrix, scales as approximately N 3. It requires

294 8 Basis Sets and Pseudopotentials in Periodic LCAO Calculations

a special consideration of the role of di use basis-set functions in crystalline systems. Very di use functions can yield numerical instabilities and risk of linear-dependent catastrophes [458]. Furthermore, due to the truncation criteria of the infinite sums, based on the overlap, decreasing the exponents of the primitive Gaussians the number of integrals to be calculated increases very rapidly. Too extended basis sets are not needed in periodic calculations because the complete basis-set limit is reached quicker than in molecular calculations. Furthermore, the risk of linear dependence problems increases. The choice of the AO basis set is a compromise between accuracy and costs. As the accuracy must be the main goal of ab-initio calculations, the good-quality basis sets should always be used in spite of their computational cost to avoid producing meaningless numbers.

The choice of the basis set (BS) is one of the critical points, due to the large variety of chemical bonding that can be found in a periodic system. For example, carbon can be involved in covalent bonds (polyacetylene, diamond) as well as in strongly ionic situations (Be2C, where the Mulliken charge of carbon is close to –4). Some general principles of the basis-set choice for periodic systems are formulated [23, 458] and we briefly reproduce them here. Much experience gained in the molecular computational chemistry can be used in the selection of basis sets for studies of crystalline solids. However, the molecular quantum chemists do not, in general, optimize basis sets by varying the exponents or contraction coe cients to minimize the energy. Rather, there is a hierarchy of basis sets with perceived qualities, and for a di cult problem where accuracy is important, one would use a “good-quality” standard basis set from a library without modification. In crystalline systems by contrast, basis-set optimization is usually necessary, essentially for two reasons. First, there is a much larger variety of binding than in molecules and basis sets are thus less transferable. Secondly, hierarchical libraries of basis sets comparable to those available for molecules do not really exist for solids. For certain types of compounds, such as molecular crystals or many covalent materials, the molecular sets can sometimes be used largely unmodified, but this has to be done carefully. However, for strongly ionic crystals and metals the basis sets, particularly the valence states, need to be redefined completely. In essentially all cases, the core states may be described using the solutions of atomic calculations, as even in the presence of strong crystal fields the core states are barely perturbed and may be described by the linear variation parameters in the SCF calculation.

Redefining basis sets in this way is obviously time consuming and even more obviously rather boring, and so over recent years various workers involved with the CRYSTAL code have contributed to an e ort to develop libraries of basis sets to be made available on the Internet. The URL of the o cial site is: www.crystal.unito.it /Basis sets /Ptable.html. We recommend to the reader to use the information given on this site. The site shows a periodic table. Clicking on the symbol for the required element will reveal a text file containing various di erent basis sets that may have been used in di erent materials containing that element type. Accompanying each basis set is a list of authors, a list of materials where the set has been used, references to publications and hints on optimization where relevant. This table is obviously not complete (in particular lanthanides and actinides are not practically represented as the calculations with f -electrons are not included even in the CRYSTAL06 code). Additional information can be found on the site of Cambridge basis set library: www.tcm.phy.cam.ac.uk/ mdt26/crystal.html.

8.1 Basis Sets in the Electron-structure Calculations of Crystals

295

We follow [458] in discussion of optimization strategies for the basis sets given in the libraries mentioned and also the adaptation of molecular bases for various types of solids. First, a number of general principles are given that should be taken into account when choosing a basis set for a periodic system.

The di use basis functions are used in atoms and molecules to describe the tails of the wavefunction, which are poorly described by the long-range decay of the Gaussian function. In periodic systems the cost of HF/DFT calculations can more essentially increase when the di use basis functions are included (in the silicon and diamond crystals, for example, the number of bielectronic integrals, can increase by a factor of 10 simply by changing the exponent of the most di use single Gaussian from 0.168 to 0.078 (Si) and from 0.296 to 0.176 (C), [458]. Fortunately, in crystalline compounds in contrast to molecules, particular in nonmetallic systems, the large overlap between neighbors in all directions drastically reduces the contribution of low-exponent Gaussians to the wavefunction. This has the consequence that a small “split-valence” basis set such as 6-21G is closer to the Hartree–Fock limit in crystals than in molecules.

The number of primitives is an important feature of the basis set. A typical basis set for all electron calculations includes “core functions” with higher exponents and a relatively large number of primitives – these will have a large weight in the expansion of the core states. The “valence functions” with a large weight in the outer orbitals will have lower exponents and contractions of only a very few primitives. It is possible to get away with putting a lot of primitives in the core since core states have very little overlap with neighboring atoms. The use of many primitives in the valence shells would add significantly to the cost of a calculation, but in many cases it is necessary.

There are several ways to improve a basis set [458]: 1) reoptimize the more diffuse exponents (and contraction coe cients if necessary; 2) decontract i.e. convert the more di use contractions into single Gaussian primitives; 3) convert sp functions into separate s and p functions; 4) add polarization functions if not already present; 5) add more primitives (watch out for linear dependence problems); 6) use a better starting point for the basis set. Optimization in this sense means varying an appropriate subset of the basis-set parameters until the energy is minimized. In principle, this is a reasonably complex multidimensional minimization, but there are various standard shell scripts available. Two of them can be downloaded from the Internet sites: BILLY code www.tcm.phy.cam.ac.uk/ mdt26/downloads/billy.tar.gz) and LoptCG code (www.crystal.unito.it/LoptCG/LoptCG.html).

The adequacy of the starting molecular basis sets depends on the type of crystalline compound.

To describe covalent crystals the small molecular split-valence basis sets can be used with confidence and essentially without modification. It is enough to reoptimize the exponent of the most di use shell, which produces a slightly improved basis, while reducing the cost of the calculation. That said, 6-21G* is not really all that good and a larger better basis set with more variation freedom is quite easy to make for these cases (see web libraries).

For fully-ionic crystals (like alkali halides, LiH or MgO) with an almost completely empty cation valence shell it often proves convenient to use a basis set containing only “core” functions plus an additional sp shell with a relatively high exponent. It is usually di cult, and often impossible, to optimize the exponents of functions that only have appreciable weight in almost empty orbitals. Anions present a di erent

296 8 Basis Sets and Pseudopotentials in Periodic LCAO Calculations

problem. Reference to isolated ion solutions is possible only for halides, because in such cases the ions are stable even at the Hartree–Fock level. For other anions, which are stabilized by the crystalline field, the basis set must be redesigned with reference to the crystalline environment. Consider, for example, the optimization of the oxygen basis set in Li2O [458]. The di culty is to allow the valence distribution to relax in the presence of two more electrons. We can begin from a standard STO-6G basis set, i.e. six contracted primitive Gaussians for the 1s shell, and six more to describe the 2sp shell. First , two more Gaussians were introduced into the 1s contraction, in order to improve the virial coe cient and total energy. The two outer Gaussians of the valence sp shell were then removed from the contraction and allowed to vary independently. The exponents of the two outer independent Gaussians and the coe cients of the four contracted ones were optimized in Li2O. The best outer exponents of the ion were found to be 0.45 and 0.15 and are therefore considerably more di use than the neutral isolated atom, where the best exponents are 0.54 and 0.24. The rest of the oxygen valence shell is unchanged with respect to the atomic situation. The introduction of d functions in the oxygen basis set gives only a minor improvement in the energy, with a population of 0.02 electrons/atom/cell (d functions may be important in the calculation of certain properties, however). Thus, for anions, reoptimization of the most di use valence shells is mandatory when starting from a standard basis set.

The majority of crystals can be classified as semi-ionic (with chemical bonding being intermediate between covalent and purely ionic limits). For such crystals the adequacy of selected basis sets must be carefully tested as is discussed in [458], for example, for semi-ionic compounds SiO2 (quartz) and Al2O3 (corundum). The exponents of the outer shell for the two cations (Si and Al) used in molecular calculations prove to be too di use. For the Si atom in quartz, reoptimization in the bulk gives

α= 0.15 (instead of the molecular value 0.09). Corundum is more ionic than quartz, and in this case it is better to eliminate the most di use valence shell of Al, and to use two Gaussians of the inner valence shells as independent functions (α = 0.94 and

α= 0.3, respectively).

For metals very di use Gaussians are required to reproduce the nearly uniform density so that it was often stated that the plane waves are a more appropriate basis for these systems. It is generally impossible to optimize atomic-like basis set in Hartree–Fock calculations of metallic systems. However, Gaussian DFT studies indicate that GTFs are able to provide a reliable and e cient description of simple metallic systems (see, for example, metallic lithium calculations [464]). Thus, the e ects of basis set and Hamiltonian were separated.

UHF and hybrid DFT calculations for strongly correlated transition-element magnetic compounds require the reasonably good basis sets for the transition elements. These are not that widely available even to molecular quantum chemists since most of the e ort in developing molecular GTF basis sets has been for firstand second-row atoms. One reason for this may be that molecules containing transition metal atoms tend to be very badly described at the Hartree–Fock level. Molecular bonds tend to have a fairly high degree of “covalency” and the existence of partially-occupied d states leads to a great many nearly degenerate levels, and thus to a large “static correlation” (i.e. the weight of the HF determinant in a CI expansion would be small, and a multideterminant treatment is more appropriate). Basis sets to describe the correlation using quantum chemistry correlated wavefunction techniques need to be

8.1 Basis Sets in the Electron-structure Calculations of Crystals

297

much richer than those for systems well described at the HF level since they need to treat all of the unoccupied levels. It may seem surprising that single-determinant HF could be so successful in periodic crystalline magnetic insulators containing transition elements, but this is an important characteristic of these ionic materials. The highly symmetric environment and long-range Coulomb forces tend to separate the orbitals into well-defined subsets with a significant gap between occupied and unoccupied states. Hence, the ground state of NiO (for example) is rather well described by a single determinant. In this sense, a strongly correlated magnetic insulator is in many ways a “simpler system” than many molecules. The success of UHF calculations in these materials (and also hybrid DFT schemes) is now well known. We consider some examples of such calculations in Chap. 9. Molecular GTF basis sets for transition elements have been reoptimized in the solid state and are available on the Torino and Cambridge Gaussian basis set library web sites referred to earlier.

In the final part of this section we shall discuss a rather serious problem, associated with Gaussian basis sets – basis-set superposition error (BSSE). A common response to this problem is to ignore it, since it will go away in the limit of a complete basis but to achieve this limit one needs calculations that are seldom performed. The problem of BSSE is a simple one: in a system comprising interacting fragments A and B, the fact that in practice the basis sets on A and B are incomplete means that the fragment energy of A will necessarily be improved by the basis functions on B, irrespective of whether there is any genuine binding interaction in the compound system or not. The improvement in the fragment energies will lower the energy of the combined system giving a spurious increase in the binding energy. It is often stated that BSSE is an e ect that one needs to worried about only in calculations on very weakly interacting systems. This is not really true [458]. BSSE is an ever-present phenomenon and accurate calculations should always include an investigation of BSSE. Examples of areas in which one should be particularly worried are the study of the binding energy of molecules adsorbed on surfaces (see, for example [465] for an interesting discussion) or the calculation of defect-formation energies.

The approach most commonly taken to estimate the e ect of BSSE is the counterpoise correction of Boys [466]: the separated fragment energies are computed not in the individual fragment basis sets, but in the total basis set for the system including “ghost functions” for the fragment that is not present. These energies are then used to define a counterpoise-corrected (CPC) interaction energy, which by comparison with perturbation theory, has been shown to converge to the BSSE-free correct value [467].

Linear dependencies of Gaussian-type orbital basis sets employed in the framework of the HF SCF method for periodic structures, which occur when di use basis functions are included in a basis set in an uncontrolled manner, were investigated [468].The basis sets constructed avoid numerical linear dependences and were optimized for a number of periodic structures. The numerical AO basis sets for solids were generated in [469] by confining atoms within spheres and smoothing the orbitals so that the first and second derivatives go to zero at the boundary. This forms small atomic-like basis sets that can be applied to solid-state problems and are e cient for treating large systems.

We considered the basis sets for all-electron calculations in which all the core electrons are involved explicitly. Both in molecular and solid-state quantum chemistry the core states can be treated implicitly as creating an e ective core potential (pseudopo-

298 8 Basis Sets and Pseudopotentials in Periodic LCAO Calculations

tential). The pseudopotential approximation becomes the most e cient for crystalline compounds of heavy elements. Simultaneously, the relativistic e ects on the electronic structure play an important role. In the next sections we consider nonrelativistic and relativistic pseudopotentials used in modern LCAO calculations of periodic systems. The choice of the corresponding valence basis sets is also discussed.

8.2 Nonrelativistic E ective Core Potentials and Valence Basis

Sets

8.2.1 E ective Core Potentials: Theoretical Grounds

It is well established from chemical experience that most chemical properties of molecules and solids are determined by the valence electrons of the constituent atoms. The core states are weakly a ected by changes in chemical bonding. The e ect of core electrons is principally to shield the nuclear charges and to provide an e ective potential for the valence electrons. The main reason for the limited role of the core electrons is the spatial separation of the core and valence shells that originates from the comparatively strong binding of the core electrons to the nucleus. This e ect is illustrated in Fig. 8.2, [470], where the radial orbitals of the oxygen atom are plotted. The spatial separation of the 1s core state from the valence 2s, 2p states and

( )r

nl rR

2.0

				Oxygen
1.5					1s
					2s
					2p
1.0					3s
0.5
0.0
-0.5
0		1	2 r, a.u. 3			4	5

Fig. 8.2. Radial orbitals of oxygen atom in the ground state, [470]

unoccupied (excited) 3s state is seen. The HF eigenvalues of 1s, 2s and 2p states (for the ground state of free oxygen atom) are –20.6687 a.u., –1.2443 a.u., –0.6319 a.u., respectively [65], i.e. the valence and core eigenvalues di er by an order of magnitude. The precise decomposition of electrons into core and valence becomes more di cult

8.2 Nonrelativistic E ective Core Potentials and Valence Basis Sets

299

when 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 and the latter are not well separated from 3s, 3p states (for example, in the case of Ti atom, HF one-electron energies are –0.2208 a.u. (4s), –0.4407 a.u (3d), –1.7952 a.u (3p), –2.8735 a.u (3s), [65]). As a result, the overlap of the 3d-orbitals (and even 3s, 3p orbitals) with the orbitals from neighboring atoms is too large to be ignored. In this situation large core and small core e ective potentials can be introduced; the former di ers from the latter by exclusion of the 3s, 3p electrons from the core.

The idea behind e ective core potentials (ECPs), also called pseudopotentials (PP), is to treat the core electrons as creating e ective averaged potentials rather than actual particles. E ective core potentials are based on the frozen-core approximation and serve to represent the potential generated by core electrons, also incorporating relativistic e ects. ECP application can introduce significant computational e ciencies as it allows formulation of a theoretical method for dealing only with the valence electrons, while retaining the accuracy of all-electron methods. Fock, Veselow and Petrashen were the first to formulate the approach for treating a subset of the electrons, such as the valence electrons, in the field of the remaining electrons. They considered the special case of atoms with two valence electrons outside a closed shell [101].

As we already noted, the pseudopotentials are essentially mandatory in planewave calculations of solids since the core orbitals have very sharp features in the region close to the nucleus and too many plane waves would be required to expand them if they were included. In atomic-like basis-set calculations pseudopotentials are formally not mandatory and have di erent characteristics from those designed for plane waves since localized basis functions actually have necessary sharp features in the core region. Nevertheless, ECP methods are also used in LCAO calculations of molecules and solids since the di culty of the standard LCAO methods rises rapidly with the number of electrons. If the CPU time in LCAO calculations were dominated by the integrals calculation, PP introduction would not give very much since the number of integrals is controlled by more di use functions that overlap strongly with neighboring atoms. These di use functions are introduced mainly to describe the change of atomic valence states. However, the use of ECP will decrease the number of coe cients in the one-electron wavefunctions and might give significant savings in the SCF part. It is also quite easy to incorporate relativistic e ects into pseudopotentials, which is increasingly important for heavy atoms. All electron relativistic calculations are very expensive and in many cases are practically di cult.

The parameters and the underlying basis set of so-called energy-consistent ECPs can be adjusted in accordance with representative experimental data, not only for the ground state, but also for excited states, electron a nities, ionization potentials and so on. Being of semiempirical origin, they can perform remarkably well for a given system, but their transferability to other environments can be poorer than that of shape-consistent ECPs. So-called shape-consistent ECPs (known in computational condensed matter as the norm conserving ECPs) are rather easy to derive and contain no adjustable parameters, i.e. these are ab-initio ECPs. The construction of shapeconsistent ECPs for molecules is based on the original proposal of Christiansen et al. [471] where shape consistency was introduced. Simultaneously, norm conserving ECP were introduced in computational condensed-matter physics by Hamann et al.

300 8 Basis Sets and Pseudopotentials in Periodic LCAO Calculations

[472]. This example demonstrates practically independent development of ECPs for molecules and crystals when for the same ECP property di erent terms are used. The work by Phillips and Kleinman (PK) [473] is an important step in the ECP applications for solids. PK developed the pseudopotential formalism as a rigorous formulation of the earlier “empirical potential” approach. They showed that ECP that has the plane-wave pseudo wavefunctions as its eigenstates could be derived from the all-electron potential and the core-state wavefunctions and energies. Thus a nonempirical approach to finding ECP was introduced.

The PK pseudopotential shortcomings are well known [474]: it depends explicitly on the one-electron eigenvalue and outside the core region the normalized pseudo orbital (PO) is proportional but not equal to the true orbital. Typically, generation of a pseudopotential proceeds as follows. First, a cuto distance rc for the core is chosen (in the case of our example, the oxygen atom, rc can be taken as 0.5 a.u., see Fig. 8.2). In a PP approach all radial Rnl(r) 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 of such a nodeless Rnl(r) that, nevertheless, agrees with the corresponding true AO in the relevant part of space, is shown in Fig. 8.3a.

It is, however, impossible to produce two (or more) nodeless orbitals in the same energy range with only a single spherical PP, as for fixed PP only the angular momentum term can generate di erences. Since PP replaces the potential of a nucleus and the core electrons, it is spherically symmetric and each angular momentum l can be treated separately, which leads to nonlocal l-dependent PP Vl(r). Consequently, the total atomic PP usually consists of several components, one for each angular momentum present in the valence space. The PP dependence upon l means that, in general, PP is a nonlocal operator, that can be written in semilocal (SL) form

ˆ P S	=		(θ, ϕ)|	(8.14)
V		|Ylm(θ, ϕ) Vl(r) Ylm

lm

where Ylm(θ, ϕ) are spherical harmonics and Vl(r) is the pseudopotential for the lth angular-momentum component. This is termed semilocal because it is nonlocal in the angular variables, but local in the radial variable: when operating on a function

ˆ P S

f (r, θ , ϕ ), V has the e ect [10]

a

b

( )r

nl rR

1.0
0.8			Oxygen 2s
				pseudo
0.6
				all-electron
0.4
0.2
0.0
-0.2
-0.4
0	1	2		3	4
		r, a.u.

( )r

nl rR

1.0

Oxygen 2p

0.8					pseudo
0.6					all-electron
0.4
0.2
0.0
	0	1	2		3	4
			r, a.u.

Fig. 8.3. Radial pseudoversus all-electron orbitals of oxygen: (a) 2s (b) 2p [470]