When calculating the spin multiplicity at either the second-order Møller-Plesset (MP2) or the iterative second-order approximate coupled-cluster singles and doubles (CC2) levels of theory using the same strategy for the calculation of the expectation value as in regular CC theory together with the usual definitions of the MP2 and CC2 density matrices, artificial spin contamination occurs in closed-shell molecules. Non-intuitively, for open-shell systems, results at the MP2 or CC2 levels of theory based on this procedure even suggest stronger contamination at the correlated level than for the Hartree–Fock reference, although treatment of electron correlation should lower spin contamination. In this Communication, the reasons behind this inconsistency are investigated and a solution is proposed, which removes spin contamination for closed-shell molecules and leads to physically meaningful results for open-shell cases. Additionally, we show that CC2 significantly outperforms MP2 in describing systems with a strongly spin-contaminated reference with a performance similar to that of full coupled-cluster with singles and doubles substitutions (CCSD).

## I. INTRODUCTION

Coupled-Cluster (CC) theory is one of the most accurate and widely used methods for *ab initio* calculations. CCSD(T), i.e., CC with singles and doubles substitutions together with a perturbative treatment of triple excitations, is used today as a gold standard.^{1} In standard CC theory, a Hartree–Fock (HF) determinant is used as the reference wavefunction. For open-shell systems, a CC treatment based on either a restricted open-shell HF (ROHF) or an unrestricted HF (UHF) reference, respectively, leads to spin-contaminated wavefunctions.^{2} In order to obtain a sound wavefunction, spin contamination needs to be monitored using appropriate tools. Usually, the treatment of electron correlation, in particular when using CC theory, significantly lowers the extent of spin contamination as compared to the HF reference.^{3–6}

The first tool to be developed to calculate the spin multiplicity and hence to study the extent of spin contamination in CC wavefunctions was introduced by Purvis, Sekino, and Bartlett.^{3} In that work, the $\u015c2$ operator was used, with $\u015c$ being the total spin operator $\u015c=\u015cxi+\u015cyj+\u015czk$ and ** i**,

**, and**

*j***being the Cartesian unit vectors. For an eigenfunction to $\u015c2$, |**

*k**S*,

*M*

_{S}⟩, the eigenvalue is given in atomic units by

with *S* being the total spin quantum number, *M*_{S} being the spin projection quantum number, and 2*S* + 1 being the spin multiplicity. To investigate the extent of spin contamination, Purvis, Sekino, and Bartlett suggested, in analogy to the CC energy expression, to project the action of this operator on the CC wavefunction onto the reference wavefunction $\u30080|\u015c2|\Psi CC\u2009\u3009=\u27e8\u015c2\u27e9proj$.^{3}

It has been pointed out that the projected expression leads only to an approximation of the spin multiplicity of the CC wavefunction.^{3–5} For this reason, Stanton proposed the calculation of the spin multiplicity using derivative theory.^{4} In this formalism, first-order properties are calculated as derivatives of the CC Lagrangian,^{7} which is equivalent to the expectation value in bivariational theory.^{8,9} Based on results for open-shell systems, Stanton suggested that there is no significant advantage in using an ROHF reference over a UHF reference since the respective CCSD energies seem to be rather invariant to the choice of reference and also the calculated spin contamination based on an ROHF vs a UHF reference is similar. Additionally, even for systems with a strongly spin-contaminated reference, the powerful treatment of correlation and effective orbital relaxation at the CCSD level brings the CC spin-expectation value very close to the exact eigenvalue *S*_{e.v.}.^{4} Note, however, that CCSD(T) or excited-state wavefunctions at the Equation of Motion (EOM)-CCSD level are rather sensitive to the spin contamination of the reference.^{10,11}

For large systems, CCSD with its *N*^{6} scaling, where *N* is the number of electrons, is often computationally very expensive and second-order Møller–Plesset perturbation (MP2) theory or the iterative second-order approximate coupled-cluster singles and doubles (CC2) method may be used instead.^{12–16}

Because of the intrinsic relation between Many-Body Peturbation Theory (MBPT) and CC theory, CC models can be viewed as an infinite-order summation of selected excitations of the MBPT theory.^{17} Hence, the projection approach for the calculation of the $\u015c2$ expectation value $\u27e8\u015c2\u27e9proj$ has already been introduced in the work of Purvis, Sekino, and Bartlett and is also used in the work of Chen and Schlegel for MP2.^{3,5} An extension to CC2 theory is readily available.

However, applying the expectation-value approach to MP2 and CC2 requires caution. In this Communication, we show that the results obtained using this formulation together with the usual definitions of the one- and, in particular, the two-particle reduced density matrices (RDMs) introduce artificial spin contamination to closed-shell systems, which is obviously not meaningful. In addition, at the MP2 and CC2 levels of theory, the $\u27e8\u015c2\u27e9$ results obtained in this manner give unexpectedly strong contamination for open-shell systems. After discussing the origin of this inconsistency, we propose a modification to this approach in Sec. II. We then describe the implementation in Sec. III, apply the modified approach to our test systems, and present the findings in Sec. IV.

## II. THEORY

The total spin-squared operator $\u015c2$ can be expressed (in atomic units) as the following sum:

where $\u015c\xb1$ are the spin ladder operators and $\u015cz$ is the spin projector operator onto the *z* axis. In second quantization, these operators are given as

The creation (*p*^{†}) and annihilation operators (*p*) refer to *α* electrons and *β* electrons when carrying overbars, respectively. The letters *p*, *q*, *r*, *s* are used for the combined hole and particle space; *a*, *b*, *c*, *d* for the particle space; and *i*, *j*, *k*, *l* for the hole space. The symbols $n^\sigma $, with *σ* = *α*, *β*, are the electron number operators for *α* and *β* electrons, respectively, and $\Delta n^=n^\alpha \u2212n^\beta $. The *S*_{±} elements of the ladder operators correspond to orbital-overlap integrals between *α* and *β* orbitals with

The norm of these elements is given as

The $\u015c2$ expectation value for the reference determinant can then be expressed by

In the UHF and ROHF approaches, the reference wavefunction is a spin eigenfunction to $\u015cz$ as it preserves the number of *α* and *β* electrons,

but only the ROHF wavefunction is an eigenfunction to $\u015c2$. For this reason, for UHF based calculations, one often calculates the deviation of the expectation value $\u27e8\u015c2\u27e9$ from the eigenvalue *S*_{e.v.} evaluated with the exact spin quantum number *S*, i.e., *S*(*S* + 1), as a diagnostic tool,^{5,6}

### A. Expressions using projection techniques

According to the derivation of Purvis, Sekino, and Bartlett, the CC correlation contribution for the projected expression $\u27e8\u015cN2\u27e9proj$ is given by^{3}

with $tia$ being the CC single amplitudes and $tijab$ being the CC double amplitudes.

Note that the projected expression yields no spin contamination with an ROHF reference |Φ_{ROHF}⟩ due to

even though the correlated wavefunction is indeed spin contaminated.^{4} As discussed by Stanton and Purvis, Sekino, and Bartlett,^{3,4}$\u27e8\u015cN2\u27e9proj$ can be considered an approximation since ⟨0| is not the left-side CC eigenstate and it does not take the response of the CC amplitudes into account. Chen and Schlegel showed with their first-order approximation to the expectation-value approach that |Δ*S*_{CC}| < |Δ*S*_{proj}| since higher-order correlation contributions are present in $\u27e8\u015cN2\u27e9CC$.^{5} On practical consideration, $\u27e8\u015cN2\u27e9proj$ is readily available after a CC energy calculation, while $\u27e8\u015cN2\u27e9CC$ requires the additional calculation of the *λ* amplitudes as usual for properties within CC theory. Furthermore, it has been pointed out by Stanton that, since in the expectation-value approach for CC theory the ket wavefunction $|\Psi CC\u2009\u3009=eT^|0\u2009\u3009$ is not the Hermitian conjugate of the bra wavefunction $\u3008\Psi \u0303CC|=\u30080|(1+\Lambda ^)e\u2212T^$, there is no lower bound for $\u27e8\u015c2\u27e9CC$ and spin contamination may even be negative Δ*S*_{CC} < 0.^{4}

Applying the projected variant to MP2 and CC2 theory is straightforward. For the MP2 case, the $tijab$ amplitudes are replaced by the first-order response of the wavefunction $tijab(1)=\u2212\u27e8ab|ij\u27e9\u2212\u27e8ab|\u2009ji\u27e9\epsilon ijab$, with $\epsilon ijab=\epsilon a+\epsilon b\u2212\epsilon i\u2212\epsilon j$ as the difference of the orbital energies *ϵ*_{p}, while the $tia(1)$ are zero. This is included in the work of Purvis, Sekino, and Bartlett, which contains a detailed description of $\u27e8\u015cN2\u27e9proj$ for the MBPT series.^{3} For the CC2 case, and generally for the CC*n* series, the CC amplitudes are replaced by the respective CC*n* solutions.

### B. Expressions based on the (bivariational) expectation value

A description improving upon the projected—and hence only approximate—expressions discussed before was given by Stanton via the CC expectation value of $\u015c2$.^{4} In this formulation, the CC Lagrangian is defined as^{17}

with $\Lambda ^$ being the de-excitation operator containing the *λ* amplitudes $\Lambda ^=\lambda aii\u2020a+14\lambda abiji\u2020aj\u2020b+\cdots \u2009$, $T^$ being the excitation operator containing the CC amplitudes $T^=tiaa\u2020i+14tijaba\u2020ib\u2020j+\cdots \u2009$, and $\u0124$ being the Hamiltonian. The $\u015c2$ expectation value is calculated as a first derivative of the Lagrangian using the modified Hamiltonian

where $\u01240$ is the electronic Hamiltonian and *μ* is a scalar prefactor that can be varied. The latter is formally set to 0 for the calculation of the energy, and its mere purpose is to introduce a perturbative linear dependency of the Lagrangian on the $\u015c2$ operator. The derivative is then given as

### C. Derivatives and density matrices in MP*n* and CC*n*

In the general case, and also here for *μ* = 0, the Lagrangian for the correlation contribution can be brought to the form^{7,8,17}

where *f*_{pq} are the elements of the Fock matrix and *W*_{pqrs} = ⟨*pq*|*rs*⟩ − ⟨*pq*|*sr*⟩ are the elements of the two-electron interaction part of the Hamiltonian, i.e., the antisymmetrized two-electron integrals. This means that the normal-ordered one- and two-particle RDMs are defined such that they reproduce the correlation energy when inserted into Eq. (3). In the case of CC theory, they are given as^{17}

with the curly bracket notation {⋯} signifying normal-ordered products of quasi-particle operators with respect to the Fermi vacuum. Essentially, the RDMs are defined by all parts of the expectation value apart from the integral itself.

In many cases, a first-order property *A*, depending on a parameter *κ*, can be calculated via

which also reflects the equivalence of generalized derivative theory and linear response apart from orbital relaxation.^{14,18,19} MP2 properties and CC2 properties are reported in Refs. 18 and 14 using derivative theory and response theory,^{19} respectively. Within a Lagrangian formulation, we may define

and take the corresponding derivative with respect to *κ*. Here, $H\u0303=e\u2212T^\u0124eT^$ and $\u01241=e\u2212T^1\u0124eT^1$. In MP2 theory, the *λ* amplitudes are simply given by the complex conjugate of the first-order response amplitudes $\lambda abij(1)=(tijab(1))*=\u2212\u27e8ij|ab\u27e9\u2212\u27e8ij|ba\u27e9\epsilon ijab$. In CC2 theory, the *λ* amplitudes are obtained iteratively.^{14} For perturbative methods such as MP*n* and CC*n* that rely on the MP partitioning of the Hamiltonian, the Fock operator $F^$ is part of the unperturbed system, while the two-electron part $\u0174$ is the perturbation. Hence, in Eq. (3), the elements *f*_{pq} are zeroth-order contributions and are therefore combined with elements $(\gamma N)pq$ of order *n*. Conversely, the elements *W*_{pqrs} comprise the perturbation and are of first order. Hence, they are combined with elements $(\Gamma N)pqrs$ of order *n* − 1. Many properties, such as geometrical gradients, can be expressed as in Eq. (4). Since they contain the derivatives of *f*_{pq} and *W*_{pqrs}, the order in perturbation theory of the RDMs is the same as for the energy.

### D. Expectation-value expressions in MP*n* and CC*n*

The derivation of the $\u27e8\u015c2\u27e9$ expectation value in CC theory is summarized in Sec. II B. Since the derivative of the Lagrangian [see Eq. (2)] is taken with respect to *μ* after modification of the Hamiltonian according to Eq. (1), it does actually *not* lead to perturbed Fock matrices and two-electron integrals. Instead, as discussed in Ref. 4, it is found that

Note here that the two-particle RDM is contracted with combinations of overlap integrals instead of (differentiated) two-electron integrals. Since the perturbation $\u0174$ is not part of the equation, to calculate $\u27e8\u015cN2\u27e9$ within MP*n* or CC*n* theory, *both* the one- and the two-particle RDMs must be of order *n*. In the case of CC2 and MP2, they should hence both be of second order. The issue stems from the fact that the additional term $\mu \u015c2$ in Eq. (1) is neither part of the unperturbed system nor does it belong to the perturbation $\u0174$. In other words, the two-body RDM that would reproduce the MP*n* or CC*n* energy is of the wrong order in perturbation theory when calculating $\u27e8\u015cN2\u27e9$. The inconsistency that two-particle RDMs of different order are required for reproducing the energy vs for calculating the $\u27e8\u015cN2\u27e9$ expectation value is solely an issue for perturbative methods such as MP*n* and CC*n* and poses no problem for regular CC methods.

To fix the problem and in order to obtain physically meaningful results in the evaluation of Eq. (5), the two-body RDM hence needs additional contributions of order *n*. For the specific cases of MP2 and CC2, the unmodified one- and two-body RDMs for MP2 and CC2 are given in Figs. 1 and 2. Compared to CCSD, the second- or higher-order contributions $\lambda 2t2k$, with *k* ≥ 1, where *λ*_{2} and *t*_{2} are the sets of *λ* and *t* double amplitudes, respectively, are missing in CC2. To construct the modified two-particle RDM for the calculation of $\u27e8\u015c2\u27e9CC2mod$, all second-order contributions need to be added to the usual definition of the CC2 two-body RDM. These modifications are denoted as $(\Gamma N)pqrsmod$ in Fig. 3 and include all *λ*_{2}*t*_{2} terms. $(\Gamma N)pqrsmod$ for CC2 is identical to the CCSD two-particle RDM, except for the *ijab* block, which in CCSD contains additional $\lambda 2t22$ contributions. $(\Gamma N)(CC2)ijabmod$ can thus be written more compactly as

Because of the additional second-order terms, the calculation of the modified two-body RDM scales as ∼*N*^{6}. In the case of MP2, the modifications lead to the MP3 two-particle RDM.

## III. IMPLEMENTATION AND COMPUTATIONAL DETAILS

The modified two-body RDMs for the calculation of the $\u015c2$ expectation value in CC2 and MP2 have been implemented in the QCUMBRE program suite.^{20} In the case of CC2, since the scaling of the calculation cannot be lowered to ∼*N*^{5} anyways, the procedure simply uses the code already available for the two-particle RDMs in CCSD^{20–22} and removes the $\lambda 2t22$ contributions from the *ijab* block.

The necessary integrals and the UHF reference wavefunction are provided via an interface to the LONDON program package.^{23,24}

We investigate a set of molecules that consists of closed-shell cases (LiH, HF, H_{2}O, and NH_{3}) as well as open-shell cases (CH, OH, CH_{3}, CN, NO, and CH_{2}O^{+}). The calculations have been carried out using the Cartesian cc-pVTZ and cc-pVQZ basis sets.^{25} For every molecule, a corresponding geometry optimization at the same level of theory has been performed with the CFOUR program package,^{26,27} using the same basis set as for the following calculation of the $\u015c2$ expectation value.

## IV. RESULTS AND DISCUSSION

### A. Closed-shell systems

The results for the molecules LiH, HF, H_{2}O, and NH_{3} are shown in Table I. For closed-shell systems, no spin contamination should occur, neither for the reference nor at correlated levels of theory. This holds in the case of $\u27e8\u015c2\u27e9proj$. However, using the unmodified expectation-value approach $\u27e8\u015c2\u27e9$, the results are non-zero at the MP2 and CC2 levels of theory, which is unphysical.

. | LiH . | HF . | H_{2}O
. | NH_{3}
. |
---|---|---|---|---|

MP2 | ||||

E_{tot} | −8.045 149 | −100.404 641 | −76.388 141 | −56.515 938 |

$\u27e8\u015c2\u27e9unmod$ | 0.035 894 | 0.104 758 | 0.127 882 | 0.135 181 |

C_{1} | 0.035 894 | 0.104 758 | 0.127 882 | 0.135 181 |

$C2unmod$ | 0 | 0 | 0 | 0 |

$\u27e8\u015c2\u27e9mod$ | 0 | 0 | 0 | 0 |

$C2mod$ | −0.035 894 | −0.104 758 | −0.127 882 | −0.135 181 |

CC2 | ||||

E_{tot} | −8.045275 | −100.406 880 | −76.390 367 | −56.517 721 |

$\u27e8\u015c2\u27e9unmod$ | 0.035995 | 0.108 436 | 0.131 697 | 0.137 845 |

C_{1} | 0.036167 | 0.109 726 | 0.133 068 | 0.138 939 |

$C2unmod$ | −0.000172 | −0.001290 | −0.001371 | −0.001 094 |

$\u27e8\u015c2\u27e9mod$ | 0 | 0 | 0 | 0 |

$C2mod$ | −0.036 167 | −0.109 726 | −0.133 068 | −0.138 939 |

. | LiH . | HF . | H_{2}O
. | NH_{3}
. |
---|---|---|---|---|

MP2 | ||||

E_{tot} | −8.045 149 | −100.404 641 | −76.388 141 | −56.515 938 |

$\u27e8\u015c2\u27e9unmod$ | 0.035 894 | 0.104 758 | 0.127 882 | 0.135 181 |

C_{1} | 0.035 894 | 0.104 758 | 0.127 882 | 0.135 181 |

$C2unmod$ | 0 | 0 | 0 | 0 |

$\u27e8\u015c2\u27e9mod$ | 0 | 0 | 0 | 0 |

$C2mod$ | −0.035 894 | −0.104 758 | −0.127 882 | −0.135 181 |

CC2 | ||||

E_{tot} | −8.045275 | −100.406 880 | −76.390 367 | −56.517 721 |

$\u27e8\u015c2\u27e9unmod$ | 0.035995 | 0.108 436 | 0.131 697 | 0.137 845 |

C_{1} | 0.036167 | 0.109 726 | 0.133 068 | 0.138 939 |

$C2unmod$ | −0.000172 | −0.001290 | −0.001371 | −0.001 094 |

$\u27e8\u015c2\u27e9mod$ | 0 | 0 | 0 | 0 |

$C2mod$ | −0.036 167 | −0.109 726 | −0.133 068 | −0.138 939 |

To be more specific, when dividing Eq. (5) into one- and two-body contribution parts *C*_{1} and *C*_{2},

*C*_{2} should exactly cancel out *C*_{1} for the theory to be physically sound. Instead, as seen in Table I, when the unmodified densities are used to calculate these contributions at the MP2 level, *C*_{2} is zero. This can be understood because for closed-shell systems in MP2, the only non-vanishing blocks in the unmodified (Γ_{N}) are $(\Gamma N)ijab$ and $(\Gamma N)abij$. Therefore, according to Eq. (6), the overlap integrals that contribute to *C*_{2} are $S+i\u0101$, $S\u2212\u012ba$, $S+a\u012b$, and $S+\u0101i$, which for the closed-shell RHF case are all zero. The one-electron contributions *C*_{1} on the other hand are non-zero and depend strongly on the system at hand. A similar situation is observed at the CC2 level of theory. Even though *C*_{2} has non-zero contributions, it is still by two to three orders of magnitude smaller than *C*_{1}.

Since for MP2 and CC2, $(\Gamma N)pqrs$ is defined only up to the first order of the perturbation, but $(\gamma N)pq$ has contributions up to the second order, it follows that *C*_{2} cannot cancel out *C*_{1}. When this inconsistency is corrected using $(\Gamma N)pqrsmod$ defined up to the second order for MP2 and CC2, *C*_{1} and *C*_{2} correctly add up to 0.

### B. Open-shell radicals

Next, we investigate the spin contamination for the MP2 and CC2 methods in open-shell systems, specifically, in monoradical systems, for which *S*_{e.v.} = 0.75 *ℏ*^{2}. We first discuss results for CH, OH, and CH_{3}, which have a weakly spin-contaminated reference. Next, we discuss systems with a strongly spin-contaminated reference: CN, NO, and CH_{2}O^{+}. The respective results are collected in Table II.

. | CH(^{2}Π)
. | OH(^{2}Π)
. | CH_{3}($\u2009\u20092A2\u2032\u2032$)
. | CN(^{2}Σ^{+})
. | NO(^{2}Π)
. | CH_{2}O^{+}(^{2}B_{1})
. |
---|---|---|---|---|---|---|

MP2 | ||||||

E_{tot} | −38.426 400 | −75.681 812 | −39.786 070 | −92.616 357 | −129.809 600 | −114.004 055 |

$\u27e8\u015c2\u27e9ref$ | 0.759 442 | 0.756 657 | 0.761 483 | 0.985 386 | 0.774 143 | 0.783 873 |

$\u27e8\u015c2\u27e9proj$ | 0.755 549 | 0.753 838 | 0.757 175 | 0.956 069 | 0.767 243 | 0.775 764 |

$\u27e8\u015c2\u27e9unmod$ | 0.811 902 | 0.827 698 | 0.835 976 | 1.081 046 | 0.943 159 | 0.917 911 |

C_{1} | 0.060 246 | 0.076 679 | 0.083 108 | 0.154 294 | 0.182 816 | 0.150 255 |

$C2unmod$ | −0.007 786 | −0.005 639 | −0.008 615 | −0.058 634 | −0.013 800 | −0.016 217 |

$\u27e8\u015c2\u27e9mod$ | 0.753 313 | 0.752 369 | 0.754 604 | 0.925 210 | 0.762 798 | 0.770 137 |

$C2mod$ | −0.066 375 | −0.080 968 | −0.089 987 | −0.214 470 | −0.194 161 | −0.163 990 |

CC2 | ||||||

E_{tot} | −38.427 290 | −75.683 383 | −39.787 259 | −92.644 981 | −129.821 743 | −114.018 357 |

$\u27e8\u015c2\u27e9ref$ | 0.759 449 | 0.756 679 | 0.761 516 | 1.090 092 | 0.795 790 | 0.784 701 |

$\u27e8\u015c2\u27e9proj$ | 0.753 269 | 0.752 521 | 0.754 245 | 0.841 376 | 0.755 196 | 0.762 381 |

$\u27e8\u015c2\u27e9unmod$ | 0.810 700 | 0.828 373 | 0.833 967 | 0.948 336 | 0.947 794 | 0.920 700 |

C_{1} | 0.057 378 | 0.076 581 | 0.078 820 | −0.228 600 | 0.146 364 | 0.156 188 |

$C2unmod$ | −0.006 127 | −0.004 887 | −0.006 369 | 0.086844 | 0.005640 | −0.020189 |

$\u27e8\u015c2\u27e9mod$ | 0.751 332 | 0.751 190 | 0.751 708 | 0.757 992 | 0.751 318 | 0.754 574 |

$C2mod$ | −0.065 495 | −0.082 070 | −0.088 627 | −0.103 500 | −0.190 837 | −0.186 315 |

CCSD | ||||||

E_{tot} | −38.450 330 | −75.694 211 | −39.806 134 | −92.645 390 | −129.812 062 | −114.030 822 |

$\u27e8\u015c2\u27e9ref$ | 0.759 479 | 0.756 669 | 0.761 614 | 1.110 029 | 0.779 116 | 0.784 357 |

$\u27e8\u015c2\u27e9proj$ | 0.750 385 | 0.750 422 | 0.750 678 | 0.791 755 | 0.752 663 | 0.755 164 |

$\u27e8\u015c2\u27e9CC$ | 0.750 197 | 0.750 161 | 0.750 261 | 0.754 323 | 0.750 603 | 0.751 500 |

. | CH(^{2}Π)
. | OH(^{2}Π)
. | CH_{3}($\u2009\u20092A2\u2032\u2032$)
. | CN(^{2}Σ^{+})
. | NO(^{2}Π)
. | CH_{2}O^{+}(^{2}B_{1})
. |
---|---|---|---|---|---|---|

MP2 | ||||||

E_{tot} | −38.426 400 | −75.681 812 | −39.786 070 | −92.616 357 | −129.809 600 | −114.004 055 |

$\u27e8\u015c2\u27e9ref$ | 0.759 442 | 0.756 657 | 0.761 483 | 0.985 386 | 0.774 143 | 0.783 873 |

$\u27e8\u015c2\u27e9proj$ | 0.755 549 | 0.753 838 | 0.757 175 | 0.956 069 | 0.767 243 | 0.775 764 |

$\u27e8\u015c2\u27e9unmod$ | 0.811 902 | 0.827 698 | 0.835 976 | 1.081 046 | 0.943 159 | 0.917 911 |

C_{1} | 0.060 246 | 0.076 679 | 0.083 108 | 0.154 294 | 0.182 816 | 0.150 255 |

$C2unmod$ | −0.007 786 | −0.005 639 | −0.008 615 | −0.058 634 | −0.013 800 | −0.016 217 |

$\u27e8\u015c2\u27e9mod$ | 0.753 313 | 0.752 369 | 0.754 604 | 0.925 210 | 0.762 798 | 0.770 137 |

$C2mod$ | −0.066 375 | −0.080 968 | −0.089 987 | −0.214 470 | −0.194 161 | −0.163 990 |

CC2 | ||||||

E_{tot} | −38.427 290 | −75.683 383 | −39.787 259 | −92.644 981 | −129.821 743 | −114.018 357 |

$\u27e8\u015c2\u27e9ref$ | 0.759 449 | 0.756 679 | 0.761 516 | 1.090 092 | 0.795 790 | 0.784 701 |

$\u27e8\u015c2\u27e9proj$ | 0.753 269 | 0.752 521 | 0.754 245 | 0.841 376 | 0.755 196 | 0.762 381 |

$\u27e8\u015c2\u27e9unmod$ | 0.810 700 | 0.828 373 | 0.833 967 | 0.948 336 | 0.947 794 | 0.920 700 |

C_{1} | 0.057 378 | 0.076 581 | 0.078 820 | −0.228 600 | 0.146 364 | 0.156 188 |

$C2unmod$ | −0.006 127 | −0.004 887 | −0.006 369 | 0.086844 | 0.005640 | −0.020189 |

$\u27e8\u015c2\u27e9mod$ | 0.751 332 | 0.751 190 | 0.751 708 | 0.757 992 | 0.751 318 | 0.754 574 |

$C2mod$ | −0.065 495 | −0.082 070 | −0.088 627 | −0.103 500 | −0.190 837 | −0.186 315 |

CCSD | ||||||

E_{tot} | −38.450 330 | −75.694 211 | −39.806 134 | −92.645 390 | −129.812 062 | −114.030 822 |

$\u27e8\u015c2\u27e9ref$ | 0.759 479 | 0.756 669 | 0.761 614 | 1.110 029 | 0.779 116 | 0.784 357 |

$\u27e8\u015c2\u27e9proj$ | 0.750 385 | 0.750 422 | 0.750 678 | 0.791 755 | 0.752 663 | 0.755 164 |

$\u27e8\u015c2\u27e9CC$ | 0.750 197 | 0.750 161 | 0.750 261 | 0.754 323 | 0.750 603 | 0.751 500 |

For CH, OH, and CH_{3}, |Δ*S*^{ref}| is smaller than 0.02 *ℏ*^{2}. $\u27e8\u015c2\u27e9proj$ largely reduces spin contamination at both the MP2 and the CC2 levels of theory, i.e., $\u27e8\u015c2\u27e9projMP2$ is ∼0.004 *ℏ*^{2} lower than $\u27e8\u015c2\u27e9ref$, while $\u27e8\u015c2\u27e9projCC2$ is lower by ∼0.006 *ℏ*^{2}. In these cases, it is found that the stronger the contamination in the reference, the greater the improvement by the treatment of correlation.

Turning to the expectation-value approach, as seen in Table II and discussed by Chen and Schlegel,^{5} results at the CCSD level of theory show $|\Delta SprojCCSD|>|\Delta SCCSD|$. This trend, however, does not hold at the MP2 and CC2 levels of theory when using the unmodified expectation value. Specifically, the spin contamination is |Δ*S*^{unmod}| ∼ 0.1 *ℏ*^{2} > |Δ*S*_{proj}| for CH, OH, and CH_{3}. At both MP2 and CC2 levels of theory, the absolute value of the unmodified two-electron contribution $C2unmod$ is also one order of magnitude smaller than the one-electron contribution *C*_{1}. Instead, for $\u27e8\u015c2\u27e9mod$, the expected trend |Δ*S*^{mod}| < |Δ*S*_{proj}| re-emerges and the extent of spin contamination is indeed lower by ∼0.006 *ℏ*^{2} at the MP2 level and ∼0.008 *ℏ*^{2} at the CC2 level. The modified two-electron contributions $C2mod$ are now of the same order of magnitude as the one-electron contributions. For these three cases at the CC2 level, $\u27e8\u015c2\u27e9CC2mod$ ranges between 0.751 and 0.752 *ℏ*^{2}.

All three radicals CN, NO, and CH_{2}O^{+} have a strongly spin-contaminated reference. For CN, it even reaches $\u27e8\u015c2\u27e9ref=1.090\u2009\u210f2$ for a cc-pVQZ/CC2 geometry. As discussed by Krylov, spin contamination can be understood as an indicator of the failure of HF theory and suggests a multi-configurational character for the wavefunction.^{6,11} Non-dynamical electron correlation is thus expected to play an important role for these cases.

As shown by the results in Table II, MP2 does not show an adequate improvement for such strongly contaminated systems. $\u27e8\u015c2\u27e9projMP2$ calculations lower the spin contamination by about 0.01 *ℏ*^{2}. For CN, $\u27e8\u015c2\u27e9projMP2$ is 0.96 *ℏ*^{2}, while for the other two systems, it ranges between 0.77 and 0.78 *ℏ*^{2}, which is still very far from *S*_{e.v.} = 0.75 *ℏ*^{2}. The unmodified $\u27e8\u015c2\u27e9MP2unmod$ ranges around 0.9–1.1 *ℏ*^{2}. In contrast, the results for $\u27e8\u015c2\u27e9MP2mod$ are closer to *S*_{e.v.} by ∼10^{−3} *ℏ*^{2} for NO and CH_{2}O^{+} and by 0.03 *ℏ*^{2} for CN as compared to $\u27e8\u015c2\u27e9projMP2$ but are still off by >0.01 *ℏ*^{2}. These findings indicate that MP2 theory does not efficiently remove spin contamination from the wavefunction as it does not adequately account for non-dynamical correlation.

Similar to MP2, results for the unmodified expectation values $\u27e8\u015c2\u27e9CC2unmod$ show large deviations from 0.75 *ℏ*^{2}, with $|\Delta SCC2unmod|\u223c0.2\u2009\u210f2$. CN at the CC2 level of theory is the only case in which $\u27e8\u015c2\u27e9CC2unmod$ improves upon $\u27e8\u015c2\u27e9ref$. Compared to MP2, although, at the CC2 level of theory, spin contamination of the reference is largely reduced when considering $\u27e8\u015c2\u27e9projCC2$, which is 0.841 *ℏ*^{2} for CN, 0.762 *ℏ*^{2} for CH_{2}O^{+}, and 0.755 *ℏ*^{2} for NO. Again, the $\u27e8\u015c2\u27e9CC2unmod$ values seem to indicate higher spin contamination as compared to $\u27e8\u015c2\u27e9projCC2$. This behavior is corrected for $\u27e8\u015c2\u27e9CC2mod$. Comparing the two-body contribution $C2mod$ to $C2unmod$, one observes that they change sign for CN and NO. Also, $C2mod$ becomes more negative for CH_{2}O^{+}. The second-order contributions for these three cases are hence negative. The final results for $\u27e8\u015c2\u27e9CC2mod$ are promising. Treatment of electron correlation at the CC2 level manages to lower the spin-expectation value of CN from 1.050 *ℏ*^{2} for the reference to 0.758 *ℏ*^{2}. The same is true for the other two “difficult” systems, for which $\u27e8\u015c2\u27e9CC2mod$ values are closer to 0.75 *ℏ*^{2}. These observations suggest an adequate ability of CC2 to address strongly spin-contaminated systems.

Finally, we compare the MP2 and CC2 results with those at the CCSD level. As expected, the projected $\u27e8\u015c2\u27e9projCC2$ results are closer to $\u27e8\u015c2\u27e9projCCSD$ than the corresponding $\u27e8\u015c2\u27e9projMP2$ values. Specifically for CN, $\u27e8\u015c2\u27e9projCCSD$ with a value of 0.792 *ℏ*^{2} indicates a strong contamination even at the CCSD level of theory. $\u27e8\u015c2\u27e9projCC2$ is 0.841 *ℏ*^{2}, while $\u27e8\u015c2\u27e9projMP2$ is 0.956 *ℏ*^{2}, both suggesting a failure of the methods to describe the system adequately. However, focusing on the modified expectation value, one observes that the CC2 result (0.758 *ℏ*^{2}) is similar to that of CCSD (0.754 *ℏ*^{2}), while the MP2 result at 0.925 *ℏ*^{2} is still far off. The mean spin contamination of the “difficult” systems at the CCSD level $|\Delta S|\xafCCSD$ is 2.1 · 10^{−3} *ℏ*^{2}, less than half that of $|\Delta S|\xafCC2mod=4.6\u2009\u22c5\u200910\u22123\u2009\u210f2$. These results show the importance of including the *T*_{1} amplitudes at zeroth order as an approximate treatment of orbital relaxation, which are missing in MP2.

Comparing in general the different basis sets used in this study (for calculations using the cc-pVTZ basis set, see the Appendix), one notices that the differences are rather small and at the order of ∼10^{−4} *ℏ*^{2}. The spin contamination for both MP2 and CC2 is, in general, lower with the larger basis set cc-pVQZ for the systems with stronger spin contamination, but it is the choice of method for the treatment of correlation that has the deciding impact.

## V. CONCLUSIONS

In this Communication, a recipe for the calculation of the $\u015c2$ expectation value at the MP2 and CC2 levels of theory has been introduced and the basis for the generalization to the MP*n* and CC*n* series has been given. The main point here is that within perturbative schemes such as MP*n* and CC*n*, the one- and two-body RDMs used for the $\u015c2$ expectation value must be defined up to the same order of perturbation for the results to be physically meaningful. Conversely, the two-particle RDM used to compute $\u27e8\u015c2\u27e9$ is different from the one that would reproduce the energy. Usually, for methods based on the MP partitioning of the Hamiltonian such as MP*n* and CC*n*, the two-particle RDM is of order *n* − 1, while it needs to be of order *n* to compute the $\u015c2$ expectation value. There is no such inconsistency in full non-perturbative bivariational CC truncations such as CCD, CCSD, and so on. The results show that after modifying the RDMs, MP2 removes spin contamination from systems that have a weakly spin-contaminated reference but fails for strongly contaminated cases. CC2 has a superior performance over MP2 and is able to bring the spin contamination of even strongly contaminated systems down to ∼10^{−3}, which is comparable to CCSD.

## DEDICATION

S. Stopkowicz dedicates this contribution to Professor Dr. Anna Krylov. Since we met, she has been an inspiration and has become an encouraging mentor and a friend with whom to discuss both science and personal growth. I also thank Anna for actively promoting women in theoretical chemistry by starting a list of female scientists in the field and by informing about and fighting against biases.

## ACKNOWLEDGMENTS

The authors thank Professor Dr. Jürgen Gauss and Professor Dr. Anna Krylov for valuable discussions. This work was supported by the Deutsche Forschungsgmeinschaft under Grant No. STO 1239/1-1.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its Appendix.

### APPENDIX: RESULTS USING THE cc-pVTZ BASIS SET

Table III reports total energies in *E*_{h} as well as modified and unmodified $\u27e8\u015c2\u27e9$ expectation values in *ℏ*^{2} for closed-shell systems using the cc-pVTZ basis set at the MP2 and CC2 levels of theory. Table IV reports total energies in *E*_{h} as well as modified and unmodified $\u27e8\u015c2\u27e9$ expectation values in *ℏ*^{2} for open-shell systems using the cc-pVTZ basis set at the MP2, CC2, and CCSD levels of theory.

. | LiH . | HF . | H_{2}O
. | NH_{3}
. |
---|---|---|---|---|

MP2 | ||||

E_{tot} | −8.029 213 | −100.348 385 | −76.336 700 | −56.472 144 |

$\u27e8\u015c2\u27e9unmod$ | 0.032 637 | 0.098 444 | 0.121 535 | 0.129 509 |

C_{1} | 0.032 637 | 0.098 444 | 0.121 535 | 0.129 509 |

$C2unmod$ | 0 | 0 | 0 | 0 |

$\u27e8\u015c2\u27e9mod$ | 0 | 0 | 0 | 0 |

$C2mod$ | −0.032 637 | −0.098 444 | −0.121 535 | −0.129 509 |

CC2 | ||||

E_{tot} | −8.029 357 | −100.350 263 | −76.338 533 | −56.473 588 |

$\u27e8\u015c2\u27e9unmod$ | 0.032 839 | 0.101 328 | 0.124 579 | 0.131 740 |

C_{1} | 0.033 011 | 0.102 224 | 0.125 526 | 0.132 521 |

$C2unmod$ | −0.000 172 | −0.000 896 | −0.000 947 | −0.000 781 |

$\u27e8\u015c2\u27e9mod$ | 0 | 0 | 0 | 0 |

$C2mod$ | −0.033 011 | −0.102 224 | −0.125 526 | −0.132 521 |

. | LiH . | HF . | H_{2}O
. | NH_{3}
. |
---|---|---|---|---|

MP2 | ||||

E_{tot} | −8.029 213 | −100.348 385 | −76.336 700 | −56.472 144 |

$\u27e8\u015c2\u27e9unmod$ | 0.032 637 | 0.098 444 | 0.121 535 | 0.129 509 |

C_{1} | 0.032 637 | 0.098 444 | 0.121 535 | 0.129 509 |

$C2unmod$ | 0 | 0 | 0 | 0 |

$\u27e8\u015c2\u27e9mod$ | 0 | 0 | 0 | 0 |

$C2mod$ | −0.032 637 | −0.098 444 | −0.121 535 | −0.129 509 |

CC2 | ||||

E_{tot} | −8.029 357 | −100.350 263 | −76.338 533 | −56.473 588 |

$\u27e8\u015c2\u27e9unmod$ | 0.032 839 | 0.101 328 | 0.124 579 | 0.131 740 |

C_{1} | 0.033 011 | 0.102 224 | 0.125 526 | 0.132 521 |

$C2unmod$ | −0.000 172 | −0.000 896 | −0.000 947 | −0.000 781 |

$\u27e8\u015c2\u27e9mod$ | 0 | 0 | 0 | 0 |

$C2mod$ | −0.033 011 | −0.102 224 | −0.125 526 | −0.132 521 |

. | CH(^{2}Π)
. | OH(^{2}Π)
. | CH_{3}($\u2009\u20092A2\u2033$)
. | CN(^{2}Σ^{+})
. | NO(^{2}Π)
. | CH_{2}O^{+}(^{2}B_{1})
. |
---|---|---|---|---|---|---|

MP2 | ||||||

E_{tot} | −38.395 620 | −75.635 320 | −39.753 806 | −92.547 788 | −129.726 703 | −113.929 682 |

$\u27e8\u015c2\u27e9ref$ | 0.758 896 | 0.756 097 | 0.761 309 | 0.990 504 | 0.774 387 | 0.783 808 |

$\u27e8\u015c2\u27e9proj$ | 0.755 130 | 0.753 442 | 0.757 065 | 0.960 725 | 0.767 567 | 0.775 824 |

$\u27e8\u015c2\u27e9unmod$ | 0.809 024 | 0.823 065 | 0.832 958 | 1.079 052 | 0.936 750 | 0.913 148 |

C_{1} | 0.057 659 | 0.072 279 | 0.080 137 | 0.148 106 | 0.176 002 | 0.145 308 |

$C2unmod$ | −0.007 531 | −0.005 310 | −0.008 489 | −0.059 558 | −0.013 639 | −0.015 968 |

$\u27e8\u015c2\u27e9mod$ | 0.752 993 | 0.752 080 | 0.754 537 | 0.929 339 | 0.763 128 | 0.770 274 |

$C2mod$ | −0.063 562 | −0.076 296 | −0.086 910 | −0.209 271 | −0.187 260 | −0.158 842 |

CC2 | ||||||

E_{tot} | −38.396 347 | −75.636 610 | −39.754 757 | −92.575 768 | −129.738 411 | −113.943 405 |

$\u27e8\u015c2\u27e9ref$ | 0.758 902 | 0.756 116 | 0.761 337 | 1.094 079 | 0.797 966 | 0.784 529 |

$\u27e8\u015c2\u27e9proj$ | 0.753 146 | 0.752 317 | 0.754 355 | 0.845 969 | 0.755 687 | 0.762 851 |

$\u27e8\u015c2\u27e9unmod$ | 0.808 057 | 0.823 580 | 0.831 012 | 0.942 594 | 0.940 596 | 0.915 995 |

C_{1} | 0.055 191 | 0.072 027 | 0.076 021 | −0.236 940 | 0.135 475 | 0.151 898 |

$C2unmod$ | −0.006 036 | −0.004 563 | −0.006 345 | −0.098 460 | 0.007 155 | −0.020 431 |

$\u27e8\u015c2\u27e9mod$ | 0.751 283 | 0.751 079 | 0.751 811 | 0.758 678 | 0.751 246 | 0.754 858 |

$C2mod$ | −0.062 810 | −0.077 064 | −0.085 547 | −0.098 460 | −0.182 194 | −0.181 568 |

CCSD | ||||||

E_{tot} | −38.420 163 | −75.648 995 | −39.774 790 | −92.578 381 | −129.731 349 | −113.957 979 |

$\u27e8\u015c2\u27e9ref$ | 0.758 933 | 0.756 115 | 0.761 427 | 1.120 451 | 0.782 809 | 0.784 282 |

$\u27e8\u015c2\u27e9proj$ | 0.750 306 | 0 750 330 | 0.750 632 | 0.791 040 | 0.752 874 | 0.755 110 |

$\u27e8\u015c2\u27e9CC$ | 0.750 167 | 0.750 128 | 0.750 245 | 0.754 095 | 0.750 579 | 0.751 491 |

. | CH(^{2}Π)
. | OH(^{2}Π)
. | CH_{3}($\u2009\u20092A2\u2033$)
. | CN(^{2}Σ^{+})
. | NO(^{2}Π)
. | CH_{2}O^{+}(^{2}B_{1})
. |
---|---|---|---|---|---|---|

MP2 | ||||||

E_{tot} | −38.395 620 | −75.635 320 | −39.753 806 | −92.547 788 | −129.726 703 | −113.929 682 |

$\u27e8\u015c2\u27e9ref$ | 0.758 896 | 0.756 097 | 0.761 309 | 0.990 504 | 0.774 387 | 0.783 808 |

$\u27e8\u015c2\u27e9proj$ | 0.755 130 | 0.753 442 | 0.757 065 | 0.960 725 | 0.767 567 | 0.775 824 |

$\u27e8\u015c2\u27e9unmod$ | 0.809 024 | 0.823 065 | 0.832 958 | 1.079 052 | 0.936 750 | 0.913 148 |

C_{1} | 0.057 659 | 0.072 279 | 0.080 137 | 0.148 106 | 0.176 002 | 0.145 308 |

$C2unmod$ | −0.007 531 | −0.005 310 | −0.008 489 | −0.059 558 | −0.013 639 | −0.015 968 |

$\u27e8\u015c2\u27e9mod$ | 0.752 993 | 0.752 080 | 0.754 537 | 0.929 339 | 0.763 128 | 0.770 274 |

$C2mod$ | −0.063 562 | −0.076 296 | −0.086 910 | −0.209 271 | −0.187 260 | −0.158 842 |

CC2 | ||||||

E_{tot} | −38.396 347 | −75.636 610 | −39.754 757 | −92.575 768 | −129.738 411 | −113.943 405 |

$\u27e8\u015c2\u27e9ref$ | 0.758 902 | 0.756 116 | 0.761 337 | 1.094 079 | 0.797 966 | 0.784 529 |

$\u27e8\u015c2\u27e9proj$ | 0.753 146 | 0.752 317 | 0.754 355 | 0.845 969 | 0.755 687 | 0.762 851 |

$\u27e8\u015c2\u27e9unmod$ | 0.808 057 | 0.823 580 | 0.831 012 | 0.942 594 | 0.940 596 | 0.915 995 |

C_{1} | 0.055 191 | 0.072 027 | 0.076 021 | −0.236 940 | 0.135 475 | 0.151 898 |

$C2unmod$ | −0.006 036 | −0.004 563 | −0.006 345 | −0.098 460 | 0.007 155 | −0.020 431 |

$\u27e8\u015c2\u27e9mod$ | 0.751 283 | 0.751 079 | 0.751 811 | 0.758 678 | 0.751 246 | 0.754 858 |

$C2mod$ | −0.062 810 | −0.077 064 | −0.085 547 | −0.098 460 | −0.182 194 | −0.181 568 |

CCSD | ||||||

E_{tot} | −38.420 163 | −75.648 995 | −39.774 790 | −92.578 381 | −129.731 349 | −113.957 979 |

$\u27e8\u015c2\u27e9ref$ | 0.758 933 | 0.756 115 | 0.761 427 | 1.120 451 | 0.782 809 | 0.784 282 |

$\u27e8\u015c2\u27e9proj$ | 0.750 306 | 0 750 330 | 0.750 632 | 0.791 040 | 0.752 874 | 0.755 110 |

$\u27e8\u015c2\u27e9CC$ | 0.750 167 | 0.750 128 | 0.750 245 | 0.754 095 | 0.750 579 | 0.751 491 |