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

7.2 Density-functional LCAO Methods for Solids

271

DFT functionals must fulfill the equations, which describe their values in the limit of

a single-electron system (hydrogen atom, for example) with density ρ (r), where σ

/ σ may be α, ρα(r) = 1, and ρβ (r) = 0, [440]. The first equation

EJ [ρ] + EX [ρα, 0] = 0

(7.66)

expresses that a single electron does not interact with itself, i.e. the self-repulsion energy (Hartree energy) EJ [ρ] is canceled by the self-exchange energy covered by the exchange functional EX [ρ]. The second equation

Ec[ρσ , 0] = 0

(7.67)

clarifies that a single electron does not possess any correlation energy (self-correlation is zero). The next two equations

vJ ([ρ]; r) + vxσ ([ρσ , 0]; r) = const			(7.68)
vσ([ρ	σ	, 0]; r) = 0	(7.69)
c

make sure that the single electron moves under the influence of the external potential v(r) rather than the Coulomb potential vJ , the exchange potential vx or the correlation potential vc. An approximate exchange-correlation functional may violate all or some of the equations (7.66)–(7.69) and, therefore, has to be corrected. The common LDAs, GGAs or Meta-GGAs have a remaining self-interaction error in their exchange part [436]. This self-interaction error is particularly critical for localized oneelectron states in molecules and crystals. The delocalized electrons move fast, thus experiencing mainly the LDA (LSDA) mean-field potential; this is why LDA methods are so fruitful for metals. The localized electrons reside on each atomic site for so long that the surroundings must respond to their presence [437]. Therefore, the LDA theory problems are the most severe in systems where the electrons tend to be localized and strongly interacting (strongly correlated systems). As was already noted, the well-known examples of such systems are the crystalline transition-metal oxides (with localized d-electron states) and rare earth elements and compounds (with localized f - electrons). These systems exhibit phenomena associated with strong correlation e ects (metal–insulator transitions, high temperature superconductivity, etc.). The already discussed hybrid DFT methods are most popular both in molecules and solids. The hybrid exchange-correlation functionals include HF exchange with some mixing coefficient and in this way partly compensate the self-interaction error. Other possibilities of DFT extension to the strongly correlated systems are known: the self-interaction- corrected (SIC) DFT method by Perdew–Zunger [352] (PZ-SIC approach) and the LDA+U method by Anisimov–Zaanen–Anderson [153]. The SIC-LSDA and LDA+U methods are described in many publications and applied to di erent molecules and solids. We refer the reader to reviews [438–441].

Let the many-electron system consist of N = Nα +Nβ electrons, where Nα and Nβ are the number of electrons with spin α and spin β, respectively. The total density of electrons with spin σ is a sum of orbital densities ρiσ (r) = |ψiσ(r)|2. The total density

Nσ

of the N -electron system is ρ(r) = ρiσ(r). In the KS method the energy of a

σ i=1

many-electron system is expressed as

272 7 Kohn–Sham LCAO Method for Periodic Systems

Nσ
	1
EKS [ρα, ρβ ] = σ=α,β i=1 ψiσ | −	2	∆|ψiσ +		ρ(r)v(r)dr + J[ρ] + Exc[ρα, ρβ ] (7.70)

The first two terms in (7.70) are the kinetic energy of a system of noninteracting electrons and the interaction between the electron density ρ(r) and external potential v(r). The Hartree energy J[ρ] = 1 / / ρ(r)ρ(r ) drdr is the Coulomb interaction of an

2|r−r |

electron density with itself, Exc is the exchange-correlation density-functional. The SIC corrected density-functional

	Nσ
ESIC KS = EKS −	σ
	(J[ρiσ ] + Exc[ρiσ, 0]) = EKS + ESIC	(7.71)

=α,β i=1

includes the self-interaction correction ESIC for all occupied spin-orbitals, which eliminates all spurious self-interaction terms orbital by orbital. In the SIC-DFT we have Ec[ρα, 0] = 0 and Ex[ρα, 0] = −EJ [ρ], i.e. an approximate exchange-correlation functional fulfills the exact conditions (7.67)–(7.68). The energy of equation (7.71) must be made stationary with regard to a mixing of occupied with occupied and occupied with virtual orbitals, which is accomplished by solving self-consistent KS equations for the SIC-DFT method. The SIC method gives significant improvements over the LSDA (SGGA) results [438], but there are also significant di culties in applying the method. The SIC exchange-correlation potential becomes orbital dependent, and since the KS orbitals thereby are not solutions to the exact same Hamiltonian, one can not in general be sure that the orbitals become orthogonal. Due to the mixing of occupied orbitals in the SCF-SIC-DFT procedure, the energy functional is no longer invariant with regard to unitary transformations of the orbitals. To circumvent these problems the Optimized Potential Method was suggested [151, 339] allowing the local singleparticle potential in the scheme that is both self-interaction free and orbital independent to be found. The simplest approximate implementation of the SIC-DFT method in molecular codes is the SIC correction applied after a self-consistent KS calculation is done [442]. The references to other SIC DFT method applications for molecules can be found in [443], where the SIC-DFT method was implemented self-consistently using a direct minimization approach. This approach was applied to calculate ionization potentials and electron a nities of atoms and small molecules [444]. It was concluded that this method works for molecules worse than for atoms, in particular it significantly overcorrects the one-electron energy for the highest-occupied MO.

In condensed matter studies the SIC DFT calculations have been extensively used. A review of solid-state SIC DFT techniques can be found in [438, 445, 446]. The wellknown most serious failure of LSDA is an underestimation of energy gaps for semiconductors and insulators and the suppression of magnetic ordering for antiferromagnetic insulators. For transition-metal monoxide insulators with the sodium chloride structure (MnO, CoO, FeO, and NiO) the relative positions of transition-metal d bands and oxygen p bands are important to clarify whether they are Mott–Hubbard insulators or charge-transfer insulators. In Mott–Hubbard insulators the occupied bands of transition-metal atoms appear at higher energy than the occupied p bands of oxygen atoms, while the relative positions are reversed in charge-transfer insulators. The lowest-unoccupied band is, in both cases, transition-metal d states. In the LSDA, the transition-metal d bands appear at higher energy than oxygen p bands and therefore

7.2 Density-functional LCAO Methods for Solids

273

they were considered as Mott–Hubbard-type materials. Such LSDA results do not agree with those by experiments and with the cluster calculations including correlation e ects by the configuration interaction (CI) method. It was shown that the SIC removes unphysical self-interaction for occupied orbitals and decreases their energies [438,447]. As in the case of molecules, the application of the LSDA SIC to solids has severe problems. Since the LSDA SIC energy functional is not invariant under the unitary transformation of the occupied orbitals one can construct many solutions in the LSDA SIC. If one chooses delocalized Bloch orbitals, the orbital charge densities vanish in the infinite-volume limit so that the SIC energy is exactly zero for such orbitals. But the localized Wannier orbitals constructed from Bloch orbitals have finite SIC energies. In many calculations for solids, the SIC was adapted to the localized orbitals, which are selected under some physical assumption. The criterion to choose these orbitals was carefully examined in [447], trying both solutions with localized and extended oxygen 2p orbitals. In solutions of types 1 (LSDA) and 2 all valencestate orbitals were extended or localized, respectively. In the solution of type 3 the transition-metal d-orbitals were localized and the oxygen p-orbitals were Bloch-type functions from the SIC potentials. In Table 7.6 are given bandgaps and magnetic moments of transition-metal monoxides for all the three types of LDA-SIC solution. The LDA+U results are discussed later. The type-1 solution is just the LSD result where all orbitals are extended, narrow d-bands appear at higher energy than the oxygen p-bands for all compounds. The bandgap of MnO is about 30% of the experimental value, there is no gap for other compounds. For the type 2 solution the transition-metal orbitals are more localized than those of oxygen p-orbitals so that the SIC potentials are larger in transition-metal d orbitals and the energies of d orbitals are pulled down significantly. The d- and p-bands become strongly hybridized, in good agreement with the observed results. As is seen from Table 7.6, the energy gaps are overestimated but the magnetic moments are well estimated. In type-3 solutions the transition-metal d-orbitals are localized and the oxygen p-orbitals are Bloch-type functions free from the SIC potential, so that the p-orbitals energies do not shift from the position of the LSDA results. On the contrary, the transition-metal d-orbitals energies are pulled down below those of the oxygen p-orbitals. Therefore, the highest-occupied band has the character of the oxygen p-band and the lowest-unoccupied band keeps mainly a character of the transition-metal d-orbitals. Then, the transition-metal monoxides are the charge-transfer type. As is seen from Table 7.6 the LDA-SIC values of the bandgaps and magnetic moments are in good agreement with the experiment. Unfortunately, the hybridization of p- and d-bands for occupied states, which appeared in the type-2 solution, is not reproduced in the type-3 solution. At the same time the type-2 solution overestimates the bandgaps. It was cleared up [438] that this overestimation is due to the overestimate in LSDSIC the Coulomb repulsion U of localized d-electron states in the solid (the SIC splits the occupied and unoccupied d-bands by the Coulomb repulsion U or slightly less than that). In condensed systems, the electron–electron Coulomb interaction is screened and the Coulomb repulsion should be much reduced. The “bare” (HF) value of the Coulomb interaction parameter U is 15–20 eV, while the screening in a solid to a much smaller value 8 eV or less [438]. This means that the LSDA-SIC functional form should be reformulated to include this screening. This is confirmed by investigations of f -electron bands in rare-earth element compounds [437].

274 7 Kohn–Sham LCAO Method for Periodic Systems

Table 7.6. Bandgaps and magnetic moments of transition-metal monoxides in SIC LSDA, [447] and LSDA+U, [439] calculations

)

B µ(

momentMagnetic

Bandgap (eV)

expt.		4.79,4.58	3.32	3.35,3.38	1.77,1.64,1.93
LSD	+U	1.67	3.62	2.63	1.59
Type	3	4.7	3.6	2.6	1.5
Type	2	4.7	3.7	2.7	1.7
Type1	(LSD)	4.5	3.2	2.0	0.7
Expt.		3.6–3.8	2.4	2.4	4.0,4.3
LSD	+U	3.5	3.2	3.2	3.1
Type	3	3.4	3.4	2.7	2.8
Type	2	6.5	6.1	5.3	5.6
Type1	(LSD)	1.1	0.0	0.0	0.0
Compound		MnO	Fe0	CoO	NiO

7.2 Density-functional LCAO Methods for Solids

275

The unscreened nature of Coulomb interaction is a serious problem of the Hartree– Fock approximation also – due to the neglect of screening, the HF energy gap values are a factor of 2–3 larger than the experimental values. A combination of the LDA and HF-type approximation called the LDA+U method was proposed [153]. The DFT+U solution can be obtained either at the LDA or GGA levels, giving rise to what have been called the LDA+U and the GGA+U formulations, respectively. In this approach the contribution of the interaction to U is added to the LDA energy. Instead of the energy-dependent interaction the averaged static Coulomb interaction between localized orbitals is used. This approach can be viewed as a way of connecting the theory of strongly correlated systems, which is based on such models, with the abinitio electronic-structure calculation methods like the HF method [151]. In contrast to the GW approximation, where the screened Coulomb interaction is used in the form of a nonlocal energy-dependent operator, in the DFT+U method a set of sitecentered atomic-type orbitals is introduced and the Coulomb interaction is present only between electrons on such orbitals. This can be called the “orbital-restricted” solution because the DFT-potential is the same for all orbitals (for example, for all five 3d-orbitals of a transition-metal ion). The DFT+U method can reproduce the splitting of d- (or f -) bands into occupied lower and unoccupied upper bands, which was the main problem of LDA. This leads to the possibility of investigating the influence of correlation e ects on structural properties, such as Jahn–Teller distortions. The orbital dependence of the DFT+U potential allows one to treat orbital and charge ordering in transition-metal compounds, which is very important for modern materials with anomalous magnetic and electronic properties.

In the LDA+U method electrons are separated into two subsystems: localized d- or

f -electrons for which Coulomb d–d (f –f ) interaction is taken into account by a term

12 U ninj (ni are d- or f -orbital occupancies) as in a mean-field (HF) approximation

i =j

and delocalized s, p-electrons that could be described by using an orbital-independent one-electron LDA potential.

In the LDA+U method the total Coulomb-interaction energy expression E =

U N (N − 1)/2(N = ni) is substracted from the LDA total energy and the orbital-

i

	dependent Coulomb-interaction energy is added:
	ELDA+U = ELDA − U N (N − 1)/2 +	1			(7.72)
			U	ninj
		2
				i =j

The orbital energies ε are derivatives of (7.72) with respect to the orbital occupations

εi =	∂E	= εLDA + U (	1		− ni)	(7.73)
	∂ni		2

This simple formula shifts the LDA orbital energy by −U/2 for occupied orbitals (ni = 1) and by +U/2 for unoccupied orbitals (ni = 0). A similar formula is found for

the orbital-dependent potential (Vi(r) = δE[ρ] ) where a variation is taken not on the

δρi(r)

total charge density ρ(r) but on the charge density of a particular ith orbital ρi(r)

Vi(r) = VLDA(r) + U (	1	− ni)	(7.74)
	2

276 7 Kohn–Sham LCAO Method for Periodic Systems

Such an approach introduces, in fact, the self-interaction correction for these orbitals. Therefore DFT-SIC and DFT+U methods give close results in many cases.

The LDA+U orbital-dependent potential (7.74) gives the energy separation between the upper valence and lower conduction bands equal to the Coulomb parameter U, thus reproducing qualitatively the correct physics for Mott–Hubbard insulators. To construct a calculation in the LDA+U scheme one needs to define an orbital basis set and to take into account properly the direct and exchange Coulomb interactions inside a partially filled d-(f -) electron subsystem [439]. To realize the LDA+U method one needs the identification of regions in a space where the atomic characteristics of the electronic states have largely survived (“atomic spheres”). The most straightforward would be to use an atomic-orbital-type basis set such as LMTO [448].

As is seen from Table 7.6 the SIC LSDA (type-3) and LSDA+U calculations give close results both for bandgaps and magnetic moments in transition-metal crystalline oxides. In both cases the agreement with the experiment appears to be essentially better than in pure DFT LSDA calculations.

The deficiency of the DFT+U method is the necessity to define explicitly the set of localized orbitals for interacting electrons. While for rare-earth and transition-metal ions the good approximation is to use atomic-type f - or d-orbitals and LMTO basis, for more extended systems like, for example, semiconductors, some more complicated Wannier-type orbitals are needed. Furthermore, the DFT+U method does not necessarily provide a reliable way of treating the U term, while, in a sense, all the important e ects of strong correlation are defined by this term. Another problem of this method is that the estimate of U may depend on the choice of the basis [449]. Di erent choices of local orbitals such as the LMTO and LCAO give di erent U because the definition of the orbital becomes ambiguous whenever there is strong hybridization. In this situation the U value fitted to experimental data can also be used. Recently, [450] the DFT+U method was implemented in the computer code VASP [451], which uses the pseudopotential description for core electrons (see Chap. 8) and plane waves as the basis in the Bloch-functions calculations. Use of this modified code allows the systematic DFT+U study of correlation e ects in di erent solids. In particular, application of DFT+U method for transition-metal sulfides yielded improved predictions for volume, magnetic moment, exchange splitting and bandgap [450]. In these calculations the dependence of results on the U parameter is studied. It was shown that the observed width of the semiconducting gap is reached at U=5 eV in the LDA +U method and U=3 eV in the GGA+U method. The comparative DFT study with B3LYP (LCAO basis) and GGA+U (PW basis) of electronic-structure and magnetic coupling in FeSbO4 crystal [452] allows the di erences due to the di erent choice of exchange-correlation functional and basis to be estimated. In Table 7.7 we give the electronic bandgap values obtained in [452].

These results correspond to the experimentally found antiferromagnetic ordering of Fe ions in two dimensions. Such an ordering is correctly reproduced giving for the B3LYP functional the magnetic coupling constant value –21 K (experimental value is –25 K). From experimental studies, FeSbO4 is known to be a charge-transfer semiconductor (with an O2p-Fe3d bandgap), with an activation energy of 0.75 eV for n-type conduction. This value represents a lower limit for the bandgap, as is the energy required to excite electrons from the donor levels to the conduction band. In DFT+U calculations the parameter Uef f is defined as the di erence U–J (J is

7.2 Density-functional LCAO Methods for Solids

277

Table 7.7. Electronic bandgap values (in eV) in FeSbO4 crystal, obtained for the antiferromagnetic solutions using di erent methods , [452]

LCAO	Bandgap	Plane Waves	Bandgap
GGA(BLYP)	0.2	GGA(PW91)	0.2
B3LYP	3.0	GGA+U (Uef f =4eV)	1.6
HF	12.8	GGA+U (Uef f =9eV)	1.7

Experimental value is estimated as being greater than 0.75 eV (see explanations to this table in the text)

a parameter representing the screened exchange energy, being almost constant at 1 eV). It was concluded [452] that the bandgap and the relative positions of the bands are a ected in similar ways by the inclusion of 20% Hartree–Fock exchange and by the GGA+U correction with Uef f = 4 eV. In particular, the occupied Fe-3d levels shift down in the valence-band region for both methods, in comparison with the GGA solutions. However, increasing the HF exchange up to 100% and using higher values of Uef f do not have equivalent e ects on the band structure, due to the local character of the GGA+U correction.

As a main conclusion of these calculations we note the consistency between B3LYP and GGA+U results. This means that LCAO calculations with hybrid functionals can play the same role in quantum chemistry of solids as DFT+U and SIC DFT calculations in solid-state physics.

Recent LDA+U calculations were made on di erent solids: ZnX( X=O, S, Se, Te), [453], NiO [454] and other transition-metal oxides (Cu2O [455], ReO2 [456], T iO2 [457]). The results of these calculations confirm the e ciency of the DFT+U approach for strongly correlated systems. We hope that future LCAO calculations with the hybrid exchange-correlation functionals will discover the same e ciency.

This chapter concludes Part I (Theory) of our book. In Part II (Applications) we consider the applications in the quantum chemistry of solids of the theoretical methods described in Part I.