Gradient-dependent upper bound for the exchange-correlation energy and application to density functional theory
aa r X i v : . [ c ond - m a t . o t h e r] J a n Gradient-dependent upper bound for the exchange-correlation energy and applicationto density functional theory
Lucian A. Constantin and Aleksandrs Terentjevs
Center for Biomolecular Nanotechnologies @UNILE,Istituto Italiano di Tecnologia, Via Barsanti, I-73010 Arnesano, Italy
Fabio Della Sala and Eduardo Fabiano
Istituto Nanoscienze-CNR, Via per Arnesano 16, I-73100 Lecce, Italy andCenter for Biomolecular Nanotechnologies @UNILE,Istituto Italiano di Tecnologia, Via Barsanti, I-73010 Arnesano, Italy
We propose a simple gradient-dependent bound for the exchange-correlation energy (sLL), basedon the recent non-local bound derived by Lewin and Lieb. We show that sLL is equivalent tothe original Lieb-Oxford bound in rapidly-varying density cases but it is tighter for slowly-varyingdensity systems. To show the utility of the sLL bound we apply it to the construction of simplesemilocal and non-local exchange and correlation functionals. In both cases improved results, withrespect to the use of Lieb-Oxford bound, are obtained showing the power of the sLL bound.
The Lieb-Oxford (LO) bound , E xc [ n λ ] ≥ LO [ n λ ] with n λ ( r ) = λ n ( λ r ) being the uniformly scaled density and LO [ n λ ] = − . Z n / λ ( r ) d r , (1)is one of the main exact constraints for the exchange-correlation (XC) energy and is used in the construc-tion of many semilocal XC approximations . Severalworks have concerned the tightening of this bound. Chanand Handy proposed to reduce the constant in LO to1.6358. Odashima and Capelle used numerical analysisto further reduce it to 1.444. Conversely, considering thelow-density limit of the spin-unpolarized uniform electrongas, it is found that the constant cannot be smaller than1.437 . Consequently the value 1.4 has been employedto build some recent XC functionals .Very recently Lewin and Lieb derived a differentbound, having a density dependence, i.e. E xc [ n λ ] ≥ LL [ n λ ] where LL [ n λ ] = − . Z n / λ ( r ) d r − (2) − . (cid:20)Z |∇ n λ | d r (cid:21) / (cid:20)Z n / λ d r (cid:21) / . Note that both the LO and LL bounds scale as the ex-change energy under a uniform scaling of the density.However, the LL bound, in the present form, is not verypractical for development in density functional theory.Thus, in this communication we consider a simplifiedform for such a bound which, at the same time, tight-ens it, such as to provide a useful tool for XC functionaldevelopment.To start we note that the LL bound is tighter than theLO one only for slowly-varying densities, where |∇ n λ | issmall. Let first consider the XC energy in this regime(spin-unpolarized case), in the low-density limit ( λ → becomes exact E xc [ n λ ] = − . Z n / λ d r − . Z n / λ s d r + Z n λ βt λ d r ≥ (3) ≥ − . Z n / λ d r − . Z n / λ s d r ≥≥ − . Z n / λ d r − . Z n / λ s α d r , (4)for any α ≤
2, and where s = |∇ n λ | / [2(3 π ) / n / λ ] and t λ = |∇ n λ | / [2(3 /π ) / n / λ ] are the reduced gradients forexchange and correlation respectively (note that s is in-variant under the uniform scaling while t λ ∝ λ / ) and β ≥ β is known only in the high-densitylimit, and in the metallic range (i.e. 2 ≤ r s ≤
6, with r s = [3 / (4 πn )] / being the bulk parameter). The lastinequality uses the fact that s ≤ E xc [ n λ ] ≥ − . Z n / λ d r − . Z n / λ s / d r . (5)This provides a slightly tighter bound than LL (notmathematically rigorous but in the spirit of Ref. 13, asexplained below), and it can be used together with Eq.(4) to obtain a more useful formula for the bound.In fact, a comparison of Eqs. (4) and (5) suggests thatthe correct bound could be sLL [ n λ ] = − . Z n / λ d r − C Z n / λ s / d r , (6)with 0 . ≤ C ≤ .
52 (to be fixed later). A similaridea is also proposed in Ref. 20. This ansatz is alsosupported by a further analysis of the XC energy underthe uniform scaling to the low-density limit at any valueof s (note that instead Eq. (4) holds only for s < E xc [ n λ ] = E LDAxc [ n λ ] + E GE xc [ n λ ] + δE x [ n λ ] + δE c [ n λ ],where LDA denotes the local density approximation, GE δE x and δE c indicate higher-order exchange and cor-relation contributions. Because in this limit t → δE c ≈
0. Moreover, following Ref. 8, we can assumethat E GE xc [ n λ ] ≈ λ →
0. Hence, E xc ≈ E LDAxc [ n λ ] + δE x [ n λ ] = (7)= − . Z n / λ d r − C X Z s ≥ n / λ [ f exactx ( s ) − d r , with C x = (3 / /π ) / and f exactx is the exact ex-change enhancement factor. Noting that at large gra-dients f exactx − ≥
0, Eq. (7) shows that a tighter LObound may fail in the low-density regime (where bothLDA correlation and exchange scale similarly, as ∼ /r s )and an s -dependent term is instead required.To fix the parameter C , we consider a tough case forthe semilocal sLL bound: the beryllium isoelectronic se-ries up to the nuclear charge Z = 140. This is in fact arapidly-varying high-density limit case with strong cor-relation ( E LDAc ≥ E exactc → −∞ for highly chargedions ). While the original LO bound gives LO/E xc ≈ .
9, the LL bound is much looser, having
LL/E xc ≈ . ≈ LOfor such a case (that is representative for rapidly-varyingdensity regime). In this way we obtain C = 0 . which was verifiedfor these cases. Hence, we consider here a few additionalpossibilities, based on interesting model systems.For slowly-varying density cases the sLL bound issurely verified in the low-density limit, because C =0 . ≥ . ≤ r s ≤
20, as well astwo interacting jellium slabs of thickness and separationdistance equal to the Fermi wave length (similar resultsare obtained also for other thicknesses and distances) and1 ≤ r s ≤
20. The latter describe several density regimes,ranging from slowly- (for thick and close slabs) to veryrapidly-varying ones (thin slabs).In all cases the exchange energy is computed exactly,whereas the correlation energy is estimated using the JSfunctional (similar results are given also by PBE andTPSS ). The LO, sLL, and LL have been calculated us-ing exact exchange densities. These systems are relevantto understand the significance on the bound of quantumoscillations in ordinary matter. Figure 1 shows that inall cases the sLL bound is respected ( sLL/E xc ≥
1) and s R a t i o sLL/E xc LO/E xc LL/E xc LO/E xc UEGLL/E xc UEG s FIG. 1: The ratios
LO/E xc , LL/E xc , and sLL/E xc versus1 /r s , for 40 e − jellium spheres (left panel) and for two in-teracting jellium slabs (right panel). Also shown are LO/E xc and LL/E xc = sLL/E xc , for the uniform electron gas (UEG). it is tighter than LO. Instead, the non-local LL boundis always rather loose. Finally, comparison with uniformelectron gas results shows that quantum oscillations andsurface effects have a negligible effect on the LO bound,which does not depend on the gradient, while they givea little weakening of the semilocal sLL bound.As a different example we consider the one-electrondensities n H ( r ) = e − r π ; n G ( r ) = e − r π / . (8)These are models for rapidly-varying density systems andhave been shown to be relevant for bonding propertiesof molecules . The calculations are analytical and wereadily find: LO/E XC = 1 . LL/E XC = 1 .
82, and sLL/E xc = 1 .
56 for H;
LO/E XC = 1 . LL/E XC =1 .
80, and sLL/E xc = 1 .
54 for G. Thus, the LO and sLLbounds have the same quality, both outperforming theLL bound. The coincidence of LO and sLL for thesemodel systems is reminiscent of the fact that, to fix Cin Eq. (6), we constrained sLL to reproduce LO for theBe series, which presents rapidly-varying densities. How-ever, notably in the present case there is no contributionfrom correlation.To conclude our analysis we consider finally the op-posite case of slowly-varying strong-correlated systems.One example of such systems is given by the point-charge-plus-continuum (PC) model . At small re-duced gradients s ≤
1, the XC energy is E P Cxc [ n ] = − . Z n / d r + 0 . Z n / s d r . (9)We observe that the local term is exactly the one enteringin the LL and sLL bounds and is very close to the numer-ical estimation of Odashima and Capelle , although thislatter bound is formally broken for s →
0. This showsthat this model is the true limit for constant densities.On the other hand, because the second-order XC contri-bution is positive, all the examined bounds are valid forthe PC model for 0 < s ≤ s F x PBEsolPBEsol(sLL)sLLLO
FIG. 2: Exchange enhancement factors versus s , for PBEsoland PBEsol(sLL) Also shown, are the asymptotes of sLL andLO bounds. Another interesting example of slowly-varying stronglycorrelated system can be obtained by considering thestrong-correlation scaling of the density , n λ ( r ) = λ − n ( λ − r ), in the limit λ → ∞ . Under this scalingthe reduced gradients behave as s λ ( r ) = λ − / s ( λ − r )and t λ ( r ) = λ − / t ( λ − r ). The XC energy is dominatedby the local terms, while the second-order contributionsdecay as λ − / . On the other hand, we have sLL = − . λ / R n / d r − . λ / R n / s / d r . Therefore,for λ → ∞ sLL provides a good bound to the XC energy(note that the local part of sLL is almost the exact limitfor the local XC energy). Moreover, we note that in thislimit sLL ≤ LO , being a tighter bound for the XC energyof slowly-varying strongly-correlated systems.To show examples of the practical utility of the sLLbound we consider its application in the construction ofa GGA and an hyper-GGA functional. GGA . We consider the new sLL bound in the PBEsol XC functional, which is a non-empirical GGA with thecorrect second-order exchange coefficient. The value of κ = 0 .
804 in the exchange enhancement factor F x =1 + κ − κ/ [1 + ( µ/κ ) s ] is replaced by κ = 0 .
559 + 0 . s / . (10)in order to respect the sLL bound at any point of thespace. For the correlation part we use β = 0 .
045 (notethat in PBEsol β = 0 . .In Fig. 2, we show a comparison of the exchange en-hancement factors. Note that any exchange functionalthat recovers the local LO bound (e.g. PBEsol, TPSSmeta-GGA ), satisfies also the local sLL bound. Instead,PBEsol with the sLL bound clearly violates the local LObound.To test the effect of the sLL bound we performed dif-ferent tests for solids and molecules. In Table I we reportthe mean absolute errors (MAEs) for solid-state testson 24 bulk solids (see computational details). We seethat using the sLL bound improves the results: latticeconstants and cohesive energies are improved with re-spect PBEsol, while bulk moduli are almost comparable.We recall that PBEsol is one of the most used semilocalfunctionals in solid-state physics and, one of the best forlattice constants and bulk moduli . AM05 performsslightly worse than PBEsol (MAEs on this test sets are32 m˚A and 8.6 Gpa, respectively). A PBEsol functional TABLE I: Mean absolute error (MAE) for the lattice con-stants (m˚A), bulk moduli (GPa), and cohesive energies(ev/atom) as obtained from the PBEsol functional with theoriginal LO bound, and with the new sLL bound, for varioustypes of solids. Reference data from Ref. 29.PBEsolLO sLLLattice constants (m˚A)simple metals 37.6 31.8transition metals 24.9 23.2semiconductors 29.6 31.8ionic solids 20.4 22.2insulators 6.5 6.3total MAE 24.5 23.8bulk moduli (GPa)simple metals 1.0 0.9transition metals 18.7 18.5semiconductors 7.8 8.1ionic solids 4.1 5.0insulators 6.7 6.7total MAE 7.7 7.9cohesive energies (eV)simple metals 0.15 0.13transition metals 0.71 0.64semiconductors 0.31 0.26ionic solids 0.10 0.06insulators 0.65 0.60total MAE 0.37 0.33 with κ = 0 .
559 is definitely worse , having MAEs of35 m˚A, 8.1 GPa and 0.51 eV/atom. This shows the im-portance of the non-local term of the sLL bound.In Table II we report the MAEs for several tests onmolecular properties that are relevant for semilocal func-tionals. In this case PBEsol with the new sLL bound isbetter than PBEsol for all tests (the only exception isthe OMRE test) with an overall relative MAE of 0.94.The improvement traces back to the larger non-localityof the sLL bound with respect to the LO one, whichprovides a better description when large reduced gradi-ents are involved . Note that a further improvementfor molecular properties can be obtained by replacingthe PBEsol correlation with the zPBEsol one , whichis not modifying the sLL bound behavior but improvesmost of the energetic properties . Hyper-GGA.
We consider the hyper-GGA correla-tion functional E c [ n ] = Z e T P SSc ( r ) X ( r ) − e T P SSx ( r ) ( X ( r ) − e x ( r )) d r , (11)where e T P SSc and e T P SSx are the TPSS correlation andexchange energy densities respectively, e x ( r ) is the ex- TABLE II: Mean absolute errors (MAE) in kcal/mol for en-ergy tests, and in m˚A for geometry tests, for several moleculesbenchmark, as obtained from the PBEsol functional with theoriginal LO bound, and with the new sLL bound.PBEsoltest LO sLLatomiz. energy (AE6) 34.5 32.1atomiz. energy (W4) 21.4 20.0atomiz. energy (W4-MR) 35.5 33.93 d -metals AE (TM10AE) 18.3 17.6Au clusters AE (AUnAE) 3.6 3.1reaction ener. (OMRE) 8.1 10.7reaction ener. (BH76RC) 6.3 6.0barrier heights (BH76) 12.2 11.83 d -metals RE (TMRE) 9.9 9.1ionization pot. (IP13) 2.3 2.1proton aff. (PA12) 1.5 1.5difficult cases (DC9/12) 82.9 78.3isomerization ener. (ISOL6) 1.7 1.7interfaces (SI12) 3.8 3.5hydrogen bonds (HB6) 1.6 1.2dihydrogen bonds (DHB23) 1.8 1.5dipole-dipole (DI6) 1.0 0.7bond lengths (MGBL19) 10.0 9.9total energy (AE17) 250.3 246.6 P i MAE i / MAE
PBEsoli act exchange energy density (we use only the conven-tional gauge ), and X ( r ) is a local upper bound forthe e T P SSx . We consider two possibilities for this bound: R X ( r ) d r = LO and R X ( r ) d r = sLL . The correspond-ing hyper-GGA functionals are labeled hyper(LO) andhyper(sLL). These correlation functionals are exact forone-electron systems, are size-consistent, and even if theyare not so realistic , they represent very good models tounderstand local hybrids, and to compare the quality ofthe LO and sLL bounds. The correlation energies of theatoms from He to Ar as well as of jellium spheres (with2, 8, 18, 20, 34, 40, 58, 92 , and 106 electrons and var-ious r s ), computed with these functionals, are reportedin Table III. The hyper(sLL) functional systematicallyimproves over hyper(LO), showing that the sLL boundmay be an interesting tool in the construction of non-local functionals.In conclusion, we showed a new simple gradient-dependent bound (sLL), which has the same quality asthe Lieb-Oxford (LO) bound for systems where the den-sity varies rapidly but is significantly tighter for slowly-varying density cases. It is also the exact limit for slowly-varying strongly correlated systems. To indicate the util- ity of this new bound we applied it to both semilocal andnon-local DFT, showing how it can be employed to con- TABLE III: Mean absolute relative errors (%) for correlationenergies of atoms from He to Ar and of magic jellium spheres(js) (average for spheres of 2, 8, 18, 20, 34, 40, 58, 92 and106 atoms) with different values of r s . The reference datafor jellium spheres are taken from Ref. 37. Accurate atomiccorrelation energies are taken from Ref. 38.PBE TPSS hyper(LO) hyper(sLL)atoms 6.3 4.9 5.5 5.4js r s = 1 6.8 6.0 5.6 5.5js r s = 2 7.9 6.6 6.1 5.9js r s = 3 .
25 7.4 5.9 5.4 5.2js r s = 4 7.1 5.4 5.0 4.8js r s = 5 .
62 7.8 5.7 5.4 5.3Average js 7.4 5.9 5.5 5.3 struct improved exchange and correlation functionals. Inparticular, we found that the sLL bound improves the ac-curacy of the PBEsol functional for both solid-state andmolecules properties. In addition, we showed that thesLL bound can be fruitfully employed in the constructionof hyper-GGAs. Even if the improvements of these func-tionals are in general small they are rather systematic.This proves the quality of the sLL bound and suggestsits importance for development of new functional forms.
Computational details.
Equilibrium lattice con-stants, bulk moduli, and cohesive energies have been cal-culated for 24 solids, including Al, Ca, K, Li, Na (simplemetals); Ag, Cu, Pd, Rh, V (transition metals); LiCl,LiF, MgO, NaCl, NaF (ionic solids); AlN, BN, BP, C(insulators); GaAs, GaP, GaN, Si, SiC (semiconductors).Reference data to construct this set were taken from Ref.29. These calculations have been performed with theVASP program using PBE-PAW pseudopotentials. AllBrillouin zone integrations were performed on Γ-centeredsymmetry-reduced Monkhorst-Pack k -point meshes, us-ing the tetrahedron method with Bl¨och corrections. Forall the calculations a 24 × × k -mesh grid was appliedand the plane-wave cutoff was chosen to be 30% largerthan maximum cutoff defined for the pseudopotential ofeach considered atom.Calculations for Table II have been performed with theTURBOMOLE program package using a def2-TZVPPbasis set . For details about the molecular test sets,see Refs. 43.Hyper-GGA calculations are non-self-consistent, usingaccurate exact exchange orbitals and densities. Acknowledgments.
We thank Prof. Andreas Savinfor useful discussions. E.H. Lieb and S. Oxford, Int. J. Quantum Chem. , 427(1981). J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996). J.P. Perdew, A. Ruzsinszky, G.I. Csonka, O.A. Vydrov,G.E. Scuseria, L.A. Constantin, X. Zhou, and K. Burke,Phys. Rev. Lett. , 136406 (2008); ibid. , 039902(2009) (E). L. A. Constantin, E. Fabiano, S. Laricchia, F. Della Sala,Phys. Rev. Lett. , 186406 (2011). E. Fabiano, L.A. Constantin, and F. Della Sala, Phys. Rev.B. , 113104 (2010). Z. Wu and R.E. Cohen, Phys. Rev. B , 235116 (2006). J. Tao , J.P. Perdew, V.N. Staroverov, and G.E. Scuseria,Phys. Rev. Lett. , 146401 (2003). J.P. Perdew, A. Ruzsinszky, G.I. Csonka, L.A. Constantin,and J. Sun, Phys. Rev. Lett. , 026403 (2009). L.A. Constantin, E. Fabiano, F. Della Sala, J. Chem. The-ory and Comput. , 2256 (2013). L.A. Constantin, E. Fabiano, F. Della Sala, Phys. Rev. B , 125112 (2013). J.M. del Campo, J.L. G´azquez, S. Trickey, A. Vela, Chem.Phys. Lett. , 179 (2012). G.K.-L. Chan and N.C. Handy, Phys. Rev. A , 3075(1999). M.M. Odashima and K. Capelle, J. Chem. Phys. ,054106 (2007). J.P. Perdew, in
Electronic structure of solids ’91 , edited byP. Ziesche and H. Eschrig (Akademie Verlag, Berlin, 1991). Y. Zhao and D.G. Truhlar, J. Chem. Phys. , 184109(2008). J. Sun, B. Xiao, and A. Ruzsinszky, J. Chem. Phys. ,051101 (2012). R. Haunschild, M.M. Odashima, G.E. Scuseria, J.P.Perdew, and K. Capelle, J. Chem. Phys. , 184102(2012). M. Lewin and E.H. Lieb, arXiv:1408.3358 [math-ph]. C.D. Hu and D.C. Langreth, Phys. Rev. B , 943 (1986). D. V. Feinblum, J. Kenison, and K. Burke,arXiv:1410.7838 [physics.chem-ph]. R.C. Morrison, J. Chem. Phys. , 10506 (2002). L.A. Constantin, L. Chiodo, E. Fabiano, I. Bodrenko, F.Della Sala, Phys. Rev. B , 045126 (2011). L.A. Constantin, E. Fabiano, F. Della Sala, Phys. Rev. B , 233103 (2011). M. Seidl and J. P. Perdew, Phys. Rev. B , 5744 (1994). M. Seidl, J.P. Perdew, and M. Levy, Phys. Rev. A , 51(1999). M. Seidl, J.P. Perdew, and S. Kurth, Phys. Rev. Lett. ,5070 (2000). M. Seidl, J.P. Perdew, and S. Kurth, Phys. Rev. A ,012502 (2000). E. Fabiano and L.A. Constantin, Phys. Rev. A , 012511(2013). J. Sun, M. Marsman, G.I. Csonka, A. Ruzsinszky, P. Hao,Y.-S. Kim, G. Kresse, and J.P. Perdew, Phys. Rev. B ,035117 (2011); L. Schimka, J. Harl, and G. Kresse, J.Chem. Phys. , 024116 (2011); F. Tran, R. Laskowski, P.Blaha, and K. Schwarz, Phys. Rev. B , 115131 (2007); J.Paier, M. Marsman, K. Hummer, G. Kresse, I. C. Gerber,and J. G. ´Angy´an, J. Chem. Phys. , 154709 (2006); C.Kittel, Introduction to Solid State Physics, 8-th ed.; JohnWiley & Sons Inc., New York, 2005. R. Armiento and A. E. Mattsson, Phys. Rev. B , 085108(2005); A. E. Mattsson and R. Armiento, Phys. Rev. B ,155101 (2009). P. Haas, F. Tran, P. Blaha, L. S. Pedroza, A. J. R. daSilva, M. M. Odashima, and K. Capelle, Phys. Rev. B ,125136 (2010). E. Fabiano, L. A. Constantin, and F. Della Sala, J. Chem.Theory Comput. , 3548 (2011). P. Haas, F. Tran, P. Blaha, and K. Schwarz, Phys. Rev. B , 205117 (2011). Y. Zhang and W. Yang, Phys. Rev. Lett. , 890 (1998). L. A. Constantin, E. Fabiano, and F. Della Sala, J. Chem.Phys. , 194105 (2012). L. A. Constantin, E. Fabiano, A. Terentjevs, and F. DellaSala, arXiv:1411.1579 [cond-mat.other]. J. Tao, J.P. Perdew, L.M. Almeida, C. Fiolhais, and S.K´ummel, Phys. Rev. B , 245107 (2008). S.J. Chakravorty, S.R. Gwaltney, and E.R. Davidson, F. A.Parpia, and C. F. Fischer, Phys. Rev. A , 3649 (1993). G. Kresse, J. Furthmuller, Phys. Rev. B 54, 11169 (1996). F. Weigend, F. Furche, R. Ahlrichs, J. Chem. Phys. 119,12753 12763 (2003). F. Weigend, R.; Ahlrichs, Phys. Chem. Chem. Phys. 7,3297 (2005). J. B. Lynch and D. G. Truhlar, J. Phys. Chem. A , 8996(2003); L. Goerigk and S. Grimme, Phys. Chem. Chem.Phys. , 6670 (2011); E. Fabiano, L. A. Constantin, andF. Della Sala, J. Chem. Phys. , 194112 (2011); E. Fabi-ano, L. A. Constantin, and F. Della Sala, Int. J. QuantumChem. , 673 (2013); R. Peverati and D. G. Truhlar, J.Chem. Theory Comput. , 2310 (2012); S. Luo, Y. Zhao,and D. G. Truhlar, Phys. Chem. Chem. Phys. 13 , 13683(2011); Y. Zhao and D. G. Truhlar, J. Chem. Theory Com-put. , 415 (2005); E. Fabiano, L. A. Constantin, and F.Della Sala, J. Chem. Theory Comput. , 3151 (2014);Y. Zhao and D. G. Truhlar, J. Chem. Phys.,125