Spin contamination in MP2 and CC2, a surprising issue

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.¹ 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.² 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.³ In that work, the $Ŝ2$ operator was used, with $Ŝ$ being the total spin operator $Ŝ=Ŝxi+Ŝyj+Ŝzk$ and i, j, and k being the Cartesian unit vectors. For an eigenfunction to $Ŝ2$⁠, |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 2S + 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 $〈0|Ŝ2|ΨCC〉=⟨Ŝ2⟩proj$⁠.³ 

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.⁴ In this formalism, first-order properties are calculated as derivatives of the CC Lagrangian,⁷ 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..⁴ 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⁶ 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.¹⁷ Hence, the projection approach for the calculation of the $Ŝ2$ expectation value $⟨Ŝ2⟩proj$ 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 $⟨Ŝ2⟩$ 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.

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

where $Ŝ±$ are the spin ladder operators and $Ŝz$ 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^σ$⁠, with σ = α, β, are the electron number operators for α and β electrons, respectively, and $Δn^=n^α−n^β$⁠. 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 $Ŝ2$ 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 $Ŝz$ as it preserves the number of α and β electrons,

but only the ROHF wavefunction is an eigenfunction to $Ŝ2$⁠. For this reason, for UHF based calculations, one often calculates the deviation of the expectation value $⟨Ŝ2⟩$ 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

According to the derivation of Purvis, Sekino, and Bartlett, the CC correlation contribution for the projected expression $⟨ŜN2⟩proj$ is given by³ 

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.⁴ As discussed by Stanton and Purvis, Sekino, and Bartlett,^3,4$⟨ŜN2⟩proj$ 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 $⟨ŜN2⟩CC$⁠.⁵ On practical consideration, $⟨ŜN2⟩proj$ is readily available after a CC energy calculation, while $⟨ŜN2⟩CC$ 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 $|ΨCC〉=eT^|0〉$ is not the Hermitian conjugate of the bra wavefunction $〈Ψ̃CC|=〈0|(1+Λ^)e−T^$⁠, there is no lower bound for $⟨Ŝ2⟩CC$ and spin contamination may even be negative ΔS_CC < 0.⁴ 

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)=−⟨ab|ij⟩−⟨ab|ji⟩εijab$⁠, with $εijab=εa+εb−εi−ε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 $⟨ŜN2⟩proj$ for the MBPT series.³ For the CC2 case, and generally for the CCn series, the CC amplitudes are replaced by the respective CCn solutions.

A description improving upon the projected—and hence only approximate—expressions discussed before was given by Stanton via the CC expectation value of $Ŝ2$⁠.⁴ In this formulation, the CC Lagrangian is defined as¹⁷ 

with $Λ^$ being the de-excitation operator containing the λ amplitudes $Λ^=λaii†a+14λabiji†aj†b+⋯$⁠, $T^$ being the excitation operator containing the CC amplitudes $T^=tiaa†i+14tijaba†ib†j+⋯$⁠, and $Ĥ$ being the Hamiltonian. The $Ŝ2$ expectation value is calculated as a first derivative of the Lagrangian using the modified Hamiltonian

where $Ĥ0$ 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 $Ŝ2$ operator. The derivative is then given as

which is identical to the definition of the expectation value in bivariational theory.^8,9 To obtain the correlation contribution $LCCcorr$⁠, the HF energy is subtracted from Eq. (2).

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¹⁷ 

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,¹⁹ respectively. Within a Lagrangian formulation, we may define

and take the corresponding derivative with respect to κ. Here, $H̃=e−T^ĤeT^$ and $Ĥ1=e−T^1ĤeT^1$⁠. In MP2 theory, the λ amplitudes are simply given by the complex conjugate of the first-order response amplitudes $λabij(1)=(tijab(1))*=−⟨ij|ab⟩−⟨ij|ba⟩εijab$⁠. In CC2 theory, the λ amplitudes are obtained iteratively.¹⁴ For perturbative methods such as MPn and CCn that rely on the MP partitioning of the Hamiltonian, the Fock operator $F^$ is part of the unperturbed system, while the two-electron part $Ŵ$ is the perturbation. Hence, in Eq. (3), the elements f_pq are zeroth-order contributions and are therefore combined with elements $(γN)pq$ of order n. Conversely, the elements W_pqrs comprise the perturbation and are of first order. Hence, they are combined with elements $(Γ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.

However, as seen in Sec. IV, MP2 and CC2 results based on Eq. (2) together with RDMs that fulfill Eq. (3) lead to artificial spin contamination, which is discussed in the following paragraph.

The derivation of the $⟨Ŝ2⟩$ 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 $Ŵ$ is not part of the equation, to calculate $⟨ŜN2⟩$ within MPn or CCn 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 $μŜ2$ in Eq. (1) is neither part of the unperturbed system nor does it belong to the perturbation $Ŵ$⁠. In other words, the two-body RDM that would reproduce the MPn or CCn energy is of the wrong order in perturbation theory when calculating $⟨ŜN2⟩$⁠. The inconsistency that two-particle RDMs of different order are required for reproducing the energy vs for calculating the $⟨ŜN2⟩$ expectation value is solely an issue for perturbative methods such as MPn and CCn 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 $λ2t2k$⁠, with k ≥ 1, where λ₂ and t₂ are the sets of λ and t double amplitudes, respectively, are missing in CC2. To construct the modified two-particle RDM for the calculation of $⟨Ŝ2⟩CC2mod$⁠, all second-order contributions need to be added to the usual definition of the CC2 two-body RDM. These modifications are denoted as $(ΓN)pqrsmod$ in Fig. 3 and include all λ₂t₂ terms. $(ΓN)pqrsmod$ for CC2 is identical to the CCSD two-particle RDM, except for the ijab block, which in CCSD contains additional $λ2t22$ contributions. $(Γ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⁶. In the case of MP2, the modifications lead to the MP3 two-particle RDM.

The modified two-body RDMs for the calculation of the $Ŝ2$ expectation value in CC2 and MP2 have been implemented in the QCUMBRE program suite.²⁰ In the case of CC2, since the scaling of the calculation cannot be lowered to ∼N⁵ anyways, the procedure simply uses the code already available for the two-particle RDMs in CCSD^20–22 and removes the $λ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₂O, and NH₃) as well as open-shell cases (CH, OH, CH₃, CN, NO, and CH₂O⁺). The calculations have been carried out using the Cartesian cc-pVTZ and cc-pVQZ basis sets.²⁵ 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 $Ŝ2$ expectation value.

The results for the molecules LiH, HF, H₂O, and NH₃ 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 $⟨Ŝ2⟩proj$⁠. However, using the unmodified expectation-value approach $⟨Ŝ2⟩$⁠, the results are non-zero at the MP2 and CC2 levels of theory, which is unphysical.

To be more specific, when dividing Eq. (5) into one- and two-body contribution parts C₁ and C₂,

C₂ should exactly cancel out C₁ 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₂ is zero. This can be understood because for closed-shell systems in MP2, the only non-vanishing blocks in the unmodified (Γ_N) are $(ΓN)ijab$ and $(ΓN)abij$⁠. Therefore, according to Eq. (6), the overlap integrals that contribute to C₂ are $S+iā$⁠, $S−īa$⁠, $S+aī$⁠, and $S+āi$⁠, which for the closed-shell RHF case are all zero. The one-electron contributions C₁ 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₂ has non-zero contributions, it is still by two to three orders of magnitude smaller than C₁.

Since for MP2 and CC2, $(ΓN)pqrs$ is defined only up to the first order of the perturbation, but $(γN)pq$ has contributions up to the second order, it follows that C₂ cannot cancel out C₁. When this inconsistency is corrected using $(ΓN)pqrsmod$ defined up to the second order for MP2 and CC2, C₁ and C₂ correctly add up to 0.

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 ℏ². We first discuss results for CH, OH, and CH₃, which have a weakly spin-contaminated reference. Next, we discuss systems with a strongly spin-contaminated reference: CN, NO, and CH₂O⁺. The respective results are collected in Table II.

For CH, OH, and CH₃, |ΔS^ref| is smaller than 0.02 ℏ². $⟨Ŝ2⟩proj$ largely reduces spin contamination at both the MP2 and the CC2 levels of theory, i.e., $⟨Ŝ2⟩projMP2$ is ∼0.004 ℏ² lower than $⟨Ŝ2⟩ref$⁠, while $⟨Ŝ2⟩projCC2$ is lower by ∼0.006 ℏ². 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,⁵ results at the CCSD level of theory show $|ΔSprojCCSD|>|Δ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 ℏ² > |ΔS_proj| for CH, OH, and CH₃. 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₁. Instead, for $⟨Ŝ2⟩mod$⁠, the expected trend |ΔS^mod| < |ΔS_proj| re-emerges and the extent of spin contamination is indeed lower by ∼0.006 ℏ² at the MP2 level and ∼0.008 ℏ² 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, $⟨Ŝ2⟩CC2mod$ ranges between 0.751 and 0.752 ℏ².

All three radicals CN, NO, and CH₂O⁺ have a strongly spin-contaminated reference. For CN, it even reaches $⟨Ŝ2⟩ref=1.090ℏ2$ 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. $⟨Ŝ2⟩projMP2$ calculations lower the spin contamination by about 0.01 ℏ². For CN, $⟨Ŝ2⟩projMP2$ is 0.96 ℏ², while for the other two systems, it ranges between 0.77 and 0.78 ℏ², which is still very far from S_e.v. = 0.75 ℏ². The unmodified $⟨Ŝ2⟩MP2unmod$ ranges around 0.9–1.1 ℏ². In contrast, the results for $⟨Ŝ2⟩MP2mod$ are closer to S_e.v. by ∼10⁻³ ℏ² for NO and CH₂O⁺ and by 0.03 ℏ² for CN as compared to $⟨Ŝ2⟩projMP2$ but are still off by >0.01 ℏ². 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 $⟨Ŝ2⟩CC2unmod$ show large deviations from 0.75 ℏ², with $|ΔSCC2unmod|∼0.2ℏ2$⁠. CN at the CC2 level of theory is the only case in which $⟨Ŝ2⟩CC2unmod$ improves upon $⟨Ŝ2⟩ref$⁠. Compared to MP2, although, at the CC2 level of theory, spin contamination of the reference is largely reduced when considering $⟨Ŝ2⟩projCC2$⁠, which is 0.841 ℏ² for CN, 0.762 ℏ² for CH₂O⁺, and 0.755 ℏ² for NO. Again, the $⟨Ŝ2⟩CC2unmod$ values seem to indicate higher spin contamination as compared to $⟨Ŝ2⟩projCC2$⁠. This behavior is corrected for $⟨Ŝ2⟩CC2mod$⁠. 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₂O⁺. The second-order contributions for these three cases are hence negative. The final results for $⟨Ŝ2⟩CC2mod$ are promising. Treatment of electron correlation at the CC2 level manages to lower the spin-expectation value of CN from 1.050 ℏ² for the reference to 0.758 ℏ². The same is true for the other two “difficult” systems, for which $⟨Ŝ2⟩CC2mod$ values are closer to 0.75 ℏ². 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 $⟨Ŝ2⟩projCC2$ results are closer to $⟨Ŝ2⟩projCCSD$ than the corresponding $⟨Ŝ2⟩projMP2$ values. Specifically for CN, $⟨Ŝ2⟩projCCSD$ with a value of 0.792 ℏ² indicates a strong contamination even at the CCSD level of theory. $⟨Ŝ2⟩projCC2$ is 0.841 ℏ², while $⟨Ŝ2⟩projMP2$ is 0.956 ℏ², 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 ℏ²) is similar to that of CCSD (0.754 ℏ²), while the MP2 result at 0.925 ℏ² is still far off. The mean spin contamination of the “difficult” systems at the CCSD level $|ΔS|¯CCSD$ is 2.1 · 10⁻³ ℏ², less than half that of $|ΔS|¯CC2mod=4.6⋅10−3ℏ2$⁠. These results show the importance of including the T₁ 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⁻⁴ ℏ². 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.

In this Communication, a recipe for the calculation of the $Ŝ2$ expectation value at the MP2 and CC2 levels of theory has been introduced and the basis for the generalization to the MPn and CCn series has been given. The main point here is that within perturbative schemes such as MPn and CCn, the one- and two-body RDMs used for the $Ŝ2$ 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 $⟨Ŝ2⟩$ is different from the one that would reproduce the energy. Usually, for methods based on the MP partitioning of the Hamiltonian such as MPn and CCn, the two-particle RDM is of order n − 1, while it needs to be of order n to compute the $Ŝ2$ 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⁻³, which is comparable to CCSD.

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.

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.

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

Table III reports total energies in E_h as well as modified and unmodified $⟨Ŝ2⟩$ expectation values in ℏ² 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 $⟨Ŝ2⟩$ expectation values in ℏ² for open-shell systems using the cc-pVTZ basis set at the MP2, CC2, and CCSD levels of theory.