Quantum Chemistry of Solids / 17-Kohn-Sham LCAO Method for Periodic Systems

7.1 Foundations of the Density-functional Theory

241

made the LDA whole physically meaningful [337]. Higher-order corrections, on the other hand, are exceedingly di cult to calculate, and little is known about them.

It was a major breakthrough when it was realized, in the early 1980s, that instead of power-series-like systematic gradient expansions one could experiment with more general functions of ρ(r) and ρ(r), which need not proceed order by order. Such functionals, of the general form

EXCGGA[ρα, ρβ ] = d3rf (ρα, ρβ , ρα, ρβ ) (7.19)

have become known as generalized-gradient approximations (GGAs), [355]. GGA functionals are the workhorses of current density-functional theory. Di erent GGAs di er in the choice of the function f (ρ, ρ). Note that this makes di erent GGAs much more di erent from each other than the di erent parametrizations of the LDA: essentially there is only one correct expression for εhomXC (ρ) and the various parametrizations of the LDA are merely di erent ways of writing it [330]. On the other hand, depending on the method of construction employed for obtaining f (ρ, ρ) one can obtain very di erent GGAs. In particular, GGAs used in molecular quantum chemistry typically proceed by fitting parameters to test sets of selected molecules. On the other hand, the GGAs used in physics tend to emphasize exact constraints on the density-functional for the exchange-correlation energy [357]. In this approach, the density-functional approximations are assigned to various rungs according to the number and kind of their local ingredients [357].

The lowest rung is the local spin-density approximation LSDA, the second rung is the generalized gradient approximation GGA (semilocal functionals). More accurate functionals of higher rungs can be called “beyond-GGA functionals” or nonlocal functionals. Third-rung Meta-GGA (MGGA) functionals depend, in addition to the density and its derivatives, also on the kinetic-energy density [362]. Fourth-rung hybrid functionals mix a fraction of the Hartree–Fock exchange into the DFT exchange functional (see Sect. 4.2). MGGA and hybrid functionals can be called orbital functionals because they are represented not only in terms of the electron density, but also contain parts represented in single-particle Kohn–Sham orbitals ϕi(r). In MGGA functionals kinetic energy density τ (r) = 1/2 i|ϕi(r)|2 is included. In hybrid functionals the nonlocal HF exchange energy is also orbital functional. Still another type of orbital functional is the self-interaction correction (SIC) (see Sect. 7.4). Higher-rung density-functionals are increasingly more complex.

The semiempirical functionals are fitted to selected data from experiment or from the ab-initio calculations. The higher the rung of the functional the larger is the number of parameters (functionals with as many as 21 fit parameters are popular in chemistry). Is DFT ab-initio or semiempirical? As was suggested in [357] it can fall in between as a nonempirical theory when the functionals are constructed without empirical fitting.

The best nonempirical functional for a given rung is constructed to satisfy as many exact theoretical constraints as possible while providing satisfactory numerical predictions for real systems [357]. Once a rung has been selected, there remains little choice about which constraints to satisfy (but greater freedom in how to satisfy them), [357]. Accuracy is expected to increase up the ladder of rungs as additional local ingredients enable the satisfaction of additional constraints. A short summary of

242 7 Kohn–Sham LCAO Method for Periodic Systems

exact constraints on EXC [ρ] can be found in [357]. In this paper some useful recommendations for users of DFT are also given. Users should not randomly mix and match functionals, but should use exchange and correlation pieces designed to work together, with their designer-recommended local parts. They should not shop indiscriminately for the functional that “works best”. Users should always say which functional they used, with its proper name and literature reference, and why they chose it. Statements like “we used density-functional theory” or “we used the generalized gradient approximation” are almost useless to a reader or listener who wants to reproduce the results.

Nowadays, the most popular (and most reliable) GGA functionals are PBE (denoting the functional proposed in 1996 by Perdew, Burke and Ernzerhof [356] in physics, and BLYP (denoting the combination of Becke’s 1988 exchange functional [358] with the 1988 correlation functional of Lee, Yang and Parr [359]) in chemistry. PWGGA denotes the GGA functional, suggested by Perdew and Wang [354, 355] and should be allowed to retire gradually. Many other GGA-type functionals are also available, and new ones continue to appear. The new meta-GGA (MGGA) functional TPSS (Tao–Perdew–Staroverov–Scuseria) [360] supersedes an older one the PKZB (Perdew– Kurth–Zupan–Blaha) functional [363]. The known functionals are modified as has been done recently for the PBE functional to improve its accuracy for thermodynamic and electronic properties of molecules [361, 362].

We close this subsection with more detailed consideration of the PBE functional [363]. It has two nonempirical derivations based on exact properties of the XC hole [364] and energy [356]. In discussion of the PBE functional we follow [365]. The

exchange-correlation (XC) energy in the LSD approximation is given by

EXCLSD[ρα, ρβ ] =

d3rρ(r)[εX (ρ(r))f (ζ, r) + εC (rs(r), ζ(r))]

(7.20)

where ζ =

ρα −ρβ

is the relative spin polarization and f (ζ) is given by

ρα +ρβ

f (ζ) =

1

(1 + ζ)4/3 + (1 − ζ)4/3

(7.21)

2

The exchange energy per electron εX (ρ) of the unpolarized uniform electron gas

3

1/3

C

of density ρ(r) depends on rs =

4π

ρ

as εX (ρ) = −rs and the correlation energy

ε

per electron in a spin-polarized uniform electron gas depends on r

and ζ [354].

C

P BE

P BE

s

The PBE functional EXC

is a sum of PBE exchange EX

and PBE correlation

EP BE functionals. The exchange PBE functional is written in the form

C

EXP BE (ρ) =

d3rρεX (ρ)FX (s)

(7.22)

with

k

FX (s) = 1 + k −

(7.23)

1 + µs2/k

and k = 0.804, µ = 0.21951. This functional is found from the exact spin-scaling

relation	1		1
EX [ρα, ρβ ] =		EX [2ρα] +		EX [2ρβ ]	(7.24)
	2		2

7.1 Foundations of the Density-functional Theory		243
The PBE correlation energy functional is
ECP BE [ρα, ρβ ] =	d3rρ[εC (rs, ζ) + H(rs, ζ, t)]	(7.25)

where the nonlocal part H(rs, ζ, t) depends on the parameter t including the density

gradient:

H = γΦ3 ln 1 +

β

1 + At2

t2

(7.26)

γ

1 + At2 + A2t4

where β = 0.066725, γ = 0.03191, and

t =

π 1/2 9π

1/6

s

(7.27)

4

Φrs1/2

4

| ρ|

3

9

1/3

s =

=

r

(7.28)

2(3π2)1/3ρ4/3

2 4π

s|

|

Φ(ζ) =

1

(1 + ζ)2/3 + (1 − ζ)2/3

(7.29)

2

A =

β

1

(7.30)

γ

exp[−εC (rs, ζ)/γΦ3] − 1

The PBE GGA correctly reduces to LSD for uniform electron densities, its correlation component recovers the slowly varying (t → 0) and rapidly varying (t → ∞) limits of numerical GGA. Under uniform scaling (ρ(r) → λ3ρ(λr)), the PBE exchange energy scales like λ (as does the exact exchange functional) and the PBE correlation energy correctly scales to a constant as λ → ∞. For small-amplitude density variations around a uniform density, LSD is a very good approximation. The PBE functional recovers this limit. The PBE functional satisfies the Lieb–Oxford bound [366]:

EXP BE [ρα, ρβ ] ≥ EXCP BE [ρα, ρβ ] ≥ 2.273EXLDA[ρ] (7.31)

A useful form to compare GGA functionals is to write

EXCGGA[ρα, ρβ ] ≈ d3rρεX (ρ)FXC (rs, ζ, s)

(7.32)

where the enhancement factor FXC (rs, ζ, s) over local exchange depends upon the local radius rs(rs 1 for core electrons and rs 1 for valence electrons), spin polarization ζ and inhomogeneity s. The s-dependence of FXC is the nonlocality of the GGA. We see that the nonempirical PBE functional best fulfils many of the physical and mathematical requirements of DFT. We will consider some applications of the PBE functional in LCAO calculations of crystals in Chapters 9–11. At the same time, there are some failures of the PBE functional essential for extended systems. For example, the exact exchange-correlation hole in crystals may display a di use long-range tail that is not properly captured by GGA (such a di use hole arises in the calculation of the surface energy of a metal). Neither LSD nor GGA describe correctly the long-range tail of the van der Waals interaction. Corrections of these failures might be possible by using orbital-dependent hybrid functionals considered in the next section.

244 7 Kohn–Sham LCAO Method for Periodic Systems

7.1.5 The Pair Density. Orbital-dependent Exchange-correlation

Functionals

Density-functional theory, even with rather crude approximations such as LDA and GGA, is often better than Hartree–Fock: LDA is remarkably accurate, for instance, for geometries and frequencies, and GGA has also made bond energies quite reliable. Therefore, “the aura of mystery” appeared around DFT (see discussion of this by Baerends and Gritsenko [367]). The simple truth is not that LDA/GGA is particularly good, but that Hartree–Fock is rather poor in the two-electron chemical-bond description. This becomes clear when one considers the statistical two-electron distribution, which is usually cast in terms of the exchange-correlation hole: the decrease in probability to find other electrons in the neighborhood of a reference electron, compared to the (unconditional) one-electron probability distribution [337].

The concept of electron density (7.1), which provides an answer to the question “how likely is it to find one electron of arbitrary spin within a particular volume element while all other electrons may be anywhere” can be extended to the probability of finding not one but a pair of two electrons with spins σ1 and σ2 simultaneously within two volume elements dr1 and dr2, while the remaining N − 2 electrons have arbitrary positions and spins. The quantity that contains this information is the pairelectron density ρ2(x1, x2) which is defined as

ρ2(x1, x2) = N (N − 1)	. . . Ψ (x1, . . . , xN )2dx3 . . . dxN	(7.33)

The pair-electron density actually contains all the information about electron correlation. The pair density is a nonnegative quantity, symmetric in the coordinates and normalized to the total number of nondistinct electron pairs N (N − 1). Obviously, if electrons would not interact, the probability of finding one electron at a particular point of coordinate-spin space would be completely independent of the position and spin of the second electron and it would be possible that both electrons are simultaneously found in the same volume element. In this case, pair density would reduce to a simple product of the individual probabilities, i.e.

	ρ2(x1, x2) =	N − 1	ρ(x1)ρ(x2)	(7.34)
		N

The N (N − 1) factor enters because the particles are identical. From the antisymmetry of the many-electron wavefunction Ψ (x1, . . . , xN ) it follows that ρ2(x1, x2) = −ρ2(x2, x1), i. e. the probability of finding two electrons with the same spin at the same point in space is zero.

This e ect is known as exchange or Fermi correlation and is a direct consequence of the Pauli principle. The Fermi hole is in no way connected to the charge of electrons and applies equally to neutral fermions. This kind of correlation is included in the HF approach due to the antisymmetry of the Slater determinant [337]. The electrostatic repulsion of electrons (the 1/r12 term in the Hamiltonian) prevents the electrons from coming too close to each other and is known as Coulomb correlation. This e ect is independent of the spin and is called simply electron correlation, and is completely neglected in the HF method.

7.1 Foundations of the Density-functional Theory

245

The conditional probability Ω(x1, x2) = ρ(x1, x2)/ρ(x1) defines the probability of finding any electron at position 2 in coordinate-spin space if there is one already known to be at position 1.

The exchange-correlation (xc) hole is defined as the di erence between Ω(x1, x2) and the uncorrelated probability of finding an electron at x2:

		hxc(x1, x2) =	ρ2(x1, x2)		− ρ(x2) = ρ(x2)f (x1, x2)		(7.35)
			ρ(x1)
and contains exactly the charge of one electron						hxc(x1, x2)dx2 = −1, [337]. The
function f (x		, x2) is the correlation factor, being				zero for the completely uncorrelated
	1					/

case.

The exchange-correlation hole describes the change in conditional probability caused by the correction for self-interaction, exchange and Coulomb correlation, compared to the completely uncorrelated situation. It can formally be split into the Fermi

hole, hxσ1=σ2 (r1, r2) and the Coulomb hole hcσ1,σ2 (r1, r2)
hxc(x1, x2) = hxσ1=σ2 (r1, r2) + hcσ1,σ2 (r1, r2)	(7.36)

The exchange hole hx applies only to electrons with the same spin, the correlation hole hc has contributions for electrons of either spin and is the hole resulting from the 1/r12 electrostatic interaction. This separation is convenient but only the total xc hole has a real physical meaning. In the HF method the Fermi hole is accounted for through the use of a single Slater determinant whereas the Coulomb hole is neglected. Like the

total hole, the Fermi hole contains exactly the charge of one electron		hx(r1, r2)dr2		=
−1 and takes care of the self-interaction correction (SIC). The	exchange hole h			is
	/		x

negative everywhere hx(r1, r2) < 0 and its actual shape depends not only on the Fermi correlation factor but also on the density at r2. As a consequence, hx is not spherically symmetric.

The Coulomb hole must be normalized to zero and this result is independent of

/

the positions of electrons with σ =σ: hc(r1, r2)dr2 = 0.

The essential error of Hartree–Fock arises from the fact that in the HF model the exchange-correlation hole in a two-electron bond is not centered around the reference electron, but is too delocalized, having considerable amplitude on both atoms involved in the bond. The LDA and GGA models incorporate, simply by using an electroncentered hole, an important part of the e ect of the interelectron Coulomb repulsion.

To improve the LDA and GGA models the orbital-dependent functionals are introduced. The introduction of orbital dependence, not only density and gradient dependence, into the functionals can be realized in di erent ways [330]. Hybrid HF-DFT functionals are the most popular beyond-GGA functionals. These functionals mix a fraction of HF exchange into the DFT exchange functional and use the DFT correlation part. But why is the exact HF exchange mixed with the approximate DFT exchange part? How can one find the weights of this mixing? May be it would be possible to use the exact HF exchange and rely on approximate functionals only for the part missing in the HF model, i.e. the electron correlation: Exc = ExHF + EcKS (EXX – orbital-dependent exact exchange functional). The HF exchange is calculated using KS orbitals so that EXX energy will be slightly higher than the HF energy. The orbital-dependent exact-exchange (EXX) group of methods has received much attention in the literature. Two advantages of the EXX functional are often mentioned:

246 7 Kohn–Sham LCAO Method for Periodic Systems

(a) the self-interaction correction is incorporated in this functional; (b) it is natural to break up the total problem of electron–electron interaction in the large (due to self-interaction correction) exchange part and the small correlation part. The hope is of course that it will be simpler to find an accurate density-functional for the small correlation energy.

A critique of the EXX method in DFT can be found, for example, in [337,367]. It is clear that the EXX method does not o er any advantage over HF for energies. On the contrary, starting with EXX actually undermines the very basis of the success of DFT, which has been the breaking away from the quantum-chemistry paradigm of dividing up the full electron–electron interaction problem into two steps: first Hartree–Fock (or exchange only) and next the remaining correlation problem. That this is indeed the case is masked by the DFT terminology that still uses the names “exchange” (e.g., Becke exchange, Perdew exchange) and “correlation” (e.g., Perdew–Wang correlation, Lee–Yang–Parr (LYP) correlation). One should be aware that, while correlation is indeed much smaller than exchange in for instance (heavy) atoms, this is not true in the electron-pair bond. Both Hartree–Fock and the use of the exact-exchange functional in Kohn–Sham DFT are very deficient models for the description of chemical bonding. They both use a two-particle probability distribution that is wrong, as can be seen from the shape of the “hole” a reference electron creates around itself: the hole is not centered around the reference electron, but is delocalized over the atoms forming the bond. From this wrong hole shape one can understand the deficiencies of this model, such as too di use orbitals and density, hence wrong kinetic energy and electron–nuclear attraction energy. The true xc hole is substantially more localized around the reference electron. That is why rough localized model holes, like those of GGA, which approximate this total hole, are so successful.

This can be illustrated for the prototype electron-pair bond, the H2 molecule, at the equilibrium H–H distance [337, 367]: in the H2 case HF and EXX are identical. The Hartree–Fock (EXX) model is a poor zero-order approximation for the electronpair bond. It is also extremely size inconsistent, in the sense that when we introduce a second H atom in the universe, remote from a given H atom, the HF/EXX model totally fails to deliver the sum of the energies of two H atoms. However, when we add the Coulomb hole (or correlation hole) to the Fermi (exchange) hole, we obtain a total hole that is localized at the nucleus where the reference electron is located. This will provide the right potential, and therefore the right density and kinetic energy. The localized hole does not require any knowledge of where the other H atom is located, it is the unphysical breaking up of this total hole into a delocalized exchange part and an equally delocalized Coulomb part that complicates the electron–electron interaction treatment and requires knowledge of the position of the second nucleus in order to build accurate holes. Actually, this example illustrates that the division in exchange and correlation holes creates one of the largest problems of DFT. The problem basically is that the total density ρ(r) around each H nucleus is identical at large H–H distances to that in a free H atom. However, the exchange and correlation energies in the stretched H2 are very di erent from those in two isolated H atoms. In particular, the correlation energy is zero in individual H atoms, in contrast to the strong left-right nondynamical correlation in the dissociating H2 molecule. It is virtually impossible to devise functionals that use only local information (local density and derivatives of the density) and still recognize the position of the other H atom and

7.1 Foundations of the Density-functional Theory

247

build correctly the individual exchange and correlation holes, the shape of which is determined by the positions of the nuclei. The LDA and GGA total holes are localized around the reference electron, which is a definite advantage, the EXX approximation can not be regarded to be the next improvement beyond GGA. There are other possibilities of such an improvement: 1) to develop orbital-dependent functionals, which represent exchange and correlation simultaneously; 2) use HF-DFT exchange mixing.

The first possibility was realized recently in [367–370] where the virtual orbitaldependent xc functional ExcBB (BB, Buijse–Baerends) is suggested:

		M	M					ϕ (r )ϕ (r )ϕ (r )ϕ (r )
	1							i 1 i 2 j 1	j	2
ExcBB [{ϕj }, {ϕa}] = −				√wiwj			dr1dr2				(7.37)
	2	i	j					r12

and applied in calculations of hydrogen chains Hn and small molecules. Here wi[ρ] is the xc orbital weight, which governs the involvement of the occupied/virtual orbital ϕi in the xc functional. It is important that the summation in (7.37) is made over M = Nocc + Nv orbitals where Nocc and Nv are the numbers of occupied and virtual KS orbitals. By giving weights to the virtual orbitals one can incorporate the e ect of correlation. In [368] the functional form of wi has been approximated with the Fermi-type distribution

wi =	2	(7.38)
	1 + exp(f (εi − εF ))

where εi are the KS orbital energies and εF is the Fermi-level parameter. We do not consider the ExcBB functional in more detail as its possibilities for crystals are not well studied.

The second possibility of GGA improvement – mixing of HF and DFT exchange – is used in the well-known B3LYP, B3PW and PBE0 hybrid functionals. The performance of hybrid density-functionals (HDF) in solid-state chemistry is examined in a review [371] and many original publications. It is demonstrated that the HDF methods allow calculation of di erent properties of solids in good agreement with the experiment. We return to this point in Sect. 7.2. In Kohn–Sham density-functional theory, the

exchange-correlation energy is rigorously given by
1			ρ(r)	hxc(r, r )
Exc =					d3r d3r	(7.39)
	2			|r − r|
hxc(r, r ) = 1 hxcλ (r, r )dλ						(7.40)

0

known as the “adiabatic connection” or “coupling strength integration”. In this most fundamental of DFT formulas, hxc is an e ective exchange-correlation hole, and λ is a coupling-strength parameter that switches on interelectronic 1/r12 repulsion, subject to a fixed electronic density (achieved, in principle, by suitably adjusting the external potential as a function of λ). The Kohn–Sham hxc is therefore a coupling strength average of hλxc as it evolves from λ = 0 through λ = 1 [372].

In an atom, the size of the hole is relatively insensitive to λ and remains of roughly atomic size. In a molecule, the changes in the character and size of the hole occur as λ

248 7 Kohn–Sham LCAO Method for Periodic Systems

varies from its noninteracting λ = 0 to its fully interacting λ = 1 limits. At λ = 0 (pure exchange with no correlation whatsoever) delocalization of the hole over two or more centers is characteristic. Accurate DFTs must recognize this λ = 0 nonlocality. As λ increases, the hole is localized by the long-range, nondynamical, left–right correlations that are absent in atoms but operative in molecules. At the fully interacting limit (λ = 1), the hole has been localized to roughly atomic size. To incorporate the necessary nonlocality Becke [373] proposed so-called hybrid functionals

Exchyb = ExcKS + amix ExHF − ExKS (7.41)

The di erence from the EXX functional is the following: the exact HF exchange ExHF is mixed with the DFT (LDA,GGA) exchange-correlation energy in hybrid functionals and with the DFT correlation energy in EXX functionals. The mixing parameter amix is used to patch an appropriate amount of HF exchange into the exchange-correlation energy. Even the simplest “half-half” hybrid functional (amix = 0.5) greatly improves the calculated properties, compared with the pure DFT results. Becke developed the 3-parameter functional expressions known as B3PW and B3LYP hybrid functionals:

ExcB3P W	= ExcLSDA	+ a ExHF		ExLSDA		+ b∆ExBecke	+ c∆ExP W	(7.42)
B3LY P	LSDA	HF	−	LSDA		Becke	LY P	(7.43)
Exc	= Exc	+ a Ex	− Ex			+ b∆Ex	+ c∆Ex
The B3PW functional (7.42) uses in Exc the					Becke [358] exchange and Perdew–Wang

exchange-correlation [355], while in the B3LYP functional (7.43) the correlation part is that suggested by Lee–Yang–Parr [359]. The a, b, c parameters were optimized to fit the experimental data and do not depend on the molecule under consideration. It is clear that the choice of hybrid functionals is motivated by reasonable physical

arguments. The term a ExHF − ExLSDA		replaces some electron-gas exchange with
the exact exchange to capture the	proper small-λ limit in (7.39). The coe cient a

reflects the rate of onset correlation as it increases from zero. The other terms allow optimum admixtures of exchange and correlation-type gradient corrections. These hybrid functionals are the simplest mixture of exact exchange, the LSDA for exchangecorrelation, and gradient corrections of exchange and correlation type, that exactly recovers the uniform electron-gas limit. The B3LYP functional is the most popular in molecular quantum chemistry. Nevertheless, the correct amount of HF exchange included in any hybrid functional cannot be a constant over all species or even all geometries of a single species [374]. The rationale for mixing exact exchange with DFT approximation is discussed in [375]. The authors write the hybrid functional in

the form
1		Ex − ExKS
Exc = ExcKS +			(7.44)
	n

where the optimum integer n can be found from the perturbation theory and n = 4 for atomization energies of typical molecules. Such a formally parameter-free hybrid functional with a PBE exchange-correlation part is known as the PBE0 hybrid functional. The simplification of hybrid exchange-correlation functionals was suggested by Becke [372] based on simulation of delocalized exact exchange by local densityfunctionals. A simple model was introduced that detects exchange hole delocalization in molecules through a local variable related to kinetic-energy density.

7.2 Density-functional LCAO Methods for Solids

249

Is it possible to use for extended systems (crystalline solids) the same hybrid functionals that are used in molecular DFT calculations? This problem is discussed in the next section.

7.2 Density-functional LCAO Methods for Solids

7.2.1 Implementation of Kohn–Sham LCAO Method

in Crystals Calculations

Although the same fundamental theory underlies the density-functional approach in both solid-state physics and molecular quantum chemistry, and the same exchangecorrelation functionals can be used in both disciplines (see Sect. 7.2.4), the detailed implementation is usually quite di erent [376]. First, this di erence concerns the basis-set choice. Since the electronic charge density of an isolated molecule is necessarily localized in a finite region of space, the traditional quantum-chemical method is LCGTO (MOs are expanded in a basis of localized Gaussian-type orbitals, centered on the atomic nuclei). In contrast, the tradition in solid-state physics is to use plane waves (PW) as a basis set for expanding the one-particle Kohn–Sham wavefunctions (crystalline orbitals). These can be used in their pure form in the first-principles pseudopotential method [377] or modified near the atomic cores (augmented plane waves [378]). The all-electron (termed also full potential – FP) PW calculations are practically difficult as the core-states description requires a huge number of plane waves. The plane waves are a reasonable first approximation to conduction-band eigenfunctions and permit many formal simplifications and computational economies. Nevertheless, PW methods have the disadvantage that very localized or inhomogeneous systems may require excessive numbers of waves for their representation. This problem can be especially serious in surface calculations, where the requirement of periodicity in three dimensions leads to a model with slabs repeating periodically along the surface normal. If the space between the slabs is wide enough to make interactions between them negligible, the “lattice constant” normal to the surface has to be quite large and the corresponding reciprocal lattice vector quite short, resulting in a large number of plane waves within a given kinetic-energy cuto and hence a relatively expensive calculation. This problem is discussed in Chap. 11 in more detail.

The distance between the slabs needs to be even larger if adsorbate molecules are to be added to the surface, since these molecules also should not interact across the space between the slabs. Localized-basis approaches to periodic systems do not su er from this disadvantage. Very localized states, including core states if desired, can be represented by a suitable choice of atom-centered basis functions. The electronic charge density goes naturally to zero in regions of space where there are no atoms, and the lack of basis functions in such regions is not a problem unless a detailed description of energetic excitations, scattering states, or tunneling is required [376]. For surface modeling there is no need for periodicity normal to the surface; a single slab of su cient finite thickness can be used, since these methods can be formulated equally easily for systems periodic in one, two or three dimensions. This facilitates the accurate treatment of molecule–surface interactions, especially when the molecules are relatively far from the surface. Also, if the basis functions are carefully constructed and optimized for the crystalline environment (we discuss this problem in Chap.

250 7 Kohn–Sham LCAO Method for Periodic Systems

8) it is possible to represent the valence-electron eigenstates with a relatively small number of basis functions: just enough to accommodate all the electrons, plus a few more functions to give the Kohn–Sham orbitals su cient variational freedom. Thus, only a few tens of basis functions per atom are typically needed, versus hundreds per atom in typical plane-wave approaches. The auxiliary basis sets needed to represent the charge density and, if desired, the exchange-correlation potential are also fairly modest in size. Once the basis set has been constructed and the required overlap and Hamiltonian matrices formed, the small (compared to plane wave) size of the basis can greatly reduce the computational cost of solving the Kohn–Sham equations.

Also, especially in systems with large unit cells, one can exploit sparcity in these matrices, since basis functions on distant centers will have negligible overlap and interaction (except for certain long-range Coulomb multipole interactions that can be summed according to Ewald’s convention or screened). Advanced techniques can be applied to calculate multicenter integrals in LCGTO methods for solids. For example, in the implementation of periodic boundary conditions in the MO LCAO program [379, 380] all terms contributing to the KS Hamiltonian are evaluated in real space, including the infinite Coulomb summations, which are calculated with the aid of the fast multipole method [381, 382]. In LCGTO methods with PBC the O(N ) linearscaling DFT calculations for large and complex systems are possible [379, 383]. For example, carbon nanostructures up to C540 were calculated [384] using DFT LCAO with the numerical atomic basis. Below, we present the basic formalism for DFT KS equations for periodic systems using a basis set of atomic orbitals. For more details of DFT LCAO methods the reader is referred to the original publications [376,385–389].

The HF LCAO method for periodic systems was considered in Sections 4.1.5 and 4.1.6. We discuss here the KS LCAO method for crystals in comparison with the HF LCAO approach. The electronic energy of the crystal (per primitive cell) as calculated within the HF approximation (EHF ) and DFT (EDF T ) can be expressed in terms of the one-electron density matrix (DM) of the crystal defined as P (k) in terms of Bloch sums of AOs, (4.125) or as ρ(R, R ) in coordinate space, (4.126). These expressions are:

EHF [ρ] = E0[ρ] + EH [ρ] + EX [ρ]

(7.45)

EDF T [ρ] = E0[ρ] + EH [ρ] + EXC [ρ]

(7.46)

ˆ

where E0[ρ] is defined as the expectation value of the one-electron operator h(R)

1

d3R[hˆ(R)ρ(R, R )]R =R

E0[ρ] =

(7.47)

N

VN

EH [ρ] is the Coulomb (Hartree) energy,

R

EH [ρ] = N d3R

d3R

R

(7.48)

1

ρ(R, R)ρ(R

, R

)

VN

VN

|

−

|

EX [ρ] is the HF exchange energy,

d3R

d3R |ρ|R − R ||

EX [ρ] = −2N

(7.49)

1

(R, R )

2

VN VN