Gradient Expansions for the Large-Coupling Strength Limit of the Møller–Plesset Adiabatic Connection

The adiabatic connection that has, as weak-interaction expansion, the Møller–Plesset perturbation series has been recently shown to have a large coupling-strength expansion, in terms of functionals of the Hartree–Fock density with a clear physical meaning. In this work, we accurately evaluate these density functionals and we extract second-order gradient coefficients from the data for neutral atoms, following ideas similar to the ones used in the literature for exchange, with some modifications. These new gradient expansions will be the key ingredient for performing interpolations that have already been shown to reduce dramatically MP2 errors for large noncovalent complexes. As a byproduct, our investigation of neutral atoms with large number of electrons N indicates that the second-order gradient expansion for exchange grows as N log(N) rather than as N, as often reported in the literature.


INTRODUCTION
Adiabatic connections (ACs) between an easy-to-solve Hamiltonian and the physical, many-electron one, have always played a crucial role in building approximations in electronic structure theory. In density functional theory (DFT), the standard AC connects the Kohn−Sham (KS) Hamiltonian with the physical one by turning on, via a parameter λ, the electron−electron interaction while keeping the one-electron density ρ(r) fixed 1 (central column of Table 1). In this case, the series expansion of the correlation energy at small coupling strengths (λ → 0) is given by the Gorling−Levy (GL) perturbation theory. 2 In the opposite limit of large-coupling strengths (λ → ∞), the correlation energy is determined by the strictly correlated-electrons (SCE) physics, 3−6 which yields the leading term. The next order is given by zero-point (ZP) oscillations 7−10 around the SCE manifold. A possible strategy to build approximations for the correlation energy is to interpolate between these two opposite limits, generalizing to any nonuniform density 7,11−16 the idea that Wigner 17 used for jellium. The advantage of such an approach is that it is not biased toward the weakly correlated regime. The lack of size consistency of these interpolations can be easily corrected at zero computational cost. 14 More recently, 18 the same interpolation idea has been applied to the AC that has the Møller−Plesset (MP) series as perturbation expansion at small coupling strengths λ (right panel of Table 1), connecting the Hartree−Fock (HF) Hamiltonian with the physical one. The λ → ∞ expansion of this MP AC is given by functionals of the Hartree−Fock density ρ HF (r), with a clear physical meaning. 19,20 The strongcoupling functionals of the DFT and the MP ACs are essentially electrostatic energies, whose exact evaluation for large particle numbers is demanding, but while for the DFT AC there are rather accurate second-order gradient expansion approximations (GEA2) 4,7,21 and, more recently, also generalized gradient approximations 16 (GGA), for the MP AC these approximations are not yet available. For this reason, in a very recent work 18 the λ → ∞ functionals of the MP AC have been modeled in an empirical way, starting from the GEA2 of the DFT ones. Quite remarkably, interpolations combined with this simple empirical model already provide very accurate results for noncovalent interactions (NCI), reducing the MP2 error by up to a factor 10 in the L7 dataset, 22 without spoiling MP2 energies for the cases in which they are accurate. 18 These interpolations work very well for diverse NCI's such as charge transfer and dipolar interactions, and they are able to correct MP2 both when it overbinds and when it underbinds, as they are able to take into account the change from concave to convex behavior of the MP AC. 18 Their appealing feature is that they use 100% of HF exchange and MP2 correlation energy, and it is the interpolation that decides for each system how much to correct with respect to MP2. This way, dispersion corrections are not needed at all to get accurate NCIs. 18 The purpose of this work is to derive the missing GEA2 for the strong-coupling functionals of the MP AC, in order to reduce empiricism and hopefully increase the accuracy of the interpolations along the MP AC. To this purpose, we use the ideas derived from the semiclassical limit of neutral atoms, which have been used in recent years in DFT for the analysis of the exchange and correlation functionals, 23−28 yielding new approximations such as PBEsol. 29 As we shall see, the functionals we need to approximate allow us to probe these ideas more extensively, revealing several interesting features that could be used more generally to build DFT approximations. We also notice that an additional term with respect to refs23, 24, 29 should be present in the second-order gradient expansion for exchange.
The paper is organized as follows. In section 2, we quickly review the large-coupling-strength functionals of the MP AC, discussing their physical meaning and the crucial differences with those of the DFT AC. Then, in section 3, we focus on the gradient expansion of the leading term at large coupling strengths: we perform an extensive analysis by filling more and more particles in a given density profile, and also by considering closed-shell neutral atoms and ions, up to the Bohr atom densities, which provide the limit of highly ionized atoms. We compute the functional in a numerically accurate way and determine a second-order gradient coefficient for the neutral-atoms case. We also discuss differences with the work of refs 23 and 24, providing an analysis that should also be relevant for the exchange and correlation functionals of DFT. In section 4, along similar lines, we extract the GEA2 coefficient for the next leading term of the MP AC largecoupling-strength expansion. The computational details are described in section 5, and the last section (section6) is devoted to conclusions and perspectives. More technical details, a curious behavior of N = 2 ions, and the discussion of an electrostatic model similar to the one used to derive the GEA2 coefficient of DFT are reported in the Appendix. Hartree atomic units will be used throughout this work.

THE LARGE COUPLING STRENGTH FUNCTIONALS
OF THE MøLLER−PLESSET ADIABATIC CONNECTION 2.1. The Møller−Plesset Adiabatic Connection. To start, we must introduce the Møller−Plesset Adiabatic Connection (MP AC), which has the following Hamiltonian: with T̂the kinetic energy, V̂e e the electron−electron repulsion, and V̂e xt the external potential due to the nuclei. The operators Ĵ= J[ρ HF ] and K̂= K̂[{ϕ i HF }] are the standard Hartree−Fock (HF) Coulomb and exchange operators in terms of the HF density ρ HF and the corresponding occupied orbitals ϕ i HF , respectively, which are determined in the initial HF calculation and are not dependent on λ. Notice that, in our notation, K̂is positive definite. This Hamiltonian links the Hartree−Fock system (λ = 0) to the physical system (λ = 1). The HF (or standard wave function theory) correlation energy, using the Hellmann−Feynman theorem, is given by with W c,λ HF being the MP AC integrand, HF ee HF HF HF (3) and Ψ λ the wave function that minimizes the expectation value of the Hamiltonian of eq 1.
The leading order, eq 6, contains the electrostatic-energy density functional E el [ρ], which entails a classical electrostatic minimization, and The density functional E el [ρ] can be understood as the total electrostatic energy of a distribution of N negative point charges and continuous "positive" charges with density ρ(r). In other words, the λ → ∞ limit of the MP AC is a crystal of classical electrons bound by minus the Hartree potential generated by the HF density. 19,20 The resulting minimizing positions {r 1 min ... r N min } in eq 9, in turn, determine the next leading term, for which eq 7 provides a rigorous variational estimate for closed-shell systems. 20 This term is given by zeropoint oscillations around the minimizing positions enhanced by the exchange operator K̂, which mixes in excited harmonic oscillator states. 20 Finally, the sum in eq 8 only runs over those minimizing positions of eq 9 that happen to be at a nucleus, and it is also a variational estimate. 20 These first three leading terms provide a rigorous framework to link MP perturbation theory with DFT, in terms of functionals of the HF density. In practice, we do not want to perform each time the classical minimization of eq 9, which is known to have many local minima and whose cost increases rapidly with N. We rather wish to find good gradient expansion approximations for the first two terms in the expansion described by eq 5. The third term, W 3/4 HF of eq 8, could instead be approximated by making the assumption that, in a large system, there is one minimizing position at each nucleus, transforming it into a functional of the HF density at the nuclei.
2.3. Comparison with the λ → ∞ Expansion of the DFT AC. In a recent work where an interpolation for W c,λ HF between MP2 and the λ → ∞ limit has been built and tested, 18 the functional W c,∞ HF of eq 6 has been approximated in terms of the strong interaction limit of the DFT AC, using the following inequality: 19 and Ψ λ DFT [ρ] is the wave function that minimizes the expectation value of eq 13. Although the λ → ∞ expansion of the DFT AC has a similar form as the MP AC one of eq 5, there are important differences between the two. The first one is the lack of the λ −3/4 term in the DFT AC, which has the following large-coupling expansion: 4,7 The reason why the MP AC can have a λ −3/4 term is that in this case there is no constraint on the density, and the electrons thus localize around the minimizing positions {r 1 min ... r N min }. The density approaches asymptotically, as λ → ∞, a sum of Dirac delta functions centered around these minimizing positions. If one of the r i min happens to be at a nucleus, the nonanalyticity of the Coulomb nuclear attraction and of the cusp in the HF orbitals and density give rise to this term. 20 In the DFT AC case, the density constraint enforces Ψ ∞ DFT to be a superposition of infinitely many classical configurations, 4 so the one with an electron at a nucleus has infinitesimal weight.
The inequality described by eq 12 can be understood on simple physical terms: the functional W ∞ DFT [ρ] can be reformulated as 21 ) where we have simply used the fact that the expectation of J[ρ] on any wave function with density ρ(r) is 2U[ρ]. Then, we can interpret 21 W ∞

DFT
[ρ] as the electrostatic energy of a system of classical electrons forced to have density ρ immersed in a classical background of charge density ρ of opposite sign. Notice that Ψ ∞ DFT [ρ] does not minimize this electrostatic energy, but it is given by The functional E el [ρ] of eq 9, in contrast, is obtained by letting Ψ relax to its minimum in eq 17, which directly implies Adding E x HF to both sides of this inequality yields eq 12. 2.4. Semilocal Functionals for the λ → ∞ Expansion of the DFT AC. In ref 18. parameters were added to both terms on the right-hand-side of eq 12 to be fitted to the S22 dataset. 30,31 This inequality was used due to the lack of approximations for E el [ρ], but also to allow the functional to be more flexible to approximate the missing but very large secondorder term. Although the exact evaluation of W ∞

DFT
[ρ] is even more expensive than the one of W ∞ HF [ρ], an inexpensive 21 but accurate 4,7 approximation called the Point Charge Plus Continuum (PC) model exists, which is a GEA2 functional:   (21) with C PC = 1/2(3π) 1/2 ≈ 1.535 and, from ref 7, D PC = −0.028957, where D PC is fixed to reproduce the helium-atom exact result. 7 In newer work by Constantin, 16 GGA functionals for both terms were derived to fix, among other things, the diverging asymptotics of the XC potentials. 32 However, these GGA's have larger errors than the original PC model when compared with accurate SCE values for small atoms. Notice that, in contrast to the DFT AC (where self-consistent calculations should in principle be carried out), we here do not need the functional derivatives of these quantities, as the MP AC is designed to directly give the HF correlation energy in a post-self-consistent-field manner.

SECOND-ORDER GRADIENT EXPANSION FOR E el [ρ]
In this section, we wish to derive a gradient expansion for E el [ρ] of eq 9.
our functional E el [ρ] displays the same scaling behavior as exchange, In practice, we wish to determine whether, for slowly varying densities, E el [ρ] is well-approximated by a second-order gradient expansion: The powers ρ(r) 4/3 in the two terms of this expression are a necessary consequence of the exact scaling law of eq 23. Defining the usual reduced gradient x of the density ρ, which essentially gives the relative change of the density on the scale of the average interparticle distance, eq 24 can also be written as The GEA2 expression should become more and more accurate as x → 0, and our goal is to determine the values of A HF and B HF . As we will discuss later, while A HF is universal, the coefficient B HF seems to be dependent on how the slowly varying limit is approached, similarly to what happens with the DFT exchange functional. 23 is equivalent to the jellium case 33 Notice that A HF = A DFT , where the latter is slightly different than the PC value A PC , which replaces 0.8959 ... with 0.9. Note the rigorous proof in refs 33 and 35 that A DFT is in fact exactly given by the Wigner crystal.
3.2. Particle-Number Scalings. As discussed in refs 23 and 24, the slowly varying limit can be approached in different ways. An extended system with uniform density can be perturbed with a slowly varying density distortion, but the resulting GEA2 coefficient might not be the one useful for chemistry. 23 More generally 23,36 for any functional that scales as eq 23, we can reach the slowly varying limit by simply putting more and more electrons in a density profile ρ ̅ (r) with ∫ dr ρ ̅ (r) = 1, by generating a discrete sequence of densities with increasing particle numbers N = 1, 2, 3, ..., using the scaling 36 With growing N, for all these densities the reduced gradient of eq 25 becomes increasingly weak, provided that For functionals that do not display a simple scaling behavior, like correlation in DFT, different values of p lead to different interesting regimes, as discussed in refs 23 and 36. When N grows, because of eq 30, we expect the gradient expansion of eq 24 to become more and more accurate. Then, by inserting the density ρ ̅ N,p into eq 24, one gets the following large-N expansion should not increase linearly with N, as expected from eq 34 with p = 1/3, since TF theory should not give exact information at this order. Nonetheless, we find numerical evidence (see Figure 1) that the GEA2 integral I GEA2 [ρ N = Z HF ] for Hartree−Fock densities of neutral atoms increases as N log(N) rather than as N.
A case for which it is even simpler to make a detailed numerical analysis of I GEA2 is the Bohr atoms, 27,28,38 which have densities constructed by occupying hydrogenic orbitals: and can be considered 27,38 as a limiting case for ions with Z ≫ N. The latter have densities that, as Z → ∞, approach those of the Bohr atom scaled as in eq 22 with γ = Z, As N → ∞, the densities ρ N Bohr (r) of eq 38 approach the Bohr atom TF profile 27,28,38 ρ ̅ TFBohr with p = −2/3, Again, I LDA [ρ ̅ TFBohr ] is finite while I GEA2 [ρ ̅ TFBohr ] diverges. If eq 34 would hold with p = −2/3, I GEA2 [ρ N Bohr ] should have a tendency toward a constant when N → ∞. Instead, we clearly see (Figure 2) that it grows as log(N). For this case, everything is analytic and it is easy to reach very large N, evaluating the GEA2 integral to high accuracy.
A detailed derivation of the behavior of I GEA2 [ρ], as a function of N for neutral atoms and for Bohr atoms, confirming the numerical evidence reported here, is also being performed independently by Argaman et al. 44 3.3. Extracting the GEA2 Coefficient B HF . The analysis in the previous section suggests that extraction of the GEA2 coefficient should not be done by using values of E el [ρ], as a function of N and fitting coefficients from eq 34, as this seems to be safe only for a scaled known profile (as in eq 29), but not for atomic densities. For this reason, we follow a route slightly different than the one used for exchange in ref 24. Namely, we directly compute The idea is that if a GEA2 expansion for E el [ρ] exists, we should observe that B̃(N → ∞) tends to a constant, which will be the sought B HF . However, such constant might not be the same for different profiles ρ ̅ or when we use the neutral atoms or the Bohr atom densities. Indeed, this seems to be the case: in Figure 3, we show for different particle numbers N: (1) Our numerical values B̃(N) for the exponential profile (2) Our numerical values B̃(N) for the Gaussian profile ρ N Bohr (r) of eq 38, including some cases in which we did not completely fill all the values for a given principal quantum number n. Notice that these latter cases cannot always be seen as the limit of highly ionized atoms, as degeneracy must be taken into account more carefully. The computational details behind the evaluation of E el [ρ] for each case are described in section 5. We see that these four sequences of data for B̃(N) seem to approach four different limits as N grows. Regarding the Bohr atoms, the cases for which the value of B̃(N) suddenly decreases to a value much closer to the one of neutral atoms are those in which we added an extra pair of s electrons to a completely filled shell. For example, N = 12 is obtained by adding 3s 2 to the filled n = 2 shell, and similarly for N = 30 and N = 62. The case N = 25 is realized by filling the orbitals as in the Mn atom. From Figure  3 we can conclude that there exists no unique GEA2 and that we should choose one of these B HF 's for our new GEA2 functional. As for the case of the exchange functional, 23,24 the most useful value for chemistry should be the one of neutral atoms.
We noticed that if we fix B HF to make the GEA2 exact for the spin-unpolarized H atom 20 (with 1/2 spin-up and 1/2 spin-down electrons 45 (44) we recover the large N limit of closed-shell neutral atoms and closed-shell ions with charges +1, +2, and −1 quite closely, as shown in Figure 4. We thus fix the GEA2 coefficient B HF to this value, which seems to be as good as a fitted one, although we lack at this point a theoretical justification of why the H[ 1 / 2 , 1 / 2 ] should provide such a good number. In Figure 5, we show the relative error of the GEA2 expansion, which, as expected, goes to zero for large neutral atoms and slightly charged ions.
However, we should stress that the GEA2 with B HF of eq 44 misses the other 27 slowly varying limit of Bohr atoms with large N, which can be regarded as the limit 27 Z ≫ N ≫ 1. To better illustrate the issue, we show in Figure 6 the values B̃(N) only for the closed-shell neutral atoms, the Bohr atoms and for selected noble-gas isoelectronic series: we then see how B̃(N) goes from one limit to the other as the nuclear charge Z is increased at a fixed electron number N. An approximation able to cover this entire range of values could possibly be designed as a metaGGA, which is a route that will be investigated in future work.   Once the minimization to obtain E el [ρ] is performed, we automatically get the functional W 1/2 HF [ρ] of eq 7 by evaluating the HF density in the minimizing positions r i min . We should still stress that, while the leading term of eq 6 is exact, eq 7 is only a variational estimate valid for closed-shell systems within restricted HF. 20 Nevertheless, we can repeat the analysis of the previous section to obtain a GEA2, which, because of the fact that W 1/2 HF [ρ] satisfies eq 32 with m = 3/2, must have the same form as the one for the DFT case of eq 21: 4.1. LDA Coefficient C HF . Within the variational expression of eq 7, the LDA coefficient C HF is readily evaluated 20 and equal to 2.8687. Notice that this is not the exact value for a uniform HF density, which should be evaluated by computing the normal modes around the bcc positions in the Wigner crystal and minimizing the total energy in the presence of the nonlocal operator K̂, which will mix in excited modes. This analysis, using the techniques recently introduced by Alves et al., 34 is the object of a work in progress.
4.2. Extraction of the GEA2 Coefficient D HF . We focus only on the relevant case of closed-shell neutral atoms and slightly charged ions, and, in analogy with eq 41, we compute and analyze the function The results are shown in Figure 7, where we see that D̃(N) gets rather flat already at N ≳ 30 around the value of ∼0.11. However, we also see a step to a slightly higher value, ∼0.13, for the largest N. We do not know whether this step is really there or whether it is due to the numerical minimization being trapped in a local minimum. The issue is that, as N increases, there are many local minima with very close values of E el [ρ], which therefore remains rather insensitive if the true global minimum is not reached. However, the functional W 1/2 HF [ρ] is dependent on the minimizing configuration and it changes more from one local minimum to another. We illustrate this in Appendix B for the case N = 2, which undergoes a transition from a symmetric to an asymmetric minimum as the nuclear charge Z varies from 2 to 1. From the data of Figure 7, we can get a rough estimate D HF ≈ 0.12 ± 0.01.  (49) where a very sharply peaked Gaussian G was used to approximate the point charge, which allows for a more efficient computation of the matrix elements using PySCF. To allow for minimization using a quasi-Newton method, we also obtained the gradient of the Hartree potential, Figure 6. B̃(N) of eq 41 for neutral atoms, Bohr atoms and for selected noble-gas isoelectronic series. We see how B̃(N) goes from one limit to the other as the nuclear charge Z is increased at fixed N.

CONCLUSIONS AND PERSPECTIVES
We have built second-order gradient expansions for the functionals of the large-coupling-strength limit (see last line of the right column of Table 1) of the adiabatic connection that has the Møller−Plesset perturbation series as smallcoupling-strength expansion (see eqs 24 and 45). For this purpose, we have used ideas from the literature based on the semiclassical limit of neutral and highly ionized atoms. 23,24,28 During our study, we have also found numerical evidence (see section 3.2 and Figures 1 and 2) that suggests that the way this semiclassical limit has been used to extract second-order gradient coefficients for exchange should be revised. 23,24,29 In future work, we will design and test new formulas for the adiabatic connection (AC) of the right-hand side of Table 1 that interpolate between MP2 and these new semilocal functionals at large coupling, including the term proportional to λ −3/4 , which can be approximated as a functional of the HF density at the nuclei. Previous work 18 showed that such functionals can be very accurate for noncovalent interactions, correcting the MP2 error for relatively large systems without using dispersion corrections. We will also analyze in the same way the functionals at strong coupling of the DFT AC (last line of the central column of Table 1), although, in this case, obtaining accurate results for large neutral atoms is numerically challenging.

■ APPENDIX A: BASIS SETS
For He−Zn 2+ we used aug-cc-pVQZ basis sets from ref 52. For the heavier atoms ranging from Zn to Xe, we used Jorge AQZP, 53,54 except for Br − , where we used a standard aug-cc-pVQZ basis set. We used Jorge ATZP for Cs to No, except for Ba and Ba 2+ , for which we used the Jorge TZP basis set. For the helium iso-electronic series, we used an aug-cc-pV6Z specifically designed for helium, 55 whereas for the isoelectronic series of neon and argon we used a standard augcc-pV5Z basisset. For the krypton iso-electronic series, we used an aug-cc-pVQZ basis set instead.
Here, we report a curiosity that we have observed about the minimizing positions of eq 9. Naively, one would expect the minimizing positions for a N = 2 atom to be symmetrically distributed on both sides of the well created by −v H (r), which is the case for the helium atom (orange curve) in Figure 8, with minimizing positions in blue. However, when the nuclear charge Z is decreased to 1 (H − ), the minimizing positions (red dots) are now asymmetrically distributed. In Figure 9, we show the difference between the distances from the nucleus of the two minimizing positions as Z varies: we see that the change from a symmetric minimum to an asymmetric minimum happens at Z = 1.16. In the case of H − , a symmetric minimum is still present, but it is only a local one. In Table 2   [ρ] (eq 20) is built by starting from eq 17 for a uniform density, in which the Wigner crystal of strictly correlated electrons is approximated by having each electron surrounded by a sphere of background density exactly integrating to 1 (PC cell). Such an approximation amounts to replace 0.9 with the value of the bcc Madelung constant 0.8959 of eq 28. The GEA2 coefficient B PC is then derived by applying a small gradient Γ = |∇ρ(r)| by requiring that the PC cell plus its electron have exactly zero dipole moment, a new modified PC cell is constructed, whose center and size are slightly changed. This way, a given density can be thought of as being composed by cells (the PC cells plus their electron) that are weakly interacting, and the total energy can be obtained as a sum of the electrostatic energy of each cell (equal to the background-background interaction plus the electron-background interaction).
In the case of E el [ρ], instead, we want to minimize the total electrostatic energy. When the density is uniform, there is no other choice than creating the Wigner crystal bcc arrangement, which could be again approximated with PC spherical cells integrating to 1. Now we apply a small gradient and, in analogy with the DFT PC model, we focus on the energy of only one cell. The electron will move from the center of the sphere in the position that minimizes the electrostatic energy, i.e., the minimum of the Hartree potential of the PC cell with dipole. An estimate (with serious limitations discussed at the end of the derivation) of the resulting B HFPC can then be easily computed as follows.
We consider the charge density which is zero outside a sphere with radius R (centered at the origin r = 0) and inside that sphere has a uniform gradient of magnitude |∇ρ(r)| = Γ in the z-direction, and The condition ρ(r) ≥ 0 implies The Hartree potential (inside the sphere) and the Hartree energy are given by 21 ρ κ ρ