Local and global interpolations along the adiabatic connection of DFT: a study at different correlation regimes

Interpolating the exchange–correlation energy along the density-fixed adiabatic connection of density functional theory is a promising way to build approximations that are not biased toward the weakly correlated regime. These interpolations can be performed at the global (integrated over all spaces) or at the local level, using energy densities. Many features of the relevant energy densities as well as several different ways to construct these interpolations, including comparisons between global and local variants, are investigated here for the analytically solvable Hooke’s atom series, which allows for an exploration of different correlation regimes. We also analyze different ways to define the correlation kinetic energy density, focusing on the peak in the kinetic correlation potential.


Introduction
The density-fixed adiabatic connection [1] of Kohn-Sham (KS) density functional theory (DFT) is a powerful theoretical tool for the construction of approximate exchange-correlation (XC) functionals: for example, hybrid [2] and double-hybrid functionals [3] can be constructed from simple models of the adiabatic connection integrand [4][5][6]. These approximations, however, use exact ingredients only for the limit of small coupling strength and are thus biased toward the weakly correlated regime.
A class of approximations that removes this bias is based on the idea of Seidl et al. [7][8][9] to interpolate the adiabatic connection integrand between its weak-and strong-interaction limits. This way, information from both extreme correlation regimes is taken into account on a similar footing. These interpolations can be performed on the global [7][8][9][10][11][12] (i.e., integrated over all spaces) ingredients, or in each point of space, using energy densities [13][14][15]. As well known, energy densities are not uniquely defined and one should be sure, when doing an interpolation between weak coupling and strong coupling in each point of space, that all the input local quantities are defined in the same way [13][14][15][16], which makes the use of semilocal approximations very difficult, a problem shared with local hybrids [17][18][19][20]. Non-local functionals for the strong-interaction limit [21,22] or the physical regime [23] are needed in this context, as full compatibility with the exact exchange energy density is required.
Interpolations constructed from the global ingredients are in general computationally cheaper than their local counterpart, not only because they can use semilocal approximations for the strong-interaction functionals, but also because they do not need energy densities from exact exchange and from secondorder perturbation theory, but only their global values. These global interpolations are in principle not size consistent, but it has been recently shown that their size consistency error can be fully corrected at no additional computational cost [12], allowing for the calculation of meaningful interaction energies This paper is dedicated to János Ángyán: his search for new understanding and passion for research continues to be a great source of inspiration.
Published as part of the special collection of articles In Memoriam of János Ángyán.
with T the electronic kinetic energy operator, Ŵ the Coulomb electron-electron interaction operator, and " → " indicating all fermionic wavefunctions yielding the oneelectron density ( ) , one obtains an exact formula [1] for the XC energy functional of KS DFT, where [ ] is the minimizing wavefunction in Eq. (1) and U[ ] is the Hartree repulsion energy. The real parameter is a knob that controls the interaction strength, defining an infinite set of systems all with the same one-electron density ( ) = =1 ( ) , but with different correlation. The global adiabatic connection integrand has the known expansions at small and large ,

Energy densities
Equation (2) can also be written in terms of real-space energy densities w ( ;[ ]), which are, of course, not uniquely defined. For the purpose of building -interpolation models on energy densities, the choice of the gauge of the electrostatic potential of the exchange-correlation hole h xc ( 1 , 2 ) seems so far to be the most suitable [16], where h xc ( 1 , 2 ) is defined in terms of the pair density P 2 ( 1 , 2 ) and the density (see also [34]), Energy density at = 0 . At = 0 , we have the Kohn-Sham or exchange hole, which yields in the case of a closed-shell singlet considered in this work (with real orbitals) where i ( ) are the occupied KS spatial orbitals.
Slope of the energy density at = 0 . The slope w � 0 ( ) of the energy density at = 0 in the gauge of Eq. (7) is given, again for a closed-shell singlet with real orbitals, by [14] (6) where a and b are unoccupied and i and j are occupied Kohn-Sham orbitals, ⟨ij�ab⟩ denotes the Coulomb integral over the spatial orbitals, and the i are the Kohn-Sham orbital energies. For systems with N > 2 , there should be also a term with single excitations [31], which we do not consider here as we focus on N = 2. Energy density at = ∞ . In the → ∞ limit we obtain a system of strictly correlated electrons (SCE), for which it has been shown [13] that where v H ( ) is the Hartree potential and i ( ) are co-motion functions that determine the position of the ith electron given the position of a chosen reference electron (as the i ( ) satisfy cyclic group properties it does not matter which electron is chosen as reference), and are non-local functionals of the density ( ) [32,35].
There is at present no local expression in the gauge of Eq. (7) for the next leading term W � ∞ [ ] in the → ∞ asymptotic expansion. In fact, the functional W � ∞ [ ] can be computed from an integral on position-dependent zero-point energies [33], which, however, do not provide an energy density within the definition of Eq. (7).

Global and local interpolations
The original idea of Seidl et al. [7][8][9] was to build an approximate adiabatic connection integrand W ISI [ ] by interpolating between the two limits of Eqs. (4) and (5). These interaction strength interpolation (ISI) functionals typically use as input the four ingredients (or a subset thereof) appearing in Eqs. (4) and (5): , denoted in short. The XC energy functional E ISI xc [ ] is then obtained from Eq. (2), by integrating W ISI [ ] over , which will result in a nonlinear function of the input ingredients . Because of this nonlinear dependence, the ISI-type functionals are not size consistent when a system dissociates into unequal fragments, even when the input ingredients are size consistent themselves. However, in this latter case, size consistency can be easily restored with a very simple correction [12]. The ISI-type functionals are, instead, automatically size extensive [12]. Several formulas for interpolating between the two limits of Eqs. (4) and (5) have been proposed in the literature and are reported in "Appendix".
More recently, these same interpolation formulas have been used to build, in each point of space, a model energy density w ISI ( ;[ ]) , with Eqs. (10)-(12) as input ingredients [14,15]. This way, by integrating w ISI ( ;[ ]) over between 0 and 1, one obtains an exchange-correlation energy density in the gauge of the coupling constant averaged

Hooke's atom series
The Hooke's atom series consists of two electrons bound by an harmonic external potential, with hamiltonian with r i = | i | and r 12 = | 1 − 2 | . At large the system has high density and is in the weakly correlated regime, which can be fully described by using the scaled coordinates i ≡ √ i , while as → 0 the system becomes more and more correlated [25], and the relevant scaled variables are ̃ i ≡ 2∕3 i . As well known, there is an infinite set of special values of for which the hamiltonian (13) is analytically solvable [24] once rewritten in terms of center of mass and relative coordinates. These analytic solutions have the center of mass in the ground state of an harmonic oscillator with mass m = 2 and frequency √ 2 , and the relative coordinate in an s-wave with the radial part described by a gaussian times a polynomial [24]. We denote here the various analytic solutions with the degree n − 1 of the polynomial in r 12 . At n = 1 we have the non-interacting system, and as n increases the system becomes more and more correlated, with smaller and smaller [24]. The values of corresponding to the different values of n considered here are reported in Table 1.

Computation of exact energy densities
Given the analytic solutions [24] (r 1 , r 2 , r 12 ) of the hamiltonian (13) for n = 2, … , 6 , we have computed the corresponding densities (r) , which are also analytic. Although leading to cumbersome expressions, these densities allowed us to obtain analytic Kohn-Sham potentials v s (r) = between the physical states with two and one electrons.  We have obtained these energy densities analytically from the exact densities. They are shown in Fig. 1 for the different analytic solutions considered here.

Energy densities for the slope at = 0
The analytic exact Kohn-Sham potentials were used to obtain the virtual Kohn-Sham orbitals needed for the evaluation of Eq. (11). We used an isotropic spherical Gaussian basis with as the width parameter. Angular momentum values were included from l = 0 to l = 9 , with 5-30 basis states for every value of l. All matrix elements were obtained analytically in this basis, including the Coulomb integrals.
We first analyze the convergence of the global slope of the coupling constant integrand, W � 0 = 2 E GL2 c , with increasing basis set size n basis in the first panel in Fig. 2. The number of basis states is that per angular momentum quantum number, with all l up to l = 9 included. As decreases (the quantum number n increases), the l = 0 contribution becomes less important, with the l > 0 contributions gaining more weight, as shown in the second panel in Fig. 2, where the result from each channel l with n basis = 30 is reported.
For the local slope w � 0 ( ) , only 10 basis states are used. In the present case of a two-electron system, w � 0 ( ) can also be simplified, as there is only one occupied Kohn-Sham spatial orbital. Additional utilization of the spherical symmetry then yields the following expression, by using the spherical harmonic expansion of the Coulomb potential, where the functions R l n (r) are the radial functions of the spatial orbitals and occ is the occupied Kohn-Sham orbital.  Table 1.
In the second panel, the energy density has been multiplied by the density and by the volume element. The high-density scaling has been used The full local slope is shown in the first panel in Fig. 3. Numerical issues appear at around the scaled variable values s ≳ 4.5 , but this is of no relevance to the integrated energy as it is clear upon multiplication by the volume element and the density (second panel in Fig. 3).

Energy densities at = ∞
The energy density w ∞ ( ) of Eq. (12) in the case of N = 2 electrons in a spherical density is known to be determined by the radial co-motion function f(r), which gives the full ( ) via ( ) = − f (r) r [13,32,38,39], yielding In turn, f(r) is a fully non-local functional of the density (r) , given in terms of the cumulant N e (r) of Eq. (15) and .
In Fig. 4, we report the energy densities w ∞ (r) for the analytical solutions corresponding to the values in Table 1.

Energy densities at = 1
Since we have exact analytic wavefunctions we can also compute the exact energy densities at physical coupling strength = 1 , which can be used to test the accuracy of local interpolations between = 0 and = ∞ , as well to study features of the energy densities as the interaction strength is changed. The exact w 1 (r) are reported in Fig. 5. We see that the physical energy densities w 1 (r) for the Hooke's atom series differ more among each other at large r, unlike w 0 (r) and w ∞ (r) . This is clearer if we look at the correlation energy density w c (r) = w 1 (r) − w 0 (r) , which is reported in Fig. 6. The correlation energy density w c (r) decays ∝ − 1 r 3 , but with different coefficients for different values of .
A comparison of the three energy densities w 0 , w 1 and w ∞ is given in Fig. 7 for the Hooke's atom with n = 6 . An interesting  Energy densities corresponding to = 1 (first panel), and energy densities corresponding to = 1 multiplied by the density and the volume element (second panel). The coordinates and energy densities are scaled according to the large limit feature of these energy densities, already observed in Ref. [13], is that for large r it can be seen that w 1 (r) < w ∞ (r) , while for the corresponding global quantities we have the strict ine- However, taking w 1 (r) ≈ w ∞ (r) for large r only has a small effect on the energy even for the most strongly correlated Hooke's atom considered here ( n = 6 ), as it becomes clear once the energy densities are multiplied by the density and the volume element (second panel in Fig. 7), which is what ultimately determines the correlation energy. This crossing of energy densities has never been observed, so far, in systems with the Coulomb external potential. , which in this case is given by [33] with
We have used the interpolation formulas reported in "Appendix", namely SPL [7], LB [40], ISI [9] and revISI [33]. The first two, SPL and LB, use only three ingredients (they do not include W � ∞ [ ] ), while ISI and revISI use all the four ingredients of Eqs. (4)- (5). Additionally, we have also used a Padé approximant (see "Appendix") which uses W 0 [ ], W � 0 [ ] and the exact W 1 [ ] , to generate plausible reference adiabatic connection curves, which are shown in Fig. 8. As expected, as the Hooke's atoms get more correlated, the AC integrand displays a stronger curvature.
The error resulting in the correlation energy E c [ ] with the different global interpolations is shown in Fig. 9. We consider only the correlation energy, since all the methods utilize 100% exact exchange. The Padé method performs best as expected, since it uses the exact W 1 , which in practical situations is unavailable. The LB interpolation formula performs second best, while SPL, containing the same ingredients, performs much worse. The ISI and revISI methods improve slightly the SPL formula, but are still outperformed by LB, despite containing more exact information in the form of W � ∞ [ ]. For comparison with traditional Density Functional Approximations (DFAs), such as the local density approximation (LDA) [41] and the PBE GGA [42], we show the error in the exchange-correlation energy E xc [ ] in the Fig. 6 Correlation energy densities (first panel) and correlation energy densities multiplied by the density and the volume element (second panel). The coordinate and energy density are scaled according to the large limit Fig. 7 Energy densities for the most strongly correlated Hooke's atom considered here ( n = 6 ), at different values of (first panel). In the second panel the energy densities have been multiplied by the density and the volume element first panel in Fig. 10. It is clear that the adiabatic connection interpolation methods outperform the PBE method, however at the increased computational cost of a double hybrid. In the second panel in Fig. 10 if we include also the second term. The LDA performs poorly already for the first Hooke's atom and its performance worsens as correlation increases. The GL2 method works well for the first Hooke's atom, which is expected since its adiabatic connection integrand resembles a straight line in Fig. 8, but it is way too negative for the exchange-correlation energy in the more correlated Hooke's atoms. The → ∞ expansion alone performs better as the Hooke's atoms become more correlated, but with the first term only is still too negative by about 15% in the strongest correlated Hooke's atom. Adding the second term contribution reduces the error for n > 3 , and the resulting XC energy becomes now less negative than the exact one.

Interpolations on energy densities
As already mentioned at the end of Sec. 2, an expression for the energy density corresponding to W � ∞ [ ] in the gauge of Eq. (7) is not available. For this reason, we can only test local interpolations using the LB and SPL interpolation formulas, which do not use the information from W � ∞ [ ] . We first compare the resulting w c (r) = w 1 (r) − w 0 (r) from the two interpolation formulas in the first panels in Figs. 11 (LB) and 12 (SPL) with the exact result obtained from the analytic wavefunctions. The errors are small on an absolute scale, so we show in both figures w c ( ) = w c, exact ( ) − w c, model ( ) and include the volume element and density. Notice that w c ( ) = w 1 ( ) since we use the exact w 0 ( ) in the construction of both the LB and SPL approximations. In order to assess the coupling constant integrated energy density w c , which is not known exactly for any of the Hooke's atoms, we compare it with the one obtained from the Padé interpolation, which includes the exact w 0 (r) , w � 0 (r) and w 1 (r). We see that in the case of LB there is an overestimation of the coupling constant averaged energy density at small r, which cancels quite well with an underestimation at large r, achieving almost perfect error cancelation. In the case of SPL, there is a smaller overestimation of the correlation at small r, coupled with a stronger underestimation of the correlation energy at large r, which worsens its performance.

Comparison between global and local interpolations
Of interest is then comparing the performance of the global and local variants of the Padé, LB and SPL interpolations. In Fig. 13 the relative error on the correlation energy obtained from the local and global interpolation is shown, where in this case we use for both 10 basis states per angular momentum quantum number for the slope. In the case of the Padé interpolation the performance worsens only slightly going from the global to the local interpolation, while for the SPL interpolation there is a dramatic worsening. In the case of the LB interpolation the error switches sign for n ≥ 3 and in general worsens. This is somehow surprising as, instead, for small chemical systems the local interpolations have been found to be of similar quality (for N = 2 ) or to significantly outperform (for N > 2 ) their global counterparts [14,15].

Kinetic correlation energy densities
The coupling constant integration is one possible way to recover the correlation part due to the difference between the true, interacting, kinetic energy Another kinetic correlation energy density that has been defined [27] and studied [43][44][45] in the literature, and that has been found to display very interesting features for strongly correlated systems [28,29,46,47], arises from the work of Baerends and coworkers [27,[43][44][45], where (2, ..., N| ) is a conditional amplitude defined in terms of a wavefunction and its density , 1, ...N denote the spatial and spin coordinates of the N electrons, and in Eq. (23) we consider the conditional amplitude from the exact wavefunction (denoted with ) and for the KS determinant (denoted with s ). Equation (23) can also be rewritten in several different interesting and more practical forms, for example in terms of first order density matrices, or in terms of natural orbitals, or with Dyson orbitals (see, e.g., [43][44][45][48][49][50][51][52][53]). In the present case of N = 2 electrons, Eq. (23) takes the simple form where ( 1 , 2 ) is the exact ground state wavefunction of the interacting system.
Both w( ) − w 1 ( ) and v c,kin ( ) integrate to T c [ ] when multiplied by the density ( ) , but they describe the kinetic correlation energy locally in a different way. Here we compare the features of these two definitions, as the kinetic correlation energy is important to capture strong correlation. Also, very recently, it has been proposed to use the correlated kinetic energy density as an additional variable in an extended KS DFT theory for lattice hamiltonians [54], and it is thus important to understand which definition is the most suitable to generalize this theory to the continuum.
In Fig. 14, we show the two different kinetic correlation energy densities, where for w(r) we have used the integration over of the Padé model, which uses the exact w 0 , w ′ 0 and w 1 as input. We see that the two are rather different: v c,kin (r) displays a peak in the center of the harmonic trap, reminiscent of the one appearing in a stretched bond [27][28][29]47], while w( ) − w 1 ( ) displays a weaker peak, which is not located at the center.

Analysis of the peak of v c,kin ( )
In the case of a stretched bond, it has been shown that the height of the peak of v c,kin ( ) at the midbond saturates as the bond is stretched [28], displaying an anomalous scaling [29], which is the way in which exact KS DFT can describe Mott insulator physics [29], and which is not captured by any approximate XC functional. In the low density (small or large n) Hooke's atom, the system forms a "Wigner molecule," with the maximum of the density located away from the center of the harmonic trap. It is interesting to analyze how the height v c,kin (0) of the peak scales when the system becomes very correlated ( → 0 ); as shown in Fig. 14 it seems to saturate when one uses the high-density scaling. For any two-electron wavefunction of the form (r 1 , r 2 , r 12 ) = e − 2 (r 2 1 +r 2 2 ) p(r 12 ) , the peak's height is given by the simple expression We have used up to the second order of the small-(strong correlation) expansion of the exact wavefunction [25], finding that in the scaling used in Fig. 14 the peak does not saturate, but eventually will decrease and then go to zero very slowly, Fig. 14 Two kinetic correlation energy densities of Eqs. (22) and (23) for the different Hooke's atoms considered here. The high-density scaling is applied Fig. 15 Peak v c,kin (0) as a function of . The first three orders in the small-(strong correlation) expansion are compared with the values (dots) from the exact wavefunctions of Taut [24], and with the large-(weak correlation) expansion as 1∕6 . In Fig. 15, we show the peak's height as a function of for the analytic solutions, compared to the first three orders in the small-(strong correlation) expansion (Eq. (32) of [25]) and with the large-(weak correlation) expansion (Eq. (22) of [25]). We see that the strong-correlation expansion for the peak is much more accurate than ordinary perturbation theory from the weak correlation limit even for very moderate correlation (the Hooke's atom with = 1∕2 resembles the He atom as far as the degree of correlation is concerned).

Conclusions
We have analyzed the performances of exchange-correlation functionals built from global and local interpolations between the weak-and the strong-interaction limits of DFT for the Hooke's atom series. This case study allows for the use of exact analytical input ingredients, thus disentangling the errors coming from the interpolation itself from those on the input quantities. Surprisingly, we have found that for these systems the global interpolations always outperform their local counterparts, in striking contrast with what had been observed so far for small chemical systems [14,15].
We have also compared two different definitions of the kinetic correlation energy density, which plays a crucial role for strongly correlated systems [28,29], and that can help in understanding how to extend to the continuum a KS theory that recovers the exact kinetic energy density recently proposed for lattice models [54].
(37) W LB = W ∞ + (y + y 4 ) , Using Eq. (2), the LB XC functional is found to be Also the LB functional does not use the information from W � ∞ .