Simple Fully Nonlocal Density Functionals for Electronic Repulsion Energy.

From a simplified version of the mathematical structure of the strong coupling limit of the exact exchange-correlation functional, we construct an approximation for the electronic repulsion energy at physical coupling strength, which is fully nonlocal. This functional is self-interaction free and yields energy densities within the definition of the electrostatic potential of the exchange-correlation hole that are locally accurate and have the correct asymptotic behavior. The model is able to capture strong correlation effects that arise from chemical bond dissociation, without relying on error cancellation. These features, which are usually missed by standard density functional theory (DFT) functionals, are captured by the highly nonlocal structure, which goes beyond the "Jacob's ladder" framework for functional construction, by using integrals of the density as the key ingredient. Possible routes for obtaining the full exchange-correlation functional by recovering the missing kinetic component of the correlation energy are also implemented and discussed.

The widespread success of Kohn–Sham density-functional theory (KS DFT)[1−4] across various chemical and physical disciplines has been also accompanied by spectacular failures,[2] reflecting fundamental issues in the present density functional approximations (DFAs) for the exchange–correlation (XC) functional. Well-known examples are the paradigmatic case of the dissociation curves of the H2 and H2+ molecules.[2,5] The usual DFA approach to construct XC functionals consists in making an ansatz in terms of “Jacob’s ladder” ingredients:[4,6−8] the local density, its gradient, its laplacian, and/or KS kinetic energy density, up to occupied and virtual KS orbitals. While this strategy has been very successful for moderately correlated systems (see, e.g., refs (3, 4, and 9−11).), it has failed so far when correlation effects become important (e.g., in stretched bonds, but also at equilibrium gemoetries when partially filled d and f subshells are present). This fact suggests that a different approach to DFAs is needed to address the problem of strong correlation.[2,3,12−14]

The strong-interaction limit of DFT[15−18] provides information on how the exact XC functional depends on the density in a well-defined mathematical limit, which is relevant for strong correlation. The thorough explorations of this limit reveal a mathematical structure totally different from that of Jacob’s ladder ingredients. Instead of the local density, density derivatives, or KS orbitals, in this limit we see that certain integrals of the density play a crucial role, encoding highly nonlocal information,[15−17] embodied in the so-called strictly correlated electrons (SCE) functional.[15−17] This functional appears to be well-equipped for solving long-standing DFAs problems: it is self-interaction free, it captures the physics of charge localization due to strong correlation without resorting to symmetry breaking,[19−21] and its functional derivative displays (in the low-density asymptotic limit) a discontinuity on the onset of fractional particle number.[22] Despite these appealing features, there are two main obstacles to the routine use of the SCE functional: its availability is restricted to small systems[16,23] and its energies are way too low for most of physical and chemical systems.[19,23−25] The nonlocal radius (NLR) functional,[26] and the newer shell model[27] are inspired to the SCE functional form and retain only some of its nonlocality. They are readily available,[27] but, being approximations to the SCE functional, their energies are also too low with respect to those of chemical systems.[26,27] The information encoded in the SCE functional or its approximations can be combined with the complementary information from the weak coupling limit. This has been recently used for constructing XC functionals from a local interpolation along the adiabatic connection.[14,27−29] Although this approach is promising for treating strong correlation within the realm of DFT,[14] it can still easily overcorrelate (for example for stretched bonds it overcorrelates the fragments), again because the SCE (exact or approximate) quantities are often far from the physical ones.[29]

Nonetheless, the way in which the information encoded in the density is transfomed into an electron–electron repulsion energy in the SCE functional is very intriguing, with many physical appealing features.[16,30−32] Motivated by this observation, in this letter we use the SCE mathematical structure to devise a new way to design fully nonlocal approximate density functionals for the electronic interaction energy at the physical coupling strength. Capturing the main structural motives of the SCE functional, we preserve many of its appealing features, but with repulsion energies that are much closer to those of physical systems. Moreover, besides accurate total repulsion energies, our model provides energy densities within the definition of the electrostratic potential of the XC hole that are also locally very close to exact ones, making it an ideal tool for the development of functionals that use the exact exchange energy density, like hyperGGA’s[33−36] or local hybrids.[37,38] In other words, it is known[33,35,39,40] that in order to use the exact exchange energy density, we need a fully nonlocal correlation functional compatible with it. It is the purpose of this work to provide a new strategy to build this fully nonlocal functional at a computational cost similar to the one of the exact exchange energy density.

In order to explain our construction, we have to first quickly review some basic DFT equations. An exact expression for the XC energy can be obtained from the density-fixed adiabatic connection formalism (AC):[41,42]where Wλ[ρ] is the global (i.e., integrated over all space) AC integrand:The wave function Ψλ[ρ] depends on the positive coupling constant λ and minimizes ⟨T̂ + λV̂ee⟩, while integrating to ρ(r), the density of the physical system (λ = 1). This way, AC links the KS non-interacting state described by Ψ0[ρ] and the physical state described by Ψ1[ρ]. It also further connects the physical and the SCE state, i.e., the state of perfect electron correlation, corresponding to the limit λ → ∞. The XC energy densities along the adiabatic connection (i.e., position-dependent quantities wλ(r) that integrate to Wλ[ρ] when multiplied by the density) are not uniquely defined and therefore we have to be specific on their gauge.[14,29,43−45] A physically sound gauge often considered in DFT is the one of the electrostatic potential of the XC hole.[14,34,35,46] Within this gauge we can express the λ-dependent energy density in terms of the corresponding spherically averaged XC hole[14,34,35,47]hxcλ(r, u) obtained from |Ψλ[ρ]|2,where u = |r – r′| is the distance from a reference electron at r. In the SCE (λ→ ∞) limit the energy density wλ(r) in the gauge of eq has the exact form[48]where vH(r) is the Hartree potential and the comotion functionsf([ρ]; r) are nonlocal functionals of the density that give the positions of the remaining N – 1 electrons when one electron is at position r.[15,16,32] From eq we see that in the λ → ∞ limit the energy density w∞(r) is fully determined by the distances RSCE([ρ]; r) = |r–f([ρ]; r)| between a reference electron in r and the remaining N – 1 ones. For example, in the case of one-dimensional systems, the distances RSCE([ρ]; x) can be constructed exactly[15,49] from the equations (with i = 2, ..., N)which can be solved in terms of the function N1D(x) = ∫–∞ρ(x′) dx′ and its inverse N1D–1(y).[15,19,30] We see that in this limit each electron is separated by the closest one by a piece of density that integrates exactly to 1 (in other words, fluctuations are totally suppressed in the limit of extreme correlation), with the key ingredient being the amount of expected electrons between two electronic positions.

In this work we propose a way to generalize the SCE form of eq by using λ-dependent distances (or “radii”) Rλ([ρ]; r) that will take into account the effect of fluctuations, which are not as suppressed as in the extreme SCE case. Thus, our “multiple-radii functional” (MRF) energy density reads asAs we shall see, we will determine the Rλ([ρ]; r) by using a simplified version of the same kind of integrals of the density that appear in the SCE limit, introducing the average effect of fluctuations by reducing the amount of expected charge between two electronic positions. Before coming to the details of the Rλ([ρ]; r) construction, we remark that eq can be also derived from the following model for the spherically averaged pair-densitywhere δ is the Dirac delta function. Given that the model of eq is properly normalized, the corresponding λ-dependent XC hole satisfies the sum rule, integrating to −1 electron.

We now turn to the construction the radii Rλ([ρ]; r). Motivated by the structure of eq , and similarly to the recent nonlocal approximations for the SCE functional for three-dimensional systems,[26,27] we introduce the spherically averaged density ρ̃(r, u) around a position r,and the function Ne(r, u),These functions have been studied and efficiently implemented by Ernzerhof and co-workers.[27,28,51] We now want to find a physical approximation for the crucial quantitieswhich give the expected number of electrons in a sphere of radius Rλ(r) centered at the reference electron in r. We notice at this point that in an inhomogeneous system, even one-dimensional, it is not possible to write the exact SCE radii explicitly in terms of the function Ne(r, u) obtained by spherically averaging the density around a reference electron in r as in eq . An exception is an homogeneous 1D system, in which ν2SCE = ν3SCE = 2, ν4SCE = ν5SCE = 4, etc.

At the physical interaction strength λ = 1, we expect a situation in which this extreme correlation is reduced, with all the ν close to i – 1. To illustrate this fact, we consider first an N = 2 system, for which we have only one radius, R2λ, which from eq will be equal toshowing that, for N = 2, R2λ(r) is the screening length associated with the Hartree-exchange-correlation potential when its response part is removed.[52−54] Given that for two-electron systems highly accurate wλ(r) have been computed,[14,50] we can use eq to obtain the “exact” R2λ(r). In Figure we show the corresponding ν2λ(r) for the hydride ion at λ = 0, λ = 1, and λ → ∞. We clearly see that, in the physical system, ν2λ = 1(r) is much closer to 1 than in the SCE extreme case, which allows for much larger numbers ν2∞(r), becoming equal to 2 at the nucleus, as in the homogeneous 1D solution.

Figure 1

ν2λ(r) quantity of eq and the corresponding R2λ(r) radius at different coupling strengths for the hydride ion obtained from eq using accurate xc energy densities wλ(r) from refs (14 and 50)

It is clear that there are several ways to define approximations for νλ(r), using different ingredients. Here, our aim is to show that already very simple approximations can yield rather accurate results, and we focus on the physical λ = 1 case. As said, we expect that ν2λ = 1(r) ≈ 1, ν3λ = 1(r) ≈ 2,..., and we writeyielding for the radii Rλ = 1(r) the equationswith σ(r) being the fluctuation function, which can push away or bring closer the ith electron to the reference one with respect to the expected distance a(r) = N–1(r, i – 1). In this first model, we consider only the case in which the ith electron is pushed further, because for this case we can use again the mathematical structure of the SCE functional as a guide. More general models will be explored in future works. From the SCE theory for spherically symmetric systems,[16,32] we know that the derivative of the radial comotion function f(r) at point r is inversely proportional to 4πf(r)2ρ(f(r)). We thus introduce the quantity S(r)which, in analogy to the SCE structure, provides information on the derivative of the R1(r) at σ = 0. When S(r) is small, the derivative of the R(r) will be very large, and we expect the electron to be pushed further, with σ approaching the average value 1/2 (which is exactly in between two expected positions). When S(r) is large, the derivative of R(r) is very small, and we expect it to stay close to σ = 0 (or even become slightly negative, a possibility not considered here). Thus, for constructing the MRF functional at the full coupling strength, hereinafter the MRF-1 functional, W1MRF[ρ] = ∫ρ(r)w1MRF(r)dr, we use a simple Gaussian ansatzwhere b = 5 has been chosen to optimize the He atom W1[ρ]. Equations , 8, 9, and 12–15 completely define w1MRF(r).

In Table we compare W1[ρ] obtained with the MRF-1 model with corresponding reference values (full-CI/CCSD), PBE, and SCE ones (W∞[ρ]) for several closed-shell atomic (ionic) systems evaluated on accurate (CCSD/full-CI) densities. The PBE values have been obtained by using the scaling relation,[58,59] with ργ(r) = γ3ρ(γ r),From Table we can see that our model, even with the very simple ansatz for σ(r) of eq , gives repulsion energy much closer to the physical ones with respect to SCE. Their quality is comparable to that of PBE, with MAE somewhat smaller (0.17 au vs 0.24 au). The purpose here is not to reach high accuracy (which requires optimization and further studies of the R), but to show that functional approximations based on modeling the quantity σ(r) is a very promising strategy, because already a primitive nonoptimized model performs very well. Even more interesting than the global W1MRF[ρ] values are the energy densities: in Figure we compare w1MRF(r) with the reference w1(r) for the neon atom (top panel) and the hydride ion (bottom panel). We also show w1(r) obtained with the LDA functional from the PW92 parametrization,[57,60]. We see that w1MRF(r) is in good agreement with the reference w1(r) in the case of Ne, but also in the more correlated case[14,48] of H–, again improving dramatically with respect to SCE. From the insets of the same figure, we can see that the local error of our model is very small, vanishing for large r due to the correct asymptotic behavior, arising from the proper normalization of eq . The availability of DFAs energy densities in this gauge is rather limited, and beyond LDA it is restricted to few approximations[61,62] to the exchange energy density (ϵ(r) = w0(r)). For instance, the gauge incompatibility[63] of the generalized gradient approximation (GGA) exchange energy densities and the exact ϵ(r), which is in the gauge of eq , has been a major hurdle for the development of local hybrid DFAs.[37,38]

Table 1

Atomic (Ionic) Repulsion Energies W1[ρ] Obtained by the MRF-1 Model and PBE Compared to Reference W1[ρ], Obtained with the Gamess-US Package[55] Using Full-CI (for the First Four Systems) and CCSD Wavefunctions (Other Systems)a

atom/ion	reference	MRF-1	PBE	SCE
He	–1.1029	–1.1844	–1.1047	–1.4982
H–	–0.4532	–0.4681	–0.4413	–0.5689
Be	–2.8341	–2.8044	–2.8430	–4.0195
Li–	–1.9462	–2.1170	–1.9617	–2.7308
F–	–10.889	–10.741	–10.997	–16.940
Ne	–12.765	–12.823	–12.876	–20.041
Mg	–16.701	–16.365	–16.913	–26.709
Cl–	–28.89	–28.48	–29.19	–47.26
Ar	–31.35	–31.19	–31.68	–51.49
Ca	–35.60	–35.92	–36.85	–60.34
MAE	-	0.17	0.24	-

The aug-cc-pCVXZ basis set of Dunning[56] has been used (X = 6 for He and H–, X = 5 for F– and Ne, X = T for Be and Li– and X = Q for the other atoms). The SCE values W∞[ρ] computed from the same densities are also reported.

Figure 2

Energy densities at full coupling strength w1(r) as a function of distance from the nucleus, r/a.u., obtained from the present model (MRF), from the local-density approximation (LDA),[57] and from the strictly correlated electrons functional (SCE), all evaluated on accurate densities, for Ne (upper panel) and H– (lower panel). The reference w1(r) are obtained at the full-CI and CCSD level of theory, as in refs (14 and 48), by using the aug-cc-pCVTZ and aug-cc-pV6Z basis sets[56] for Ne and H–, respectively. Insets show the absolute error of approximate energy densities, δw1(r) = w1 – w1apx(r).

The main point of introducing the full nonlocal dependence is of course to treat static and strong correlation. In Figure we show the energy density w1MRF(r) for the H2 molecule at different bond lengths Lb along the internuclear axis, compared with accurate ones from ref (14). For comparison, we also show w1(r) obtained from the interpolation model of Liu and Burke (LB)[64] applied to energy densities, with the exact w0(r), w0′(r), and w∞(r) as input ingredients.[14] As we can see from Figure , in the equilibrium region the MRF-1 energy densities are still somewhat lower than the reference ones, whereas the LB is highly accurate, as it is also the case with atoms.[14,29] However, we can also see that in the stretched case (Lb = 10 au) the MRF-1 energy densities are very accurate, even more accurate than the LB interpolated ones, whose error is already small (note again that they use the exact w0(r), w0′(r) and w∞(r) as input for the interpolation). While in the stretched H2 molecule the static correlation effects are dominant, at intermediate bond lengths, around Lb ∼ 5.0 au, there is a subtle interplay between dynamic and static correlation effects.[29] This region can be even more challenging for DFAs than the stretched case, given that certain DFAs that dissociate H2 correctly fail in this scenario yielding a positive “bump” (see, e.g., refs (14 and 65−67)). We can see that the MRF-1 energy densities are very accurate at Lb = 5.0 au, hardly distinguishable from the reference ones.

Figure 3

Energy densities at full coupling strength as a function of the distance from the bond midpoint z along the internuclear axis for the H2 molecule at different bond-lengths Lb obtained with MRF-1 using accurate FCI/aug-cc-pCVTZ densities. Reference energy densities and those obtained with the LB local interpolation scheme are from ref (14).

In Figure we show the dissociation curve for the H2 molecule obtained using the MRF-1 functional evaluated on the accurate FCI/aug-cc-pCVTZ densities. We can see that around equilibrium MRF-1 underestimates the total energy, because it misses the positive kinetic correlation component Tc[ρ] and slightly underestimates the exact W1[ρ], as already shown in Figure . Despite missing Tc[ρ], MRF-1 dissociates H2 correctly, because Tc[ρ] vanishes as the H2 dissociates into atoms. To recover the missing Tc[ρ] component, one can combine W1MRF[ρ] with the quantities from the weak coupling limit, namely, W0[ρ] and W0′[ρ], to interpolate Wλ[ρ] and thus obtain Exc[ρ]. For this purpose, we employ a very simple interpolation form, the two-legged representation,[14,29,68] which has recently been used to construct a tight lower bound to correlation energies.[29] This form reads asAs in this work, we use W1MRF[ρ] as an approximation to W1[ρ], we call this approach the “2-leg MRF”, and from Figure we can see that it substantially improves the MRF-1 energies. In this case very similar results are obtained if we do the interpolation on the energy densities rather than on integrated quantities. Besides dissociating correctly H2, the dissociation of H2+ is also correctly described within the MRF-1 and 2-leg MRF approaches, because our model of eq is equal to 0 for all N = 1 systems.

Figure 4

H2 molecule dissociation curve as a function of the internuclear distance Lb/a.u. obtained with the MRF-1 and 2-leg MRF approaches presented in this work, compared to restricted PBE and FCI. All the curves have been obtained using the aug-cc-pCVTZ basis set.[56]

Finally, one may wonder if the MRF would encounter problems for extended systems. As a paradigmatic example, we consider the uniform electron gas (UEG) with density , for which Ne(r, u) = Ne(u) = u3/r3 and Ne–1(i – 1) = rs(i – 1)1/3. Then, by using the model pair density of eq , we obtain w1MRF(rs) = w̃(rs)/rs withwhere σ(rs) is given by eq and, from eq , S = 3(i – 1)2/3/rs. This expression has a fast N → ∞ convergence and when σ = 0 can be evaluated in closed form. It yields reasonable values for the UEG, with a maximum relative error of 23%. The function w̃(s) of eq has the same qualitative behavior of the exact one: it is monotonically decreasing with rs, bounded between the two limting values w̃(0) = −0.487 and w̃(∞) = −0.7564, not so far (considering the simplicity of the model for σ(rs)) from the exact ones, and ≈ – 0.876, respectively (this latter value is currently a matter of discussion, see refs (69 and 70)). The UEG also illustrates the physics of the model: when rs → 0 (weak correlation) σ → 0, while when rs → ∞ (strong correlation), σ ≠ 0 for i larger and larger (long-range fluctuations become more and more important). It also suggests that for extended systems the explicit functional can be confined to a set i < imax, and the rest can be resummed. The value imax is determined by correlation (for example, in the UEG it is automatically determined by σ(rs)).

In summary, we have proposed a strategy to build fully nonlocal DFAs inspired by the mathematical structure of the exact XC functional in the strong coupling limit, reducing the problem to the construction of the fluctation function σ(r) in terms of S(r) of eq . Already an extremely simple model such as the one of eq is locally accurate, it is able to dissociate correctly the H2 and H2+ molecules, and gives very reasonable results for the uniform electron gas. We thus believe that the nonlocal structure of our functional, which goes beyond the Jacob’s ladder framework, opens up new perspectives for the development of XC functionals able to tackle strong correlation. Although the functional is highly nonlocal, it can be obtained at a computational cost comparable to that of the NLR and shell functionals, which have been recently implemented in a very efficient way.[27] Many strategies to improve the accuracy can be pursued: the inclusion of kinetic correlation through interpolation along the adiabatic connection (as in Figure ); trying to model directly the λ-dependence of σ; the generalization to noninteger number of electrons[22] and spin densities; improving the accuracy for the UEG, and adding the dependence on the gradient of S(r). The functional can also be readily applied to other dimensionalities, e.g., electrons confined in quasi-1D and quasi-2D geometries, for which the SCE approach has already proven very useful.[19,20] It can be also applied to other isotropic interactions, such as the error function used in range separation[5] but also effective interactions for ultracold quantum gases.[21]