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 ). In this case, the series expansion of the correlation energy at small coupling strengths (λ → 0) is given by the Görling–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]

Table 1

Two Adiabatic Connections (ACs) Considered in This Worka

	DFT adiabatic connection	Møller–Plesset adiabatic connection
Ĥλ	T̂ + V̂ext + λV̂ee + V̂λ[ρ]	T̂ + V̂ext + V̂HF + λ (V̂ee – V̂HF)
		V̂HF = Ĵ[ρHF] – K̂[{ϕiHF}]
ρλ=0	ρ	ρHF
ρλ=1	ρ	ρ
ρλ	ρ	ρλ
Wc,λ	⟨Ψλ|V̂ee|Ψλ⟩ – ⟨Ψ0|V̂ee|Ψ0⟩	⟨Ψλ|V̂ee – V̂HF|Ψλ⟩ – ⟨Ψ0|V̂ee – V̂HF|Ψ0⟩
Ec	∫01Wc,λDFT dλ	∫01Wc,λHF dλ
Wc,λ→0	∑n=2∞nEcGLnλn–1	∑n=2∞nEcMPnλn–1
Wc,λ→∞	Wc,∞DFT + W1/2DFTλ–1/2 + O(λ–5/4)	Wc,∞HF + W1/2HFλ–1/2 + W3/4HFλ–3/4 + O(λ–5/4)

Middle column: the standard density-fixed DFT AC that starts at λ = 0 with the Kohn–Sham determinant. Right column: the AC that has the Møller–Plesset series as expansion for small coupling strengths λ and starts at λ = 0 with the Hartree–Fock Slater determinant.

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 ), 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 strong-coupling 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 refs[23, 24, 29] should be present in the second-order gradient expansion for exchange.

The paper is organized as follows. In section , 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 , 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 , along similar lines, we extract the GEA2 coefficient for the next leading term of the MP AC large-coupling-strength expansion. The computational details are described in section , and the last section (section) 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

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̂ee the electron–electron repulsion, and V̂ext the external potential due to the nuclei. The operators Ĵ = Ĵ[ρHF] and K̂ = K̂[{ϕHF}] are the standard Hartree–Fock (HF) Coulomb and exchange operators in terms of the HF density ρHF and the corresponding occupied orbitals ϕ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 bywith WHF being the MP AC integrand,and Ψλ the wave function that minimizes the expectation value of the Hamiltonian of eq . The last two terms, U[ρHF] and E[{ϕHF}], are the Hartree energy and the HF exchange energy, respectively, whose sum gives minus the expectation value of V̂ – Ĵ + K̂ on the HF Slater determinant (see right column of Table ). The small-λ expansion of WHF is the familiar MP perturbation series,

The λ → ∞ Expansion of the Møller–Plesset Adiabatic Connection

The large-λ expansion of the MP AC has recently been uncovered[19,20] for closed-shell systems as follows:The leading order, eq , contains the electrostatic-energy density functional Eel[ρ], which entails a classical electrostatic minimization,withandThe density functional Eel[ρ] 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 {r1min... rmin} in eq , in turn, determine the next leading term, for which eq provides a rigorous variational estimate for closed-shell systems.[20] This term is given by zero-point oscillations around the minimizing positions enhanced by the exchange operator K̂, which mixes in excited harmonic oscillator states.[20] Finally, the sum in eq only runs over those minimizing positions of eq 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 , 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, W3/4HF of eq , 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.

Comparison with the λ → ∞ Expansion of the DFT AC

In a recent work where an interpolation for WHF between MP2 and the λ → ∞ limit has been built and tested,[18] the functional WHF of eq has been approximated in terms of the strong interaction limit of the DFT AC, using the following inequality:[19]The DFT AC of the central column of Table uses the Hamiltonian,with V̂λ[ρ] = ∑vλ(r, [ρ]) being the one-body potential that forces the density to be equal to the physical one for all values of λ. With this Hamiltonian, the KS exchange-correlation (XC) energy is given bywhere the DFT coupling constant integrand isand ΨλDFT[ρ] is the wave function that minimizes the expectation value of eq . Although the λ → ∞ expansion of the DFT AC has a similar form as the MP AC one of eq , 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 {r1min... rmin}. The density approaches asymptotically, as λ → ∞, a sum of Dirac delta functions centered around these minimizing positions. If one of the rmin 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 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 Ĵ[ρ] 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 byThe functional Eel[ρ] of eq , in contrast, is obtained by letting Ψ relax to its minimum in eq , which directly impliesAdding EHF to both sides of this inequality yields eq .

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 to be fitted to the S22 dataset.[30,31] This inequality was used due to the lack of approximations for Eel[ρ], but also to allow the functional to be more flexible to approximate the missing but very large second-order 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:with APC = and BPC = ≈ 0.005317. The PC model was built from the physical interpretation of W∞DFT[ρ] provided by eq : perfectly correlated electrons that need to minimize their interaction while giving the same density ρ of the classical positive background will tend to neutralize the classical charge distribution ρ (which is different than minimize the total electrostatic energy, as in Eel[ρ]). Along similar lines, by considering zero-point oscillations around the PC positions, a GEA2 functional for the second-order term, was constructed:[21]with CPC = 1/2(3π)1/2 ≈ 1.535 and, from ref (7), DPC = −0.028957, where DPC 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 Eel[ρ]

In this section, we wish to derive a gradient expansion for Eel[ρ] of eq . As detailed in Appendix C, we cannot proceed along lines similar to the derivation of the PC model used for the DFT AC, because the charge distribution with the electrostatic energy Eel[ρ] cannot easily be divided into weakly interacting cells. Moreover, we are only interested in Eel[ρ] for ρ(r) that are HF densities of atoms and molecules. For this reason, we follow the procedure used for the DFT exchange functional EDFT[ρ] in refs (23 and 24) with some modifications. This procedure extracts the GEA2 coefficient from accurate data, and is very suitable because, under uniform coordinate scaling at fixed particle number N,our functional Eel[ρ] displays the same scaling behavior as exchange,since vH([ργ], r) = γvH([ρ], γr) and U[ργ] = γU[ρ].

In practice, we wish to determine whether, for slowly varying densities, Eel[ρ] 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 . 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 can also be written asThe GEA2 expression should become more and more accurate as x → 0, and our goal is to determine the values of AHF and BHF. As we will discuss later, while AHF is universal, the coefficient BHF seems to be dependent on how the slowly varying limit is approached, similarly to what happens with the DFT exchange functional.[23,27]

LDA Coefficient AHF

The uniform density limit N → ∞ of a constant droplet density with radius R = N1/3r, taken per particle:is equivalent to the jellium case[33] and has been analyzed already in ref (20). The result is the Wigner crystal energy per particle,[34], leading toNotice that AHF = ADFT, where the latter is slightly different than the PC value APC, which replaces 0.8959 ... with 0.9. Note the rigorous proof in refs (33 and 35) that ADFT is in fact exactly given by the Wigner crystal.

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 , 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 becomes increasingly weak,provided thatExamples of relevant values of p areFor any density functional G[ρ] that, under the uniform coordinate scaling of eq , behaves asfor a fixed profile ρ̅, all the different choices of p in eq are equivalent and simply related to the case p = 0,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).

p = 1/3: the Thomas–Fermi scaling of neutral atoms[23,28,37]N = Z;

p = −2/3: the Thomas–Fermi scaling of the Bohr atoms;[27,28,38,39]

p = 0: the scaling used in refs (40−42) to analyze the Lieb–Oxford bound;

p = −1/3: the scaling used in ref (43) to analyze the asymptotic exactness of the local density approximation.

When N grows, because of eq , we expect the gradient expansion of eq to become more and more accurate. Then, by inserting the density ρ̅ into eq , one gets the following large-N expansionwhereClearly, eq holds only as long as the integrals ILDA[ρ̅] and IGEA2[ρ̅] are finite for the given density profile ρ̅. In refs (23, 24) the fact that the neutral atoms densities for large N asymptotically satisfy the Thomas–Fermi (TF) scaling with p = 1/3,with ρ̅TFna(r) the TF profile (integrating to 1) for neutral atoms,[25−28,37] led to the conclusion that their exchange energy, as a function of N = Z should have, to leading orders, the large-N expansion aN5/3+bN. Extracting b from exchange energies of neutral atoms allowed to fix the GEA2 coefficient for exchange in ref (24). However, while the GEA2 integral IGEA2[ρ] for neutral atoms is finite, the integral IGEA2[ρ̅TFna] for the asymptotic TF profile diverges (while ILDA[ρ̅TFna] is also finite). This does not automatically imply that IGEA2[ρ] should not increase linearly with N, as expected from eq with p = 1/3, since TF theory should not give exact information at this order. Nonetheless, we find numerical evidence (see Figure ) that the GEA2 integral IGEA2[ρHF] for Hartree–Fock densities of neutral atoms increases as N log(N) rather than as N.

Figure 1

GEA2 integral of eq for the Hartree–Fock densities of neutral atoms, IGEA2[ρ], divided by the number of electrons N (log scale on the x-axis). Numerical values (red dots) are compared with a logarithmic fit (blue line).

A case for which it is even simpler to make a detailed numerical analysis of IGEA2 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 with γ = Z,As N → ∞, the densities ρBohr(r) of eq approach the Bohr atom TF profile[27,28,38] ρ̅TFBohr with p = −2/3,Again, ILDA[ρ̅TFBohr] is finite while IGEA2[ρ̅TFBohr] diverges.

If eq would hold with p = −2/3, IGEA2[ρBohr] should have a tendency toward a constant when N → ∞. Instead, we clearly see (Figure ) 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.

Figure 2

GEA2 integral of eq for the Bohr atom densities, IGEA2[ρBohr] (log scale on the x-axis). Numerical values (red dots) are compared with a logarithmic fit (blue line).

A detailed derivation of the behavior of IGEA2[ρ], 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]

Extracting the GEA2 Coefficient BHF

The analysis in the previous section suggests that extraction of the GEA2 coefficient should not be done by using values of Eel[ρ], as a function of N and fitting coefficients from eq , as this seems to be safe only for a scaled known profile (as in eq ), 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 computeThe idea is that if a GEA2 expansion for Eel[ρ] exists, we should observe that B̃(N → ∞) tends to a constant, which will be the sought BHF. 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 , we show for different particle numbers N:The computational details behind the evaluation of Eel[ρ] for each case are described in section . 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 3s2 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 we can conclude that there exists no unique GEA2 and that we should choose one of these BHF’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.

Figure 3

Plot of eq for the four cases described in section .

Our numerical values B̃(N) for the exponential profile

Our numerical values B̃(N) for the Gaussian profile

Our numerical values B̃(N) for the Hartree–Fock densities ρHF(r) of neutral atoms.

Our numerical values B̃(N) for the Bohr atom densities ρBohr(r) of eq , 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.

We noticed that if we fix BHF to make the GEA2 exact for the spin-unpolarized H atom[20] (with 1/2 spin-up and 1/2 spin-down electrons[45]),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 . We thus fix the GEA2 coefficient BHF 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 , we show the relative error of the GEA2 expansion, which, as expected, goes to zero for large neutral atoms and slightly charged ions.

Figure 4

Value BH[HF = −0.0150578 that makes the GEA2 exact for the spin-unpolarized H atom accurately recovers the large N limit of B̃(N) of eq for closed-shell neutral atoms (q = 0) and slightly charged closed shell ions, with q = +1, +2, and −1.

Figure 5

Relative error of the GEA2 expansion for the functional Eel[ρ] with BHF = BHHF[ = −0.0150578 for closed-shell neutral atoms (q = 0) and ions with q = +1, +2, and −1.

However, we should stress that the GEA2 with BHF of eq 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 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.

Figure 6

B̃(N) of eq 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.

Second-Order Gradient Expansion for W1/2HF[ρ]

Once the minimization to obtain Eel[ρ] is performed, we automatically get the functional W1/2HF[ρ] of eq by evaluating the HF density in the minimizing positions rmin. We should still stress that, while the leading term of eq is exact, eq 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 W1/2HF[ρ] satisfies eq with m = 3/2, must have the same form as the one for the DFT case of eq :

LDA Coefficient CHF

Within the variational expression of eq , the LDA coefficient CHF 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.

Extraction of the GEA2 Coefficient DHF

We focus only on the relevant case of closed-shell neutral atoms and slightly charged ions, and, in analogy with eq , we compute and analyze the functionThe results are shown in Figure , 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 Eel[ρ], which therefore remains rather insensitive if the true global minimum is not reached. However, the functional W1/2HF[ρ] 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 , we can get a rough estimate DHF ≈ 0.12 ± 0.01.

Figure 7

Values of D̃(N) of eq for closed-shell neutral atoms (q = 0) and slightly charged closed shell ions, with q = +1, +2, and −1.

Computational Details

To obtain reference values for Eel[ρ] for closed-shell neutral atoms and slightly charged ions, we first performed RHF calculations with PySCF 1.7.6,[46] with the basis sets specified in Appendix A. For a given set of positions {r1, r2, ..., r}, we calculated the value of Vel, whereWe computed the value of vH[ρHF](r) by contracting,where γHF is the Hartree–Fock 1-body Reduced Density Matrix (1-RDM) and the matrix element is given bywhere 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,The total gradient then isFinally, Eel[ρHF] was obtained by minimizing Vel[ρHF] using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm[47−50] as in the scipy.optimize.minimize function of scipy.[51]

For selected cases, such as Ne and Ar, we have also double-checked the minimum by using Mathematica 12.3.1, experimenting with different minimizers. For the scaled densities and the Bohr atoms, we have used both Mathematica and Python with the same scipy.optimize.minimize function used for the HF densities.

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 ) of the adiabatic connection that has the Møller–Plesset perturbation series as small-coupling-strength expansion (see eqs 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 and Figures 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 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 ), although, in this case, obtaining accurate results for large neutral atoms is numerically challenging.

Table 2

Distances from the Nucleus of the Minimizing Positions r1 and r2 (Hartree atomic units), Eel[ρ] and W1/2HF[ρ] of the Symmetric Local Minimum and Asymmetric Global Minimum of H–

	H–[Sym]	H–[Asym]
r1	0.8477	1.2515
r2	0.8477	0.5116
Eel[ρ]	–0.9219	–0.9228
W1/2HF[ρ]	1.4545	1.5003