Asymptotic correction schemes for semilocal exchange-correlation functionals
Chi-Ruei Pan,1Po-Tung Fang,1and Jeng-Da Chai1,2,*
1Department of Physics, National Taiwan University, Taipei 10617, Taiwan
2Center for Theoretical Sciences and Center for Quantum Science and Engineering, National Taiwan University, Taipei 10617, Taiwan (Received 30 August 2012; published 16 May 2013)
Aiming to remedy the incorrect asymptotic behavior of conventional semilocal exchange-correlation (XC) density functionals for finite systems, we propose an asymptotic correction scheme, wherein an exchange density functional whose functional derivative has the correct (−1/r) asymptote can be directly added to any semilocal density functional. In contrast to semilocal approximations, our resulting exchange kernel in reciprocal space exhibits the desirable singularity of the type O(−1/q2) as q → 0, which is a necessary feature for describing the excitonic effects in nonmetallic solids. By applying this scheme to a popular semilocal density functional, PBE [Perdew, Burke, and Ernzerhof,Phys. Rev. Lett.77, 3865 (1996)], the predictions of the properties that are sensitive to the asymptote are significantly improved, while the predictions of the properties that are insensitive to the asymptote remain essentially the same as PBE. Relative to the popular model XC potential scheme, our scheme is significantly superior for ground-state energies and related properties. In addition, without loss of accuracy, two closely related schemes are developed for the efficient treatment of large systems.
DOI:10.1103/PhysRevA.87.052510 PACS number(s): 31.15.es, 31.15.ee, 71.15.Mb
I. INTRODUCTION
Over the past two decades, Kohn-Sham density-functional theory (KS-DFT) [1,2] has been one of the most power- ful theoretical methods for the ground-state properties of large electronic systems. Its time-dependent extension, time- dependent density-functional theory (TDDFT) [3–5], has gradually become popular for the study of excited-state and time-dependent properties.
In KS-DFT, the exact exchange-correlation (XC) den- sity functional Exc[ρ] remains unknown and needs to be approximated. Accurate density-functional approximations to Exc[ρ] have been successively developed to extend the applicability of KS-DFT to a wide variety of systems. Despite the recent advances in the orbital-dependent density-functional approach [6], semilocal density functionals remain popular due to their computational efficiency for large systems and reason- able accuracy for applications governed by short-range XC effects [7]. However, due to the associated several qualitative failures, semilocal functionals can produce erroneous results in situations where the accurate treatment of nonlocality of the XC hole is important [8–10].
One of the important and long-standing subjects in KS-DFT is the asymptotic behavior of the XC potential vxc(r) = δExc[ρ]/δρ(r). For finite systems, the exact vxc(r) exhibits the Coulombic (−1/r) decay as r → ∞ [11–14]. However, due to the severe self-interaction error (SIE) [15], the XC potential of semilocal functionals fails to describe the correct asymptotic behavior, yielding qualitatively incorrect predictions for the properties sensitive to the asymptote, such as the vertical ionization potentials and high-lying (Rydberg) excitation energies of atoms and molecules [16,17].
Currently, perhaps the most successful density-functional methods in practice to improve the asymptote of the XC potential are provided by the long-range corrected (LC) hybrid
*Author to whom correspondence should be addressed:
scheme [18–28] and asymptotically corrected (AC) model potential scheme [29–39]. For the LC hybrid scheme, the nonlocal Hartree-Fock (HF) exchange for the long-range electron-electron interactions is added to a semilocal func- tional. Therefore, the LC hybrid scheme can be impractical for very large systems due to the inclusion of the long-range HF exchange (which significantly increases the computational cost relative to the semilocal functional). By contrast, for the AC model potential scheme, an AC XC potential is directly modeled and added to a semilocal functional, maintaining the similar cost as the semilocal functional. In principle, a model XC potential should be a functional derivative of some Exc[ρ]. However, as a number of popular model potentials are found not to be functional derivatives, several necessary conditions for a functional derivative can be violated [40].
Besides, as these model potentials are not variationally stable, the associated XC energies are not uniquely defined, and properties obtained from these model potentials need to be carefully interpreted [40,41]. Recently, we have examined the performance of the LC hybrid scheme and AC model potential scheme on a very wide range of applications [42]. Despite its computational efficiency, the popular model potential scheme can exhibit severe errors in the calculated ground-state energies and related properties, due to the lack of Exc[ρ].
On the other hand, for a system of N electrons, the Fermi- Amaldi (FA) XC functional [43],
ExcFA[ρ] = − 1 2N
! ! ρ(r)ρ(r$)
|r − r$| dr dr$, (1) which is simply (−1/N) times the Hartree energy functional, appears to be the simplest XC functional whose functional derivative has the correct (−1/r) asymptote. However, there are several problems with the FA model [44]. While the FA XC potential is correct in the asymptotic region, it is inaccurate elsewhere. Besides, due to its delocalized XC hole, the FA model is not size consistent, where the energy of a system composed of two or more well-separated subsystems
is not identical to the sum of the energies of the separate subsystems [45].
In this work, we propose an AC scheme for any system composed of atoms (e.g., atoms, molecules, and solids), wherein a modified FA XC functional, which is size consistent in the calculated energy and whose functional derivative has the correct (−1/r) asymptote, can be directly added to any semilocal functional. Without loss of accuracy, two related efficient schemes are also developed for large systems.
II. THEORETICAL METHODS A. LFA scheme
By partitioning and localizing a modified FA XC hole to the atoms in a system, we propose the “localized”FA (LFA) exchange functional,
ExLFA[ρα,ρβ] = − "
σ=α,β
"
A
1 2NA,σ
×
! !
ρA,σ(r)ρA,σ(r$)erf(ω|r − r$|)
|r − r$| dr dr$, (2) to resolve the size-inconsistency issue associated with the FA model (see Supplemental Material [46]). Here, the second sum is over all the atoms in the system, ρA,σ(r) is the σ -spin (σ = α for spin up or β for spin down) electron density associated with the atom A,
ρA,σ(r) = wA(r)ρσ(r), (3) and the weight function wA(r), ranging between 0 and 1, is of the Hirshfeld type [47,48]:
wA(r) = ρA0(r)
#
BρB0(r), (4)
where ρ0A(r) is the spherically averaged electron density computed for the isolated atom A. NA,σ is the number of the σ -spin electrons associated with the atom A,
NA,σ =
!
ρA,σ(r)dr, (5)
and the long-range interelectron repulsion operator erf(ω|r − r$|)/|r − r$| is to retain the correct asymptotic behavior without the (unneeded) energy contribution from the complementary short-range operator, where ω is a parameter defining the range of the operators. Due to the sum rule of#
AwA(r) = 1,
#
AρA,σ(r) = ρσ(r), and #ANA,σ = Nσ (the number of σ-spin electrons).
By taking the functional derivative of ExLFA[ρα,ρβ], the LFA exchange potential for σ -spin electrons is
vLFAx,σ(r) = δELFAx [ρα,ρβ] δρσ(r)
= −"
A
wA(r) NA,σ
!
ρA,σ(r$)erf(ω|r − r$|)
|r − r$| dr$. (6)
If the functional derivative of NA,σis also taken, an additional constant term
vLFAx,σ(∞) ≡"
A
1 2NA,σ2
×
! !
ρA,σ(r)ρA,σ(r$)erf(ω|r − r$|)
|r − r$| dr dr$ (7) should be added to Eq.(6). However, as will be shown later, this constant is of no consequence.
In the asymptotic limit, vx,σLFA(r) has the correct asymptotic form,
rlim→∞vx,σLFA(r) = −"
A
wA(r) NA,σ
!
ρA,σ(r$) 1
|r|dr$
= −1 r
"
A
wA(r) = −1
r. (8)
From Eqs.(2)and(6), we have ExLFA[ρα,ρβ] = 1
2
"
σ=α,β
!
ρσ(r)vx,σLFA(r)dr, (9) showing that the LFA exchange energy density per electron also has the correct (−1/2r) asymptote [29,49].
For the calculation of excitation energies using adiabatic linear-response TDDFT [3–5], the functional derivative of vx,σLFA(r) yields the LFA exchange kernel for σ -spin electrons,
fx,σLFA(r,r$) = δvLFAx,σ(r) δρσ(r$) =
δ2ExLFA[ρα,ρβ] δρσ(r)δρσ(r$)
= −erf(ω|r − r$|)
|r − r$|
"
A
wA(r)wA(r$) NA,σ . (10) In contrast to semilocal approximations, fx,σLFA(r,r$) in recip- rocal space has the correct long-wavelength O(−1/q2) diver- gence as q → 0, which is crucially important for the proper description of excitonic effects in nonmetallic solids [50–53].
We emphasize that this striking feature appears naturally from our fully nonlocal ELFAx [ρα,ρβ].
To improve its description of short-range XC effects, ELFAx [ρα,ρβ] is combined with a popular semilocal functional, PBE (Perdew, Burke, and Ernzerhof) [54]. However, this will produce a double-counting (DC) energy EDC, which needs to be removed. As the significant fraction of ExLFA[ρα,ρβ] should be from the core regions of the atoms, we presuppose that ρA,σ(r) is strictly localized at RA(the position of the atom A) in Eq.(2), which gives ρA,σ(r) ≈ NA,σδ(r − RA) for satisfying Eq.(5), to estimate EDC,
EDC= − "
σ=α,β
"
A
1 2NA,σ
×
! !
{NA,σδ(r − RA)}{NA,σδ(r$− RA)}
×erf(ω|r − r$|)
|r − r$| dr dr$
= − "
σ=α,β
"
A
1 2NA,σ
(NA,σ)2
$
r→RlimA
erf(ω|r − RA|)
|r − RA|
%
= − "
σ=α,β
"
A
NA,σ
2
$ 2ω√π
%
= − ω
√π
"
σ=α,β
"
A
NA,σ = − ω
√πN. (11)
This estimate is very accurate for systems with highly localized charges (e.g., HCl), and less accurate for systems with delocalized charges (e.g., benzene) (see Supplemental Material [46]). Our resulting LFA-PBE functional is given by ExcLFA-PBE= ExcPBE+ ExLFA− EDC. (12) Note that vDCx,σ(r) = δEDC/δρσ(r) = −ω/√
π is simply a constant, which can be absorbed into the constant vx,σLFA(∞).
As the KS potential is only defined within an arbitrary constant, without loss of generality, we require the KS potential to vanish asymptotically, which sets vx,σLFA(∞) + ω/√
π= 0.
Unlike the FA model, LFA-PBE is size consistent for any system composed of atoms (see Supplemental Material [46]).
Note that LFA-PBE (with ω = 0) reduces to PBE.
B. RILFA scheme
For systems composed of many atoms, LFA-PBE can be computationally unfavorable due to the numerical integration of many Hartree-like potentials in vx,σLFA(r). To resolve this com- putational bottleneck without loss of much accuracy, Eq.(6) can be efficiently evaluated by the resolution-of-identity (RI) approximation [55,56]. Following Ref. [56], ρA,σ(r) is expanded with an auxiliary basis set {gp(r)}, i.e., ρA,σ(r) ≈
˜ρA,σ(r) =#
papgp(r), where the expansion coefficients {ap} are given by Eq. (8) of Ref. [56] [with ρ(r) being replaced by ρA,σ(r)]. The RILFA exchange potential is evaluated by
vx,σRILFA(r) = −"
A
wA(r) NA,σ
!
˜ρA,σ(r$)erf(ω|r − r$|)
|r − r$| dr$. (13) From Eq. (10) of Ref. [56], the RILFA exchange energy is given by
ERILFAx = "
σ=α,β
"
A
$
− 1
NA,σ
×
! !
ρA,σ(r) ˜ρA,σ(r$)erf(ω|r − r$|)
|r − r$| dr dr$
+ 1
2NA,σ
! !
˜ρA,σ(r) ˜ρA,σ(r$)erf(ω|r − r$|)
|r − r$| dr dr$
% . (14) For a sufficiently large {gp(r)}, the RILFA scheme approaches the LFA scheme. Here, RILFA-PBE is defined by Eq.(12), with ExLFAbeing replaced by ExRILFA.
C. LFAs scheme
For very large systems, both the LFA and RILFA schemes may be impractical, compared to the efficient semilocal density-functional approach. Aiming to retain the correct (−1/r) asymptote with essentially no added computational cost relative to semilocal functionals, thestrict localization of
ρA,σ(r) at RA[i.e., ρA,σ(r) ≈ NA,σδ(r − RA)] is presupposed in Eq.(6), to define the LFAs exchange potential,
vLFAsx,σ (r) = −"
A
wA(r) NA,σ
!
{NA,σδ(r$− RA)}erf(ω|r − r$|)
|r − r$| dr$
= −"
A
wA(r)erf(ω|r − RA|)
|r − RA| . (15)
The asymptote of vx,σLFAs(r) remains correct,
rlim→∞vx,σLFAs(r) = −"
A
wA(r)1
|r| = −1
r. (16)
Based on Eq.(9), the LFAs exchange energy is given by ELFAsx = 1
2
"
σ=α,β
!
ρσ(r)vx,σLFAs(r)dr (17) to retain the correct (−1/2r) asymptote of the LFAs exchange energy density per electron. Although vLFAsx,σ (r) differs from the functional derivative of ELFAsx by a factor of 2, the prescribed LFAs scheme approaches the LFA scheme for a sufficiently small ω value, where vx,σLFAs(r) becomes an excellent approximation of vx,σLFA(r). Similarly, LFAs-PBE is defined by Eq.(12), with ExLFAbeing replaced by ExLFAs.
III. DEFINITION OF AN OPTIMAL ω VALUE For the exact KS-DFT, the minus highest occupied molecular orbital (HOMO) energy of a molecule should be the same as the vertical ionization potential (IP) of the molecule [11,12,57–59]. Therefore, the optimal ω values for LFA-PBE, RILFA-PBE, and LFAs-PBE are determined by fitting the predicted IPs (calculated by the minus HOMO energies) of 18 atoms and 113 molecules in the IP131 database to the corresponding experimental IPs [60]. All calculations are performed with a development version ofQ-CHEM3.2[61], using the 6-311 + + G(3df,3pd) basis set (and sufficiently large auxiliary basis sets for the RILFA scheme), unless noted otherwise. The error for each entry is defined as (error = theoretical value − reference value).
As shown in Fig.1, the minimum root-mean-square (rms) errors of LFA-PBE, RILFA-PBE, and LFAs-PBE for the IP131 database, which all occur at ω = 0.15 bohr−1, are more than three times smaller than the rms error of PBE (the ω = 0 case), reflecting the importance of the correct asymptote of XC potential here (see Supplemental Material [46]).
Adopting ω = 0.15 bohr−1 for all the LFA-corrected PBE functionals, the calculated IPs are plotted against the experimental values in Fig.2. As can be seen, the differences between the IPs calculated by RILFA-PBE and LFA-PBE are within 0.005 eV, and the differences between the IPs calculated by LFAs-PBE and LFA-PBE are within 0.1 eV.
Therefore, all the LFA-corrected PBE functionals yield very similar results, indicating that their XC potentials should be very similar (see Supplemental Material [46]). By contrast, the IPs calculated by PBE are seriously underestimated due to its incorrect asymptote.
Similar results are found when our LFA-related schemes are combined with LDA [62,63]. As both the LDA and PBE XC potentials decay exponentially in the asymptotic
rms
FIG. 1. (Color online) Root-mean-square (rms) errors of LFA- PBE, RILFA-PBE, and LFAs-PBE for the IP131 database [60]. The ω= 0 case corresponds to PBE.
region, their predicted IPs are similar, requiring essentially the same corrections from the LFA-related schemes (i.e., with the same optimal ω) (see Supplemental Material [46]). Therefore, ω= 0.15 bohr−1can be recommended as the optimal ω value, when the LFA-related schemes are combined with a local or semilocal functional whose functional derivative has the (incorrect) exponential asymptote.
IV. RESULTS AND DISCUSSION
Here, we examine the performance of the PBE and LFA- corrected PBE functionals (ω = 0.15 bohr−1) on various test sets, involving the reaction energies of 30 chemical reactions
FIG. 2. (Color online) Calculated versus experimental ionization potentials (IPs) for the IP131 database [60]. Inset I shows the differences between the IPs calculated by RILFA-PBE and LFA-PBE, while inset II shows the differences between the IPs calculated by LFAs-PBE and LFA-PBE. ω = 0.15 bohr−1 is adopted for all the LFA-corrected PBE functionals.
TABLE I. Mean absolute errors (in kcal/mol) of various test sets (see the text for details). The LB94 results are taken from Ref. [42].
(1 kcal/mol = 0.0434 eV.)
System PBE LFA-PBE RILFA-PBE LFAs-PBE LB94
Reaction (30) 4.38 4.48 4.47 4.42
G3/99 (223) 21.51 27.36 27.36 24.58 484.91
NHTBH (38) 8.62 8.71 8.71 8.66 93.94
HTBH (38) 9.67 9.69 9.69 9.69 44.31
S22 (22) 2.72 2.37 2.37 2.52 51.70
(a test set described in Ref. [23]), the 223 atomization energies of the G3/99 set [64–66], the 76 barrier heights (BHs) of the NHTBH38/04 and HTBH38/04 sets [67,68], the 22 noncovalent interactions of the S22 set [69], 19 valence excitation energies, and 23 Rydberg excitation energies. There are in total 393 pieces of data in our test sets, which are quite large and diverse. Unspecified detailed information of the test sets is given in Ref. [60]. For comparison, the results calculated by the van Leeuwen and Baerends 1994 (LB94) potential (a popular AC model XC potential) [29] are taken from Ref. [42]. Note that the LB94 potential is a linear combination of the LDA exchange potential, the LDA correlation potential, and a gradient-dependent exchange potential [e.g., see Eq. (55) of Ref. [29]]. Due to the inclusion of the gradient-dependent exchange potential, the LB94 potential is not a functional derivative [40,41]. In Ref. [42], the exchange energy from the LB94 exchange potential was evaluated by the popular Levy- Perdew virial relation [70] [e.g., see Eq. (1) of Ref. [42]], while the correlation energy from the LB94 correlation potential was directly evaluated by the LDA correlation energy functional.
As shown in TableI, the performance of the LFA-corrected PBE functionals is similar to that of PBE (see Supplemental Material [46]). As these properties are rather insensitive to the asymptote of the XC potential, our schemes do not affect the already good performance of PBE. By contrast, due to the lack of Exc[ρ], LB94 performs the worst. Therefore, one should avoid using the AC model potential scheme for the calculation of total energies and related properties.
For the valence and Rydberg excitation energies, we per- form adiabatic linear-response TDDFT calculations, using the 6-311(2+,2+)G** basis set, on five molecules: nitrogen gas (N2), carbon monoxide (CO), water (H2O), ethylene (C2H4), and formaldehyde (CH2O) on the experimental geometries taken from Ref. [71]. For the TDDFT calculations using the LFA-corrected PBE functionals, both the PBE XC kernel and the LFA exchange kernel should be adopted for a consistent approximation on Exc[ρ]. However, in this work, we only adopt the PBE XC kernel, and neglect the LFA exchange kernel for computational simplicity. Note that similar tricks have been constantly used in the AC model potential approach (e.g., the XC kernel of a local or semilocal functional is adopted) [16,30,33,42,72–74]. For example, the LDA XC kernel is frequently adopted for the TDDFT calculations using the LB94 potential [42]. For finite systems, this approximation should not make much difference in the prediction of valence and Rydberg excitation energies. As shown in TableII, all the LFA-corrected PBE functionals and LB94 perform well for both the valence and Rydberg excitations, while PBE severely
TABLE II. Mean absolute errors (in eV) of the 19 valence and 23 Rydberg excitation energies of five molecules [71]. The LB94 results are taken from Ref. [42].
System PBE LFA-PBE RILFA-PBE LFAs-PBE LB94
Valence (19) 0.32 0.29 0.29 0.29 0.36
Rydberg (23) 1.30 0.46 0.46 0.49 0.73
underestimates Rydberg excitation energies due to its incorrect asymptote (see Supplemental Material [46]).
V. CONCLUSIONS
In conclusion, we have developed the LFA scheme, wherein an exchange density functional whose functional derivative has the correct (−1/r) asymptote can be directly added to any semilocal density functional. In contrast to semilocal approximations, the LFA exchange kernel in reciprocal space exhibits the desirable singularity of the type O(−1/q2), which is an important feature for the description of excitonic effects in nonmetallic solids. Applying the LFA scheme to PBE, the resulting LFA-PBE (ω = 0.15 bohr−1) has yielded accurate IPs and Rydberg excitation energies for a wide range of atoms and molecules, while performing similarly to PBE for various properties that are insensitive to the asymptote. Without loss of accuracy, two closely related schemes (RILFA and LFAs) have been developed for the efficient treatment of large systems.
Relative to the popular model XC potential scheme, LFA-PBE is significantly superior for ground-state energies and related properties. It remains to be seen if the LFA-corrected PBE
functionals will perform well for properties sensitive to the details of the XC potential (not just to the asymptote), such as quantum defects [75].
As with all pure density-functional methods (e.g., semilo- cal functionals and model XC potentials), some limitations remain. Due to the lack of HF exchange, the LFA-corrected PBE functionals may suffer from the SIE problems, energy- gap problems, and charge-transfer problems (e.g., see the discussions in Ref. [42]). Nevertheless, the energy-gap prob- lems may be circumvented by the perturbation approach recently developed in Ref. [76]. Although the LC hybrid scheme, which has remedied several qualitative failures of pure density-functional methods, could be reliably accurate for a very wide range of applications [42], it can be impractical for very large systems due to the expensive computational cost. By contrast, our LFAs-PBE, which has the correct (−1/r) asymptote with essentially no added computational cost relative to PBE, is potentially very useful for the study of the ground-state energies and related properties, frontier orbital energies, valence and Rydberg excitation energies, and time-dependent properties of very large systems.
ACKNOWLEDGMENTS
This work was supported by the National Science Coun- cil of Taiwan (Grant No. NSC101-2112-M-002-017-MY3), National Taiwan University (Grants No. 99R70304, No.
101R891401, and No. 101R891403), and the National Center for Theoretical Sciences of Taiwan. We thank Dr. Yihan Shao (Q-Chem, Inc.) for helpful discussions.
[1] P. Hohenberg and W. Kohn,Phys. Rev.136, B864 (1964).
[2] W. Kohn and L. J. Sham,Phys. Rev.140, A1133 (1965).
[3] E. Runge and E. K. U. Gross,Phys. Rev. Lett.52, 997 (1984).
[4] M. Petersilka, U. J. Gossmann, and E. K. U. Gross,Phys. Rev.
Lett.76, 1212 (1996).
[5] C. A. Ullrich, Time-Dependent Density-Functional Theory:
Concepts and Applications (Oxford University, New York, 2012).
[6] S. K¨ummel and L. Kronik,Rev. Mod. Phys.80, 3 (2008).
[7] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, L. A. Constantin, and J. Sun,Phys. Rev. Lett.103, 026403 (2009).
[8] A. J. Cohen, P. Mori-S´anchez, and W. Yang,Science321, 792 (2008).
[9] A. J. Cohen, P. Mori-S´anchez, and W. Yang,Chem. Rev.112, 289 (2011).
[10] J.-D. Chai,J. Chem. Phys.136, 154104 (2012).
[11] M. Levy, J. P. Perdew, and V. Sahni,Phys. Rev. A 30, 2745 (1984).
[12] C.-O. Almbladh and U. von Barth,Phys. Rev. B31, 3231 (1985).
[13] A. G¨orling,Phys. Rev. Lett.83, 5459 (1999).
[14] P. W. Ayers and M. Levy,J. Chem. Phys.115, 4438 (2001).
[15] J. P. Perdew and A. Zunger,Phys. Rev. B23, 5048 (1981).
[16] M. E. Casida, C. Jamorski, K. C. Casida, and D. R. Salahub, J. Chem. Phys.108, 4439 (1998).
[17] D. J. Tozer and N. C. Handy,J. Chem. Phys.109, 10180 (1998).
[18] A. Savin, inRecent Developments and Applications of Modern Density Functional Theory, edited by J. M. Seminario (Elsevier, Amsterdam, 1996), p. 327.
[19] H. Iikura, T. Tsuneda, T. Yanai, and K. Hirao,J. Chem. Phys.
115, 3540 (2001).
[20] T. Yanai, D. P. Tew, and N. C. Handy,Chem. Phys. Lett.393, 51 (2004).
[21] O. A. Vydrov, J. Heyd, A. V. Krukau, and G. E. Scuseria, J. Chem. Phys.125, 074106 (2006).
[22] E. Livshits and R. Baer, Phys. Chem. Chem. Phys.9, 2932 (2007).
[23] J.-D. Chai and M. Head-Gordon,J. Chem. Phys.128, 084106 (2008).
[24] J.-D. Chai and M. Head-Gordon,Phys. Chem. Chem. Phys.10, 6615 (2008).
[25] J.-D. Chai and M. Head-Gordon,Chem. Phys. Lett.467, 176 (2008).
[26] J.-D. Chai and M. Head-Gordon,J. Chem. Phys.131, 174105 (2009).
[27] Y.-S. Lin, C.-W. Tsai, G.-D. Li, and J.-D. Chai,J. Chem. Phys.
136, 154109 (2012).
[28] Y.-S. Lin, G.-D. Li, S.-P. Mao, and J.-D. Chai,J. Chem. Theory Comput.9, 263 (2013).
[29] R. van Leeuwen and E. J. Baerends,Phys. Rev. A 49, 2421 (1994).
[30] D. J. Tozer and N. C. Handy, J. Chem. Phys. 109, 10180 (1998).
[31] D. J. Tozer,J. Chem. Phys.112, 3507 (2000).
[32] P. R. T. Schipper, O. V. Gritsenko, S. J. A. van Gisbergen, and E. J. Baerends,J. Chem. Phys.112, 1344 (2000).
[33] M. E. Casida and D. R. Salahub, J. Chem. Phys.113, 8918 (2000).
[34] M. Gr¨uning, O. V. Gritsenko, S. J. A. van Gisbergen, and E. J.
Baerends,J. Chem. Phys.114, 652 (2001).
[35] Q. Wu, P. W. Ayers, and W. T. Yang,J. Chem. Phys.119, 2978 (2003).
[36] A. D. Becke and E. R. Johnson,J. Chem. Phys.124, 221101 (2006).
[37] F. Tran and P. Blaha,Phys. Rev. Lett.102, 226401 (2009).
[38] X. Andrade and A. Aspuru-Guzik,Phys. Rev. Lett.107, 183002 (2011).
[39] A. P. Gaiduk, D. S. Firaha, and V. N. Staroverov, Phys. Rev.
Lett.108, 253005 (2012).
[40] A. P. Gaiduk and V. N. Staroverov,J. Chem. Phys.131, 044107 (2009).
[41] A. P. Gaiduk and V. N. Staroverov, Phys. Rev. A83, 012509 (2011).
[42] C.-W. Tsai, Y.-C. Su, G.-D. Li, and J.-D. Chai, Phys. Chem.
Chem. Phys.15, 8352 (2013).
[43] E. Fermi and E. Amaldi, Mem. R. Accad. Ital.6, 117 (1934).
[44] P. W. Ayers, R. C. Morrison, and R. G. Parr,Mol. Phys.103, 2061 (2005).
[45] J. P. Perdew,Adv. Quantum Chem.21, 113 (1990).
[46] See Supplemental Material athttp://link.aps.org/supplemental/
10.1103/PhysRevA.87.052510for additional information.
[47] R. F. W. Bader, Atoms in Molecules: A Quantum Theory (Clarendon, Oxford, 1990).
[48] F. L. Hirshfeld,Theor. Chim. Acta44, 129 (1977).
[49] N. H. March,Phys. Rev. A36, 5077 (1987).
[50] Ph. Ghosez, X. Gonze, and R. W. Godby,Phys. Rev. B56, 12811 (1997).
[51] L. Reining, V. Olevano, A. Rubio, and G. Onida,Phys. Rev. Lett.
88, 066404 (2002).
[52] Y.-H. Kim and A. G¨orling,Phys. Rev. Lett.89, 096402 (2002).
[53] V. U. Nazarov and G. Vignale,Phys. Rev. Lett.107, 216402 (2011).
[54] J. P. Perdew, K. Burke, and M. Ernzerhof,Phys. Rev. Lett.77, 3865 (1996).
[55] R. A. Kendall and H. A. Fr¨uchtl,Theor. Chem. Acc.97, 158 (1997).
[56] J. E. Jaffe and A. C. Hess,J. Chem. Phys.105, 10983 (1996).
[57] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Jr.,Phys.
Rev. Lett.49, 1691 (1982).
[58] J. P. Perdew and M. Levy,Phys. Rev. B56, 16021 (1997).
[59] M. E. Casida,Phys. Rev. B59, 4694 (1999).
[60] Y.-S. Lin, C.-W. Tsai, G.-D. Li, and J.-D. Chai,J. Chem. Phys.
136, 154109 (2012).
[61] Y. Shaoet al.,Phys. Chem. Chem. Phys.8, 3172 (2006).
[62] P. A. M. Dirac,Proc. Cambridge Philos. Soc.26, 376 (1930).
[63] J. P. Perdew and Y. Wang,Phys. Rev. B45, 13244 (1992).
[64] L. A. Curtiss, K. Raghavachari, P. C. Redfern, and J. A. Pople, J. Chem. Phys.106, 1063 (1997).
[65] L. A. Curtiss, P. C. Redfern, K. Raghavachari, and J. A. Pople, J. Chem. Phys.109, 42 (1998).
[66] L. A. Curtiss, K. Raghavachari, P. C. Redfern, and J. A. Pople, J. Chem. Phys.112, 7374 (2000).
[67] Y. Zhao, B. J. Lynch, and D. G. Truhlar,J. Phys. Chem. A108, 2715 (2004).
[68] Y. Zhao, N. Gonz´alez-Garc´ıa, and D. G. Truhlar,J. Phys. Chem.
A109, 2012 (2005); 110, 4942(E) (2006).
[69] P. Jureˇcka, J. ˇSponer, J. ˇCern´y, and P. Hobza,Phys. Chem. Chem.
Phys.8, 1985 (2006).
[70] M. Levy and J. P. Perdew,Phys. Rev. A32, 2010 (1985).
[71] S. Hirata and M. Head-Gordon,Chem. Phys. Lett. 314, 291 (1999).
[72] S. J. A. van Gisbergen, J. G. Snijders, and E. J. Baerends,Phys.
Rev. Lett.78, 3097 (1997).
[73] I. Vasiliev and R. M. Martin,Phys. Rev. A69, 052508 (2004).
[74] I. Vasiliev and J. R. Chelikowsky,Phys. Rev. A 82, 012502 (2010).
[75] M. van Faassen and K. Burke, J. Chem. Phys. 124, 094102 (2006).
[76] J.-D. Chai and P.-T. Chen,Phys. Rev. Lett.110, 033002 (2013).
