Exchange-Correlation Functionals via Local Interpolation along the Adiabatic Connection.

The construction of density-functional approximations is explored by modeling the adiabatic connection locally, using energy densities defined in terms of the electrostatic potential of the exchange-correlation hole. These local models are more amenable to the construction of size-consistent approximations than their global counterparts. In this work we use accurate input local ingredients to assess the accuracy of a range of local interpolation models against accurate exchange-correlation energy densities. The importance of the strictly correlated electrons (SCE) functional describing the strong coupling limit is emphasized, enabling the corresponding interpolated functionals to treat strong correlation effects. In addition to exploring the performance of such models numerically for the helium and beryllium isoelectronic series and the dissociation of the hydrogen molecule, an approximate analytic model is presented for the initial slope of the local adiabatic connection. Comparisons are made with approaches based on global models, and prospects for future approximations based on the local adiabatic connection are discussed.

Introduction

Kohn–Sham density-functional theory (KS DFT)[1] is the method most widely used in electronic structure calculations, due to its modest computational cost combined with an accuracy that is often competitive with much more expensive ab initio methods. The accuracy of the method is limited by the quality of approximations required to describe the quantum mechanical exchange and correlation (XC) interactions of electrons. A large number of density-functional approximations (DFAs) for the XC energy have been developed in recent decades.

The simplest DFAs are based on the local density approximation (LDA), as proposed by KS in their 1965 paper,[1] in which the XC energy is approximated as a functional of the density at a given point in space. The generalized gradient approximations (GGAs)[2−6] go beyond the LDA by modeling the XC energy as a functional of the local density and its first derivative. The meta-GGAs[7−9] are closely related but their functional forms are also dependent on the KS kinetic energy density and/or, less commonly, the Laplacian of the electron density. Further developments led to the introduction of the occupied KS orbitals as ingredients for the XC energy (hybrid functionals[10−12] and self-interaction corrections,[13−16]), and more recently also the virtual KS orbitals (double-hybrid functionals[17−19] and random-phase approximations[20−22]). Local hybrid functionals[23−26] are also an interesting alternative approach to construct hybrid methods that are pertinent to the context of this work.

The inclusion of additional dependencies in XC functionals has often resulted in significant improvements in their accuracy for general calculations. However, these improvements cannot be described as systematic in the same way that the accuracy of an ab initio calculation may be systematically improved by considering a larger number of excited determinants; some DFAs give excellent results for particular systems but perform very poorly otherwise, and vice versa. There also remain many problems that none of the currently available DFAs can accurately model. An important example of this, which is pertinent to this work, are strong correlation effects, commonly found in systems with near-degenerate orbitals such as the d- and f-block elements, but also in systems where chemical bonds are being broken or formed.

In the present work, the problem of constructing DFAs accurately for systems with and without strong correlation is examined by considering the adiabatic connection (AC)[27,28] at the local level, i.e., at each point of space.[29] The AC, discussed in subsection , provides an exact expression for the exchange and correlation energies by considering the changes that occur as the strength of electron interaction is smoothly increased from zero. This formalism has provided the basis for the development of several DFAs,[10,30−32] which attempt to interpolate the AC between the non-interacting and physical systems in order to estimate the XC energy. An advantage of the AC formalism crucial to our construction is that it allows the problem of strong correlation to be addressed in a more direct way, by creating interpolation models that are explicitly dependent on the strongly interacting limit, in addition to the non-interacting limit, of the AC.

The strongly interacting limit of the AC has recently become the subject of much interest.[29,33−37] The properties of the AC integrand in this limit reveal a highly nonlocal density dependence of correlation effects[33,35,38−40] that cannot be obtained from the standard (semi)local or hybrid functionals. Study of the strongly interacting limit in DFT has focused on the strictly correlated electrons (SCE) functional, where the electrons have an infinite interaction strength. This limit is of particular interest from a theoretical point of view as it can be studied exactly for one-dimensional systems[41] and may be closely approximated in systems with spherical or cylindrical symmetry.[33,42] These studies show that in the limit of infinite interaction strength certain integrals of the density appear in the exchange–correlation functionals, revealing a mathematical structure very different from the one of the usual semilocal or orbital-dependent approximations. The nonlocal radius (NLR) functional proposed in ref (43) approximates the SCE functional with a model that retains some of the SCE nonlocality, introducing integrals of the spherically averaged density around a reference electron. The inclusion of the NLR functional into global and local interpolations along the adiabatic connection has been very recently explored by Zhou et al.[44] In another recent work, Kong and Proynov constructed a functional combining the information from the Becke’13 model[45] and approximating local quantities along the AC.[46]

The aim of the present work is to start a systematic study of local interpolation models along the adiabatic connection, using at a first stage exact input ingredients, thus disentangling the errors due to the interpolation models from those due to the approximate ingredients. The local AC for several closed-shell atoms has been recently computed[47] to high accuracy between the non-interacting and physical systems using the Legendre–Fenchel formulation of DFT due to Lieb,[48] and the Lieb maximization method of refs (49−52). Local information pertaining to the strongly interacting limit is calculated using the SCE functional, and together these quantities are used to calculate local analogues of some established global AC interpolation functionals. We also discuss how to approximate crucial local ingredients such as the initial slope of the local adiabatic curve.

In section , relevant theoretical background is given including an overview of the AC formalism and the construction of DFAs from both global and local variants of the AC. Techniques for computing quantities along the AC are discussed, including the determination of the local AC as introduced in ref (47). In section the construction of a local model for the AC is discussed, considering the non-interacting and strong-interaction limits carefully in this context. The role of the SCE in constructing local AC interpolation models is examined. Finally the forms of some local interpolation models, taken from successful existing global models, are introduced. In section the performance of these models is assessed for the helium and beryllium isoelectronic series and for dissociation of the H2 molecule, a system that typifies the failure of present DFAs to properly account for strong correlation. Directions for future work are outlined in section .

Theoretical Background

Adiabatic Connection

The AC was proposed in a series of reports,[27,28,53,54] which suggested that further insight into electronic correlations in DFT may be gained by considering a system at constant electron density as the interaction strength is smoothly scaled between zero, i.e., the KS auxiliary system, and the full physical interaction strength. This scaling of the interaction strength is achieved by the introduction of a simple coupling-constant coefficient, λ, such that the Hamiltonian for any given λ is written aswhere T̂ is the kinetic energy operator, Ŵ is the physical electron-interaction operator, and V̂λ is the operator representing an external potential vλ that binds the electron density at λ, such that it is always equal to the density of the physically interacting system (ρλ = ρ1, ∀ λ). As the value of λ is smoothly increased from zero to one, the system evolves adiabatically through a family of λ-dependent wave functions Ψλ to the physical system described by Ψ1.

Given a Hamiltonian Ĥλ, one can define the corresponding λ-dependent universal density functional aswhere eq 2b follows from the application of the Hellmann–Feynman theorem to eq 2a. This allows the well-known AC formula to be derived, yielding the following exact expression for the XC energy of an electronic system,where is the (global) AC integrand, given byΨλ[ρ] is the ground-state wave function of Ĥλ in eq , and U[ρ] is the Hartree (Coulomb) energy.

The AC integrand may be characterized by several features that can be exactly defined: the expansion of in the non-interacting limit is given by[55]while its expansion in the strongly interacting limit can be expressed as[33,34,38]Here, the non-interacting terms and are the exchange energy and twice the second-order correlation energy given by Görling–Levy perturbation theory (GL2)[55,56] (see section ), respectively. Their analogues at the strongly interacting limit, and , have been studied in refs (33, 34, and 38) and will be discussed further in section . In addition to these asymptotic limits, the behavior of under uniform coordinate scaling is also well-defined, as discussed in ref (57).

DFAs Based on the Global Adiabatic Connection

To construct practical DFAs one could consider modeling the integrand of eq using a function that interpolates between the limits of eq and eq . The SCE limit is of particular importance in the present work; however, one could also consider models that intercept any other known point on the adiabatic connection for λ > 0. Several attempts to develop DFAs based on these ideas have been put forward in the literature; see, e.g., refs (30, 31, 33, 34, and 58−62). Each form makes a choice of a simple model function and the parameters on which to base the model. These parameters often include the known exact expressions for the parameters in eq , and , since these may be computed from the set of Kohn–Sham orbitals ({ϕ }) and orbital energies ({ϵ }).

The calculation of requires the GL2 correlation energy,[55,56] which leads to a computational cost similar to the second-order Møller–Plesset (MP2) model used in ab initio quantum chemistry.[63] The parameters in the SCE limit entering eq are clearly also of special interest in this context, and they can be computed numerically for atomic systems and molecules with cylindrical symmetry.[33,42,64] More frequently DFAs are derived for points along the AC with λ > 0, often by employing scaling relations to derive forms from existing DFAs. A similar strategy can also be used to approximate by a DFA; see for example ref (60).

In tandem with choosing a set of exact or approximate values to parametrize a model for the AC one must also choose an appropriate model function for the AC integrand. A number of these have been suggested and many have been benchmarked in practical applications. One of the simplest (and most often used) is that of a [1/1] Padé, as suggested by Ernzerhof.[30] A range of forms were suggested by Cohen et al. and tested using approximate parametrizations,[60] leading to the MCY1 functional in which a [1/1] Padé model is employed. Peach et al.[65,66] attempted to disentangle approximations in the choice of parameters from those in the choice of model AC function by utilizing nearly exact KS orbitals and orbital energies derived from full configuration-interaction data to calculate and and the corresponding interacting wave functions to evaluate via eq . Our present study follows a similar philosophy, but applied to local, rather than global, interpolations.

Seidl and co-workers[58,59] were the first to make use of the strong-interaction limit (although approximated at a semilocal level, using the so-called point-charge-plus-continuum, or PC, functional) in constructing a global AC model, known as the interaction strength interpolation (ISI) functional. The revISI model[34] and the models of Liu and Burke[61] were later constructed to take account of the λ–1/2 dependence of the second term of eq , which was not correctly described by the ISI approach. Teale et al. also proposed forms for the AC integrand based on the structure of traditional ab initio methodologies[62] and parametrized these forms to intercept values of the AC at any λ > 0.

The majority of these models for the global AC suffer from the fact they are not size consistent in practice. This deficiency arises from a nonlinear dependence on the parameters , , and a chosen approximation to . When these global parameters enter in a nonlinear fashion (often as ratios), then size consistency is difficult to achieve. One route forward is to construct local AC models, which can replace these global parameters with local values defined at each point in space and may be more amenable to the construction of models that recover size consistency (at least in the usual density-functional sense[67,68]).

Constructing a Local Adiabatic Connection

The AC expression for the XC energy of a system as given in eq describes a global quantity, integrated over the coupling constant λ. However, it may equally be written as the spatial integral of a local quantity analogous to the local value of an XC functional. To this effect, eq may be rewritten aswhere wλ(r) is the energy density at a given λ. It is well-known that the energy density cannot be uniquely defined;[69−71] an arbitrary number of terms may be added to wλ(r), yet an identical recovered if their spatial integral is zero. Thus, any such energy densities are only defined within a particular gauge, and only energy densities defined in the same gauge may be meaningfully compared.

In the context of the present work, it is both convenient and physically meaningful to define wxc,λ(r) in the gauge of the electrostatic potential of the exchange–correlation hole,where hxcλ(r,r′) is the exchange–correlation hole,and P2λ(r,r′) is the pair density obtained from the wave function Ψλ[ρ]The definition of energy densities in the gauge of the XC hole is well-established in the literature, and further discussion may be found in refs (29, 72, and 73). The coupling-constant-averaged (λ-averaged) XC energy density is defined asSince the spatial integral of the product of this quantity and the density yields the XC energy, the same quantity may be considered as a target to be modeled by XC functionals,[47] although GGAs and metaGGAs often aim at energy densities within different definitions.[5,7,72] Given the invariance of the exchange energy density to electron-interaction strength, eq may be trivially resolved into separate exchange and correlation terms asThe aim of the local interpolation schemes examined in this work is to approximate w̅xc(r) and w̅c(r) through interpolating the local AC. In principle, this approach is analogous to that of the global AC interpolation schemes previously discussed, but rather than depending on quantities pertaining to the global AC, they are instead constructed from their local equivalents, wλ(r). Obviously, a local interpolation will only be meaningful if all of the local terms are defined in the same gauge. It is again both convenient and logical to define all local quantities in the gauge of eq , as in which highly accurate energy densities wλ(r) in the range 0 ≤ λ ≤ 1 have previously been calculated,[47] and additionally can be computed for small systems in the limit λ → ∞.[29,37]

At λ = 0, the energy density in the gauge of eq is the exchange energy density w0(r) = w(r), often denoted ϵ(r) in the literature (also equal to 1/2 the nonlocal Slater potential[74]), which is the crucial ingredient of local hybrid functionals. Accurate and efficient computational schemes for this quantity have become available in recent years.[26,75] In a way, local interpolation models can be viewed as local hybrids that carefully address the gauge problem.

The local equivalent of is not as simple to define, yet is an essential component of AC interpolation schemes as it provides a measure of the departure from exchange-only behavior, in other words provides the information from which the correlation energy is approximated. While many global models use GL2 theory for this purpose, its dependence on global quantities makes it unclear how it could be applied to a local interpolation scheme. This is discussed in detail in section .

Lieb Maximization

In order to assess the quality of our local interpolation functionals, it is necessary to have accurate data of energy densities, defined in the gauge of the XC hole. These may be acquired by the method of Lieb maximization, described in refs (50−52).

The Lieb maximization is an optimization algorithm developed using the convex–conjugate functional defined by Lieb in ref (48) as the Legendre–Fenchel transform to the energy,in which the density ρ and potential v are conjugate variables, belonging to the dual vector spacesand Eλ[v] is the energy yielded by a given electronic structure calculation at potential v(r). This convex–conjugate formulation follows from the concavity of variationally determined energy Eλ[v] in v, from which Lieb showed that Fλ[ρ] must be convex in ρ. Furthermore, the conjugate functional to a nonconcave energy, such as that which may result from a nonvariational calculation, remains well-defined as it is necessarily convex. A subsequent Legendre–Fenchel transform of Fλ[ρ] yields the concave envelope (least concave upper bound) to Eλ[v]; hence unique solutions to Fλ[ρ] can always be obtained.

In the Lieb maximization, the optimized density ρ(r) is obtained by maximizing Fλ[ρ] with respect to variations in the potential v(r), rather than by minimizing Eλ[v] with respect to ρ(r) as is usually the case. Therefore, at convergence, the potential v(r) in eq 13b is that which yields ρ(r). In the present work, Lieb maximizations have been carried out at a number of points along the AC in the range 0 ≤ λ ≤ 1; hence the density is constrained such that ρλ = ρλ=1, resulting in a λ-dependent optimizing potential.

In order to effectively optimize with respect to the potential, we parametrize it by using the method of Wu and Yang (WY)[49] aswhere vext(r) is the external potential due to nuclei, vref(r) is a reference potential chosen to ensure that v(r) has the correct asymptotic behavior, and {g} are a set of Gaussian functions with coefficients {b}. In all calculations in this work we choose the potential expansion basis set to be identical to the primary orbital basis set. The reference potential used in this work is the Fermi–Amaldi potential,[76] and Fλ[ρ] is optimized with respect to the coefficients of the potential basis {b}. Additionally, convergence is accelerated through the use of the Newton method described in refs (50−52). The relaxed-Lagrangian formulation of Helgaker and Jørgensen[77] is used to obtain relaxed densities for nonvariational wave functions, which serve as input to the Lieb functional and are used in the determination of the derivatives required for its optimization.

In this work, Lieb maximization calculations are performed using the implementation of refs (50−52) in a development version of the Dalton quantum chemistry software package,[78] in which Eλ[v] is computed by using coupled-cluster singles and doubles (CCSD)[79] and full configuration-interaction (FCI) wave functions. At convergence, where the optimizing potential is such that ρλ = ρλ=1, the relaxed λ-interacting one- and two-particle density matrices are computed, with which the λ-dependent XC energy densities may be obtained as

Modeling the Local AC

Local Slope in the Non-interacting Limit

As described in section , the initial slope of the AC is an important part of many global AC models, in which it may be calculated directly by GL2 perturbation theory; however there is no analogous expression that yields the local equivalent, and we give such an expression in section . Here, the local initial slope of the XC energy density that is given in eq is defined asand is related to the global slope, hence the GL2 correlation energy, by

Numerical Calculation of the Local Slope

In this study w0′(r) is numerically approximated by the method of finite difference, with a series of wλ(r) for λ ≪ 1.

The resulting local slopes in the H2 molecule with bond lengths of 1.4 and 6.0 au are plotted along the H–H bond in Figure , along with the densities from which they are calculated, at the FCI level of theory and in the uncontracted aug-cc-pCVTZ basis set.[80] In both cases, the local slope is greatest in magnitude at the nuclei, as has been seen previously in atoms.[47] It can be seen that the magnitude of the local slope is significantly larger in the stretched H2 molecule, mirroring observations previously made of the global AC in the dissociating hydrogen molecule.[51]

Figure 1

Plots comparing the values of −ρ(r) and w0′(r), with respect to the distance from the bond midpoint, z/a.u., along the principal axis of the H2 molecule with bond lengths of 1.4 au (upper panel) and 6.0 au (lower panel).

The local slopes in the He and Be isoelectronic series are plotted in Figures and 3, respectively. It is clear that, with increasing nuclear charge, the charge densities in both series become increaslingly contracted. The x-axis in both plots has been scaled by nuclear charge, highlighting a contrast in their behavior with respect to the uniform scaling condition,which holds for nondegenerate KS systems.[81] In Figure , it can be seen that the slope of the AC for the He series becomes less negative with increasing Z, tending to an asymptotic value as Z → ∞, consistent with the scaling relation of eq . However, the slope of the AC in the Be isoelectronic series becomes more negative with increasing Z, indicating that the scaling relation is not satisfied by this series.[82]

Figure 2

Plots of w0′(r) for the helium isoelectronic series, with nuclear charges 1 ≤ Z ≤ 10, and with radial distance from the nucleus r/a.u. scaled by nuclear charge.

Figure 3

Plots of w0′(r) for the beryllium isoelectronic series, with nuclear charges 4 ≤ Z ≤ 10, and with radial distance from the nucleus r/a.u. scaled by nuclear charge.

Functional Approximation for the Local Slope

While it is useful to numerically approximate the local slope for the purposes of evaluating local interpolation schemes, such functionals would only be viable for mainstream use in DFT calculations if they can be described by simple functional forms.

In global models, the initial slope can be calculated directly from the occupied and virtual KS orbitals according to GL2 theory,where the indices i, j and a, b pertain to occupied and virtual KS orbitals, respectively, v̂xKS is the local KS potential, and v̂xHF the nonlocal Hartree–Fock (HF) exchange potential. The first term in eq is analogous to the correlation energy given by MP2 theory, in which ϕ and ϵ are canonical HF orbitals and eigenvalues rather than KS ones. The second term accounts for the difference between the KS and HF exchange potentials and has a form similar to a singles term in many-body perturbation theory. Previous studies of GL2 theory have found that the second term, although non-negligible, is small in magnitude relative to the MP2-like term evaluated with the KS orbitals.[83] Therefore, it follows that an approximation to the GL2 correlation energy may be obtained by evaluating the MP2 correlation energy[17] with the KS orbitals and eigenvalues, EGL2 ≈ EMP2.

Given that an approximation to the global AC slope may be obtained from an MP2-like calculation, it follows that an approximation to the local AC slope may be obtained by deriving a local form of this expression. While MP2 theory treats perturbations of the wave function, the analysis may be extended to energy densities in the gauge of the XC hole by means of eq , as the substitution of eq into eq yields the following,where P2′(r, r′) is the derivative of the pair density at λ = 0,Notice that eq ensures that w0′(r) is in the gauge of the electrostatic potential of the XC hole. Given a non-interacting ground-state wave function Ψ(0), the perturbed wave function Ψλ for |λ| ≪ |Ψ(1) – Ψ(0)|2 can be appproximated by the series expansionIf one assumes that Ψ(0) is nondegenerate and has the form of a single Slater determinant, the first-order correction to the wave function is given byRestricting the space of Ψ(0) to doubly-excited determinants reduces this expression toIn MP2 theory, contributions to the correlation energy from singly-excited determinants are necessarily zero due to Brillouin’s theorem. However, this is not strictly true in GL2 theory as the singles term in eq makes a small, but nonzero, contribution to the GL2 correlation energy.[55] As such, considering only double excitations in the model for the local slope can only yield approximations to the local slope; spatial integration of this quantity will not return the exact GL2 correlation energy.

Application of the Slater–Condon rules to eq allows it to be rewritten aswhere the coefficient to Ψ may be identified as an MP2 doubles amplitude t,To obtain P2′(r, r′), it is necessary to take the derivative of the pair density corresponding to the perturbed wave function, Ψλ ≈ Ψ(0) + λΨ(1). Substituting this into eq and rearranging the resulting expressions yields the following,where we assume that Ψ(0) and Ψ(1) are real and P̂2(r, r′) = N(N – 1) ∑δ(r–r) δ(r–r) is the pair-density operator. Substituting eq into this expression giveswhich may then be resolved into the following orbital-explicit expression,where δσ is the Kronecker δ over two spin functions: ∫σ*(m) σ(m) dm = δσ. Substituting eq into eq finally results in an expression for the local slope,where v(r) is the antisymmetrized orbital potential,Multiplying the right-hand side of eq by the density and integrating over all space, we recover twice the MP2-like expression. This is not an exact expression for the local slope, as the second term of eq is not accounted for. However, the omitted term is generally small relative to the MP2-like term and vanishes entirely for two-electron systems; hence the expression for the local slope in eq should, in principle, be a fair approximation of the exact local slope.

In future work we will implement and test eq against the numerical results in section . The doubles amplitudes t are readily obtainable from standard quantum chemical packages, and the potential v(r) can also be readily calculated; however, it would likely be computationally expensive to evaluate on a numerical grid. To reduce this cost, a range of techniques, commonly used to accelerate the calculation of integrals in linear-scaling software packages, may be employed.[26,84−86]

We note that the behavior of the local slopes presented in Figures , 2, and 3 may be rationalized in a manner similar to that commonly discussed for global models in terms of eq . This is because of the key role of the doubles amplitude t in eq . The doubles amplitude has a dependence on the orbitals and orbital energies that is similar to that of the GL2 energy in eq . We see in Figure that the local slope of the hydrogen molecule displays the minima at the nuclei. Equation , which is exact for two-electron systems, can be used to rationalize this observation. For closed shell two-electron systems with only one virtual orbital, eq is simplified as follows:Even if we used a minimal orbital basis for the evaluation of the expression given in eq for the hydrogen molecule, we would see that the local slope is most negative at the two nuclei, for any bond length. While this effect is captured with the minimal basis, the same minimal basis model would incorrectly describe the slope at the bond midpoint. For example, in the top panel of Figure we see that w0′(r) is less than 0 at the bond midpoint of H2 at R = 1.4. Within the minimal basis, the local slope would be exactly 0 for any R, as the antibonding ϕ2(r) orbital which enters eq has a node at the bond midpoint.

We also see in Figure that the correlation energy density for the He isoelectronic series scales quickly toward an asymptotic constant as Z increases. Furthermore, the local slope decays smoothly with distance from the nucleus, reflecting the behavior of v(r). The behavior for the Be isoelectronic series in Figure is more complex. The KS HOMO–LUMO gap is known to increase[87] as Z increases from 4 to 10, from which one would expect the correlation energy to become less negative according to the behavior of t. In the core region this behavior holds; however in the valence region the trend is opposite, with the correlation energy density becoming more negative with increasing Z. This suggests that the numerator of t and the spatial dependence of v(r) due to the form of the KS orbitals are dominant in this region, provided that eq is sufficiently accurate for the Be isoelectronic series.

SCE Model and the Strong Interaction Limit

In recent years, the exact strong-coupling limit of the AC has been intensively studied.[29,33−37] This limit reveals a new structure for the XC functional: instead of the traditional ingredients of DFAs (local density, density gradients, KS kinetic energy density, and occupied and unoccupied KS orbitals), it is observed that certain integrals of the density appear in this limit, encoding highly nonlocal information.[33,35,38−40]

Tests on model physical and chemical systems (electrons confined in low-dimensional geometries and low-density, ultracold dipolar systems, simple stretched bonds and anions) have shown[35,37,39,40,88−90] that taking into account this exact behavior can pave the way for the solution of the strong correlation problem in DFT. However, the exact information encoded in the infinite coupling limit, described by the SCE functional, does not come for free: the SCE problem is ultra-nonlocal, and, although sparse in principle, its nonlinearity makes its exact evaluation for general three-dimensional geometry a complex task. A possible route to find suitable algorithms relies on the fact that constructing the exact SCE functional for a given density is equivalent to solving an optimal transport (or mass transportation theory) problem with a cost function given by the Coulomb interaction.[91,92] This equivalence has triggered interest from the applied mathematics community working on optimal transport problems, which has led to the suggestion of several algorithms,[89,93−95] together with very interesting exact results.[96−98]

So far, the SCE solution is known exactly for one-dimensional systems.[41] For spherically symmetric systems, a conjectured solution[33] that is very close to the exact one[64] (and it is in many cases, but not always,[98] exact) has been proposed and used to address interesting physical problems.[90,99] Using algorithms and ideas from optimal transport, the SCE problem for the hydrogen molecule along the dissociation curve has just recently been solved and both the global[37,89] and local[37] SCE quantities have been computed. A more practical way to proceed is to build approximations for the SCE functional inspired by its exact form, as it was done in the construction of the NLR functional.[43,44]

The SCE system complements the KS system.[33,34,38] It corresponds to the wave function that minimizes the Hamiltonian of eq when λ → ∞. One can argue that the SCE system is a better starting point than the Kohn–Sham system for the description of very strongly correlated systems.[37,39,40,89]

The SCE functional is defined as[29,33,35]The XC part can be easily extracted from , as . The KS SCE approximation, proposed in ref (35), uses the SCE functional to approximate the Hartree and exchange–correlation energy, and it is equivalent to setting for all λ. It has been shown that KS SCE yields good energies for systems where correlation plays a dominant role, such as electrons confined in low-density nanodevices or extremely stretched bonds.[37,39,40,89] On the other hand, KS SCE treats moderately and weakly correlated systems very poorly, giving energies that are unacceptably too low.[37,88] A less drastic approximation is to construct a model, in such a way that its λ → ∞ limit is given by the exact or approximate value of , as done in the pioneering work of Seidl et al.[58,59] Analogously, one can also model wλ(r), imposing that its λ → ∞ limit is given by w∞(r). This latter approach is the main object of the following sections.

In the SCE limit, the electrons are infinitely or perfectly correlated and their positions are given by an infinite superposition of classical configurations. The basic idea is that the electronic positions are all determined by a collective variable r, a feature that is captured by the so-called comotion functionsf(r).[33] If a reference electron is at r, then the position of all of the other electrons in the system will be given by r = f(r).[33] Since the electrons are perfectly correlated, the probability of finding the reference electron at r has to be the same as the probability of finding the ith electron at f(r). Therefore, the comotion functions have to satisfy the following differential equation:[33]For more details on the comotion functions, including their group properties, see refs (29, 91, 99), and (33).

In terms of the comotion functions, the SCE functional is given by[29]Despite the high nonlocality of the SCE functional, evident from eq , we can easily compute its functional derivative from the following expression[35,91]Equation suggests the following energy density in the SCE limit:where vH(r) is the Hartree potential. This expression is indeed in the gauge of the XC hole potential of eq , as proven in ref (29). Being derived from a wave function, the w∞(r) energy density decays like , similar to the physical (λ = 1) and the exchange (λ = 0) energy densities of eq . Its functional derivative, eq , has also the correct asymptotic behavior.

To solve the SCE problem for spherically symmetric systems (the He and Be isoelectronic series considered in this work), we have used the approach presented in ref (33), which is exact if N = 2. For atomic densities with N > 2 it could be either a very good approximation for the true minimum of eq or, again, the exact result.[64,98] For the H2 molecule we have used the results of ref (37), where the SCE energy density has been computed by obtaining the comotion function from the dual Kantorovich formulation[91,100] of the SCE problem.

Local Interpolation Models

The local interpolation models tested in this work are largely simple translations of the well-established global interpolation models into a local form. This was done for the model of Seidl, Perdew, and Levy (SPL),[58] the “simplified” model of Liu and Burke,[61] which will be referred to here as the LB model, and the Padé[1/1] model.[30,101] Each of the energy densities resulting from the three mentioned models is constructed from three local parameters, a, b, and c, which are defined in the gauge of the XC hole. The functional forms of these three models are summarized in Table .

Table 1

Mathematical Forms of the Local AC Interpolation Models (for the Padé[1/1] Model, p > 0, )

	wλ(r)	a(r)	b(r)	c(r)	refs
SPL		w∞(r)	w0(r) – w∞(r)		(33, 58)
LB		w∞(r)	(w0(r) – w∞(r))/2		(61)
Padé[1/1]		w0(r)	w0′(r)		(30, 101)

In addition to these, we constructed a local form of the two-legged representation[31] which, given some value of w1(r), takes the form

Whenever we used the two-legged representation to model the local AC in this work, we did it by incorporating the interpolated w1(r) of the LB model: w1(r) ≈ w1LB(r). By doing the local interpolation this way, we use the following three input quantities: w0(r), w0′(r), and w∞(r) and circumvent the direct utilization of the full interacting energy density, w1(r). In each of these four models, integration of wλ(r) with respect to coupling constant gives the λ-averaged energy density w̅xc(r), which, if spatially integrated according to eq , yields the XC energy Exc[ρ].

An important observation in the translation of global to local models is that while the following global inequalities are always satisfied,their local counterparts do not necessarily satisfy these same inequalities. It has previously been observed for the Hooke’s atom series that, in the tail regions of the density, w∞(r) can be less negative than w1(r).[29] In this work, the crossing of w∞(r) with w̅xc(r), w1(r), and w0(r) has only been observed in the tail regions of the density and is thought to be an artifact of the numerical instability that occurs where the density is very small.

Results

Helium Isoelectronic Series

Although the helium isoelectronic series is a set of only two-electron systems, it is a useful series to consider in evaluating the local interpolation models as most standard DFAs incorrectly characterize the hydride ion (H–), failing to predict its existence as a bound electronic system.[66,102] Here, local interpolation models are constructed from energy densities acquired by the Lieb maximization at the FCI level, as described in section , in the range 0 ≤ λ ≤ 1 and at λ = ∞ by evaluating the SCE functional on the λ = 1 density, also at the FCI level of theory.

In Table , the correlation energies given by local forms of the SPL, LB, two-legged representation (the column labeled “two-leg”) and Padé[1/1] models (the latter parametrized using the accurate values for w1(r), in order to compare with models that, instead, use the λ → ∞ information) are given, along with that given by the global SPL model and the FCI correlation energy for comparison. These data show that the local interpolation correlation energies are in close agreement with the FCI reference values; the mean absolute errors (MAE) of the local interpolation models are 2.0, 1.5, 0.5, and 0.1 mH, for the two-legged representation, SPL, LB, and Padé[1/1] models, respectively.

Table 2

Reference and Interpolated Ec Values, in hartree, for the He Isoelectronic Series

Z	FCI	local SPL	global SPL	local LB	Padé[1/1]	local two-leg
1	–0.0409	–0.0367	–0.0368	–0.0398	–0.0401	–0.0477
2	–0.0400	–0.0378	–0.0380	–0.0394	–0.0399	–0.0435
3	–0.0410	–0.0393	–0.0395	–0.0404	–0.0409	–0.0431
4	–0.0416	–0.0402	–0.0404	–0.0411	–0.0415	–0.0433
5	–0.0418	–0.0408	–0.0409	–0.0415	–0.0418	–0.0433
6	–0.0419	–0.0410	–0.0411	–0.0416	–0.0418	–0.0431
7	–0.0414	–0.0407	–0.0408	–0.0412	–0.0414	–0.0423
8	–0.0412	–0.0405	–0.0406	–0.0410	–0.0412	–0.0420
9	–0.0411	–0.0405	–0.0406	–0.0409	–0.0411	–0.0419
10	–0.0411	–0.0405	–0.0407	–0.0408	–0.0411	–0.0418

As would be expected, the local Padé[1/1] is the most accurate of the models, given that it is derived from the full interacting energy density. These data further suggest that the local LB model is marginally superior to the local SPL and the two-legged representation. However, comparing the global and local models shows a slightly lower error for the global model; the local SPL model has an MAE of 1.5 mH, compared to 1.3 mH for the global model.

Figure compares the FCI w̅c(r) with that of the local LB and SPL models, for the helium atom. This reflects the numerical data in Table , both being very close to the FCI energy density but with slightly lower error in the LB model.

Figure 4

Plots comparing the FCI, local LB, and local SPL λ-averaged correlation energy density in the helium atom.

Beryllium Isoelectronic Series

The changes in correlation energy across the beryllium isoelectronic series are somewhat more complicated than those in the helium isoelectronic series, and its explanation involves the interplay of several effects. With increasing nuclear charge, the density becomes increasingly contracted, suggesting that the correlation energy should approach the high-density limit for very large Z. However, this is accompanied by a changing KS HOMO–LUMO gap, here the energy difference between 2s and 2p orbitals, which increases from Z = 4 → 13 before decreasing with higher Z values.[87]

Table shows the reference and interpolated E results for the Be isoelectronic series, with Z in the range 4–10. The wave function for λ values between 0 and 1 has been computed in the same way as that for the He isoelectronic series, however at the CCSD level of theory rather than FCI. As for the helium series, the local Padé[1/1] that uses w1(r) is the most accurate of the local interpolation models. However, in contrast to the findings for He isoelectronic series, the local interpolation models are much more accurate than the global models. For example, in the case of F5+ the global SPL model has an MAE of 47.4 mH, whereas the error for the local SPL model is 3.5 mH. The local two-legged interpolation underestimates the correlation energies of the elements of the given series. We discuss in more detail this model in the next subsection.

Table 3

Reference and Interpolated Ec Values, in hartree, for the Be Isoelectronic Series

Z	CCSD	local SPL	global SPL	local LB	Padé[1/1]	local two-leg
4	–0.0920	–0.0876	–0.1049	–0.0925	–0.0911	–0.1046
5	–0.1089	–0.1041	–0.1250	–0.1100	–0.1076	–0.1246
6	–0.1244	–0.1202	–0.1455	–0.1271	–0.1229	–0.1444
7	–0.1389	–0.1363	–0.1668	–0.1443	–0.1373	–0.1645
8	–0.1534	–0.1532	–0.1898	–0.1626	–0.1517	–0.1859
9	–0.1683	–0.1717	–0.2157	–0.1826	–0.1666	–0.2098
10	–0.1833	–0.1920	–0.2447	–0.2046	–0.1817	–0.2361

Figure shows the λ-averaged correlation energy densities for the beryllium atom. The shape of w̅c(r) reflects the shell structure of the Be atom.[47,82] The local SPL and LB interpolation models appear to qualitatively capture the shell structure of w̅c(r); however, in some regions it overestimates the reference value while in other regions the converse is the case. The error cancellation that results from this is the most likely explanation for the superior accuracy of the local models in comparison to the global models.

Figure 5

Plots comparing the CCSD, local LB, and local SPL λ-averaged correlation energy density in the beryllium atom.

Hydrogen Molecule

Despite the development of DFT into the most widely applied electronic structure method, and the wealth of XC DFAs that have been developed, there are some systems for which no combination of DFAs provide an accurate description. A well-known example of such a system is the dissociating H2 molecule.[42,103] Standard DFAs become increasingly inaccurate with greater H–H bond length, reflecting a fundamental flaw of DFAs in their inability to properly treat strong correlation.

It has been seen previously[37,89] that KS SCE correctly predicts the dissociation of H2 in a spin-restricted formalism; however at equilibrium geometry the energies it predicts are extremely low and the bond lengths predicted are overly short. The overall accuracy of KS-SCE for H2 dissociation can be substantially improved by the addition of nonlocal corrections.[37]

Figure shows the dissociation curves for H2 given by the local interpolation models, along with those given by HF, FCI, and the PBE functional[5] for comparison. The computational details are the same as those of the He isoelectronic series, and the PBE, HF, and FCI curves have been obtained from the Dalton quantum chemistry package[78] all within the uncontracted aug-cc-pCVTZ basis set.[80] The SCE energy density has been computed by using the dual Kantorovich method.[37]

Figure 6

Potential energy curves for the H2 molecule with internuclear distance R/a.u., which are obtained using the local interpolation methods: SPL, Liu–Burke, two-legged representation combined with the Liu–Burke model, Padé[1/1] with w1(r). Restricted HF, PBE, and FCI curves are also shown for comparison.

It can be seen in Figure that all of the interpolation models correctly predict the dissociation of H2, which follows from their inclusion of w∞(r). In global AC models, at infinite separation the initial slope diverges as a result of the vanishing HOMO–LUMO gap and the SPL and LB models reduce to , yielding the exact energies. However, the dissociation curves produced by the local models approach the FCI curve slowly, resulting in an unphysical “bump”-like feature. This is a well-known failing of DFT, having been observed with other functionals, such as the random-phase approximation[103] and even the global Padé[1/1] model with .[65] It can be seen in Figure that this is not remedied by the local interpolation approach, as the curve obtained by the local Padé[1/1] also exhibits this unphysical bump, as does that given by the local SPL model and, to a lesser extent, the local LB model.

To analyze why the intermediate region is less accurately described by the local interpolation methods than the equilibrium and stretched region, we show in Figure the correlation component of the local AC at one of the nuclei of the hydrogen molecule at different bond lengths: R = 1.4 au (at equilibrium), R = 5.0 au (the intermediate region), and R = 13.0 au (stretched bond). The structure of the three local AC curves at one of the nuclei is very similar to the structure of the corresponding global AC curves.[50] From the given figure we see that at equilibrium the local AC is almost linear, so we can expect that even a single line segment approximation to the local AC: wλ(r) ≈ w0(r) + λw0′(r) would properly capture the shape of the given local AC curve. The local AC curve at the nuclei of the stretched H2 exhibits the characteristic “L-shape”, which was also observed in the case of the corresponding global AC curve.[50] We would expect that the two-legged representation would capture the given local AC very well, but even a single line segment approximation, wλ(r) ≈ w∞(r), this time coming from the strong coupling limit, would be highly accurate for the stretched H2.[37] In contrast to the local AC curves of the stretched and H2 at equilibrium, the curvature of the local AC curve at the intermediate bond length is highly pronounced. The shapes of the local AC curves at the nuclei mirror the difference in correlation regimes present in the hydrogen molecule at different bond lengths. While in H2 at equilibrium and at very stretched bond length, correlation is almost purely dynamical and almost purely static, respectively, in the intermediate dissociation region there is a subtle interplay between the dynamical and static correlation.

Figure 7

FCI local correlation AC curves at one of the nuclei of H2 for different internuclear separations, R.

In the intermediate region of the dissociation curve, where the unphysical bump is present, the local two-legged representation model is more accurate than the local Padé[1/1] which we always use here with w1(r). This may be understood by comparing the exact local AC data with the interpolated quantities. The top panel of Figure shows the difference between w̅FCI(r) and that of each of the local interpolation models, along the H–H bond at the 5.0 au geometry, as a function of the distance from the bond midpoint z. This difference δw(r) = w̅FCI(r) – w̅model(r), is multiplied by the density to represent an energy per volume element. It shows that the local SPL energy density is the one that most overestimates w̅(r). The error is smaller for the LB model and even more so for the local Padé[1/1] model. The error is smallest in the two-legged model, obtained using the w1(r) of the local LB. Furthermore, there is the error cancellation in the two-legged model, as there are regions where the w̅(r) of this model underestimates w̅FCI(r).

Figure 8

Plots of the difference between FCI and interpolated λ-averaged energy densities, δw(z) = w̅FCI(z) – w̅model(z), with respect to the distance from the bond midpoint, z/a.u. (upper panel), and the local AC curves at one of the nuclei of the FCI and local interpolation models (lower panel), both in H2 with a 5.0 au bond length.

It can also be seen that the curves shown in the top panel of Figure have a maximum at the nucleus (z = 2.5 a.u.). Focusing on this region, it appears that the FCI curve meets that of the Padé[1/1] at λ = 1 and that the two-legged representation curve meets that of the LB model also at λ = 1. This follows from the construction of the Padé[1/1] and two-legged curves from w1FCI(r) and w1LB(r), respectively. All curves, except for that of the two-legged model, lie above the FCI curve. In the case of the two-legged interpolation model, the first line segment is below the FCI curve, as a result of eq 39a and the convexity of the given local AC curve. The second line segment that starts at xλ ∼ 0.1 is given by w1LB(r) and lies above the FCI curve. The resulting error cancellation makes it clear why the two-legged representation appears more accurate than the other models.

Conclusion and Perspectives

In this work we have studied local interpolations along the adiabatic connection for the He and Be isoelectronic series and the hydrogen molecule, by using accurate input local quantities computed in the gauge of the electrostatic potential of the XC hole and comparing the results with nearly exact energy densities defined in the same way. In order to obtain approximations to the local AC over the physical regime (0 ≤ λ ≤ 1), we constructed interpolation models between the weak and strong coupling limits of DFT. The weak coupling energy densities were obtained using the Lieb variation principle, while the strong coupling limit energy densities were obtained using the strictly correlated electrons (SCE) approach. The inclusion of the SCE information in density-functional approximations helps to ensure their ability to capture the strong correlation effects.

Unlike previous attempts in this direction that used global (integrated over all space) input quantities to model the AC, the local approach is more amenable to the construction of approximations that do not violate size consistency, at least in the usual DFT sense.[67,68] Since the aim here is to work in a restricted formalism, avoiding mimicking strong correlation with symmetry breaking, some care must be taken when discussing size consistency. In fact, strictly speaking, in a restricted framework the energy densities of the second-order perturbation theory and exact exchange are not intensive quantities in the presence of near degeneracy,[67,68] which is the main challenge of capturing strong correlation within DFT.[104−106]

In future work we will test different approximations for the SCE energy densities and the local slope. The development of algorithms for solution of the SCE problem is a very active research field. In spite of the recent improvements, we still lack an algorithm that will solve the SCE problem for general 3D molecular geometries at low computational cost. However, a good candidate to approximate the SCE energy density in the gauge of the XC hole potential is the nonlocal radius functional (NLR),[43] which has been already implemented and used in ref (44).

In addition to numerically exploring the local AC we have also reported the local weak-coupling slope of the adiabatic connection and derived an approximate expression for it in terms of occupied and unoccupied orbitals. This quantity is very important to signal the amount of correlation at each point of space. In our future work we will implement this expression and test it against the results reported here.