PHYSICAL REVIEW B VOLUME 54, NUMBER 23 15 DECEMBER 1996-I

Generalized gradient approximation for the exchange-correlation hole of a many-electron system

John P. Perdew
Department of Physics and Quantum Theory Group, Tulane University, New Orleans, Louisiana 70118

Kieron Burke
Department of Chemistry, Rutgers University, Camden, New Jersey 08102

Yue Wang
Department of Biostatistics, School of Public Health, University of North Carolina, Chapel Hill, North Carolina 27599-7400
~Received 14 February 1996; revised manuscript received 29 July 1996!
We construct a generalized gradient approximation ~GGA! for the density n xc(r,r1u) at position r1u of
the exchange-correlation hole surrounding an electron at r, or more precisely for its system and spherical
average ^ n xc(u) & 5(4 p ) 21 * dV u N 21 * d 3 r n(r)n xc(r,r1u). Starting from the second-order density gradient
expansion, which involves the local spin densities n (r),n (r) and their gradients n (r),n (r), we cut off
the spurious large-u contributions to restore those exact conditions on the hole that the local spin density ~LSD!
approximation respects. Our GGA hole recovers the Perdew-Wang 1991 and Perdew-Burke-Ernzerhof GGAs
for the exchange-correlation energy, which therefore respect the same powerful hole constraints as LSD. When
applied to real systems, our hole model provides a more detailed test of these energy functionals, and also
predicts the observable electron-electron structure factor. @S0163-1829~96!04847-3#

I. INTRODUCTION AND SUMMARY See Refs. 14 and 15 for reviews. LSD models the hole den-
OF CONCLUSIONS sity as

xc ~ r,r1u ! 5n xc n ~ r ! ,n ~ r ! ;u,
n LSD ~4!
unif
Kohn-Sham density functional theory14 would yield the
exact ground-state energy E and spin densities n (r),n (r) where n unif
xc (n ,n ;u) is the hole density in an electron gas
of a many-electron system, if the exact exchange-correlation
of uniform spin densities n ,n , for which an accurate ana-
energy functional E xc@ n ,n # were known. Good results for
lytic model now exists.15 Thus LSD respects the sum rules
solids are often found from the local spin density ~LSD!
approximation1
E d 3 u n x~ r,r1u! 521, ~5!

xc @ n ,n # 5
E LSD E d 3 r n ~ r! e xcn ~ r! ,n ~ r! , ~1!

where e xc(n ,n ) is the known5 exchange-correlation energy

E d 3 u n c~ r,r1u! 50, ~6!

per particle of an electron gas with uniform spin densities which constrain the energy integral of Eq. ~3!. Moreover, Eq.
n , n , and n5n 1n . Equation ~1! is clearly valid when ~3! involves only the system and spherically averaged hole
the spin densities vary slowly over space, but this condition density
is violated by real atoms, molecules, and solids. Indeed, the
next systematic correction in the slowly varying limit, the
second-order gradient expansion approximation611 ~GEA!, ^ n xc~ u ! & 5 E dV u 1
4p N
E d 3 r n ~ r! n xc~ r,r1u! , ~7!

xc @ n ,n # 5E xc @ n ,n # 1
E GEA LSD
(
s,s8
E 3
d r C ss 8
xc ~ n ,n !
for which LSD provides a fairly reasonable description.16
The LSD hole has other desirable features, such as the
negativity condition on the exchange hole17
2/3
3n s n s 8 /n s2/3n s 8 , ~2!
n x~ r,r1u! <0, ~8!
is less accurate than LSD.6,12,13 an accurate ~but inexact ! on-top value n xc(r,r), and the
18 19
Gunnarsson and Lundqvist14 explained the success of correct electron-electron cusp condition at u50.20 GEA im-
LSD as follows: The exchange-correlation energy is the elec- proves upon LSD at small u, but displays spurious large-u
trostatic interaction of each electron at r with the density behavior,21,17,12,22 which is sampled by the long-range Cou-
n xc(r,r1u)5n x1n c at r1u of the exchange-correlation lomb interaction 1/u. Langreth and Perdew21 showed that the
hole which surrounds it. In atomic units (\5e 2 5m51), GEA hole, which is a truncated expansion and not the exact
hole of any system, violates Eq. ~6!.
E xc@ n ,n # 5
1
2
E d 3 r n ~ r! E d 3 u n xc~ r,r1u! /u. ~3!
Nonempirical generalized gradient approximations
~GGAs!,12,23

0163-1829/96/54~23!/16533~7!/$10.00 54 16 533 1996 The American Physical Society

16 534 JOHN P. PERDEW, KIERON BURKE, AND YUE WANG 54

where m 50.219 51 ~to preserve the good LSD description of

the exchange-correlation energy in the linear response of the
uniform gas! and k 50.804. In Sec. II, we construct a GGA
hole that accurately reproduces Eq. ~15!, by applying a
damping factor to the PW86 exchange hole. The damping
factor, used only for exchange, reflects the more pathological
large-u behavior of the n GEA
x , and the double nature of its
GGA cutoffs, which enforce both Eqs. ~5! and ~8!.
In Sec. III, we present our model for the GGA correlation
hole. The first such model,23 which led to the PW91 corre-
lation energy functional, was based upon sharp cutoffs of
crude approximations for both the LSD and gradient contri-
butions to the hole. We refine these crude approximations,
FIG. 1. Spherically averaged exchange hole density n X for but find essentially the same correlation energy, which can
s51 in LSD ~circles!, GEA ~crosses!, and GGA ~solid line!. be accurately represented by the PBE functional

xc @ n ,n # 5
E GGA E d 3 r f ~ n ,n ,n ,n ! , ~9! c @ n ,n # 5
E PBE E d 3 r n $ e c~ r s , z ! 1H PBE~ r s , z ,t ! % , ~16!

where
often start from the GEA for the hole n GEA
xc and cut off its
large-u contributions to restore exact conditions such as Eqs. r s 5 ~ 3/4p n ! 1/3, ~17!
~5!, ~6!, and ~8!. Since only the system average of Eq. ~7! is
needed, 2 n contributions to the GEA are first transformed z 5 ~ n 2n ! /n, ~18!
via integration by parts on r. GGAs may be applied directly
or hybridized with exact exchange.24,25 t5 u n u /2k s f n, ~19!
In Sec. II, we present our GGA model for the exchange
hole. The first such model was that of Perdew and Wang12 in f 5 21 @~ 11 z ! 2/31 ~ 12 z ! 2/3# , ~20!
1986 ~PW86!, who used sharp cutoffs on n GEA x to enforce
Eqs. ~5! and ~8!, yielding k s 5 ~ 4k F / p ! 1/2, ~21!

E GGA
1 GGA
x @ n ,n # 5 E x
2
1
@ 2n # 1 E GGA
2 x
@ 2n # , ~10! H PBE5 gf 3 ln 11 H b 2
t F
11At 2
g 11At 2 1A 2 t 4 GJ , ~22!

x @n#5
E GGA E d 3 r n e unif
x ~ n ! F x~ s ! , ~11! A5
b
g
@ exp$ 2 e unif
c / gf % 21 #
3 21
~23!

where and g 50.031 091, b 50.066 725. The reduced gradients s

and t measure how fast n(r) is varying on the scales of the
e unif
x ~ n ! 523k F /4p , ~12! local Fermi wavelength 2 p /k F and the local Thomas-Fermi
screening length 1/k s , respectively.
k F 5 ~ 3 p 2 n ! 1/3, ~13! In Ref. 30, Eq. ~22! was derived from three limits:

s5 u n u /2k F n. ~14! H PBE bf 3 t 2 ~ t0 ! , ~24!

The real-space cutoff gave a numerical function F x(s) ~see H PBE2 e unif ~ t` ! , ~25!
C
Fig. 1 of Ref. 12!, which was fitted to an analytic form,
F PW86(s). In the later work of Perdew and Wang23,26 in 1991 and
x
~PW91!, Beckes27,28 semiempirical refinements plus addi-
c @ n g ,n g # const
E PBE ~ g ` ! , ~26!
tional theoretical constraints were included in F PW91
x (s), al-
though F PW91
x (s) was a worse fit to the numerical function where n s g (r)5 g 3 n s ( g r) is a uniformly scaled density.31
than was F PW86
x (s). Both the PW86 and PW91 parametriza- These limits also emerge naturally from the real-space cutoff
tions were contorted29 at small s to recover the expected construction of Sec. III, as shown in Ref. 32. The high-
GEA of Eq. ~2!. density limit of Eq. ~26! is violated by both LSD and
Recently, Perdew, Burke, and Ernzerhof ~PBE!30 pre- PW91.33
sented a simplified construction of a simplified GGA for ex- Thus the PBE correlation energy functional of Eq. ~16!
change and correlation, in which all parameters ~other than can be derived either from various limits, as in Ref. 30, or
those in LSD! are fundamental constants. Although indepen- from a real-space construction of the GEA correlation hole,
dent of PW91 or any model for the hole, the PBE functional as in Sec. III. The PBE exchange energy functional of Eq.
is numerically equivalent to PW91 for most purposes, and ~15! is derived from its limits in Ref. 30, and is then used to
improve the real-space cutoff of the GEA exchange hole in
x ~ s ! 511 k 2 k / ~ 11 m s / k ! ,
F PBE ~15!
2
Sec. II.
54 GENERALIZED GRADIENT APPROXIMATION FOR THE . . . 16 535

Sharp cutoffs produce a choppy n GGAxc (r,r1u) ~e.g., is the reduced electron-electron separation on the scale of the
see the GGA exchange hole in Ref. 34!, but are smoothed in Fermi wavelength. The functions J(z), L(z), M (z), and
the system and spherical average of Eq. ~7!. In Sec. IV, we N(z) are known, oscillating functions of z, given in Eqs.
describe the general features of the GGA holes. ~8!~11! of Ref. 12, respectively. D(z) is a damping factor,
We conclude the following: ~1! The PW91 and PBE which equals 1 at the pure GEA level. To first order in u,
exchange-correlation energy functionals correspond to a n GEA
x 52n(r1u)/2, so that the GEA hole is deeper on the
GGA exchange-correlation hole, with the known correct fea- high-density side of the electron.
tures of the LSD hole, plus additional correct inhomogeneity In our generalized gradient approximation ~GGA!, the
effects. These functionals should therefore perform as reli- hole is represented as
ably as LSD, except under special conditions ~see below!. ~2!
x ~ r,r1u ! 52 2 n ~ r ! y u ~ y ! u u x~ r ! 2u,
n GGA ~32!
1
The GGA system-averaged exchange-correlation hole of Eq.
~7! can be constructed for any system via the formulas given
where u (x)51 for x.0 and 0 for x,0. The first step func-
in this work. The coupling-constant integration
tion on the right enforces the negativity condition of Eq. ~8!,

^ n xc~ u ! & 5 E0
1
dl ^ n xc,l~ u ! & ~27!
while the second involves a cutoff separation u x chosen to
enforce the normalization condition of Eq. ~5!, which be-
comes

E
can be undone, as in Refs. 35 and 15, to extract the physical
1 zx
exchange-correlation hole at full coupling strength (l51), 2 dz z 2 y sph av~ z,s ! 521, ~33!
for comparison with the results of accurate wave-function 12p 0
calculations. ~3! The system-averaged hole at full coupling
where z x52k F u x and
strength is an observable, since its Fourier transform is es-
sentially the electron-electron structure factor S(k) measured
in quasielastic scattering processes. We expect the GGA to
make a useful prediction for this observable.
y sph av~ z,s ! 5
1
4p
E dV uyz,s,uu y ~ z,s,u! . ~34!

Finally, the PW91 and PBE functionals yield great im- Equation ~33! determines z X as a function of s, and Eqs. ~3!
provements over LSD for the total energies of atoms ~and and ~11! yield
their separate exchange and correlation contributions! or at-
omization energies of molecules,26,3640 but have a mixed
history of successes and failures for solids.4151 This may be F x~ s ! 5
1
9
E
0
zX
dz zy sph av~ z,s ! . ~35!
because the exchange-correlation hole can have a diffuse
large-u tail in a solid, but not in an atom or small molecule, The angular integration over V u in Eq. ~34! is performed
where the density n(r) itself is well localized. As we shall analytically ~as explained at the end of this section! and the
see in Sec. II, a sharp radial cutoff corresponds to z integrations of Eqs. ~33! and ~35! are performed numeri-
k 50.804 in Eq. ~15!, while a more diffuse cutoff leads to a cally.
smaller value of k . This uncertainty is also reflected in the As discussed in Refs. 17 and 12, the GEA exchange hole
PBE derivation of Eq. ~15!, in which k is set to the maxi- displays an undamped cos(2k F u) oscillation as u`. We
mum value allowed by the Lieb-Oxford bound23,52 on E xc . damp this oscillation by taking

D ~ z ! 51/@ 11 ~ bz ! q # ~36!
II. EXCHANGE HOLE
in Eq. ~30!. To preserve the GEA hole at small u, q>2.
The exact exchange hole arises from the Kohn-Sham non- Since the LSD hole density is essentially confined to the
interacting wave function ~a Slater determinant! and satisfies region 0&z&2 p , and the GEA breaks down when the GEA
the spin-scaling relation hole density is much greater than the LSD density, we expect
b;1/2p . The choices q52.5 and b51/2p provide a good
n s ~ r!
n x@ n ,n #~ r,r1u! 5 (s n ~ r! x
n @ 2n s #~ r,r1u! . ~28! x (s) of Eq. ~15!.
fit to F PBE
Figure 1 is a plot of 2y sph av(z,s)z 2 /12p , appearing in
Thus we need only model n x@ n # (r,r1u), the exchange hole Eq. ~33!, as a function of z for s51.0, for the LSD hole and
as a functional of the density for a spin-unpolarized system. the damped GEA and GGA holes. For large z, the damped
We write GEA hole displays unphysical oscillations about n x50. The
GGA hole is sharply cut off at z x;10.5, to satisfy Eq. ~33!,
and is identically zero for 8&z x&9, due to the negativity
x ~ r,r1u ! 52 2 n ~ r ! y,
n GEA ~29!
1

cutoff in Eq. ~34!. Elsewhere, the GGA hole is not equal to

where the damped GEA hole because of the step function inside the
spherical average in Eq. ~34!, which produces derivative dis-
y ~ z,s,u! 5J ~ z ! 1D ~ z ! $ 4L ~ z ! us/3216M ~ z !~ us! 2 /27 continuities as a function of z, e.g., at z;7.
Figure 2 shows numerical results for the reduced cutoff
216N ~ z ! s 2 /3% . ~30!
radius z x as a function of s. As s0, the cutoff radius moves
Here u5u/u, s5n/2k F n, and out to `, and the GGA hole reduces to the damped GEA
hole. As s`, the cutoff radius slowly approaches zero, and
z52k F u ~31! the hole becomes highly localized around the electron. For
16 536 JOHN P. PERDEW, KIERON BURKE, AND YUE WANG 54

where A, B, and C are independent of n , and given by Eq.

~30!. The value of y sph av can be stated for several different
cases: If T,0, where T5B 2 24AC, or u n 6 u .1, where
n 6 52B/2A6 AT/2u A u , then
y sph av5 @ F ~ 1 ! 2F ~ 21 !# u ~ V ! , ~38!
with V5A2B1C and F( n )5A n /61B n /41C n /2. If 3 2

u n 1 u .1 but u n 2 u ,1, then

y sph av5 @ F ~ n 2 ! 2F ~ 21 !# u ~ V ! 1 @ F ~ 1 ! 2F ~ n 2 !# u ~ 2V ! .
~39!
If u n 2 u .1 but u n 1 u ,1, then

FIG. 2. Reduced cutoff separation for the GGA exchange hole y sph av5 @ F ~ n 1 ! 2F ~ 21 !# u ~ V ! 1 @ F ~ 1 ! 2F ~ n 1 !# u ~ 2V ! .
as a function of reduced density gradient. For s*3, we find ~40!
u x&r s . Finally, if u n 6 u ,1, then

intermediate s, the steps in z x(s) occur where the normaliza- y sph av5 @ F ~ n 2 ! 2F ~ 21 ! 1F ~ 1 ! 2F ~ n 1 !# u ~ V !
tion cutoff passes through a negativity cutoff. 1 @ F ~ n 1 ! 2F ~ n 2 !# u ~ 2V ! . ~41!
Figure 3 shows various numerical results for the enhance-
ment factor F x(s) over local exchange as a function of s in
the physical range53 0&s&3, as defined by Eq. ~35!, as well III. CORRELATION HOLE
as the PBE enhancement factor of Eq. ~15!, with which it The local density n(r) sets only one length scale for ex-
agrees closely. We also present the numerical GGA enhance- change ~the Fermi wavelength! but it sets a second for cor-
ment factor resulting from ~1! a diffuse radial cutoff factor relation, the Thomas-Fermi screening length. Since the cor-
@ 11(u/u x) 2 # exp@2(u/u x) 2 # , with u x fixed by Eq. ~33!, relation hole is not required to satisfy a negativity constraint
which leads to a smaller F x(s) for s*1, and ~2! with sharp like Eq. ~8!, we do not need the non-spherical component of
radial cutoffs but with no damping @i.e., D51 in Eq. ~36!#, its GEA density. We write the spherically averaged GEA
which leads to a slightly larger F x(s), and a linear depen- correlation hole as
dence on s as s0. Thus the real-space cutoff procedure
c ~ r,u ! 5n c ~ r s , z , v ! 1t d n c~ r s , z , v ! ,
n GEA ~42!
LSD 2
determines the general features of F x(s), but not its exact
behavior, which must then be fixed by other constraints @as
where
has been done in the construction of the PW86,12 PW91,23
and PBE ~Ref. 30! exchange functionals#. v 5 f k su ~43!
We close this section with a technical point: the analytic
result for the angular integral in Eq. ~34!. Let n be the cosine is the reduced electron-electron separation on the scale of the
of the angle between u and s. Then screening length.
c (r s , z , v ) is accu-
The LSD correlation hole function n LSD

E
1 1 rately known,15 and has been confirmed by recent quantum
y sph av5 d n ~ A n 2 1B n 1C ! u ~ A n 2 1B n 1C ! , ~37! Monte Carlo calculations.54 We write
2 21

c ~ r s , z , v ! 5 f k s A c~ r s , z , v ! ,
n LSD ~44!
5 2

where ~for r s <10)

4 p v 2 A c~ r s , z , v ! 5 f 1 ~ v ! 1 f 2 ~ r s , z , v ! . ~45!
Here f 1 ( v ) is a nonoscillatory long-range contribution,
known from the random phase approximation ~RPA!. @Cor-
rection to Eq. ~22! of Ref. 15: a 3 50.002 431 7#. The short-
range contribution f 2 ( v ) vanishes rapidly for v @1/Ap,
where

p ~ r s , z ! 5 p k F d ~ z ! /4 f 4 ~46!
and d( z )50.30520.136z . 2

Similarly, we write the GEA correlation hole as

FIG. 3. Enhancement factor over local exchange as a function of d n c~ r s , z , v ! 5 f 5 k 2s B c~ r s , z , v ! , ~47!

reduced density gradient. The solid line is the damped numerical
GGA of Eq. ~35!, the open circles are the PBE of Eq. ~15!, the B c~ r s , z , v ! 5B LM
c ~ v !@ 12exp~ 2 p v !#
2
pluses are the numerical results with no damping, and the crosses
are the numerical results with a diffuse radial cutoff. 1 b ~ r s , z !v 2 exp~ 2 p v 2 ! , ~48!
54 GENERALIZED GRADIENT APPROXIMATION FOR THE . . . 16 537

FIG. 4. The function B c(r s , z , v ), which defines the shape of the FIG. 5. Spherically averaged correlation hole density n c for
gradient correction to the correlation hole via Eq. ~48!, for r s 52 r s 52 and z 50. GEA holes are shown for four values of the re-
with z 50 and z 51. Also shown is the entire Langreth-Mehl duced density gradient, t5 u n u /(2k z n). The vertical lines indicate
(r s 5 z 50) curve. where the numerical GGA cuts off the GEA hole to make
v
* 0 Cd v 4 p v 2 n c( v )50.
where
c ~ r s , z ,t, v ! 5 f k s @ A c~ r s , z , v ! 1t B c~ r s , z , v !#
n GGA 5 2 2

2 21
c ~ v ! 5 @ 18p ~ 11 v /12 ! #
B LM ~49!
3 2
3 u ~ v c2 v ! , ~52!

is the nonoscillating long-range contribution, the Fourier where v c(r s , z ,t) is the largest root satisfying the normaliza-
transform of the Langreth-Mehl55 ~LM! exponential model tion condition
for the wave-vector decomposition of the gradient contribu-

E
tion to the correlation hole in the RPA.
vc
Equation ~48! contains no constant or linear term in v , so d v 4 p v 2 @ A c~ r s , z , v ! 1t 2 B c~ r s , z , v !# 50. ~53!
it does not alter the highly accurate LSD on-top hole or its 0
cusp. The coefficient b (r s , z ) is determined from the known
GEA correction to the energy, for which we use the Figure 5 is a plot of the spherically averaged GGA correla-
Langreth-Mehl approximation to the high-density (r s 0) tion hole for r s 52 and z 50 for several different values of
limit, with the spin dependence of Ref. 56: t. We see that, for a small value of t (t50.5), the gradient
correction to LSD is small and so v c is large ('10.0), tend-

DE GEA
c 5C LM
c E d 3 r f ~ z ! u n u 2 /n 4/3
ing to its GEA value (`) as t0. On the other hand, for
t51.5, the gradient correction is 9 times larger, causing the
cutoff to occur at a much smaller value of v c ( v c52.3). In
516~ 3/p ! 1/3C LM
c E d 3r n f 3t 2, ~50!
the limit t`, v c0 ~Fig. 6!, turning off the correlation
contribution altogether @as in Eq. ~25!#. In all cases, the GGA
correlation hole is more localized than either the LSD or
GEA holes.
where C LM c 5( p /3) /(24p ) differs only slightly from the
1/3 2

exact high-density limit of Ma and Brueckner,6 In Fig. 6, we follow v c as a function of t for r s 52, for
b ( p /3) 1/3/16. Requiring our GEA correlation hole to recover both the spin-unpolarized ( z 50) and the fully spin-polarized
this energy fixes

2p2
b~ rs ,z !5 @ 12E 1 ~ 12p !# , ~51!
3p2

where E 1 (x)5x exp(x)*`x dt exp(2t)/t. Figure 4 shows both

the Langreth-Mehl B c and how the short-range contribution
changes it for r s 52 and z 50,1.
A key fact is that Eq. ~47! and its normalization integral
are well defined but positive. While the LSD correlation hole
properly integrates to zero, the GEA correlation hole does
not.21 A simplification is that A c of Eq. ~45! and B c of Eq.
~48!, which are functions of three variables, correctly reduce
to functions of v alone in either the high-density (r s 0,
p`) or the long-range ( v `) limits, where RPA is valid. FIG. 6. Reduced cutoff separation for the GGA correlation hole,
With the GEA correlation hole fully defined, we construct for r s 52 and z 50 or z 51. For z 51, v c50 beyond t'2.5. For
the spherically averaged GGA hole the case r s 52 and z 50, u c&r s means t*1.4 and s*1.6.
16 538 JOHN P. PERDEW, KIERON BURKE, AND YUE WANG 54

FIG. 7. The function H5 e GGAc 2 e LSD

c for several values of r s for FIG. 8. Same as Fig. 7, but for the fully spin-polarized case
the spin-unpolarized case ( z 50). The solid lines are the numerical ( z 51).
result of the real-space cutoff procedure, while the open circles are
from the PBE parametrization of Eq. ~22!.
is diffused over a much greater length scale } Ar s , and each
( z 51) cases. Note that, for small t, v c is large, and the displays a simple scaling behavior. Thus exchange domi-
21
curves ( z 50 or 1) merge because the long-range contribu- nates correlation: e unif
x }r s and e unif
c }ln(rs). As r s increases,
tion to the hole is independent of z . The same qualitative the exchange hole expands more rapidly than the correlation
behavior occurs for all densities, although for higher densi- hole, and ultimately engulfs it. In the low-density (r s `)
ties (r s 0) the short-range contribution becomes negligible limit, the exchange and correlation holes each scale with r s
for all t and the two curves become everywhere identical. and are of comparable size. The extreme long-range part of
The GGA correlation energy per particle at position r is the correlation hole is always given by RPA, as is the
now r s 0 limit.
The gradient corrections to the hole are known less reli-

E
ably than the LSD terms. Second-order gradient corrections
vc 4pv2
e GGA
c ~ r s , z ;t ! 5 f
3
dv @ A c~ r s , z , v ! 1t 2 B c~ r s , z , v !# to the LSD on-top hole and cusp are small18 and neglected
0 2v here. In fact, self-consistent LSD or GGA calculations pro-
5 e unif vide accurate predictions57,58 for n(r) and n xc(r,r), even
c ~ r s , z ! 1H ~ r s , z ,t ! . ~54!
when symmetry breaking leads to serious errors in z (r) of
Eq. ~18!. The gradient corrections display the same small-
In Fig. 7, we plot the difference between the GGA and LSD
and large-r s scalings as do the LSD holes.
correlation energies for different values of r s , as a function
of t for the spin-unpolarized case ( z 50). For t 2 0, this A nonzero reduced gradient s at the position of the elec-
figure recovers the t 2 behavior of GEA. Since the GGA cor- tron affects the holes at small u, deepening the exchange part
relation energy vanishes at large gradients, the limit as and raising the correlation part. The large-u contributions to
t` in this figure is precisely 2 e unif the GGA holes are chopped off. This leads to the scaling
c . Figure 8 shows the
same for z 51. Figures 7 and 8 also show the fair agreement limit of Eq. ~26!, since r s (r) g 21 r s ( g r), s(r)s( g r), and
between the analytic PBE expression of Eq. ~22! and the t(r) g 1/2t( g r) under uniform scaling.
numerical GGA. The results for PW91 are very similar, ex- As the reduced gradient s increases, the exchange hole
cept that as r s 0, H PW91(t) becomes a simple parabola for @constrained by Eqs. ~5! and ~8!# becomes deeper and more
large t, due to the H 1 term in that functional.23 short ranged in u, so the negative exchange energy turns on
Finally, compare Figs. 7 and 8 to deduce an approximate more strongly. But the correlation hole @constrained by Eq.
spin-scaling relationship for r s &6: ~6!# is gradually cut down to zero, so the correlation energy
turns off. All these effects may be seen in the GGA exchange
and correlation energies, as depicted in Figs. 3, 7, and 8, and
e GGA
c ~ r s , z ;t ! ' f ~ z ! e c r s ,0;t/ f ~ z ! .
3 GGA
~55! in Fig. 1 of Ref. 30; see also Refs. 52 and 59.

IV. QUALITATIVE FEATURES OF THE GGA HOLES

ACKNOWLEDGMENTS
At the LSD level, our models are essentially exact. The
only missing ingredient is the large-u oscillation of the cor- This work has been supported by NSF Grants No.
relation hole, which is energetically unimportant15 and is DMR92-13755 and No. DMR95-21353. We thank Matthias
washed away in the average of Eq. ~7!. When the reduced Ernzerhof, Jingsong He, Mel Levy, and Cyrus Umrigar for
gradient s of Eq. ~14! vanishes at the position of the electron, useful discussions. K.B. also acknowledges the generous
our GGA holes reduce to LSD, and so are fixed by the local hospitality of both the Physics and Mathematics Departments
value of r s . In the high-density (r s 0) limit, the exchange of Trinity College, Dublin, where part of this work was per-
hole exists on a length scale }r s , while the correlation hole formed.
54 GENERALIZED GRADIENT APPROXIMATION FOR THE . . . 16 539

1
W. Kohn and L.J. Sham, Phys. Rev. 140, A1133 ~1965!. 31
M. Levy, Int. J. Quantum Chem. S23, 617 ~1989!.
2
R.O. Jones and O. Gunnarsson, Rev. Mod. Phys. 61, 689 ~1989!. 32
K. Burke, J. P. Perdew, and Y. Wang ~unpublished!.
3 33
R.G. Parr and W. Yang, Density Functional Theory of Atoms and C.J. Umrigar and X. Gonze, in High Performance Computing and
Molecules ~Oxford, New York, 1989!. its Application to the Physical Sciences, Proceedings of the
4
R.M. Dreizler and E.K.U. Gross, Density Functional Theory Mardi Gras 1993 Conference, edited by D.A. Browne et al.
~Springer-Verlag, Berlin, 1990!. ~World Scientific, Singapore, 1993!.
5
J.P. Perdew and Y. Wang, Phys. Rev. B 45, 13 244 ~1992!. 34
M. Slamet and V. Sahni, Phys. Rev. B 44, 10 921 ~1991!.
6
S.-K. Ma and K.A. Brueckner, Phys. Rev. 165, 18 ~1968!. 35
K. Burke, J.P. Perdew, and M. Levy, in Modern Density Func-
7
L.J. Sham, in Computational Methods in Band Theory, edited by tional Theory: A Tool for Chemistry, edited by J.M. Seminario
P.M. Marcus, J.F. Janak, and A.R. Williams ~Plenum, New and P. Politzer ~Elsevier, Amsterdam, 1995!.
York, 1971!, p. 458. 36
A.D. Becke, J. Chem. Phys. 97, 9173 ~1992!.
8
M. Rasolt and H.L. Davis, Phys. Lett. 86A, 45 ~1981!. 37
J.C. Grossman, L. Mitas, and K. Raghavachari, Phys. Rev. Lett.
9
M. Rasolt and D.J.W. Geldart, Phys. Rev. B 34, 1325 ~1986!.
75, 3870 ~1995!; 76, 1006~E! ~1996!.
10
L. Kleinman and S. Lee, Phys. Rev. B 37, 4634 ~1988!.
11
38
D. Porezag and M.R. Pederson, J. Chem. Phys. 102, 9345 ~1995!.
D.C. Langreth and S.H. Vosko, Adv. Quantum Chem. 21, 175 39
E.I. Proynov, E. Ruiz, A. Vela, and D.R. Salahub, Int. J. Quantum
~1990!.
12
J.P. Perdew and Y. Wang, Phys. Rev. B 33, 8800 ~1986!; 40, Chem. S29, 61 ~1995!.
3399~E! ~1989!.
40
A.C. Scheiner, J. Baker, and J.W. Andzelm, J. Comput. Chem. ~to
13
J.P. Perdew, Phys. Rev. B 33, 8822 ~1986!; 34, 7406~E! ~1986!. be published!.
41
14
O. Gunnarsson and B.I. Lundqvist, Phys. Rev. B 13, 4274 ~1976!. A. Garcia, C. Elsasser, J. Zhu, S.G. Louie, and M.L. Cohen, Phys.
15
J.P. Perdew and Y. Wang, Phys. Rev. B 46, 12 947 ~1992!. Rev. B 46, 9829 ~1992!; 47, 4130~E! ~1993!.
16
M. Ernzerhof, J.P. Perdew, and K. Burke, in Density Functional
42
D.J. Singh and J. Ashkenazi, Phys. Rev. B 46, 11 570 ~1992!.
Theory, edited by R. Nalewajski ~Spinger-Verlag, Berlin, 1996!.
43
V. Ozolins and M. Korling, Phys. Rev. B 48, 18 304 ~1993!.
44
17
J.P. Perdew, Phys. Rev. Lett. 55, 1665 ~1985!; 55, 2370~E! C. Filippi, D.J. Singh, and C. Umrigar, Phys. Rev. B 50, 14 947
~1985!. ~1994!.
18
K. Burke, M. Ernzerhof, and J.P. Perdew ~unpublished!. 45
P. Dufek, P. Blaha, and K. Schwarz, Phys. Rev. B 50, 7279
19
K. Burke, J.P. Perdew, and D.C. Langreth, Phys. Rev. Lett. 73, ~1994!.
1283 ~1994!. 46
N.A.W. Holzwarth and Y. Zeng, Phys. Rev. B 49, 2351 ~1994!.
20
K. Burke and J.P. Perdew, Int. J. Quantum Chem. 56, 199 ~1995!. 47
P. Soderlind, O. Eriksson, J.M. Willis, and B. Johannson, Phys.
21
D.C. Langreth and J.P. Perdew, Phys. Rev. B 21, 5469 ~1980!. Rev. B 50, 7291 ~1994!.
22
Y. Wang, J.P. Perdew, J.A. Chevary, L.D. MacDonald, and S.H. 48
N. Moll, M. Bockstedte, M. Fuchs, E. Pehlke, and M. Scheffler,
Vosko, Phys. Rev. A 41, 78 ~1990!. Phys. Rev. B 52, 2550 ~1995!.
23
J.P. Perdew, in Electronic Structure of Solids 91, edited by P. 49
L. Stixrude and R.E. Cohen, Science 267, 1972 ~1995!.
Ziesche and H. Eschrig ~Akademie Verlag, Berlin, 1991!, p. 11. 50
B. Hammer and J.K. No rskov, Nature 376, 238 ~1995!.
24
A.D. Becke, J. Chem. Phys. 98, 5648 ~1993!. 51
D.R. Hamann, Phys. Rev. Lett. 76, 660 ~1996!.
25
M. Ernzerhof, Chem. Phys. Lett. ~to be published!; K. Burke, M. 52
M. Levy and J.P. Perdew, Phys. Rev. B 48, 11 638 ~1993!.
Ernzerhof, and J. P. Perdew, ibid. ~to be published!; J. P. Per- 53
A. Zupan, J.P. Perdew, K. Burke, and M. Causa, Int. J. Quantum
dew, M. Ernzerhof, and K. Burke, J. Chem. Phys. ~to be pub-
Chem. Symp. ~to be published!.
lished!. 54
26
J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson, M.R. Ped- W.E. Pickett and J.Q. Broughton, Phys. Rev. B 48, 14 859
erson, D.J. Singh, and C. Fiolhais, Phys. Rev. B 46, 6671 ~1993!.
~1992!; 48, 4978~E! ~1993!.
55
D.C. Langreth and M.J. Mehl, Phys. Rev. B 28, 1809 ~1983!.
27
A.D. Becke, Phys. Rev. A 38, 3098 ~1988!.
56
Y. Wang and J.P. Perdew, Phys. Rev. B 43, 8911 ~1991!.
57
28
B.G. Johnson, P.M.W. Gill, and J.A. Pople, J. Chem. Phys. 98, J.P. Perdew, A. Savin, and K. Burke, Phys. Rev. A 51, 4531
5612 ~1993!. ~1995!.
58
29
E. Engel and S.H. Vosko, Phys. Rev. B 47, 13 164 ~1993!. J.P. Perdew, M. Ernzerhof, K. Burke, and A. Savin, Int. J. Quan-
30
J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. ~to be tum Chem. ~to be published!.
published!. 59
J.P. Perdew and K. Burke, Int. J. Quantum Chem. 57, 309 ~1996!.