The shell model for the exchange-correlation hole in the strong-correlation limit Hilke Bahmann, Yongxi Zhou, and Matthias Ernzerhof

Citation: The Journal of Chemical Physics 145, 124104 (2016); doi: 10.1063/1.4962738 View online: http://dx.doi.org/10.1063/1.4962738 View Table of Contents: http://aip.scitation.org/toc/jcp/145/12 Published by the American Institute of Physics

Articles you may be interested in Construction of exchange-correlation functionals through interpolation between the non-interacting and the strong-correlation limit The Journal of Chemical Physics 143, 124103 (2015); 10.1063/1.4931160 Design of exchange-correlation functionals through the correlation factor approach The Journal of Chemical Physics 143, 144102 (2015); 10.1063/1.4932074 Perspective: Kohn-Sham density functional theory descending a staircase The Journal of Chemical Physics 145, 130901 (2016); 10.1063/1.4963168 Long-range interactions from the many-pair expansion: A different avenue to dispersion in DFT The Journal of Chemical Physics 146, 024111 (2017); 10.1063/1.4973728 The impact of long-range electron-hole interaction on the charge separation yield of molecular photocells The Journal of Chemical Physics 146, 034103 (2017); 10.1063/1.4973984 The Perdew–Burke–Ernzerhof exchange-correlation functional applied to the G2-1 test set using a planewave basis set The Journal of Chemical Physics 122, 234102 (2005); 10.1063/1.1926272

THE JOURNAL OF CHEMICAL PHYSICS 145, 124104 (2016)

The shell model for the exchange-correlation hole in the strong-correlation limit Hilke Bahmann,1,a) Yongxi Zhou,2 and Matthias Ernzerhof2,a) 1 2

Department of Chemistry, Technische Universität Berlin, Strasse des 17 Juni 135, 10623 Berlin, Germany Département de Chimie, Université de Montréal, C.P. 6128 Succursale A, Montréal, Québec H3C 3J7, Canada

(Received 23 June 2016; accepted 1 September 2016; published online 23 September 2016) We present a model for the exchange-correlation hole and the exchange-correlation energy in the strong-correlation (SC) limit of density functional theory. The SC limit is useful in the construction of exchange-correlation functionals through interpolation of the adiabatic connection. The new approximation (referred to as shell model) is an improvement of the non-local radius (NLR) model recently proposed by Wagner and Gori-Giorgi [Phys. Rev. A 90, 052512 (2014)]. The NLR model does not correctly reproduce the limit of the strongly correlated homogeneous electron gas and this shortcoming is remedied by the shell As in the case of the NLR model, the spherically  model. u ρ(r + u) is the starting point for the construction of the averaged electron density ρ(r,u) = dΩ 4π shell model and it is also its computational bottleneck. We show how ρ(r,u), the NLR, and the shell model can be implemented efficiently. For this purpose, analytical integrals for the normalization and the energy density of the underlying holes are provided. Employing the shell model, we illustrate how improved adiabatic connection interpolations can be constructed. Published by AIP Publishing. [http://dx.doi.org/10.1063/1.4962738]

I. INTRODUCTION

Over the past two decades, Kohn-Sham (KS) density functional theory (DFT)1–5 has developed into the most widely used electronic structure theory method. This is the result of the advances made in the construction of approximate exchange-correlation (E XC) functionals that evolved into very sophisticated and complicated expressions which proved to be accurate in numerous applications. Various of the popular E XC functionals rely on physical constraints, others draw on numerical fitting to data sets, and most combine both, first-principles and empirical elements. Despite the remarkable successes of modern exchangecorrelation functionals, there are several notorious problem areas that are not accurately described by the popular functionals. One such problem area is characterized through strong correlation. In the strong-correlation limit (SCL), the Coulomb repulsion among the electrons dominates their kinetic energy. This is achieved by scaling the repulsion among the electrons to infinity, while keeping the electron density fixed at its physical value. The exchange-correlation functional in the SCL is a focus of interest in DFT and exact solutions and models for the SCL are available.6–11 However, in actual applications the exact solutions are computationally expensive and thus efficient and simple approximations are required. Various models have been proposed in the past, starting from the point charge plus continuum model of Seidl, Perdew, and Levy.6 Recently Wagner and Gori-Giorgi developed the non-local radius (NLR) model11 to the exchange-correlation energy which is based on a model to the exchange- correlation a)Authors to whom correspondence should be addressed. Electronic

addresses: [email protected] UMontreal.ca.

and

0021-9606/2016/145(12)/124104/10/$30.00

Matthias.Ernzerhof@

hole in the SCL. Assuming an electron located at a reference point r, the electron density is removed completely around r up to a given radius set by the normalization condition to the exchange-correlation hole.2 There are various problems where the interaction among the electrons resembles the strong-correlated limit— an example being the stretched H2 molecule. Another area where strong-correlation is important is the description of open-shell atoms12 that are characterized by degeneracies of the ground-states. Apart from these two examples, the main motivation for the study of the SCL is provided by the adiabatic connection (AC).13 The AC involves a couplingconstant (λ) decomposition of E XC and of the underlying exchange-correlation hole. It connects the noninteracting, exchange-only limit (λ = 0), and the strong correlation limit (λ = ∞). The physical E XC is obtained by integration over the coupling constant from zero to one. Numerous hybrid schemes are based on interpolations of the AC. Some interpolations (e.g., Refs. 6–11) use the exchange-correlation energy in the SCL as well. The size-consistency problem as well as other difficulties of some global hybrid schemes are circumvented14 by local interpolations along the AC where energy densities are mixed at each point in space.15,16 Another approach that also relies on the AC has recently been developed by Becke.12 In the resulting functional for static, dynamical, and strong correlation, a model for the exchange-correlation energy at λ = 1 is proposed which is then extended to an exchangecorrelation functional through an implicit interpolation of the adiabatic connection. In Section II, we provide relevant notions and concepts employed for the development of our model. The shell model and its construction are detailed in Section III. In Section IV, we illustrate the shell model through applications to atoms

145, 124104-1

Published by AIP Publishing.

124104-2

Bahmann, Zhou, and Ernzerhof

J. Chem. Phys. 145, 124104 (2016)

and molecules. We employ it in constructions of the adiabatic connection and provide evidence that the shell model describes the SCL more accurately than the NLR model. In Section V we provide details of the efficient implementation of both the shell and the NLR model.

II. PRELIMINARIES

To prepare the derivation of the shell model and its combination with exact exchange in a local hybrid functional, the underlying concepts of the NLR model and the adiabatic connection are briefly reviewed. The SCL can be analyzed through the key quantity in the description of exchange-correlation effects, the exchangecorrelation hole ρXC(r, u). It is commonly expressed as a function of two coordinates representing the position of the reference electron r and the distance u = r′ − r to a second electron at r′. ρXC(r, u) describes the reduction in the electron density at position r′ = r + u due to the presence of the reference electron at point r. Semi-local models to ρXC(r, u) that are used to derive functionals within the generalizedgradient approximation (GGA) are based on the electron density ρ(r) and its gradient ∇ρ(r) at the reference point. In the SCL ρXC(r, u) can be quite non-local and strongly dependent on the non-local electron density ρ(r + u). By the non-local electron density we refer to the electron density as a function of u, ρ(r) → ρ(r + u).

(1)

The sum of ρXC(r, u) and ρ(r + u), i.e., ρXC(r, u) + ρ(r + u) = ρint(r, u), yields the interacting charge density seen at r + u by the reference electron located at r. Consequently the total electron-electron interaction, as a sum of the Hartree (J) and E XC contribution, is given in terms of the density and hole density by   ρ(r)ρint(r, u) 3 3 1 d rd u J + E XC = 2 u ( )   1 ρ(r)ρ(r + u) ρ(r)ρXC(r, u) 3 3 = + d r d u. 2 u u (2) Since the Coulomb operator depends only on the absolute value of the electron-electron distance u = |u|, virtually all models for the exchange-correlation hole are functions of u only, effectively eliminating two variables through an integration over the orientation of u (i.e., Ωu). This motivates the introduction of the spherically averaged exchangecorrelation hole and spherically averaged electron density  dΩu ρXC(r,u) = ρXC(r, r + u) (3) 4π and  dΩu ρ(r,u) = ρ(r + u), (4) 4π respectively. A problem that we characterize as local is determined by the density for vanishing u, whereas a nonlocal problem depends on the electron density for larger interelectronic distances. The notion of small and large values

of u is determined on the basis of the non-local Wigner-Seitz radius uc introduced within the NLR model,11  uc 1 = 4π du u2 ρ(r,u). (5) 0

It defines the radius of a sphere around the reference electron that contains exactly one electron. In the homogeneous limit, uc reduces to the well-known local Wigner-Seitz radius that characterizes the HEG in terms of the density, (

3 r S (r) = 4π ρ(r)

) 13

.

(6)

For non-homogeneous densities, uc is a function of the position r. If the density varies on the scale of uc , we refer to it as a problem of non-local character, otherwise we refer to the problem as being local. In particular for systems containing only a few electrons, uc can extend over a significant portion of the system leading to highly non-local exchange-correlation holes. For one-electron systems, ρXC(r,u) is highly non-local since uc becomes infinite and then ρ(r,u) varies of course strongly on the scale of uc . In the present work, ρ(r,u) is the starting point for the approximations to exchange-correlation holes in the SCL. Therefore, the functionals discussed below are best qualified as completely non-local. They are closely related to the branch of approximate exchange-correlation functionals referred to as weighted density approximations (WDA).17–22 To elaborate on the NLR model, we write its ρXC(r,u) in terms of an exchange-correlation factor ( f XC(r,u))23–25 that is applied to the spherically averaged density, NLR ρNLR XC (r,u) = f XC (r,u)ρ(r,u),

(7)

where f NLR XC (r,u) is given by f NLR XC (r,u) = −θ(uc (r) − u).

(8)

The Heaviside step function θ(uc (r) − u) drops from one to zero if u exceeds the non-local radius uc (r). ρNLR XC (r,u) is then employed to calculate the exchange-correlation energy according to the conventional expression,   NLR 3 E XC = d r ρ(r)2π du u ρNLR (9) XC (r,u). This functional has been analyzed in detail by Wagner and Gori-Giorgi, and it has various desirable properties: it is intuitive and accessible through efficient algorithms that have been discussed before16 and are further advanced in this work. However, the NLR model has the disadvantage that it does not reproduce the SCL of the homogeneous electron gas. The reason for this error is that the hole in Eq. (7) has negative values only while the exact exchangecorrelation hole in the SCL, in general, has positive values as well. We address this problem below and provide an improved approximation that correctly reproduces the HEG limit. As mentioned in the introduction, the SCL of the exchange-correlation energy or approximations thereof can be combined with existing approximations to E XC through hybridization by invoking the AC. In detail, the AC involves a coupling-constant decomposition of E XC and of the underlying

124104-3

Bahmann, Zhou, and Ernzerhof

J. Chem. Phys. 145, 124104 (2016)

exchange-correlation hole,  1  1   3 λ E XC = dλW XC(λ) = dλ d r ρ(r)2π du u ρXC (r,u), 0

0

(10) where W XC(λ) is given by W XC(λ) = ⟨Ψ(λ)|Vee |Ψ(λ)⟩ − J.

(11) λ ρXC (r,u)

The λ-dependent exchange-correlation hole is defined2 through the pair density of Ψ(λ). Ψ(λ) is the coupling-constant-dependent ground-state wave function of the Hamiltonian, H(λ) = T + v(λ) + λVee .

(12)

The λ-dependent external potential v(λ) ensures that the ground-state density of H(λ) is equal to the physical density of the system obtained for λ = 1. For λ = 1, the external potential v(λ) reduces to the true external potential generated by the nuclei. Approximate exchange-correlation functionals can be obtained by performing the λ-average in Eq. (10) using interpolation schemes. The interpolations typically employ λ various values of W XC(λ) or ρXC (r,u) along the adiabatic connection (λ-discretization). For instance, in the λ = 0 λ=0 or noninteracting case, both W XC(λ = 0) and ρXC (r,u) = ρX(r,u) are readily accessible. Similarly, fairly accurate λ=1 approximations to W XC(λ = 1) and ρXC (r,u) are avail12,26,27 able which are then used as another interpolation point along the adiabatic connection. The resulting interpolations give rise to various hybrid schemes (see, e.g., Refs. 5, 12, 26, 28, and 29). Additionally considering an approximate slope of W XC(λ) at λ = 0 results in double hybrid functionals.30–33 Various interpolation schemes use the exact or approximate 6–11,16 W XC(λ = ∞) or ρ∞ In the case where XC(r,u) as well. λ (r,u) is approximated for different λ values instead of ρXC W XC(λ), the corresponding exchange-correlation energies per particle are available,  λ w XC(λ, r) = 2π du u ρXC (r,u), (13) and r-dependent local AC interpolations can be formulated. The coupling constant integral of w XC(λ, r) yields the exchange-correlation energy per particle,  1 ϵ XC(r) = dλ w XC(λ, r), (14) 0

such that

 E XC =

d r ρ(r)ϵ XC(r). 3

(15)

Using accurate ingredients for small and infinite coupling strengths (w XC(λ = 0, r) and w XC(λ = ∞, r)), several interpolation schemes have been compared for the construction of exchange-correlation functionals in a recent study.15 To illustrate the performance of the shell model, we use it within the ACLSD approximation16 to calculate exchangecorrelation energies of atoms and atomization energies of molecules. The ACLSD approximation represents a particular simple example of an AC interpolation scheme that results in

the expression for ϵ XC(r,u),
LSD ϵ ACLSD (r) = ϵ SC (r) ϵ X(r) − ϵ SC XC XC(r) + a XC(r) .

(16)

In this formula, ϵ X(r) is the exact-exchange energy density LSD and ϵ SC (r) depends XC(r) = w XC(λ = ∞, r). The coefficient a on r through the local spin density. It is determined by replacing all energies per particle in Eq. (16), on the left and on the right hand side, by their counterparts for the HEG. This results in an equation for aLSD(r). In actual calculations, an approximation to the LSD correlation energy is required. Here we employ the Perdew-Wang parametrization34 for all the LSD calculations. Note that we do not employ the SCL obtained from a given parametrization (in our case the Perdew-Wang parametrization). Instead we use the available estimates of the SCL ϵ SC XC of the HEG that will be discussed in more detail in Section III. For the HEG, the mixing parameter should go to zero in the SCL reducing the ACLSD functional to strong correlation only. With a parametrization for the LSD correlation energy that does not exhibit the correct SCL, small absolute values of aLSD might occur in the SCL. For finite systems, this happens only in low-density regions where the correlation energy density is small anyway. While more sophisticated local mixing coefficients a(r) can be constructed (see, e.g., Ref. 15), ACLSD is suitable to illustrate the shell model. A useful property of the ACLSD approximation is that it is exact in the one-electron limit: In this case, ϵ X(r) and ϵ SC XC(r) become identical and both exactly cancel the spurious self-interaction. It is therefore irrelevant what function a(r) is used. Another important property of the ACLSD model is that it correctly reproduces the high density limit in which the exchange-correlation energy is dominated by exact exchange. Similarly, in the low density limit, the ACLSD correctly reduces to ϵ SC XC(r). By construction, ACLSD yields the right answer if applied to the HEG. Now in practical applications, ϵ SC XC(r) has to be approximated and in case that we employ the NLR model, we do not obtain the correct result in the SCL (or equivalently the low density limit) of the HEG. However, this limit is one of the essential columns in the construction of non-empirical functionals and we restore it in the shell model.

III. THE SHELL MODEL FOR THE EXCHANGE-CORRELATION HOLE

For uniform densities, uc of the NLR model reduces to the Wigner-Seitz radius (uc → r S ) and the corresponding energy density becomes 0.75 , (17) rS where r S depends on the density through Eq. (6). Recently a renewed discussion emerged35,36 about the correct strongcorrelation limit of the HEG. In the low density limit, which is a strong-correlation limit, the exchange-correlation energy of the HEG scales like a pure exchange energy, ϵ NLR–HEG (ρ) = − XC

C , (18) rS where the value of C = 0.895 93 has been obtained from a Wigner-crystal-based calculation.8 This value has been put ϵ SC–HEG (ρ) = − XC

124104-4

Bahmann, Zhou, and Ernzerhof

J. Chem. Phys. 145, 124104 (2016)

into question by the work of Lewin and Lieb.35 More recently, Seidl, Vuckovic, and Gori-Giorgi36 provided a new estimate for the upper bound of ϵ SC XC. The bound obtained by these authors indicates that C ≥ 0.875 90. According to Ref. 36 C = 0.875 90 provides the best estimate of the SCL of the HEG. The correct value of C seems to be a matter of ongoing investigations; therefore, we employ the original value of C = 0.895 93 to generate the tabulated data. For comparison, we also use the lower estimate of 0.875 90 and report the corresponding MAE values. In any case the energy density of the NLR is substantially in error and above the exact strongcorrelation limit of the HEG. This shortcoming originates from the fact that the NLR exchange-correlation hole has only negative values while the exchange-correlation hole in the SCL in general can have positive as well as negative values. This means that the range of NLR hole is reduced to a minimum. In fact the NLR model provides the most short-ranged hole possible satisfying the exact condition, ρXC(r,u) ≥ −ρ(r,u).

(19)

If the hole crosses to positive values, the negative area has to be extended to comply with the normalization of ρXC(r,u),  4π

du u2 ρXC(r,u) = −1.

(20)

As a consequence, the range of the exchange-correlation hole is enlarged compared to the NLR model, allowing for lower exchange-correlation energy densities. This reasoning leads to the introduction of the shell model for the exchangecorrelation hole in which in addition to a contribution equal to −ρ(r,u), there is now also one equal to ρ(r,u). Setting the value of the second shell to ρ(r,u) is a choice and not dictated by a constraint or condition. A corresponding scaling parameter can in principle be used to accommodate other exact conditions or it can be adjusted to empirical data. In Fig. 1 the shell model is compared to the NLR model for the Be atom. The interacting charge density ρint(r,u) = ρ(r,u) + ρXC(r,u) illustrates the modification that the hole makes to the electron density. For the NLR model ρint(r,u) reduces to ρ(r,u) beyond uc . The increased range of the our new model can be seen as well as the additional shell given by ρ(r,u) for u ∈ [u s ,uc ]. In the limit of uniform densities, the NLR and shell-model holes are composed of one or two step functions, respectively.

FIG. 2. Co-motion function for a spin-unpolarized two-electron system defined through an electron density given by ρ(r) = Nα e −αr with α = 3. For this system, the co-motion function f2(r) can be expressed (f2(r) = − f (r )r/r ) through f (r ), a scalar function of r . If the electron at the reference position r approaches the nucleus, the second electron at f (r ) moves out to infinity.

To further motivate the shell model for ρSC XC(r,u), we discuss the exact pair density PSC(r1, r2) in the SCL, N   ρ(s) δ(r1 − fi (s))δ(r2 − f j (s)). (21) PSC(r1, r2) = d3 s N i, j=1 i, j

In this expression, fi (r) represents a co-motion function of electron i, and its determination is described in Ref. 9. In Fig. 2 we provide a simple illustration of the co-motion function for a spherically symmetric two-electron system defined through an exponential density.7 In this case, f1(r) = r, and f2(r) = − f (r)r/r, where f (r) is plotted as a function of r. Using this expression for PSC(r1, r2), we focus on two-electron problems and calculate the exchange-correlation hole in the SCL, ρSC XC(r,u) = −ρ(r,u) + δ(u − ( f (r) + r)).

(22)

To obtain this expression, we used the fact that the electrons are located on a line passing through the nucleus. The exchangecorrelation hole in Eq. (22) clearly shows that the hole in the SCL has negative as well as δ-function-type positive parts. The δ function contributions are a characteristic of the strongcorrelation limit, exclusively. For finite coupling strength, the electrons are not localized in δ-function-type orbits. Our objective is thus to develop a model hole that accounts for the oscillatory behavior of the exchange-correlation hole in the SCL, describing the fact that it assumes positive as well as negative values, without using the extreme electron localization. To this end, we introduce the simple shell model

FIG. 1. Illustration of the non-local radius (left) and the shell model (right) for the Be atom. The blue curves show the electron density and the red dashed ones the interacting charge density ρ int(r, u) = ρ(r, u) + ρ XC(r, u) as a function of the interelectronic distance u with r = 0.42a 0.

124104-5

Bahmann, Zhou, and Ernzerhof

J. Chem. Phys. 145, 124104 (2016)

(Fig. 1) in which, in addition to a shell of negative density around the reference electron, a positive shell of electron density is added at long range. In this model, we represent the N − 1 electrons that interact with the reference electron simply through the electron density instead of representing it through N − 1 δ functions. This is of course an approximation and it appears to be appropriate in view of the final objective of modelling the exchange-correlation hole for values of the coupling constant λ between 0 and 1. ρSC XC(r,u) of Eq. (22), and more generally, the pair density of Eq. (21) implies that in the strong-correlation limit, the exchange-correlation hole always has positive contributions. This result seems to be in contradiction to the known exact exchange-correlation hole of the infinitely stretched H2 molecule37 which has negative contributions only. However, in the strong-correlation limit, the coupling constant λ is taken to infinity and in the stretched H2 the bond length is increased to infinity at finite coupling constant. In general there is no reason to assume that the same physical solution should result since both systems have different Hamiltonians. The two approaches nonetheless yield the same energy for the infinitely stretched H2 since at increasing separation between the hydrogen atoms, in the strong-correlation limit, the hole contribution of the form of a δ function moves out to an infinite distance from the reference electron and becomes irrelevant for the exchange-correlation energy. Interestingly, the same happens within the shell model, upon stretching H2 the positive shell moves out to infinity and does not contribute to energetics. So while the stretched H2 and the stretched H2 in the strong-correlation limit are not identical, stretching H2 to infinity in the strong-correlated limit yields the E XC obtained by just stretching H2. In the sense described, the strong-correlation limit is capable of dealing with this problem of stretched H2 as shown in detail in Ref. 38. The mathematical form of the shell model is obtained by multiplying the spherically averaged density with an exchange-correlation factor consisting of two step functions, shell ρshell (r,u) XC (r,u) = ρ(r,u) f

= ρ(r,u)θ(uc (r) − u)sign(u − u s (r)),

(23)

where the two parameters denote the radius of the hole (uc ) and the zero crossing (u s ) respectively. To eliminate u s and to ensure that the HEG limit is recovered, we impose the condition that if the shell model is applied to uniform densities, the corresponding strong-correlation limit of the exchangecorrelation energy per particle (ϵ SC–HEG (ρ)) is obtained XC  ϵ shell u ρshell XC (ρ(r)) = 2π XC (r,u) du  (24) HEG 2π ρ u θ(uc − u)sign(u − u s ) du

determined for any system, homogeneous or inhomogeneous, through the normalization condition (Eq. (19)). The resulting values for u s and uc are subsequently employed to calculate the exchange correlation energy per particle ϵ shell XC (r) of the shell model,  ϵ shell (r) = 2π du u ρshell (26) XC XC (r,u).

IV. APPLICATIONS OF THE SHELL MODEL

All energies reported here are obtained in a post-SCF manner with the Turbomole program39,40 package using PBE41 densities and 6-311g** basis sets.42,43 The exception being the results reported in Table I where CCSD(T) densities and aug-cc-pVQZ basis sets44,45 have been employed. The exactexchange energy density per particle used in the ACLSD functional (cf. Eq. (16)) is defined in terms of the molecular orbitals {φi (r)} as  occ φi (r′)φ j (r′) 3 ′ 1  φ (r)φ (r) dr. (27) ϵ exact (r) = − i j X 2ρ(r) i j |r − r′| It is calculated analytically at each grid point in real space as explained in the context of the semi-numerical implementation of local hybrid functionals.46 In Ref. 11 the exact exchange-correlation energy in the SCL E SC XC has been calculated for various atoms. To these values we add the ones obtained with the shell model (Table I). The shell model consistently yields lower XC energies compared to NLR. Furthermore, the results obtained with the shell model are a considerable improvement over NLR. If we use the more recently calculated (and slightly less negative value) for the SCL of the HEG in Eq. (18) the exchange-correlation energies increase and the MAE becomes 0.246 hartree. The energies obtained with the shell model seem to be quite sensitive to the underlying value of the SC-HEG exchange-correlation energy. Note that we employ CCSD(T) densities, whereas the results of Ref. 11 were calculated with other highly accurate densities. Since we have the non-local radius of the NLR model at our disposition, we employ it to replace the local r S in the SCL of the local spin density approximations (LSDA). In the homogeneous limit, both radii coincide, but they are different, in general, for inhomogeneous systems. As can be seen from Table I, the LSDA results significantly underestimate E SC XC TABLE I. Exchange-correlation energies in hartree obtained from the nonlocal radius model (NLR) and the shell model (SM) using CCSD(T)/aug-ccpVQZ densities.

=

−0.895 93 . rS Based on this equation we find that the ratio between u s and uc is independent of the density and given by us = 0.849 488. (25) uc For the revised SCL of the HEG, this ratio would be 0.912 17. This equation eliminates one unknown (u s ) and uc is then = ϵ SC–HEG (ρ) = XC

H− He Li Be Ne MAE a Exact

SCEa

NLR

SM

SC-LSDA(r S )

SC-LSDA(u c )

−0.569 −1.498 −2.596 −4.021 −19.99

−0.556 −1.427 −2.496 −3.819 −18.29

−0.566 −1.455 −2.553 −3.924 −19.91

−0.674 −1.729 −2.969 −4.522 −21.57

−0.718 −2.398 −4.378 −6.404 −26.87

0.42

0.053

values for the SCL taken from Ref. 11.

0.56

2.42

124104-6

Bahmann, Zhou, and Ernzerhof

J. Chem. Phys. 145, 124104 (2016) TABLE II. Atomization energies in kcal/mol obtained with ACLSD16 using the non-local radius model (NLR) and the shell model (SM). Mean absolute and mean signed errors (MAE and MSE, respectively) with respect to experimental energies (zero-point energy corrections subtracted as in Ref. 41) are shown for the subsets as well. ACLSD

FIG. 3. r S and u c in the Be atom as a function of the distance to the nucleus. The underlying density has been obtained through CCSD(T) calculations employing the aug-cc-pVQZ basis set.

which is even more pronounced when the Wigner-Seitz radius is replaced by the non-local radius uc . The latter is overall smaller than r S as confirmed by the plot in Fig. 3 where r S and uc are compared for the Be atom. Close to the nucleus, however, uc is larger than r S . While r S exhibits some shell structure, uc varies much less. At large distances from the nucleus, uc becomes a linear function of r. This is not surprising since uc measures the radius of a sphere occupied by one electron and this radius starts to coincide with the distance from the nucleus if the latter becomes large. Now we turn to the use of the shell model in AC interpolations. In Ref. 16 we put forward the ACLSD approach that is reviewed in Sec. II. We found that ACLSD can compete with PBE41 and PBE hybrid28,30,47,48 for total energies of atoms while atomization energies are significantly underestimated. The MAEs for PBE and PBE hybrid are 2.19 kcal/mol and 4.62 kcal/mol (single-bonded molecules) and 10.74 kcal/mol and 3.63 kcal/mol (congested systems, characterized by multiple overlapping occupied orbitals in the bond region). The resulting total MAEs are 6.69 kcal/mol (PBE) and 4.09 kcal/mol (PBE hybrid). We concluded that an additional measure for non-dynamical correlation is required for improved atomization energies. It is, however, also desirable to use a more accurate model for the energy density in the SCL. In particular, it should be useful to have a strongcorrelation model that recovers the homogeneous limit. This is so because the ACLSD model relies on this limit to determine the real-space function aLSD(r) (cf. Eq. (16)). Consequently, we calculate the atomization energies for selected systems (Table II) using the ACLSD functional with either the NLR or the shell model for the strong correlation limit of the AC. Upon substituting the NLR by the shell model in the ACLSD functional (cf. Eq. (16)) all atomization energies are increased, in many cases even leading to overbinding which is reflected in the overall mean signed error of only 2.78 kcal/mol and a mean absolute error of 18.36 kcal/mol. The MAE for congested systems is improved by a factor 2 and the MAE for the atomization energies of single-bonded molecules is reduced from 14.82 kcal/mol with the NLR to 10.69 kcal/mol with the shell model. As mentioned above, the exchange-correlation energy of the shell model increases when the more recently calculated upper bound for the XC energy of the HEG in the SCL is used. Consequently, the overbinding tendency is less pronounced and the MAE for the single-bonded molecules is slightly improved (7.83 kcal/mol)

H2 LiH CH4 NH3 OH H2O HF Li2 LiF

Expt.

SM

NLR

109 58 419 297 107 232 141 24 139

113.07 58.69 460.56 320.84 110.26 240.69 141.36 21.27 127.99

110.69 51.66 409.57 275.85 94.29 207.01 122.66 17.00 107.28

10.69 7.64

14.82 −14.44

410.96 601.77 297.85 236.00 195.47 117.59 82.94 −3.95 129.01 66.77

363.57 530.28 259.02 207.26 164.99 89.30 57.26 −21.82 74.76 30.79

25.26 −12.16

50.06 −50.06

18.36 −2.78

33.36 −33.19

MAE MSE C2H2 C2H4 HCN CO N2 NO O2 F2 P2 Cl2

405 563 312 259 229 153 121 39 117 58

MAE MSE MAE total MSE total

while the congested systems are again all underbound (MAE of 33.83 kcal). The total energies for selected atoms are shown in Table III. For the ACLSD with the NLR model, the MAE of 0.063 hartree is comparable to the PBE hybrid functional. The energies are systematically underestimated. Due to the TABLE III. Total energies in hartree obtained with ACLSD16 using the non-local radius model (NLR) and the shell model (SM). Experimental energies49–51 and mean absolute errors are shown for comparison. ACLSD

H Li C N O F Cl P MAE

Expt.

SM

NLR

−0.5 −7.478 −37.845 −54.589 −75.067 −99.734 −460.148 −341.259

−0.499 −7.520 −37.933 −54.703 −75.207 −99.903 −460.819 −341.781

−0.499 −7.512 −37.892 −54.641 −75.117 −99.777 −460.294 −341.395

0.218

0.063

124104-7

Bahmann, Zhou, and Ernzerhof

J. Chem. Phys. 145, 124104 (2016)

more negative exchange-correlation energies obtained with the shell model, this underestimation is even more pronounced when the shell model is used in the interpolation scheme. The MAE amounts to 0.218 hartree. With the less negative exchange-correlation energies from the shell model based on the more recent SC-HEG limit, the total energies are improved, yielding a MAE of 0.145 hartree.

γ µ , and a set of allowed magnetic quantum numbers i, j, k such that l = i + j + k, χ µ (r A) = Nµ (x − X A)i ( y − YA) j (z − Z A)k
× exp γ µ |r − R A |2 .

(32)

Introducing the auxiliary vectors S A = S A x Sˆ x + S A y Sˆ y + S A z Sˆ z, S B = SB x Sˆ x + SB y Sˆ y + SB z Sˆ z,

V. IMPLEMENTATION

Both the NLR and the shell model were implemented into a developers version of the TURBOMOLE program package. For the convenience of the interested reader samples of the new integral routines can be retrieved from this website.52 First, we will show how ρ(r,u), underlying the NLR and the shell model, can be calculated efficiently within a Gaussian basis set. Subsequently the analytical integrals for the normalization and the energy density as defined in Eqs. (20) and (26) are derived and the numerical stability of their implementation will be addressed. The basic integrals are essentially the same for both, the NLR and the shell model. As will be explained in more detail below, the former requires only integrals to the upper limit uc , while the shell model can be expressed as a linear combination of integrals up to the limits of u s and uc .

In the context of the exchange-factor model to approximate exchange energies, it has been shown how the integration for the spherically averaged electron density can be done analytically in programs based on Gaussian basis sets.23 This section is mostly an elaboration of these results which are then subsequently used to derive analytical integrals over the interelectronic distance u. ρ(r,u) is defined in terms of the density matrix P (in the atomic orbital basis χ µ ) and integrals over pairs of basis functions  ρ(r,u) = Pµν Aµν (r,u), (28) µν

(34)

as well as the function G µν (S A, S B) = exp(−γ µr 2A − γν r 2B) exp −(γ µ + γν )u2 × exp(S A.r A + S B .r B)

sinh(w u) , wu


(35)

where w = |S A + S B − 2γ µ r A − 2γν r B |,

(36)

the spherical average over a pair of primitive Cartesian GTOs with azimuthal quantum numbers l and l˜ is defined by the partial derivative of G µν (S A, S B) with respect to the components of the auxiliary vectors in the limit of S A = S B = 0, Aµν (r,u) ˜

A. The spherically averaged density

(33)

˜

˜

∂i ∂ j ∂k ∂i ∂ j ∂k G µν (S A, S B)+/ . = *. i ∂S A x ∂S Aj ∂S Ak z ∂SBi˜ ∂S j˜ ∂SBk˜ x z B y , -S A=S B =0 y (37) Since all Cartesian GTOs belonging to the same shell share the same exponent and origin, all integrals that belong to a given shell pair with quantum numbers l and l˜ are calculated simultaneously. To simplify the partial derivative of G(Sa , S B) in Eq. (35), we rewrite it as
G µν (S A, S B) = K µν exp −γu2 × exp(S A.r A + S B .r B)F(w,u),

(38)

with the constant prefactor

where A denotes the integral over angular coordinates of u of a pair of basis functions at r + u,  dΩu χ µ (r A + u) χν (r B + u). (29) Aµν (r,u) = 4π The basis functions are located at atoms A and B with coordinates R A and R B, respectively, and r A = r − R A = (x − X A)ˆx + ( y − YA)ˆy + (z − Z A)ˆz,

(30)

r B = r − R B = (x − X B)ˆx + ( y − YB)ˆy + (z − Z B)ˆz.

(31)

For simplicity, an uncontracted Gaussian type orbital (GTO) basis is assumed in the following. Contraction coefficients can be included as a prefactor in the matrix elements defined in Eq. (29) and the extension to contracted basis sets is straightforward. A primitive Cartesian GTO that belongs to a shell with the azimuthal quantum number l is defined by its exponent

K µν = exp(−γ µr 2A − γν r 2B),

(39)

the sum of exponents γ = γ µ + γν ,

(40)

and the auxiliary function F(w,u) =

sinh(w u) . wu

(41)

Note that F depends through w on the components of S A and S B and the partial derivatives (Eq. (37)) apply only to the last two factors of Eq. (38). Using the chain rule for derivatives of F, the rhs of Eq. (37) can be simplified to a sum over derivatives of F(w,u) with respect to w in the limit of S A = S B = 0 and coefficients that are independent of the interelectronic distance u

124104-8

Bahmann, Zhou, and Ernzerhof

J. Chem. Phys. 145, 124104 (2016)

˜˜˜

jki jk G iµν

j˜

i˜

k˜

∂ ∂ ∂ ∂ ∂ ∂ = *. i G(S A, S B)+/ ˜ ˜ j k i j ∂S A x ∂S A ∂S A z ∂SB ∂S ∂SBk˜ x z By y , -S A=S B=0 i

j

k

= K µν exp(−γu2)

l+l˜ 

˜˜˜

Cni j k i j k F (n)(w0,u),

(42)

n=0

where the limiting case of w is abbreviated as w0 = w(S A = 0, S B = 0) = | − 2γ µ r A − 2γν r B |.

(43)

In Eq. (42) F (n)(w0,u) denotes the nth derivative of F with respect to w0 which is equivalent to taking the derivative with respect to w and setting the auxiliary vectors to zero, F (n)(w0,u) =

∂n ∂n F(w,u)|S A=S B =0 = F(w0,u). n ∂w ∂w0n

(44)

The coefficients contain derivatives of the exponential function exp(S A.r A + S B .r B) and w with respect to the components of the auxiliary vectors in the limit S A = S B = 0. As an example, ˜ the zero order and the maximum (l + l)th-order coefficients ˜ for F (0) and F (l+l), respectively, in Eq. (42) are given by ˜˜˜ C0i j k i j k

=

˜ ˜ ˜ r iA x r Aj y r kA z r Bi x r Bj y r Bk z

˜˜˜

˜

(45) ˜

˜

i jki jk = (w 100)i+i (w 010) j+ j (w 001)k+k . Cl+ l˜

(46)

For the other coefficients, explicit expressions for given ˜ j, ˜ k˜ are obtained with a symbolic values of i, j, k and i, math program.53 In the quantum chemistry program, they are precalculated as functions of r A and r B and derivatives of w before evaluating F(w0,u) for different values of u or integrals over the u-dependent terms in Eq. (42). From the definition of w (Eq. (36)) follows that ∂w = ∂S∂wB , ∂S∂w = ∂S∂wB , and ∂S∂w = ∂S∂wB and thus ∂S Ax

w

x

Ay

y

Az

z

˜ ˜ j+ j)(k+ ˜ (i+i)( k)

B. Recursion scheme for the normalization and energy integrals

The normalization of the exchange-correlation hole is used in both models to determine the cut-off radius uc . For the NLR hole it is given in terms of the integrals over spherically averaged basis-function pairs introduced above  uc NLR N (r,uc ) = −4π u2 ρ(r,u)du 0  uc  = −4π Pµν u2 Aµν (r,u)du. (48) 0

µν

The same type of integrals over two different ranges (u s and uc ) are required for the shell model and the normalization integral is just a linear combination thereof, ( u c )  us shell 2 2 N (r,u s ,uc ) = 4π u ρ(r,u)du − 2 u ρ(r,u)du 0 0 ( uc  = 4π Pµν u2 Aµν (r,u)du 0

µν



−2

us

) u2 Aµν (r,u)du .

(49)

0

We will thus derive the recursion scheme for the uc -dependent integral only. j k i˜ j˜k˜ Inserting the definition for Aµν (r,u) in terms of G iµν and the simplified expression for the latter from the previous section (Eqs. (37) and (42)) into the normalization integral over a pair of Cartesian GTOS gives  4π

uc

u Aµν (r,u) du = K µν 2

0

l+l˜ 

˜˜˜

Cni j k i j k ifnn (w0, γ,uc ),

(50)

n=0

where the integrals depending on the exponent γ, w0 (cf. Eqs. (40) and (43)), and uc are defined as  uc ifnn (w0, γ,uc ) = 4π u2 exp(−γu2)F (n)(w0,u)du. (51) 0

˜

˜

˜

∂i ∂ j ∂k ∂i ∂ j ∂k + = *. i w/ . ˜ j˜ ∂S ∂S j ∂S k ∂S i˜ ∂S k , A x A y A z B x ∂SB y B z -S A=S B =0

(47)

That is, only partial derivatives with respect to 3 instead of 6 variables are required for the subsequent computation of the coefficients. Closer inspection of Eq. (42) reveals that for a pair of ˜ only l + l˜ basis functions with quantum numbers l and l, derivatives of F(w0,u) need to be calculated additionally to the function itself, i.e., for a pair of two s-functions only one, for a ps-pair two (F and its first derivative), and for a pp-pair three (F and up to the second derivative of F) functions of u need to be computed, etc. This reduces the computational effort particularly if the spherically averaged density is needed for many values of u. Also, the spherical average over Cartesian GTOs belonging to a given shell pair with quantum numbers l and l˜ differs only in the coefficients in Eq. (42). It is thus recommended to calculate all respective integrals simultaneously as is commonly done for other kinds of two- and one-electron integrals.

Using the Leibniz identity to evaluate F (n)(w0,u), the integrals are expressed as a sum over auxiliary integrals, ifnn (w0, γ,uc ) n 
−
+  (−1)n−k n!  1 + (−1)k Ik+1 + 1 + (−1)k+1 Ik+1 , = 2π n−k+1 k!w0 k=0 (52) which are given by  uc Ik− = uk exp(−γu2) sinh(w0 · u) du, 0  uc Ik+ = uk exp(−γu2) cosh(w0 · u) du,

(53) (54)

0

and calculated recursively starting from I0± for k ≥ 2, Ik− =

uk−1 w0 + k −1 − Ik−1 + Ik−2 − c exp(−γu2) sinh w0u, 2γ 2γ 2γ

(55)

Ik+ =

uk−1 w0 − k −1 + Ik−1 + Ik−2 − c exp(−γu2) cosh w0u. 2γ 2γ 2γ

(56)

124104-9

Bahmann, Zhou, and Ernzerhof

J. Chem. Phys. 145, 124104 (2016)

FIG. 4. Algorithm for the calculation of the normalization integral (energy density) for the non-local radius and shell model at each grid point. For clarity an implicit loop over n = 0, . . ., l + l˜ is assumed at each occurrence of the index n. If the normalization condition is not fulfilled, the algorithm is repeated with an improved u c value.

The special cases for k = 1 are w0 + exp(−γu ) sinh w0u I − , (57) 2γ 0 2γ w0 − exp(−γu2) cosh w0u − 1 I1+ = I − . (58) 2γ 0 2γ From Eq. (52), it follows that for the normalization integral only Ik− with impair k and Ik+ with pair k are needed. Once the range uc of the hole is determined, the exchangecorrelation energy density for the NLR model is calculated through  uc NLR ϵ XC (r) = −2π u ρ(r,u) du 0  uc  = −2π Pµν u Aµν (r,u) du. (59) 2

I1− =

0

µν

Analogously to the normalization integral, an explicit expression in terms of integrals over pairs of Cartesian GTOs is derived using Eqs. (37) and (42),  2π

uc

uAµν (r,u)du = K pq

0

l+l˜ 

˜˜˜

Cni j k i j k ifen (w0, γ,uc ).

(60)

This procedure is rather cheap since, after calculation of the coefficients in Eq. (50), the ifnn (Eq. (51)) integrals can be calculated efficiently for many values of uc . Due to numerical instabilities for small w0 inherent to the normalization and energy integrals ifnn and ifen (as well as the spherically averaged density), a Taylor series is used in these cases. To avoid jumps or oscillations and thus convergence problems in the iterative procedure to find uc , the Taylor expansion is used for all values of uc if w0 is below a given threshold. This happens for basis functions with small exponents at grid points close to the nucleus. Additionally the normalization integral becomes unstable for some combinations of w0, γ, and uc that lead to very small values of the first auxiliary integral I0+ in the recursive formula (Eq. (56)). Due to numerical limitations, the integral becomes zero and spuriously large negative integral values of ifnn ensue (cf. Eq. (52)). For values of I0+ below a certain threshold, all ifnn integrals are thus calculated numerically using a 11-point Gauss-Legendrequadrature with precalculated weights and abscissae. The algorithm to calculate the normalization and energy integrals defined at each grid point in Eqs. (48) and (59) is outlined in Fig. 4.

n=0

Finally, the energy integrals are obtained from the previously defined auxiliary integrals Ik±, (Eq. (54)) ifen (w0, γ,uc ) = π

n  (−1)n−k n!

k!w0n−k+1 

 × 1 + (−1)k Ik− + 1 + (−1)k+1 Ik+ . k=0

(61)

In contrast to the normalization, for the energy integral, Ik− with pair k and Ik+ with impair k are required. C. Efficiency and numerical stability

Within the NLR and the shell model, the range of the exchange-correlation hole is obtained from the normalization condition (cf. Eq. (20)) and uc is determined iteratively using a bisection algorithm. To speed up the calculation, a scan over the whole range of possible uc values is performed and the bisection starts at the closest values bracketing −1.

VI. CONCLUSIONS AND OUTLOOK

At the center of the present work is the strong-correlation limit of the exchange-correlation hole. We employ a nonempirical approach to the construction of ρSC XC(r,u) in which the spherically averaged density is multiplied by an exchangecorrelation factor that turns the density into an exchangecorrelation hole. This procedure resembles an earlier approach to the construction of the exchange and exchange-correlation hole through the exchange and correlation factor ansatz, respectively.23–25 The shell model improves upon the NLR model by respecting the strong correlation limit of the HEG while all advantageous aspects of NLR are retained, including exactness of the corresponding energy density for one-electron systems as well as for infinitely stretched H2. The homogeneous limit is a corner stone in the first-principles construction of exchange-correlation functionals and virtually all existing first-principles functionals aim at reproducing it.

124104-10

Bahmann, Zhou, and Ernzerhof

While the NLR model is a very short-ranged hole with maximum depth that removes the electron density around a reference point, the shell model extends the range of the exchange-correlation hole by adding a positive part to it. Other than the additional formal conditions satisfied by the shell model, it also provides improvements in accuracy. For a selection of atoms, the shell model yields better strongcorrelation exchange-correlation energies compared to NLR. Employing the shell model to the SCL, we construct the ACLSD adiabatic connection interpolation and compare it to the analogous interpolations drawing on the NLR model. While atomization energies are significantly improved with the functional based on the shell model, the total energies for atomic systems are further underestimated. It has already been noted that the simple ACLSD interpolation has various limitations that are addressed in ongoing work. For instance, strong-correlation in the core region, which is the main source of errors for total energies, can be suppressed additionally. Similarly, the atomization energies obtained with ACLSD in conjunction with the shell model can be improved through the construction of a suitable local mixing coefficient a(r) that accounts for static correlation. ACLSD building on the shell model satisfies various desirable conditions such as exactness in the high and low density limit, exact for one-electron systems as well as for the HEG. Therefore it should provide a better starting point for the construction of first principles approximations than LSD itself. Important for applications, we show how the shell model (and the parent NLR model) can be implemented efficiently for Gaussian basis sets. The integrals for the normalization of the exchange-correlation hole and the exchange-correlation energy density are calculated analytically using a recursion scheme. Numerical limitations (e.g., for grid points near a nucleus) are circumvented by using a series expansion or by resorting to numerical integration. ACKNOWLEDGMENTS

We would like to acknowledge the financial support provided by The Natural Sciences and Engineering Research Council of Canada (NSERC) and the Deutsche Forschungsgemeinschaft (DFG). H.B. thanks the Turbomole consortium and Martin Kaupp for access to the source code of the Turbomole program package. 1W.

Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).

2R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules

(Oxford University Press, 1989). Koch and M. C. Holthausen, A Chemist’s Guide to Density Functional Theory, 2nd ed. (Wiley-Vch, Weinheim, 2004). 4R. Martin, Electronic Structure Basic Theory and Practical Methods (Cambridge University Press, 2004). 5A. D. Becke, J. Chem. Phys. 140, 18A301 (2014). 3W.

J. Chem. Phys. 145, 124104 (2016) 6M.

Seidl, J. P. Perdew, and M. Levy, Phys. Rev. A 59, 51 (1999). Seidl, Phys. Rev. A 60, 4387 (1999). 8M. Seidl, J. P. Perdew, and S. Kurth, Phys. Rev. A 62, 012502 (2000). 9A. Mirtschink, M. Seidl, and P. Gori-Giorgi, J. Chem. Theory Comput. 8, 3097 (2012). 10F. Malet, A. Mirtschink, K. J. H. Giesbertz, L. O. Wagner, and P. Gori-Giorgi, Phys. Chem. Chem. Phys. 16, 14551 (2014). 11L. O. Wagner and P. Gori-Giorgi, Phys. Rev. A 90, 052512 (2014). 12A. D. Becke, J. Chem. Phys. 138, 074109 (2013). 13D. Langreth and J. Perdew, Solid State Commun. 17, 1425 (1975). 14P. Gori-Giorgi and M. Seidl, Phys. Chem. Chem. Phys. 12, 14405 (2010). 15S. Vuckovic, T. J. P. Irons, A. Savin, A. M. Teale, and P. Gori-Giorgi, J. Chem. Theory Comput. 12, 2598 (2016). 16Y. Zhou, H. Bahmann, and M. Ernzerhof, J. Chem. Phys. 143, 124103 (2015). 17J. A. Alonso and L. A. Girifalco, Phys. Rev. B 17, 3735 (1978). 18J. Alonso and L. Girifalco, Solid State Commun. 24, 135 (1977). 19O. Gunnarsson, M. Jonson, and B. Lundqvist, Solid State Commun. 24, 765 (1977). 20J. Charlesworth, Phys. Rev. B 53, 12666 (1996). 21Z. Guang, S. Clark, S. Brand, and R. Abram, Chin. Phys. Lett. 24, 807 (2007). 22Z. Wu, R. E. Cohen, and D. J. Singh, Phys. Rev. B 70, 104112 (2004). 23H. Antaya, Y. Zhou, and M. Ernzerhof, Phys. Rev. A 90, 032513 (2014). 24J. Precechtelova, H. Bahmann, M. Kaupp, and M. Ernzerhof, J. Chem. Phys. 141, 111102 (2014). 25J. Precechtelova, H. Bahmann, M. Kaupp, and M. Ernzerhof, J. Chem. Phys. 143, 144102 (2015). 26M. Ernzerhof, Chem. Phys. Lett. 263, 499 (1996). 27M. Ernzerhof, J. P. Perdew, and K. Burke, Int. J. Quantum Chem. 64, 285 (1997). 28A. D. Becke, J. Chem. Phys. 98, 1372 (1993). 29P. Mori-Sánchez, A. J. Cohen, and W. Yang, J. Chem. Phys. 124, 091102 (2006). 30M. Ernzerhof, in Lecture Notes in Physics, edited by D. Joubert (Springer, Berlin, Heidelberg, 1998), Vol. 500, pp. 60–90. 31S. Grimme, J. Chem. Phys. 124, 034108 (2006). 32K. Sharkas, J. Toulouse, and A. Savin, J. Chem. Phys. 134, 064113 (2011). 33J. C. Sancho-Garcia and C. Adamo, Phys. Chem. Chem. Phys. 15, 14581 (2013). 34J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). 35M. Lewin and E. H. Lieb, Phys. Rev. A 91, 022507 (2015). 36M. Seidl, S. Vuckovic, and P. Gori-Giorgi, Mol. Phys. 114, 1076 (2016). 37E. J. Baerends, Phys. Rev. Lett. 87, 133004 (2001). 38S. Vuckovic, L. O. Wagner, A. Mirtschink, and P. Gori-Giorgi, J. Chem. Theory Comput. 11, 3153 (2015). 39TURBOMOLE V6.6 2014, A Development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007, available from http://www.turbomole.com. 40R. Ahlrichs, M. Bär, M. Häser, H. Horn, and C. Kölmel, Chem. Phys. Lett. 162, 165 (1989). 41J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996). 42R. Krishnan, J. S. Binkley, R. Seeger, and J. A. Pople, J. Chem. Phys. 72, 650 (1980). 43A. D. McLean and G. S. Chandler, J. Chem. Phys. 72, 5639 (1980). 44T. H. Dunning, Jr., J. Chem. Phys. 90, 1007 (1989). 45D. Woon and T. H. Dunning, Jr., J. Chem. Phys. 100, 2975 (1994). 46H. Bahmann and M. Kaupp, J. Chem. Theory Comput. 11, 1540 (2015). 47J. P. Perdew, M. Ernzerhof, and K. Burke, J. Chem. Phys. 105, 9982 (1996). 48M. Ernzerhof and G. E. Scuseria, J. Chem. Phys. 110, 5029 (1999). 49C. L. Pekeris, Phys. Rev. 126, 1470 (1962). 50C. Schwartz, Int. J. Mod. Phys. E 15, 877 (2006). 51S. J. Chakravorty, S. R. Gwaltney, E. R. Davidson, F. A. Parpia, and C. F. Fischer, Phys. Rev. A 47, 3649 (1993). 52See http://www.chemie.tu-berlin.de/dr_hilke_bahmann/home/ for sample routines for the integrals described in Sec. V. 53Wolfram Research, Inc., Mathematica 10.2, Champaign, Illinois, 2015. 7M.