Ring coupled cluster doubles at the multireference level

Margócsy, Á.; Szabados, Á.

doi:10.1063/5.0005075

A ring approximation within an internally contracted multireference (MR) Coupled Cluster (CC) framework is worked out and tested. Derivation of equations utilizes MR based, generalized normal ordering and the corresponding generalized Wick-theorem (MR-GWT). Contractions among cluster operators are avoided by adopting a normal ordered exponential ansatz. The original version of the MR ring CC doubles (MR-rCCD) equations [Á. Szabados and Á. Margócsy, Mol. Phys. 115, 2731 (2017)] is rectified in two aspects. On the one hand, over-completeness of double excitations is treated by relying on the concept of frames. On the other hand, restriction on the maximal cumulant rank is lifted from two to four. This is found essential for obtaining reliable correlation corrections to the energy. The MR function underlying the approach is provided by the Generalized Valence Bond (GVB) model. The pair structure of the reference ensures a fragment structure of GVB cumulants. This represents a benefit when evaluating cumulant contractions appearing as a consequence of MR-GWT. In particular, cumulant involving terms remain less expensive than their traditional, pair-contracted counterpart, facilitating an $O (N^{6})$ eventual scaling of the proposed MR-rCCD method. Pilot applications are presented for covalent bond breaking, deprotonation energies, and torsional potentials.

I. INTRODUCTION

Coupled cluster (CC) theory is remarkably successful in capturing dynamical correlation effects when describing the electronic structure of molecular systems.¹ The singles and doubles CC scheme with perturbative triples, CCSD(T),² is an established “gold standard” in situations where the Hartree–Fock (HF) method represents an appropriate starting point. Research efforts invested in efficient implementation have recently allowed us to push the range of application of CCSD(T) to incredibly large systems.^3,4 When formulated properly, single reference based CC schemes are applicable also in situations where dynamical and static correlation are intertwined: methods based on the ansatz of Oliphant and Adamowicz,^5,6 Ivanov and Adamowicz,^7,8 or those harnessing the so-called CC moments⁹ are typical examples. Efficient implementation of these approaches has also been addressed.^10–13 The idea put forward by Rolik and Kállay¹⁴ represents an alternative single reference type technique, incorporating multireference (MR) effects via formulating the theory in terms of quasi-particles.

Adopting a single reference based methodology for an inherently MR problem is often accompanied by undesirable pivot dependence, rationalizing the development of genuine MR-CC techniques. Considerations based on the Bloch equation^15,16 have been traditionally pursued in this field, applying a formulation either in the Hilbert-space^17–20 or in the Fock-space.^21–23 Jezioski–Monkhorst (JM) parametrization¹⁷ of the wave operator has been prevalent in Hilbert-space approaches, involving an individual cluster operator assigned to each model space function. At difference with this, Fock-space theories assume a single cluster operator when composing the valence universal wave operator. The normal ordered exponential form for the wave operator

Ω = : \exp (T) :

(1)

proposed by Lindgren²⁴ has been prevalent in Fock-space theories, and application of Eq. (1) in the Hilbert-space context is relatively rare. In the above and throughout this work, colons indicate MR based, generalized normal order (MR-GNO).^25–27

An example for a Hilbert-space formulation exploiting Eq. (1) is the theory put forward by Mukherjee²⁵ and developed by Mukherjee and co-workers,²⁸ currently categorized as an internally contracted (ic), Hilbert-space MR-CC theory. The ic concept applied within MR-CC (icMR-CC) circumvents several difficulties associated with the JM ansatz.^29–31 The issue of sufficiency conditions, the unfavorable scaling of the theory with the number of determinants contributing to the reference as well as the “proper residual” problem,^32,33 is mitigated. In return, the wave operator of Eq. (1) introduces difficulties in MR generalized Wick-theorem (MR-GWT)^25,27,34,35 based matrix element evaluation, leading to a tremendous blow-up in the number and complexity of the contractions to derive and implement. At the time of its introduction, decay of the cumulant norm was anticipated with an increase in the cumulant rank, conceivably allowing one to put MR-GWT in practice.²⁷ It has been demonstrated later that high rank cumulants remaining comparable in magnitude to their low rank counterparts hinder the cumulant rank based cutoff, e.g., in covalent bond breaking situations.^36,37 A further problem associated with the wave operator of Eq. (1) is that its inverse is nontrivial to devise, impeding similarity transformation¹⁶ based manipulations.

Due to the complications associated with Eq. (1), icMR-CC formulations put in practice to date have mostly been employing the traditional

Ω = \exp (T)

(2)

form of the wave operator. The similarity transformed Hamiltonian is straightforward to introduce in this framework; the Baker–Campbell–Hausdorff (BCH) expansion, however, includes terms beyond four nested commutators when the cluster operator involves non-commuting excitations. Truncation of the BCH expansion is a common practice, reasoned with the anticipated smallness of cluster amplitudes, provided that the reference function is of sufficiently good quality. Matrix elements are typically evaluated based on the traditional Wick-theorem, with the closed shell core portion of the reference serving as Fermi-vacuum. The first icMR-CC implementations based on Eq. (2) were made by Evangelista and Gauss,³⁸ followed closely by Hanauer and Köhn.³⁹ A cluster operator in MR-GNO was later introduced in the argument of the exponential by Hanauer and Köhn in order to achieve size extensivity.⁴⁰ Among alternative routes explored within the domain of ordinary exponential parametrization, sequential transformation⁴¹ or the assumption of a unitary wave operator implying an anti-Hermitian cluster operator^42–44 may be mentioned.

One of the difficulties common to icMR-CC formulations irrespective of the parametrization of the wave operator is the redundancy among excited functions. Orthogonalization procedures involving redundancy filtering are commonly applied as a remedy.^28,42,45,46 Various ways of redundancy filtering have been extensively explored by Hanauer and Köhn.^39,40 Alternatively, Nooijen^47–49 advocated so-called many-body residuals instead of projections with redundancy filtered excited functions. The large number and complexity of the terms needing implementation are other adverse features of icMR-CC techniques, which are usually addressed by applying automated derivation tools, possibly invoking symbolic algebra.^{12,45,50–54}

Among the scarce use of the normal ordered exponential ansatz of Eq. (1) out of the scope of Fock-space theory, one may mention works elated by Nooijen^48,55 and Mukherjee.^56–59 Applications of MR-GWT in the context of Hilbert-space CC theories are more numerous, encompassing Canonical Transformation (CT) theory⁶⁰ and several studies by Nooijen and co-workers^45,49,55,61 or Sinha and co-workers.⁶² The challenge presented by MR-GWT may be met by truncation based on the cumulant rank^45,54 or based on the operator rank.^55,61 Original CT works with cumulants up to rank 2, while its quadratic extension involves cumulants of rank 3 in addition.^42,54 Besides neglecting cumulants beyond rank 2, terms quadratic in rank 2 cumulants are also discarded in state specific equation of motion theory.⁴⁵ Many-body residuals are in fact beneficial in this respect, cumulant involvement of the contributing terms being less complex.⁴⁹ Cumulant based truncation is occasionally avoided by transcribing the MR-GNO form of operators into normal order with respect to a genuine vacuum, following MR-GWT based manipulations, just for the purpose of matrix element evaluation by the ordinary Wick-theorem.^49,55,61

This work presents the first results on utilizing the ansatz of Eq. (1) and MR-GWT in the context of internally contracted, Hilbert-space MR-CC, as suggested by Mukherjee et al.^25,28 In particular, we report the ring approximation introduced by Čížek,⁶³ extended for this framework. The theory presented in Ref. 64 was largely motivated by the known correspondence between the random phase approximation (RPA) and ring CC doubles (rCCD) when based on the HF reference^65–67 together with the extended RPA (ERPA) put forward by Pernal.^68–70 The reason for opting for the icMR-CC framework of Mukherjee et al.^25,28 is that it facilitates an ic ring MR-CC doubles (icMR-rCCD) formulation paralleling ERPA at a large extent.

Our original icMR-rCCD theory⁶⁴ preceded numerical realization that has subsequently revealed two essential weaknesses, which are presently amended. One deficiency constitutes lack of redundancy handling when solving the projected residual equations. Presently, a treatment based on the concept of frames^71,72 is introduced to resolve this issue. Application of frames in MR based electronic structure methodology is rare though not unprecedented.⁷³ The second amendment affects the cutoff in MR-GWT based evaluation of matrix elements. Contrary to previous experience, we find unacceptably large errors in correlation energies when neglecting cumulants beyond rank two. This may be reasoned with the lack of commutators in the present theory. The cumulant rank based cutoff is consequently lifted from rank 2 to rank 4, implying python assisted derivation and implementation of the terms additional to those reported in Ref. 64.

The ring approximation applied in this work represents one of the means of keeping at bay the incredibly large number of terms generated by MR-GWT. The pair structure of the reference wavefunction is a second ingredient that contributes to achieving an icMR-rCCD method applicable in practice. The antisymmetrized product of strongly orthogonal geminals (APSG) function^74,75 assumed as reference lends a beneficial fragment structure to cumulants,^76,77 which is further simplified when assigning not more than two orbitals to an electron pair, resulting in the function known as Generalized Valence Bond (GVB).⁷⁸ The fact that nonzero cumulants can appear solely with all indices assigned to the same electron pair ensures that cumulant involving contractions are less expensive than their pair-contracted counterpart. As a consequence, the $O (N^{6})$ scaling characteristic of HF based rCCD remains valid for the icMR-rCCD worked out presently.

A common trait of the here proposed icMR-rCCD method and previously reported geminal based correlation schemes is the inherent utilization of the two-electron fragment structure of the reference. Excited functions contributing to the correction typically inherit a geminal structure characteristic of the ground state.^73,79–87 In the ic framework proposed here, the geminal structure of the reference naturally gets harnessed at the level of cumulants.

This paper is organized as follows: Basic notations introduced in Sec. II A are followed by the derivation of the icMR-rCCD equations in Sec. II B, with some of the formulas deferred to Appendix A for transparency. Details on the cumulant structure and evaluation are presented in Sec. II C and Appendix B. Parametrization of the cluster operator is developed in Sec. II D, and redundancy treatment is briefly outlined in Sec. II E with more details provided in Appendix C. Formal considerations are closed by some remarks on size consistency and extensivity in Sec. II F. Pilot numerical applications, given in Sec. III, help assess icMR-rCCD compared to ERPA and to an appropriate benchmark. The role of the cumulant based cutoff as well as various parametrizations of the cluster operator is also illustrated numerically.

Since we shall be concerned solely with internally contracted theory, the “ic” designation is omitted from acronyms further on.

II. THEORY

A. Notations

The tensor notation of Harris et al.⁸⁸ is adopted, $a^{p σ} = a_{p σ}^{†}$ standing for the creation operator of spin-orbital pσ that is assumed to be the product of spatial orbital p and spin function σ. In general, upper and lower indexed quantities are related to each other by Hermitian conjugation. Spin summation is implied in excitation operators of rank K, denoted as

E_{q_{1} \dots q_{K}}^{p_{1} \dots p_{K}} = \sum_{σ_{1}, \dots, σ_{K}} a_{q_{1} σ_{1} \dots q_{K} σ_{K}}^{p_{1} σ_{1} \dots p_{K} σ_{K}},

ensuring commutation with total spin, $[S^{2}, E_{q_{1} \dots q_{K}}^{p_{1} \dots p_{K}}] = 0$ ⁠. Throughout this paper, we apply MR-GNO^25,27,35 associated with arbitrary reference function Φ. The normal product form of the Hamiltonian, H_N, reads as

H_{N} = H - ⟨) Φ | H | Φ (⟩ .

Based on MR-GWT,^25,27,35 the normal product form, H_N, involves one-body and two-body terms expressed with MR-GNO fermion products,

H_{N} = : H_{1} : + : H_{2} : .

(3)

Terms of H_N can be written with the aid of spin-free excitation operators as

: H_{1} : = \sum_{p, q} f_{p}^{q} : E_{q}^{p} :,

(4a)

: H_{2} : = \frac{1}{2} \sum_{p, q, r, s} v_{p r}^{q s} : E_{q s}^{p r} : .

(4b)

In the above equations, colons (:) indicate MR-GNO, $v_{p r}^{q s} = \int φ_{p}^{*} (1) φ_{r}^{*} (2) r_{12}^{- 1} φ_{q} (1) φ_{s} (2) d r_{1} d r_{2}$ is the two-electron integral, and $f_{p}^{q}$ stands for the element of the generalized Fockian, given by

f_{p}^{q} = h_{p}^{q} + \sum_{r} {\bar{v}}_{p r}^{q r} n_{r},

(5)

where $h_{p}^{q}$ is the one-electron integral including kinetic energy and nuclear-electron repulsion terms, and the antisymmetrized two-electron integral is defined as

{\bar{v}}_{p r}^{q s} = 2 v_{p r}^{q s} - v_{r p}^{q s} .

The expression of the Fockian in Eq. (5) implies natural orbitals of Φ; moreover, the spin density is assumed to be zero, allowing us to write

⟨) Φ | a_{q α}^{p α} | Φ (⟩ = ⟨) Φ | a_{q β}^{p β} | Φ (⟩ = δ_{p q} n_{p} .

The hole occupation number is introduced as

{\bar{n}}_{p} = 1 - n_{p},

with occupation number n_p and hole occupation number ${\bar{n}}_{p}$ both falling in the interval [0, 1].

We shall primarily be concerned with an APSG reference function. Strong orthogonality being equivalent to an expansion of geminals in mutually disjoint subsets of orthonormal orbitals,⁸⁹ the APSG wavefunction takes the form

| Φ^{APSG} (⟩ = \prod_{P = 1}^{N / 2} \sum_{p \in P} c_{p} a^{p α} a^{p β} |vac⟩,

(6)

where P is the geminal index, c_p is the geminal coefficient, $|vac⟩$ denotes the bare vacuum, and N stands for the even number of electrons in the system. Note that Eq. (6) involves singlet geminals with s_z = 0 (so-called perfect pairing) and orbitals that are natural orbitals of Φ^APSG, c_p being related to n_p as |c_p|² = n_p. The wavefunction widely known as GVB⁷⁸ represents a special case of APSG, geminals accommodating two orbitals at most. While some of the results below concern an arbitrary reference, others resort to an APSG construction. Track shall be kept by using notation Φ or Φ^APSG, respectively.

B. Ring coupled cluster doubles equation

We start by presenting the internally connected MR-CCD equations obtained from the work of Mukherjee et al.^25,28 The procedure parallels the textbook derivation of HF based CCD,¹⁶ with the differences that (i) the expectation values are evaluated using MR-GWT^27,35 and (ii) a normal ordered exponential wave operator [Eq. (1)] is used instead of the exponential of MR-GNO operators. The latter step is taken in order to avoid contractions between terms of the cluster operator. Introduction of the similarity transformed Hamiltonian is avoided for the complications associated with the inverse of the wave operator of Eq. (1).

The normal ordered exponential ansatz, put forward by Lindgren,²⁴ is written as

| Ψ^{CCD} (⟩ = : \exp (T) : | Φ (⟩

(7)

for an arbitrary reference function Φ. The cluster operator is expressed with spin-free double excitations as

T = \frac{1}{2} \sum_{\begin{array}{c} i, j \\ a, b \end{array}} t_{a b}^{i j} E_{i j}^{a b},

(8)

indices i, j, a, b referring to spatial orbitals. While no particle/hole character is associated with the orbitals, a distinction between indices i, j and a, b is still implied. Letters i, j, … are used for orbitals that are bound to contribute an occupation number in the leading (i.e., pair-contracted) term when evaluating matrix elements of the CCD equations by MR-GWT. Similarly, orbitals that are bound to contribute a hole occupation number in the leading term are denoted by a, b, …. Letters p, q stand for arbitrary indices in this respect. At this point, no restriction is imposed on indices i, j, a, b in Eq. (8). This shall be amended in Sec. II D, in knowledge of the eventual structure of the ring CCD equations. The cluster amplitude $t_{a b}^{i j}$ is not anti-symmetric with respect to permutation of either its upper or lower indices. At the same time, $t_{a b}^{i j}$ is invariant to permuting both its upper and lower indices, i.e., $t_{a b}^{i j} = t_{b a}^{j i}$ ⁠.

Substituting Eq. (7) into the

H_{N} | Ψ (⟩ = Δ E | Ψ (⟩

(9)

Schrödinger equation with ΔE = E −⟨Φ|H|Φ⟩ and E denoting the eigenvalue corresponding to Ψ, projection of Eq. (9) with ⟨Φ| yields the MR-CCD energy correction

\begin{matrix} Δ E & = & ⟨ H_{N} : \exp (T) : ⟩ = \frac{1}{2} \sum_{\begin{array}{c} i, j \\ a, b \end{array}} ⟨ H_{N} : E_{i j}^{a b} : ⟩ t_{a b}^{i j} + O (t^{2}) \\ = & \sum_{\begin{array}{c} i, j \\ a, b \end{array}} n_{i} n_{j} {\bar{n}}^{a} {\bar{n}}^{b} B_{i j}^{a b} t_{a b}^{i j} + O (t^{2}), \end{matrix}

(10)

making use of ⟨Φ|: exp(T): Φ⟩ = 1, a consequence of the definition of MR-GNO. Brackets here and further on stand for the expectation value taken with Φ. In the last line of Eq. (10), we have introduced

n_{i} n_{j} {\bar{n}}^{a} {\bar{n}}^{b} B_{i j}^{a b} = \frac{1}{2} ⟨ H_{N} : E_{i j}^{a b} : ⟩,

(11)

many-body expression of $B_{i j}^{a b}$ being given in Appendix A. The leading term of $B_{i j}^{a b}$ is ${\bar{v}}_{i j}^{a b}$ ⁠, followed by cumulant involving terms that give zero upon Φ becoming the HF determinant. Indices of particle/hole occupation numbers appearing in the subscript or superscript follow indexing of the associated tensor element, e.g., $B_{i j}^{a b}$ in Eq. (11). It has no meaning beyond that, in particular, n_i = nⁱ and ${\bar{n}}_{a} = {\bar{n}}^{a}$ ⁠.

The CCD amplitude equation is obtained by projecting the Schrödinger equation with $⟨) Φ | : E_{a b}^{i j} :$ ⁠, leading to

⟨ : E_{a b}^{i j} : H_{N} : \exp (T) : ⟩ = Δ E ⟨ : E_{a b}^{i j} : : \exp (T) : ⟩ .

(12)

Expanding both sides and resorting to terms of ΔE linear in t, the amplitude equation takes the form

\begin{matrix} 2 n^{i} n^{j} {\bar{n}}_{a} {\bar{n}}_{b} B_{a b}^{i j} + ⟨ : E_{a b}^{i j} : H_{N} : T : ⟩ + \frac{1}{2} {⟨ : E_{a b}^{i j} : H_{N} : T^{2} : ⟩}_{C} + O (t^{3}) = 0, \end{matrix}

(13)

where the subscript C refers to the connected part of the expectation value (the disconnected part canceling with the linear term of the energy correction).

The MR-CCD equations of the present formalism are given by Eq. (13) in the approximation where ΔE is linear in T. Main differences between Eq. (13) and its HF based counterpart stem from MR-GNO and the associated MR-GWT. A first observation is that matrix element evaluation by MR-GWT generates a plethora of cumulant involving terms beyond the leading, pair-contracted expressions. In Eq. (13), cumulants up to rank 6 appear in the matrix element linear in T, and the quadratic term involves cumulants up to rank 8 when evaluated without any neglect. A further remark concerns the expansion of the left hand side of Eq. (13) in powers of T. This is a non-terminating series due to the appearance of cumulants in MR-GWT and due to the fact that cumulants do not vanish with their rank increased.

The terms of Eq. (13) are consequently evaluated by introducing approximations. It is already apparent at this stage that cluster amplitudes are assumed small enough to allow for neglection based on powers of t. A further step in the direction of breaking down the number of terms is to evoke the ring approximation introduced by Čížek.⁶³ The ring approximation generalized for the MR case in Ref. 64 resorted to terms linear in cumulants and involved a truncation beyond cumulant rank Λ₂. Both of these restrictions are abandoned presently, keeping all terms conforming with the Riccati type amplitude equation

B^{†} + 2 A^{T} t + 2 t A + 4 t M B M t = 0

(14)

characteristic of ring CCD. Tensor M in the above equation includes the product of occupation numbers along its diagonal,

M_{i j}^{a b} = δ_{i j} δ_{a b} n_{i} {\bar{n}}^{a} .

Elements of tensors A and B are given in Appendix A. Inspection of Appendix A reveals that the ring approximation retains just a handful of terms in Eq. (13). The rank of cumulants involved in Eq. (13) is also effectively reduced, the highest cumulant rank appearing in Eq. (14) being Λ₄.

Generalization of the ring approximation as described in Ref. 64 is not invariant to rotation among orbitals of degenerate occupation number at the level of APSG, i.e., within the core (n_i = 1) and virtual (⁠ ${\bar{n}}_{a} = 1$ ⁠) subset. While this is a common feature with HF based ring CCD, orbital invariance is restored in the present formulation by allowing the full tensor F to enter tensor A, cf. the first two terms of Eq. (A3) in Appendix A. Rotation among orbitals of fractional occupation number is allowed only within geminal subspaces indexed by P in Eq. (6). Even this latter freedom is eliminated in the present MR-rCCD by requiring the APSG 1-RDM to be diagonal.

C. Cumulants

Cumulants λ_K are introduced via their relation to reduced density matrices (RDMs), γ_K.⁷⁶ The corresponding spin-orbital expressions, up to rank K = 4, necessary in the present approach, read

γ_{1} = λ_{1},

(15a)

γ_{2} = λ_{2} + \frac{1}{2} λ_{1} \otimes λ_{1},

(15b)

γ_{3} = λ_{3} + λ_{2} \otimes λ_{1} + \frac{1}{6} λ_{1} \otimes λ_{1} \otimes λ_{1},

(15c)

\begin{matrix} γ_{4} = λ_{4} + λ_{3} \otimes λ_{1} + \frac{1}{2} λ_{2} \otimes λ_{2} + \frac{1}{2} λ_{2} \otimes λ_{1} \otimes λ_{1} + \frac{1}{24} λ_{1} \otimes λ_{1} \otimes λ_{1} \otimes λ_{1} . \end{matrix}

(15d)

The notation of Ref. 37 is adopted in the above equation, with ⊗ standing for the antisymmetrized tensor product with index permutations performed only between terms of the product and never within one term. Factors 1/m! compensate m!-tuple counting where m equivalent terms are present in the product.

The structure of the APSG wavefunction together with cumulant decomposition of the RDMs [Eq. (15)] ensures that cumulants of rank 2 and above exhibit what is termed here “geminal connectedness.” Under “geminal connectedness,” it is understood that nonzero elements occur for all orbital indices of λ belonging to the same geminal.⁷⁶ This is a consequence of the APSG wavefunction being an antisymmetrized direct product of complete active space functions. The fact that these active spaces are furnished with two electrons brings further simplification: RDM elements of rank γ₃ and above are zero when indexed by orbitals of the same geminal. As a consequence, $γ_{K}^{APSG}$ is not contributing to the nonzero elements of the $λ_{K}^{APSG}$ for K = 3, … Focusing on intrageminal elements, the APSG cumulants satisfy the simplified relations⁹⁰

λ_{3}^{APSG} = - λ_{2} \otimes λ_{1} - \frac{1}{6} λ_{1} \otimes λ_{1} \otimes λ_{1},

(16a)

\begin{align} λ_{4}^{APSG} = & - λ_{3} \otimes λ_{1} - \frac{1}{2} λ_{2} \otimes λ_{2} - \frac{1}{2} λ_{2} \otimes λ_{1} \otimes λ_{1} \\ - \frac{1}{24} λ_{1} \otimes λ_{1} \otimes λ_{1} \otimes λ_{1} \\ = & - \frac{1}{2} λ_{2} \otimes λ_{2} + \frac{1}{2} λ_{2} \otimes λ_{1} \otimes λ_{1} + \frac{1}{8} λ_{1} \otimes λ_{1} \otimes λ_{1} \otimes λ_{1} . \end{align}

(16b)

While being geminal connected, cumulants of Eqs. (16a) and (16b) are manifestly disconnected. Their role of canceling the disconnected part of γ_K is lost due to the vanishing intrageminal elements of the corresponding RDMs. The situation is even more peculiar for a two-electron system, since cumulants of Eqs. (16a) and (16b) do not vanish, unlike the RDMs of rank 3 and higher. (The ring approximation presently limits highest cumulant at Λ₄, and it therefore can only occur for two-electron systems that cumulant rank, K, exceeds the total number of electrons, N.) Although a disconnected cumulant contradicts the usual physical interpretation, Λ_K up to K = 4 are essential for the fulfillment of MR-GWT used to evaluate ΔE, and cumulants up to rank four are therefore retained. An illustration on the role of the cumulant rank in the energy correction, provided in Sec. III A 1, supports the claim that these cumulants are indispensable.

Cumulant decomposition of RDMs and application of the MR-GWT are in principle straightforward for spin-orbitals; it is, however, highly nontrivial to express the results in terms of spin-summed cumulants, Λ_K, defined as

Λ_{p_{1}^{'} \dots p_{K}^{'}}^{p_{1} \dots p_{K}} = \sum_{σ_{1}, \dots, σ_{K}} λ_{p_{1}^{'} σ_{1} \dots p_{K}^{'} σ_{K}}^{p_{1} σ_{1} \dots p_{K} σ_{K}} .

(17)

Difficulty in introducing Λ_K is associated with the appearance of products of cumulants indexed by spin-orbitals and summed for spin. Expressing the products in terms of individual, spin-summed cumulants requires the inverse of relation (17), i.e., elements of λ_K expressed in terms of Λ_K. This can be obtained with the aid of additional spin symmetry relations exhibited by λ_K; the task, however, gradually becomes more complicated with the increase in the cumulant rank.^34,91 Derivation of the working formulas reported here has been carried out in terms of spin-orbitals, with spin summation performed at the final stage, e.g., by applying Eqs. (B3) and (B4).

Elements of the APSG cumulants Λ_K up to rank K = 3 are given in Appendix B. A glance at the expressions in Appendix B reveals that APSG cumulants are highly factorized. Each nonzero term of Λ₂ is proportional to a product of Kronecker deltas, yielding nonzero elements with at most 2 nonidentical orbital indices. This property carries through to Λ₃ and Λ₄, their nonzero elements involving at most 3 and 4 nonidentical orbital indices, respectively. The reason behind vanishing cumulant elements can be quite general. For example, cumulants of odd rank are zero for wavefunctions exhibiting particle-hole symmetry.³⁷ The Pauli principle has the consequence that three spatial indices repeated in the subscript or superscript results in a zero cumulant element

Λ_{\dots}^{p p p \dots} = Λ_{p p p \dots}^{\dots} = 0 .

Cumulants Λ_K consequently vanish with APSG for rank 2d_max < K, where d_max is the size of the biggest geminal subspace.⁷⁶ Considering GVB where d_max = 2, cumulants Λ_K vanish for 4 < K.⁹² This provides further justification for the neglect of terms involving cumulants Λ_K, 4 < K ≤ N when evaluating Eq. (13), for the case of GVB reference. An upper limit N2^2K−1 for the nonzero elements of GVB cumulants can also be deduced from the above for K ≤ 4, giving not more than, e.g., 8N nonzero elements for Λ₄.

D. Pruning the excitation manifold

In lack of particle-hole categorization of orbitals, no restriction on indices entering Eq. (8) has been assumed so far. This picture needs refinement for the following reason. Consider the first approximation of the amplitudes, $t_{a b}^{i j (1)}$ obtained at the first Jacobi iteration step when solving Eq. (14),

t_{a b}^{i j (1)} = - \frac{1}{2} {(B_{i j}^{a b})}^{*} {(δ_{a}^{c} δ_{k}^{i} A_{k a}^{c i} + δ_{b}^{c} δ_{k}^{j} A_{k b}^{c j})}^{- 1},

(18)

and inspect the case where the reference Φ becomes HF and orbital b becomes occupied, i.e., ${\bar{n}}_{b} \to 0$ ⁠. The contribution of amplitudes involving excitation to orbital b is expected to be small when ${\bar{n}}_{b}$ is close to 0, and the contribution should vanish in the limit. This is ensured by the factor ${\bar{n}}^{b}$ in Eq. (10), provided that the amplitude is finite. Unfortunately, Eq. (18) does not conform to this expectation since terms of the numerator may tend to infinity. Examine, e.g., the third term in the fourth line of Eq. (A1). Performing complex conjugation and substituting the expression of Λ₂ from Eq. (B1), one obtains

\begin{matrix} - \frac{1}{2 n^{i} {\bar{n}}_{b}} \sum_{p, q} v_{a q}^{p j} Λ_{p b}^{i q} & = & - δ_{b}^{i} \sum_{p \in B} \frac{n_{p}}{{\bar{n}}_{b}} v_{a p}^{p j} + 2 δ_{B}^{I} \frac{n_{b}}{{\bar{n}}_{b}} v_{a b}^{i j} - δ_{B}^{I} \frac{c^{i} c_{b}}{n^{i} {\bar{n}}_{b}} v_{a i}^{b j} \\ = & - δ_{b}^{i} (v_{a b}^{b j} + \sum_{\begin{array}{c} p \neq b \\ p \in B \end{array}} \frac{n_{p}}{{\bar{n}}_{b}} v_{a p}^{p j}) + (1 - δ_{b}^{i}) δ_{B}^{I} \\ \times (2 \frac{n_{b}}{{\bar{n}}_{b}} v_{a b}^{i j} - \frac{c^{i} c_{b}}{n^{i} {\bar{n}}_{b}} v_{a i}^{b j}) . \end{matrix}

(19)

Analyzing Eq. (19), it is apparent that the nature of the outer index pair ib occurring on Λ is decisive on the magnitude of the term. When i and b belong to different geminals, the term is zero due to $Λ_{i q}^{p b}$ being zero. When indices i and b match, the first term on the right hand side of Eq. (19) is, in effect, finite despite ${\bar{n}}_{b}$ figuring in the denominator. To see this, one has to take into account that p resorts to a single term for a GVB reference and $n_{p} = {\bar{n}}_{b}$ as a consequence of n_p + n_b = 1. The quotient $n_{p} / {\bar{n}}_{b}$ therefore simplifies to 1, avoiding divergence as ${\bar{n}}_{b} \to 0$ ⁠. In the APSG case, $n_{p} / {\bar{n}}_{b} \leq 1$ holds for all p ∈ B, p ≠ b, again ensuring that divergence is avoided. When i and b belong to the same geminal and i ≠ b, the second term on the right hand side of Eq. (19) is to be examined and found divergent for ${\bar{n}}_{b} \to 0$ ⁠. Excitations entering Eq. (8) are pruned in order to avoid such divergence. Pursuing the analysis, analogous terms of ${(B_{i j}^{a b})}^{*}$ can be identified that behave similarly and call for pruning for index pairs ja, ia, and jb.

Regarding the index quartet i, j, a, b, there are two more index pairings left, ab and ij. They can both cause trouble, harmful terms being the first and the second of Eq. (A1) in line three, reading as

\begin{matrix} \frac{1}{2 {\bar{n}}_{a} {\bar{n}}_{b}} \sum_{p, q} v_{p q}^{i j} Λ_{a b}^{p q} & = & δ_{a b} (\frac{n_{a}}{{\bar{n}}_{a}} v_{a a}^{i j} + \frac{c^{a}}{{\bar{n}}_{a}^{2}} \sum_{\begin{array}{c} p \neq a \\ p \in A \end{array}} v_{p p}^{i j} c^{p}) \\ - (1 - δ_{a b}) δ_{A B} \frac{n_{a}}{{\bar{n}}_{b}} \frac{n_{b}}{{\bar{n}}_{a}} {\bar{v}}_{a b}^{i j} \end{matrix}

(20)

and

\frac{1}{2 n^{i} n^{j}} \sum_{p, q} v_{a b}^{p q} Λ_{p q}^{i j} = δ_{i j} (\frac{{\bar{n}}^{i}}{n^{i}} v_{a b}^{i i} + \frac{c^{i}}{{(n^{i})}^{2}} \sum_{\begin{array}{c} p \neq i \\ p \in I \end{array}} v_{a b}^{p p} c_{p}) - (1 - δ_{i j}) δ_{I J} {\bar{v}}_{a b}^{i j},

(21)

respectively. Similar to the case of Eq. (19), indices of Λ occurring on the same geminal are problematic, since the cumulant is zero otherwise. At difference with Eq. (19), it is index matching that causes divergence of the expressions. The first term on the right hand side of Eq. (20) diverges for δ_ab in the limit ${\bar{n}}_{a} \to 0$ ⁠. The first term on the right hand side of Eq. (21) similarly tends to infinity for δ_ij in the limit nⁱ → 0. The i ≠ j term of Eq. (21) is apparently insensitive to the value of the occupation numbers. As mentioned before, $n_{a} / {\bar{n}}_{b} \leq 1$ as well as $n_{b} / {\bar{n}}_{a} \leq 1$ for a ≠ b; consequently, the second term on the right hand side of Eq. (20) is also well behaving for ${\bar{n}}_{a} \to 0$ or ${\bar{n}}_{b} \to 0$ ⁠. Divergence is again avoided by introducing pruning in Eq. (8) for i = j and a = b.

In order to complete the picture, it remains to be seen that the denominator of Eq. (18) is finite when the numerator diverges. This can be checked by finding the dependence of the diagonal terms of tensor A on occupation numbers and geminal coefficients. Explicit evaluation has been performed for the contributions of tensor F, i.e., the first two terms on the right hand side of Eq. (A3). These involve

\begin{matrix} \frac{1}{2} δ_{i k} n_{i} F_{k}^{i} ({\bar{n}}^{i} \to - n^{i}) \\ = (2 n_{i} - 1) (f_{i}^{i} - \sum_{\begin{array}{c} p \in I \\ p \neq i \end{array}} {\bar{v}}_{i p}^{i p} n_{p}) + 2 \sum_{\begin{array}{c} p \in I \\ p \neq i \end{array}} f_{p}^{p} n_{p} + 2 n_{i} {\bar{n}}_{i} v_{i i}^{i i} \\ + \sum_{\begin{array}{c} p \in I \\ p \neq i \end{array}} (v_{p p}^{i i} c_{i} c^{p} + v_{i i}^{p p} c^{i} c_{p} - 2 {\bar{v}}_{i p}^{i p} n_{i} n_{p}) \\ - \sum_{\begin{array}{c} p, q \in I \\ p, q \neq i \end{array}} (2 {\bar{v}}_{p q}^{p q} n_{p} n_{q} - v_{q q}^{p p} c_{p} c^{q}) \end{matrix}

(22)

and

\begin{matrix} \frac{1}{2} δ_{a c} {\bar{n}}_{a} F_{a}^{c} (n_{a} \to - {\bar{n}}_{a}) & = & - (2 n_{a} - 1) (f_{a}^{a} + n_{a} \sum_{\begin{array}{c} p \in A \\ p \neq a \end{array}} {\bar{v}}_{a p}^{a p} \frac{n_{p}}{{\bar{n}}_{a}}) \\ + 2 n_{a} \sum_{\begin{array}{c} p \in A \\ p \neq a \end{array}} f_{p}^{p} \frac{n_{p}}{{\bar{n}}_{a}} - 2 n_{a} {\bar{n}}_{a} v_{a a}^{a a} \\ - \sum_{\begin{array}{c} p \in A \\ p \neq a \end{array}} (v_{p p}^{a a} c_{a} c^{p} + v_{a a}^{p p} c^{a} c_{p} - 2 {\bar{v}}_{a p}^{a p} n_{a} n_{p}) \\ - n_{a} \sum_{\begin{array}{c} p, q \in A \\ p, q \neq a \end{array}} (2 {\bar{v}}_{p q}^{p q} \frac{n_{p} n_{q}}{{\bar{n}}_{a}} - v_{q q}^{p p} \frac{c_{p} c^{q}}{{\bar{n}}_{a}}) \end{matrix}

(23)

with the notation of Appendix A. Apparently, Eq. (22) shows no ill-effect as n_i → 0. Although ${\bar{n}}_{a}$ appears in the denominator of Eq. (23), the expression involves $n_{p} / {\bar{n}}_{a}$ and $\sqrt{n_{p} / {\bar{n}}_{a}}$ for a ≠ p, both being smaller than one for ${\bar{n}}_{a} \to 0$ ⁠. Examination of the contribution of the third term on the right hand side of Eq. (A3) to Eq. (18) similarly reveals no divergence with nⁱ, n^j, ${\bar{n}}_{a}$ ⁠, or ${\bar{n}}_{b}$ tending to zero.

In view of the above, restriction on indices i, j and a, b entering Eq. (8) is set as

$n^{i} {\bar{n}}_{a} > 1 / 4$ for i ≠ a and I = A;
$n^{i} {\bar{n}}_{b} > 1 / 4$ for i ≠ b and I = B;
$n^{j} {\bar{n}}_{a} > 1 / 4$ for j ≠ a and J = A;
$n^{j} {\bar{n}}_{b} > 1 / 4$ for j ≠ b and J = B;
nⁱ > 1/2 for i = j; and
${\bar{n}}_{a} > 1 / 2$ for a = b.

Pruning is an obvious disadvantage since a change in the set of amplitudes allowed in Eq. (8) may lead to discontinuity on the potential energy surface (PES). The above thresholds have been set to avoid such problems for covalent bond dissociation described by GVB. Focusing on a dissociating geminal, the process is characterized by one occupation number evolving from 1 − η to 1/2 + η and the other changing from η to 1/2 − η, the two adding up to one and η being some small positive number. No combination of the particle and hole occupation numbers of these two orbitals can cross the thresholds in (i)–(iv) during dissociation. A new situation (e.g., APSG reference instead of GVB) may call for revision of the restrictions affecting Eq. (8), the principal aim being avoidance of discontinuities. Ultimately, it would be beneficial to eliminate the need of pruning. It, however, appears unavoidable at the present stage of the theory.

It may be interesting to observe that restrictions (i)–(vi) allow for spectator excitations. [The case i = j = a = b is excluded by points (v) and (vi).] It is also admissible that n_i < n_a, provided that I ≠ A. This case deserves a note since the excitation operator of the companion ERPA correction^64,69,70 involves such deexcitation-type transitions, parametrized by the amplitude matrix Y. A numerical example on the role of spectator and deexcitation type transitions is given in Sec. III A 1.

E. Redundancy handling

Internally contracted excitations bring along over-completeness of the set of excited functions, $: E_{i j}^{a b} : | Φ (⟩$ ⁠. Redundancy is usually handled by elimination of certain amplitudes or linear combinations thereof. An unfortunate consequence of amplitude elimination is the possible appearance of discontinuities on the PES, as alluded to in Sec. II D. Interference with orbital invariance and size consistency has been reported for some of the redundancy treatments investigated in Ref. 39, with size extensivity restored in a follow-up work.⁴⁰ To avoid such pitfalls, Nooijen and co-workers advocate the use of many-body residual equations.

The present work applies a hitherto unexplored scheme for handling redundancy when CC amplitude equations are obtained as ordinary projections. We rely on the concept of frames, which facilitates keeping all the amplitudes corresponding to the redundant set of excited functions. Frames are generalizations of basis sets, relaxing the requirement of linear independence on the constituting vectors. Frames were originally introduced in connection with signal analysis and have been finding ever increasing applications in physics, engineering, and computer science.^71,72 In this section, we give a recipe-style presentation of the procedure followed when solving Eq. (12). For rationalization, we refer to Appendix C.

Let us introduce shorthand

| ψ_{μ} (⟩ = : E_{i j}^{a b} : | Φ (⟩

(24)

for bookkeeping. Whenever the overlap matrix S, composed of the elements

S_{μ ν} = ⟨) ψ_{μ} | ψ_{ν} (⟩,

possesses a zero eigenvalue, amplitude equations [Eq. (12)] contain less information than the number of parameters. Among such circumstances, the amplitude vector t, collecting elements t_μ, is ill-defined. With the aim of defining amplitudes in a unique manner, write the spectral decomposition of S as

S = V Σ V^{†},

(25)

where the matrix Σ contains Σ_μ along its diagonal, among which $μ = 1, \dots, M$ are nonzero and Σ_μ = 0 for $μ = M + 1, \dots, N$ ⁠. For further convenience, we denote as Σ′ = diag(Σ_μ) the $M \times M$ matrix collecting nonzero Σ_μ and $σ^{'} = \sqrt{Σ^{'}}$ ⁠.

As detailed in Appendix C, the amplitude update boils down to determining t^e corresponding to the auxiliary Löwdin-basis of Eq. (C1). The equation determining t^e reads

α t^{e} = - β,

(26)

where the matrix α and vector β are given by Eqs. (C19) and (C21), respectively. Once t^e is determined, the amplitude vector corresponding to the set of over-complete functions ${ψ_{μ}}_{μ = 1}^{N}$ is generated as

t = V (\begin{matrix} σ^{' - 1} \\ 0 \end{matrix}) t^{e} .

(27)

With the overlap matrix exhibiting a block-diagonal structure, the procedure can be performed separately for each block of S of dimension larger than one.

In order to be more specific, we briefly work out the update for the ring CC equation of Eq. (12), using the example of a two-dimensional block of S, formed of exchange spectator excitations $: E_{i 2}^{2 a} : | Φ (⟩ = | ψ_{1} (⟩$ and $: E_{i 3}^{3 a} : | Φ (⟩ = | ψ_{2} (⟩$ ⁠. We assume that i is a core orbital (n_i = 1), a is a virtual orbital (n_a = 0), and orbitals 2, 3 belong to a two-dimensional geminal (n₂ + n₃ = 1). Form Eq. (12) of the amplitude equations allows us to directly identify constituents of the linear system of equations obtained with ψ, denoted α^ψ and β^ψ and read as

α^{ψ} (\begin{matrix} 2 A_{22}^{i i} + 2 A_{22}^{a a} & 0 \\ 0 & 2 A_{33}^{i i} + 2 A_{33}^{a a} \end{matrix})

and

β^{ψ} (\begin{matrix} {(B_{2 a}^{i 2})}^{*} + 2 {(A^{T} t)}_{2 a}^{i 2} + 2 {(t A)}_{2 a}^{i 2} + 4 {(t M B M t)}_{2 a}^{i 2} - (2 A_{22}^{i i} + 2 A_{22}^{a a}) t_{2 a}^{i 2} \\ {(B_{3 a}^{i 3})}^{*} + 2 {(A^{T} t)}_{3 a}^{i 3} + 2 {(t A)}_{3 a}^{i 3} + 4 {(t M B M t)}_{3 a}^{i 3} - (2 A_{33}^{i i} + 2 A_{33}^{a a}) t_{3 a}^{i 3} () \end{matrix}) .

Constituents of Eq. (26) can be obtained based on Eq. (C1) by the transformations

α = (σ^{' - 1} 0) V^{†} α^{ψ} V (\begin{matrix} σ^{' - 1} \\ 0 \end{matrix})

and

β = (σ^{' - 1} 0) V^{†} β^{ψ} .

Working out the 2 × 2 block of the overlap matrix of this example, one arrives at

S = 2 n_{2} n_{3} (\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}),

giving rise to

σ^{'} = 2 \sqrt{n_{2} n_{3}}

and

V = \frac{1}{\sqrt{2}} (\begin{matrix} 1 & 1 \\ - 1 & 1 \end{matrix}) .

The matrix α and vector β are consequently built with a single component, which reads as

α = \frac{1}{4 n_{2} n_{3}} (A_{22}^{i i} + A_{22}^{a a} + A_{33}^{i i} + A_{33}^{a a})

and

β = \frac{1}{2 \sqrt{2 n_{2} n_{3}}} (β_{1}^{ψ} - β_{2}^{ψ}) .

The single amplitude on the Löwdin-basis is updated as

t^{e} = \sqrt{2 n_{2} n_{3}} \frac{β_{2}^{ψ} - β_{1}^{ψ}}{A_{22}^{i i} + A_{22}^{a a} + A_{33}^{i i} + A_{33}^{a a}},

yielding the update for the 2-component amplitude vector

(\begin{matrix} t_{2 a}^{i 2} \\ t_{3 a}^{i 3} \end{matrix}) = \frac{β_{2}^{ψ} - β_{1}^{ψ}}{2 (A_{22}^{i i} + A_{22}^{a a} + A_{33}^{i i} + A_{33}^{a a})} (\begin{matrix} 1 \\ - 1 \end{matrix}) .

(28)

The above 2 × 2 example is special in the sense that the linear term of Eq. (12) involves no coupling between the two selected excitations. Whenever such a coupling is present, the corresponding term is to be retained in α^ψ, which consequently becomes non-diagonal.

The dimension of the redundant blocks obviously depends on the pruning affecting Eq. (8). With the pruning strategy described in Sec. II D, the redundant blocks do not get larger than 8 × 8 in the examples studied in Sec. III. Blocks are identified based on the expression of S_μν given in terms of cumulants according to MR-GWT. Terms of S_μν linear in Λ₂, Λ₃, or Λ₄ can yield a nonzero S_μν matrix element when some orbital indices constituting μ and ν match while indices not matching between μ and ν belong to the same geminal (allowing for a nonzero Λ element). In addition, the MR-GWT expression of S_μν involves terms quadratic in Λ₂, which necessitates watching for a bunch of four orbital indices constituting μ and ν belonging to one geminal and the other four again belonging to one geminal.

Determining zero eigenvalues of the overlap requires a threshold, which was set at 10⁻⁴ by Kong et al.⁴⁵ and by Hanauer and Köhn.³⁹ A second threshold of 10⁻² was applied by Kong et al. while system dependent thresholds on the order of 10⁻²–10⁻¹ were reported for CT.⁴³ In the calculations reported in Sec. III, a threshold of 10⁻¹⁰ is used uniformly. Despite the zero overlap threshold being comparatively small, no ill-effect has been detected on the potential surfaces explored in Sec. III. A numerical example on the role of redundancy affected amplitudes is given in Sec. III A 1.

The redundancy treatment described above applies even in the ultimate case where no overlap eigenvalue falls above the zero threshold either for a block of S or for an individual excitation. These excitations are eventually eliminated as the Moore–Penrose inverse of the matrix composed of $⟨) {\tilde{ψ}}_{ν} | H_{N} | ψ_{μ} (⟩$ is zero (cf. Appendix C). For computational economy, it is still practical to discard such excitations in advance when they can be easily identified. For this reason, core (n_p = 1) orbitals are excluded as a, b and virtual (n_p = 0) orbitals are not allowed as i, j in Eq. (8).

F. Size extensivity

Size extensivity follows directly from the connected construction of the theory, assuming a separation of the system for subsystems exhibiting zero inter-system cumulants, e.g., breaking up a system for two closed shell fragments. Excitations affected by redundancy are necessarily connected by cumulants (cf. the MR-GWT expansion of overlap matrix elements); redundancy treatment is therefore not expected to interfere with separability in lack of inter-system cumulants.

Separation of the system for subsystems that remain spin coupled is a more difficult question due to inter-system spin cumulants connecting the fragments.^36,37 This question has been addressed by Mukherjee⁹³ in connection with a JM based, unitary group adapted MRCC. The analysis being non-trivial, we wish to not address this topic here. At the same time, it appears worth to remark that taking APSG as reference, it is single covalent bond dissociation that deserves inspection from this point of view. Since multiple covalent bonds are not separated correctly by APSG,⁹⁴ such situations are better analyzed assuming, e.g., an unrestricted geminal reference.^95,96

III. NUMERICAL ILLUSTRATIONS

Numerical calculations assuming a GVB reference function are reported below. Starting orbitals of the GVB optimization process are generated by the Foster–Boys localization procedure.⁹⁷ Valence geminals corresponding to bonding and lone electron pairs are assigned two orbitals. Core orbitals belonging to inner shells are kept at the HF level, the two electrons being assigned a single spatial orbital. The orbitals underlying the calculations are those diagonalizing the GVB 1-RDM, exhibiting fractional occupation number when belonging to a two-dimensional geminal. Orbitals in the subspace of degenerate occupation numbers (i.e., cores characterized by n_i = 1 and virtual orbitals with n_a = 0) do no need further specification since the theory is unitary invariant in these two subspaces.

Numerical results obtained with GVB-rCCD are shown in parallel with the GVB-ERPA method of Pernal,⁹⁸ based on the same reference. This comparison is prompted by the fact that the two theories are related⁶⁴ and they are of comparable computational cost. Thresholds associated with pruning and redundancy treatment in rCCD are according to Secs. II D and II E, respectively. In the case of the GVB-ERPA, the categorization of orbitals as core, active, and virtual is performed based on the occupation number thresholds 0.995 and 0.005 separating the core/active and the active/virtual set, respectively. Single excitations allowed are of the type core → virt, core → act, act → virt, as well as act → act. The latter type requires further attention. Excitation p → q is in effect if n_q/n_p ≤ 0.99 holds when p and q belong to different geminals, while n_q ≤ n_p is required when p and q belong to the same geminal. In addition, ERPA amplitude vectors corresponding to complex excitation energy are discarded, with a threshold of 10⁻¹⁵ set on the imaginary part. These technical details are essentially in line with the procedure reported in Refs. 70, 98, and 99.

The computational cost of iterating the rCCD amplitude equation [Eq. (14)] roughly scales as $O ({(N_{c} + N_{a})}^{3} {(N_{a} + N_{v})}^{3})$ ⁠, where N_c is the number of core orbitals (with n_i = 1), N_a is the number of orbitals of fractional occupation number, and $N_{v}$ is the number of virtual orbitals (with n_a = 0). The cost of solving the ERPA equations necessary for ERPA-APSG similarly scales with the third power of the number of single excitations allowed to enter the excitation operator ansatz [a number that can be approximated from above by (N_c + N_a) (N_a + $N_{v}$ ⁠)].

Working equations of GVB-rCCD are Eqs. (10) and (14) with the tensor elements A and B given in Appendix A. In some of the examples, neglect of terms of F and B based on the cumulant rank is investigated. The cumulants retained are indicated in the acronym of these results, i.e., GVB-rCCD(Λ₁, Λ₂) refer to the exclusion of Λ₃ and Λ₄ involving terms from Eq. (10) and from the underlying amplitude equation (14).

Whenever the basis set permits, Full CI (FCI) results are reported as benchmarks. For larger basis sets, the CCSD(T) benchmark is used when strong correlation is negligible, while in the genuine multireference situation of N₂ dissociation, the semi-stochastic heat-bath CI (shCI) algorithm of Sharma et al.^100,101 serves as a benchmark. An in-house implementation of Olsen’s algorithm¹⁰² was used for obtaining FCI results; CCSD(T) energies were calculated with the Gaussian 09 software.¹⁰³

A. Covalent bond breaking and formation

Breaking and forming covalent bonds represent test cases of the genuine multireference situation. Among the examples shown, dissociation of a single covalent bond is captured properly by GVB, providing ideal circumstances for the assessment of MR based correlation correction schemes. Examples where different regions of the PES can be characterized by different Lewis structures are less favorable from the point of view of GVB. Deficiencies are therefore anticipated that might not be attributed exclusively to the correlation correction scheme in these situations. Such difficult cases are nevertheless of interest and are investigated.

As examples of the latter, difficult types dissociating multiple covalent bonds and elongation of two or more single covalent bonds attached to a common atom are provided. Spin recoupling, taking place during these processes, is not described properly by perfect pairing.⁹⁴ While the flaw is usually not apparent on the GVB PES, it may become an issue when devising correlation corrections built upon GVB.^99,104 The performance of GVB-rCCD is examined in two examples of this sort: symmetric dissociation of the H₂O molecule and the dissociation of the N₂ molecule.

Square to rectangle transformation of the H₄ system represents another difficult example. The problem here is apparent on the GVB PES, since two solutions corresponding to two Lewis structures cross at the square geometry. This produces an incorrect cusp at square geometry on the minimal energy PES. Insertion of Be into H₂ leading to a BeH₂ molecule is another example for a change in the characteristic Lewis structure taking place along the process. Since the work of Purvis et al.,¹⁰⁵ this has been a standard test system of CC methods.³⁰

1. H₂O molecule

The results on the elongation of a single OH bond of the H₂O molecule displayed in Fig. 1 demonstrate that both rCCD and ERPA corrected PES are significantly improved as compared with the GVB PES. The error of rCCD and ERPA is in the same order of magnitude, the former being somewhat smaller and both varying a couple of mE_h in the geometry range studied. There appears a small hump on the order of 0.1 mE_h with the maximum around 3 Å on the GVB-rCCD total energy curve, which is not discernible in Fig. 1. It may be worth to remind that such a behavior is not uncommon in the case of RPA-related methods.¹⁰⁶

FIG. 1.

View large Download slide

Single bond dissociation of the H₂O molecule in the 6-31G^* basis set at α_OHO = 104.5° bond angle and $R_{O H}^{'}$ = 1.000 Å bond length for the non-dissociating bond. All valence GVB based rCCD and ERPA corrected results are displayed; FCI serves as the benchmark. The total energy and energy difference computed as ΔE = E − E _FCI are displayed.

Results on the symmetric OH bond stretch of the H₂O molecule are shown in Figs. 2 and 3. This example is similar to the single OH bond stretch in that rCCD yields a greater portion of electron correlation than the ERPA correction. A hump is again present on the GVB-rCCD total energy curve on the order of a mE_h peaking at around 3 Å. At difference with the GVB based linearized CCD correction,¹⁰⁴ the GVB-rCCD curve shows no ill-effect in the spin recoupling region. The result obtained by neglecting Λ₃ including terms when calculating GVB-rCCD is also shown in Fig. 3. The latter approximation results in slightly smaller error than GVB-rCCD in the [1.5, 3.5] Å bond length region and a hump on the total energy curve, reduced to the order of 0.1 mE_h. The range of variation of the error is similar for the three methods displayed in Fig. 3.

FIG. 2.

View large Download slide

Symmetric dissociation of the H₂O molecule in the 6-31G^* basis set at α_OHO = 104.5° bond angle. All valence GVB based rCCD and ERPA corrected results are displayed; FCI serves as the benchmark.

FIG. 3.

View large Download slide

Energy difference taken with FCI for the methods and process shown in Fig. 2. The label GVB-rCCD(Λ₁, Λ₂, Λ₄) refers to omission of Λ₃ in the working equations.

A systematic study on the effect of cumulants involved in the expression of tensors F and B is provided on the example of the H₂O molecule in Fig. 4. The HF-like approximation where all cumulants of rank 2 and higher are neglected, GVB-rCCD(Λ₁), apparently inherits the divergent behavior of RHF-based rCCD.¹⁰⁷ Since cumulants decay fast with the increase in rank around equilibrium, it is comprehensible that inclusion of cumulants up to Λ₂ yields good results at around 1 Å in Fig. 4. This scheme, denoted GVB-rCCD(Λ₁, Λ₂), gives a finite but rather disappointing solution in the dissociating regime, where higher rank cumulants are not in general negligible compared to their lower rank counterpart. Accordingly, GVB-rCCD(Λ₁, Λ₂) produces a hump in the spin recoupling region and tends to an energy value lying off from FCI by some −50 mE_h in the dissociation limit. The scheme GVB-rCCD(Λ₁, Λ₂, Λ₃) produces a more expressed hump and tends essentially to the same limit as GVB-rCCD(Λ₁, Λ₂), since Λ₃ becomes negligible in the dissociated limit.³⁷ Allowing Λ₄ to enter, but discarding Λ₃, gives the scheme GVB-rCCD(Λ₁, Λ₂, Λ₄) that is already acceptable in Fig. 4. As apparent in Fig. 3, GVB-rCCD(Λ₁, Λ₂, Λ₄) stays somewhat more close to FCI than GVB-rCCD with the full inclusion of cumulants [denoted GVB-rCCD(Λ₁, Λ₂, Λ₃, Λ₄) in Fig. 4]. One can generally conclude that retaining Λ₄ in the energy formula is essential to obtain correct PES for covalent bond dissociation. At the same time, the success of GVB-rCCD(Λ₁, Λ₂, Λ₄) originates in a rather fortuitous than systematic error compensation taking place in this example. The scheme GVB-rCCD(Λ₁, Λ₂, Λ₄) outperforming GVB-rCCD does not occur in general.

FIG. 4.

View large Download slide

Total energy by various cumulant rank based cutoff schemes in GVB-rCCD for the process shown in Fig. 2. Cumulants retained in the working equations (10) and (14) are indicated in parentheses.

Symmetric dissociation of the H₂O molecule serves for one further test: examination of the numerical effect of various pruning strategies as well as the role of amplitudes affected by redundancy. The results are presented at two geometries in Table I, at around equilibrium and in the dissociated regime. The number of amplitudes retained in a given scheme, collected in Table I, has to be contested with 10 107, giving the number of original amplitudes when following the pruning described in Sec. II D. (Individual excitations or blocks of excitations with all overlap eigenvalues below the numerical threshold are already excluded in the original scheme, cf. the last paragraph of Sec. II D). The total number of amplitudes as well as the structure of the redundant blocks remains unchanged when stepping from R_OH = 1.0 Å to R_OH = 4.0 Å.

TABLE I.

Number of retained amplitudes and difference in total energy compared to GVB-rCCD when neglecting excitation types from Eq. (8). An example is provided by the symmetric dissociation of the H₂O molecule. See the legend of Fig. 2 for basis and geometry. The abbreviation “red” refers to redundancy affected amplitudes, “deexc” denotes deexcitation type transitions, and “exc.spec” and “dir.spec” stand for exchange and direct spectators, respectively. See text for more details.

		(E − E_rCCD) (E_h)
Amplitude type omitted	n_T	R_OH = 1.0 (Å)	R_OH = 4.0 (Å)
red	8079	9.562 × 10⁻⁴	8.284 × 10⁻⁴
deexc	7779	1.365 × 10⁻⁴	1.507 × 10⁻⁴
exc.spec	9067	7.764 × 10⁻⁵	6.154 × 10⁻⁵
dir.spec	9067	1.976 × 10⁻⁴	5.560 × 10⁻⁴
exc.spec, dir.spec	8127	3.216 × 10⁻⁴	2.429 × 10⁻⁴
exc.spec, dir.spec, deexc	6415	4.055 × 10⁻⁴	3.624 × 10⁻⁴

		(E − E_rCCD) (E_h)
Amplitude type omitted	n_T	R_OH = 1.0 (Å)	R_OH = 4.0 (Å)
red	8079	9.562 × 10⁻⁴	8.284 × 10⁻⁴
deexc	7779	1.365 × 10⁻⁴	1.507 × 10⁻⁴
exc.spec	9067	7.764 × 10⁻⁵	6.154 × 10⁻⁵
dir.spec	9067	1.976 × 10⁻⁴	5.560 × 10⁻⁴
exc.spec, dir.spec	8127	3.216 × 10⁻⁴	2.429 × 10⁻⁴
exc.spec, dir.spec, deexc	6415	4.055 × 10⁻⁴	3.624 × 10⁻⁴

Amplitude types omitted in the schemes examined in Table I involve exchange spectators, i.e., $E_{i j}^{a i}$ type transitions, abbreviated “exc.spec”; direct spectators, i.e., $E_{i j}^{i b}$ type transitions, abbreviated as “dir.spec”; and deexcitation type transitions, i.e., $E_{i j}^{a b}$ where n_i < n_a with I ≠ A or n_j < n_b with J ≠ B, abbreviated as “deexc.” In addition, elimination of amplitudes affected by redundancy (i.e., at least one overlap eigenvalue below the numerical zero in the corresponding block of S) is also examined, abbreviated as “red” in Table I. (Note that these amplitude types do not define disjoint sets of amplitudes. Some transitions, e.g., $E_{i i}^{a i}$ ⁠, may be considered both exchange and direct spectators. Moreover, spectators are heavily affected by redundancy.)

As given in Table I, roughly 20% of the amplitudes appear in redundant blocks of the matrix S in this example. Any attempt to converge the projection equations without redundancy treatment was unsuccessful. Eliminating redundancy affected amplitudes altogether introduces a change in energy, essentially on the mE_h level. Other schemes examined in Table I are accompanied by about an order of magnitude smaller energy change. Interestingly, the role of exchange spectators appears somewhat smaller than direct spectators. Spectators altogether amount to roughly 20% of the amplitudes. Deexcitation type transitions are somewhat more numerous, around 23% of the total amplitudes. The energetic role of deexcitations, however, appears smaller than the effect of all spectators. Variation of the energy change with geometry is the most pronounced for direct spectators, but even for this excitation type, it falls in the same order of magnitude for R_OH = 1.0 Å and 4.0 Å.

2. Bond dissociation of N₂

Dissociation of the N₂ molecule, presented in Fig. 5, is a second example for the multiple bond breaking process. The results are similar to those in the case of the H₂O molecule. The error of GVB, on the order of hundreds of mE_h, is reduced to the ten mE_h range by ERPA and by rCCD. The performance of rCCD is somewhat better than ERPA at the price of a hump discernible at around 2.4 Å. Parallelity of CCSD(T) is considerably better than either of the GVB based schemes in the equilibrium region, but it becomes unreliable beyond 1.8 Å.

FIG. 5.

View large Download slide

Dissociation of the N₂ molecule in the 6-31G^* basis set. All valence GVB based rCCD and ERPA corrected results are displayed along with HF based CCSD(T). The benchmark is provided by shCI.

3. H₄ system

The characteristic GVB PES with a cusp at α_HXH = 90° is displayed in Fig. 6 for the transformation of the H₄ system from rectangle to square geometry. The cusp shape is roughly conserved by ERPA as well as by rCCD, making a striking difference with the zero derivative at the square arrangement of the correct, FCI curve. While the incorrect curve shape is not alleviated much by the correction methods, the error of GVB is reduced considerably by rCCD, as shown in Fig. 6. Similar to the previous examples, rCCD and ERPA are comparable, the latter lagging behind rCCD by a couple of mE_h-s.

FIG. 6.

View large Download slide

Rectangle to square transformation of H₄ in the 6-31G^** basis set. A dummy atom at the center of mass is denoted by X, and distance XH is fixed at R_XH = 0.75 Å. The benchmark is provided by FCI.

4. BeH₂ system

Atom Be is placed at the origin, and coordinates of the hydrogen atoms in Å are (0, ±1.344, 0), (0, ±1.1005, 0.529), (0, ±0.8575, 1.058), (0, ±0.7355, 1.323), (0, ±0.6745, 1.455), (0, ±0.614, 1.588) (0, ±0.492, 1.852), (0, ±0.3705, 2.117), and (0, ±0.3705, 3.175), respectively, at points A, B, C, D, E, F, G, H, and I. In geometry points A–E, the two valence geminals can be assigned to the two Be–H bonds. Starting at point F, geminals of the lowest energy GVB solution can rather be identified as one residing on atom Be and another describing the H–H bond. Switching between two Lewis structures is not obvious from the results shown in Fig. 7 since the switching point is not explored with fine resolution for this system. Dunning’s DZ set¹⁰⁸ is taken for the hydrogen atoms. For beryllium, the basis of Purvis et al.¹⁰⁵ is used with the p function decontracted, leaving the most compact primitive (exponent 5.693 880) alone and contracting the remaining two into a second p function (exponents 1.555 630, 0.171 855 and coefficients 0.144 045, 0.949 692, respectively). A pair of single excitations with imaginary excitation energy is omitted from the ERPA correction shown in Fig. 7 at geometry points E and I. The rCCD calculation is also tinkered at point E: excitations $E_{p q}^{r r}$ and $E_{p q}^{s s}$ had to be omitted from T in order to achieve convergence. (At this point, orbitals p, r and q, s constitute one and the other BeH geminal, respectively.)

FIG. 7.

View large Download slide

Total energy of the BeH₂ system starting from a linear molecule (point A) and ending in a distant Be atom and a H₂ molecule (point I). All valence GVB based rCCD and ERPA corrected results are displayed; FCI serves as the benchmark. The basis set is of DZ quality, and geometry agrees with Ref. 105. The beryllium atom lies at the origin of the coordinate system. The coordinate z of the H atoms is plotted on axis x, labeled $d_{{Be - H}_{2}}$ ⁠. See text for more details.

Compared to the FCI benchmark, the GVB-ERPA describes a larger portion of the correlation than GVB-rCCD in this example. Figure 7 also reflects that the parallelity of GVB-rCCD is slightly better: difference from FCI lies in the range 3–29 mE_h for the GVB-ERPA and 13–30 mE_h for GVB-rCCD.

The BeH₂ system, as computed by Purvis and Bartlett,¹⁰⁵ has become a widely used test of performance in the presence of strong correlation. It offers the possibility of comparison with some MR based CC methods mentioned in the Introduction. Energy differences taken with FCI, collected in Table II, reflect that the performance of GVB-rCCD lags behind MR based CCSD schemes by roughly two orders of magnitude. This holds true regarding either of the benchmark methods displayed in Table II. Among these, CASCCSD(sw) is a single-reference based, pivot-dependent formulation of Lyakh and co-workers,¹¹ adapted for the MR situation following the basic idea of Ivanov and Adamowicz.^7,8 The abbreviation “sw” in the acronym refers to a swap of the pivot between $d_{{Be - H}_{2}} = 2.0$ bohrs and $d_{{Be - H}_{2}} = 2.5$ bohrs. Methods Mk-MRCCSD and ic-MRCCSD fall in the category of genuine MR-CC techniques formulated in the Hilbert-space. The JM parametrization is harnessed by Mk-MRCCSD developed by Mukherjee and co-workers,¹⁹ while ic-MRCCSD refers to the internally connected approach of Evangelista and Gauss.³⁸ Concerning the theoretical framework, it is the latter of the above three methods that bears the most kinship with GVB-rCCD. The uniformly huge error of GVB-rCCD as compared to either of the benchmarks is an obvious consequence of the ring approximation, which wipes out an ample number of terms of the CCSD equations and together with these a large portion of dynamical correlation. See, e.g., Refs. 109 and 110 to grab the extent of error of the ring approximation in the single reference based context. The assessment of Table II highlights that abandoning the ring approximation is essential for the present theory to become competitive with state-of-the-art MR-CC methods. Work along this line is currently in progress in our laboratory.

TABLE II.

Energy difference computed as ΔE = E − E_FCI for the BeH₂ molecule. The CASCCSD(sw) results are quoted from Ref. 111. The Mk-MRCCSD and ic-MRCCSD results are quoted from Ref. 112. The reference function underlying CASCCSD(sw), Mk-MRCCSD, and ic-MRCCSD is CAS(2,2). The basis set is of DZ quality, and geometry agrees with Ref. 112. The beryllium atom lies at the origin of the coordinate system. The coordinate z of the H atoms is labeled $d_{{Be - H}_{2}}$ ⁠. See text for more details.

$d_{{Be - H}_{2}}$	ΔE_GVB-rCCD	ΔE_CASCCSD(sw)	ΔE_Mk-MRCCSD	ΔE_ic-MRCCSD
(bohr)	(mE_h)	(mE_h)	(mE_h)	(mE_h)
2.000	13.013	0.024	0.609	0.239
2.500	13.811	0.044	0.320	0.385
2.750	29.772	0.072	−0.613	1.043
3.000	26.611	0.402	−1.140	0.534
3.250	23.241	0.275	0.058	0.303
3.500	24.675	0.206	0.257	0.222
3.750	24.559	0.158	0.227	0.169

$d_{{Be - H}_{2}}$	ΔE_GVB-rCCD	ΔE_CASCCSD(sw)	ΔE_Mk-MRCCSD	ΔE_ic-MRCCSD
(bohr)	(mE_h)	(mE_h)	(mE_h)	(mE_h)
2.000	13.013	0.024	0.609	0.239
2.500	13.811	0.044	0.320	0.385
2.750	29.772	0.072	−0.613	1.043
3.000	26.611	0.402	−1.140	0.534
3.250	23.241	0.275	0.058	0.303
3.500	24.675	0.206	0.257	0.222
3.750	24.559	0.158	0.227	0.169

B. Torsional motion of small molecules

The GVB model accounts for the interaction of chemical bonds only at the SCF level. The interest in torsional barriers described by geminal based methods^70,114 derives from this fact, since covalent bonds are essentially conserved, but their interaction changes during such processes. The three examples picked here are the bending motion of H₂O shown in Fig. 8, the torsional motion of H₂O₂ shown in Fig. 9, and the umbrella inversion of NH₃ shown in Fig. 10.

FIG. 8.

View large Download slide

Bending PES of the H₂O molecule in the Dunning-type cc-pVTZ basis set.¹¹³ Bond lengths are fixed at R_OH = 0.96 Å. All valence GVB based rCCD and ERPA corrected results are displayed; CCSD(T) serves as the benchmark. The energy difference is measured from the minimum of the respective curve in the top panel, while it is computed as ΔE = E − E_CCSD(T) in the bottom panel.

FIG. 9.

View large Download slide

Torsional PES of the H₂O₂ molecule in the Dunning-type cc-pVDZ basis set.¹¹³ Bond lengths are fixed at R_OH = 0.97 Å and R_OO = 1.45 Å; the value of the bond angle is α_OHO = 99.5°. The dihedral angle is labeled α _tors on axis x. All valence GVB based rCCD and ERPA corrected results are displayed; CCSD(T) serves as the benchmark. The energy difference is measured from the minimum of the respective curve in the top panel, while it is computed as ΔE = E − E_CCSD(T) in the bottom panel.

FIG. 10.

View large Download slide

Umbrella inversion of the NH₃ molecule in the cc-pVTZ basis set.¹¹³ Bond lengths are fixed at R_NH = 1.01 Å; X denotes a dummy atom in the center of mass. All valence GVB based rCCD and ERPA corrected results are displayed; CCSD(T) serves as the benchmark. The energy difference is measured from the minimum of the respective curve in the top panel, while it is computed as ΔE = E − E_CCSD(T) in the bottom panel.

Figures 8–10 uniformly show a considerable improvement over GVB both by rCCD and ERPA. The rCCD correction is closer to the CCSD(T) than ERPA for H₂O bending and H₂O₂ torsion. The quality of rCCD and ERPA is similar when compared with CCSD(T) on the example of NH₃ inversion; the error curves in the bottom panel of Fig. 10 fall in the kcal mol⁻¹ range. According to the top panel of Fig. 10, rCCD barriers are slightly better in this example.

C. Deprotonation energies

The deprotonation energy of the H₂O and CH₃OH molecules is calculated as the adiabatic energy difference between the neutral molecule and the anion. Geometry for both structures is optimized at the B3LYP/6-31G^* level. Comparing the results collected in Table III with the CCSD(T) benchmark, we see that the error of the dissociation energy is on the order of 10 kcal mol⁻¹ by GVB. This is reduced to a couple of kcal mol⁻¹ both by GVB-ERPA and GVB-rCCD. The accuracy of the two correction schemes does not appear remarkably different in these test cases.

TABLE III.

Deprotonation energies of the H₂O and CH₃OH molecules in the aug-cc-pVDZ basis set. All valence GVB based rCCD and ERPA corrected results are displayed; CCSD(T) serves as the benchmark.

		E_neutral (E_h)	E_anion (E_h)	ΔE (kcal mol⁻¹)
	GVB	−76.1060	−75.4627	403.62
H₂O	GVB-rCCD	−76.2719	−75.6453	393.18
	GVB-ERPA	−76.2580	−75.6245	397.51
	CCSD(T)	−76.2739	−75.6438	395.42
	GVB	−115.1717	−114.5318	401.58
CH₃OH	GVB-rCCD	−115.4435	−114.8234	389.12
	GVB-ERPA	−115.4269	−114.7986	394.25
	CCSD(T)	−115.4556	−114.8336	390.29

		E_neutral (E_h)	E_anion (E_h)	ΔE (kcal mol⁻¹)
	GVB	−76.1060	−75.4627	403.62
H₂O	GVB-rCCD	−76.2719	−75.6453	393.18
	GVB-ERPA	−76.2580	−75.6245	397.51
	CCSD(T)	−76.2739	−75.6438	395.42
	GVB	−115.1717	−114.5318	401.58
CH₃OH	GVB-rCCD	−115.4435	−114.8234	389.12
	GVB-ERPA	−115.4269	−114.7986	394.25
	CCSD(T)	−115.4556	−114.8336	390.29

IV. CONCLUSION

The ring approximation worked out and assessed with pilot numerical applications represents a new member in the family of icMR-CC methods. Assuming a GVB reference, it offers an $O (N^{6})$ way of incorporating a part of dynamical correlation into the wavefunction. The MR-rCCD correction brings a significant improvement over GVB, as illustrated on computed potential energy curves and energy differences. At the same time, MR-rCCD cannot amend the deficiencies of GVB originating in perfect pairing. The GVB model being relevant for ground states, applications resort to the ground state of the systems studied.

The performance of MR-rCCD is often comparable to the companion GVB-ERPA method. Small humps, on the order of mE_h, have been observed on the bond breaking PES obtained with MR-rCCD, similar to the single reference based RPA. A common feature of the rCCD and ERPA corrections is the necessity of pruning excitations admitted in the ansatz, which is eventually controlled by numerical thresholds. When following a PES, it has to be ensured that excitations are uniformly admitted in the geometry range considered. An advantageous feature of MR-rCCD over the GVB-ERPA is that N-representable RDMs are possible to be constructed based on the MR-rCCD corrected wavefunction.

Overlap treatment represents the source of a further numerical threshold in the procedure examined. This has, however, not been observed to cause any numerical issue in the pilot applications. Handling of the overlap based on frames, as performed here, is a generally applicable technique, irrespective of the truncation of the CC equations or the MR-CC flavor adopted. An important aspect of the overlap treatment based on frames is that no excitations are dropped, despite over-completeness. Since excitation is not eliminated, overlap treatment is not expected to undermine desirable properties, e.g., unitary invariance.

The intruder effect has been detected with the present MR-rCCD in one instance and deserves further exploration. Extensivity analysis in the case of spin-cumulants connecting the separated subsystems is likewise deferred to a follow-up work.

Contemplating on possible extensions of the theory, it has to be kept in mind that cumulants appearing as a consequence of MR-GWT represent the main obstacle. The advantage of a reference function exhibiting a pair structure cannot be overemphasized in this respect, since the majority of cumulant elements are wiped out just for this reason. Extension toward keeping the pair function character of the reference but abandoning the perfect pairing approximation is a possible way forward, since unrestriction at the geminal level (i.e., singlet-triplet mixed geminals) does not destroy the fragment structure of the cumulant nor does half-projection.¹¹⁵

Rank four is the highest rank of cumulants retained in the present scheme, essentially due to the ring approximation. Investigation in the line of systematic selection among relatively simple terms beyond the ring approximation is in progress and will be presented in a forthcoming report. Incorporation of single excitations is relatively straightforward and will be included together with beyond ring terms. We finally mention that cumulants of rank between 5 and N being zero with GVB (N standing for the number of electrons) may point to a GVB based full CCSD theory deserving exploration.

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

ACKNOWLEDGMENTS

This work was supported by the National Research, Development and Innovation Office (NKFIH), under Grant No. K115744. Support of the ÚNKP-19-3-I-ELTE-594 New National Excellence Program of the Ministry for Innovation and Technology is acknowledged by Á.M. This work was completed in the ELTE Institutional Excellence Program (No. 1783-3/2018/FEKUTSRAT) supported by the Hungarian Ministry of Human Capacities.

The FCI results were computed with a code implemented by Zoltán Rolik (Budapest University of Technology and Economics).

APPENDIX A: TENSORS OF THE SPIN-FREE rCCD EQUATIONS

Elements of tensor B arising from Eq. (11) fulfill

\begin{matrix} n_{i} n_{j} {\bar{n}}^{a} {\bar{n}}^{b} B_{i j}^{a b} & = & n_{i} n_{j} {\bar{n}}^{a} {\bar{n}}^{b} {\bar{v}}_{i j}^{a b} + \frac{1}{2} \sum_{p} ({\bar{n}}^{a} f_{p}^{a} Λ_{i j}^{p b} - n_{i} f_{i}^{p} Λ_{p j}^{a b} + {\bar{n}}^{b} f_{p}^{b} Λ_{i j}^{a p} - n_{j} f_{j}^{p} Λ_{i p}^{a b}) \\ + \frac{1}{2} \sum_{p, q} ({\bar{n}}^{a} {\bar{n}}^{b} v_{p q}^{a b} Λ_{i j}^{p q} + n_{i} n_{j} v_{i j}^{p q} Λ_{p q}^{a b} + n_{j} {\bar{n}}^{b} {\bar{v}}_{j p}^{b q} Λ_{i q}^{a p} + n_{i} {\bar{n}}^{a} {\bar{v}}_{i p}^{a q} Λ_{j q}^{b p}) - \frac{1}{2} \sum_{p, q} [n_{j} {\bar{n}}^{a} (v_{j p}^{q a} Λ_{i q}^{p b} + v_{j p}^{a q} Λ_{i q}^{b p}) + n_{i} {\bar{n}}^{b} (v_{i p}^{q b} Λ_{j q}^{p a} + v_{i p}^{b q} Λ_{j q}^{a p})] \\ + \frac{1}{4} \sum_{p, q, r, s} [{\bar{v}}_{p r}^{q s} Λ_{i q}^{a p} Λ_{j s}^{b r} - v_{p r}^{q s} (Λ_{i j}^{p b} Λ_{q s}^{a r} + Λ_{i j}^{a p} Λ_{q s}^{b r} + Λ_{i s}^{p r} Λ_{q j}^{a b} + Λ_{i q}^{a b} Λ_{j s}^{p r} + Λ_{i q}^{b t} Λ_{j s}^{a r})] \\ + \frac{1}{12} \sum_{p, q, r, s} v_{p r}^{q s} [Λ_{i j}^{p r} (2 Λ_{q s}^{a b} + Λ_{s q}^{a b}) + Λ_{i s}^{b p} (2 Λ_{j q}^{a r} + Λ_{j q}^{r a}) + Λ_{i s}^{p b} (2 Λ_{j q}^{r a} + Λ_{j q}^{a r})] \\ + \frac{1}{2} \sum_{p, q, r} ({\bar{n}}^{a} v_{p r}^{q a} Λ_{i j q}^{r b p} - n_{i} v_{i p}^{r q} Λ_{r j q}^{a b p} + {\bar{n}}^{b} v_{p r}^{q b} Λ_{i j q}^{a r p} - n_{j} v_{j p}^{r q} Λ_{i r q}^{a b p}) + \frac{1}{2} \sum_{p, q} f_{p}^{q} Λ_{i j q}^{a b p} + \frac{1}{4} \sum_{p, q, r, s} v_{p r}^{q s} Λ_{i j q s}^{a b p r} . \end{matrix}

(A1)

Comparing with the expression in Appendix B of the previous study,⁶⁴ lines 5–8 of Eq. (A1) are new. These terms appear since we presently allow for products of cumulants and cumulants beyond rank Λ₂. The terms in line 2 of Eq. (A1) were overlooked in the previous work. Concise amplitude equation, Eq. (14) assumes tensor B being a two-index quantity. Hiper index is formed of orbital indices aligned in column [e.g., ia and jb in Eq. (A1)] of B. The same holds for tensors A and t. Elements of B^T read as ${(B^{T})}_{i j}^{a b} = B_{j i}^{b a}$ ⁠. Symmetry relations applying to tensor B involve B^T = B. In accordance with

n^{i} n^{j} {\bar{n}}_{a} {\bar{n}}_{b} {(B^{†})}_{i j}^{a b} = \frac{1}{2} ⟨ : E_{b a}^{j i} : H_{N} ⟩,

we require ${(B^{†})}_{i j}^{a b} = {(B_{j i}^{b a})}^{*}$ ⁠. Due to the occupation numbers appearing in Eq. (A1), B is not self-adjoint, since ${(B_{j i}^{b a})}^{*} \neq B_{b a}^{j i}$ ⁠. A relation, however, exists between these two quantities, given by ${(B_{j i}^{b a})}^{*} = B_{b a}^{j i} (n_{b} \to - {\bar{n}}_{b}, n_{a} \to - {\bar{n}}_{a}, {\bar{n}}^{j} \to - n^{j}, {\bar{n}}^{i} \to - n^{i})$ ⁠, with $n_{b} \to - {\bar{n}}_{b}$ in parentheses referring to substituting n_b by $- {\bar{n}}_{b}$ ⁠. [Note that particle/hole occupation numbers appearing as a consequence of MR-GWT are tied to the upper/lower position of the corresponding index in Eq. (A1). Index conventions i, j and a, b are used just as a guide.] In the case of real quantities f, $v$ ⁠, and Λ, the simple relation ${(B^{†})}_{i j}^{a b} = B_{i j}^{a b}$ holds due to the structure of Eq. (A1).

In analogy with the above, elements of tensor F^† are introduced as

n^{i} {\bar{n}}_{a} {(F^{†})}_{i}^{a} = \frac{1}{2} ⟨ : E_{a}^{i} : H_{N} ⟩,

${(F^{†})}_{i}^{a}$ being given by

\begin{matrix} n^{i} {\bar{n}}_{a} {(F^{†})}_{i}^{a} & = & n^{i} {\bar{n}}_{a} f_{a}^{i} + \frac{1}{2} \sum_{p, q} f_{p}^{q} Λ_{q a}^{p i} \\ - \frac{1}{2} \sum_{p, q, r} (n^{i} v_{p r}^{i q} Λ_{a q}^{p r} - {\bar{n}}_{a} v_{a q}^{p r} Λ_{p r}^{i q}) + \frac{1}{4} \sum_{p, q, r, s} v_{p r}^{q s} Λ_{q s a}^{p r i} . \end{matrix}

(A2)

Similar to the case of tensor B, $F_{i}^{a} = {(F_{a}^{i})}^{*} ({\bar{n}}_{i} \to - n_{i}, n^{a} \to - {\bar{n}}^{a}) = {(F^{†})}_{a}^{i} ({\bar{n}}_{i} \to - n_{i}, n^{a} \to - {\bar{n}}^{a})$ holds for tensor F. As an example for the effect of occupation number substitutions, take the second term on the right hand side of Eq. (A2). It contributes $\frac{1}{2} \sum_{p, q} {(n^{a} {\bar{n}}_{i})}^{- 1} f_{p}^{q} Λ_{q i}^{p a}$ to ${(F^{†})}_{a}^{i}$ ⁠, while it becomes $\frac{1}{2} \sum_{p, q} {({\bar{n}}^{a} n_{i})}^{- 1} f_{p}^{q} Λ_{q i}^{p a}$ in $F_{i}^{a}$ ⁠.

Tensor A of the Riccati equation [Eq. (14)] is composed of the elements

\begin{matrix} A_{k b}^{c j} & = & δ_{k j} {\bar{n}}^{c} F_{b}^{c} (n_{b} \to - {\bar{n}}_{b}) - δ_{c b} n_{k} F_{k}^{j} ({\bar{n}}^{j} \to - n^{j}) \\ + {\bar{n}}^{c} n_{k} B_{k b}^{c j} ({\bar{n}}^{j} \to - n^{j}, n_{b} \to - {\bar{n}}_{b}), \end{matrix}

(A3)

with notation $n_{b} \to - {\bar{n}}_{b}$ and ${\bar{n}}^{j} \to - n^{j}$ explained above. Tensor A is not symmetric, ${(A^{T})}_{a k}^{i c} = A_{k a}^{c i} \neq A_{a k}^{i c}$ ⁠. A comparison with tensor A of Ref. 64 reveals that the Λ₃ involving last term of Eq. (A2) is new. Apart from that, new terms of A are generated by new terms of B given in Eq. (A1). A comparison with tensor A in Appendix B of Ref. 64 further reveals that δ_cb is missing from the first term and δ_kj is deleted from the second term on the right hand side of Eq. (A3). This ensures invariance of the theory to a rotation among core orbitals (n_i = 1) or among virtual orbitals (⁠ ${\bar{n}}_{a} = 1$ ⁠). These are the orbital rotations leaving the 1-RDM of the reference APSG function intact in general.

APPENDIX B: CUMULANTS OF THE APSG WAVEFUNCTION

Cumulants up to rank λ₄ were generated in the spin-orbital form by Eqs. (15a), (15b), (16a) and (16b). After spin-summation, elements up to rank 3 are given by

Λ_{q}^{p} = 2 δ_{q}^{p} n^{p},

Λ_{q s}^{p r} = 2 δ_{Q}^{P} c^{p} c_{q} δ^{p r} δ_{q s} - δ^{P R} (Λ_{q}^{p} Λ_{s}^{r} - \frac{1}{2} Λ_{s}^{p} Λ_{q}^{r}),

(B1)

\begin{matrix} Λ_{q s u}^{p r t} & = & - δ^{P R} δ^{R T} δ^{T P} (Λ_{q}^{p} Λ_{s}^{r} Λ_{u}^{t} - \frac{1}{2} Λ_{q}^{p} Λ_{u}^{r} Λ_{s}^{t} + \frac{1}{4} Λ_{u}^{p} Λ_{q}^{r} Λ_{s}^{t}) \\ - (\frac{1}{2} Λ_{u}^{p} Λ_{s}^{r} Λ_{q}^{t} + \frac{1}{4} Λ_{s}^{p} Λ_{u}^{r} Λ_{q}^{t} - \frac{1}{2} Λ_{s}^{p} Λ_{q}^{r} Λ_{u}^{t}) \\ - δ^{P R} (Λ_{q}^{p} Λ_{s u}^{r t} - \frac{1}{2} Λ_{s}^{p} Λ_{q u}^{r t} - \frac{1}{2} Λ_{u}^{p} Λ_{s q}^{r t}) \\ - δ^{R T} (Λ_{s}^{r} Λ_{u q}^{t p} - \frac{1}{2} Λ_{u}^{r} Λ_{s q}^{t p} - \frac{1}{2} Λ_{q}^{r} Λ_{u s}^{t p}) \\ - δ^{T P} (Λ_{u}^{t} Λ_{q s}^{p r} - \frac{1}{2} Λ_{q}^{t} Λ_{u s}^{p r} - \frac{1}{2} Λ_{s}^{t} Λ_{q u}^{p r}) . \end{matrix}

(B2)

The expression for Λ₄ is too lengthy to type out. It can still be manipulated rather easily with symbolic algebra programs. Spin-summations involved in Λ₄ can be performed with the help of relations⁹¹

\sum_{σ_{1}, σ_{2}, σ_{3}, σ_{4}} λ_{q σ_{1} s σ_{3}}^{p σ_{1} r σ_{2}} λ_{u σ_{2} w σ_{4}}^{t σ_{3} v σ_{4}} = \frac{1}{2} Λ_{q s}^{p r} Λ_{u w}^{t v}

(B3)

and

\begin{matrix} \sum_{σ_{1}, σ_{2}, σ_{3}, σ_{4}} λ_{q σ_{3} s σ_{4}}^{p σ_{1} r σ_{2}} λ_{u σ_{1} w σ_{2}}^{t σ_{3} v σ_{4}} \\ = \frac{1}{6} (2 Λ_{q s}^{p r} Λ_{u w}^{t v} + Λ_{q s}^{p r} Λ_{w u}^{t v} + Λ_{s q}^{p r} Λ_{u w}^{t v} + 2 Λ_{s q}^{p r} Λ_{w u}^{t v}) . \end{matrix}

(B4)

The structure of cumulants in the case of geminal dissociation merits a short comment. Taking the example of a GVB geminal with orbitals p and q assigned to it, geminal coefficients and occupation numbers in the dissociation limit take the form

c^{p} = - c^{q} = \frac{1}{\sqrt{2}}

(B5)

and

n^{p} = n^{q} = \frac{1}{2} .

(B6)

Nonzero elements of Λ₁ and Λ₂ are

Λ_{p}^{p} = Λ_{q}^{q} = 1

(B7)

and

Λ_{p p}^{p p} = Λ_{q q}^{q q} = Λ_{q p}^{p q} = \frac{1}{2},

(B8)

Λ_{p q}^{p q} = Λ_{q q}^{p p} = Λ_{p p}^{q q} = - 1 .

(B9)

Substitution of Eqs. (B7) and (B8) into Eq. (B2) reveals that all elements of Λ₃ are zero, in accordance with the general statement on the vanishing of odd rank cumulants in the limit of reaching perfect particle-hole symmetry.^37,76

Note that the requirement of Λ₁ being diagonal does not fix geminal orbitals in the dissociation limit due to the degeneracy in occupation numbers, cf. Eq. (B6). It is the diagonal form of the coefficient matrix [Eq. (B5)] that fixes the orbitals in a unique manner in the dissociated limit in the present formulation. The respective orbitals p and q are delocalized over the dissociated fragments in this situation.

APPENDIX C: DETAILS OF REDUNDANCY HANDLING

In order to justify the procedure described in Sec. II E, we introduce an auxiliary orthonormal basis obtained following Löwdin’s canonical orthogonalization scheme,¹¹⁶

| e_{μ} (⟩ = \sum_{ν} | ψ_{ν} (⟩ V_{ν μ} σ_{μ}^{- 1}, μ = 1, \dots, M,

(C1)

where V_νμ are elements of V containing the eigenvectors of S as columns, $σ_{μ} = \sqrt{Σ_{μ}}$ ⁠, and Σ_μ are eigenvalues of S, cf. Eq. (25). The matrix F (not to be mixed with tensor F of Appendix A) composed of elements F_νμ = ⟨ e_ν|ψ_μ⟩ provides the matrix of the so-called synthesis operator,⁷² its singular value decomposition (SVD) form reading as

F = U σ V^{†},

(C2)

where

σ = (\begin{matrix} σ^{'} 0 \end{matrix}),

σ′ = diag(σ_μ) is $M \times M$ ⁠, σ is $M \times N$ ⁠, and U is an $M \times M$ unit matrix. The matrix F obviously fulfills F^†F = S. The product in reverse order, FF^† = Σ′ obeying

A P_{U} \leq F F^{†} \leq B P_{U},

shows that functions ψ_μ constitute a frame. In the above equation, $A = min_{μ} Σ_{μ}^{'}$ and $B = max_{μ} Σ_{μ}^{'}$ provide the tightest bounds and P_U = UU^† is the $M \times M$ unit matrix.

Projector P_U can be expressed with F in the form

F {\tilde{F}}^{†} = P_{U},

(C3)

where

\tilde{F} = F S^{- 1} = U (\begin{matrix} σ^{' - 1} 0 \end{matrix}) V^{†}

(C4)

is the synthesis matrix of the so-called canonical dual frame. In Eq. (C4), S⁻¹ is understood as a Moore–Penrose inverse. The dual frame vectors, ${\tilde{ψ}}_{μ}$ ⁠, are associated with columns of $\tilde{F}$ ⁠. Based on Eq. (C3), ${\tilde{ψ}}_{μ}$ are used as bra vectors and ψ_μ are considered as ket vectors when writing the matrix vector form of the amplitude equations.

A justification of working with the over-complete set of frame vectors and their dual may be given by introducing a linearly independent set via extension of the matrix σ to dimension $N \times N$ ⁠, according to

\bar{σ} = (\begin{matrix} σ^{'} & 0 \\ 0 & 1 \end{matrix}) .

The SVD form of the extended synthesis matrix is constructed in analogy with Eq. (C2) as

\bar{F} = \bar{U} \bar{σ} V^{†},

(C5)

with $\bar{U}$ denoting the unit $N \times N$ matrix. Columns of $\bar{F}$ represent $N$ linearly independent vectors, which follows from the overlap matrix $\bar{S} = {\bar{F}}^{†} \bar{F}$ lacking zero eigenvalues. Inverse of $\bar{S}$ facilitates the construction of the dual set in the extended space as

\tilde{\bar{F}} = \bar{F} {\bar{S}}^{- 1} = \bar{U} {\bar{σ}}^{- 1} V^{†} .

Columns of $\bar{F}$ and $\tilde{\bar{F}}$ are biorthonormal, characterized by ${\tilde{\bar{F}}}^{†} \bar{F} = V V^{†}$ ⁠. In comparison, frame vectors and their dual are not biorthonormal, their product yielding

{\tilde{F}}^{†} F = S^{- 1} S = V (\begin{matrix} P_{U} & 0 \\ 0 & 0 \end{matrix}) V^{†} .

(C6)

The matrix of Eq. (C6), denoted as

R = V (\begin{matrix} P_{U} & 0 \\ 0 & 0 \end{matrix}) V^{†},

is obviously idempotent. It is a key quantity, connecting synthesis matrices via

(\begin{matrix} F \\ 0 \end{matrix}) = \bar{F} R,

(C7)

(\begin{matrix} \tilde{F} \\ 0 \end{matrix}) = \tilde{\bar{F}} R,

(C8)

with the bottom zero block on the lhs of the above equations being of dimension $(N - M) \times N$ ⁠.

We now proceed to the expansion of the wavefunction. With the help of frame vectors, one can write the $\bar{U}$ -component of Ψ^CCD as

P_{\bar{U}} | Ψ^{CCD} (⟩ = \sum_{μ} | ψ_{μ} (⟩ t_{μ} + O (t^{2}),

(C9)

while the extended frame vectors allow us to parametrize the same component as

P_{\bar{U}} | Ψ^{CCD} (⟩ = \sum_{μ} | {\bar{ψ}}_{μ} (⟩ {\bar{t}}_{μ} + O ({\bar{t}}^{2}),

(C10)

with $P_{\bar{U}}$ standing for the projector of the space spanned by ${{\bar{ψ}}_{μ}}_{μ = 1}^{N}$ ⁠, and ψ_μ in Eq. (C9) are associated with the columns of the matrix F padded with a zero $(N - M) \times N$ block at the bottom. Associating ${\bar{ψ}}_{ν}$ with the columns of $\bar{F}$ ⁠, the relation

| ψ_{μ} (⟩ = P_{\bar{U}} \sum_{ν} | {\bar{ψ}}_{ν} (⟩ R_{ν μ}

(C11)

holds based on Eq. (C7). Substituting Eq. (C11) into Eq. (C9) results in

P_{\bar{U}} | Ψ^{CCD} (⟩ = P_{\bar{U}} \sum_{μ} \sum_{ν} | {\bar{ψ}}_{ν} (⟩ R_{ν μ} t_{μ} + O (t^{2})

(C12a)

= P_{\bar{U}} \sum_{μ} | {\bar{ψ}}_{μ} (⟩ {(R t)}_{μ} + O (t^{2}) .

(C12b)

Comparing Eqs. (C12b) and (C10), one can deduce

R t = \bar{t},

(C13)

omitting terms quadratic in the amplitude.

Elaborating further on amplitudes, let us relate t and $\bar{t}$ to a third amplitude vector, corresponding to ${e_{μ}}_{μ = 1}^{M}$ denoted by t^e. We regard t^e a column vector of length $M$ ⁠. Based on Eq. (C1), the relation between t and t^e reads

t = V (\begin{matrix} σ^{' - 1} \\ 0 \end{matrix}) t^{e} .

(C14)

At the same time,

| e_{μ} (⟩ = \sum_{ν} | {\bar{ψ}}_{ν} (⟩ V_{ν μ} {\bar{σ}}_{μ}^{- 1}, μ = 1, \dots, N,

(C15)

facilitates to relate $\bar{t}$ and t^e as

\bar{t} = V (\begin{matrix} σ^{' - 1} & 0 \\ 0 & 1 \end{matrix}) (\begin{matrix} t^{e} \\ 0 \end{matrix}) .

(C16)

Comparison of Eqs. (C14) and (C16) reveals that not only Eq. (C13) holds but also that t of Eq. (C14) and $\bar{t}$ are the same. Note, however, that Eq. (C1) is not the only way to express ${e_{μ}}_{μ = 1}^{M}$ with ${ψ_{μ}}_{μ = 1}^{N}$ due to the over-completeness of the latter set. Consequently, Eq. (C14) is just one of the possibilities of generating t from t^e. The expression of ${e_{μ}}_{μ = 1}^{N}$ with ${{\bar{ψ}}_{μ}}_{μ = 1}^{N}$ according to Eq. (C15) is on the other hand unique. The equivalence of vectors t and $\bar{t}$ of Eqs. (C14) and (C16) can be interpreted as follows: Introduction of the linearly independent extended functions singles out one particular vector t, which is related to t^e according to Eq. (C14). This is the origin of removing the ill-defined nature of t while conserving an amplitude vector of length $N$ ⁠. Other vectors t, fulfilling Eq. (C9) and not matching $\bar{t}$ ⁠, can be related to it via Eq. (C13).

It is also instructive to project Eq. (C12b) with functions ${\tilde{\bar{ψ}}}_{ν}$ and project Eq. (C9) with the dual frame vectors ${\tilde{ψ}}_{ν}$ ⁠. In both cases, one obtains

{(R t)}_{ν} + O (t^{2}) = ⟨) {\tilde{\bar{ψ}}}_{ν} | Ψ^{CCD} (⟩ = ⟨) {\tilde{ψ}}_{ν} | Ψ^{CCD} (⟩,

based on Eqs. (C13) and (C6). This indicates that it is not necessary to explicitly construct $\tilde{\bar{F}}$ and $\bar{F}$ ⁠, since the same result is obtained when working with the set $\tilde{F}, F$ ⁠.

With this in hand, we consider the amplitude equations written in the form

\begin{matrix} \sum_{μ} ⟨) {\tilde{ψ}}_{ν} | H_{N} | ψ_{μ} (⟩ t_{μ} \\ = - ⟨) {\tilde{ψ}}_{ν} | H_{N} | Φ (⟩ - \frac{1}{2} \sum_{μ λ} ⟨) {\tilde{ψ}}_{ν} | H_{N} : E_{μ} E_{λ} : | Φ {(⟩}_{C} t_{μ} t_{λ} - O (t^{3}) . \end{matrix}

(C17)

Coupled cluster iteration usually proceeds via arranging coupling terms $\sum_{μ \neq ν} ⟨) {\tilde{ψ}}_{ν} | H_{N} | ψ_{μ} (⟩ t_{μ}$ on the rhs of Eq. (C17) and getting an update for t_μ via division by $⟨) {\tilde{ψ}}_{μ} | H_{N} | ψ_{μ} (⟩$ ⁠. Instead of this route, let us consider updating t_μ via keeping coupling terms on the lhs and inversion of the coefficient matrix given by the elements $⟨) {\tilde{ψ}}_{ν} | H_{N} | ψ_{μ} (⟩$ ⁠. In the case where the set of ${ψ_{μ}}_{μ = 1}^{N}$ is over-complete, the coefficient matrix is invertible only in the Moore–Penrose sense. To make this transparent, let us rewrite $⟨) {\tilde{ψ}}_{ν} | H_{N} | ψ_{μ} (⟩$ based on Eqs. (C2) and (C4) to get

⟨) {\tilde{ψ}}_{ν} | H_{N} | ψ_{μ} (⟩ = \sum_{γ λ}^{M} V_{ν γ} σ_{γ}^{' - 1} α_{γ λ} σ_{λ}^{'} V_{μ λ}^{*}

(C18)

with

α_{γ λ} = ⟨) e_{γ} | H_{N} | e_{λ} (⟩, γ, λ = 1, \dots, M .

(C19)

Since the set ${e_{μ}}_{μ = 1}^{N}$ is linearly independent, the matrix α given by the elements in Eq. (C19) is invertible. The coefficient matrix of Eq. (C17)

V (\begin{matrix} σ^{' - 1} \\ 0 \end{matrix}) α (\begin{matrix} σ^{'} & 0 \end{matrix}) V^{†}

can be inverted in Moore–Penrose sense with the use of α⁻¹ as

V (\begin{matrix} σ^{' - 1} \\ 0 \end{matrix}) α^{- 1} (\begin{matrix} σ^{'} & 0 \end{matrix}) V^{†} .

(C20)

Introducing

\begin{matrix} β_{λ} (t) = ⟨) e_{λ} | H_{N} | Φ (⟩ + \frac{1}{2} \sum_{μ ρ} ⟨) e_{λ} | H_{N} : E_{μ} E_{ρ} : | Φ {(⟩}_{C} t_{μ} t_{ρ} + O (t^{3}), \\ λ = 1, \dots, M, \end{matrix}

(C21)

the rhs of Eq. (C17) can be expressed as $- \sum_{λ}^{M} V_{ν λ} σ_{λ}^{' - 1} β_{λ} (t)$ ⁠, giving rise to

t = V (\begin{matrix} σ^{' - 1} \\ 0 \end{matrix}) α^{- 1} β (t)

(C22)

upon multiplying Eq. (C17) with Eq. (C20). Comparing Eqs. (C14) and (C22), one can deduce

t^{e} = α^{- 1} β (t) .

(C23)

Amplitude update based on the Moore–Penrose inverse of the coefficient matrix given by $⟨) {\tilde{ψ}}_{ν} | H_{N} | ψ_{μ} (⟩$ is hence equivalent to updating in the Löwdin-basis according to Eq. (C23) and transforming t^e with the help of Eq. (C14). The treatment is general in the sense that it is applicable to any excitation level or approximate form of the CC equations. Merely, matrices α and β(t) need to be adapted to the particular CC scheme.

REFERENCES

1.

R. J.

Bartlett

and

M.

Musiał

,

Rev. Mod. Phys.

79

,

291

(

2007

).

https://doi.org/10.1103/revmodphys.79.291

Google Scholar

Crossref

2.

K.

Raghvachari

,

G. W.

Trucks

,

J. A.

Pople

, and

M.

Head-Gordon

,

Chem. Phys. Lett.

157

,

479

(

1989

).

https://doi.org/10.1016/s0009-2614(89)87395-6

Google Scholar

Crossref

3.

Z.

Rolik

,

L.

Szegedy

,

I.

Ladjánszki

,

B.

Ladóczki

, and

M.

Kállay

,

J. Chem. Phys.

139

,

094105

(

2013

).

https://doi.org/10.1063/1.4819401

Google Scholar

Crossref

PubMed

4.

L.

Gyevi-Nagy

,

M.

Kállay

, and

P. R.

Nagy

,

J. Chem. Theory Comput.

16

,

366

(

2020

).

https://doi.org/10.1021/acs.jctc.9b00957

Google Scholar

Crossref

PubMed

5.

N.

Oliphant

and

L.

Adamowicz

,

J. Chem. Phys.

94

,

1229

(

1991

).

https://doi.org/10.1063/1.460031

Google Scholar

Crossref

6.

P.

Piecuch

,

N.

Oliphant

, and

L.

Adamowicz

,

J. Chem. Phys.

99

,

1875

(

1993

).

https://doi.org/10.1063/1.466179

Google Scholar

Crossref

7.

V. V.

Ivanov

and

L.

Adamowicz

,

J. Chem. Phys.

112

,

9258

(

2000

).

https://doi.org/10.1063/1.481547

Google Scholar

Crossref

8.

L.

Adamowicz

,

J.-P.

Malrieu

, and

V. V.

Ivanov

,

J. Chem. Phys.

112

,

10075

(

2000

).

https://doi.org/10.1063/1.481649

Google Scholar

Crossref

9.

P.

Piecuch

,

K.

Kowalski

,

I. S. O.

Pimienta

, and

M. J.

McGuire

,

Int. Rev. Phys. Chem.

21

,

527

(

2002

).

https://doi.org/10.1080/0144235021000053811

Google Scholar

Crossref

10.

Z.

Rolik

and

M.

Kállay

,

J. Chem. Phys.

148

,

124108

(

2018

).

https://doi.org/10.1063/1.5004971

Google Scholar

Crossref

PubMed

11.

D. I.

Lyakh

,

V. V.

Ivanov

, and

L.

Adamowicz

,

J. Chem. Phys.

122

,

024108

(

2005

).

https://doi.org/10.1063/1.1824897

Google Scholar

Crossref

PubMed

12.

P.

Piecuch

,

S.

Hirata

,

K.

Kowalski

,

P.-D.

Fan

, and

T. L.

Windus

,

Int. J. Quantum Chem.

106

,

79

(

2006

).

https://doi.org/10.1002/qua.20753

Google Scholar

Crossref

13.

W.

Li

,

P.

Piecuch

,

J. R.

Gour

, and

S.

Li

,

J. Chem. Phys.

131

,

114109

(

2009

).

https://doi.org/10.1063/1.3218842

Google Scholar

Crossref

PubMed

14.

Z.

Rolik

and

M.

Kállay

,

J. Chem. Phys.

141

,

134112

(

2014

).

https://doi.org/10.1063/1.4896703

Google Scholar

Crossref

PubMed

15.

I.

Lindgren

and

J.

Morrison

,

Atomic Many-Body Theory

(

Springer

,

Berlin

,

1986

).

Google Scholar

Crossref

16.

I.

Shavitt

and

R. J.

Bartlett

,

Many-Body Methods in Chemistry and Physics

(

Cambridge University Press

,

Cambridge

,

2009

).

Google Scholar

Crossref

17.

B.

Jeziorski

and

H. J.

Monkhorst

,

Phys. Rev. A

24

,

1668

(

1981

).

https://doi.org/10.1103/physreva.24.1668

Google Scholar

Crossref

18.

S. A.

Kucharski

and

R. J.

Bartlett

,

J. Chem. Phys.

95

,

8227

(

1991

).

https://doi.org/10.1063/1.461301

Google Scholar

Crossref

19.

U. S.

Mahapatra

,

B.

Datta

, and

D.

Mukherjee

,

J. Chem. Phys.

110

,

6171

(

1999

).

https://doi.org/10.1063/1.478523

Google Scholar

Crossref

20.

J.

Paldus

,

J.

Pittner

, and

P.

Čársky

, “

Multireference coupled-cluster methods: Recent developments

,” in

Recent Progress in Coupled Cluster Methods: Theory and Applications

, edited by

P.

Cársky

,

J.

Paldus

, and

J.

Pittner

(

Springer Netherlands

,

Dordrecht

,

2010

), pp.

455

–

489

.

Google Scholar

Crossref

21.

D.

Mukhopadhyay

,

B.

Datta (nee Kundu)

, and

D.

Mukherjee

,

Chem. Phys. Lett.

197

,

236

(

1992

).

https://doi.org/10.1016/0009-2614(92)85761-x

Google Scholar

Crossref

22.

L.

Meissner

,

Chem. Phys. Lett.

255

,

244

(

1996

).

https://doi.org/10.1016/0009-2614(96)00400-9

Google Scholar

Crossref

23.

L.

Meissner

and

M.

Musiał

, “

Intermediate Hamiltonian formulations of the fock-space coupled-cluster method: Details, comparisons, examples

,” in

Recent Progress in Coupled Cluster Methods: Theory and Applications

, edited by

P.

Cársky

,

J.

Paldus

, and

J.

Pittner

(

Springer Netherlands

,

Dordrecht

,

2010

), pp.

395

–

428

.

Google Scholar

Crossref

24.

I.

Lindgren

,

Int. J. Quantum Chem. Symp.

12

,

33

(

1978

).

https://doi.org/10.1002/qua.560140804

Google Scholar

Crossref

25.

D.

Mukherjee

, “

Cluster expansion from a multi-determinant reference function: Wick reduction formula

,” in

Recent Progress in Many Body Theories

(

Plenum Press

,

New York

,

1995

), Vol. 4, pp.

127

–

134

.

Google Scholar

26.

D.

Mukherjee

,

Chem. Phys. Lett.

274

,

561

(

1997

).

https://doi.org/10.1016/s0009-2614(97)00714-8

Google Scholar

Crossref

27.

W.

Kutzelnigg

and

D.

Mukherjee

,

J. Chem. Phys.

107

,

432

(

1997

).

https://doi.org/10.1063/1.474405

Google Scholar

Crossref

28.

U. S.

Mahapatra

,

B.

Datta

,

B.

Bandyopadhyay

, and

D.

Mukherjee

,

Adv. Quantum Chem.

30

,

163

(

1998

).

https://doi.org/10.1016/s0065-3276(08)60507-9

Google Scholar

Crossref

29.

D. I.

Lyakh

,

M.

Musiał

,

V. F.

Lotrich

, and

R. J.

Bartlett

,

Chem. Rev.

112

,

182

(

2012

).

https://doi.org/10.1021/cr2001417

Google Scholar

Crossref

PubMed

30.

A.

Köhn

,

M.

Hanauer

,

L. A.

Mück

,

T.-C.

Jagau

, and

J.

Gauss

,

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

3

,

176

(

2013

).

https://doi.org/10.1002/wcms.1120

Google Scholar

Crossref

31.

F. A.

Evangelista

,

J. Chem. Phys.

149

,

030901

(

2018

).

https://doi.org/10.1063/1.5039496

Google Scholar

Crossref

PubMed

32.

M.

Hanrath

,

J. Chem. Phys.

128

,

154118

(

2008

).

https://doi.org/10.1063/1.2899645

Google Scholar

Crossref

PubMed

33.

L.

Kong

,

Int. J. Quantum Chem.

109

,

441

(

2009

).

https://doi.org/10.1002/qua.21822

Google Scholar

Crossref

34.

W.

Kutzelnigg

,

K. R.

Shamasundar

, and

D.

Mukherjee

,

Mol. Phys.

108

,

433

(

2010

).

https://doi.org/10.1080/00268970903547926

Google Scholar

Crossref

35.

L.

Kong

,

M.

Nooijen

, and

D.

Mukherjee

,

J. Chem. Phys.

132

,

234107

(

2010

).

https://doi.org/10.1063/1.3439395

Google Scholar

Crossref

PubMed

36.

J. M.

Herbert

,

Int. J. Quantum Chem.

107

,

703

(

2007

).

https://doi.org/10.1002/qua.21160

Google Scholar

Crossref

37.

M.

Hanauer

and

A.

Köhn

,

Chem. Phys.

401

,

50

(

2012

).

https://doi.org/10.1016/j.chemphys.2011.09.024

Google Scholar

Crossref

38.

F. A.

Evangelista

and

J.

Gauss

,

J. Chem. Phys.

134

,

114102

(

2011

).

https://doi.org/10.1063/1.3559149

Google Scholar

Crossref

PubMed

39.

M.

Hanauer

and

A.

Köhn

,

J. Chem. Phys.

134

,

204111

(

2011

).

https://doi.org/10.1063/1.3592786

Google Scholar

Crossref

PubMed

40.

M.

Hanauer

and

A.

Köhn

,

J. Chem. Phys.

137

,

131103

(

2012

).

https://doi.org/10.1063/1.4757728

Google Scholar

Crossref

PubMed

41.

F. A.

Evangelista

,

M.

Hanauer

,

A.

Köhn

, and

J.

Gauss

,

J. Chem. Phys.

136

,

204108

(

2012

).

https://doi.org/10.1063/1.4718704

Google Scholar

Crossref

PubMed

42.

T.

Yanai

and

G. K.-L.

Chan

,

J. Chem. Phys.

124

,

194106

(

2006

).

https://doi.org/10.1063/1.2196410

Google Scholar

Crossref

PubMed

43.

E.

Neuscamman

,

T.

Yanai

, and

G. K.-L.

Chan

,

J. Chem. Phys.

132

,

024106

(

2010

).

https://doi.org/10.1063/1.3274822

Google Scholar

Crossref

PubMed

44.

Z.

Chen

and

M. R.

Hoffmann

,

J. Chem. Phys.

137

,

014108

(

2012

).

https://doi.org/10.1063/1.4731634

Google Scholar

Crossref

PubMed

45.

L.

Kong

,

K. R.

Shamasundar

,

O.

Demel

, and

M.

Nooijen

,

J. Chem. Phys.

130

,

114101

(

2009

).

https://doi.org/10.1063/1.3089302

Google Scholar

Crossref

PubMed

46.

P.

Jeszenszki

,

P. R.

Surján

, and

Á.

Szabados

,

J. Chem. Phys.

138

,

124110

(

2013

).

https://doi.org/10.1063/1.4795436

Google Scholar

Crossref

PubMed

47.

M.

Nooijen

and

R. J.

Bartlett

,

J. Chem. Phys.

104

,

2652

(

1996

).

https://doi.org/10.1063/1.471010

Google Scholar

Crossref

48.

M.

Nooijen

,

J. Chem. Phys.

104

,

2638

(

1996

).

https://doi.org/10.1063/1.470988

Google Scholar

Crossref

49.

D.

Datta

,

L.

Kong

, and

M.

Nooijen

,

J. Chem. Phys.

134

,

214116

(

2011

).

https://doi.org/10.1063/1.3592494

Google Scholar

Crossref

PubMed

50.

C. L.

Janssen

and

H. F.

Schaefer

,

Theor. Chim. Acta

79

,

1

(

1991

).

https://doi.org/10.1007/bf01113327

Google Scholar

Crossref

51.

M.

Kállay

and

P. R.

Surján

,

J. Chem. Phys.

115

,

2945

(

2001

).

https://doi.org/10.1063/1.1383290

Google Scholar

Crossref

52.

S.

Hirata

,

J. Phys. Chem. A

107

,

9887

(

2003

).

https://doi.org/10.1021/jp034596z

Google Scholar

Crossref

53.

A.

Köhn

,

J. Chem. Phys.

130

,

104104

(

2009

).

https://doi.org/10.1063/1.3079543

Google Scholar

Crossref

PubMed

54.

E.

Neuscamman

,

T.

Yanai

, and

G. K.-L.

Chan

,

J. Chem. Phys.

130

,

124102

(

2009

).

https://doi.org/10.1063/1.3086932

Google Scholar

Crossref

PubMed

55.

O.

Demel

,

D.

Datta

, and

M.

Nooijen

,

J. Chem. Phys.

138

,

134108

(

2013

).

https://doi.org/10.1063/1.4796523

Google Scholar

Crossref

PubMed

56.

D.

Jana

,

U. S.

Mahapatra

, and

D.

Mukherjee

,

Chem. Phys. Lett.

353

,

100

(

2002

).

https://doi.org/10.1016/s0009-2614(01)01500-7

Google Scholar

Crossref

57.

D.

Datta

and

D.

Mukherjee

,

J. Chem. Phys.

131

,

044124

(

2009

).

https://doi.org/10.1063/1.3185356

Google Scholar

Crossref

PubMed

58.

R.

Maitra

,

D.

Sinha

, and

D.

Mukherjee

,

J. Chem. Phys.

137

,

024105

(

2012

).

https://doi.org/10.1063/1.4731341

Google Scholar

Crossref

PubMed

59.

D.

Sinha

,

R.

Maitra

, and

D.

Mukherjee

,

J. Chem. Phys.

137

,

094104

(

2012

).

https://doi.org/10.1063/1.4747702

Google Scholar

Crossref

PubMed

60.

G.

Chan

,

J. Chem. Phys.

127

,

104107

(

2007

).

https://doi.org/10.1063/1.2761870

Google Scholar

Crossref

PubMed

61.

D.

Datta

and

M.

Nooijen

,

J. Chem. Phys.

137

,

204107

(

2012

).

https://doi.org/10.1063/1.4766361

Google Scholar

Crossref

PubMed

62.

D.

Sinha

,

R.

Maitra

, and

D.

Mukherjee

,

Comput. Theor. Chem.

1003

,

62

(

2013

).

https://doi.org/10.1016/j.comptc.2012.09.035

Google Scholar

Crossref

63.

J.

Čížek

,

J. Chem. Phys.

45

,

4256

(

1966

).

https://doi.org/10.1063/1.1727484

Google Scholar

Crossref

64.

Á.

Szabados

and

Á.

Margócsy

,

Mol. Phys.

115

,

2731

(

2017

).

https://doi.org/10.1080/00268976.2017.1317111

Google Scholar

Crossref

65.

N.

Ostlund

and

M.

Karplus

,

Chem. Phys. Lett.

11

,

450

(

1971

).

https://doi.org/10.1016/0009-2614(71)80382-2

Google Scholar

Crossref

66.

A. E.

Hansen

and

T. D.

Bouman

,

Mol. Phys.

37

,

1713

(

1979

).

https://doi.org/10.1080/00268977900101271

Google Scholar

Crossref

67.

G. E.

Scuseria

,

T. M.

Henderson

, and

D. C.

Sorensen

,

J. Chem. Phys.

129

,

231101

(

2008

).

https://doi.org/10.1063/1.3043729

Google Scholar

Crossref

PubMed

68.

K.

Chatterjee

and

K.

Pernal

,

J. Chem. Phys.

137

,

204109

(

2012

).

https://doi.org/10.1063/1.4766934

Google Scholar

Crossref

PubMed

69.

K.

Pernal

,

J. Chem. Theory Comput.

10

,

4332

(

2014

).

https://doi.org/10.1021/ct500478t

Google Scholar

Crossref

PubMed

70.

E.

Pastorczak

and

K.

Pernal

,

Phys. Chem. Chem. Phys.

17

,

8622

(

2015

).

https://doi.org/10.1039/c4cp05958a

Google Scholar

Crossref

PubMed

71.

V. I.

Morgenshtern

and

H.

Bölcskei

, in

Mathematical Foundations for Signal Processing, Communications, and Networking

, edited by

T. C. E.

Serpedin

and

D.

Rajan

(

CRC Press

,

2012

), Chap. 20, pp.

737

–

789

.

Google Scholar

72.

P. G.

Casazza

,

G.

Kutyniok

, and

F.

Philipp

, in

Finite Frames: Theory and Applications

, edited by

P. G.

Casazza

,

G.

Kutyniok

, and

B.

Basel

(

Springer Scient and Business Media

,

2012

), Chap. 1, pp.

1

–

54

.

Google Scholar

73.

D. W.

Small

and

M.

Head-Gordon

,

J. Chem. Phys.

130

,

084103

(

2009

).

https://doi.org/10.1063/1.3069296

Google Scholar

Crossref

PubMed

74.

P. R.

Surján

,

Top. Curr. Chem.

203

,

63

(

1999

).

https://doi.org/10.1007/3-540-48972-x_4

Google Scholar

Crossref

75.

P.

Jeszenszki

,

P. R.

Nagy

,

T.

Zoboki

,

Á.

Szabados

, and

P. R.

Surján

,

Int. J. Quantum Chem.

114

,

1048

(

2014

).

https://doi.org/10.1002/qua.24634

Google Scholar

Crossref

76.

W.

Kutzelnigg

and

D.

Mukherjee

,

J. Chem. Phys.

110

,

2800

(

1999

).

https://doi.org/10.1063/1.478189

Google Scholar

Crossref

77.

W.

Kutzelnigg

,

Chem. Phys.

401

,

119

(

2012

).

https://doi.org/10.1016/j.chemphys.2011.10.020

Google Scholar

Crossref

78.

F. W.

Bobrowicz

and

W. A.

Goddard

III, in

Methods of Electronic Structure Theory

, edited by

H. F.

Schaefer

III (

Plenum

,

New York

,

1977

), p.

79

.

Google Scholar

Crossref

79.

P. R.

Surján

,

Int. J. Quantum Chem.

55

,

109

(

1995

).

https://doi.org/10.1002/qua.560550205

Google Scholar

Crossref

80.

E.

Rosta

and

P. R.

Surján

,

J. Chem. Phys.

116

,

878

(

2002

).

https://doi.org/10.1063/1.1427918

Google Scholar

Crossref

81.

S.

Li

,

J.

Ma

, and

Y.

Jiang

,

J. Chem. Phys.

118

,

5736

(

2003

).

https://doi.org/10.1063/1.1555983

Google Scholar

Crossref

82.

S.

Li

,

J. Chem. Phys.

120

,

5017

(

2004

).

https://doi.org/10.1063/1.1646355

Google Scholar

Crossref

PubMed

83.

V. A.

Rassolov

,

F.

Xu

, and

S.

Garashchuk

,

J. Chem. Phys.

120

,

10385

(

2004

).

https://doi.org/10.1063/1.1738110

Google Scholar

Crossref

PubMed

84.

G. J. O.

Beran

,

M.

Head-Gordon

, and

S. R.

Gwaltney

,

J. Chem. Phys.

124

,

114107

(

2006

).

https://doi.org/10.1063/1.2176603

Google Scholar

Crossref

PubMed

85.

J. A.

Parkhill

and

M.

Head-Gordon

,

J. Chem. Phys.

133

,

124102

(

2010

).

https://doi.org/10.1063/1.3483556

Google Scholar

Crossref

PubMed

86.

T.

Zoboki

,

Á.

Szabados

, and

P. R.

Surján

,

J. Chem. Theory Comput.

9

,

2602

(

2013

).

https://doi.org/10.1021/ct400138m

Google Scholar

Crossref

PubMed

87.

E.

Xu

and

S.

Li

,

J. Chem. Phys.

139

,

174111

(

2013

).

https://doi.org/10.1063/1.4828739

Google Scholar

Crossref

PubMed

88.

F. E.

Harris

,

B.

Jeziorski

, and

H. J.

Monkhorst

,

Phys. Rev. A

23

,

1632

(

1981

).

https://doi.org/10.1103/physreva.23.1632

Google Scholar

Crossref

89.

T.

Arai

,

J. Chem. Phys.

33

,

95

(

1960

).

https://doi.org/10.1063/1.1731142

Google Scholar

Crossref

90.

In Ref. 64, it was supposed falsely that Λ_K would be zero for K = 3, … for APSG. It is only the γ_K contribution of nonzero elements that vanishes beyond Λ₂, according to Eq. (16).

91.

K. R.

Shamasundar

,

J. Chem. Phys.

131

,

174109

(

2009

).

https://doi.org/10.1063/1.3256237

Google Scholar

Crossref

PubMed

92.

This holds true as far as 4 < K ≤ N since geminal disconnected cumulants appear for 4 < N < K.

93.

R.

Maitra

,

D.

Sinha

,

S.

Sen

, and

D.

Mukherjee

,

Theor. Chem. Acc.

133

,

1522

(

2014

).

https://doi.org/10.1007/s00214-014-1522-5

Google Scholar

Crossref

94.

P.

Jeszenszki

,

V.

Rassolov

,

P. R.

Surján

, and

Á.

Szabados

,

Mol. Phys.

113

,

249

(

2015

).

https://doi.org/10.1080/00268976.2014.936919

Google Scholar

Crossref

95.

V. A.

Rassolov

and

F.

Xu

,

J. Chem. Phys.

126

,

234112

(

2007

).

https://doi.org/10.1063/1.2741266

Google Scholar

Crossref

PubMed

96.

K. V.

Lawler

,

D. W.

Small

, and

M.

Head-Gordon

,

J. Phys. Chem. A

114

,

2930

(

2010

).

https://doi.org/10.1021/jp911009f

Google Scholar

Crossref

PubMed

97.

J. M.

Foster

and

S. F.

Boys

,

Rev. Mod. Phys.

32

,

300

(

1960

).

https://doi.org/10.1103/revmodphys.32.300

Google Scholar

Crossref

98.

K.

Chatterjee

,

E.

Pastorczak

,

K.

Jawulski

, and

K.

Pernal

,

J. Chem. Phys.

144

,

244111

(

2016

).

https://doi.org/10.1063/1.4954694

Google Scholar

Crossref

PubMed

99.

Á.

Margócsy

,

P.

Kowalski

,

K.

Pernal

, and

Á.

Szabados

,

Theor. Chem. Acc.

137

,

159

(

2018

).

https://doi.org/10.1007/s00214-018-2355-4

Google Scholar

Crossref

100.

A. A.

Holmes

,

H. J.

Changlani

, and

C. J.

Umrigar

,

J. Chem. Theory Comput.

12

,

1561

(

2016

).

https://doi.org/10.1021/acs.jctc.5b01170

Google Scholar

Crossref

PubMed

101.

J.

Li

,

M.

Otten

,

A. A.

Holmes

,

S.

Sharma

, and

C. J.

Umrigar

,

J. Chem. Phys.

149

,

214110

(

2018

).

https://doi.org/10.1063/1.5055390

Google Scholar

Crossref

PubMed

102.

J.

Olsen

,

B. O.

Roos

,

P.

Jørgensen

, and

H. J. A.

Jensen

,

J. Chem. Phys.

89

,

2185

(

1988

).

https://doi.org/10.1063/1.455063

Google Scholar

Crossref

103.

M. J.

Frisch

,

G. W.

Trucks

,

H. B.

Schlegel

,

G. E.

Scuseria

,

M. A.

Robb

,

J. R.

Cheeseman

,

G.

Scalmani

,

V.

Barone

,

B.

Mennucci

,

G. A.

Petersson

,

H.

Nakatsuji

,

M.

Caricato

,

X.

Li

,

H. P.

Hratchian

,

A. F.

Izmaylov

,

J.

Bloino

,

G.

Zheng

,

J. L.

Sonnenberg

,

M.

Hada

,

M.

Ehara

,

K.

Toyota

,

R.

Fukuda

,

J.

Hasegawa

,

M.

Ishida

,

T.

Nakajima

,

Y.

Honda

,

O.

Kitao

,

H.

Nakai

,

T.

Vreven

,

J. A.

Montgomery

, Jr.,

J. E.

Peralta

,

F.

Ogliaro

,

M.

Bearpark

,

J. J.

Heyd

,

E.

Brothers

,

K. N.

Kudin

,

V. N.

Staroverov

,

R.

Kobayashi

,

J.

Normand

,

K.

Raghavachari

,

A.

Rendell

,

J. C.

Burant

,

S. S.

Iyengar

,

J.

Tomasi

,

M.

Cossi

,

N.

Rega

,

J. M.

Millam

,

M.

Klene

,

J. E.

Knox

,

J. B.

Cross

,

V.

Bakken

,

C.

Adamo

,

J.

Jaramillo

,

R.

Gomperts

,

R. E.

Stratmann

,

O.

Yazyev

,

A. J.

Austin

,

R.

Cammi

,

C.

Pomelli

,

J. W.

Ochterski

,

R. L.

Martin

,

K.

Morokuma

,

V. G.

Zakrzewski

,

G. A.

Voth

,

P.

Salvador

,

J. J.

Dannenberg

,

S.

Dapprich

,

A. D.

Daniels

,

Ö.

Farkas

,

J. B.

Foresman

,

J. V.

Ortiz

,

J.

Cioslowski

, and

D. J.

Fox

, Gaussian 09 Revision C.1,

Gaussian, Inc.

,

Wallingford, CT

,

2009

.

104.

P. R.

Surján

,

P.

Jeszenszki

, and

Á.

Szabados

,

Mol. Phys.

113

,

2960

(

2015

).

https://doi.org/10.1080/00268976.2015.1060366

Google Scholar

Crossref

105.

G. D.

Purvis

III,

R.

Shepard

,

F. B.

Brown

, and

R. J.

Bartlett

,

Int. J. Quantum Chem.

23

,

835

(

1983

).

https://doi.org/10.1002/qua.560230307

Google Scholar

Crossref

106.

X.

Ren

,

P.

Rinke

,

C.

Joas

, and

M.

Scheffler

,

J. Mater. Sci.

47

,

7447

(

2012

).

https://doi.org/10.1007/s10853-012-6570-4

Google Scholar

Crossref

107.

G. E.

Scuseria

,

T. M.

Henderson

, and

I. W.

Bulik

,

J. Chem. Phys.

139

,

104113

(

2013

).

https://doi.org/10.1063/1.4820557

Google Scholar

Crossref

PubMed

108.

T. H.

Dunning

, Jr.,

J. Chem. Phys.

53

,

2823

(

1970

).

https://doi.org/10.1063/1.1674408

Google Scholar

Crossref

109.

W.

Klopper

,

A. M.

Teale

,

S.

Coriani

,

T. B.

Pedersen

, and

T.

Helgaker

,

Chem. Phys. Lett.

510

,

147

(

2011

).

https://doi.org/10.1016/j.cplett.2011.04.101

Google Scholar

Crossref

110.

B.

Mussard

,

P.

Reinhardt

,

J. G.

Ángyán

, and

J.

Toulouse

,

J. Chem. Phys.

142

,

154123

(

2015

).

https://doi.org/10.1063/1.4918710

Google Scholar

Crossref

PubMed

111.

I. A.

Zaporozhets

,

V. V.

Ivanov

,

D. I.

Lyakh

, and

L.

Adamowicz

,

J. Chem. Phys.

143

,

024109

(

2015

).

https://doi.org/10.1063/1.4926392

Google Scholar

Crossref

PubMed

112.

F. A.

Evangelista

,

J. Chem. Phys.

134

,

224102

(

2011

).

https://doi.org/10.1063/1.3598471

Google Scholar

Crossref

PubMed

113.

T. H.

Dunning

,

J. Chem. Phys.

90

,

1007

(

1989

).

https://doi.org/10.1063/1.456153

Google Scholar

Crossref

114.

P. R.

Surján

,

I.

Mayer

, and

M.

Kertész

,

J. Chem. Phys.

77

,

2454

(

1982

).

https://doi.org/10.1063/1.444115

Google Scholar

Crossref

115.

Z. É.

Mihálka

,

P. R.

Surján

, and

Á.

Szabados

,

J. Chem. Theory Comput.

16

,

892

(

2020

).

https://doi.org/10.1021/acs.jctc.9b00858

Google Scholar

Crossref

PubMed

116.

P.-O.

Löwdin

,

Adv. Phys.

5

,

1

(

1956

).

https://doi.org/10.1080/00018735600101155

Google Scholar

Crossref

2020

Author(s)

Ring coupled cluster doubles at the multireference level

I. INTRODUCTION

II. THEORY

A. Notations

B. Ring coupled cluster doubles equation

C. Cumulants

D. Pruning the excitation manifold

E. Redundancy handling

F. Size extensivity

III. NUMERICAL ILLUSTRATIONS

A. Covalent bond breaking and formation

1. H₂O molecule

2. Bond dissociation of N₂

3. H₄ system

4. BeH₂ system

B. Torsional motion of small molecules

C. Deprotonation energies

IV. CONCLUSION

DATA AVAILABILITY

ACKNOWLEDGMENTS

APPENDIX A: TENSORS OF THE SPIN-FREE rCCD EQUATIONS

APPENDIX B: CUMULANTS OF THE APSG WAVEFUNCTION

APPENDIX C: DETAILS OF REDUNDANCY HANDLING

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Ring coupled cluster doubles at the multireference level Open Access

I. INTRODUCTION

II. THEORY

A. Notations

B. Ring coupled cluster doubles equation

C. Cumulants

D. Pruning the excitation manifold

E. Redundancy handling

F. Size extensivity

III. NUMERICAL ILLUSTRATIONS

A. Covalent bond breaking and formation

1. H2O molecule

2. Bond dissociation of N2

3. H4 system

4. BeH2 system

B. Torsional motion of small molecules

C. Deprotonation energies

IV. CONCLUSION

DATA AVAILABILITY

ACKNOWLEDGMENTS

APPENDIX A: TENSORS OF THE SPIN-FREE rCCD EQUATIONS

APPENDIX B: CUMULANTS OF THE APSG WAVEFUNCTION

APPENDIX C: DETAILS OF REDUNDANCY HANDLING

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

Ring coupled cluster doubles at the multireference level

1. H₂O molecule

2. Bond dissociation of N₂

3. H₄ system

4. BeH₂ system