Assessment of the accuracy of coupled cluster perturbation theory for open-shell systems. I. Triples expansions

The coupled cluster (CC) hierarchy of size-extensive correlated methods^1–3 features some of the most prominent ab initio wave function-based models in use today, namely, the CC singles and doubles (CCSD) model,⁴ the CC singles, doubles, and triples (CCSDT) model,⁵ the CC singles, doubles, triples, and quadruples (CCSDTQ) model,⁶ etc. In the limit of correlation between all N electrons of a molecular system, the full configuration interaction (FCI) energy is returned. While the iterative CC models offer excellent precision, the practical usefulness—in particular of the higher-level models for modest- and large-sized systems—is hampered by a polynomial increase of the computational cost up through the hierarchy. For this reason, a related branch of correlated methods, which incorporate non-iterative many-body perturbation theory³ (MBPT) energy corrections into the CC treatment, has been devised. However, while these perturbative corrections lower the overall cost of the computation, a proper iterative description of the effects of single and double excitations (at times even triple and quadruple excitations) within a CC scheme is called for, which is why methods have been devised in which CC and perturbation theory are coupled. The most renowned of such combined methods is the CCSD(T) model,⁷ in which the iterative CCSD energy is augmented by selected MBPT-rationalized corrections for triple excitations,⁸ and owing in part to fortuitous, yet remarkably consistent error cancellations and in part a favorable computational scaling, CCSD(T) has developed into the model of choice whenever high-accuracy total, reaction, and interaction energies, equilibrium geometries, and various other molecular properties are desired, cf., for instance, the calibration chapter of Ref. 2.

The authoritative position held today by the CCSD(T) model has been justified by the results of a great number of numerical studies and experimental comparisons, of which the clear majority has been based on closed-shell molecules in their equilibrium geometry despite the fact that UHF- and ROHF-based CCSD(T) implementations have been available almost from the beginning. However, being able to accurately account for triples effects in closed-shell systems will only provide access to a fraction of the manifold of chemical relevant problems encountered in, e.g., modern catalysis, biological, or even material sciences. Within the diversified world of reaction-based chemistry, for instance, radicals and bare atoms are omnipresent. Since such reactive species are notoriously difficult to characterize experimentally, computational chemistry, preferably of sufficiently high accuracy, hence has a substantial and valuable role to play in their description and characterization. Thus, taking into account today’s automatic predilection for the CCSD(T) model, it becomes vital to pose the following question: is it obvious that the striking cancellation of errors, which gives rise to the stability of the model, necessarily carries over from closed- to open-shell calculations? In answering this, it appears essential to rigorously assess to what accuracy CCSD(T) total energies of atoms and organic radicals can be obtained, and, if this level of accuracy is not qualitatively acceptable, establish new methods and protocols that enable a balanced treatment of both closed- and open-shell species. On this note, it has previously been shown how the precision of equilibrium geometries and harmonic vibrational frequencies of small- to modest-sized molecular radicals, as calculated by the CCSD(T) model, do not mirror that found for comparable closed-shell molecules.⁹ Furthermore, the notion of a CC/MBPT hierarchy of systematically improved models seems not to apply straightaway for open-shell molecules, as second-order Møller-Plesset (MP2)¹⁰ results are often no improvement over the underlying mean-field Hartree-Fock (HF) solution.^9,11

In the present work, we wish to suggest and advocate—as an alternative to the widely-used CCSD(T) model—the recently proposed CCSD(T–n) models,^12–14 which not only offer a convergent and well-mannered series of CC triples models for the accurate determination of total electronic energies, but do so, as we will demonstrate, regardless of the spin of the ground state. The CCSD(T–n) series is defined as an order expansion in the Møller-Plesset fluctuation potential of a bivariational CCSDT energy Lagrangian where, from a CCSD zeroth-order expansion point, perturbative solutions of both the exponential CCSDT cluster state and linear Λ-state are embedded into the energy corrections.¹⁵ In Ref. 14, the accuracy of the second- through fourth-order models of the CCSD(T–n) series (the CCSD(T–2), CCSD(T–3), and CCSD(T–4) models, respectively) was statistically assessed for total electronic energies of a test set of 17 small closed-shell molecules through a comparison with alternative CC triples expansions, including the CCSD(T) model. Furthermore, non-parallelity errors with respect to the native CCSDT solution were analyzed for the (symmetric) bond stretches in hydrogen fluoride and water. In summary, the fourth-order CCSD(T–4) model was found to offer a stable, albeit relatively costly alternative to, e.g., the CCSD(T) model, returning the smallest errors among all of the tested models. In the present study, we will extend the analysis in Ref. 14, not only by conducting a similar assessment for a test set of 18 atoms and small organic radicals, but also by reporting results for higher-order (fifth- and sixth-order) models of the CCSD(T–n) series, for total (equilibrium) electronic energies and for the symmetric bond stretches in the closed-shell water molecule and the open-shell amidogen radical, as examples of systems at the onset of bond-cleavage for which static correlation effects become increasingly important. In the second part of the present series,¹⁶ we conduct a related study for the CCSDT(Q–n) quadruples series,^12,17 which is defined—theoretically on par with the CCSD(T–n) series—as an order expansion of a bivariational CCSDTQ energy Lagrangian around a CCSDT zeroth-order expansion point.

For closed-shell systems, given that the ground state is dominated by a single determinant, the solution from a restricted HF (RHF) calculation is adequate as the reference for a subsequent treatment of correlation effects. However, for ground states with S≠0, or even for molecules which are most appropriately described by (biradical) resonance structures, the traditional approach has been to invoke different orbitals for different spins, which is what is done in a more general unrestricted HF (UHF) calculation.¹⁸ However, a UHF determinant will in general not correspond to a pure spin state, i.e., the resulting wave function will usually not be an eigenfunction of the $Sˆ2$ operator. Such contamination of the spin, arising from spurious and highly unphysical contributions from other spin states, is commonly considered a major drawback of any UHF-based method, although the inclusion of correlation effects will generally help to reduce the contamination.¹⁹ For this reason, a restricted open-shell HF (ROHF) trial function is often preferred over the corresponding UHF one whenever the UHF expectation value $〈Sˆ2〉$ deviates significantly from the correct value, i.e., whenever $〈Sˆ2〉−Sz(Sz+1)≠0$⁠.²⁰ On the other hand, the lack of spin polarization in ROHF solutions may occasionally manifest itself as a drawback and render UHF the better choice. For this reason, the choice between the two types of open-shell references is an intricate one, and the assumption that ROHF-based correlation methods should always be preferred over corresponding UHF-based methods thus cannot be expected to hold for all radicals.

For both parts of the present series, we have deliberately chosen to limit our focus to correlated treatments on top of the UHF and ROHF references, as these are the two trial functions which are most commonly and routinely used for open-shell systems. For the sake of completeness, however, we acknowledge that the use of orbitals from Brueckner²¹ or general orbital-optimized²² CC has previously been shown to reduce spin contamination in some cases^23,24 as well as to alleviate possible symmetry-breaking problems.²⁵ In addition, the use of Kohn-Sham orbitals has been advocated as an inexpensive alternative to ordinary UHF references.⁹ However, an investigation of the possible gain from using any of these non-standard determinants falls outside the scope of the present work.

All of the perturbative CC triples models tested here give corrections to the CCSD model, as formulated for either a UHF or an ROHF reference.²⁶ As in Ref. 14, the natural reference target energy, against which we wish to compare the various models, will be that of the CCSDT model.^5,27 This choice of reference is motivated not only by the fact that the CCSDT model is the theoretical yardstick for all of the approximative triples models—for restricted as well as unrestricted trial functions—but also the fact that any spin contamination of the wave function is practically removed at the CCSDT level of theory. In Table I, CC/cc-pVTZ expectation values of $〈Sˆ2〉$ are summarized for all the members of our open-shell test set (due to the frozen-core approximation, the results in Table I do not converge onto the theoretical values in the limit of full correlation among the valence electrons).²⁸ As may be recognized from the UHF- and ROHF-based results, the CCH and CN radicals (both doublets) as well as O₂ (triplet) differ from the rest of the test set, which we will further discuss in Sections III A and III B.

Besides the CCSD(T–2), CCSD(T–3), and CCSD(T–4) models of Ref. 14, we will also assess the accuracy of the fifth- and sixth-order models of the CCSD(T–n) triples series—the CCSD(T–5) and CCSD(T–6) models, respectively—for which closed- and open-shell results will be presented for the first time. Furthermore, the results for all three references (RHF, UHF, and ROHF) will be compared to corresponding results for the CCSD(T) model, the accuracy of which will be initially evaluated. In terms of computational cost, the CCSD(T) and CCSD(T–2) models both scale with system size, N, as $O(N7)$ in the non-iterative step, while all of the higher-order models of the CCSD(T–n) series (n > 2) scale as $O(N8)$ with an increasing number of cost-determining contractions upon moving up through the series, ultimately making the difference between a perturbative and a full CCSDT calculation negligible for orders n > 6.

On par with both of our previous analyses for perturbative CC triples and quadruples expansions,^14,17 which were carried out for a test set of 17 small closed-shell molecules at their equilibrium geometries,²⁹ we will here perform a statistical analysis for a test set of 18 atoms and small radicals, all optimized at the all-electron CCSD(T)/cc-pVQZ level of theory and listed in Table I. We note that this test set has previously been used in the calibration of the HEAT thermochemical model.³⁰ For measuring the numerical performance of each of the models, we report results in terms of (i) the relative recovery of the contribution to the CCSDT correlation energy from triple excitations (i.e., the CCSDT–CCSD correlation energy difference) as well as (ii) the actual deviation from this difference. All of the RHF-based calculations have been performed within the ncc module³¹ of the cfour quantum chemical program package,³² while the recently developed aquarius program³³ has been used for all of the UHF-/ROHF-based calculations (in both the ROHF-based CCSD(T–n) and CCSD(T) calculations, a semicanonical reference is used).^34–36 The correlation-consistent cc-pVTZ basis set³⁷ is used throughout for all of the reported valence-electron (frozen-core) results.

In this section, we start by investigating the performance of the CCSD(T) model against the CCSDT model for UHF and ROHF references in Section III A, while corresponding results for the CCSD(T–n) models are reported in Section III B. In both cases, we report statistical error measures generated from the individual results, cf. the supplementary material,³⁸ and these are, in turn, compared to RHF-based results for the closed-shell test set of Refs. 14 and 17. Finally, limited-range CCSD(T–n) and CCSD(T) potential energy surfaces for the symmetric vibrational modes in the water molecule and amidogen radical are reported in Section III C.

In Figure 1, we present CCSD(T) mean recoveries of/deviations from the CCSDT triples contribution to the correlation energy as well as standard deviations around the means for the RHF, UHF, and ROHF calculations, all compactly depicted in terms of normal distributions. Inspecting first the results in Figure 1(a), two things are immediately clear: first, the consistency of the CCSD(T) model in recovering CCSDT triples effects for closed-shell systems is clearly worsened in the transition from closed- to open-shell systems, and second, the statistical results for the two open-shell references appear very similar. Whereas the first observation agrees well with the aforementioned discrepancy between closed- and open-shell CCSD(T) results for equilibrium geometries and harmonic frequencies, the second is slightly surprising, taking into account that some of the UHF reference determinants may be significantly spin contaminated (cf. Table I). This becomes clearer upon inspecting the actual deviations in Figure 1(b). Here, a pronounced difference is observed between the UHF- and ROHF-based CCSD(T) results where the latter, i.e., the agreement of the CCSD(T) and CCSDT models for ROHF determinants, correspond more closely to the RHF-based results (RHF-CCSD(T) versus RHF-CCSDT for the comparable closed-shell test set). By comparing the individual entries of Tables S5 and S6 of the supplementary material,³⁸ in which results are presented for the deviations from UHF- and ROHF-CCSDT correlation energies, respectively, we find that the large difference between the UHF- and ROHF-based results in Figure 1(b) is primarily due to the inclusion of two radicals into the test set which are heavily plagued by spin contamination, namely, the challenging isoelectronic cases of the ethynyl and cyano radicals, cf. Table I. In fact, CCH and CN behave radically different from the rest of the set at the UHF-CCSD(T) level, as has already been pointed out in Ref. 30, and they are the only two radicals for which the ROHF-CCSDT triples contribution is larger than the corresponding UHF-CCSDT contribution (although less than 1.0 kJ/mol). If one were to exclude these two outliers from the statistical analysis (UHF-CCSD(T) deviations of 0.90 and 1.44 kcal/mol from the CCSDT–CCSD differences for CCH and CN, respectively), the UHF and ROHF distributions in Figure 1(b) would essentially overlap, while it would have no positive effect on the results for the recoveries in Figure 1(a).

As mentioned in Section I, the remarkable performance of the CCSD(T) model in chemical applications for closed-shell molecules (in particular for energy differences and properties) is often attributed to an overestimation of effects from triple excitations to such an extent that RHF-CCSD(T) results are in better agreement with more advanced CC models, e.g., the CCSDTQ model or even the FCI limit, than with its natural reference, the CCSDT model. By comparing individual RHF-based CCSD(T) and CCSDT correlation energies to corresponding RHF-CCSDTQ results in Tables S7 and S8 of the supplementary material,³⁸ this is indeed recognized as an occasional, although far from general feature of the model in the context of total energies (unlike for molecular properties). This is graphically illustrated in Figure 2, which gives CCSD(T) and CCSDT mean recoveries of/deviations from the combined CCSDTQ triples and quadruples contribution to the correlation energy as well as standard deviations around the means. For the closed-shell test set, the molecules that overestimate the CCSDT triples contribution the most (N₂, O₃, and CO₂) only manage to “recover” 18.0%, 12.7%, and 12.6%, respectively, of the total contribution to the CCSDTQ energy from quadruple excitations. On the other hand, we note how this potential asset of the model is in general absent for open-shell calculations, in comparing UHF and ROHF-based CCSD(T) and CCSDT results to corresponding UHF- and ROHF-CCSDTQ results (cf. Figure 2 and Tables S7 and S8 of the supplementary material³⁸). In fact, only for one open-shell species (the oxygen molecule, and only in a triple-ζ basis) does the CCSD(T) model give a marginal overestimation of the CCSDT triples contribution.

The results in Figure 1, in particular the rather poor performance of the CCSD(T) model for the two open-shell references, now raise an evident follow-up question: is it possible to improve upon the CCSD(T) model without having to resort to iterative triples models? For instance, if any perturbative model is to be an improvement over the CCSD(T) model for, say, simple atomization energies of CCSDT quality, it would have to challenge (if not improve) the level of accuracy that may be achieved for closed-shell species while improving upon it for open-shell species. In Ref. 14, we compared the CCSD(T) model to the second-, third-, and fourth-order models of the recently proposed CCSD(T–n) perturbation series (alongside a broad range of alternative triples models) and found that only the fourth-order model (CCSD(T–4)) was capable of improving the accuracy—albeit at an increased computational cost—as measured against the target CCSDT triples model. For this reason, we will compare the UHF- and ROHF-CCSD(T) models to the second- through sixth-order models of the CCSD(T–n) series for both of these references. In Figure 3, results are presented for all of these models, again with RHF-based results included for a full comparison (the CCSD(T) results for RHF, UHF, and ROHF references are the same of those in Figure 1). Inspecting first the recovery of the CCSDT triples contribution in Figures 3(a), 3(c), and 3(e), we note how the theoretically predicted convergence of the CCSD(T–n) series towards the CCSDT target energy is numerically confirmed through sixth order in the perturbation, and, most importantly, how this holds true not only for closed-shell species, but regardless of the spin of the ground state. The three different curves for RHF, UHF, and ROHF are remarkably similar, with the most notable differences being the RHF results for the CCSD(T–2) and CCSD(T–4) models, which are most likely slightly better than what should be expected, as may be seen by comparing these to the corresponding open-shell results. Furthermore, we recognize—as was hinted to in Ref. 14—how the ordinary staircase convergence of HF-based perturbation theory (MBPT) is preserved in the CCSD(T–n) series, which is based on a CCSD zeroth-order reference; while the even-ordered corrections to the CCSD energy are large, the odd-ordered corrections are generally much smaller in magnitude, although the fifth-order corrections often work to rectify the somewhat disproportionate, yet rather consistent results achieved at fourth order.

As far as the deviations of the CCSD(T–n) models from the CCSDT triples contribution in Figures 3(b), 3(d), and 3(f) are concerned, the picture is observed to be much the same. We note how the RHF-based CCSD(T–2) and CCSD(T–3) results are worse than the corresponding UHF- and ROHF-based ones, both in terms of mean and standard deviations, but how this discrepancy between the closed- and open-shell descriptions vanishes at fourth order. We find that this is due in particular to the ozone molecule, cf. Table S4 of the supplementary material,³⁸ which, at the CCSD(T–2) and CCSD(T–3) levels, has a contribution from triples excitations that deviates by more than 4.0 kcal/mol from the CCSDT reference value. At the CCSD(T–4) level, however, this deviation is reduced to 0.1 kcal/mol as a result of the relaxation and higher-order triples effects which enter the CCSD(T–n) series at fourth order.¹⁴ Finally, we note from Table S6 of the supplementary material³⁸ that O₂, when described by an ROHF reference determinant, is an odd case for the CCSD(T–n) models, giving a deviation at the CCSD(T–6)/cc-pVTZ level that is about six times larger than the second largest deviation. One possible reason for this is the relatively large deviation from the theoretical $〈Sˆ2〉$ value left at the ROHF-CCSDT level for this molecule (cf. Table I). As mentioned in Section III A, the oxygen molecule is the only member of the open-shell test set for which the CCSD(T) model overestimates the CCSDT contribution from triple excitations, and the large differences exhibited by the models of the CCSD(T–n) series, which otherwise account satisfactorily for triples effects, thus strongly indicate that the ROHF-CCSD(T) result for this molecule might be rather artificial.

In the course of a chemical reaction, a molecule is often distorted away from its equilibrium geometry along a given internal coordinate. In such situations, the stretching of bond(s) will introduce multiconfigurational character, that is, static correlation, into the wave function, which makes a mean-field HF description of the reference state increasingly worse. An immediate consequence of this is that the overall description of the system becomes exceedingly more difficult. Furthermore, the weights of higher-level excited determinants in the wave function will grow, making it inherently more demanding to apply perturbative models for the recovery of these pronounced effects. Thus, in order to further test the models of the CCSD(T–n) series—despite the fact that neither of these have been designed specifically to recover static correlation effects—we will here evaluate the models for the symmetric bond stretches in the closed-shell water (H₂O) molecule and the open-shell amidogen (NH₂) radical and compare them in terms of their non-parallelity errors with respect to the target CCSDT solution (the validity of which we will initially assess).

For both water (RHF reference) and amidogen (UHF reference) in a cc-pVTZ basis (C_2v point group), the error of the CCSDT solution with respect to that of the more advanced CCSDTQ model is positive and increases up to a certain bond length, at which the CCSDT solution collapses. Thus, it is therefore only of interest to test the performance of the CCSD(T–n) models within this region; for H₂O in Figure 4(a), the CCSDT solution starts to collapse at an elongation of 80%, while for NH₂ in Figure 5(a), this critical point is met at an elongation of 60%.

In Figures 4(b) and 5(b), the total deviation from the correlation energy at the CCSDT/cc-pVTZ level of theory is reported for both of these stretches, again presented alongside corresponding results for the CCSD(T) model, which—like the models of the CCSD(T–n) series—is designed with the purpose of recovering dynamical correlation only. If we define a conservative confidence interval to be spanned by a maximum error of ±0.2 kcal/mol, we note how for the water stretch in Figure 4(b), the results for the CCSD(T–2) model all fall outside the confidence interval, the results for the CCSD(T–3) model nearly so, while the CCSD(T) results stay within up until a stretch of about 60%. Only the higher-order CCSD(T–n) models are capable of staying within the interval for all tested bond lengths with, e.g., the performance of the CCSD(T–4) model only showing signs of deteriorating near the very end of the tested interval of bond lengths. For the NH₂ stretch in Figure 5(b), both the native CCSD and CCSDT models as well as the CCSD(T–n) models are observed to behave similarly to the case of the H₂O stretch in Figure 4. However, while the CCSD and CCSDT models deviate less from the CCSDTQ solution than for H₂O, the CCSD(T–n) models (as measured against the CCSDT model) deviate slightly more. For the CCSD(T) model, on the other hand, the performance for the two test systems is markedly different; whereas the CCSD(T) results for NH₂ are similar to the corresponding CCSD(T–n) results, that is, the model consistently (and increasingly) underestimates the CCSDT triples effects, this is not the case for H₂O, for which the error is positive only up to a certain point (a stretch of about 40%), at which the model suddenly starts to overestimate these effects.

Thus, whenever high-accuracy potential energy surfaces in the vicinity of the equilibrium geometry are required, as, e.g., when these form the basis for vibrational configuration interaction or coupled cluster response theory^39,40 or in the determination of activation energies, the higher-order models of the CCSD(T–n) series might offer a stable non-iterative alternative to the CCSDT model. In particular, the fact that triples relaxation effects are accounted for in the fourth- and higher-orders models of the series will make these significantly more stable, albeit more expensive, than various other non-iterative models, as seen, e.g., from the comparison with the CCSD(T) model in Figures 4 and 5 and various other triples models in Ref. 14 (in which a similar analysis of the bond stretch in hydrogen fluoride is also presented).

By evaluating frozen-core/cc-pVTZ total energies against CCSDT reference results, the performance of the second- through sixth-order models of the recently proposed CCSD(T–n) series of triples expansions for UHF and ROHF reference determinants has been statistically assessed for a test set of 18 atoms and small radicals. Furthermore, these results have been compared to (i) the performance of the RHF-based variants of the models for a comparable test set of 17 closed-shell molecules as well as (ii) corresponding CCSD(T) results for all three references. In summary, no differences are found between the treatments of closed- and open-shell species by the higher-order CCSD(T–n) models, neither for the size-intensive recovery of the CCSDT triples contribution nor the actual deviation from this. Moreover, the results presented here strongly indicate, as is the case for regular iterative CC models, that the higher-order CCSD(T–n) models are capable of significantly reducing the amount of spin contamination which may be present in underlying UHF trial functions, since no behavioural differences are observed between the UHF- and ROHF-based CCSD(T–n) results reported herein. These models, in particular the fourth-order CCSD(T–4) model, thus offer proper non-iterative, albeit more computationally expensive alternatives to the current de facto standard, the CCSD(T) model, for which the appraised numerical consistency in recovering CCSDT triples effects is worsened in the transition from closed- to open-shell systems, irrespective of the choice of reference determinant used for the latter type (UHF/ROHF).

In addition, the performance assessment of the CCSD(T–n) models has been extended to two prototypical bond stretches (the symmetric vibrational modes in the closed-shell H₂O molecule and open-shell NH₂ radical), as examples of systems for which static correlation effects become increasingly dominant. As for the calculations at equilibrium geometries, the higher-order models of the CCSD(T–n) triples series are found to offer results that are in significantly better agreement with the native CCSDT solution than the traditional CCSD(T) model for molecules at distorted geometries, as, e.g., encountered along reaction paths or at transition states. For instance, the performance of the CCSD(T–4) model is again found to be essentially independent of the reference used, and it will thus be capable of offering a balanced non-iterative treatment of closed- and open-shell species as required in, e.g., the determination of activation energies, even at intermediate bond lengths where alternative triples expansions may start to collapse.

J.J.E. and P.J. acknowledge support from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. 291371 and The Danish Council for Independent Research – Natural Sciences. J.G. acknowledges financial support from the Deutsche Forschungsgemeinschaft (Grant No. DFG GA 370/5-1). D.A.M. acknowledges support from the US National Science Foundation (NSF) under Grant No. ACI-1148125/1340293.