Nonlocal correlation functional involving the Laplacian of the density

2 December 1994

ELSEVIER

CHEMICAL PHYSICS LETTERS

Chemical Physics Letters 230 ( 1994) 419-428

Nonlocal correlation functional involving the Laplacian of the density E.I. Proynov

‘, A. Vela 2, D.R. Salahub

D@artement de Chimie et Centre d’Excellence sur la Dynamique MolPculaire et Interfaciale, C’niversitkde Montreal, MontrPal, Quebec, Canada H3C 3J7 Received 3 1 May 1994; in final form 10 August 1994

Abstract The gradient-free correlation functional derived in a previous work is extended to a nonlocal scheme beyond the gradientcorrected ones. It involves the Laplacian of the density and the self-consistent (SCF) Kohn-Sham kinetic energy density. Using two additional parameters, the correlation contribution to the binding energy is adjusted to compensate the error from a chosen exchange functional. When combined with the exchange functionals of Becke and Perdew the new scheme gives improved results for the binding energies and geometries of test molecules.

In previous works [ 1,2 ] we derived and tested a correlation functional implemented in the linear combination of Gaussian-type orbital (LCGTO) code deMon [ 31. The derivation was based on a model pair-correlation function, using the adiabatic connection formula of the Kohn-Sham (KS) density functional theory (DFT). Although derivatives of the electron density were not employed, the new correlation functional yielded more accurate binding energies and equilibrium geometries than the local density approximation (LSD, correlation of Vosko, Wilk and Nusair [4] ) when it was combined either with the standard LSD exchange, or with a modified gradientless exchange [ 21. For many of the molecules considered, the results were even comparable to those of the gradient-corrected functionals currently in wide use: Becke-Perdew (BP) [ 5,6 1, Perdew-Perdew

(PP) [6,7] and Perdew-Wang (PW91) [8,9]. However, for some ‘difficult’ molecules like FZ, O2 and C2H4, the binding energies were still considerably overestimated, albeit to a lesser extent than for LSD [ 21. Even the nonlocal schemes of the generalized gradient approximation (GGA) type, such as BP and PW9 1, meet difficulties in these cases [ 21. In the present work the correlation functional derived in Ref. [ 1 ] is extended to a nonlocal scheme involving the Laplacian of the electron density, as well as the self-consistently determined Kohn-Sham kinetic energy density. The nonlocal extension is based on some results of the local thermodynamic approach in DFT [ 10,111. Our starting point is the adiabatic connection formula [ 10,111 for the correlation energy, as defined in the Kohn-Sham scheme:

’ On leave from Bulgarian Academy of Sciences, Institute of Kinetics and Catalysis. * On Sabbatical leave from Department of Chemistry, Universidad Autonoma Metropolitana-Iztapalapa. 0009-2614/94/$07.00 0 1994 Elsevier Science B.V. All rights reserved SSDIOOO9-2614(94)01189-3

420

E.I. Proynov et al. / Chemical Physics Letters 230 (1994) 419-428

with the correlation-energy

s

x da.kk4R,

r)-11

>

density

(6)

(1)

0

where R= $ ( r1 + r2) (electron center-of-mass coordinate); r=r, -r2 (relativecoordinate); r= 1rl; ( )av denotes an averaging over the angular components of the relative coordinate r. The coupling-strength parameter 2 varying between 0 and 1 connects adiabatically the actual system of N interacting electrons (A= 1) with a fictitious (reference) system of noninteracting electrons (n=O). A key quantity in the above formula is the ;1- and spin-dependent electron pair correlation function (CF) g”,, . The following special class of Cf was considered [ 1 ] : g:,(R,r)=l-exp(-k$r2) x

[F,(~)+F*(~)l

>

F,(A)=m-@:,(R)(m+Ar),

(2) (3)

F,(A)=-exp(-/&r2)[1-@$,(m+1r) +(@?,)2(l+Ar+$12r2)],

(4)

where k, has the meaning of an ‘inverse correlation length’ with the dimensionality of a wave vector, and m is a parameter. For any real value of m the corresponding CF, Eqs. (2 )- (4)) automatically obeys the known exact cusp conditions [ 1 I, 121 for a l.-dependent antiparallel-spin CF (but not for the parallel-spin CF) in (R, Y) coordinates [ 11. Correlation between electrons of parallel spin (beyond the exchange-only level, n=O) is not taken into account in this model. The ;i dependence in Eqs. (2)- (4) is incorporated using the coordinate-scaling recipes of Levy [ 13 1. In the particular case of m = 2, ii = 1 and in the spin-unpolarized limit Eqs. (2)-(4) recover the ColleSalvetti pair CF [ 141. The function tiTJ (R) is determined from the sum rule for g:,, to zeroth order in the expansion of the quantity (n,(R+fr)n,(R-jr)),,aboutR [1,2].UsingEqs. (2)-(4) (with m=2) in formula (1) leads to the following result for the l-averaged (antiparallel-spin ) correlation energy [ 1,2] :

Q(k) = - &k

+ 2 in(v)

(7) where bj are numerical constants resulting from the integration over the relative coordinate in Eq. ( 1) [ 1,2]. For the particular case of m = 2 considered in Refs. [ 1,2] these are: b,=2.763169; b2= 1.757515; b3= 1.741397; b,=0.568985; bs= 1.572202; b6= 1.885389. Note that the way the functional (5)-( 7) is derived here does not depend on the particular form of ktL, which (apart from its dimensionality fixed by scaling arguments) still remains to be specified. The choice of k, is a key point of the present work, since via this quantity we incorporate nonlocal corrections. As in our previous works [ 1,2] (following Becke [ 15,16]), we approximate the anti-parallel spin correlation length as an average of the parallel spin ones,

,&-a2k,k, k, +k, ’

(8)

where cr is treated here as a fitting parameter. To zeroth order in the gradients of the electron density, the parallel spin k, is proportional to the corresponding local Fermi wavevector [ 17 ] : kFo= (67~~n,)“~,

(9)

and this was the way k, and k, were evaluated in Refs. [ 1,2 1. However, in inhomogeneous systems the correlation length should depend also on higher derivatives of the electron density, so as to ‘feel’ the boundaries of the system. This idea was successfully developed by Becke [ 15,16 1, especially in his recent model achieving results beyond the GGA [ 161. In this model the parallel-spin correlation lengths are determined from the Coulomb potential arising from the (nonlocal) exchange hole, while the antiparallel-spin correlation length is assumed to be proportional to

E.I. Prqynov et al. I Chemical

the sum of the parallel-spin ones (see Eq. ( 8) ). In the present work we consider a different from Ref. [ 16 ] nonlocal expression for the parallel-spin correlation lengths (k;’ in our notations), based on some results of the local thermodynamic approach in DFT [ 10,111. In the latter, the exchange hole (averaged over the angular components of the electron-distance coordinate r) is considered in the following approximate form: c&(R, r) = n,(R)

exp i --~‘/P~t~>

7:

(101

where P,(R) is the so-called ‘local temperature’ parameter, which actually plays the role of a Fermi-hole correlation length in the present context. Referring to this model enables us to approximate k2, with the inverse of /3, and use the expression for /$, derived in the local thermodynamic approach by means of a Wigner transformation of the first-order density matrix [ll]:

(11) where t,, is the noninteracting KS kinetic energy density (in the spin-polarized case). In this work t,, is evaluated self-consistently through the occupied KS orbitals {pig} [ 10,111:

1 MoPPic7~v$%o) t.w? = yj c I

IPi0 I’

(12)

Inserting Eqs. ( 11) and ( 12 ) in Eq. (7) we obtain a nonlocal generalization of the functional (6). The only derivative of the electron density entering explicitly in Eqs. (6), ( 11) and ( 12) is the Laplacian (inEq. (11)). All quantities entering Eqs, (S), ( 11) and ( 12) were readily available in the code deMon, so that almost no additional extension of the program was necessary for implementing the new nonlocal option. However, using the SCF procedure for the evaluation of the kinetic energy density makes the derivation of the corresponding nonlocal correlation potential rather obscure. In the present work we test the functional (Eqs. (5)-(7)) in a non-selfconsistent manner, using the corresponding gradientless limit of the correlation potential [ 17. We try to compensate partly for this inconsistency by using two different values of the correlation parameter LYin Eq. (8) - one for k,

Fhysics Letters 230 (I994j 419-428

421

entering the potential (a,) (in this case Eq. (9) is used instead of Eqs. ( 11) and ( 12) ), and another for the nonlocal expression of the correlation energy (a,) (Eqs. (5)- ( 12) ). Nonlocal corrections are often incorporated non-selfconsistently using the potential of the homogeneous-gas limit. From this perspective, the value of the parameter (Yeused in the (gradientless) correlation potential can be determined from the homogeneous-gas test explored previously [ 11. We found that a small reduction of (Y, from its homogeneous gas value is necessary to yield optimal bond distances. Variations of cy, appear to affect primarily the geometry. On the other hand, the parameter ay, entering the Laplacian correction, Eqs. (8)-( 12 ), mainly influences the energy values. The problem here is how to achieve the best fit for cr,? As it was shown in Refs. [ 1,2], such a tit must be synchronized with, and depends strongly on, the chosen exchange counterpart. Here we use the nonlocal exchange functionals of Becke [ 5 ] and Perdew [ 6,7] as partners of Eqs. ( 5 )- ( 7 ). The usual way of proceeding then would be to adjust o!, in Eqs. f 8) and ( 11) so as to reproduce as much as possible the experimentally derived correlation energies of a set of atoms and small moiecules [ 1,2,6-S 1. It was found, however [ 21, that such a procedure does not necessarily lead to optimal binding energies - one of the main observables evaluated in DFT. Moreover, one has to keep in mind that the correlation energy as defined and evaluated in Kohn-Sham DFT may in general be different from the conventional one, based on Hartree-Fock theory [ 18-223. In Table 1 we present further evidence that, at Ieast for some molecules, the above mentioned difference may be quite substantial. The first two rows of Table 1 contain the experimentally derived binding-energy contribution from the correlation alone (A&), as well as the remainder (exchange, kinetic and electronnuclear contributions denoted as AE,,) for several well known small molecules. The corresponding experimental binding energies are given in the third row. We have calculated all these quantities with two different nonlocal (GGA) functionals - PW9 1 and BP, which have been widely applied in the deMon code, and the results (given as deviations from the experiment) appear in the body of Table 1. (For the sake of comparabihty we have used the same Gaussian basis set (orbital and auxiliary), as in our previous

E. I. Proptov ef al. / Ch~~icu~ Physics Letters 230 (I 994) 4 19-428

422 Table 1 Correlation

contribution

bE,(exp.) A.%, (exp.

&Jew) PW91

BP

BLap 1

BLap2

PLapl

1

(talc.-exp.)

Hz0

to the binding

N2

energy (eV) for several molecules

h

-

3.1 7.07

4.73 5.18

2.53 -0.87

10.17

9.91

1.66

CN

PN

CK

C&

4.4 1 3.48

3.56 2.88

3.7 14.7

4.63 20.02

7.89

6.44

18.4

24.65

-0.49 +0.70

-2.17 i-2.54

-1.69 +2.25

-2.39 +3.12

-2.45 + 2.29

- 1.36 -t 1.43

+0.46 -0.45

+0.86 -1.19

+0.21

+0.37

i-O.56

f0.73

-0.16

+0.07

+o.oi

-0.33

-0.30 to.36

-2.09 +2.26

- 1.60 +2.07

-2.26 +2.83

-2.29 i-2.53

- 1.30 +I.24

+0.49 -0.20

+ 1.25 -0.88

+ 0.06

to.17

+ 0.47

+0.57

to.24

-0.06

+0.29

+0.37

-0.65 +0.64

-2.13 +2.31

-1.71 +2.10

-2.53 +2.84

-2.51 +2.49

-1.34 + 1.24

+0.23 -0.26

t 0.49 -0.63

-0.01

+0.18

+0.39

+0.31

-0.08

-0.10

-0.03

-0.14

-0.58 +0.59

- 2.08 •t 2.27

- 1.64 + 1.97

-2.51 +2.71

- 2.52 + 2.46

-1.27 + 1.20

+0.2s -0.20

+0.57 -0.60

+0.01

io.19

+0.33

+0.20

-0.06

-0.07

to.08

-0.03

- 0.66 f0.63

-2.14 t 2.63

-1.72 + 2.09

-2.67 +2.96

-2.58 i-2.71

- 1.35 +1.40

+0.23 -0.28

$0.49 -0.60

-0.03

+0.49

+0.37

+0.29

+0.13

+0.05

-0.05

-0.11

works [ 1,2].) It is seen that A.& obtained with the GGA is smaller than the experimental one, except for the two hydrocarbon molecules involved in the test, CH, and C2H4. In many cases the underestimation is quite pronounced - between 2 and 3 eV. However, the rest of the binding energy (A&,) (for which we use the label ‘exchange’ contribution) always turns out to deviate from its experimental counterpart in the opposite direction (compared to A&) so that the two ‘errors’ cancel each other to a large extent. As a result, the overall binding energies calculated with PW9 1 and BP are pretty good. As it is seen from Table 1) the above-mentioned ‘cancelation of errors’ is too systematic to be fortuitous. Since at the same time the correlation functionals within the PW9 1 and BP schemes give atomic correlation energies very close to the experimentally derived (conventionally defined) ones (see Tabie 2), it is clear from Table 1 that at least for the molecules shown the situation is just the opposite - the KS-DFT based correlation en-

ergy of molecules may be quite different from the conventionally defined one. This problem has been intensively discussed in the literature (see Refs. [ 2 1,221 and the references therein). A comparative analysis of the difference between the conventional correlation energy and the one evaluated in DFT was made recently by Pople et al. [22] on a large set of molecules, using the Becke (86) exchange plus LeeYang-Parr correlation (BLYP). This paper underlined some of the difficulties associated with the partitioning of &, into parallel and antiparallel spin parts in DFT. Keeping the above discussed problems in mind, we decided to use binding-energy data like those in Table 1, involving exchange and correlation together, rather than the experimental correlation energies afone, in order to fix properly the parameter LY,.Taking the exchange functionals of Becke [ 5 ] and Perdew 173, we tried to adjust the correlation term (Eqs. (S)-(7), (8)-{ll)),soastoachievethebest possible blinding energies, i.e. to compensate as much

423

E.I. Proynov et al. / Chemical Physics Letters 230 (1994) 419-428 Table 2 Correlation -

energy differences

(ET

H

--4X,” PW9 1 BP BLap 1 BLap2 PLap 1 PLap2 ’ The experimental

-E:“* ) (au) for some atoms N

0.0000 -0.0010 -0.0025 0.000 0.000 0.000 0.000 correlation

0

F 0.189

0.322

0.258

+ 0.009 +0.015 + 0.002 +0.0008 + 0.002 +0.001

- 0.004 +0.006 +0.007 +o.oos + 0.007 +0.001

-0.001 +0.007 - 0.003 -0.005 -0.003 -0.004

energies for atoms are taken from Ref.

as possible for the missing part of that energy beyond the exchange-only contribution. In doing so we found that the magnitude of AEc is very sensitive to variations of the coefficient b, multiplying the logarithmic term of the rhs of Eq. (7 ). Thus, at the expense of one extra parameter (the coefficient b3 in Eq. ( 7 ) ) , the present correlation functional can be adjusted to any exchange counterpart separately, so as to achieve optimal binding energies. Let us recall that the original values of all bi coefficients in Eq. (7) were obtained as a result of the analytical integration over the relative coordinate in Eq. ( 1) using Eqs. (2 )- (4 ). A more consistent way of governing the magnitude of b3 would be to vary the parameter m in Eqs. (2)(4) i.e. to pick up other members of this family of CF. This, however, is a rather tedious procedure and will not be considered in this work. Using a trial-anderror approach we found values of the parameters a’, and b3 for which the combination of Becke exchange [ 5 ] and Eqs. (5)-( 11) (this scheme is denoted here as BLap 1) gives improved results for the binding energies of a set of 19 test molecules. In the BLap 1 scheme the exchange is handled self-consistently using the corresponding nonlocal potential. As was mentioned above, the Laplacian corrections are included (at each SCF cycle) only in the correlation energy while the correlation potential is obtained in the gradientless limit. The same set of parameters turned out to be suitable when Eqs. (5)-( 7) are combined with the exchange functional of Perdes I:86) [ 71 (PLapl scheme). Two other schemes are also tested here (BLap2 and PLap:!), in which all nonlocal corrections (including the exchange ones) are taken into account only in the energy and in a perturbative way (after the last SCF cycle). In these

Be

C 0.157 +0.007 + 0.007 0.000 - 0.0005 0.000 - 0.0005

0.0944 +0.0006 -0.0001 +0.003 + 0.002 +0.003 +0.003

[ 91. schemes the self-consistency is first reached using the local version of the correlation functional, Eqs. (5)(7), i.e. using Eqs. (8) and (9) for k, [ 11, combined with the standard LSD exchange. To summarize, the optimal values of the parameters used in the new schemes are found to be: in BLapl and PLapl: (Y,= 1.29 (in the Laplacian correction); (~~~0.2 (in the correlation potential); b3= 1.45: in BLap2 and PLap2: ay,= 1.27; cr,=O. 196; b3= 1.48. Note that the optimal value of b3 is slightly reduced compared to its theoretical value b, = 1.74 1. The value of (x, used in the correlation potential is also slightly reduced with respect to the value found from the homogeneous-gas test [ 1 ] (a = 0.2026). As mentioned above, the re-scaling of (Y, in the correlation potential allows us to achieve improved molecular geometries. Strictly speaking, (Y, and (Y, represent one and the same parameter ((x in Eq. (8 ) ) used in different contexts. The use of these two parameters was necessary to partly compensate for the missing (unknown) precise SCF potential given by the functional derivative of the Laplacian corrected functional with respect to the spin densities. The large difference between the optimal values of cy, and o!, indicates on one hand that the SCF potential of the Laplacian-corrected functional, Eqs. ( 5 )- ( 7 ), ( 8 ) and ( 11) should be rather different from its gradientless counterpart corresponding to Eqs. (5)-( 9) (as derived in Ref. [ 1 ] ). On the other hand, we found that the parametrical dependence of the correlation energy, Eq. ( 5 ), on a! (a,) changes quantitatively when one uses in Eq. (8) the nonlocal expression ( 11) for the parallel-spin wavevectors k, instead of the local

E.I. Pro.vnovet al. / ChemrcalPhysics Letters 230 (1994) 419-428

424

one, proportional to the local Fermi wavevector, Eq. (9). The real space distribution of the nonlocal wavevector k, given by Eqs. (8) and ( I 1f is generally quite different from its local counterpart, Eqs. (8) and (9). While the rhs of Eq. (9) is a positively defined monotonic function of the electron (spin) density only, the rhs of Eq. ( 11) encounters the changes of the curvature of the electron density distribution, which is one of the important aspects of the inhomogeneity. However, at certain points of space Eq. ( 11) may lead to negative values for ki and these nonphysical values have to be avoided. At such points we use the gradientless limit of the rhs of Eq. ( 1 1 ), in which kc is propo~ional to kFOwith a constant of prop~)~ionality 1,/a [ IO,1 11. This problem resembles the ‘turning point’ problem encountered in the Kirzhnits expansion method of deriving gradient corrections in DFT [ 11, p. 85 1. The new options were tested on a set of atoms and small molecules, and the results are shown in Tables 2-6. It is worth noticing, that while the parameters of the correlation functional were optimized with respect to the binding energies of the test molecules, the so-adjusted scheme

Table 3 Binding energy differences

(Ep

-ET@)

NZ

2

$”

a

PW9 BP PP

1

BLapl BLap2 PLapl PLap2

9.9

CH

1.66

CK 18.40

3.65

NH3

GH‘l 24.65

PH3

10.77

13.05

PW91 BP PP

-0.39 -to.15 +0.20

+0.01 +0.29 +0.27

-0.33 +0.37 +0.41

+0..57 $0.38 +0.46

+O.Oi -0.09 -0.13

BLapl BLap2 PLapl PLap2

+0.14 +0.17 +0.19 +0.23

-0.03 +0.08 -0.05 +0.08

-0.14 - 0.03 -0.11 +0.03

+0.18 +0.28 +0.26 +0.37

-0.19 -0.12 -0.24 -0.15

gives correlation energies for atoms in very good agreement with the results of PW9 1 and PP, and with the experimental values (Table 2). This is in line with earlier observations that for light atoms (in contrast to the molecules shown on Table 1) the KS-DFT based correlation energy should be close to the conventional one [ 20 1. Turning to the binding energies (Tables 3 and 4), we note that the absolute error for BLapl and BLap2

feV) for several molecules

F2 1

Table 4 As in Table 3 for some other molecules

(the zero-point

energy is added) HF

02 5.21

4.41

5.08

2.51

6.12

+0.37 +0.17 +0.49

+0.56 i-o.47 +0.45

+0.73 +0.57 +0.56

+0.20 +o.os -0.03

-0.03 -0.17 -0.14

-0.04 -0.16 -0.23

+0.10 -0.01 - 0.49

+0.17 +0.19 + 0.49 +0.52

+0.39 +0.33 +0.37 + 0.28

+0.31 +0.20 +0.29 f0.28

-0.14 -0.21 -0.21 -0.21

-0.29 -0.21 -0.29 -0.25

-0.26 - 0.29 -0.32 -0.36

-0.001 i-o.04 -0.09 -0.04

50

SiO

NO

CN

PN

5s

NC1

E‘$”

10.17

PW9l BP PP

+0.21 +0.06 - 0.005

-0.14 -0.26 -0.26

+0.62 +0.43 +0.59

-0.16 +0.24 +0.46

+0.07 -0.06 -to.10

+0.17 +0.02 - 0.06

+0.12 +0.04 -0.02

BLapl BLap2 PLap L PLap2

-0.01 +0.01 -0.03 +0.0005

-0.23 -0.28 -0.22 -0.21

+0.27 +0.10 +0.43 +0.31

-0.08 -0.07 i-o.13 +0.16

-0.10 -0.08 -to.05 +o.os

-0.06 -0.03 -0.14 -0.07

-0.02 +0.02 -0.09 -0.04

a The experimental

binding energies

8.34

(corrected

6.61

for the zero-point

7.89

energy)

6.44

are taken from Refs.

[ 21.231

8.02

4.62

425

E.I. Proynov et al. / Chemical Physics Letters 230 (I 994) 419-428

Table 5 Ground state bond distances

(au) (talc. -exp.)

N2 R

a

“P

F2

2.0142

I

PW9

for several molecules O* 2.2818

2.6682

obtained

with different

functionals

NH,

Hz0

HS

CzH4

NH 1.9123

OH 1.8088

SH 2.5099

CH: 2.0333 (CC: 2.5132) + 0.007 ( + 0.005 ) +0.034 (+0.025) +0.035 (+0.027)

+0.037

+0.011

+ 0.066

+0.016

+0.017

+0.027

BP

+ 0.037

+ 0.038

+0.059

+0.021

+0.030

+0.052

PP

+0.042

+ 0.060

+0.071

+0.02s

+0.032

+0.068

BLapl

+0.022

+ 0.047

+0.055

-0.004

+0.010

+ 0.022

BLap2

+0.011

-0.055

+ 0.00s

+0.0002

+0.013

+ 0.03

+ 0.024

+0.08

+ 0.065

-0.003

+0.015

+0.026

+0.011

-0.055

+ 0.008

+ 0.0002

+0.013

+ 0.02

HF

CH

NO

CN

PN

CH,

PH,

CH 2.0502

PH 2.6834

-0.001 + 0.026 +0.026

-0.001 + 0.03 +0.02s

- 0.0004 + 0.007 +0.0001 + 0.007

-0.008 +0.002 - 0.003 + 0.002

PLap

I

PLap2

R CXP

1.7328

PW91 BP

2.1 162

+0.019 +0.03s

2.1746

+0.049 +0.043

2.2143

+0.045 + 0.046

+ 0.006 + 0.020

2.8173

PP

+ 0.042

+0.039

+0.051

+0.022

+0.032 +0.035 +0.055

BLap 1 BLap2 PLapl PLap2

+0.021 + 0.023 + 0.027 +0.023

+0.006 +0.025 + 0.004 +0.025

+ 0.030 +0.006 +0.037 + 0.006

+0.003 -0.008 +0.006 -0.008

+0.014 - 0.008 +0.021 - 0.00s

a The experimental

geometries

are taken from Refs.

b

+ 0.006 (+0.004) +0.019 (-0.013) + 0.007 (+o.oos) +0.01s (-0.014)

[ 23,241.

Table 6 As in Table 5 for some other molecules

sz R exp

a

PW9

I

P2

3.5699

Cl2

3.5779

3.1547

SiO 2.8528

NaCl 4.4611

NaF 3.6393

+0.095

+ 0.039

+0.107

+0.044

-0.109

BP PP

+0.10s +0.15

+0.057 +0.057

+0.125 +0.126

+0.059 +0.077

+ 0.065 + 0.063

-0.057 +0.065 +0.09

BLapl BLap2 PLap 1 PLap2

+0.102 +0.035 +0.119 +0.036

+ + -

+0.141 +0.019 + 0.087 +0.011

+0.042 +0.013 +0.049 +0.013

+ 0.067 -0.061 + 0.062 - 0.065

+0.073 -0.034 +0.070 -0.032

0.033 0.004 0.047 0.004

does not exceed 0.4 eV anywhere. The same holds also for PLapl and PLap2 with the exception of Nz, and NO, where the absolute error is slightly above 0.4 eV.

Within the set of homonuclear diatomics (Table 3 ), the Laplacian-corrected schemes give improved resuits for the ‘difficult’ molecules F2 and 02, while for

426

E.I. Proynov et al. /Chemical

N2 the results of BLap and PLap are almost the same as those of BP and PPK, respectively. The binding energy of the second row diatomics SZ, P2 and CIZ is somewhat underestimated with the new schemes, the errors being relatively larger than the corresponding GGA ones. For the rest of the molecules considered, the Laplacian-corrected functionals give systematically better binding energies than the GGA schemes, especially for NO (BLap), HF (PLap), NH,, CH, and CzH4. An exception is PH3, where the SCF (with respect to the exchange) options BLapl and PLapl give slightly worsened results, in contrast to their perturbative counterparts BLap2 and PLap2. Concerning the bond distances (Tables 5 and 6), the new options achieve small but nearly systematic improvement. Among the homonuclear diatomics, F2 is the only case where both SCF schemes BLap 1 and PLap I do not improve the bond distance compared to GGA. BLapl also gives slightly worsened bond distance for CIZ. The success of the perturbative schemes BLap2 and PLap2 is here even more pronounced, giving the best bond distances in all cases, except for F2, and C2H4. In the latter case the SCF options BLapl and PLapl give the best result, comparable to PW9 1. The absolute errors in the bond lengths calculated with the Laplacian corrected options show an overall decrease when moving from homonuclear diatomics to polyatomic molecules. It is worth mentioning at this point that while the GGA leads to improved energetics for many molecules (compared to LSD), the GGA exchange and correlation potentials appear to contain some inaccuracies [25] that can be ‘turned off’ by using the perturbative approach. In Figs. 1 and 2 we present a percentile statistical analysis [26] of the results in terms of error distributions. Each box shown in the figures contains 90% of the data. The solid line inside each box corresponds to the median, i.e. the value corresponding to 50% of the data sorted in increasing value of the error. The lower and upper dashed lines represent 25% and 75% of the data, respectively. The distance between these two dashed lines is called the interquartile distance (IQD). The smaller the IQD, the narrower the distribution of the data around the median. One gets better agreement with experiment when the median is closer to the zero line. For the binding energies (Fig. 1)) all Laplacian-corrected schemes

Physics Leilers 230 (1994) 419-428

II,88, \//II,, / INI% I, /,j II, /, %/I /)I,

Binding Energies Calculated - Experimental LO I ~~~ I

Fig. 1. Percentile analysis of the binding energy errors. erage absolute error for each of the functionals considered in the bottom of the figure.

Bond Calculated 0.1

I

Distances - Experimental _~

Fig. 2. The same as Fig. 1 for the equilibrium

The avis given

*II1IYi(l(.lltl1 1 ~~~~ N I tI \ cI -,

I

/// \,,a uo
geometry

errors.

have a smaller IQD than the corresponding GGA functionals. The medians of BLap2 and PLap2 data are closer to the zero line than those of BLapland PLapl, and all the new options are better in this respect than the GGA options. For the bond distances (Fig. 2), the perturbative schemes BLap2 and PLap2 have the smallest IQD and a median very close to zero. The SCF schemes BLap 1 and PLap 1 have slightly larger IQD than BP and PP, respectively, but their medians are closer to the zero line, similar to PW91. To identify the molecules whose calculated binding energy and bond distances deviate considerably from the experimental values we plotted in Fig. 2 those molecules that are outside the

EL Proynov et al. / Chemical Physics Letters 230 (1994) 419-428

box containing 90% of the data (outliers) for at least one functional. For the binding energies there are no outliers (Fig. I), the worst cases lying on the borders of the boxes. In contrast, for the bond distances we observe outliers, Fig. 2. These are among the halogen containing molecules and/or Sz. since we have used standard basis sets in these calculations, we cannot at this stage rule out the possibility that some of the outliers are associated with basis set problems. Comparing BLap 1 ( BLap2 ) with PLap 1 (PLap2 ) options, in general we find that the nonlocal exchange of Becke is, overall, a better partner to the present correlation functional, Eqs. (5)-(7), (8) and ( 11)) than the Perdew (86 ) exchange. The above analysis shows that the Laplacian-corrected options improve both the evaluated geometry and binding energy of the majority of the molecules considered, compared to GGA functionals. The perturbative schemes BLap2 and PLap2 give binding energies of a comparable quality to BLap 1 and PLapl , while being superior with respect to the geometry in most of the cases. Important exceptions are methane and ethylene, where we find the SCF schemes PLap 1 and BLapl to be particularly precise. Extended tests in that direction are currently in progress. Concerning the exchange alone, the comparison between BLap 1 ( PLap 1) and BLap2 (PLap2) shows that the CGA nonlocal corrections in the SCF exchange potential has primarily impact on the optimal geometries, while the energetics depends mainly on the nonlocal corrections in the functional itself. Note in this respect that the perturbative schemes BLap2 and PLap2 give slightly better binding energies for the ‘difficult’ molecules F2, O2 (Table 3 ) and C2H, (Table 4) than the corresponding SCF schemes BLap 1 and PLap 1. To conclude, new nonlocal exchange-correlation schemes have been derived and implemented in the LCGTO DFT code deMon, based on the Laplaciancorrected correlation functional, Eqs. ( 5 )- ( 7 ), and the GGA exchange functionals of Becke (88) and Perdew (86). The new options give improved results for the binding energy and equilibrium geometry within the set of test molecules considered. The perturbative options BLap2 and PLap2 seem to be especially promising: their accuracy in the binding energies is comparable to that of the SCF options while giving superior geometries, with less computer

427

time. Tests on further systems (transition metals, hydrogen bonds, intermolecular interactions, etc.) will be necessary to establish the realm of applicability of the proposed functionals. Such tests are planned for the near future. We thank Dr. A. Becke for the advance preprint of Ref. [ 161. We also gratefully acknowledge financial support from NSERC, and computing resources from the Services Informatiques de 1’UniversitC de Montreal. We also acknowledge DGSCA/UNAM - Mexico for providing us with computer time on the Cray Y-MP41464.

References [ 1 ] E.I. Proynov and D.R. Salahub, Phys. Rev. B 49 ( 1994) 7874. [2] E.I. Proynov, A. Vela and D.R. Salahub, submitted for publication. [3] A. St-Amant and D.R. Salahub, Chem. Phys. Letters 169 (1990) 3x7; A. St-Amant, Ph.D. Thesis, Universitt de Montreal (1992). [4] S.H. Vosko, L. Wilk and M. Nusair, Can. J. Phys. 5X (1980) 1200. [5] A.D. Becke, Phys. Rev. A 38 (1988) 3098. [6] J.P. Perdew, Phys. Rev. B 33 (1986) X822; 34 (1986) 7406. [7] J.P. Perdew and Y. Wang. Phys. Rev. B 33 (1986) 8800. [ 8 ] J.P. Perdew, Physica B 172 ( 199 1) 1; J.P. Perdew and Y. Wang, unpublished. [9] 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. lo] R.G. Parr and W. Yang, Density-functional theory of atoms and molecules (Oxford Univ. Press. Oxford, 1989). 11 ] R.M. Dreizler and E.K.U. Gross, Density functional theory. An approach to the quantum many-body problem (Springer. Berlin, 1990). 121 A. Rajagopal. J. Kimball and M. Banerjee, Phys. Rev. B 1X (197X) 2339. [ 131 M. Levy, Phys. Rev, .A 43 ( 1991) 4637. [ 141 R. Colle and 0. Salvetti, Theoret. Chim. Acta 37 (1975) 329. [ 151 A. Becke, J. Chem. Phys. XX (1988) 1053. [ 161 A. Becke, Intern. J. Quantum Chem. Symp., in press. [ 171 E.I. Proynov and D.R. Salahub, Intern. Quantum Chem. 49 (1994) 67. [1X] M. Levy and J.P. Perdew, Intern. J. Quantum Chem. 49 ( 1994) 539. [ 19 1M. Levy. Advan. Quantum Chem. 2 I ( 1990) 7 1. [20] L.C. Wilson and M. Levy, Phys. Rev. B 41 (1990) 12930. [ 2 11 E. Clementi and S.J. Chakravorty, J. Chem. Phys. 93 ( 1990) 2591.

428

E.I. Proynov et al. / Chemml

f21] P.hf.W. Gill, B.G. folmson and J-A. Popie, Intern. J.

Quantum Chem. Synrp. 25 ( I992f 319. f23] KP. Hnber and G. Her&erg, Molecular spectra and molecuIar structure. [23] KM. Hellwege and A.M. Hellwege, eds., La~do~t-Holstein, Numerical data and functional relationships in science and technology, Group II, Vol. 7. Structure data of free polyatomic molecules (Springer, Berlin, 1979).

Physics Letters 230 (1994) 419-428

[ZS] C.3. Umrigar and X. Gonze, in: Proceedings of the Conference an Concurrent fomenting in the Physical Sciences,inpress.

[as] L.L. Haviicek and R.D. &in. Practical statistics for the physical sciences, ACS professionat reference book (Am. &hem. Sot., Washington, 1958 ).