Chemical Physics Letters

The exchange energy of a uniform electron gas experiencing a new, flexible range separation

John A. Parkhill

^a,b

, Jeng-Da Chai

^a,b,c

, Anthony D. Dutoi

^d

, Martin Head-Gordon

^a,b,*

aDepartment of Chemistry, University of California, Berkeley, CA 94720, United States

bChemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, United States

cDepartment of Physics, National Taiwan University, Taipei 10617, Taiwan

dInstitute of Physical and Theoretical Chemistry, Ruprecht Karls-University Heidelberg, Im Neuenheimer Feld, 229, 69120, Germany

a r t i c l e i n f o

Article history:

Received 12 April 2009 In final form 16 July 2009 Available online 19 July 2009

a b s t r a c t

The exchange energy of a uniform electron gas which experiences a two-parameter separation of the Coulomb interaction is derived as a local functional of the electron density. The two parameter range sep- arator allows separate control of where and how rapidly the Coulomb interaction is switched off. The use- fulness of the functional is briefly assessed by combination with a recently published pair of exchange and correlation functionals. The self-interaction error of noble-gas dimer cation dissociation is found to be reduced while thermochemistry is relatively unperturbed. These results suggest that changes in attenuator shape can improve range-separated functionals.

Ó 2009 Elsevier B.V. All rights reserved.

1. Introduction

Despite rough beginnings[1], the local density approximation (LDA) has been developed through decades of work on the Kohn–

Sham (KS)[2]construction into one of the most successful approx- imations in quantum chemistry and solid state physics. Within this framework, the LDA exchange correlation functional is combined by an adiabatic connection with a non-interacting wavefunction so that an approximate kinetic energy may be extracted and there is no need to develop accurate functionals for the kinetic energy[3]

which have proven elusive. Along these lines, Becke[4]realized that the accuracy of Kohn–Sham energy functionals could be im- proved by the admixture of ‘exact’ exchange coming from the ex- plicit exchange energy of the fictitious Kohn–Sham wavefunction.

The resulting hybrid density functionals have been the most com- monly applied model chemistry for many years[5]because they have been found to be remarkably accurate with computational costs virtually equivalent to those of the Hartree–Fock (HF) method.

One of the few remaining substantial defects of the Kohn–Sham construction which has attracted theoretical effort is the so-called self-interaction problem[6–12], and it is directly related to the treatment of the exchange energy[13]. In the HF energy expression the Coulomb repulsion of a one-electron function with itself is can- celled exactly by the corresponding exchange integral. In the KS

construction with a pure, local functional the Coulomb energy is non-local, but the exchange energy is not. Considering the one-par- ticle functions provided by the KS wavefunction we might say that the electron repels itself if the particle is spread over space because the antisymmetric complement of the Coulomb interaction, non- local exchange, is missing. At equilibrium geometries the effect on predicted ground state energies is not severe, but this defect means that dissociation problems may lead to fragments which only possess a fractional number of electrons[14,15], or response properties which reflect serious artifacts if charge is significantly redistributed[16]. If globally a fraction of the exchange energy of the KS determinant is mixed with the DFT exchange energy these artifacts are partially remediated. To a stranger unfamiliar with the history of hybrid DFT’s development the situation must seem confusing, because it is not obvious why any mixture of ‘exact’ ex- change with (semi-)local exchange is advantageous. The answer is that the accuracy of most DFT functionals lies in a cancellation of errors between exchange and correlation functionals. Both are compensating for the single-reference nature of the fictitious KS determinant[17,18]. From another angle, one might say that these two non-local objects[19] are best considered together because what results is more local.

A way to preserve the local cancellation of errors yet recover correct exchange at long range has emerged in the form of range-separated hybrid functionals. The idea which goes back to the pioneering work of Gill and Savin[20–22], is to divide 1=r by multiplying it with a function which varies between 0 and 1, such that both this function and its complement are integrable. The greatest fraction of work in this very active area[23–28]has em- ployed the standard error function to achieve this separation:

0009-2614/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved.

doi:10.1016/j.cplett.2009.07.052

* Corresponding author. Address: Department of Chemistry, University of Cali- fornia, Berkeley, CA 94720, United States. Fax: +1 510 643 1255.

E-mail addresses: [email protected] (J.A. Parkhill), [email protected].

edu.tw(J.-D. Chai),[email protected](A.D. Dutoi),mhg@cchem.

berkeley.edu,[email protected](M. Head-Gordon).

Chemical Physics Letters 478 (2009) 283–286

Contents lists available atScienceDirect

Chemical Physics Letters

j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c p l e t t

1=r ¼ erfð

x

rÞ=r þ erfcð

x

rÞ=r: ð1Þ The LDA exchange functional corresponding to erfcð

x

rÞ=r, and inte- gral kernel for the exact exchange over erfð

x

rÞ=r can then be de- rived so that locally exchange is provided by the LDA and at a distance exchange is provided by the KS wavefunction. The position where the transition is smoothly made between the two treatments is determined by the adjustable parameter 1=

x

. The choice of the error function as a Coulomb attenuator is both practical (for most implementations one must be able to perform the integral of Gaus- sians over the function analytically[29]) and arbitrary because the error function is just one of many which possess this property. In some recent studies promising results have been attributed to more flexible range separation[30].

One can extend the idea of mixing ab initio and DFT strengths further by imagining range separation of the correlation part of the functional as well. In this scheme (which we will not pursue in this work beyond mention) the ab initio method is made respon- sible for static and long range dispersion effects while the LDA cor- relation functional is adjusted for the modified Coulomb interaction to avoid double counting. Savin and coworkers have experimented with the choice of another attenuated Coulomb interaction for these purposes[31], a linear combination of an error function and a Gaussian which offers a sharper separation:

v

^{ee;erfgau¼ erfð}

^x

^{rÞ=r  ð2}

^l

^{Þ=pffiffiffiffi}

^p

^{ eð1=3Þl2r2: ð2Þ}

More recently this erf–gau LDA functional was combined with a standard GGA in an attempt to surpass the accuracy of exchange- hybrid functionals based on erf[32]with a GGA correlation treat- ment. Another one-parameter attenuated Coulomb interaction, the Yukawa potential, has also been the subject of recent investiga- tions[33,34]

v

ee;yukawa¼ expð

c

rÞ=r: ð3Þ

Improved performance of the resulting functional was attributed to an increased fraction of short-range exchange[35].

Our group has recently published an analytical integral over a more general sort of Coulomb attenuator which allows for separate control of where and how rapidly the shift is made between parts of the Coulomb interaction[36]. It allows continuous variation of sharpness between the limits of erf and the Heavyside function.

The function is a linear combination of two error functions (although note that the erf integral formulas are not sufficient to describe it) and so we have adopted the name ‘terf’:

terfr0;xðrÞ ¼ ð1=2Þðerfð

x

r 

x

r0Þ þ erfð

x

r þ

x

r0ÞÞ: ð4Þ Investigations into the performance of range-separated hybrids[37]

have found that existing separations cannot simultaneously de- scribe ground state and excited state properties with a single choice of

x

. Along these lines a common area of intersection has been lo- cated amongst many optimized attenuators at roughly 0.8 Bohr[38]

in the ðr; VoptimalðrÞÞ plane. The physical implication is that cancella- tion of GGA-exchange and GGA-correlation errors in this region is balanced with the error induced by semi-local exchange, and our attenuator should be shaped similarly in this region for thermo- chemical accuracy. Yet to repair self-interaction error the attenua- tor should reach its asymptote as rapidly as possible once we leave this region. The terf functional form can pass through this point while still reaching its asymptote more rapidly than erf, and so there is reason to hope that terf may be useful in this respect. An- other nice feature of this choice of separation is that it reduces to the erf attenuator if the r0parameter is chosen to be zero. The atten- uator controlling the fraction of exact exchange to include as a func- tion of electronic separation is plotted for various choices of parameters inFig. 1: it consists of terf itself as defined above, plus

a fraction of short-range exact exchange (i.e. terfc) controlled by a third parameter, cx:

LðrÞ ¼ cxterfcr₀;xðrÞ þ terfr₀;xðrÞ: ð5Þ The specific sets of parameters plotted inFig. 1are choices whose performance for chemical problems will be explored later in this paper. They all include roughly similar fractions of exact exchange in the mid-range region, as is believed to be important[38].

2. The terf-attenuated LDA

The exchange energy of a many-fermion system, charge bal- anced by a structureless positive background is our starting point.

This matrix element is [39] (where h denotes the Heaviside function)

Ex¼ k³_f 12

p

⁴

Z 1 0

q²

v

^{eeðqÞ 1 3}₂^{x þ1}₂^x3

 

hð1  xÞdq;

where kf ¼ ð3

p

²nÞ¹⁼³;x ¼ q=ð2kfÞ: ð6Þ So we must obtain the Fourier transform of terfcr0;xðrÞ=r.

Ffterfcr0;xðrÞ=rg ¼

v

eeðqÞ ¼

4

p

1  e^4q2x2cos^qx_x^r0

q² : ð7Þ

The integration of this function is algebraically quite tedious, but can be done. Unfortunately the complex error function enters. Note that for z 2 C; erfðz^Þ ¼ erfðzÞ^and for z 2 R; erfiðzÞ 2 R. We report the exchange energy per particle of the Fermion gas experiencing this interaction,



xcwhich may be readily implemented in any KS- DFT code.

Ex¼ Z

nðRÞ



xcðnðRÞÞdR: ð8Þ

The spinless k_f can be easily replaced to obtain the spin-density functional



^x_x^;r0ðnÞ ¼ x⁴

192A⁴p^{3 8Ae}

x²r²

0^pffiffiffiffip Axr0

 3 þ A²



4x²r²₀ 6



erfi½xr0







þIm erf 1 2A ixr0

 

 

þ Re erf 1 2A ixr0

 

 

þ



ð3 þ 24A²þ 32A⁴ðx²r²₀ 1Þ



þe^4A2116A² ð1 þ 2A²ðx²r²₀ 1ÞÞ cos xr0

A



þ Axr0sin xr0

A



 o

; where A ¼ x

2kf

: ð9Þ

One may easily verify that this expression matches the known erf expressions for



x as r0! 0 [20,40]. For large values of A (>1), a Fig. 1. Various attenuators that control the fraction of exact exchange as a function of inter-electronic distance, r, plotted for comparison. The first two are equivalent to erf ðx¼ 0:3; 0:4Þ.

284 J.A. Parkhill et al. / Chemical Physics Letters 478 (2009) 283–286

series expansion in powers of 1=A is employed up to 10th order in our implementation for purposes of numerical stability.

3. Application to range-separated hybrids

Without semi-local gradient information the thermochemistry of this functional would be undoubtedly poor and it would be dif- ficult to determine if terf could improve functionals in use today.

There are several recipes for combining this LDA exchange func- tional with a GGA enhancement factor ranging in degrees of tech- nical difficulty and empiricism. Ideally the GGA factor will depend on the attenuating parameters[41], but recent results have shown

that superior accuracy [13]can be obtained even if this is only done implicitly through optimized parameters of the GGA. At the end of the day the choice of

x

parameter is quite empirical, as will be r0, even if we introduce them for physical reasons. The final measure of a range-separated hybrid is optimization over a large training set, and evaluation over an independent test set, roughly a year of computer effort. We seek some justification for such an effort and so we combine the terfc-LDA exchange energy with the GGA exchange enhancement factor and correlation functionals of

x

B97X[13]and run some basic tests to establish whether the resulting functional shows promise. To be precise, the resulting functional is obtained directly from xB97X by replacing the FðarÞ of Eq. (7) in that paper with the corresponding terfc-LDA Fða_rÞ obtained from Eq.(8). To begin from a functional as close as possible to the parent (see the previous paper to clarify the nota- tion), we also incorporate a variable fraction of short-range HF ex- change in such a way that the UEG limit is respected (Eqs.(9) and (10))

E^SR—DFA_x ¼X r

Z

e^terfc—LSDA_xr ð

q

_rÞg^xxr^B97Xðs²_rÞdr; ð10Þ Exc¼ E^LR—HFx þ cxE^SR—HF_x þ ð1  cxÞE^SR—DFAx þ E^xc^B97X: ð11Þ Aside from the many parameters associated with the GGA we must choose reasonable guesses of fr₀;

x

;c_xg. The physically moti- vated guess is to reach the asymptote as rapidly as possible while still overlapping significantly with the established attenuators in the previously mentioned critical region[38]. An initial choice of parameters fr0;

x

;cxg ¼ f1:2; 1:016; c^x_x^B97Xg was made by this physically motivated criterion, and the usefulness of the resulting functional was assessed on noble-gas dimer cation dissociation (Figs. 2 and 3). Even with only very conservative changes made to the functional form of the attenuator, terf was able to signifi- cantly increase the accuracy (relative to its predecessor xB97X) of the dissociation asymptote associated with the self-interaction problem (SIE) (Fig. 3). In the case of He^þ₂ the valence density lies so close to the neighboring atom that it seems challenging to reach a compromise between thermochemistry and exact exchange within a transferable range-separated exchange functional but the results for Ar^þ₂are encouraging. The solid thermochemical per- formance of the parent functional seemed more-or-less conserved, and so we investigated a little further. We make the approximation that the GGA parameters are unchanged between erf (xB97X) and terf attenuated Coulomb interactions. Undoubtedly this should be improved upon in future work and the literature already describes many ways this may be done.

A few more sets of attenuator parameters were obtained by maximizing the least-squared overlap of a terf attenuator with those of

x

B97 andxB97X varying r0;cxfor a given

x

. Noble-gas dissociation curves and atomization energies were calculated for the standard G2 [44] thermochemical test set in the -5.02

-5 -4.98 -4.96 -4.94 -4.92 -4.9 -4.88

0 1 2 3 4 5 6 7

Energy (Eh)

R(He-He) Angstrom Terf (1.016) wB97X Correct Asyptote

Fig. 2. Dissociation of He^þ₂.

-1054.56 -1054.54 -1054.52 -1054.50 -1054.48 -1054.46 -1054.44

2 3 4 5 6 7

Energy (Eh)

R(Ar-Ar) Angstrom Terf (1.016) wB97X

Correct Asymptote (Terf) Correct Asymptote (Erf)

Fig. 3. Dissociation of Ar^þ₂.

Table 1

Mean absolute error (kcal/mol) of G2 set atomization energies and errors of dimer cation asymptotes for various functionals.

x(a.u.) r0(a.u.) cx MAE He^þ₂ error Ne^þ₂error Ar^þ₂error Kr^þ₂error

BLYP^c – – – 87.55 80.37 48.51 49.0

0:3^a 0 0.1577 2.09 38.9 34.4 14.0 11.2

0:4^b 0 0 2.53 35.1 33.3 9.7 7.7

1 1.48 0.2395 5.71 31.7 29.3 6.7 4.2

1.016 1.2 0.1577 4.02 24.9 24.7 3.0 1.5

1.4 1.345 0.248 5.77 26.6 26.2 2.9 1.2

2 1.329 0.316 9.16 23.2 21.9 0.4 0.0

2 0.98 0.138 7.38 13.0 16.2 .03 0.0

a xB97X.

b xB97.

c Pure Becke 88 exchange[42]and LYP[43]correlation, errors in this row are upper bounds.

J.A. Parkhill et al. / Chemical Physics Letters 478 (2009) 283–286 285

6-311++G(3df,3pd) basis with a saturated quadrature grid. The purpose was not to search for parameters which would surpass xB97X, because gradient corrections are vital to thermochemical accuracy, but rather to document the balance between thermo- chemistry and correction of the self-interaction error (Table 1).

As expected the results were quite sensitive to the steepness of the attenuator and the amount of middle-range exchange, but note that the further this attenuator departs from erf the more severe becomes the GGA approximation. An accurate asymptote cannot be obtained by simply increasing the amount of exact exchange in the attenuator (because eventually the correlation part of the problem is disturbed), but by maximizing overlap with the estab- lished ones in the critical region we do obtain reduction of SIE with increasing ‘exact’ exchange. At the moment where the singly-occu- pied MO’s bond is breaking if this MO’s density around atom 2 is far enough from the bulk of this MO’s density around atom 1 so that it experiences Hartree exchange the density will localize on an atom (as it should physically). Smaller atoms cause greater dif- ficulty in general depending specifically on shell structure. The re- sults suggest that if the GGA enhancement factor of the functional were re-optimized for the modified attenuator it would possess thermochemistry much like xB97X with a significantly larger amount of exact exchange, and thus reduced self-interaction. Even in its current incarnation, the modified functional is accurate en- ough to be used in lieu of others for problems where ‘exact’ ex- change might be important.

4. Discussion and conclusions

Owing largely to the popularity of hybrid functionals, range sep- aration of exchange has become an intense area of research, and more flexible range separation may prove desirable[18]. Indeed this has already been done with the erf–gau type attenuator[32], although in this case this was done at the expense of abandoning a physically motivated choice of a parameter. An analytical formula [36]is available for the exact exchange energy with the terf atten- uation, and this expression has already been efficiently imple- mented in the publicly available release of the Q-CHEM package [45]. This paper provides the other building block, the short-range LDA exchange energy and a proof-of-concept GGA functional. Preli- minary results with unoptimized parameters indicate that the new functional may be a useful improvement and development should continue in the area of more general exchange attenuators. Further improvement over the functional developed here might realized through complete reoptimization[13]. Alternatively one could de- rive the corresponding PBE type GGA-functional [46,47,41]. In either case, the path is clear and only limited by one’s curiosity. It will be especially interesting to see if the flexibility of the new attenuator can simultaneously describe ground state electronic structure and excited states. Special attention should be paid to the size of the chromophore relative to the scale of the attenuator, and the distance over which electron density is redistributed. This direction is currently being pursued in our group.

Acknowledgements

This work was supported by the Director, Office of Energy Research, Office of Basic Energy Sciences, Chemical Sciences Division of the US Department of Energy under Contract DE-AC0376SF00098, and by a Grant from the SciDACProgram.

References

[1] J.C. Slater, Phys. Rev. 81 (1951) 385.

[2] W. Kohn, L.J. Sham, Phys. Rev. 140 (1965) A1133.

[3] P.W. Ayers, J. Math. Phys. 46 (2005) 062107.

[4] A.D. Becke, J. Chem. Phys. 98 (1993) 1372.

[5] S. Yanagisawa, T. Tsuneda, K. Hirao, J. Chem. Phys. 112 (2000) 545.

[6] J.P. Perdew, A. Zunger, Phys. Rev. B 23 (1981) 5048.

[7] Y. Zhang, W. Yang, J. Chem. Phys. 109 (1998) 2604.

[8] J.P. Perdew, R.G. Parr, M. Levy, J.L. Balduz, Phys. Rev. Lett. 49 (1982) 1691.

[9] O.A. Vydrov, G.E. Scuseria, J.P. Perdew, A. Ruzsinszky, G.I. Csonka, J. Chem. Phys.

124 (2006) 094108.

[10] R.O. Jones, O. Gunnarsson, Phys. Rev. Lett. 55 (1985) 107.

[11] M.R. Pederson, R.A. Heaton, C.C. Lin, J. Chem. Phys. 80 (1984) 1972.

[12] A. Ruzsinszky, J.P. Perdew, G.I. Csonka, O.A. Vydrov, G.E. Scuseria, J. Chem.

Phys. 126 (2007) 104102.

[13] J.-D. Chai, M. Head-Gordon, J. Chem. Phys. 128 (2008) 084106.

[14] A.D. Dutoi, M. Head-Gordon, Chem. Phys. Lett. 422 (2006) 230.

[15] O.A. Vydrov, G.E. Scuseria, J.P. Perdew, J. Chem. Phys. 126 (2007) 154109.

[16] A. Dreuw, M. Head-Gordon, Chem. Rev. 105 (2005) 4009.

[17] J.P. Perdew, V.N. Staroverov, J. Tao, G.E. Scuseria, Phys. Rev. A (Atomic Molecul.

Opt. Phys.) 78 (2008) 052513.

[18] A.V. Krukau, G.E. Scuseria, J.P. Perdew, A. Savin, J. Chem. Phys. 129 (2008) 124103.

[19] A.D. Becke, J. Chem. Phys. 112 (2000) 4020.

[20] P.M.W. Gill, R.D. Adamson, J.A. Pople, Molecul. Phys. 88 (1996) 1005.

[21] T. Leininger, H. Stoll, H.-J. Werner, A. Savin, Chem. Phys. Lett. 275 (1997) 151.

[22] R. Pollet, A. Savin, T. Leininger, H. Stoll, J. Chem. Phys. 116 (2002) 1250.

[23] H. Iikura, T. Tsuneda, T. Yanai, K. Hirao, J. Chem. Phys. 115 (2001) 3540.

[24] J.-W. Song, T. Hirosawa, T. Tsuneda, K. Hirao, J. Chem. Phys. 126 (2007) 154105.

[25] D. Jacquemin, E.A. Perpète, O.A. Vydrov, G.E. Scuseria, C. Adamo, J. Chem. Phys.

127 (2007) 094102.

[26] J.-W. Song, M.A. Watson, A. Nakata, K. Hirao, J. Chem. Phys. 129 (2008) 184113.

[27] B.M. Wong, J.G. Cordaro, J. Chem. Phys. 129 (2008) 214703.

[28] B.G. Janesko, T.M. Henderson, G.E. Scuseria, Phys. Chem. Chem. Phys. 11 (2009) 443.

[29] P.M.W. Gill, R.D. Adamson, Chem. Phys. Lett. 261 (1996) 105.

[30] T.M. Henderson, A.F. Izmaylov, G.E. Scuseria, A. Savin, J. Chem. Phys. 127 (2007) 221103.

[31] J. Toulouse, F. Colonna, A. Savin, Phys. Rev. A 70 (2004) 062505.

[32] J.-W. Song, S. Tokura, T. Sato, M.A. Watson, K. Hirao, J. Chem. Phys. 127 (2007) 154109.

[33] A. Savin, H.-J. Flad, Int. J. Quant. Chem. 56 (1995) 327.

[34] R. Baer, E. Livshits, D. Neuhauser, Chem. Phys. 329 (2006) 266.

[35] Y. Akinaga, S. Ten-no, Chem. Phys. Lett. 462 (2008) 348.

[36] A.D. Dutoi, M. Head-Gordon, J. Phys. Chem. A 112 (2008) 2110.

[37] M.A. Rohrdanz, J.M. Herbert, J. Chem. Phys. 129 (2008) 034107.

[38] J.-D. Chai, M. Head-Gordon, Chem. Phys. Lett. 467 (2008) 176.

[39] A. Fetter, D. Walecka, Quantum Theory of Many-Particle Systems, Dover, 2003.

[40] J. Toulouse, A. Savin, H. Flad, Int. J. Quant. Chem. 100 (2004) 1047.

[41] J. Toulouse, F. Colonna, A. Savin, J. Chem. Phys. 122 (2005) 014110.

[42] A.D. Becke, Phys. Rev. A 38 (1988) 3098.

[43] C. Lee, W. Yang, R.G. Parr, Phys. Rev. B 37 (1988) 785.

[44] L.A. Curtiss, K. Raghavachari, P.C. Redfern, J.A. Pople, J. Chem. Phys. 106 (1997) 1063.

[45] Y. Shao et al., Phys. Chem. Chem. Phys. 8 (2006) 3172.

[46] J. Heyd, G.E. Scuseria, M. Ernzerhof, J. Chem. Phys. 118 (2003) 8207.

[47] T.M. Henderson, B.G. Janesko, G.E. Scuseria, J. Chem. Phys. 128 (2008) 194105.

286 J.A. Parkhill et al. / Chemical Physics Letters 478 (2009) 283–286