Restricted single-reference coupled cluster theory truncated to single and double excitations accurately describes weakly correlated systems, but often breaks down in the presence of static or strong correlation. Good coupled cluster energies in the presence of degeneracies can be obtained by using a symmetry-broken reference, such as unrestricted Hartree-Fock, but at the cost of good quantum numbers. A large body of work has shown that modifying the coupled cluster *ansatz* allows for the treatment of strong correlation within a single-reference, symmetry-adapted framework. The recently introduced singlet-paired coupled cluster doubles (CCD0) method is one such model, which recovers correct behavior for strong correlation without requiring symmetry breaking in the reference. Here, we extend singlet-paired coupled cluster for application to open shells via restricted open-shell singlet-paired coupled cluster singles and doubles (ROCCSD0). The ROCCSD0 approach retains the benefits of standard coupled cluster theory and recovers correct behavior for strongly correlated, open-shell systems using a spin-preserving ROHF reference.

## I. INTRODUCTION

The electronic structure of weakly correlated systems is accurately described by the coupled cluster (CC) family of methods.^{1–4} Coupled cluster with low order excitations, such as CC with single, double, and perturbative triple excitations, CCSD(T), gives quantitatively correct results for weakly correlated systems, and the method is routinely applied to a wide range of quantum chemical problems. However, as a single-reference theory, coupled cluster usually fails to describe strongly correlated systems, where the reference wavefunction, e.g., restricted Hartree-Fock (RHF), fails to even qualitatively describe the physics in the many-body spectrum, usually reflected in the single-particle spectrum. Such multi-reference problems are characterized by degeneracies, i.e., for symmetry-adapted RHF orbitals, the HOMO-LUMO gap vanishes.

One solution is to allow symmetry breaking in the reference, e.g., using an unrestricted Hartree-Fock wavefunction. Symmetry breaking lifts the degeneracy, allowing the HOMO and LUMO to become gapped. Unrestricted coupled cluster usually has well-behaved energies, but can suffer greatly from spin contamination. Thus, though the energies may be good, properties based on the contaminated wavefunction may be unreliable. Ideally, we might solve the problem exactly in the one-particle basis via full configuration interaction (FCI), but this approach is rarely computationally feasible. Another approach is to use a multi-reference method, such as multi-reference coupled cluster, and significant progress has been made along these lines.^{5} However, it would be interesting to have a simpler alternative more in the spirit of single-reference coupled cluster.

Recent work has shown that it might be possible to describe strong correlation with a judicious modification of single-reference coupled cluster theory or its variants.^{6–19} It seems that the failure of coupled cluster to describe strong correlation can be attributed not merely to inadequacies of the reference but also to instabilities in the truncated coupled cluster *ansatz*. Eliminating or decoupling pairing channels of the *T*_{2} operator enables restricted coupled cluster to recover qualitatively correct behavior in describing strongly correlated systems.^{20–27}

With some notable exceptions,^{17} most work using a simplification of single-reference coupled cluster theory to recover qualitatively correct behavior for strong correlation has been done only for closed-shell systems. In this manuscript, we extend the search for an inexpensive method for strong correlation in the open-shell regime via restricted open-shell singlet-paired coupled cluster theory. For completeness, we first review singlet-paired coupled cluster theory for closed shells^{27} before generalizing the method for open shells. We then present some interesting benchmark calculations before offering concluding remarks.

## II. THEORY AND METHODS

### A. Closed-shell singlet-paired coupled cluster

In closed-shell coupled cluster with double excitations (CCD), we write the wavefunction as

where |0〉 is the restricted Hartree-Fock (RHF) reference. The cluster operator *T*_{2} is spin adapted^{28} and creates double excitations,

Here, orbitals *i j* (*a b*) are occupied (unoccupied) in the RHF reference, spin indices *σ* and *ξ* are summed over both up and down spin, and summation over repeated indices is implied. The amplitudes $tijab$ are found by solving

where $|ijab\u3009$ is notation for referring to the space of symmetry-adapted, doubly excited determinants.^{28} The energy is then found by evaluating

using the amplitudes from Eq. (3).

With the addition of single excitations, CCSD accurately describes weakly correlated systems, but breaks down when systems take on multireference character, as evidenced by the familiar unphysical hump in the dissociation of N_{2} in the STO-3G basis,^{29} shown in Fig. 1. Although CCSD is a good approximation of the full configuration interaction result near equilibrium, this curve erroneously predicts molecular N_{2} to be only slightly stabler than the atomic radicals and gives energies at dissociation below the exact FCI result.

It has been known for a long time that *T*_{2} can be written in terms of a particle-particle/hole-hole recoupling.^{30,31} This led Bulik *et al.*^{27} to notice that Eq. (2) can be recast as

where $T2[0]$ and $T2[1]$ are the singlet- and triplet-paired components of the full *T*_{2}, respectively, and are given by

The pair operators are given by

where

For a closed-shell reference |0〉, $Pab\u2020|0\u3009$ and *P _{ij}*|0〉 are

*N*+ 2 and

*N*− 2 electron singlets, respectively, where

*N*is the number of electrons in the reference. We also note that $Pab\u2020$ is given by

where we first substitute *i* → *a*, *j* → *b* into Eq. (8) and then transpose.

Since $Pab\u2020$ is symmetric on interchange of indices *a* and *b*, so too are the $\sigma ijab$ amplitudes. It can be shown that $\sigma ijab$ is related to the full $tijab$ by

Singlet-paired coupled cluster has a wavefunction *ansatz* from which coupled cluster equations follow. That is, one can derive fully symmetric amplitude and energy equations for singlet-paired coupled cluster doubles (CCD0) by writing the CCD0 wavefunction as

and inserting $T2[0]$ instead of the full *T*_{2} into the working equations for standard coupled cluster, given by Eqs. (3) and (4). However, rather than re-deriving the necessary equations, one can quickly obtain the desired result from a standard CCD code by symmetrizing the $tijab$ amplitudes according to Eq. (12) at the end of each coupled cluster iteration.

It is worth mentioning that the $\sigma ijab$ amplitudes obtained from CCD0 are different from the symmetric part of the $tijab$ amplitudes obtained from standard CCD. Equation (5) writes *T*_{2} in terms of singlet- and triplet-paired components of the operator. For the amplitudes, we also have

where $\sigma ijab$ and $\pi ijab$ are the singlet- and triplet-paired amplitudes, respectively. As $\sigma ijab$ is the symmetric part of $tijab$, $\pi ijab$ is the antisymmetric part, given by

Singlet-paired coupled cluster decouples the symmetric part of the amplitudes from the antisymmetric part by setting $\pi ijab$ in Eq. (14) equal to zero. In standard coupled cluster, the singlet- and triplet-pairing channels are fully coupled. Including other orders of excitations with $T2[0]$, e.g., CCSDT0, maintains correct physics while improving on the accuracy of the method.^{27} As can be seen in Fig. 1, although singlet-paired coupled cluster does not give exact energies, it gives the qualitatively correct shape of the dissociation curve.

Before moving on to the generalization of CCD0 for open shells, we note that one can use $T2[1]$ rather than $T2[0]$ in the above procedure to generate triplet-paired coupled cluster doubles (CCD1) and solve for the antisymmetric $\pi ijab$ amplitudes, rather than finding the symmetric $\sigma ijab$ amplitudes via CCD0.^{27} In the present work, we focus on the open-shell extension of singlet-paired coupled cluster. The electron correlation comes primarily from opposite spin interactions, which CCD0 describes more fully than CCD1. Triplet-paired coupled cluster is interesting in that it recovers a qualitatively correct description of strong correlation, but because its energies are generally substantially worse than those given by CCD0, we do not discuss it further here.

### B. Open-shell singlet-paired coupled cluster

Extending the singlet-pairing scheme discussed above for use in describing open-shell systems poses three primary challenges.

First, notice that the singlet-paired annihilation operator in Eq. (8) is only singlet-paired when the indices *i* and *j* run over the same set of spatial orbitals. Otherwise, *P _{ij}* is not antisymmetric on interchange of the spin indices, and neither $Pab\u2020|0\u3009$ nor

*P*|0〉 is a singlet wavefunction. Thus, the reference orbitals used to define the pairing scheme cannot come from unrestricted Hartree-Fock (UHF), in which the

_{ij}*α*and

*β*spatial orbitals are allowed to differ. Even the semicanonical orbitals from restricted open-shell Hartree Fock (ROHF) calculations will not suffice, as these only guarantee that

*α*and

*β*orbitals span the same space. We can solve this problem by using as our reference the so-called Natural Orbitals (NO), i.e., ROHF orbitals in the basis that diagonalizes the charge density matrix,

^{32}

where *γ*^{α} and *γ*^{β} are the alpha and beta density matrices, respectively. The transformation to the NO Basis is defined by unitary *U* such that

where *n* is diagonal, with eigenvalues that are occupation numbers of 1, $12$, and 0. In this non-canonical basis, the occupied-occupied and virtual-virtual blocks of the Fock matrices are no longer diagonal. In our code, we use conjugate gradients^{33} to solve the non-canonical coupled cluster equations at each iteration.

Second, the *T*_{2} operator becomes more complicated in the open-shell case. While for RHF-based calculations, *T*_{2} takes the form shown in Eq. (2), ROHF-based calculations on open shells introduce other classes of excitations, namely, excitations into and out of the singly occupied space.^{34} Such excitations do not appear in the derivation of RHF-based CCD0, and it is at first glance unclear what role they play in the extension of CCD0 to open shells. We might be tempted to do blocked-CCD0 on open shells, wherein we only allow excitations from the doubly occupied to virtual orbitals in the ROHF reference. In other words, we might freeze the open-shell orbitals, and they would not participate in the excitations as either originators or recipients. Such a procedure eliminates the need to deal with the open-shell excitations, but is clearly unsatisfactory, as these open-shell excitations are needed for an accurate description of the system. We have found that these open-shell excitations seem to play little role in the breakdown of restricted open-shell coupled cluster in the systems studied here. As such, we fully include these amplitudes while taking only the singlet-paired component of the doubly occupied to virtual (*oo-vv*) amplitude block. We refer to this method as restricted open-shell singlet-paired coupled cluster doubles (ROCCD0). In practice, the coupled-cluster description of open shells is greatly improved by the inclusion of single excitations. In singlet-paired coupled cluster for open shells, we always include single excitations as in RHF-based CCSD0 to give ROCCSD0. To recapitulate, for ROCCSD0 we take the symmetric part of $toovv$, the full amplitudes in the other blocks of *T*_{2}, and full single excitations.

Third, we write the coupled cluster equations schematically as

where *Ht* is shorthand for the contractions between the Fock matrix and amplitudes,

For RHF-based CCD0, this equation takes the form

where we solve for the symmetric *σ* amplitudes from a symmetric right hand side *G*_{sym}. In addition, for the closed-shell case, *H* is also symmetric because *f*^{α} = *f*^{β}.

For the open-shell case, however, *H* is not necessarily symmetric, because *f*^{α}≠*f*^{β}. Spin-summing the coupled cluster equations gives three sets of equations: two same-spin equations for alpha/alpha and beta/beta excitations and the alpha/beta equations describing opposite-spin excitations. For the alpha/beta equations, Eq. (19) becomes

where an overbar indicates a beta spin index, a lack of an overbar indicates an alpha spin index, and the indices run over the occupied/virtual dimensions appropriate for that spin. Because we now have *f*^{α}≠*f*^{β}, Eq. (21) is not invariant to interchange of the spin indices. For ROCCSD0, then, we must also take the symmetric part of *H* in the *oo-vv* block in addition to the symmetric part of *T* and *G*,

We construct *H*_{sym} by replacing the alpha and beta Fock matrices contracting over the doubly occupied to virtual block in Eq. (21) with Fock charge, defined in analogy with the charge density matrix,

Amplitudes outside of the doubly occupied to virtual block are contracted as usual over their respective Fock matrices, and the *oo-vv* block of Eq. (21) becomes

where now *I*, *J* run over doubly occupied orbitals, *A*, *B* run over virtual orbitals, and *X*, *Y* run over singly occupied orbitals. We emphasize that, while the Fock matrices are zero in the *oo-vv* block, the singly occupied to virtual blocks of the Fock matrices are nonzero, resulting in the last two terms in Eq. (24). We then move the last two terms in Eq. (24) to the right hand side of Eq. (22) before symmetrizing *G*.

ROCCSD0 can be implemented as follows:

Construct the ROHF reference.

Transform the integrals to the NO basis.

Iterate:

Build standard

*G*.Move excitations involving singly occupied orbitals that couple to the

*oo-vv*block over to*G*. From Eq. (24): $GIJAB\u2190GIJAB\u2212fIXtXJ\xafAB\xaf+fY\xafB\xaftIJ\xafAY\xaf$.Symmetrize the doubly occupied to doubly virtual block of

*G*: $GIJAB\u219012(GIJAB+GIJBA)$.Solve

*H*_{sym}*σ*=*G*_{sym}(*σ*) in the doubly occupied to virtual block.Solve standard

*Ht*=*G*(*t*) in all other amplitude blocks, including singles.

### C. Computational details

All Hartree-Fock and standard coupled cluster calculations were done in *Gaussian 09*,^{35} while the non-canonical ROCCSD0 calculations were performed using in-house code. Our in-house code does not support point group symmetry, so our ROHF calculations were done without point group symmetry, and we used the lowest-energy orbitals we could find as our coupled-cluster references. Where we found multiple smooth ROHF curves, we have shown ROCCSD and ROCCSD0 results on both references. Throughout, we use cartesian *d* functions and work in relatively small basis sets.^{36,37} Because dynamic correlation is basis-set dependent, while static correlation is not, static correlation dominates in small bases. Working in such basis sets allows us to emphasize the deficiencies of standard coupled cluster to describe strong correlation.

## III. RESULTS

First, we study the behavior of ROCCSD0 on the dissociation of a set of linear, equally spaced hydrogen chains with open boundary conditions. The top panel of Fig. 2 shows the dissociation of a doublet 5-hydrogen chain in the STO-3G basis set. For ROCCSD, we observe a breakdown similar to that seen in the CCSD dissociation of N_{2}. The ROCCSD curve turns over at around 3.5 bohrs before becoming difficult to converge. ROCCSD0 on the other hand is both well behaved and convergent at longer bond lengths. In addition, ROCCSD0 gives energies across the range of bond lengths nearly indistinguishable from UCCSD(T), which gives the right energy for this system, but the wrong wavefunction because of spin contamination. Looking at the bottom panel of Fig. 2, we see the dissociation of the same system in the larger cc-pVDZ basis set.^{38} The larger basis makes ROCCSD generally better behaved, pushing the point of breakdown out to around 4 bohrs. However, ROCCSD still breaks down as catastrophically in this basis as in the minimal basis case. Once again, ROCCSD0 recovers the correct behavior, is indistinguishable from UCCSD(T) at dissociation, and sacrifices only a fraction of the correlation energy at equilibrium.

Figure 3 shows the dissociation of an equally spaced 7-hydrogen chain in the cc-pVDZ basis. Here, the breakdown of ROCCSD becomes more pronounced relative to that observed for the 5-hydrogen case; the ROCCSD curve turns over sooner, the divergence is much steeper, and the equations become harder to converge at shorter bond lengths. Nevertheless, we see the same recovery of correct physics by ROCCSD0. At dissociation, ROCCSD0 is very near UCCSD(T), despite the significant breakdown of ROCCSD.

We now apply ROCCSD0 to some simple triplet hydrogen chains. In Fig. 4, we see the dissociation of an equally spaced triplet 6-hydrogen chain in the cc-pVDZ basis. Relative to the 5-hydrogen chain shown above, this system is evidently a tougher problem in that ROCCSD breaks down and stops converging sooner. In addition, although ROCCSD0 does recover the correct description of the dissociation, it also stops converging sooner. We also see a more notable difference between the energies predicted by UCCSD(T) and ROCCSD0, a result we would generally expect.

Figure 5 shows the dissociation of a triplet 8-hydrogen chain in cc-pVDZ. Once again, the larger system is more difficult to describe, as ROCCSD turns over much more sharply than for the 6-hydrogen case and stops converging more quickly. In contrast, ROCCSD0 can be converged nearly as far for the 8-hydrogen chain as for the 6-hydrogen case. We again see the expected recovery of correct dissociation behavior by ROCCSD0 and note that we now have a noticeable difference between the ROCCSD0 and UCCSD(T) energies.

While ROCCSD0 correctly describes these systems using reference determinants that are eigenfunctions of $S2$, UCCSD(T) requires severely symmetry-broken UHF references to accomplish the same feat. While $S2$ for a doublet is 0.75, at dissociation, UHF in cc-pVDZ predicts $S2$ values for H_{5} and H_{7} of 2.84 and 3.73, respectively. Similarly, for a triplet, we should have $S2=2$. However, at 6 bohrs, UHF in cc-pVDZ predicts $S2$ for triplet H_{6} and triplet H_{8} both to be around 3.99. These numbers emphasize the need to preserve wavefunction symmetries if possible. Though the energies in such cases may be good, symmetry-dependent molecular properties will be poor. Another question, however, is how spin contaminated the coupled cluster wavefunctions are. If ROCCSD is spin adapted, i.e., the expectation value of *S*^{2} is correct, so is ROCCSD0. However, both methods retain terms whose recoupling can yield non-singlets in open-shell cases. A precise quantification of this effect would require implementation of ROCCSD0 response density matrices, which we will carry out in due course.

We see slightly different behavior in the dissociation of a triplet 8-hydrogen ring in Fig. 6. Here, two ROHF curves result in pathological behavior of the coupled cluster energies. The two SCF energy surfaces and their corresponding coupled cluster solutions are nearly identical until around 4.3 bohrs, at which point ROHF-A begins to dissociate to a lower limit than ROHF-B. It is interesting that this point is also near the point where ROCCSD begins to break down. Although ROCCSD-A turns over less sharply than ROCCSD-B, both curves are disastrous. ROCCSD0-B does not completely recover correct behavior, but is much improved over ROCCSD-B. ROCCSD0-A, on the other hand, is well behaved. One of the difficulties of restricted open-shell coupled cluster is the choice of ROHF reference. In this case, the choice is clear, as the lower energy ROHF curve results in the better ROCCSD0 curve, and ROCCSD0 on this ROHF-A reference gives excellent results.

We now turn our attention to a slightly more complicated problem: the dissociation of triplet N_{2}. The breakdown of RHF-based coupled cluster in the description of the dissociation of singlet N_{2} is well known. As is apparent in Fig. 7, this breakdown also extends to the ROCCSD description of the dissociation of triplet N_{2} for the solution curves labeled “A.” As in the case of the 8-hydrogen systems, we find at least two smooth ROHF curves here. For triplet N_{2} however, the choice of which ROHF reference to use for ROCCSD is muddled because, while ROHF-A is lower at equilibrium, ROHF-B goes to a lower limit at dissociation. In the top panel of Fig. 7, we plot the energy of triplet N_{2} at each bond length relative to the quartet N and doublet N atomic energies from the same method. Here, the breakdown of ROCCSD-A is exacerbated relative to the singlet case, which at least predicted molecular N_{2} to be more stable than the atomic ions even in the STO-3G basis. Even in the larger cc-pVDZ basis^{38} calculations shown in Fig. 7, ROCCSD-A unphysically turns over to eventually predict triplet N_{2} to be unbound. Although ROCCSD0-A predicts too deep of a minimum at equilibrium, it goes to nearly the correct limit of the sum of the ROCCSD0 quartet N and doublet N radical energies at dissociation. Coupled cluster on the ROHF-B reference is well behaved, dissociating to the correct limit. Since there is no breakdown in this case to correct, ROCCSD0-B simplify shifts the energy up from ROCCSD-B across the range of bondlengths. Although ROCCSD0-B dissociates to too high of a limit, the depth of the energy minimum is comparable to that of the ROCCSD-B curve.

In the bottom panel of Fig. 7, we plot the total energies for the same system. This plot makes it easier to see how ROCCSD0 recovers correct behavior from ROCCSD, but also highlights an important feature of ROCCSD0 not as obvious as in the applications on hydrogen chains. That is, singlet-paired coupled cluster generally sacrifices a significant portion of the weak correlation in order to protect the method from failure due to strong correlation. This observation is in line with results from closed-shell CCD0,^{27} and work to remedy this deficiency is on going.^{39,40}

## IV. DISCUSSION

Bulik *et al.*^{27} showed that decoupling the singlet- and triplet-pairing channels of the *T*_{2} operator substantially ameliorates the failure of RHF-based coupled cluster to describe strong correlation. Here, we have extended their formalism to open shells via ROCCSD0, which, like RHF-based CCD0, maintains the advantages of traditional coupled cluster while recovering correct physics for describing strong correlation. It is interesting that for the systems studied here, the breakdown of ROCCSD can be attributed to the same coupling between the singlet-paired and triplet-paired channels of *T*_{2} that is apparently the culprit for the breakdown of RHF-based CCD. We are able to fully include excitations involving singly occupied orbitals in ROCCSD0, which are only indirectly affected by being coupled to the $\sigma ijab$ amplitudes in the doubly occupied to virtual block rather than the full $tijab$.

For the applications to hydrogen chains, ROCCSD0 produces curves comparable to unrestricted coupled cluster, which required massive spin contamination in the reference, while correlating a spin-preserving ROHF wavefunction. For the triplet 8-hydrogen ring, ROCCSD0 is greatly improved relative to ROCCSD, but this success appears to depend somewhat on the choice of ROHF reference. For triplet N_{2}, ROCCSD0 recovers correct behavior again when ROCCSD breaks down, dissociating to nearly the correct limit, i.e., the sum of the ROCCSD0 energies for the individual atoms, but the loss of weak correlation is more apparent. One ROHF reference for this system resulted in good ROCCSD energies. ROCCSD0 on this reference maintains correct behavior, but sacrifices some correlation. To highlight the effectiveness of ROCCSD0, we have presented examples where standard ROCCSD fails dramatically at long bond lengths, but it is worth mentioning that ROCCSD is a robust method that often does not fail as obviously for strongly correlated systems as its closed-shell counter part. We speculate that if the breakdown of ROCCSD is associated with quadratic terms in the coupled cluster equations connecting doubly occupied to virtual orbitals, then the presence of singly occupied orbitals protects ROCCSD from breakdown by widening the effective gap between the highest doubly occupied molecular orbital and lowest virtual orbital.

A useful feature of RHF-based CCD0 is that it completely eliminates explicit same-spin correlation, making it possible to combine CCD0 correlation energies with same-spin correlation from homogenous electron gas-based density functionals and the random phase approximation to obtain more quantitative results.^{39,40} This combination is not straightforward without double counting in our current formulation of ROCCSD0, because we retain same-spin excitations involving open-shell orbitals. However, since it is apparently the doubly occupied to virtual block of *T*_{2} that must be decoupled to recover correct behavior for strong correlation, there is some flexibility in how we treat the open-shell excitations. Since ROCCSD0 already removes same-spin excitations in the doubly occupied to virtual block, it might be fruitful to eliminate these excitations in the open-shell block as well, enabling the combination of ROCCSD0 with same-spin correlation from DFT. Such an ROCCSD0+DFT method might prove to be a useful tool in the study of systems such as open-shell actinide complexes, as RHF-based CCD0+DFT seems to be,^{39,40} but this is not an avenue we explore here. We are currently working on how to recouple the singlet- and triplet-pairing channels of *T*_{2} in order to obtain more accurate results for closed- and open-shell systems without reintroducing the failures of standard coupled cluster, but such a procedure is also not straightforward.

We conclude by offering a general remark on the apparent success of singlet-paired coupled cluster over traditional coupled cluster. What this work and previous work have shown is that it is not simply deficiencies in the single-determinant reference that cause correlated methods such as coupled cluster theory to fail. This is something we already knew, of course, since full configuration interaction built upon a single-reference wavefunction can provide the exact solution in the space of single-particle basis functions, even in the presence of strong correlation. Full configuration interaction tells us that it is possible to describe strong correlation at the single-reference level. Singlet-paired coupled cluster tells us that it may be possible to do so in a computationally tractable manner, provided that we have the right wavefunction *ansatz*, and it is our hope that this work, though not the final goal, will be a step in the right direction.

## Acknowledgments

This work was supported in part by the U.S. Department of Energy Heavy Element Chemistry Program (No. DE-FG02-04ER15523). G.E.S. is a Welch Foundation Chair (C-0036). J.A.G. gratefully acknowledges the support of the National Science Foundation Graduate Research Fellowship Program (No. DGE-1450681).