Elementary Concepts in Density Functional Theory

J.M. Seminario (F_Aitor) Recent Developments and Applications of Modern Density Functional Theory Theoretical and Computational Chemistry, Vol. 4 9 1996 Elsevier Science B.V. All rights reserved.

E l e m e n t a r y C o n c e p t s in D e n s i t y F u n c t i o n a l T h e o r y Mel Levy Department of Chemistry and Quantum Theory Group, Tulane University, New Orleans, Louisiana 70118 A few of the essential concepts in ground-state density functional theory are studied, and several of the basic theorems are proved in a simple manner. Then it is reviewed that coordinate scaling and the constrained-search approach provide dimensional requirements for exchange and correlation, provide formulas for the kinetic and potential components of the correlation energy, and provide a convenient path for a line integration that generates the correlation energy functional from the correlation potential. Recent research is reviewed concerning the adiabatic connection and hybrid functionals, and it is pointed out that the correction term to the two-point adiabatic connection formula consists of half the derivative of the correlation energy with respect to a coordinate scale factor in the density. This correction term is analyzed, and properties of the adiabatic connection integrand are listed for approximation purposes. 1. INTRODUCTION Density-Functional Theory (DFT) is presently enormously popular. The electron density is quite attractive to work with because it contains only three dimensions, no matter how many electrons are being considered. Superb books have recently been published that have carefully and accurately reviewed, in a systematic manner, the basic principles, theorems, and derivations in DFT. See, for instance, references [1-4]. Nevertheless, new users of DFT computer programs, who typically have little or no experience with the calculus of variations, often request to be provided, on just a few digestible pages, with quick and simple proofs of the very basic ground-state theorems in DFT, including derivations of the resultant Euler equations. Consequently, I have written the first part of this chapter with these requests in mind. For instance, I shall actually show how the corresponding Euler equation follows from the marvelous Hohenberg-Kohn variational principle [5,6], instead of just presenting the Euler equation as a consequence of the minimization process. In a sense, the first section of this contribution supplements part of the valuable earlier chapter, by Jorge Seminario [4], in a companion volume. The constrained-search formulation shall be employed throughout, even for the adiabatic connection derivation, because the constrained-search approach is simple, explicit, automatically encompasses degeneracies, and circumvents v-representability considerations.

The reader might observe that certain technical points are going to be oversimplified, even those that I myself have spent a good bit of time studying, and that I shall obviously omit many important and fascinating features of the theory in an attempt to keep the focus well-localized. Also, the reader is going to have to infer the spin-density results from the density results. In the second part of this chapter, coordinate scaling is first discussed. The coordinate scaling is then used to help quickly reveal aspects of the form of the exchange-correlation functional. Then, in the third part of the chapter, coordinate scaling is employed to display the exact formal scaling correction to the two-point approximation [7-9] in the adiabatic connection method [10-13]. Then a hybrid scheme is discussed in terms of a generalized Kohn-Sham formulation [14-16]. Finally, for approximation purposes, relations are explored [17,18] that connect the exchange-correlation potential (functional derivative) to the corresponding exchangecorrelation functional, and a line integral is investigated. In essence, the second and third parts of this chapter concern powerful techniques that are now in use for approximating exchange-correlation and its functional derivative. 2. HOHENBI~RG-KOHN VARIATIONAL T H E O R E M Suppose we want the ground-state energy and density for the following Hamiltonian of N electrons: N

)i = Zv(~) + t + 9 .

(i)

i=l

where t

-

S 1 r --_w, i-I 2

(2)

where

9.=

N-1

~

N

~

I~-~l-',

(3)

i=l j=i+l --+

and where v(r) is the electron-nuclear attraction operator of interest. Then according to the Hohenberg-Kohn theorem [5], there exists a universal variational functional, F[p], of trial electron densities, p, such that (4) p

and

Eo : fd3r vffjpoff ) + F[po ]

(5)

where E o and Po are, correspondingly, the ground-state energy E o and the groundstate density of v. What makes F[p] attractive is that p is only 3-dimensional, independent of the size of the system of interest.

2.1.

Constrained-Search

P r o o f Of T h e o r e m

It is my purpose here to review how simple the proof of equations (4) and (5) is by using the constrained-search formulation. We simply start with the familiar variational theorem:

Eo

: min <'P I~IIS'

(6)

>.

Next, the key point in the proof is the utilization of the fact that Eq. (6) may be formulated in terms of a two-step minimization. Namely [19,20] = minmin<~]I:I[~> = min<~p

Eo

p

Y~p

~p

(7)

>

p

or

Eo : man{< ~ m !~N vr p

m > + < ~'~n l t + 9 ~ l ~'.~ > }

(8)

i=l

where ~pmin is the antisymmetric wavefunction that yields p and simultaneously mlmmlzes the expectation value of H. Equation (8) concludes the proof of the theorem [expressions (4) and (5)], with 9

~

~

P

^

rain

N

(9)

rain

fd~r vCOp(O : <~L I~ vff~)l% > i--'l

with F[p]

: < ~Pp ~11:

+ 9.1

~I'p

(10)

>

and with the realization that ~mm has to be the ground-state wavefunction of I~ Po The proof is that simple! But, the following very important fact still has to be utilized. Since the lei~-hand-side of Eq. (9) is obtained as the expectation value of N

~v(fi) with any 9 that yields p, it follows that

~mm

is identified as the

P

iffil

wavefunction that yields p and minimizes the expectation value ofiust T + Vee; the v(r) term need not be included in the constrained search! Consec~uently, the F[p] given in expression (10) is truly universal because the particular electron-nuclear attraction operator of interest does not appear on the right-hand-side of expression (10). The functional in Eq. (10) follows closely the one first put down by Percus for the non-interacting problem [21]. In actual calculations, F[p] must be approximated, and good approximations to it have to take into consideration, directly or indirectly, its definition in Eq. (10). This universality is important for approximating F[p] because we could study its properties, independent of the different v's under study; the value of F depends only upon the trial p within F[p]. 2.2

Functional Derivative and Euler Equation

Functional derivatives are very important entities in density-functional theory (DFT). They appear as central characters for the minimization process in Eqs. (4-5). There are different equivalent ways to define the first functional derivative. Given

some functional called G[p], a particularly useful definition for its first functional derivative, for application in DFT, is given through aG[~ + e g] ]

I~: o

a~

= fd3r b G[p]

I g~.

8 pCO Ip=~

(11)

8G[p] I Here 8 pff) I _ is the functional derivative of G[p] at the density p=~, e is a p=p scalar, and g is an arbitrary function. By setting F = G, we are now able to derive the Euler equation corresponding to our minimization problem, expressions (4) and (5). Expressions (4) and (5) dictate that

Eo: min{fd3rv(D[poff) + eg(D] + F[Po+eg] }

(12)

for all g such that Ig(r)dar = O. Since the minimizing e is e = O, it follows that

{fdSr0e a . v~[po(-r) + eg(D] + F[Po+ eg]}l Ie=O = 0

(13)

f d3rv(f) g(f) + 0F[p o + eg] {

(14)

or

0e

: 0 .

{e:0

By expression (11), Eq. (14) becomes fd3rg(f){v(f) + 8F[p] , 8 p rp

} = 0. ~Po'

(15)

This, in turn, means that V(~

-

-

8F[p] I 8p Ip=po

(16)

because identity (15) applies for arbitrary g. Expression (16) is a fundamental Euler equation in DFT. Observe that Euler equation (16) implies that v(r) is a unique functional of the density, as was first shown by Hohenberg and Kohn via a different proof. In other words, note that Eq. (16) dictates that, except for the possible addition of a trivial constant, Po determines the one-body multiplicative attractive potential, v. Also observe that Eq. (12) requires that [22] 02F[Po+eg] I

0c2

]e=0

>0

(17)

for arbitrary g. It is anticipated that the above inequality, the convexity constraint, shall prove useful for generating approximations to F in the future.

3. DERIVATION OF KOHN-SHAM THEORY Assume that Po is simultaneously the ground-state density of the noninteracting Hamiltonian

H, :

N

t

+ ~ v~ff~)

(18)

i=l

It then follows from the Hohenberg-Kohn proof for the above corresponding interacting situation that

E~-min{fd3rv~(Dp(D+ TJp]}

: fd3rv~)po~ +TJpo]

(19)

P

where E s is the ground-state energy of I~. Here Ts[ p] = ,

(20)

and Omin .p is the antisymmetric wavefunction that ^yields p and simultaneously minimizes [21,23,19,20] just the expectation value ofT. By the previous section, the Euler equation corresponding to expression (19) is then

5Tjp] I

v, ff) -

(21)

5p(f) [P:Po

(Note that Ts[ p] is often called the non-interacting kinetic energy functional or the Kohn-Sham kinetic energy functional. Also, O p TM most often turns out to be a single determinant.) Now, partition F[p] as (22)

F[p] = Ts[ p] + U[p] + Exc[P], where 1

utpl=

ff

P(fl) P(r2) d3rld3r2

(23)

i.~1_~2 [

and where Exc[p] is the unknown exchange-correlation functional. The partitioning in Eq. (22) is a convenient one because of uniform density considerations (see, for instance, the discussion on page 18 in reference [4]), and because T s, which helps to build in the Pauli principle, can often be computed exactly. Next combine Eq. (22) with Euler equations (16) and (21), to obtain Vs(r) = v(r) + u([Po]; r) + Vxc([Po]; r),

(24)

where

po(~2)

u([po]; r-') - 58U[p]pff)lip:p0 : fd3r2 !~-~2[

and where the exchange-correlation potential, Vxc, is defined by

(25)

Vx=([Po]; 0 =

5 F_,~c[p]I 8 p('t") [p=po

(26)

This all leads to the following iterative K o h n - S h a m procedure [6]. One begins w i t h a s t a r t i n g guess for Po" C a l l this guess n r This, in turn, leads to a s t a r t i n g guess for u + Vxc, a n d hence for I-Is. T h a t is, the s t a r t i n g H s is N

I:!,[nl] = "r + P" v~([nl]; ri)

(27)

i=1

with

(28)

Vs([nl]; r) = v(r) + u([nl], r) + vxe([nl];r) where nl~2)

dSr2

(29)

and w h e r e

Vxc([nl];0 =

bE, o[p] I 5p

(30)

lp__n~

Next, solve for the ground-state of I:Is[nl]. Call this density n 2. Use n 2 for the next iteration. The process continues until n~+1 = n r In other words, the process continues until self-consistency is achieved, upon which ns+1 = Do. It is often the case that the ground-state of Hs[n k] is non-degenerate. In this situation, the ground-state wavefunction of I~[n k] emta is a single determinant. '

l~k+l

'

This single d e t e r m i n a n t consists of Kohn S h a m orbitals (pik+l t h a t are the lowest eigenfunctions of - 1/2 V 2 + vs([n k]; r). T h a t is __ ~

}

k+l,,-~

k+l

k+l~

+ Vs([nk];~) q)i tr) = ei q~i v, ;i

I~ 2 where

....

12,

N

(31)

N nk+~ = ~:

ill

I~k+ll2

Note t h a t in our notation,

(32)

bT,[p] I

~ , unless self-consistency is ~P Ip_-nk reached and Po is obtained. At self-consistency, n k = Po, and thus in our notation w e

have

v~([Po];r-3 : v~Ct) : -

v,([n~]; it) * -

~T,[P] I 8p~l

P=Po

4. C O O R D I N A T E SCALING The exchange-correlation functional must be approximated in density-functional theory. In order to help one understand its form and in order to aid with our discussion later about the adiabatic connection method, some key coordinate scaling properties of Exc[p] shall now be briefly reviewed.

4.1

Definition Of Exchange And Correlation

To discuss coordinate scaling, it is first convenient and useful to partition Exc[p] into its exchange and correlation components (33)

Exc[p] = Ex[ p] + Ec[P] with E~[p] = < r

rain

rain

19~ECp

(34)

> - U[pl

so that Eo[pl = < e 0 I t §

4.2

>-

§

(35)

Coordinate Scaling Requirements For Ex[ p] And Ec[ p]

The study of coordinate scaling is the study of the exact dimensional properties in density-functional theory. Coordinate scaling provides key requirements of E x and E c for approximation purposes. Moreover, as we shall see in section 5, comprehension of the adiabatic connection method is greatly enhanced through coordinate scaling. With these thoughts in mind, define the uniformly scaled electron density, px(r), by px(r) - px(x,y,z) - ~3p()~x,~y,~z)

(36)

With definition (36), we know that E x satisfies the following simple homogeneous scaling relationship [17]: Ex[Px] = ~.Ex[P]

(37)

Requirement (37) dictates that a dimensionless function must multiply p4/3(r) in the integrand for a viable approximation to Ex[P]. The 4/3 results from (1+3)/3, where the number 1 stems from the fact that ~. is raised to the power of unity in Eq. (37), and where the number 3 stems from the fact that space is 3 dimensional [24]. The scaling for E c is much more complicated than for exchange. This fact is reflected in the forms of present approximate functionals. For instance [17], Ec[px] > ~Ec[p]; ~. > 1,

(38)

and [14,25,26] Ec[P~] : ~c,2[P] + ~-lc~,s[P] + ~-2~,4[P] + --for sufficiently high X, so that [25]

(39)

10 lim E [ p x ]

.

(40)

Ec[px] = ~D[p] + "-"

(41)

> -oo

While [251

for ~. close enough to zero, where D[p] is an unknown functional. In other words, E c transforms differently for high and low ~.. This means that unlike in the E x case, the integrand, for a viable approximation to Ec, does not simply consist of a dimensionless function which multiplies pro, where m is some constant. It has been shown [17] that the simple scaling for E x arises from the fact that the constrained-search definition for (I)~fin, in Eq. (34), involves just the operator q?, in the search for Cram, and from the fact that V is homogeneous of degree -1, through Irir~j 1-1. Consequently ~. is raised to the power of unity in expression (37). On the other hand, the complicated scaling for E c arises from the fact that the constrained-search definition of^ tI/m]n in Eq (35) involves two operators with different homogeneous p properties; T is homogeneous of degree -2, because it contains second derivatives, while Vee is homogeneous of degree -1 as stated above. (Eq (38) shall be derived later to illustrate a scaling derivation.) Equation (39) is equivalent, at ~ = 1, to E[p] : 2] ~ [ p ]

(42)

i:2

where ec,i is homogeneous of degree i-2. That is [27-29] ec,i[px] = ~i-2 ec,i[Pl.

(43)

Moreover, e c 2[P] is known exactly as the second-order contribution [14] in a recently developed D~T perturbation theory [14,25]. In this perturbation theory, the noninteracting Kohn-Sham single determinant of p is the zero-order wavefunction and

9,=-

N

z

N

(44)

v ([pl;g~)- ~u([p];~)

i--1

i--1

is the perturbation [14] for the second-order term. 4.3

Kinetic And Potential

Contributions

To The Correlation

Energy

Coordinate scaling allows one to generate the kinetic energy, Tc, and the potential energy, Vc, contributions to the correlation energy. The definitions of these contributions are, respectively, Tr

mm

: <~f'p

mm

Itl~p > -
n~n

Itl|

~

(45)

,

and rain

rain

They are obtained through coordinate scaling by [171

(46)

ll

Tr

=-EJo]

+

OEJo3 1 aX IX=I

(47)

OE~[P~]I

(48)

and [17] V~[p] : 2E~[p] -

ax

Ix=l

Expressions (47) and (48) shall be derived in section 7.1.

5. ADIABATIC C O N N E C T I O N M E T H O D The adiabatic connection formula is indeed a powerful one because it allows one to just focus upon the electron-electron repulsion operator when approximating Exc[P]. In this section, the adiabatic connection formula [10-13] and method shall be derived and discussed. The object of the adiabatic connection method is to approximate Exc[p] by approximating the integrand in the adiabatic connection formula. The study of the adiabatic connection method is presently a very active area of research and is now beginning to provide exchange-correlation functionals that are among the best available. 5.1

Derivation

Of Adiabatic

Connection

Formula

For the purpose of quickly deriving the adiabatic connection formula [10-12] through the simple constrained-search orientation [13], let's now define tlJpmin'~^as the antisymmetric wavefunction that yields p and simultaneously minimizes
(49)

it follows that 1

E~r

+U[p] = f dX a-; O , -" It+xg,,,Iv,,

(50)

0

where utilization, in Eq. (50), has been made of the identities mi~ mi~,i ~ mi~,0 (51) ~o = ~ p and Cp = ~ o " Next, we take advantage of the minimizing (stationary) quality embodied in the definition of ~Ijmin')', concerning the variation of ~. at fixed p, to ascertain that in Eq. (50) we only l~ave to differentiate with respect to the operator, ~"§ ~'~ree ; the derivative with respect to ~min,~, at fixed operator, is zero. Consequently --p

12 1

Er

+ U[p] :

fd; < 'Fp

min,~

min,~

{'q~l'l'p

>

(52)

0 or 1

Exr

=f

<53)

Vxc[p]d;~

0

where x

YAp]

: < v

l%.lv

mi~,X> U[p] -

9

(54)

The partitioning of Exc into E x and E c in expressions (33)-(35) enables a separate adiabatic connection formula for Ec[P]. Namely, 1

(55)

EJp]= fvSp]ax 0

where Vc~[p] = <~iJpm,k e "i~re i~iJpmink_,>_

(56)

5.2

Exact P r o p e r t i e s Of Adiabatic C o n n e c t i o n I n t e g r a n d The object in the adiabatic connection method is to approximate its integrand Vx~cas a function of ~. See Eq. (53). Knowledge of V~xcas a function of ~ is equivalent to having exact knowledge of Exc[P]. Accordingly, for approximation purposes in mind, I now list exact properties of Vx~c[p]. Note that the superscript ~ represents an adiabatic connection parameter, while ~ is a scale factor when it is a subscript with p. Let's commence by partitioning Vx~c[p] as x Vxr

(57)

= Ex[p] +V#[p]. ,

where definitions (34), (54),and (56) have been employed. Hence, we need only focus upon how Vcx behaves in order to understand how VxXc behaves, because Ex[P! separates out of VxXc. With this in mind the basic uniform scaling requirement of V c is [17,25] Vc [P] =

xv '[p 1/•1

=

x Vc[Pl/• ]

,

(58)

where pl~(x,y,z) = ~,-3p(~'1x,%-ly,~.'1z).

Next, let'sfollow Vc~ along the ~ path. Begin at ~.=0. W e know the result exactly at this point. It is Vck=~

= 0.

(59)

13 Now let's look at ~ slightly greater than zero. In this region Vcx is linear [9] in ~. To see this begin with [17,14]

V~[pl - 2xEe[pl/x]

+ ~2 0Ee[Pl/x]

ax

(60)

and for Ec[Pl/x] substitute [14,25] Ec[Pl/x] = ec,2[p] + ~ec,3[p] + ...

(61)

which applies when ~ is sufficiently small (See Eq. (39)). The result is [14,16,30] VcX[p] = 2~ ec,2[p] + 3~2ec,3[p] + ...

(62)

Now let's investigate Vcx for ~. greater than the region near ~=0. For this purpose, an important general property is the fact that V c i s a monotonically decreasing function of ~. That is [17] OX

Vc [p] < 0, all X > 0.

(63)

Next, although one integrates only up to ~.=1 in expressions (53) and (55), it is useful to give a result for ~-~oo. Accordingly, we exploit the fact that < ~p

' 19.1 ~p

>>0

(64)

to imply [25] V c [p] ~ -U[p] - Ex[P]

(65)

which means that V~ is bounded as )~ ~ ~,. In fact, through implementation of the Lieb-Oxford bound [311 it has been proven that [32,221 ~-.oo

v:

f

d3r - Ex[p].

(66)

Observe that the coefficient of the linear term, 2ec,2[p] , is twice the second-order contribution to the correlation energy [14,30]. That VcX[p] is linear in the region of small enough ~ was utilized in reference [9]. Finally, expressions (59), (62), (63), and (64) suggest that perhaps the derivative of Vcx is a monotonically increasing function of $.. In other words, perhaps O2VcX[ol > 0, all X > 0 8X2

(67)

It would be worthwhile to investigate the validity of expression (67). 5.3

E x a c t C o r r e c t i o n T e r m To L i n e a r A p p r o x i m a t i o n In A d i a b a t i c Connection Method: It is the purpose here to deduce the exact expression for Exc[ p] in terms of the linear approximation to the adiabatic connection integrand, Vxcx, and a scaling correction. Start with

14

1

Er

0Er

1

= ~Vctp] + ~

I

(68)

0------~lX=l

which is a r e a r r a n g e m e n t of Eq. (48). Then simply add Ex[ p] to both sides of Eq. (68) to obtain the desired result which is [16] 1

o

1

1

1

E~o[p] = ~ V~J0] + ~ V,o[0]. + -

2

OEr

I

0X

(69)

I~.=l

or, equivalently, Exr

1 E [p]

=

2

x

+

1 2

~ Vxr

1 OEr +- ~ 2

O•

I

(70)

[~.=1

because vxO[p] = Ex[p ]

(71)

by Eq. (54), and Vx'c[p] = Ex[p] + V[p]

(72)

by Eq. (57). The first two terms on the right-hand-side of Eq. (70) constitute the linear approximation for Exc within the adiabatic connection method. Finally, an alternative way of arriving at expressions (69) and (70) is given in section 9. The last (scaling) term in equations (69-70) constitutes the correction to the linear approximation (the assumption t h a t Vxcx is linear in ~..). Since we have derived Eqs. (69-70) by starting with relation 68, it follows that the linear approximation would be exact if the density-functional definition of the correlation energy, Ec[P], and the density-functional definition of the potential energy contribution to the correlation energy, Vc[P] , were to satisfy the same virial relation, for atoms, as are satisfied by their more familiar standard q u a n t u m chemistry counterparts. In other words, the error in the linear approximation would be zero if

E[p] = ~1 Vr

.

(See also

references [7] and [9].) Alternatively, the linear approximation would be exact if the perturbation expansion for Ec[Px] , Eq. (39), where to contain only the second order term, because 0Eo[p~] I

0X

[X=I

= - er

p] - 2 er

....

(73)

by expressions (39), (42), and (43). Consequently, the two-point (linear) formula with the added higher-order perturbation correction is [16]

15

= 5.4

1

+

1

~

V;c[p] -

An Analysis Of Considerations

The

1

[p] - e 4 [ p ] . . . . Adiabatic

Connection

(74) Method

And

Future

With expression (70) we can understand why immediately the linear (two-point) approximation offers improvement over the local density approximation (LDA) for the exchange-correlation energy when the LDA repulsion contribution to exchangecorrelation is employed at a coupling-constant of unity but not at the couplingconstant of zero. Although Eq. (70) was not published at the time of Becke's work [8], Becke, in effect approximates the right-hand-side of Eq. (70) by neglecting

1 OEo[Px] I ~ and by using 2 Ox Iz--1

-

1 E,[p] + I ,,,I v,; LDArtpl

(75)

where E x is, of course, the exact exchange energy and where Vxlc'LDA is the local density approximation for Vxlc. Analogous to Eqs. (72) and (48), VxI'LDA is obtained from the LDA exchange and correlation functionals by

v,,~ ~ A [p] =F. xLDA[p] § 2 E ~LDA[p ] -

oEf'[p~]

ox

t [X:l

(76)

Equation (76) is equivalent to the LDA expression in reference [33]. Expression (75) typically improves upon the ExLcDA because, relative to the other terms, the neglected exact value of

1 OEJp~] I 2 ox Ix:l

-

is generally small in magnitude and

because

_1 --xcV 1,LDA is generally closer to _1 V:c than is -xc ~' LDA to Exc. Indeed, it has 2 2 recently been argued that the LDA describes the exchange-correlation hole better at the coupling-constant of unity than at the coupling-constant of zero. Exact values for

1 aEo[p~] ]

-2 ~ ak

Ik=l

have recently been obtained.

The obtained

values, in atomic units, for this derivative correction to the two-point formula are small in magnitude. They are -0.003, -0.012, and -0.032, for the He, Be, and Ne atoms [34-36], respectively, while the corresponding total values for Exc[p] are approximately -1.068, -2.76, and -12.50. (Incidentally, the corresponding LDA values for 1/2 the above scaling derivative are much too negative. They are, respectively, -0.045, -0.086, and -0.247. The Perdew-Wang (91) values [37], however, are very good here. They are, respectively, -0.004, -0.011, and-0.034).

16 1 When the exact Vxcand E x are employed, the two-point approximation,

1 Ex + 1 1 2 2 vx~

(77)

will yield atomic estimates for the exchange-correlation energies that are generally

1 too positive because the signs of ~

aEjo ] I OX

[)~=1 have always been found to be

negative. That is 1 Ex[p] +

1

1 Vxjp ]

(78)

is to be expected. Expression (78) is consistent with inequality in conjecture (67). Even if one were to know the exact vie[p], then the above numbers for the exact scaling derivative correction indicate that the linear (two-point) approximation might not really be quite accurate enough. So, it looks like one has to employ, directly or indirectly, the scaling derivative term. With this in mind, I think that one of many possibilities to try would be to use mod

modl 1 E , , + 1 Vx~ ' + - 10E~ [p~]] 2 2 2 0~. I~= 1

(79a)

where the superscript "mod" indicates that an approximate (model) functional is used for the term. It would even be worthwhile to try an approximation for the scaling derivative term which comes from a model that is different from the one used for Vxcm~ Along these lines, it has been found that the accuracy becomes spectacular when the exact functional is used for the third term as well as for the first term in Eq. (79a), and the PW-91 is used for the second term [38]. There are obviously countless possiblities and opportunities for investigation [9,38]. Finally, it is interesting that if the model functionals used in the last two terms of expression (79a) are the same, then the expression becomes 1 Ex[p] + 1 Exmod[p] + Errood [p]

(79b)

2

and a scaling derivative thus disappears. It is consequently noteworth that until quite recently, all the hybrid functionals constructed after the publication of reference [8] have not contained a scaling derivative term.

5.5

G e n e r a t i n g A p p r o x i m a t i o n s For Vxc~ The object of course is to appproximate as best as possible the adiabatic connection integrand, Vxc~, from ~=0 to ~,=1, in expression (53). For this purpose the exact properties listed in section 5.2, via equations (57) through (67), should be of great help. For instance, we know that Vxc~ is E x at ~=0 and we know through Eq. (62) that its slope is 2e c 2 at ~=0. Both E x and e c 2 can be evaluated through the Kohn-Sham orbitals and orbital energies. (In fact encouraging implementation of this knowledge of ec,2 has recently been made by Ernzerhof [30]). The rest of the constraints listed in the preceding section deal with the shape of the exact Vxc~ curve and its value at ~=oo. By employing their knowledge of the known properties of the shape of V xc, x Burke, Ernzerhof, and Perdew have fruitfully been exploring new approximate

17 functionals [30, 39-41]. In approximating the Vxcx integrand one needs, of course, to approximate Vxcx for at least one point other than for ;L=0. Accordingly, it is emphasized here that one may generate a model Vxcx from a model Exc through use of expressions analogous to (57) and (60). That is [17,14,9] mod

V m~

- g m~

+ 2 X grrood[~ ll;t] + ~2

aEr OX[P~/x]

(80)

In fact, for rood = LDA, Becke [8] employed an expression from the literature [88] that is equivalent to the fight-hand-side in the above equation. In Becke's second paper, however, a Vxcm~ was not employed when he investigated gradient functionals, perhaps because expression (80) is not all that well known. I feel that it is worthwhile to pursue the study of Vxcx with respect to gradient funetionals beyond the LDA. In fact, one could actually employ a different approximate functional (model) for each term on the fight-hand-side of Eq. (80). Essentially all of the published prevalent hybrid schemes employ neither a discussion of Vx~x nor an approximation to it in the hybrid formula. Again, p e r h a p s this is because expression (80) is not well known. Instead, a hybrid functional of the following form is usually presented [42]" hybrid

E~

.=mod

[p] = aE~[p] +(1-a)~ x

--mod

(81)

+Er

where, consistent with our notation throughout this chapter, E x is the exact exchange -rood and Er -rood are approximations. The bar on each functional energy, but both E~ indicates that it may consist of a linear combination of present approximations.

6. H Y B R I D M E T H O D S A S G E N E R A L I Z E D K O H N - S H A M T H E O R I E S

In a class of present hybrid schemes, a part of the exact Fock potential is employed, and the rest of the functional derivative is approximated by some linear combination of a variety of present approximations for exchange and correlation. In order to improve these hybrid approximations, it is necessary, through formal construction, to define what has to be approximated. With this in mind, I now review here the proof of this formal construction [14-16]. Define F~IF by

F~[p] = < r

I~ +

~r162

>

(82)

where q)~F[P] is the single determinant that yields p and minimizes . Here HF signifies Hartree-Fock. Observe that (I)~tF[PHF] is the familiar Hartree-Fock single determinant for the Hamiltonian in Eq. (1), where PIlE is the corresponding Hartree-Fock density. Next define v a such that

min{fd~rv~p~

+ F~[p]} =

F rpj

s3)

P

where Po, our physical interacting density of interest, is simultaneously the groundstate density of H in Eq. (1). (Notice that Po ~ PHF) The Euler equation

18 corresponding to (83) is v~(f) :

8F~[P]I 8p(f) Ip:po

(84)

Now utilize F~F as a component of F. functional G a such that

In other words, define an unknown (85)

F[p] : F~F[ p] + Ga[p]. That is, Ga[P] : -
IT + aVeelr

>

(86)

Euler expression (16) thereby becomes vff)

:

- 8P~[p]

- g~,([po];r-'),

[

(87)

8pff) I p:po where

8O.[p] I

g~([Po]; r-') :

8pff)

[p:po

.

(88)

By comparing (84) and (87), we have va(r) = v(r) + ga([Po]; r).

4(89)

Hence, Po, the physical interacting ground-state density of v(r), may be obtained through the following Hartree-Fock calculation:

N rain < CSDIZ {vff~) + g~([pfl;ft)} + t + ag.~ [ r

@sD

i=l

> '

(90)

where the CSD are single determinants. In other words, the ground-state density is obtained as the Hartree-Fock density for the attractive potential v(F) + ga([Po]; F) and the repulsion operator a~;ee. The minimizing single determinant, r is thus obtained from those orbitals, ~?, that solve the following corresponding system of Hartree-Fock equations"

I_ 1V2+v(f)+ga([Po];~+au([Po];~+aVx^ l'IF,,f al r-')1 {~k(~ i3@i/; 2

(91)

= %~k(f); k=l,2,...,N The operator in expression (91) is the usual Hartree-Fock operator, except that it contains the coupling parameter a and the multiplication potential ga- The operator VxHF is the familiar non-local Fock exchange operator. Here it is composed of the

19 orbitals r The density in ga([p];f) and in u([p];r), and the orbtials in ^v xHF , are changed during each interation until self-consistency is reached. Upon selfconsistency, q)~F[Po] is generated. The physical ground-state energy is then given by Eo =
11:I+ (a-1)Vee[r

+ Ga[Po]"

(92)

or

Eo = ~ d3r v(F)Po(F) + H~[Po]> + Ga[Po]

(93)

In a sense expression (93) is a reasonably direct justification for the class of approximations commonly employed on the right-hand-side of Eq. (81). To help best appreciate this section, let's investigate the a=O and a=l situations. At a = 0, v + ga becomes the Kohn-Sham effective potential v s and ~ F becomes the Kohn-Sham single determinant. At a=l, Ga[Po] becomes a correlation functional that should be typically close in value to Ec[Po], with the inequality Gl[Po ] > Ec[Po]

(94)

dictated by the definitions of E c and G 1. For al__!a, we have e~ = -I,

(95)

where e~ is the highest-occupied orbital energy and I is the exact ionization energy. The value of e~ is independent of the value of a because the fact that the groundstate density, Po, is independent of a implies that the asymptotic decay of Po is independent of a, and both I and e~ depend only upon the asymptotic decay of the density (see, for instance, references [1] and [2] and the works cited within).

7. C O O R D I N A T E SCALING P R O O F S

The reader will gain some understanding of the simplicity of the constrainedsearch orientation and coordinate scaling if we now go through some of the scaling proofs [17]. They are indeed straightforward. Let's begin by immediately observing that

F[px] _<<)~3N/2 ~pmm " (~rl,-..~rN)IT + ~Tee]~3N/2~Ijpmin(~rl,'')~rN)>

(96)

because )~3N/2~IJpmin()~l,...,~,r N) --->pk(r)

(97)

and because of the <~? + ~Tee> minimizing quality of ~_in. Since T is homogeneous of degree -2, since Vee is homogeneous of degree -1, and ~since it is known that Ts[Px] = ~2Ts[P]

(98a)

20 and U[px] = ~.U[p],

(98b)

Eq. (96) becomes ~.2Ts[p] + )~U[p] + )~Ex[p] + Ec[Px] < ~.2T[p] + kVee[ p] ,

(99)

where T[p] : : Ts[P ] + Tc[O],

(100)

and where

Vee[ p] = = V[p] + Ex[ p] + Ec[ p] - Tc[Pl. Next, employ expressions (100) and (101) in combination with Eq. (99). number of cancellations, the result is Ec[Ox] + )~(1-~.) Tc[ p] - kEc[ p] < 0,

(101) After a

(102)

which implies inequality (38) because To[p] > 0. S c a l i n g D e r i v a t i o n O f T c A n d Vxcx To get expression (47) for Tc, utilize the fact that the left-hand-side of Eq. (102) achieves its maximum at )~ = 1. Consequently, Eq. (47) is generated [17] by setting the derivative equal to zero, at )~ = 1, of the lei~-hand-side of Eq. (102). Further, since the T c expression implies the expression for Vc, Eq. (48), we have also derived the latter. Finally, combine Eq. (48) with 7.1

to arrive at [17,9] Vr [p] = 2XEr

+ X2 0Er 0X or, by adding E x to both sides,

(104)

Vxc[p] = 2~.Exc[Pllx] + ~.2 0Exc[Plp.] O~

(105)

7.2

D e r i v a t i o n VcX[p] -- ~.Vc[pm] Let's now prove the important identity (103), which is a restatement of relation (58). To achieve this, let's first review that it has been shown that ~.3N/2tI~pmin(~l...,~r~N) yields px(r') and minimizes [17]

(106)

21 Consequently, it follows that 9

"'"

(107)

)~3N/2~min'l(~'f "'" )~fN)" 1,

,

In other words, the right-hand-side of Eq. (107) yields p(x,y,z) = ~3pyx()~x,~.y,~z)

(108)

and minimizes the operator in expression (106). This infers that = < x3N/2~min, l(X~p,/x x l""'XfN) 19~(fl .....fs) IX3N/Z'l'm ' 'm' ' pZ,(Xfl /x

V c][~Pk

., XfN)> -C[p] -E~[p] (109)

or Vc~[p] = ~,Vee[Pl/X] - V [ p ] - Ex[ p]

(il0)

where Vee[Pl/~.] = <

~ptl~(-~l,...,fN)[gee(T1,...,~N)] ~pllx(-~l,...,fN)

>.

(111)

As earlier in this chapter, use has been made here that Vee is homogeneous of degree -1. Finally, identity (58) is arrived at [17,25] when U[px] = ~U[p] is employed, etc.

8. Exc FROM Vxc BY LINE INTEGRATION THROUGH SCALING PATH Knowledge of the exact relationship between Exc and Vxc is especially important now that an active area of research involves the direct approximation of Vxc, in order to obtain more accurate orbital energies and densitites. The relation for exchange is simply [17] Ex[ p]

= -Sd3rp(~)~.Vvx([p]; ~)

(ii2)

The situation is more complicated for correlation. Instead, for correlation Eq. (112) transforms into the following inequality (113)

Ec[P] < -~p(~)~'Vvc([p]; ~)

However, all is not lost. We could get E c from Vc, provided that we take a line integral along some specified path. For instance [18]

Ec[P] : foX dXfd3rvr

rD ap~ff)

(114)

where

I vc([P~'];r-'): 8 Ejp] 8p Ip:px

(115)

Eq. (114) above, of van Leeuwen and Baerends, results from a line integration of a formula of Ghosh and Parr [43] or a formula of Levy and Perdew [17]. It is

22 instructive to now derive the formula. Commence with Ec[Px,] - Ec[P x] = S d3r Vc([px]; r~)[px,(~) - px(~)] + higher order terms in Px' " Px,

(116)

where )J is very close to )~. Next divide both sides by 2~' - 2~ and let )J ~ 2~. Obtain aE~[p~]

0x

0 p~ (t")

fd3rvr

:

(117) '

which resembles more closely the representation of the corresponding expression in reference [43] t h a n in reference [17]. Next integrate both sides of Eq. (117) from 2~=0 to ~=1 to generate our desired result, expression (114), by incorporating the crucial end point identity [25] Ec[Px=0] = 0

(118)

which is consistent with Eq. (41). I feel t h a t relations (112) and (114) will see ever increasing roles in the future. The object is to model v x and v c directly by utilizing important aspects of their known properties. Then one uses relations (112) and (114) to obtain, respectively, E x and E c from v x and v c.

9. C L O S I N G R E M A R K S Various parts of the development within are closely connected to other studies. See, for instance, references [44] and [45]. Moreover, the exact correction term (the scaling derivative) to the linear approximation in the adiabatic connection method can be obtained by a straightforward manipulation of Eq. (3) in reference [38]. This latter equation is E~c[p] = ~

pff)f'Vv, c([p];~]

(119)

where our present notation is used. Specifically, to obtain equations (69) and (70), in Eq. (119) partition Vxc as vr

= v,([p];~ + vr

(120)

utilize Eq. (112), which relates E x to Vx, and utilize [17] 0Er 0X

[ = - fd3r pCf)f'Wr IX=l

(121)

which is an alternative way of writing the corresponding expression in reference [43]. Equation (119) was employed to obtain extremely accurate values for Exc[ p] from experimental densities [38] with the generalized gradient approximation [32,37]. Advantage was taken of the fact t h a t the exact exchange-correlation energy functional is more local for full-coupling strength t h a n for the coupling-constant average [46].

23 A number of valuable scaling relations are known for Ec, as reviewed in this chapter. Nevertheless, it would be absolutely wonderfull if one were to find a general scaling identity for E c that is comparable to the one for Ex, but would of course have to be more complicated. Such an identity would indeed impose a powerful constraint for approximating E c. Towards this goal, combine expressions (39), (40), and (42), to obtain E[p] = limEc[p~ ] ~

ax

I~=1

- 2ec4[P] . . . . . '

(122)

Thus, Eq. (122) pro~des this ]dnd of identity, through third order, for those densities where the Taylor series is valid. Further study is of course needed.

REFERENCES 1. R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules (Oxford University Press, Oxford, 1989). 2. R.M. Dreizler and E. K. U. Gross, Density Functional Theory (Springer Verlag, Berlin, 1990). 3. Density Functional Theory, edited by E. K. U. Gross and R. M. Dreizler (Plenum, New York, 1995). 4. Modern Density Functional Theory: A Tool for Chemistry, edited by J. M. Seminario and P. Politzer (Elsevier, Amsterdam, 1995). 5. P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 6. W. Kohn and L. J. Sham, Phys. Rev. 140, A 1133 (1965). 7. M. Levy in Density Functional Theory, edited by J. Keller and J. L. Gasquez (Springer, New York, 1983). 8. A.D. Becke, J. Chem. Phys. 98, 1372 (1993). 9. M. Levy, N. H. March, and N. C. Handy, J. Chem. Phys. 104, 1989 (1996). 10. D. C. Langreth and J. P. Perdew, Solid State Commun. ~ 1425 (1975). 11. O. Gunnarsson and B. I. Lundqvist, Phys. Rev. B13, 4274 (1976). 12. D. C. Langreth and J. P. Perdew, Phys. Rev. B15, 2884 (1976). 13. A. G5rling, M. Levy, and J. P. Perdew, Phys. Rev. B47, 1167 (1993). 14. A. GSrling and M. Levy, Phys. Rev. B47, 13105 (1993). 15. A. Seidl, A. GSrling, P. Vogl, and J. A. Majewski, and M. Levy, Phys. Rev. B5__33, 3764 (1996). 16. A. GSrling and M. Levy, J. Chem. Phys., in press. 17. M. Levy and J. P. Perdew, Phys. Rev. A32, 2010 (1985). 18. R. van Leeuwen and E. J. Baerends, Phys. Rev. A51, 170 (1995). 19. M. Levy, Proc. Natl. Acad. Sci. USA ~ 6062 (1979). 20. M. Levy, Phys. Rev. A26, 1200 (1982). 21. J. K. Percus, Int. J. Quantum Chem. 13, 89 (1978). 22. M. Levy and J. P. Perdew, Phys. Rev. B48, 11638 (1993). 23. M. Levy, Bull. Am. Phys. Soc. ~ 626 (1979). 24. L. J. Sham, Phys. Rev. AI_, 969 (1970). 25. M. Levy, Phys. Rev. A4__3~34637(1991). 26. A. GSrling and M. Levy, Phys. Rev. A45, 1509 (1992). 27. A. Gfrling and M. Levy, Int. J. Quantum Chem. Symposium ~ 93 (1995).

24 28. 29. 30. 31. 32. 33. 34.

35. 36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46.

A. GSrling and M. Levy, Phys. Rev. A5__2_2,4493 (1995). S. Liu and R. G. Parr, Phys. Rev. A, ~ 2211 (1996). M. Ernzerhof, submitted to Chem. Phys. Letters. E. H. Lieb and S. Oxford, Int. J. Quantum Chem. 19, 427 (1981). J. P. Perdew, in Electronic Structure of Solids '91, edited by P. Ziesche and H. Eschrig (Akademie Verlag, Berlin, 1991). J. P. Perdew and W. Yang, Phys. Rev. B45, 13244 (1992), and references within. C. J. Umrigar and X. Gonze, in High Performance Computing and its Application to the Physical Sciences, Proceedings of the Mardi Gras '93 Conferences, edited by D. A. Browne et.al. (World Scientific, Singapore, 1993). C. J. Umrigar, private communication. C. J. Umrigar and X. Gonze, Phys. Rev. A5_0~03827 (1994). J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Phys. Rev. B ~ 6671 (1992); 48, 4978(E) 1993. K. Burke, J. P. Perdew, and M. Levy, Phys. Rev. A53, R2915 (1996). M. Ernzerhof, K. Burke, and J. P. Perdew, submitted to the International Journal of Quantum Chemistry Symposium Issue. K. Burke, M. Ernzerhof, and J. P. Perdew, submitted to Chem. Phys. Letters. J. P. Perdew, K. Burke, and M. Ernzerhof, submitted to J. Chem. Phys. A. D. Becke, J. Chem. Phys. 98, 5648 (1993). S. K. Ghosh and R. G. Parr, J. Chem. Phys. 82, 3307 (1985). E. I. Proynov, E. Ruiz, A. Vela, and D. R. Salahub, Int. J. Quantum Chem. $29, 61 (1995). O. V. Gritsenko, R. van Leeuwen, and E. J. Baerends, Int. J. Quantum Chem., to be published. K. Burke, J. P. Perdew, and M. Ernzerhof, Int. J. Quantum Chem., to be published.