Møller–Plesset Adiabatic Connection at Large Coupling Strengths for Open-Shell Systems

We study the adiabatic connection that has as weak-coupling expansion the Møller-Plesset perturbation series, generalizing to the open-shell case previous closed-shell results for the large-coupling limit. We first focus on the hydrogen atom with fractional spins, providing results along the adiabatic connection from small to large coupling strengths. We reveal an intriguing phase diagram and an equation for the large-coupling leading order that has closed-form solutions for specific choices of its relevant quantum numbers. We then show that the hydrogen atom results provide variational estimates for the large-coupling leading terms for the general many-electron open-shell case in terms of functionals of the Hartree–Fock α-spin and β-spin densities.


INTRODUCTION
−20 The so-called flat plane conditions 12,13,[15][16][17][18]20 that can guide the construction of approximate functionals 21−23 were derived by using the prototypical example of the spin dependence in the H atom.
−29 All three of these terms are functionals of the HF density, ρ HF , only.Using these results, functionals that interpolate between the weak and strong coupling limits of the Møller−Plesset adiabatic connection (MPAC) have been introduced using the exact HF exchange and MP2 correlation energies combined with the strong-coupling limit. 30These interpolated functionals have been shown to massively improve MP2 interaction energies for a wide variety of noncovalent interactions, ranging from small charge-transfer complexes to larger π−π-bonded complexes. 30Generalizations to include other variants of MP2, such as opposite-spin only and regularization, have been found to be even more accurate at a lower computational cost. 31−45 A fundamental difference between the two ACs is that in the DFT one, the density remains fixed as the coupling constant λ is turned on, whereas in the MPAC, the density can roam freely.In the DFT AC, the role of the spin state has been found to enter in the large coupling limit only at orders 36,41,46 e , which means that it can be ignored in the two leading λ → ∞ terms used in the interpolating functionals.
However, in the MPAC case, the spin state already affects the second leading term at strong coupling due to the role of the exchange operator at this order and the lack of the density constraint. 20,26This spin-dependence becomes easy to study in the closed-shell case, where it has been shown 26 ) provides a variational estimate for the general many-electron case.However, this unnecessarily restricts the chemical space for which the new MPAC functionals can be used.Studying how the spin affects the MPAC strong coupling limit and also its role along the whole AC path beyond the closed-shell case is the gap that we fill in this work.
The paper starts with an introduction of the MPAC in Section 2, including a summary of previous results for its strong coupling limit.Since the closed-shell many-electron case was obtained by generalizing the result for H , , we start by studying the H atom MPAC beyond the spin-unpolarized case in Section 3, where we derive and solve numerically the relevant equations, revealing an interesting phase diagram along the AC path.The spin-dependence of the strongcoupling limit coefficients is then extracted in Section 4. As we shall see, this limit defines an equation that has closed form solutions only for some special values of its parameters (orbital angular momentum and spin).We then show in Section 5 that the results for the H atom with fractional spins provide a variational estimate for the strong-coupling MPAC functionals for the general many-electron case in terms of the HF α-spin and β-spin densities.Conclusions and perspectives are discussed in Section 6.

MøLLER−PLESSET AC
We keep the notation general, such that the following equations apply both to restricted (the spatial part of the α and β orbitals is forced to stay the same) and unrestricted openshell HF.
Fixed in terms of these spin−orbitals (which are determined in the initial HF calculation), the standard Hartree and exchange operators = [ ] J J HF and

HF
(that appear in the initial HF equations) do not change along the AC defined below.Subsequently, treating J ̂and K ̂as fixed (λindependent) one-body operators, the MPAC for this Nelectron system is represented by the generalized (λ-dependent) Hamiltonian where T ̂, V ext , and The Hellmann−Feynman theorem implies

While
ext is the physical Hamiltonian (including the two-body interaction V ee ), the one-body operator is the original HF Hamiltonian whose ground state |Ψ 0 ⟩ = |Ψ λ=0 ⟩ is the Slater determinant made of the occupied HF spin−orbitals {ϕ i HF }.The HF energy of the N-electron system is defined as

HF
where the Hartree energy while is the usual HF exchange energy.The difference between the physical ground-state energy HF and E HF is the HF correlation energy Here we have introduced the MPAC integrand where we have eq 3 for E d d

HF
. The Taylor expansion of W c,λ for λ → 0 is the MP perturbation series This expansion holds for closed systems.The first part of this work, however, addresses open fragments of larger systems (for example, an H atom within an infinitely stretched H 2 molecule).In such cases, eq 6 applies to the whole system, while for the subsystem, we may find W c,λ=0 ≠ 0. An explicit example is reported in Figure 10 of ref 26, where it is shown that the MPAC result for the stretched H 2 molecule tends to twice the result for the H atom with 1 2 spin-up and 1 2 spin-down as the distance R between the two H atoms increases, except at λ = 0, where the order of limits, R → ∞ and λ → 0, matters.
2.1.Previous Results on the λ → ∞ Expansion.The counterpart of eq 6 is the large coupling strength (λ → ∞) expansion 25,26 The Journal of Physical Chemistry A For general N-electron systems (atoms or molecules with M fixed nuclear positions R k with charge numbers Z k , where k = 1, ... , M), the leading term is 25 Here, E el [ρ] is the classical electrostatic energy of N negative point charges (classical electrons) sitting at equilibrium positions in a rigid continuous positive background charge distribution with given density ρ(r) with the electrostatic (Hartree) potential due to the charge distribution ρ(r) While the minimizing set {r 1 min ,...,r N min } of equilibrium positions in eq 9 is typically not unique (depending on the symmetry group of the molecule), the set of N density values ρ HF (r i min ), for i = 1, ..., N, should be unique.Independently, we expect a certain number I nuc ≤ N of these positions r i min to coincide with some of the fixed nuclear positions R k (with k = 1, ..., M), implying that I nuc ≤ M.Then, after relabeling the r i min if necessary, we have In terms of these values ρ HF (r i min ), the coefficients of the remaining terms in eq 7 have been shown 26 to have the variational estimate for closed-shell N-electron systems Here, = 1.6185 ( ) ) For the hydrogen atom , and the only minimizing position r 1 min = 0 coincides with the only nuclear position R 1 = 0, where minus the Hartree potential has its minimum.In ref 26, only the two cases s = 1 of a spin-polarized regular atom H[1,0] and = s 1 2 of a spin-unpolarized ensemble H , were studied.The values of the three quantities ρ HF (r), 1/2 , and 1/4 are s dependent, with the latter two for reported again in Table 1 for completeness.
As a particular feature, the expansion ( 7) has a term O(λ −3/4 ).Such a term is absent in the corresponding λ → ∞ expansion for the density-fixed AC in DFT.According to eq 12, the term O(λ −3/4 ) in eq 7 occurs only in molecules with I nuc > 0. For a short explanation, 26 we note that the spatial probability distribution of the N electrons in the state HF as λ → ∞ and therefore O(λ −3/4 ) to W c,λ .Moreover, due to the Kato cusps of ρ HF (r) at nuclear positions, the term case was proven to yield variational estimates for the many-electron closed-shell case. 26

MøLLER−PLESSET ADIABATIC CONNECTION FOR THE H ATOM WITH FRACTIONAL SPIN
In this section, we generalize and compute the MPAC for the hydrogen atom at arbitrary values of s 1 1 2 of the weight parameter s.This parameter s, defined in eq 17 below, must not be confused with the spin quantum number.In Section 5, we show that in the large-λ limit of the MPAC for a general open-shell system, the parameter s is linked to the local spin polarization, eq 54.
Here, we obtain results of the MPAC for λ ∈ [0, ∞) for the H atom with general s, which, as we shall see, reveal an interesting phase diagram.

HF Orbital ϕ s (r). We consider a one-electron system
Instead of being in a pure quantum state, however, this system is described by an ensemble of a spin-up state ϕ α (r)|α⟩ and a spin-down state ⟨ϕ(i)|α⟩ with weights of 1 − w and w, respectively.Statistically, our system has a spin In this study, we stay in a restricted open-shell HF (ROHF) framework, forcing both spin states to have the same real spatial orbital ϕ(r) = ϕ s (r), which will depend on the weight w via 20 the parameter s The influence of the restricted open-shell choice made here for the generalization to the many electron open shell case ofthe unrestricted case is discussed in Section 5.
Then, ϕ s (r) 2 = ρ s HF (r) will be the HF density of eq 1, with the Hartree energy U[ρ s HF ].The exchange functional for this ensemble system is explicitly weight-dependent 20

[ ] =
[ ] E sU s x, 2 In the pure-state cases (w = 0 or w = 1), E x,s [ϕ] exactly compensates for the spurious Hartree interaction U[ϕ 2 ].In the ensemble case (0 < w < 1), this compensation is incomplete. 15,20This one-electron RHF functional of ϕ explicitly reads In cases with s ≠ 1, the spherically symmetric (real-valued) minimizer (19)   will be different from the hydrogen ground state ψ 1s (r).However, the radial wave function R s (r) is still finite at r = 0 and satisfies Kato's cusp condition R s ′(0) = −ZR s (0). 48he nonlinear self-consistent field Euler−Lagrange equation for ϕ s (r) corresponding to the minimization in eq 19 is solved using a basis set expansion in terms of Slater-type orbitals, as in ref 26

Møller−Plesset Adiabatic Connection: Equations and Numerical Solutions.
Employing the HF orbital ϕ s (r), fixed by eq 19 for a given value of s , 1 , we now consider the λ-dependent Hamiltonian of eq 2 for one-electron systems 3.2.1.Operators J ̂and K s .For the present case of oneelectron systems with fractional spin, the Hartree operator J ̂[ϕ] and the (explicitly weight-dependent) exchange operator [ ] K s are defined by their action on a general (single-particle) spin orbital while [ ] K s is nonlocal and, in addition, acts on the spin variable σ in Ψ(x) = Ψ(r, σ).Equivalently, writing

Wave Function for Arbitrary
HF ext is the physical Hamiltonian of the H atom, which has a spin-degenerate ground state.Any superposition of the kind where ψ 1s (r) is the hydrogenic 1s orbital with q ∈ [0, 1] is a valid ground state.Alternatively, instead of the superposition, we can consider again a statistical ensemble.Generalizing eq 25 to cases with arbitrary λ > 0, we write the ground state Ψ s,λ (x) of H s, where we have chosen q = w, forcing the spin expectation S z in this pure state to equal the ensemble average = m w s 1 2 . In other words, we suppress spin flip, and we stay on the restricted open-shell curve by enforcing the spatial orbital to be the same for both spins.We also notice that we could alternatively use an ensemble in eq 26 instead of a superposition.This does not change the result for the energy along the AC since the exchange kernel of eq 24, being diagonal in the spin part, yields the same expectation value for a superposition or an ensemble.All the equations reported below, in which the spin-dependence is explicitly transformed into a weight dependence in the MPAC Hamiltonian, are thus the same whether for the wave function at λ > 0, we use a superposition or an ensemble.The only constraint that matters is forbidding spin flips, which we enforce to keep the AC curve smooth.In fact, if we allow the spin to relax, the MPAC has a discontinuity 20,26 as we cross λ = 1 (except in the case = s 1 2 ).

The Journal of Physical Chemistry A
Since our aim is to build interpolations by using the information at large λ, we want to follow the AC that connects smoothly the λ ∈ [0, 1] region with the λ → ∞ limit.With eq 26, the expectation of the Hamiltonian ( 21) can be written as Here, we have performed the spin summation, turning [ ] K s into the simpler operator K ̂[ϕ], which no longer acts on spin but instead needs the weight parameter s from eq 17 as a prefactor in eq 27.As mentioned before, this holds regardless of the choice of using for Ψ(x) a superposition or an ensemble.

Evaluation of the MPAC Integrand W c,λ .
To evaluate the weight-dependent MPAC integrand of eq 5 for the present system we need the ground state energy E s,λ HF of H s, Since our (fixed) HF orbital and u(0) = 0, in Then, the minimum ( 30) is obtained as For any fixed value of s, the λ-dependent minimizer = s (λ) can jump between different integers , as λ continuously grows from λ = 0 to λ = ∞.Plotted versus λ, E s,λ HF will be continuous with possible kinks (and corresponding jumps in the derivative The Euler−Lagrange equation for the minimization (31) is solved for a given {s, } pair with the spectral renormalization method, 49−52 following the algorithm described in ref 26.HF starts at = 0 for small λ, then a first crossing of states from = 0 to = 1 occurs at λ = 2.3, followed by a second crossing back to = 0 around λ = 11.5, as shown here again on the bottom left panel of Figure 1.Instead, for the case = s 1 2 , the = 0 state was found to be the lowest for all λ ≥ 0. This means that the integrand W c,λ has discontinuities in the case s = 1 (see the top left panel Figure 1), while it is continuous for = s 1 2 .The first question we address is thus whether the two crossings of states persist as we lower s from 1 to 1 2 and how they eventually disappear at = s 1 2 .We find that, as we lower s starting from s = 1, the region of λ values for which = 1 is the ground state shrinks, until it disappears entirely at s = 0.810.
In Figure 2, we report the resulting phase diagram in the λ, s plane, showing the regions in which each gives the lowest energy.We observe two distinct regions: the first, with 0.5 ≤ s ≤ 0.810 (0.1063 ≤ w ≤ 0.8937), has no crossing of states, with = 0 being the ground state at all λ ≥ 0, which means that there are also no discontinuities in W c,λ .The second region, s > 0.810, has two crossings of states, with an intermediate λ range in which the = 1 channel gives the lowest energy.The > 1 channels have been found to always have higher energy in the physical range s , 1 Similarly to constraining the spin, one can also follow the AC along the constrained = 0 channel.This removes the crossings of states that introduce discontinuities in W c,λ ; see, as an example, the right panels of Figure 1.Since the crossings always happen for λ > 1, interpolating between the large and small λ limits along the = 0 curve allows us to make approximations of W c,λ without affecting the resulting correlation energy.The only exception remains s = 1, which still has a discontinuity, this time at λ = 2.5 (see the left panels of Figure 1), due to a crossing between the flat 1s curve E s=1,λ HF = −0.5 and the second = 0 state, with a radial node.However, in all the other cases, the transitions are smooth, meaning that continuous curves for 0.5 ≤ s < 1 can be obtained when constraining = 0.

λ → ∞ COEFFICIENTS FOR THE H ATOM
In this section, we compute the spin-dependence of the first three leading terms of the large-λ expansion of the H atom MPAC studied in the previous Section 3. We will then show in the next Section 5 that, similarly to the closed-shell case, 26 these coefficients can be used in the general many-electron case.
We start from the Euler−Lagrange equation in eq 33.As λ → ∞, the effect of the term −λv s H (r)u(r) becomes dominant: 26 despite the presence of the quantum-mechanical kinetic-energy term and of the exchange operator, the solution u(r) will for λ → ∞ concentrate indefinitely at the minimum of −v s H (r) at r = 0, implying (34)   To this explicitly, eq 37 below, we expand 26 in eq 33 where we have used the cusp condition R s ′(0) = −ZR s (0) (see Section 3.1 above).We emphasize that v s H (r), in the hypothetic case R s ′(0) = 0, would have no term O(r 3 ).In terms of the scaled coordinate , eq 33 takes the form 1/4 0 (37)   we see that ϵ 1/2 is the zero-order eigenvalue ϵ (0) in the perturbation expansion for the operator while ϵ 1/4 is the corresponding first-order correction ϵ (1) Via the operators H 1/2 and H 1/4 , the values ϵ 1/2 and ϵ 1/4 depend on = R (0) 4 (0)

HF
. This dependence is easily revealed when we write = = ( ) , with . Then, eq 38 becomes a universal equation 26 with two parameters, and s Numerically, the eigenvalue = s ( , ) is always lowest for = 0, independently of s, confirming the results of the phase diagram of Figure 2. Therefore, for the determination of the large-λ coefficients, we can set = 0 everywhere.Then eq 39 reads We then see that 1/2 and 1/4 for = 0 are pure functions of s, see Figures 3 and 4.
In summary, eq 37 yields for for large λ → ∞ (when the minimizer in eq 32 is always = 0), the expansion of eq 7 where the coefficients now depend on the weight parameter , we set = 0 in eq 40 and expand u 1/2 (p) on the basis of the quantum isotropic harmonic oscillator (IHO) problem that arises if we set s = 0 in eq 40, which has frequency = 1 3 , and energies given by , with L n 1/2 being generalized Laguerre polynomials.Numerical solutions of eq 40 with = 0 for different values of s , 1 have thus been obtained both by using the IHO basis set expansion of eq 42, which converges very fast, and by using the spectral renormalization method, finding perfect agreement for s ( ) 1/2 and s ( ) , which are shown, respectively, in Figures 3 and 4.
In both figures, s are accompanied by a quadratic fit, which proves to accurately interpolate the results and can be used for all practical purposes with maximum absolute deviations of 0.0026 and 0.0032, respectively.
In appendix A, we report an intriguing curiosity regarding eq 40: it has closed-form solutions for special pairs of s and values.It was already noticed in ref 26 that the case s = 1, = 0 has a simple closed-form solution, and in appendix A, we investigate the structure of eq 40 further, finding an infinite set of such solutions.Unfortunately, they all appear at s > 1 and thus have no physical meaning in the context we are analyzing.

GENERAL MANY-ELECTRON CASE
We show in this section that the λ → ∞ results for the H atom at s , 1 provide a variational estimate for the strongcoupling leading terms of the general many-electron case in terms of functionals of the HF α-spin and β-spin densities.
The derivation is a generalization of the one for the spinunpolarized (closed-shell) case: 26 we start from a variational ansatz for the wave function Ψ λ→∞ that minimizes the Hamiltonian (2) when λ → ∞.We use a simple Hartree product of localized orbitals around each minimizing position r i min of eq 9, with each spin in a superposition (again, for all the expectation values we need to consider in the derivation, it is

The Journal of Physical Chemistry A
the same if we use an ensemble for the spin part instead of a superposition).The antisymmetry of the wave function can be neglected as it contributes to orders e to the energy. 25,26,36The Hartree product then reads which needs to be determined variationally, 0 ≤ q i ≤ 1, and we will set later = n 1 4 . 25,26In other words, we know that the wave function squared |Ψ λ→∞ | 2 tends to be a product of delta functions centered around the minimizing positions r i min , and with our ansatz, we seek the best variational spherical representation of the delta function that minimizes the next leading term.This approach does not take into account the coupling between the localized states and their anisotropy.As such, it can only provide a variational upper bound for the λ → ∞ functionals.The same kind of approximation was used by Wigner 53 to compute the zero-point energy in the low-density electron gas, yielding an error of ∼12% with respect to the full coupled exact solution, 54 which could provide an indication on the tightness of the upper bound we provide.
We thus evaluate the expectation of the Hamiltonian H HF of eq 2 on Ψ λ h , and we retain only the leading orders at λ → ∞.The kinetic energy and V J ee are spin-independent: their expectation values at large λ are the same as for the closed-shell case considered in ref 26, which we report here again for completeness and and its expectation on Ψ λ h at large λ is then  where we have expanded the HF spin−orbitals ϕ a,σ in scaled coordinates t i = λ n (r i -r i min ) at large λ and we have introduced the usual spin-polarization parameter for the HF α-spin and β-spin densities Now we see that the local spin polarization at the minimizing positions plays exactly the same role as w in our derivation for the H atom, with = w With this choice, the leading (λ −2n ) term of eq 49 becomes We thus see that if we set (which varies between , and introducing the scaled variable p = (4π ρ) 1/4 t, we obtain eq 40 with = 0.This means that the best possible spherical variational ansatz for i is the same as the one we found for the s-dependent H atom in Section 3, around each equilibrium position r i min , yielding the variational estimate where we can use for s ( ) the quadratic fit of Figure 3, and s(r) is a functional of the HF α-spin and β-spin densities via eqs 51 and 54.
The next leading order works exactly as in the closed-shell case. 26Its generalization to the open-shell case is then where the sum runs only over the nuclear positions that coincide with minimizing electronic positions r i min , as in eq 10.For the function s , we can use the quadratic fit of Figure 4.
We notice that these variational estimates are strictly valid for the restricted open-shell HF case only.For the unrestricted case, there would be an additional dependence on the α-spin and β-spin densities entering when we solve the s-dependent equation for the H atom, as the pair ϕ s,α (0) and ϕ s,β (0) appear in those equations, and the resulting functions s ( ) 1/2 and s ( ) 1/4 will be slightly different.However, we may expect these effects to be much smaller than the main s-dependence studied here (see, e.g., Figure 1 of ref 20).

CONCLUSIONS AND PERSPECTIVES
We have extended the results 26 for the large-coupling limit of the AC that has as small-coupling expansion the MP series to the open-shell case.We first studied the paradigmatic case of the H atom, revealing an interesting phase diagram (Figure 2), and then showed that the results for the H atom at large coupling strength can be used for the general many-electron open shell case (Section 5), yielding functionals of the HF αspin and β-spin densities.
Using these results, we plan to extend the construction of MPAC functionals 30,31 to open shell systems, either by developing generalized gradient approximations for the leading term functionals, as done for the closed-shell case, 27 or by using inequalities and relationships with the DFT AC case. 30,31APPENDIX A Closed-Form Solutions for Special Values of s at = 0 We start by noticing that the energy, 1/2 , for s = 1 and = 0 is exactly equal to the energy of the first excited state of the isotropic harmonic oscillator, which arises by setting s = 0 in eq 40 = ( 1/ 3 ).Furthermore, the first excited state of s = 1 is again equal to the energy of the second excited state of s = 0 and has as a closed-form solution a linear combination of the ground state and the first two excited states.In summary, the energies of s = 1 and s = 0 are related, and their structure is shown in the first three columns of Table 2. Notice that in order to yield an energy of n (7 4 ) 2 for s = 1, one needs to mix in all IHO's orbitals until n.
One can find other values of s that give closed-form solutions to eq 40 at = 0.The next one, = s 10 3 , is degenerate with the second excited state of the IHO (s = 0), = 1/2 11 2 3 , and is again a linear combination of the ground state and the first two excited states at s = 0. Closed-form ground-and excited-state energies for s = 7, s = 12, and = s 55 3 can also be obtained, see Table 2.One can find a formula for the values of these special s m and the corresponding energies s ( ) The Journal of Physical Chemistry A

For a system
with N α spin-up and N β spin-down electrons in a given external potential v ext (r), a HF calculation amounts to minimizing the expectation value of the physical Hamiltonian over single Slater determinants only, yielding N = N α + N β occupied HF spin orbitals ϕ i HF (r,σ) = ϕ i HF (x), with the HF electron density and the α-spin and β-spin densities

s
is spherically symmetric, the same is true for the Hartree potential in eq 23, v s H (r) = v s H (r), and the problem becomes block-diagonal in the orbital angular momentum .Therefore, the minimization (30) is performed separately for each ,

Figure 1 .
Figure 1.W c,λ graph of s = 1 (w = {0, 1}) (top left) and s = 0.82 (w = {0.1,0.9}) (top right) when is unconstrained and when we constrain = 0, with the blue shaded area being the correlation energy.The bottom two panels contain the corresponding = 0 and = 1 of the E sλ HF curves, showing the crossings of states between the two channels.

Figure 2 .
Figure 2. Phase diagram of E HF , showing in each region the angular momentum corresponding to the lowest energy.

Figure 3 .
Figure 3. Dependence of 1/2 on s between 0.5 and 1 and a quadratic fit.

Figure 4 .
Figure 4. Dependence of 1/4 on s 0.5 and 1 and a quadratic fit.
whereC = E el [ρ HF ] − U[ρ HF ].The expectation of K ̂is the only part that changes with respect to the closed-shell case.The kernel of K ̂for a general open-shell system with N α spin-up electrons and N β = N − N α spin-down electrons reads we want to use the λ → ∞ expansion to build interpolations, we need to forbid spin flip (to keep W c,λ continuous) and thus set

Table 2 .
n increasing the excitation rank.The corresponding wave functions for these values of s are Structure of the Closed-Form Solutions of eq 40, Where the c n 's of Every Combination of s

Table 1 .
Value of .