Calculation of dispersion interactions with the geminal-based ring Coupled Cluster Doubles method

Abstract

The performance of the recently developed multi-reference extension of ring coupled cluster doubles is investigated for dispersion energy calculations, applied to the generalized valence bond wave function. The leading-order contribution to the dispersion energy is shown to have the correct asymptotic behaviour. Illustrative calculations on noble gas dimers are presented.

1 Introduction

Calculation of dispersion energy between atoms or molecules is a long-standing problem of quantum chemistry. Physically, dispersion interaction arises from the charge density fluctuations of the subsystems, leading to weakly bound states between, for example, closed-shell partners even if they lack permanent electric moments. The importance of dispersion interaction in both chemical and biological systems cannot be overemphasized.

For two subsystems at a large separation distance R, the asymptotic dispersion energy is manifested as [1]

$$\begin{aligned} E_{{\rm disp.}}\approx -\frac{C_6}{R^6} , \end{aligned}$$

(1)

provided that R is not large enough for retardation effects to become noticeable. The dispersion coefficient $C_6$ can be expressed in terms of the frequency-dependent polarizability (FDP) of each subsystem. Providing an equally accurate description of dispersion both in the equilibrium and in the asymptotic regime has remained a challenge for electronic structure methods to this day.

Mean-field methods, such as Hartree–Fock (HF), cannot even qualitatively account for dispersion; such phenomena are usually introduced via correlation corrections, obtained, for example by perturbation theory (PT). Szabo and Ostlund showed [2] that a supermolecular Møller–Plesset PT2 (MP2) calculation gives rise to dispersion energy exhibiting the correct asymptotic behaviour [i.e. conforming (1)], but with $C_6$ calculated only at the level of uncoupled HF [3]. An alternative of the supermolecule approach is the symmetry-adapted perturbation theory (SAPT) of Jeziorski et al. (see Ref. [4] and references cited therein), in which the intersystem interaction is treated as perturbation. The relationship between supermolecular MP and SAPT has been investigated in detail [5].

The quality of the supermolecular approach can be improved by applying corrections more involved than MP2, e.g. the all-order inclusion of direct ring-type Goldstone diagrams (and possibly certain kinds of exchange ring diagrams) [6]. In the random phase approximation (RPA), the effect of these diagrams can be easily taken into account. Different variants of RPA (distinguished by the different degrees of including exchange diagrams) improve quantitatively on the MP2 results, but not all of them lead to a correct expression for the dispersion coefficient [2, 7].

Another way to include higher-order diagrams in the energy is through coupled cluster (CC) [8, 9]. The CC method with single, double and perturbative triple excitations (CCSD(T)) [10] yields essentially exact results for intermolecular interactions and is often used as benchmark in such cases. The CC model also provides a way to sum up only ring diagrams. The ring CCD (rCCD) approach, originally introduced by Čížek [11], is a computationally less demanding alternative to CCSD(T), scaling as ${\mathcal{O}}(n^6)$ (in some cases ${\mathcal{O}}(n^5)$ [12]) instead of ${\mathcal{O}}(n^7)$ (n being the number of basis functions). The formal equivalence of certain kinds of RPA and rCCD methods has been noticed by many [6, 11,12,13].

RPA and rCCD have recently seen a renewed interest in the context of density functional theory [14,15,16], since correlation energy functionals derived from RPA can serve to construct post-Kohn–Sham models [17] and can account for dispersion effects. Development of new methods (e.g. RPA+SOSEX [18, 19] or renormalized PT2 [20]) and refinement of existing ones (e.g. inclusion of explicitly correlated terms in rCCD [21]) have been reported. The combination of range-separated hybrid functionals with RPA/rCCD has also proved to be a successful method for calculating dispersion energy [22, 23].

Dispersion interaction is a result of the dynamic (weak) part of electron correlation. It is still interesting to investigate the performance of multi-reference (MR) methods (developed to give reasonable results for strongly correlated systems) when it comes to calculating weak interactions. This is the purpose of the present study.

MR methods assume a multi-determinantal reference function. The reference function we are concerned with in this work is the antisymmetrized product of strongly orthogonal geminals (APSG) [24]. The APSG wave function for an N-electron singlet system (N being even) is constructed using two-electron functions (geminals) with a handful of orbitals assigned to each of them. The strong orthogonality condition implies the mutually disjoint distribution of orbitals among the geminals [25]. The wave function takes a product form:

$$\begin{aligned} |{\text{APSG}}\rangle =\prod _{\mathcal{{P}}=1}^{N/2}\psi _{\mathcal{{P}}}^{\dagger }|{\text{vac}}\rangle , \end{aligned}$$

with

$$\begin{aligned} \psi _{\mathcal{{P}}}^{\dagger }=\sum _{p \in {\mathcal{P}}}c_{p}a_{p\alpha }^{\dagger }a_{p\beta }^{\dagger } . \end{aligned}$$

Here, $|{\text{vac}}\rangle$ is the true vacuum, $a_{p\sigma }^{\dagger }$ is the creation operator of spatial orbital $\phi _p$ with spin $\sigma$ and $\mathcal{{P}}$ is the index of the geminal $\phi _p$ belongs to. We arrive at a particularly simple wave function if each geminal consists of at most two orbitals; this choice defines the generalized valence bond (GVB) method [26]. We work on the basis of natural orbitals diagonalizing the one-particle reduced density matrix. The geminal coefficients obey the relation $|c_p|^2=n_p$, where $n_p \in [0,1]$ is the occupation number. For later use, we also introduce the hole occupation number ${\bar{n}}_p=1-n_p$. During an APSG calculation, geminal coefficients and orbitals are optimized simultaneously. The intrageminal part of electron correlation is taken into account, while the intergeminal interaction is treated in a mean-field manner. While APSG usually provides a good starting point for the description of left–right correlation, it is no better at calculating dispersion energy than HF. Various correction schemes applied to APSG can, however, tackle this problem. The APSG-based linearized CC of Zoboki et al. was shown to give very accurate results for a He dimer [27]. The extended RPA for APSG (APSG-ERPA) developed by Pernal et al. [28, 29] was also shown to account for a great portion of dispersion energy; later this approach was further refined by electron-groups embedding [30, 31].

We presented our approach to GVB-based correction methods in Ref. [32]: an MR extension of rCCD, formulated in the internally contracted (ic) MR-CC framework of Mukherjee et al. [33]. The method greatly improves on GVB results in various scenarios (bond breaking, conformational barriers, etc.) and gives results somewhat similar to ERPA, although the formal RPA-rCCD correspondence does not hold anymore [34]. In this paper, we assess GVB-rCCD on the example of calculating interaction energy between noble gas atoms.

The rest of the paper is organized as follows: in Sect. 2, we demonstrate that in leading order the spin-orbital GVB-rCCD energy shows the correct asymptotics, according to (1), and argue that the full rCCD energy correction shows a similar decay. In Sect. 3, a brief review of spin-free GVB-rCCD is given, and numerical results are presented for the interaction energy of He$_2$ and Ne$_2$. Conclusions are drawn in Sect. 4.

Atomic units are used throughout the paper.

2 Asymptotic behaviour of GVB-rCCD

Our ic-MR-CC formalism is based on the MR extension of the generalized normal ordering (denoted by  ::) and the corresponding generalized Wick theorem (MR-GWT) of Kutzelnigg and Mukherjee [35]. The correlated ground-state wave function is parametrized by a normal-ordered exponential wave operator, originally proposed by Lindgren [36]:

$$\begin{aligned} |\Psi \rangle =:\exp (T):|\Phi \rangle , \end{aligned}$$

(2)

where $|\Phi \rangle$ can be—in principle—any MR wave function, and

$$\begin{aligned} T=\frac{1}{4}\sum _{\substack{ I,A \\ J,B }}t^{IJ}_{AB} \, a^{AB}_{IJ} , \end{aligned}$$

with $a^{AB}_{IJ}=a_{A}^{\dagger }a_{B}^{\dagger }a_{J}a_{I}$. Indices I, J, K, L and A, B, C, D correspond to arbitrary core/active and active/virtual spin orbitals, respectively (excluding a few ill-defined cases; see Sect. 3). The working equations are obtained by plugging (2) in the Schrödinger equation and projecting with $\langle \Phi |$, and $\langle \Phi |:a^{IJ}_{AB}:$. The resulting formulae are to be evaluated using the MR-GWT, leading to an expansion in terms of density cumulants. Due to the incredibly large number of terms generated by the MR-GWT, it is virtually impossible to proceed without further approximations. A way of arriving at a manageable set of equations is given by the MR extension of the ring approximation. Following the same procedure as in Ref. [32], one arrives at the spin-orbital variant of the MR-rCCD equations:

$$\begin{aligned} \Delta E&=\frac{1}{4}\sum _{\substack{ I,A \\ J,B }}V_{IJ}^{AB} \, t^{IJ}_{AB} \, n_{I}n_{J}{\bar{n}}^{A}{\bar{n}}^{B} , \end{aligned}$$

(3)

$$\begin{aligned} 0&=V^{{I \,}{J \,}}_{{A \,} {B \,}} -\sum _{K}{\mathcal{A}}^{IJ}\Big (F^{{J}}_{K} \, t^{IK}_{AB}\Big )n_{K} +\sum _{C}{\mathcal{A}}_{AB}\Big (F^{C}_{{B}} \, t^{IJ}_{AC} \Big ){\bar{n}}^{C} \nonumber \\&\quad +\sum _{K,C} \Big (V^{{I}C}_{{A}K} \, t^{KJ}_{CB} +V^{{J}C}_{{B}K} \, t^{KI}_{CA}\Big ) n_{K}{\bar{n}}^{C} \nonumber \\&\quad +\sum _{\substack{ K,C \\ L,D }} t^{IK}_{AC} \, V^{CD}_{KL} \, t^{LJ}_{DA} \, n_{K}n_{L}{\bar{n}}^{C}{\bar{n}}^{D} , \end{aligned}$$

(4)

where ${\mathcal{A}}^{IJ}=I-P^{IJ}$ is the antisymmetrizing operator. The MR-rCCD equations (3) and (4) formally look just like their single-reference (SR) counterparts, but contain “dressed” one- and two-particle vertices F and V instead of “bare” vertices f (the Fockian) and ${\tilde{v}}$ (the antisymmetrized two-electron integral). The explicit form of these vertices is given in “Appendix 1”. The dressed vertices are functions of density cumulants and therefore contain a zeroth-order description of static electron correlation ($|\Phi \rangle$ being MR). If $|\Phi \rangle$ happens to be a single Slater determinant, then all cumulants vanish, the dressed vertices fall back to bare ones, and the usual SR-rCCD formulae are retrieved. Setting $|\Phi \rangle =|{\text{GVB}}\rangle$ defines GVB-rCCD.

Using the SR-like energy denominator

$$\begin{aligned} \epsilon ^{A}_{I}={\bar{n}}^{A}F^{A}_{{A}}-n_{I}F^{{I}}_{I} \end{aligned}$$

(5)

to solve (4) iteratively and plugging back the iterations into (3) leads to a PT-like expansion for the energy correction:

$$\begin{aligned} \Delta E=-\frac{1}{4}\sum _{\substack{ I,A \\ J,B }} \frac{V_{IJ}^{AB} \, V^{{I \,}{J \,}}_{{A \,}{B \,}}}{\epsilon ^{A}_{I}+\epsilon ^{B}_{J}} n_{I}{\bar{n}}^{A}n_{J}{\bar{n}}^{B} + \cdots \end{aligned}$$

(6)

Because of its formal similarity to a second-order PT energy, we shall refer to the first term on the rhs of (6) as $E^{(2)}$. The connectedness of the formulae and the vertices ensures the extensivity of the method, provided that intersystem cumulants vanish (which is always the case for closed-shell subsystems [37]).

The spin-orbital SR-rCCD energy correction is known to be somewhat “incomplete”: in its diagrammatic expansion, ring-type antisymmetrized Goldstone diagrams beyond second order contain an extra factor of 1/2. This factor would be cancelled by the crossed-ring terms (with, for example, amplitudes $t^{IK}_{CB}$) missing in rCCD. An easy way to make the diagrams complete (without double counting $E^{(2)}$) is through the modified energy formula [6, 38] $\Delta {\tilde{E}}=2\Delta E -E^{(2)}$. The same reasoning applies to our spin-orbital MR-rCCD energy; however, we shall not pursue this question any further.

Most RPA/rCCD-like methods do not give a satisfactory description of long-range dispersion interactions [2, 23]: while they produce the correct $\sim 1/R^6$ decay, the dispersion coefficient $C_6$ can rarely be expressed in terms of FDPs. We show here that the leading-order contribution of the GVB-rCCD energy exhibits the correct asymptotics and its dispersion coefficient can be expressed in terms of FDP-like quantities. We believe that this result gives valuable insight into how cancellations occur in the dressed vertices. We also argue that the $\sim 1/R^6$ decay holds for the complete GVB-rCCD energy as well, but without the correct form of $C_6$. The derivation follows closely the MP2 analysis of Szabo and Ostlund [2, 7] and can be thought of as a possible MR extension.

We start from $E^{(2)}$, the leading-order energy in (6), and consider singlet systems. After breaking up the spin-orbital indices to a spatial and a spin part (e.g. $A=a\sigma$), and using relations $n_i :=n_{i\alpha }=n_{i\beta }$ and (27), one finds that the energy denominator is spin-independent: $\epsilon ^{a}_{i} :=\epsilon ^{a\alpha }_{i\alpha }=\epsilon ^{a\beta }_{i\beta }$. The summation over the remaining spin degrees of freedom can be then easily performed with the help of (28). Utilizing the spin-free two-particle vertex ${\mathcal{V}}^{ab}_{ij}$ (given in “Appendix 1”), $E^{(2)}$ takes the form:

$$\begin{aligned} E^{(2)}=-\frac{1}{12}\sum _{i,j,a,b}\frac{{\mathcal{V}}_{ij}^{ab} \left( 2{\mathcal{V}}^{ij}_{ab} +{\mathcal{V}}^{ij}_{ba}\right) }{\epsilon ^a_i+\epsilon ^b_j}n_i{\bar{n}}^an_j{\bar{n}}^b . \end{aligned}$$

(7)

We now consider a supersystem composed of two closed shell, electrically neutral atoms ${\mathbb {A}}$ and ${\mathbb {B}}$ separated at a large distance R. Orbitals are assumed to be localized on either ${\mathbb {A}}$ or ${\mathbb {B}}$. Based on this, indices in (7) can be categorized, resulting in two intrasystem terms ($i,j,a,b \in {\mathbb {A}}$ or $i,j,a,b \in {\mathbb {B}}$) and the remaining intersystem term:

$$\begin{aligned} E^{(2)}=E^{(2)}_{{\mathbb {A}}}+E^{(2)}_{{\mathbb {B}}}+E^{(2)}_{disp} . \end{aligned}$$

Turning to the dispersion energy, the differential overlap of the subsystems is exponentially small, making most of the intersystem integrals negligible. Intersystem cumulants also vanish due to the strongly orthogonal geminals being localized on subsystems. Therefore, cases like $(i,j \in {\mathbb {A}} ; \ a,b \in {\mathbb {B}})$ or $(i,j,a \in {\mathbb {A}} ; \ b \in {\mathbb {B}})$ can be safely neglected. The only remaining cases are of the form

$$\begin{aligned}&(i,a \in {\mathbb {A}} ; \ j,b \in {\mathbb {B}}) \ \rightarrow \ {\mathcal{V}}^{ab}_{ij} \approx 4W^{ab}_{ij} , \end{aligned}$$

(8)

$$\begin{aligned}&(j,a \in {\mathbb {A}} ; \ i,b \in {\mathbb {B}}) \ \rightarrow \ {\mathcal{V}}^{ab}_{ij} \approx -2W^{ab}_{ji} . \end{aligned}$$

(9)

The expression of vertex $W^{ab}_{ij}$ is given as

$$\begin{aligned} W^{ab}_{ij}&=v^{ab}_{ij}+\frac{1}{2}\sum _{p,q}\left( \frac{1}{n_i{\bar{n}}^a} \ v^{bq}_{jp}\Lambda ^{ap}_{iq} +\frac{1}{n_j{\bar{n}}^b} \ v^{aq}_{ip}\Lambda ^{bp}_{jq}\right) \nonumber \\&\quad +\frac{1}{4}\frac{1}{n_i{\bar{n}}^an_j{\bar{n}}^b}\sum _{p,q,r,s} v^{qs}_{pr} \ \Lambda ^{ap}_{iq}\Lambda ^{br}_{js} \nonumber \\&=\int {\mathrm {d}}^3r_1\int {\mathrm {d}}^3r_2\frac{\Omega ^a_i({\varvec{r}}_1)\Omega ^b_j({\varvec{r}}_2)}{|{\varvec{r}}_1-{\varvec{r}}_2|} , \end{aligned}$$

(10)

and

$$\begin{aligned} \Omega ^a_i({\varvec{r}})=\phi _i^{*}({\varvec{r}})\phi _a({\varvec{r}}) +\frac{1}{2n_i{\bar{n}}^a}\sum _{p,q}\Lambda ^{ap}_{iq} \phi _p^{*}({\varvec{r}})\phi _q({\varvec{r}}) . \end{aligned}$$

Using the partial trace relation [39, 40]

$$\begin{aligned} \sum _{p}\Lambda ^{ap}_{ip}=-2\delta ^{a}_{i}n_i{\bar{n}}^a , \end{aligned}$$

it is easy to see that $\Omega ^a_i({\varvec{r}})$ integrates to zero for all i, a:

$$\begin{aligned} \int {\mathrm {d}}^3r \, \Omega ^a_i({\varvec{r}})=\delta ^a_i+\frac{1}{2n_i{\bar{n}}^a} \sum _{p}\Lambda ^{ap}_{ip}=0 . \end{aligned}$$

(11)

Quantities ${W_{ab}^{ij}}$, ${\Omega _a^i({\varvec{r}})}$ appearing in the asymptotic limit of ${\mathcal{V}}_{ab}^{ij}$ can be obtained analogously to the above. Based on “Appendix 1”, the relation of ${\mathcal{V}}_{ab}^{ij}$ and ${\mathcal{V}}^{ab}_{ij}$ is Hermitian conjugation, and the same holds for ${W_{ab}^{ij}}$ and ${\Omega _a^i({\varvec{r}})}$.

Substituting (8), (9) and their Hermitian conjugate into (7) gives

$$\begin{aligned} E^{(2)}_{disp} \approx -4\sum _{\substack{ i,a \in {\mathbb {A}} \\ j,b \in {\mathbb {B}} }} \frac{W_{ij}^{ab}W^{{i \,} {j \,}}_{{a \,} {b \,}}}{\epsilon ^a_i+\epsilon ^b_j}n_i{\bar{n}}^an_j{\bar{n}}^b . \end{aligned}$$

(12)

We shall now make use of the multipole expansion of the Coulomb potential:

$$\begin{aligned} \frac{1}{|{\varvec{r}}_1-{\varvec{r}}_2|}= \frac{1}{|{\varvec{r}}-{\varvec{R}}|} =\frac{1}{R}+\frac{{\varvec{n}}\cdot {\varvec{r}}}{R^2} +\frac{3({\varvec{n}}\cdot {\varvec{r}})^2-r^2}{2R^3} +{\mathcal{O}}(R^{-4}) , \end{aligned}$$

(13)

with ${\varvec{n}}$ being a unit vector parallel to ${\varvec{R}}={\varvec{R}}_{{\mathbb {B}}}-{\varvec{R}}_{{\mathbb {A}}}$ and ${\varvec{r}}={\varvec{r}}_{{\mathbb {A}}1}-{\varvec{r}}_{{\mathbb {B}}2}$ being the difference of electronic coordinates measured from the nuclei.

Inserting (13) into (10), the appearance of anomalous terms proportional to 1/R or $1/R^2$ in $W^{ab}_{ij}$ is prevented by (11). Keeping only the first non-vanishing term yields

$$\begin{aligned} W^{ab}_{ij} \approx \frac{1}{R^3}\varvec{\mu }^a_i\varvec{\tau }\varvec{\mu }^b_j , \end{aligned}$$

(14)

where $\varvec{\tau }={\varvec{I}}-3{\varvec{n}}\otimes {\varvec{n}}$ is a symmetric matrix ($\otimes$ being the dyadic product of the 3-vectors) and $\varvec{\mu }^a_i$ can be interpreted as the transition dipole moment with a cumulant correction:

$$\begin{aligned} \varvec{\mu }^a_i=-\int {\mathrm {d}}^3r' \, \Omega ^a_i({\varvec{r}}'){\varvec{r}}'= \frac{1}{2n_i{\bar{n}}^a}\left\langle {\text{GVB}}|\varvec{\mu }:E_{i}^{a}:|{\text{GVB}}\right\rangle . \end{aligned}$$

(15)

In the above, $E^a_i=a_{i\alpha }^{a\alpha }+a_{i\beta }^{a\beta }$ is a spin-free excitation operator. As before, the cumulant-corrected transition dipole ${\varvec{\mu }_b^j}$ is the Hermitian conjugate of ${\varvec{\mu }^b_j}$.

Substituting (14) and its Hermitian conjugate into (12), we have ${E^{(2)}_{{\rm disp}} \approx -C_6/R^6}$, with the coefficient $C_6$ given by

$$\begin{aligned} C_6=4\sum _{\substack{ i,a, \in {\mathbb {A}} \\ j,b \in {\mathbb {B}} }} \frac{n_i{\bar{n}}^an_j{\bar{n}}^b}{\epsilon ^a_i+\epsilon ^b_j} (\varvec{\mu }^a_i\varvec{\tau }\varvec{\mu }^b_j) (\varvec{\mu }^{{i}}_{{a}}\varvec{\tau }\varvec{\mu }^{{j}}_{{b}}) . \end{aligned}$$

Assuming the positivity of $\epsilon _i^a$ (as argued in “Appendix 3”), we can use the identity

$$\begin{aligned} \frac{1}{\epsilon +\epsilon '}=\frac{2}{\pi }\int \limits _0^{\infty }{\mathrm {d}}\omega \frac{\epsilon \epsilon '}{(\epsilon ^2+\omega ^2)(\epsilon '^2+\omega ^2)} \ \ \ \ (\epsilon ,\epsilon '>0) , \end{aligned}$$

to obtain the final form of $C_6$, reading as

$$\begin{aligned} C_6=\int \limits _0^{\infty }\frac{{\mathrm {d}}\omega }{2\pi } \mathrm {Tr}\left[ \varvec{\alpha }_{\mathbb {A}}(i\omega )\varvec{\tau } \varvec{\alpha }_{\mathbb {B}}(i\omega )\varvec{\tau }\right] , \end{aligned}$$

(16)

where orbitals were taken to be real for simplicity, and we introduced the FDP-like quantities:

$$\begin{aligned} \varvec{\alpha }_{\mathbb {A}}(i\omega )=4\sum _{i,a \in {\mathbb {A}}} \frac{\epsilon ^a_in_i{\bar{n}}^a}{(\epsilon ^a_i)^2-(i\omega )^2} \varvec{\mu }^a_i\otimes \varvec{\mu }^{{i}}_{{a}} . \end{aligned}$$

For isotropic subsystems, $\varvec{\alpha }_{\mathbb {A}}(i\omega )$ must be proportional to the unit tensor: $\varvec{\alpha }_{\mathbb {A}}(i\omega )=\alpha _{\mathbb {A}}(i\omega ){\varvec{I}}$. Using this together with $\mathrm {Tr}[\varvec{\tau }^2]=6$ allows us to write (16) as

$$\begin{aligned} C_6=\frac{3}{\pi }\int \limits _0^{\infty }{\mathrm {d}}\omega \, \alpha _{\mathbb {A}}(i\omega )\alpha _{\mathbb {B}}(i\omega ) . \end{aligned}$$

The obtained expressions are formally similar to the PT-based formulae of Szabo and Ostlund [2], but involve cumulant-dependent quantities, such as $\varvec{\mu }^a_i$ or $\varvec{\alpha }_{\mathbb {A}}(i\omega )$.

As it was demonstrated in this section, dressed vertices show the proper asymptotic behaviour, just like their bare counterparts. Based on this, one could continue the analysis and investigate the asymptotics of the full rCCD equations (3) and (4); while the $\sim 1/R^6$ decay could be argued this way, the coefficient $C_6$ would not exhibit the form of (16) involving FDPs.

3 Numerical results for noble gas dimers in spin-free formalism

Numerical applications presented in this work adopt a spin-free formulation of GVB-rCCD facilitated by the Unitary Group Approach (UGA) [41], as described in Ref. [32]. The SR ring approximation is not invariant to the various alternative treatments of spin during the derivation of equations. The situation is analogous in the MR case, spin-orbital and UGA parametrizations of GVB-rCCD being non-equivalent. The major findings of Sect. 2, however, hold for the UGA case as well, the only difference being that the second-order approximation of the dispersion energy is 5/4 times that of (12). In the following, we give a brief account of GVB-rCCD(UGA) used for obtaining the pilot results. The cluster operator is parametrized as

$$\begin{aligned} T=\frac{1}{2}\sum _{i,j,a,b}t^{ij}_{ab} \ \sum _{\sigma ,\sigma '}a^{a\sigma b\sigma '}_{i\sigma j\sigma '} =\frac{1}{2}\sum _{i,j,a,b}t^{ij}_{ab}E_{ij}^{ab} , \end{aligned}$$

(17)

with $t^{ij}_{ab}=t^{ji}_{ba}$. Indices i, j, k, l and a, b, c, d belong to core/active and active/virtual spatial orbitals, respectively. The energy formula reads:

$$\begin{aligned} \Delta E=\sum _{i,j,a,b}B_{ij}^{ab} \ t^{ij}_{ab} \, n_in_j{\bar{n}}^a{\bar{n}}^b , \end{aligned}$$

(18)

and the amplitude equation takes the form:

$$\begin{aligned} 0 =B^{{i \,} {j \,}}_{{a \,} {b \,}} +2\sum _{k,c}\left( t^{ik}_{ac}A^{cj}_{kb}+A^{ic}_{ak}t^{kj}_{cb}\right) +4\sum _{ k,c, l,d }t^{ik}_{ac}B^{cd}_{kl}t^{lj}_{db} \, n_kn_l{\bar{n}}^c{\bar{n}}^d , \end{aligned}$$

(19)

where

$$\begin{aligned} B^{ab}_{ij}&= \frac{1}{2}{\mathcal{V}}_{ij}^{ab} , \\ A^{aj}_{ib}&= \frac{1}{2}\Big ({\mathcal{F}}_{{\, b}}^{\, a} \ \delta ^j_i \, {\bar{n}}^a -{\mathcal{F}}^{{\, j}}_{\, i} \ \delta ^a_b \, n_i +{\mathcal{V}}^{a{j}}_{i{b}} \, n_i{\bar{n}}^a\Big ) . \end{aligned}$$

The algebraic expressions of spin-free vertices are given in “Appendix 1”.

There are two technical issues to address when solving (19). The first is the redundancy of excited states used for projection, which is presently handled based on the concept of frames [42, 43]. The associated numerical threshold for identifying zero eigenvalues of the overlap matrix is set to $s_{\min}=10^{-10}$. The second issue concerns the removal of certain kinds of excitations which could otherwise make the vertices and thus the amplitudes singular. To prevent this, an excitation selection scheme has been set up: excitation $E_{ij}^{ab}$ is allowed in (17) only if

• (i) $n_i{\bar{n}}^a >1/4 \,$ for $i \ne a$ and ${\mathcal{I}}={\mathcal{A}}$ ,

• (ii) $n_i{\bar{n}}^b >1/4 \,$ for $i \ne b$ and ${\mathcal{I}}=\mathcal{{B}}$ ,

• (iii) $n_j{\bar{n}}^a >1/4 \,$ for $j \ne a$ and ${\mathcal{J}}={\mathcal{A}}$ ,

• (iv) $n_j{\bar{n}}^b >1/4 \,$ for $j \ne b$ and ${\mathcal{J}}=\mathcal{{B}}$ ,

• (v) $n_i >1/2 \,$ for $i = j$ ,

• (vi) ${\bar{n}}^a >1/2 \,$ for $a = b$ ,

calligraphic letters denoting geminal indices. All other excitations (including spectators, such as $E^{aj}_{ij}$) are allowed. The above constraints can be introduced in the results of Sect. 2 without any problem. Restricting the excitations this way does not violate size consistency.

GVB-rCCD(UGA) dispersion energy curves are presented for He$_2$ and Ne$_2$. The results are compared to the GVB-ERPA method of Pernal et al. [28, 29], and HF-CCSD(T) is used as benchmark. Calculations based on the GVB reference were performed by the Budapest version of the MUNGAUSS program package [44], and HF-CCSD(T) was calculated with the Gaussian 09 software [45]. Calculations were carried out using the Dunning aug-cc-pVDZ atomic basis set [46] (with five d-orbitals).

All calculations assume GVB natural orbitals, with geminals consisting of two orbitals, except for core geminals (having only one orbital).

In reporting the results, $E_{{\rm disp}}$ stands for the total energy of the supersystem minus total energies of the subsystems, computed according to the Boys–Bernardi scheme [47]. Total energy is calculated as the sum of the GVB energy and $\Delta E$ of (18).

Equilibrium atomic distances and the corresponding dispersion energies are given in Table 1, while Figs. 1 and 2 show the dispersion energy curves for He$_2$ and Ne$_2$, respectively. Just like HF, GVB cannot describe dispersion interactions, resulting in non-bonding energy curves. GVB-ERPA and GVB-rCCD(UGA) improve significantly on this situation, although the results are far from perfect. While GVB-ERPA underestimates the dispersion energy in both cases, GVB-rCCD(UGA) overestimates it, although to a somewhat lesser extent. A more significant difference shows up in the equilibrium atomic distances: GVB-ERPA overshoots by 0.4 bohr in each case, while GVB-rCCD(UGA) reproduces the benchmark CCSD(T) value to the digits displayed in Table 1.

Table 1 Equilibrium atomic distances and interaction energies $E_{{\rm int}}=E_{{\rm disp}}(R_{\min})$ for He$_2$ and Ne$_2$, calculated on the aug-cc-pVDZ basis set

Fig. 1

Dispersion energy of He$_2$, calculated on the aug-cc-pVDZ basis set

Fig. 2

Dispersion energy of Ne$_2$, calculated on the aug-cc-pVDZ basis set

Deficiencies of GVB-ERPA have been noted before in describing weak interaction, and electron-groups embedding has been suggested as a cure (GVB-EERPA) [30, 31]. It was also demonstrated that GVB-EERPA represents a significant improvement over GVB-ERPA in such situations. Comparison with GVB-ERPA is warranted here, since this unamended version represents the theory paralleling GVB-rCCD which is of prime interest in this study.

4 Conclusions

A recently proposed GVB-rCCD method built with dressed vertices (20) and (21) has been investigated concerning the calculation of dispersion interactions of noble gas dimers. An analytical inspection has revealed that the GVB-rCCD dispersion energy follows the correct $\sim 1/R^6$ asymptotic decay, and in the leading order of the energy, the dispersion coefficient can be expressed in terms uncoupled FDPs that involve a dressing taking into account the MR character of the reference function. Pilot numerical tests of full GVB-rCCD result in dispersion energies and equilibrium atomic distances in fairly good agreement with benchmark data.

Analysis of the $R\rightarrow \infty$ asymptotics of the full GVB-rCCD scheme is warranted as a continuation of this study since an expression of the form of (16) may arise with improved FDPs [2, 15], depending on the incorporation of exchange terms in the equations. Our current formulation is not favourable in this respect. Investigation in the line of GVB-rCCD variants yielding a $C_6$ expressed with coupled FDPs is deferred to a follow-up study.

Improving upon the correlation content of the equations, i.e. stepping beyond the ring approximation, is a further line of research for the future. This can be done by a systematic inclusion of a couple of types of non-ring diagrams, expressed with dressed vertices. Derivation and testing of these more involved MR-CC schemes shall be the subject of a forthcoming paper.

Appendices

Appendix 1: Vertices appearing in the MR-rCCD equations

The bare vertices are $f_{p\sigma }^{q\sigma '}=f^{p}_{q}\delta ^{\sigma '}_{\sigma }$ and ${\tilde{v}}^{q\sigma ' s\sigma '''}_{p\sigma r\sigma ''}= v^{qs}_{pr}\delta ^{\sigma '}_{\sigma }\delta ^{\sigma '''}_{\sigma ''}- v^{qs}_{rp}\delta ^{\sigma '}_{\sigma ''}\delta ^{\sigma '''}_{\sigma }$, with

$$\begin{aligned} f^q_p&= h^q_p+\sum _r{\bar{v}}^{qs}_{pr} \, n_r , \\ h^q_p&= \int {\mathrm {d}}^3r \, \phi _p^{*}({\varvec{r}})\left( -\frac{1}{2}\nabla ^2 -\sum _{{\mathbb {A}}}\frac{Z_{{\mathbb {A}}}}{|{\varvec{r}}-{\varvec{R}}_{{\mathbb {A}}}|}\right) \phi _q({\varvec{r}}) , \\ v^{qs}_{pr}&= \int {\mathrm {d}}^3r_1\int {\mathrm {d}}^3r_2 \, \phi _p^{*}({\varvec{r}}_1)\phi _r^{*}({\varvec{r}}_2)\frac{1}{|{\varvec{r}}_1-{\varvec{r}}_2|} \phi _s({\varvec{r}}_2)\phi _q({\varvec{r}}_1) , \\ {\bar{v}}^{qs}_{pr}&= 2v^{qs}_{pr}-v^{qs}_{rp} . \end{aligned}$$

Indices p, q, r, s correspond to arbitrary spatial orbitals.

The spin-orbital dressed vertices read at the one-particle level:

$$\begin{aligned}\zeta _{U}\zeta ^{X}F_{U}^{X}& =\zeta _{U}\zeta ^{X}f^{X}_{U} -\sum _{P,Q}f^{Q}_{P} \ \lambda ^{PX}_{QU}\nonumber \\& -\frac{1}{2}\sum _{P,Q,R}\Big ( \zeta ^{X}{\tilde{v}}^{XQ}_{PR} \ \lambda ^{PR}_{UQ} +\zeta _{U}{\tilde{v}}_{UQ}^{PR} \ \lambda _{PR}^{XQ}\Big ) \nonumber \\& -\frac{1}{4}\sum _{P,Q,R,S} {\tilde{v}}^{QS}_{PR} \ \lambda ^{PRX}_{QSU} , \end{aligned}$$

(20)

and at the two-particle level:

$$\begin{aligned}\zeta _{U}{\zeta }^{X}\zeta _{W}{\zeta }^{Y}V_{UW}^{XY}& =\zeta _{U}{\zeta }^{X}\zeta _{W}{\zeta }^{Y}{\tilde{v}}_{UW}^{XY} +\sum_{P}\Big [{\mathcal{A}}^{XY}\Big ({\zeta }^{X}f_{P}^{X} \lambda _{UW}^{PY}\Big ) +{\mathcal{A}}_{UW}\Big (\zeta _{U}f_{U}^{P}\lambda _{PW}^{XY}\Big )\Big ] \nonumber \\& +\frac{1}{2}\sum_{P,Q}\Big ({\zeta }^{X}{\zeta }^{Y} {\tilde{v}}_{PQ}^{XY}\lambda _{UW}^{PQ} +\zeta _{U}\zeta _{W}{\tilde{v}}_{UW}^{PQ}\lambda _{PQ}^{XY}\Big ) \nonumber \\& -\sum _{P,Q}{\mathcal{A}}_{UW}{\mathcal{A}}^{XY}\Big (\zeta _{U}{\zeta }^{X} {\tilde{v}}_{UP}^{XQ}\lambda _{WQ}^{YP}\Big ) \nonumber \\& +\frac{1}{4}\sum _{P,Q,R,S} {\tilde{v}}_{PR}^{QS}\lambda _{UW}^{PR}\lambda _{QS}^{XY} +\sum _{P,Q,R,S}{\tilde{v}}_{PR}^{QS} {\mathcal{A}}_{UW}\Big (\lambda _{QU}^{PX}\lambda _{SW}^{RY}\Big ) \nonumber \\& -\frac{1}{2}\sum _{P,Q,R,S}{\tilde{v}}_{PR}^{QS} \Big [{\mathcal{A}}^{XY}\Big (\lambda _{QS}^{PX}\lambda _{UW}^{RY}\Big ) +{\mathcal{A}}_{UW}\Big (\lambda _{QU}^{PR}\lambda _{SW}^{XY}\Big )\Big ] \nonumber \\& +\sum _{P,Q}f_{P}^{Q}\lambda _{UWQ}^{XYP} +\frac{1}{2}\sum _{P,Q,R}\Big [{\mathcal{A}}^{XY}\Big ({\zeta }^{X} {\tilde{v}}_{PR}^{QX}\lambda _{UWQ}^{RYP}\Big ) +{\mathcal{A}}_{UW}\Big (\zeta _{U}{\tilde{v}}_{UP}^{RQ} \lambda _{RWQ}^{XYP}\Big )\Big ] \nonumber \\& +\frac{1}{4}\sum _{P,Q,R,S}{\tilde{v}}_{PR}^{QS} \lambda _{UWQS}^{XYPR} . \end{aligned}$$

(21)

Indices P, Q, R, S correspond to arbitrary spin orbitals, and $\lambda$ denotes the spin-orbital cumulants [35, 39]. In this paper, we distinguish between indices corresponding to core/active (I, J, K, L) and active/virtual (A, B, C, D) orbitals. Different dressed vertices belong to the different combination of indices, distinguished by the values $\zeta$ takes:

$$\begin{aligned} \zeta _U=\left\{ \begin{array}{l} -n_U \ \ \ \mbox{if} \ \ \ U \in \{I,J,K,L\} \\ \ \ {\bar{n}}_U \ \ \ \ \mbox{if} \ \ \ U \in \{A,B,C,D\} . \\ \end{array} \right. \end{aligned}$$

(22)

The spin summation of (20) and (21) yields the spin-free dressed vertices:

$$\begin{aligned} \zeta _u{\zeta }^x{\mathcal{F}}_u^x&= \zeta _u{\zeta }^x\sum _{\sigma }F_{u\sigma }^{x\sigma } \nonumber \\&=2\zeta _u{\zeta }^xf^x_u-\sum _{p,q}f^q_p\Lambda ^{px}_{qu}\nonumber \\&\quad -\sum _{p,q,r}\Big ({\zeta }^xv^{xq}_{pr}\Lambda ^{pr}_{uq}+\zeta _uv_{uq}^{pr}\Lambda _{pr}^{xq}\Big ) -\frac{1}{2}\sum _{p,q,r,s}v^{qs}_{pr}\Lambda ^{prx}_{qsu} , \end{aligned}$$

(23)

and

$$\begin{aligned}&\zeta _u{\zeta }^x\zeta _w{\zeta }^y{\mathcal{V}}_{uw}^{xy}\nonumber \\&\quad = \zeta _u{\zeta }^x\zeta _w{\zeta }^y\sum _{\sigma ,\sigma '}V_{u\sigma w\sigma '}^{x\sigma y\sigma '} \nonumber \\&\quad =2\zeta _u{\zeta }^x\zeta _w{\zeta }^y{\bar{v}}_{uw}^{xy} +\sum _{p}\Big ({\zeta }^xf_{p}^{x}\Lambda _{uw}^{py}+\zeta _uf_{u}^{p}\Lambda _{pw}^{xy} +{\zeta }^yf_{p}^{y}\Lambda _{uw}^{xp}\nonumber \\&\qquad +\zeta _wf_{w}^{p}\Lambda _{up}^{xy}\Big ) \nonumber \\&\qquad +\sum _{p,q}\Big ({\zeta }^x{\zeta }^yv_{pq}^{xy}\Lambda _{uw}^{pq}+\zeta _u\zeta _wv_{uw}^{pq}\Lambda _{pq}^{xy} -\zeta _w{\zeta }^y{\bar{v}}_{wp}^{yq}\Lambda _{uq}^{xp}\nonumber \\&\qquad -\zeta _u{\zeta }^x{\bar{v}}_{up}^{xq}\Lambda _{wq}^{yp}\Big ) \nonumber \\&\qquad +\sum _{p,q}\Big [\zeta _w{\zeta }^x\Big (v_{wp}^{qx}\Lambda _{uq}^{py}+v_{wp}^{xq}\Lambda _{uq}^{yp}\Big ) \nonumber \\&\qquad +\zeta _u{\zeta }^y\Big (v_{up}^{qy}\Lambda _{wq}^{px}+v_{up}^{yq}\Lambda _{wq}^{xp}\Big )\Big ] \nonumber \\&\qquad +\frac{1}{2}\sum _{p,q,r,s}\Big [{\bar{v}}_{pr}^{qs} \Lambda _{uq}^{xp}\Lambda _{ws}^{yr} -v_{pr}^{qs}\Big ( \Lambda _{uw}^{py}\Lambda _{qs}^{xr} +\Lambda _{uw}^{xp}\Lambda _{qs}^{yr} +\Lambda _{us}^{pr}\Lambda _{qw}^{xy} \nonumber \\&\qquad +\Lambda _{uq}^{xy}\Lambda _{ws}^{pr} +\Lambda _{uq}^{yp}\Lambda _{ws}^{xr} \Big )\Big ] \nonumber \\&\qquad +\frac{1}{6}\sum _{p,q,r,s}v_{pr}^{qs}\Big [ \Lambda _{uw}^{pr}\Big (2\Lambda _{qs}^{xy}+\Lambda _{sq}^{xy}\Big ) +\Lambda _{us}^{yp}\Big (2\Lambda _{wq}^{xr}+\Lambda _{wq}^{rx}\Big ) \nonumber \\&\qquad +\Lambda _{us}^{py}\Big (2\Lambda _{wq}^{rx}+\Lambda _{wq}^{xr}\Big ) \Big ] \nonumber \\&\qquad +\sum _{p,q}f_{p}^{q}\Lambda _{uwq}^{xyp} +\sum _{p,q,r}\Big ( {\zeta }^xv_{pr}^{qx}\Lambda _{uwq}^{ryp} +\zeta _uv_{up}^{rq}\Lambda _{rwq}^{xyp} \nonumber \\&\qquad +{\zeta }^yv_{pr}^{qy}\Lambda _{uwq}^{xrp} +\zeta _wv_{wp}^{rq}\Lambda _{urq}^{xyp} \Big ) \nonumber \\&\qquad +\frac{1}{2}\sum _{p,q,r,s}v_{pr}^{qs}\Lambda _{uwqs}^{xypr} . \end{aligned}$$

(24)

Function $\zeta$ in the above is the analogue of (22), with all orbital indices in lower case. Spin-free cumulants are denoted by $\Lambda$; see “Appendix 2” for the corresponding formulae.

In order to make the arguments of “Appendix 3” easier to follow, we give the explicit formulae for ${\mathcal{F}}$ with both indices being either core/active or active/virtual type:

$$\begin{aligned} n_in^j{\mathcal{F}}_{ \, i}^{{\, j}} &=2n_in^jf^j_i-\sum _{p,q}f^q_p\Lambda ^{pj}_{qi} +\sum _{p,q,r}\Big (n^jv^{jq}_{pr}\Lambda ^{pr}_{iq}+n_iv_{iq}^{pr}\Lambda _{pr}^{jq}\Big )\nonumber \\&\quad -\frac{1}{2}\sum _{p,q,r,s}v^{qs}_{pr}\Lambda ^{prj}_{qsi} , \end{aligned}$$

(25)

$$\begin{aligned} {\bar{n}}_b{\bar{n}}^a{\mathcal{F}}_{{b}}^a &=2{\bar{n}}_b{\bar{n}}^af^a_b-\sum _{p,q}f^q_p\Lambda ^{pa}_{qb} -\sum _{p,q,r}\Big ({\bar{n}}^av^{aq}_{pr}\Lambda ^{pr}_{bq}+{\bar{n}}_bv_{bq}^{pr}\Lambda _{pr}^{aq}\Big )\nonumber \\&\quad -\frac{1}{2}\sum _{p,q,r,s}v^{qs}_{pr}\Lambda ^{pra}_{qsb} . \end{aligned}$$

(26)

For singlet systems, the relation between spin-orbital tensors invariant under spin rotations and their spin-free counterparts was established by Shamasundar (originally discussed in the context of density matrices and cumulants) [40, 48]. The same kind of formulae is applied to the dressed vertices as well:

$$\begin{aligned} F^{a\alpha }_{i\alpha }=\frac{1}{2}{\mathcal{F}}^a_i , \end{aligned}$$

(27)

and

$$\begin{aligned} V^{a\alpha b\alpha }_{i\alpha j\alpha }=\frac{1}{6}\left( {\mathcal{V}}^{ab}_{ij}-{\mathcal{V}}^{ab}_{ji}\right) , \ \ V^{a\alpha b\beta }_{i\alpha j\beta }=\frac{1}{6}\left( {2\mathcal {V}}^{ab}_{ij}+{\mathcal{V}}^{ab}_{ji}\right) . \end{aligned}$$

(28)

The vertices are invariant to spin-flip, and they vanish in all spin non-conserving cases.

Appendix 2: The structure of spin-free GVB cumulants

The calculation of dressed vertices (23)–(26) requires cumulants up to rank four. Computation and storage of such a cumulant (having eight indices running through all active orbitals) would be computationally too demanding, for example, for a CAS-type reference function. Strongly orthogonal geminal functions, on the other hand, offer a relatively simple and economic way to generate cumulants. Strong orthogonality implies the vanishing of all cumulant elements with indices belonging to different geminals [39]: for GVB, this condition renders the maximal number of nonzero elements in $\Lambda _K$ to only $N2^{2K-1}$. Because of strong orthogonality, cumulant terms in (23)–(26) are actually less expensive to calculate than SR terms, leading to GVB-rCCD exhibiting the same computational scaling as SR-rCCD.

Here, we give the spin-free expressions of GVB cumulants up to rank three:

$$\begin{aligned}&\Lambda ^{p}_{q}=2\delta ^{p}_{q} n_{p} , \nonumber \\&\Lambda ^{pr}_{qs} =2\delta ^{{\mathcal{P}}}_{{\mathcal{Q}}}c^pc_q\delta ^{pr}\delta _{qs} -\delta ^{{\mathcal{P}}{\mathcal{R}}}\left( \Lambda ^p_q\Lambda ^r_s -\frac{1}{2}\Lambda ^p_s\Lambda ^r_q\right) , \end{aligned}$$

(29)

$$\begin{aligned}\Lambda ^{prt}_{qsu}&=-\delta ^{\mathcal{{PR}}}\delta ^{\mathcal{{RT}}}\delta ^{\mathcal{{TP}}} \Big (\Lambda ^p_q\Lambda ^r_s\Lambda ^t_u -\frac{1}{2}\Lambda ^p_q\Lambda ^r_u\Lambda ^t_s +\frac{1}{4}\Lambda ^p_u\Lambda ^r_q\Lambda ^t_s\nonumber \\&\quad -\frac{1}{2}\Lambda ^p_u\Lambda ^r_s\Lambda ^t_q +\frac{1}{4}\Lambda ^p_s\Lambda ^r_u\Lambda ^t_q -\frac{1}{2}\Lambda ^p_s\Lambda ^r_q\Lambda ^t_u\Big )- \nonumber \\&\quad -\delta ^{\mathcal{{PR}}}\Big (\Lambda ^p_q\Lambda ^{rt}_{su} -\frac{1}{2}\Lambda ^p_s\Lambda ^{rt}_{qu}-\frac{1}{2}\Lambda ^p_u\Lambda ^{rt}_{sq}\Big )\nonumber \\&\quad -\delta ^{\mathcal{{RT}}}\Big (\Lambda ^r_s\Lambda ^{tp}_{uq} -\frac{1}{2}\Lambda ^r_u\Lambda ^{tp}_{sq}-\frac{1}{2}\Lambda ^r_q\Lambda ^{tp}_{us}\Big ) \nonumber \\&\quad -\delta ^{\mathcal{{TP}}}\Big (\Lambda ^t_u\Lambda ^{pr}_{qs} -\frac{1}{2}\Lambda ^t_q\Lambda ^{pr}_{us}-\frac{1}{2}\Lambda ^t_s\Lambda ^{pr}_{qu}\Big ) . \end{aligned}$$

(30)

The expression of the rank four cumulant ($\Lambda _4$) is too lengthy to be presented; it can be, however, easily manipulated with the aid of symbolic algebra codes.

Appendix 3: On the positivity of $\epsilon ^a_i$

In Sect. 2, when expressing $C_6$ in terms of FDP-like quantities (see (16)) we assumed the positivity of $\epsilon ^a_i$ (defined in (5)), a statement supported by numerical results for GVB, tested in various (e.g. bond breaking) scenarios. While we cannot prove this statement in an exact way, we can provide some formal rationalization based on certain reasonable assumptions concerning the diagonal elements of the GVB Fockian: $f^p_p$ is assumed to be negative if orbital p is strongly occupied ($n_p \lesssim 1$) and positive (or negative with a very small absolute value) if p is weakly occupied ($n_p \gtrsim 0$). Furthermore, GVB geminals consist of two orbitals, here denoted as $i,i'\in {\mathcal{I}}$ and $a,a'\in {\mathcal{A}}$, so that $n_i+n_{i'}=n_a+n_{a'}=1$. With the help of (27), $\epsilon ^a_i$ of (5) can be written equivalently as

$$\begin{aligned} \epsilon ^a_i=\frac{1}{2}\left( {\bar{n}}^a{\mathcal{F}}^a_{{a}}-n_i{\mathcal{F}}_i^{{i}}\right) . \end{aligned}$$

Using the explicit form of the GVB $\varvec{\Lambda }_2$ and $\varvec{\Lambda }_3$ given by (29) and (30), one can express (25) and (26) as

$$\begin{aligned}\frac{1}{2}n_i{\mathcal{F}}_i^{{i}} &=(2n_i-1)f_i^i+2{\bar{n}}_if_{i'}^{i'}+2n_i{\bar{n}}_iv_{ii}^{ii}\\&\quad -2\sqrt{n_i{\bar{n}}_i}v^{ii}_{i'i'}-2n_i{\bar{n}}_i{\bar{v}}^{ii'}_{ii'} -{\bar{n}}_i(2n_i-1)({\bar{v}}^{ii'}_{ii'}-v^{i'i'}_{i'i'}) , \\ \frac{1}{2}{\bar{n}}^a{\mathcal{F}}^a_{{a}} &=-(2n^a-1)f_a^a+2n_af_{a'}^{a'}-2n_a{\bar{n}}_av^{aa}_{aa}\\&\quad +2\sqrt{n_a{\bar{n}}_a}v^{aa}_{a'a'}+2n_a{\bar{n}}_a{\bar{v}}^{aa'}_{aa'} -n_a(2n_a-1)({\bar{v}}^{aa'}_{aa'}-v^{a'a'}_{a'a'}) . \end{aligned}$$

In the above, we made use of the fact that the coefficients of orbitals i and $i'$ have opposite signs (the same applying to a and $a'$). Real orbitals are assumed.

We shall analyse these expressions in two limiting cases: the HF-like case (one orbital in the geminal is strongly occupied, and the other is weakly occupied) and the case of the dissociating geminal (both orbitals in the geminal having the same occupation). We have:

$$\begin{aligned}&n_i \approx 1 \ \rightarrow \ \frac{1}{2}n_i{\mathcal{F}}_i^{{i}} \approx f_i^i , \end{aligned}$$

(31)

$$\begin{aligned}&n_i \approx 0 \ \rightarrow \ \frac{1}{2}n_i{\mathcal{F}}_i^{{i}} \approx -f_i^i + 2f_{i'}^{i'} -\big (v^{i'i'}_{i'i'}+v^{ii}_{i'i'}-2v^{ii'}_{ii'}\big ) , \end{aligned}$$

(32)

$$\begin{aligned}n_i & \approx \frac{1}{2} \ \rightarrow \ \frac{1}{2}n_i{\mathcal{F}}_i^{{i}} \approx f_{i'}^{i'}\nonumber \\&\quad +\frac{1}{2}\big (v_{ii}^{ii}-v^{ii}_{i'i'}-2v_{ii'}^{ii'}\big ) \approx f_{i}^{i} -\frac{1}{2}\big (v_{ii}^{ii}+v^{ii}_{i'i'}\big ) , \end{aligned}$$

(33)

$$\begin{aligned}&{\bar{n}}^a \approx 1 \ \rightarrow \ \frac{1}{2}{\bar{n}}^a{\mathcal{F}}^a_{{a}} \approx f_a^a , \end{aligned}$$

(34)

$$\begin{aligned}&{\bar{n}}^a \approx 0 \ \rightarrow \ \frac{1}{2}{\bar{n}}^a{\mathcal{F}}^a_{{a}} \approx -f_a^a + 2f_{a'}^{a'} +\big (v^{a'a'}_{a'a'}+v^{aa}_{a'a'}-2v^{aa'}_{aa'}\big ) , \end{aligned}$$

(35)

$$\begin{aligned}{\bar{n}}^a & \approx \frac{1}{2} \ \rightarrow \ \frac{1}{2}{\bar{n}}^a{\mathcal{F}}^a_{{a}} \approx f_{a'}^{a'}\nonumber \\&\quad -\frac{1}{2}\big (v_{aa}^{aa}-v^{aa}_{a'a'}-2v_{aa'}^{aa'}\big ) \approx f_{a}^{a}+\frac{1}{2}\big (v_{aa}^{aa}+v^{aa}_{a'a'}\big ) . \end{aligned}$$

(36)

We used the fact that in the dissociating regime $f_{i}^{i}\approx f_{i'}^{i'}$ and $v^{ii}_{ii}\approx v^{i'i'}_{i'i'}\approx v^{ii'}_{ii'}$ (and similar relations for a). Note that in the weakly interacting case considered in this paper, geminals are localized on subsystems. A dissociating geminal is assumed to be part of the subsystem, e.g. a molecule with stretched geometry. In this paper, the case $n_i\approx 1/2$ does not appear, all subsystems being closed-shell atoms.

All two-electron integrals appearing in (31)–(36) are positive. The two-electron contribution is small compared to the Fockian part around equilibrium geometry, while it is the other way around in the dissociating regime. Based on these equations and our initial assumption concerning the sign of $f_p^p$, one can argue the positivity of $\epsilon _i^a$. In, for example, (32), $2f_{i'}^{i'}$ is expected to give the (large negative) dominant contribution. The positivity of the dominant contribution in (35) can be argued similarly, leading to the positivity of $\epsilon ^a_i$, despite $\phi _i$ and $\phi _a$ being weakly and strongly occupied, respectively. As a further example, let us assume a dissociating geminal with ${\mathcal{I}}={\mathcal{A}}$, $i \ne a$ and $n_i \approx n_a \approx 1/2$. Then $i'=a$ and $a'=i$. Using (33) and (36), $\epsilon ^a_i$ reads:

$$\begin{aligned} \epsilon ^a_i \approx v^{ii}_{ii}+v^{aa}_{ii} \ge 0 . \end{aligned}$$