Construction of the adiabatic connection

13 December 1996

CHEMICAL PHYSICS LETTERS ELSEVIER

Chemical Physics Letters 263 (1996) 499-506

Construction of the adiabatic connection Matthias Ernzerhof Department of Physics and Quantum Theory Group, Tulane University, New Orleans, LA 70118, USA Received 15 August 1996; in final form 11 October 1996

Abstract olecule Two models for the change in the ,t-dependent exchange-correlation energy upon atomization AExc.a= ~xxt°~s _ Et~_~ xc,a are proposed, where Ex~,a = (¢AlfZ~ICA>- fd3r dSr' p(r)p(r')/21r- r'[. The wavefunction tpa yields the ground-state density p and minimizes l?+ Af'~. These models (A Exx¢.~*t) make use of the exact Ex and generalized gradient approximations (GGAs) to Ex~. The construction of the simplest model is verified by calculating the exact dAExc.A/dAIA_--0from density functional perturbation theory and comparing it to dAE~xc'~a~/dhla_-o. For systems with strong static correlation, explicit inclusion of dAEx~.a/dAIA=ofurther improves the approximation to AExc.a. Atomization energies calculated from AEm~~l show a significant improvement over GGA.

1. Introduction Hybrid schemes [1-7] which mix a fraction of exact exchange with density functional approximations (DFAs) to exchange and correlation are widely used in quantum chemical calculations. Recently, progress has been made to provide a rationale [6] for the accurate empirical hybrid scheme of Becke [1] and to construct non-empirical hybrids [6,7]. In this work, which is based on the analysis of the adiabatic connection of Ref. [ 5 ], we develop models for the change in the A-dependent exchange-correlation olecule . energy upon atomization AExc.~ = ]h-"atoms "~xc,a - Era.. xc,A The A-dependent exchange-correlation energy Exc,a is related to the exchange-correlation energy Exc by the adiabatic-connection formula

I

Exc = / d A Exc,a-

(1)

The simplest model proposed here is similar in its inputs to the non-analytic model of Ref. [ 7 ]. First we explain the physical background of our construction of the adiabatic connection. First-principles local and semi-local density functionals used in this work, i.e. the local spin-density approximation (LSD) [8] and the Perdew-Wang 91 (PW91) [9-12] generalized gradient approximation (GGA), are derived from the electron gas of slowly-varying density. In the electron gas of slowly-varying density, there is no equivalent phenomenon to the strong static correlation observed in systems such as F2 or O3. If we consider for instance atomization processes, we expect more static correlation in the molecule than in the separated atoms. An approximation which does not account for static correlation is likely to underestimate the correlation contribution to the atomization energy [ 13]. On the other hand, it has been argued [ 1-3,13-17] that an overestimation of the exchange energy contribution A E x = Ex at°ms - Ex m°lecule by local and semilocal approximations mimics static correlation. Indeed it has

0 0009-2614/96/$12.00 Copyright (~) 1996 Elsevier Science B.V. All rights reserved. PII S 0 0 0 9 - 2 6 1 4 ( 9 6 ) 0 1 2 2 5 - 0

M. Ernzerhof/Chemical Physics Letters 263 (1996) 499-506

500

been verified [ 5 ] that LSD or PW91 severely overestimate AEx for molecules with strong static correlation. This overestimation of AEx often accurately compensates for the missing static correlation contribution, so that atomization energies of systems such as CH4, NH3, and HF are well approximated by PW91 [5,13]. The overestimation of AEx by LSD or PW91 is particularly pronounced in molecules where bond formation leads to a strong spatial overlap between different molecular bonds as in CO or N2, so that these systems are overbound. In this contribution we propose the first analytic model for the A-dependent change in the exchangecorrelation energy AExc,x which is based on the above observations and which respects known constraints on AExc,x [ 18,19]. 2. The adiabatic connection formula

In Kohn-Sham theory [ 8] the ground-state energy is decomposed as E = T~ + Vext+ U + Ex¢.

(2)

Vextaccounts for the interaction energy of the electrons with the external potential and Ts denotes the kinetic energy of a non-interacting ground-state wavefunction which yields the exact density p(r). U is the Hartree energy, i.e. U = ~1 f d3r

d3rt p(r)p(r')/[r - r'[,

d3r

Tc = f da Exc,a - Exc,a=l,

(5)

o where Exc,x = Vee,x- V and Vee,,~= (,tkalf"ee[~'a). The wavefunction 0a yields the ground-state density and minimizes ~ + af,'~. Note that Exc,a---o= Ex. In practice Ec and Ex are approximated by local or semilocal density functionals, and Exc.a can be extracted from these approximations via the expression (Refs. [23,24] and Eq. ( 4 9 ) o f R e f . [25]):

Exc,x[PT,Pl] Ex[PT,Pl] + Ec,X[PT,Pl] =

= d [X2Exc[pT(r/X)/X3,pl(r/X)/k3]] " dA

(6)

For local and semilocal density functionals Ex°cFA is obtained from an expression of the form EDcFA / d3r f ( P T (r), P~ (r) ' ~7/OZ(r), ~7pj. ( r ) )

(7) EOxcF, A[pT,p~]

(3)

Vee stands for the electron-electron repulsion energy and Tc = T-Ts is the difference between the interacting and non-interacting kinetic energies. The exchangecorrelation energy is usually decomposed into an exchange and a correlation contribution Exc = Ex + Ec, where the exchange energy is given by the Fock integral Ex = - 2 ~

1

therefore

and Exc is the exchange-correlation energy, Exc = Vee - U + Tc.

and OKs denotes a Kohn-Sham orbital. The correlation energy Ec = Vee - U - Ex + T~ has a potential energy contribution and a kinetic energy contribution Tc. An important result of density functional theory shows that Tc can be transformed into a potential energy term via the coupling-constant integration or adiabatic connection formula [20-22]

d3r'

= d
(r/a), a-3~Tp~(r/a) ) }

= d{A5 fd3r

f(a-3pT(r),k-3pl(r),

A-4~TpT (r), a-4•pl ( r ) ) }.

(8)

(Jr

x °~f &/KS*(r',o-)o/KS(r,~r).

2/[r'--r I ,

(4)

The right hand side ofF_x]. (8) is obtained by variable substitution and the differentiation on the right hand side can be carried out analytically or numerically.

M. Ernzerhof/Chemical Physics Letters 263 (1996) 499-506

The importance of the coupling-constant decomposition of Exc is in part due to the fact that local and semi-local density functional approximations to Exc,a are expected to work better at a = 1 than at small A-values [ 1-3,26 ] and the coupling-constant decomposition facilitates the elimination of the error in ExOFA c,a=0 = -F. x DFA • In light of what has been said in the introduction, the success of this approach is limited to cases, where the error in AExDFA compensates for missing static correlation in A~DFA ~'~c,a=l [ 13].

3. Construction of the adiabatic connection for the atomization process Before we turn to the construction of the adiabatic connection, we discuss qualitative features of Exc,a. In finite systems such as atoms or molecules, Exc,a becomes linear in a when the electron-electron interaction can be treated as a small perturbation [27] (i.e. in the high-density limit), so that Exc,a = Ex + aEc,a=!. Thus I

Tc = / d a

Exc,a - Exc,a=! = -Ec,a=l/2.

(9)

, J

0

If correlation effects become more important (i.e. in the low-density limit), we need higher-order terms in a to account for the correlation energy which saturates as a ---+ oo [18,19]. The Ec,A curve is bent upward [5], and this implies Tc < -Ec,a=!/2.

(10)

These observations suggest that the a-dependence of Exc,a for the core electrons ( 'high-density electrons' ) is qualitatively different from the behavior of Exc,a for the valence electrons ('low-density electrons'). Thus we separate off the core electrons by modeling toms the atomization energy difference AExc,a = Ea...xc,a E~_.. °lecule The contributions from the core electrons are XC,)t * canceled in AE~c,a, Since lima--.oo Exe,a = const [ 18,19], we find lim AExc,a =const.

(11)

We introduce a [1/1]-Pad6 ansatz [28] for AExc,a which respects this constraint,

501

I1 + ab]

A iT[l/l].Pad6

~xc.a

= a L1--7-~c j .

(12)

Note that previous models for AExc,a [7,6] do not obey Eq. ( 11 ). In the following we employ this Pad6 form to represent various approximations to AExc.a, which differ from each other by the conditions used to determine the parameters a, b, and c. The parameters a and b are always fixed by the conditions AE~:I~ Pad~ AEx and AExc,a=l I1/llPad~ = AExc,:t=l,where the approximations to AEx and AExc,a=! are to be specified. This leads to [I/1]-Pad6 _ A E x + { ( 1

AExc'a

+ c)AExc,A= 1 -

-

AEx} a

1 + ac (13)

The parameter c determines the upward bending of Ab-,ll/I]-Pad6 For small c, Ab"[t/ll-Pad4 becomes --~XC,~. " --~XC,~. a straight line which joins -AEx and -AExc,a=l, providing a realistic description for atomization processes of molecules with little static correlation. For large c, [ 1 / I ] -Pad~ -AExc,a drops rapidly from its a = 0 value towards -AExc,a=l, which is a realistic behavior for atomization of molecules with strong static correlation. We can use Eq. (13) to represent aExw9'. In this case AEx = AEPW91, AExc,a=l = AEPW911 and we determine the parameter c by the condition !

da

A~PW9I

-- A F PW9I

~

--

~gC,~

--~XC

"

(14)

0

For the systems we studied so far, AExPwgl evaluated according to Eq. (6) can be accurately reproduced by the approximation Eq. (13). Note that ExPw91, as a function of A, has a derivative singularity at a = 0. This derivative singularity is smoothed in the representation o f AExPw91 via Eq. ( 13 ). In what follows A ExPw91 refers to the described Pad6 representation. As explained in the introduction, the overestimation of AExPwgl often overcompensates for the underestimation of AEe°w91 . For the coupling-constant resolved -AExPw91, that means that the a = 0 end of the curve is too negative since at a = 0 the static correlation is zero, but - - A E Pw91, which is independent of a, is still underestimated. At a ~ 1 the error i n AEP,WA91approximately cancels the error in A----x ~ PW91 , SO that ~xc.A~l A K'PW91 is expected to be accurate.

502

M. Ernzerhof/Chemical Physics Letters 263 (1996) 499-506

To eliminate the error in AExPw91 at h = 0 we replace AExPw91 by AEx evaluated from the Kohn-Sham orbitals according to Eq. (4). Since we do not know the exact AExe, to fix the parameter c we need to impose a condition different from Eq. (14) • Aib"exact is ~XC,A assumed to be identical to AEPxw91 for h ~ I, so a reasonable constraint to fix the parameter c is

- A Exc,.~

-0.

°: !t

:~

co

\

O3

"\\

-0.~

[ I/l]-Pad~

dAEx~c'a

a--J

dh

rlA ~'VW91 ~ ' a ' "t"XC, t

a=l'

(15)

-0.;

dh 0._~

which ensures that Ab-'PW91 and A~'[l/ll'Pad~ become ~XC,A ~XC,a as 'similar' as possible at h = 1, where density functionals often work best. From Eq. (15) it follows that c=

A ~-,PW91 _ A Ex LaXC,A=I

PW91 dAExrc,a /dhla=,

- 1.

(16)

With the above-defined model for the adiabatic connection, we obtain a new approximation AExl~/1]Pad~ p[~/l]Pad~ ' which to AExc, i.e. --Ab'~[l/1]-Pad~-xc= fo dh A -~xc,a [ I/l]-Pad~ facilitates a check of the quality of AExtc,a . Table 1 shows that the atomization energies (details of the calculations are compiled in the Appendix) of systems, where the bond formation leads to strongly overlapping molecular bonds, are greatly improved by the above construction. The atomization energies obtained from the [ 1/1 I-Pad6 model for the effective twoelectron bonds such as H2 and Li2 are unchanged or slightly worse compared to PW91. The LiF molecule, which is one of the most difficult molecules for hybrid methods (cf. Ref. [ 1 ] ), shows the biggest error in AEx[cI/ll-Pad6. To analyze these results we plot the ,,l-dependent AExc.a in the various approximations. The solid lines in the left and right half of the figure show -AEPxw9t for CO and 03 respectively. The dotted lines represent the proposed [ 1 / 1 ]-Pad6 model which gives accurate atomization energies, so that this curve is probably very close to the exact - A Exc,a. As already mentioned, the error in --AExpw91 compensates for the missing static correlation so that - A E ~ w91 is accurate for h 1. For A ---, 0 the static correlation goes to zero, but the error in - A E Pw9l persists. This error is eliminated by the present [ 1/ 1 ]-Pad4 model and by other hybrid schemes. In cases where the error in AEPxw91 is small, such as in H2 or Li2 (see Table 2), the PW91 approximation

-0.1 -0.4

0.2

0.4

0.6

0.8

0.2

0.4

0.6

0.8

l Fig. 1. Adiabatic connection for the atomization of CO and 03. Shown are -AEPxw91 (solid line), the [ 1/1 ]-Pad6 model (dotted line), the [ 2/21 -Pad6 model (long dashes) and - AExc,a obtained from second-order density functional perturbation theory (dashed line). Note that for CO the dotted and the long-dashed line are on-top of each other. (Energy in hartree.)

is likely to predict too small atomization energies, and we expect that AEPxW91= 1 is too small compared to the exact result• This has been verified [ 13] for the H2 molecule, where aexW321 is too small by about 14 mhartree. These findings are in disagreement with the justification of the hybrid scheme of Ref. [ 1 ] for the H2 molecule and similar systems• As already mentioned the parameter c in Eq. (13) is a measure for the static-correlation contribution to the atomization energy. For H2, N2, and O3 we get c-values from Eq. (16) of 0.6, 1.8, and 3.4 respectively. This is in agreement with our qualitative expectation. The c-values of the [ 1/ 1 I-Pad6 representation of AEPxw91 fall into a much closer range, for H2, N2, and O3 we obtain 0.6, 0.6, and 0.7 respectively• Note that the error in AEx°FA is an error in the description of the interaction between electrons with parallel spin [ 5 ]. The underestimation of the static correlation contribution to AEc is mainly due to an inadequate description of the interaction between antiparallel spin electrons which occupy the same bond orbital. Same-spin electrons are in different bond orbitals and therefore parallel-electron correlation is often less pronounced. In cases such as H2, where we

M. Ernzerhof /Chemical Physics Letters 263 (1996) 499-506

503

Table I Atomization energies (De) in the LSD, the PW91 and the [1/1 I-Pad6 approximation are compared to experimental values. The mean absolute errors (m.a.e.) are given separately for systems where the bond formation leads to weakly or strongly overlapping molecular bonds. Unrestricted Hartree-Fock (UHF) and second-order perturbation theory calculations (UMP2) based on the UHF reference states are reported. Also shown are atomization energies obtained from second-order density functional perturbation theory, either directly (GL2) or via the 12/21-Pad6 model. (All energies are in kcal/mol)

System

UHF

UM P2

GL2

LSD

PW91

[ I / 1] - Pad6

[ 2 / 2 ] -Pad6

Expt.

H2 LiH Li2 LiF Be2 CH4 NH3 OH H20 HF

84 33 3 89 7 328 200 68 154 97

104 54 16 142 -1 415 290 104 230 143

114 70 39 193 22 454 340 128 274 173

113 60 23 153 13 462 337 124 267 162

105 53 20 137 10 422 301 II0 233 143

105 52 20 131 8 419 294 106 226 138

105 54 20 133 7 418 294 106 226 138

109 58 24 139 3 419 297 107 232 141

ma.e.

48

4

28

19

3

4

4

54 90 174 115 53 33 -27 -10

64 161 266 231 152 126 96 70

190 335 355 342 265 230 407 134

89 216 299 267 199 175 242 78

77 197 269 242 171 143 185 54

70 181 257 227 155 128 150 41

69 181 257 226 155 126 146 38

91

15

133

46

18

3

2

B2

CN CO N2 NO 02 03 F2 m.a.e.

do not have overlapping molecular bonds, there is no large error in AE DFA which could compensate for the missing static correlation among the anti-parallel spin electrons.

4. Perturbation expansion of AExc,~

^ KS KS 2 I(~biKS ~jKS IV~oIG ~ >L

1

E~cL2 - - -

(jail

- ~ ia

7 6--;: 6-7:6-7

I(~b~Sl~x - 371~s)12 ' E a -- 6 i

where the indices i and j run over the occupied and a and fl over the unoccupied Kohn-Sham orbitals, ei denotes a Kohn-Sham eigenvalue and 37 is the nonlocal exchange potential of the unrestricted Fock operator formed from the Kohn-Sham orbitals. The GL2 approximation to the total energy is defined by E GL2 -- Ts

In the coupling-constant perturbation expansion derived by GiSrling and Levy (GL) [25], Ec is written as a Taylor series expansion in powers of A. The second-order term becomes

(17)

71 179 259 229 153 121 146 39

+ Vext + U "1"-E x + E~3cL2 .

(18)

Atomization energies in the GL2 approximation are given in Table 1. For comparison, second-order M¢ller-Plesset perturbation theory (UMP2) results using unrestricted Hartree-Fock reference determinants are also listed. For the systems in the upper part of the table, we observe that UMP2 gives good results and that GL2 overbinds. Note that GL2 and UMP2 both noticeably overestimate De for LiF, which indicates strong static correlation. In the lower part of Table 1, GL2 shows a pronounced overestimation of De and UMP2 still gives realistic atomization energies,

M. Ernzerhof/Chemical Physics Letters 263 (1996) 499-506

504

Table 2 Comparison of AE'xc,a__o = dAExc' a/dh]a=o obtained from AExc,a Pw9t and from AExc,a ll/ll-Pad6 with GL2 results. The first two columns show the exchange contribution to the atomization energy AEx in the PW91 approximation and exact. (Atomic units are used.) System

AEPw9!

AE~xxact

Alc,pw91 ~ ~xc,A=O

H2 LiH Li2 LiF Be2 CH4 NH3 OH H20 HF

0.043 0.049 0.009 0.195 0.028 0.307 0.200 0.071 0.174 0. I 12

0.042 0.043 0.002 0.148 0.006 0.289 0.156 0.043 0.129 0.078

0.089 0.074 0.080 0.092 0.038 0.389 0.339 0.113 0.226 0.109

0.093 0.094 0.116 0.282 0.200 0.456 0.506 0.220 0.398 0.244

0.097 0.120 0.118 0.345 0.105 0.407 0.453 0.195 0.386 0.245

B2 CN CO N2 NO 02 03 F2

0.095 0.093 0.151 0.064 0.055 0.068 0.032 0.001

0.020 -0.032 0.063 -0.054 -0.069 -0.060 -0.278 -0.112

0.105 0.224 0.200 0.265 0.221 0.179 0.319 0.091

0.681 0.854 0.596 0.790 0.820 0.881 2.216 0.895

0.570 0.844 0.592 0.742 0.701 0.654 1.548 0.566

except for 03 where static correlation is extremely strong. By taking the derivative d /daIa_-0 = d(Ex + 2aE~ L2)/da[a=0 we obtain the exact initial slope of the adiabatic connection curve: d~X,~t/dAla_--o = dE~LaZ/dAla_-o= 2E~3~L2.

(19)

Eq. (17) makes use of detailed information about the Kohn-Sham orbitals and Vx, both of which we evaluate in the PW91 approximation. Table 2 shows dAExc,a/dA[a--0, exact and in various approximations. For systems with non-overlapping molecular bonds (e.g., H2, Li2, and LiH), our model in general shows an underestimation of the slope obtained from GL2. Consistent with the underestimation of the slope, the model curve yields underestimated atomization energies. Systems with strongly overlapping molecular bonds (lower part of Table 2) show strong static correlation, and the slope predicted by the [ 1/1 ]-Pad6 model is too big. The improvement of the slope at A = 0 obtained from the proposed model compared to PW91 is significant and supports the [ 1/ 1 ]-Pad6 model. The straight lines in the Figure represent AE~xL2. For CO and to a lesser extent for 03 the

AE II/l]-padt~ t xc,a--O

Ab.,GL2 t --~xc,A~)

GL2 curves nicely approximate the tangents to the [ I / 1 ]-Pad6 curves at A = 0.

5. A model for AExc,a making use of E~cL2 In cases where we find a very big error in AExPW91 (see Table 2) and where dAExe,a/d,~la_-~ is small, the bending of the A~XC,~. tZll/ll-Pad~ curve becomes very strong and the interpolation becomes uncertain for small A. An example of such a case is 03, illustrated in the figure. In 03, which is not bound at the UHF level, the static correlation is particularly pronounced. The inclusion of the initial slope 2AE~¢L2 in our model promises an improved description of AExc,a and hence of AExc. Following the lines of the construction of the first model, we make the [2/2]-Pad6 ansatz [2/2]-Parle

AExc,a

[ I -'~ a b + ,)t2c 1

=a

1 + h d + a2eJ '

(20)

where the 5 parameters are determined by the conditions used in the [ 1/ 1 ]-Pad6 model plus the additional conditions

M. Ernzerhof/Chemical Physics Letters 263 (1996) 499-506

12/21-Pad~ dAExc.a la=0 = 2AE~cL2, dA

(21)

A2A K,[2/2]-Pad6-

0 ~< - ~ -~xc,a dA2

= min, a--0

Art2/21-pao~ is analytic for 0 ~< A ~< 1 ~XC,/[

(22)

The first condition is obvious, and the second one ensures that --A~12/2]-Pad~ follows --AE~.., L2 for small ~XC,A XC,A A without dropping below it. c12A iC,[2/2]-Pad~-

0>

--~

a.~XC,h

dA2

A---0

would imply that third-order density functional perturbation theory gives worse atomization energies than GL2. Note that the ansatz of Eq. (20) respects the constraint of Eq. ( 11 ). The long dashed curves in the Figure represent the [2/2]-Pad6 model. As expected we observe a comparatively big difference between the [1/l]-Pad6 and [2/2J-Pad6 model for the 03 molecule, where the slope at A = 0 obtained from the [ 1/ 1] -Pad6 model has the biggest error. Atomization energies obtained from AE~c2fa21"pad' (Table 1) show a small overall improvement compared to AExt~/llpad6, validating the construction of Ab-,12/2]-Pad~ Note the ~XC,A excellent result for 03, which is a known trouble case for previous hybrid schemes [ 29 ] and conventional ab initio methods (see, e.g., Ref. [30] ). The remaining error in AEI --XC2/2J-Pad6 is again attributed to the inaccuracy of AExPW921.

6. Conclusion Starting from the general features of the adiabatic connection and physical insight into the A-dependence of AExc.a, simple approximations to AExc,a have been proposed. Our models for AExc,a are based on the assumption that A E DFA mimics static correlation and on the accuracy of AED~A=I. No improvement over AEx°cFA is obtained if these criteria are not met. Recent empirical [ l ] and non-empirical [6,7] hybrid schemes are subject to similar limitations. An independent check of the present models is provided by comparison of predicted atomization energies with experimental results, and in the case of the [ 1/1 ]-Pad6 model by

505

comparison of the initial slope with the GL2 prediction for this quantity. The inclusion of E~cL2 into the model for AExc,a leads to a further improvement for systems with very strong static correlation. The non-empirical construction of AExc.a, via calculation of the GL perturbation expansion, pursued in the present contribution can be generalized to higher than second order, and represents a practical way to systematically investigate various energy differences as well as Exc,a. Studies comparing the present schemes with previous empirical [ 1] and non-empirical [6,7] hybrid schemes are in progress. The constructions of the adiabatic connection proposed here should prove useful in conjunction with any GGA functional.

Acknowledgements Helpful discussions with Professor John P. Perdew and with Dr. Michael Seidl are acknowledged. This work was supported by the Deutsche Forschungsgemeinschaft and the National Science Foundation under grant DMR95-21353.

Appendix. Technical details of the calculations The spin-unrestricted calculations are performed with a modified version of the CADPAC program [ 31 ]. Nonspherical densities and Kohn-Sham potentials have been used for open-shell atoms [32]. The experimental geometries employed in this work are listed in Ref. [33]. The De values are obtained from the experimental atomization energies and the zero point energies given in Ref. [ 34]. The Gaussian basis sets used are of triple zeta quality with two p- and one d-type polarization functions for the hydrogen and two d- and one f-type polarization functions for the second-row elements.

References [1] A.D. Becke,J. Chem. Phys. 104 (1996) 1040. 12] A.D. Becke,J. Chem. Phys. 98 (1993) 5648. [31 A.D. Becke,J. Chem. Phys. 98 (1993) 1372. 14] A. G6rling and M. Levy,unpublished. [5] M. Emzerhof,J.R Perdew and K. Burke, Int. J. Quantum Chem. Symp., acceptedfor publication.

506

M. Ernzerhof/Chemical Physics Letters 263 (1996) 499-506

[6] J.P. Perdew, M. Ernzerhof and K. Burke, J. Chem. Phys., to be published. [7] K. Burke, M. Emzerhof and J.E Perdew, submitted to Chem. Phys. Lett. 18] W. Kohn and L.J. Sham, Phys. Rev. 140 (1965) A1133. [9] J.E Perdew and Y. Wang, Phys. Rev. B 33 (1986) 8800. [10] J.P. Perdew, in: Electronic structure of solids 91, eds. P. Ziesche and H. Eschrig (Akademie Verlag, Berlin, 1991 ). [ I 1] J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson, M.R. Pederson, D.J. Singh and C. Fiolhais, Phys. Rev. B 46 (1992) 6671; 48 (1993) 4978 (E). [12] J.P. Perdew, K. Burke and Y. Wang, Phys. Rev. B, accepted for publication. [13] M. Ernzerhof, K. Burke and J.P. Perdew, in: Recent developments in density functional theory, ed. J.M. Seminario (Elsevier, Amsterdam, 1997), to be published. 1141 V. Tschinke and T. Ziegler, J. Chem. Phys. 93 (1990) 8051. 115] J.P. Perdew, A. Savin and K. Burke, Phys. Rev. A 51 (1995) 4531. [ 16] J.P. Perdew, M. Ernzerhof, K. Burke and A. Savin, Int. J. Quantum Chem. Symp., to appear. [ 171 C.W. Murray, N.C. Handy and R.D. Amos, J. Chem. Phys. 98 (1993) 7145. 118] E.H. Lieb and S. Oxford, Int. J. Quantum Chem. 19 (1981) 427. [ 19] M. Levy and J.P. Perdew, Phys. Rev. B 48 (1993) 11638. [20] D.C. Langreth and J.P. Perdew, Solid State Commun. 17 (1975) 1425. [21] D.C. Langreth and J.P. Perdew, Phys. Rev. B 15 (1977) 2884. [22] O. Gunnarsson and B.I. Lundqvist, Phys. Rev. B 13 (1976) 4274.

[23] M. Levy and J.P. Perdew, Phys. Rev. A 32 (1985) 2010. [24] M. Levy, N.H. March and N.C. Handy J. Chem. Phys. 104 (1996) 1989. [25] A. G6rling and M. Levy, Phys. Rev. B 47 (1993) 13105. [26] M. Ernzerhof, J.P. Perdew and K. Burke, in: Density functional theory, ed. R. Nalewajski ( Springer, Berlin, 1996). [27] A. G~rling and M. Levy, Phys. Rev. A 50 (1994) 196. 1281 G.A. Baker, Essentials of Pad6 approximants, (Academic Press, New York, 1975). [29] D.J. Tozer, J. Chem. Phys. 104 (1996) 4166. [301 J. Jenderek and C.M. Marian. Theor. Chim. Acta 88 (1994) 13. [311 CADPAC6: The Cambridge Analytical Derivatives Package Issue 6.0 Cambridge (1995) A suite for quantum chemistry programs developed by, R.D. Amos, with contributions from 11. Alberts, J.S. Andrews, S.M. Colwell, N.C. Handy, D. Jayatilaka, P.J. Knowles, R. Kobayashi, G.J. Laming, A.M. Lee, P.E. Moslan, C.W. Murray, P. Palmieri, J.E. Rice, E.D. Simandiras, A.J. Stone, M.-D. Su and D.J. Tozer. [32] EW. Kutzler and G.S. Painter, Phys. Rev. Lett. 59 (1987) 1285. [33] B.G. Johnson, P.M.W. Gill and J.A. Pople, J. Chem. Phys. 98 (1993) 5612. Geometry of Be: V.E. Bondybey and J.H. English, J. Chem. Phys. 80 (1984) 568; 03: M.V. Rama Krishna and K.D. Jordan, Chem. Phys. Lett. 115 (1987) 423; B2: K.P. Huber and G. Herzberg, Molecular spectra and molecular structure, Vol. IV. Constants of diatomic molecules (Van Nostrand Reinhold, New York, 1979). [34] J.A. Pople, M. Head-Gordon, D.J. Fox, K. Raghavachari and L.A. Curtiss, J. Chem. Phys. 90 (1989) 5622; For Be, 03 and BE see Ref. [33].