We compare the numerical performance of various non-iterative coupled cluster (CC) quadruples models. The results collectively show how approaches that attempt to correct the CC singles and doubles energy for the combined effect of triple and quadruple excitations all fail at recovering the correlation energy of the full CC singles, doubles, triples, and quadruples (CCSDTQ) model to within sufficient accuracy. Such a level of accuracy is only achieved by models that make corrections to the full CC singles, doubles, and triples (CCSDT) energy for the isolated effect of quadruple excitations of which the CCSDT(Q–3) and CCSDT(Q–4) models of the Lagrangian-based CCSDT(Q–*n*) perturbation series are found to outperform alternative models that add either of the established [Q] and (Q) corrections to the CCSDT energy.

In computational quantum chemistry, the coupled cluster (CC) family of methods^{1} and many-body perturbation theory (MBPT)^{2} both offer a systematic approach toward the full configuration interaction (FCI)^{3} limit—the exact solution of the time-independent, non-relativistic electronic Schrödinger equation—within a given one-electron basis set. The non-iterative second-order Møller-Plesset (MP2) model^{4} gives the lowest-order perturbative correction to the Hartree-Fock (HF) energy for the isolated effect of connected double excitations. As an improvement, the iterative CC singles and doubles (CCSD) model^{5} accounts not only for connected doubles but also does so to infinite order in the space of all single and double excitations out of the reference, which further introduces disconnected double, triple, etc., excitations in the wave function as a direct consequence of the exponential CC ansatz. In most cases, however, the accuracy of the CCSD model is not adequate for the study of chemical phenomena, and higher-level connected excitations (triples, quadruples, etc.) also have to be accounted for. In the CC singles, doubles, and triples (CCSDT) model,^{6} triple excitations enter the CC cluster operator. However, while the CCSDT model offers a clear improvement over the MP2 and CCSD models in terms of accuracy, the difference in computational cost between the CCSD model ($O ( N 6 ) $) and the CCSDT model ($O ( N 8 ) $) is extremely steep (where *N* is a measure of the total system size).

For this reason, a plethora of approximate triples models have been devised over the years, scaling either iteratively or non-iteratively as $O ( N 7 ) $ or non-iteratively as $O ( N 8 ) $. For a recent review and a numerical comparison of various non-iterative CC triples models, cf. Refs. 7 and 8, respectively. Importantly, most of these approximate methods combine coupled cluster and perturbation theory, either by forming higher-order corrections to the CCSD energy from non-iteratively constructed coupled cluster amplitudes, or by modifying the CCSDT equations on the basis of a perturbation analysis, and approaches based on either the HF^{9–11} or coupled cluster^{7,12–14} reference states have been proposed, many of which are capable of approximating the CCSDT energy to within a reasonable degree of accuracy.

If one, however, is aiming for an increasingly more fine-grained level of detail, which is a crucial requirement in, for instance, the field of computational thermochemistry^{15–19} and for the precise determination of molecular equilibrium geometries,^{20–23} even CC triples models (including the full CCSDT model) will often fail to provide the level of accuracy required, which is often in the sub-kJ/mol range. In these cases, including the effect of connected quadruple excitations proves necessary. Alas, the full iterative CC singles, doubles, triples, and quadruples (CCSDTQ) model^{24} exhibits a staggering asymptotic $O ( N 10 ) $ scaling, which is bound to prevent it from routine applications in anything but the niche of atoms and modest-sized molecules. For this reason, CC quadruples models have been devised that scale non-iteratively as anything from $O ( N 7 ) $ to $O ( N 10 ) $. In the present communication, we wish to investigate to what extent non-iterative CC quadruples models rationalized from either HF or CC state-based perturbation theory manage to recover the energy of the iterative CCSDTQ model. In the course of doing so, we will also present the first results obtained with the lowest-order models of the CCSD(TQ–*n*) and CCSDT(Q–*n*) series,^{7} which are derived from Lagrangian-based CC perturbation theory for a CCSD or CCSDT zeroth-order state, respectively.

In order to conduct a proper investigation on the performance of various perturbative quadruples models, numerical comparisons will be made between a broad spectrum of established triples and quadruples models from the literature as well as the recently proposed CCSD(TQ–*n*) and CCSDT(Q–*n*) models. The various models that make corrections to the converged CCSDT energy for the isolated effect of connected quadruple excitations will be compared to one another (the CCSDT(Q–*n*), CCSDT[Q],^{10} CCSDT(Q),^{13} ΛCCSDT[Q], and ΛCCSDT(Q)^{25} models), while the models that correct the converged CCSD energy for the combined effect of triple and quadruples excitations (the CCSD(TQ–*n*), CCSD(TQ_{f}),^{10} and CCSD+TQ^{*}^{11} models) will not only be compared against one another but also against two triples-only models, namely, the CCSD(T)^{26} and CCSDT^{6} models.^{27}

We begin by giving a brief introduction to each of the tested models. (i) By truncating the CC energy Lagrangian^{3} at a given excitation level in the cluster operator, we have recently shown how corrections to any CC state of interest can be determined that formally converge toward a chosen target CC limit in terms of term-wise size-extensive energy contributions.^{7,14} For example, the energy corrections of the models that form the CCSD(T–*n*) perturbation series^{8} are derived by parametrizing the CCSDT energy Lagrangian around the CCSD energy as the reference point. Likewise, the models of the CCSD(TQ–*n*) series are derived by parametrizing the CCSDTQ energy Lagrangian around the CCSD energy, while the CCSDT(Q–*n*) models are derived by parametrizing the CCSDTQ energy Lagrangian around the CCSDT energy. (ii) The CCSDT[Q]^{10} and CCSDT(Q)^{13} models are straightforward quadruples extensions of the traditional CCSD[T]^{9} and CCSD(T)^{12,26} triples models, which are rationalized from MBPT and therefore based on the underlying HF state.^{28} This is true also for the ΛCCSDT[Q] and ΛCCSDT(Q)^{25} models, which differ from the CCSDT[Q] and CCSDT(Q) models by expressing the quadruples energy corrections in terms of both converged CCSDT amplitudes and multipliers (or Λ state parameters), which correctly account for relaxation of the CC amplitudes. These models are thus distinguished from the CCSDT[Q] and CCSDT(Q) models in the same manner as the ΛCCSD[T]^{8} and ΛCCSD(T)^{29,30} models are distinguished from the CCSD[T] and CCSD(T) models. (iii) The CCSD+TQ^{*} model^{11} extends the CCSD[T] correction to quadruple excitations by including all terms which enter the CCSD energy expression at fifth order in MBPT. This makes the quadruples contribution factorizable into terms of at most $O ( N 7 ) $ cost, as in the fourth- and fifth-order Møller-Plesset schemes, while the fifth-order triples corrections scale as $O ( N 8 ) $. Similarly, the CCSD(TQ_{f}) model^{10} constructs a CCSDT[Q]-like correction on top of the second-order triples amplitudes as defined in the CCSD[T] and CCSD(T) models. However, this approach differs from the CCSD+TQ^{*} approach by not including any fifth-order triples corrections beyond the ones already present in the CCSD(T) model, and an enforced factorization of the quadruples terms is thus required (although in practice, this is often a very good approximation).

In numerically comparing the various quadruples models, we perform calculations on the same test set of 17 small closed-shell molecules at their equilibrium geometries—all optimized at the all-electron CCSD(T)/cc-pCVQZ level of theory—as was previously used in Refs. 3 and 8 for benchmarking the performance of the CCSD(T) and CCSD(T–*n*) models, respectively. For measuring the recovery of the contribution to the CCSDTQ correlation energy from quadruple excitations or the combined effect of triple and quadruple excitations, we use two size-intensive performance measures, $ \Delta i Q $ and $ \Delta i TQ $, for the recovery of the CCSDTQ/CCSDT and CCSDTQ/CCSD differences, respectively. Similarly, we present results in terms of actual differences from the CCSDTQ reference (*δ _{i}*) where this is either the total CCSDTQ correlation energy or atomization energy.

In all of the valence-electron (frozen-core) calculations of the present communication, we use the Dunning correlation-consistent cc-pVDZ and cc-pVTZ basis sets.^{31} Statistical results in terms of mean recoveries ($ \Delta \xaf $) and mean differences ($ \delta \xaf $) as well as standard deviations around the means (Δ_{std} and *δ*_{std}, respectively) will be compactly depicted in terms of normal distributions as in Ref. 8. For individual results, cf. the supplementary material.^{32} All of the closed-shell calculations have been performed within the NCC module^{33} of the cfour quantum chemical program package,^{34} while all of the atomic open-shell calculations have been performed within the recently developed Aquarius program.^{35}

In Figure 1, we present normal distributions of the recovery of CCSDTQ/CCSDT correlation energy differences for the perturbation models that make quadruples corrections to the CCSDT energy. As is clear from comparing the performance of the models in Figure 1(a) with those in Figure 1(b), none of the (Λ)CCSDT(Q)/[Q] models are capable of recovering the CCSDTQ correlation energy with the same consistent accuracy as the CCSDT(Q–3) and CCSDT(Q–4) models. The poor performance of the CCSDT[Q] model previously reported by Kállay and Gauss^{25} is confirmed herein, and as may be recognized by comparing individual cc-pVDZ and cc-pVTZ results (cf. Tables S1 and S3 of the supplementary material^{32}), the severe basis set dependence of the [Q] correction^{13} is replicated by the more advanced ΛCCSDT[Q] model. This behavior is less pronounced for the CCSDT(Q–2) model, even less so for the CCSDT(Q) and ΛCCSDT(Q) models, and completely absent for the CCSDT(Q–3) and CCSDT(Q–4) models. Finally, we note how the (Q) correction is always too negative (except for the difficult case of methylene, cf. the discussion in Ref. 8), while the CCSDT(Q–2) correction is consistently not negative enough, the CCSDT(Q–3) correction lies centred around the CCSDTQ reference, and the CCSDT(Q–4) correction is the best with only a minor positive error.

For the models in Figure 2, which are all designed to recover the combined effect of triple and quadruple excitations on top of the CCSD energy, we note how the most advanced of these, the CCSD(TQ–4) model, is the only one that shows any improvement over the CCSDT model in Figure 2(a). However, since the energy gap between CCSDTQ and CCSD is always larger than the corresponding difference in energy between CCSDTQ and CCSDT, the results in Figures 1 and 2 cannot be compared directly to one another. For this reason, we have also compared the various models based on absolute rather than relative energy differences. These results are given in Figure 3 for the best among the CCSDT and CCSD state-based models and again reported in full in the supplementary material (Tables S3 and S4^{32}). As is clear from Figure 3, the CCSDT(Q–3) model gives results in better agreement with the CCSDTQ reference than the CCSDT(Q–4) model, despite being less consistent with respect to relative errors (cf. Figure 1(b)). Furthermore, we note that even the CCSD(TQ–4) model, being the best of the models in Figure 2, cannot be used for targeting the CCSDTQ model, as the absolute errors are too large and erratic, and we thus conclude that triples and quadruples effects cannot generally be treated at the same level, i.e., a full iterative account of triples effects is integral before quadruples effects can be addressed. This is also exemplified by comparing the performance of any of the quadruples models in Figure 2(a) to that of the CCSD(T) and CCSDT models, which only make a correction for triples effects. Finally, by comparing the CCSD(T) and CCSDT normal distributions in Figure 2(a), we conclude that the beneficial cancellation of errors, which at times brings the CCSD(T) model in closer agreement with the CCSDTQ than the CCSDT model, is far from general.

Although the relative and absolute recoveries of the quadruples contribution to the total CCSDTQ correlation energy in Figures 1–3 are fair error measures on which we can compare the various models, neither of them will necessarily reflect the actual consistency of any of the models in potential chemical applications, as this most often relies on the ability of a given model to describe energy differences rather than absolute energies accurately, and neither will the measures say anything about how the individual models might advantage from possible general cancellations of errors for such applications. For this reason, we have compared those models that exhibit the smallest differences from the total CCSDTQ correlation energies on the basis of a third measure, which is the capability of the models to recover atomization energies calculated at the frozen-core CCSDTQ/cc-pVTZ level of theory. While atomization energies calculated using the CCSDTQ model in themselves will not reflect reality, as various other effects such as core correlation, zero-point vibrational motion, and spin-orbit corrections^{16–19} are not accounted for, the remaining contributions from higher-level excitations (pentuples, hextuples, etc.) may be assumed to be minor in comparison to the contributions from excitations up to the level of quadruples.

By comparing individual CCSDT and CCSDTQ atomization energies (Table S5 of the supplementay material^{32}), it is obvious that the case of ozone stands out with a total contribution from quadruple excitations of 14 kJ/mol, which is more than 10 kJ/mol in excess of the results for all of the other tested molecules. For this reason, we have chosen to exclude ozone from the statistical test set, which has a pronounced positive effect on the results for the CCSDT(Q) and ΛCCSDT(Q) models, while less so for the models of the CCSDT(Q–*n*) series. Plotting the final results for the atomization energies in Figure 4, however, reveals how the most consistent and precise results are still those obtained with the higher-order CCSDT(Q–*n*) models; while the CCSDT(Q–3) model as in Figure 3 shows the lowest standard error both with and without ozone included in the test set, only the CCSDT(Q–4) model exhibits a consistent sign of the error with respect to the CCSDTQ results across the test set.

Finally, the question remains of how much one *gains* by using any of the non-iterative quadruples models in Figure 4 over the full CCSDTQ model? The answer to this question depends on both the level of accuracy obtained as well as the computational cost of the method. For instance, while thermochemical quantities such as atomization energies most often have to be determined to within an error margin of 1 kJ/mol, larger errors might be tolerated for other application purposes, which could possibly allow for a less-accurate, but more cost-efficient model to be used. While all of the results and discussions up to this point have revolved around the potential of each of the tested quadruples models in terms of numerical accuracy, we will now comment on the computational cost associated with using any of the most accurate of these for three representative cases of varying size and difficulty; H_{2}O, C_{2}H_{4}, and O_{3}.

In Table I, timings for the individual models are presented relative to the time required for the iterative CCSDTQ solution. Generally, while the absolute cost grows with the level of complexity of the quadruples models, we note how the relative cost drops notably whenever the cost of the CCSDTQ solution increases (from water over ethene to ozone as well as upon increasing the basis set). In both basis sets, the cost of computing the (Q) correction on top of the CCSDT energy, which scales non-iteratively as $O ( N 9 ) $, is marginal. For example, the increase in compute time over the CCSDT solution for ozone in the cc-pVDZ basis is no more than a factor of 1.6 for the CCSDT(Q) model and 3.8 for the more advanced ΛCCSDT(Q) model, in which also the CCSDT Λ equations are converged, while these ratios trivially grow in the larger cc-pVTZ basis. Turning to the third- and fourth-order CCSDT(Q–*n*) models, the additional accuracy of these is found to come at a price, as the cost of computing any of the corrections is significantly larger than that of the (Q) correction; while the CCSDT(Q–2) correction has a formal cost that is slightly less than a single CCSDTQ iteration (but still scales non-iteratively as $O ( N 10 ) $), the CCSDT(Q–3) and CCSDT(Q–4) corrections add approximately the cost of one and two full CCSDTQ iterations, respectively. Comparing their cost to that of the ΛCCSDT(Q) model, however, reveals how the CCSDT(Q–3) and CCSDT(Q–4) corrections only increase the total time-to-solution by a factor of 5 and 10, respectively, in the cc-pVTZ basis, and less so in the smaller cc-pVDZ basis.

Basis set . | CCSDT . | CCSDT(Q) . | ΛCCSDT(Q) . | CCSDT(Q–2) . | CCSDT(Q–3) . | CCSDT(Q–4) . |
---|---|---|---|---|---|---|

cc-pVDZ | 10.3/1.7/0.9 | 10.6/2.4/1.4 | 19.8/5.5/3.3 | 22.8/7.9/5.2 | 28.2/11.1/7.8 | 42.1/21.6/15.3 |

cc-pVTZ | 0.9/0.3/0.2 | 1.9/0.9/0.7 | 3.8/1.6/1.2 | 5.9/3.0/2.6 | 11.6/7.8/5.8 | 23.2/15.4/11.0 |

Basis set . | CCSDT . | CCSDT(Q) . | ΛCCSDT(Q) . | CCSDT(Q–2) . | CCSDT(Q–3) . | CCSDT(Q–4) . |
---|---|---|---|---|---|---|

cc-pVDZ | 10.3/1.7/0.9 | 10.6/2.4/1.4 | 19.8/5.5/3.3 | 22.8/7.9/5.2 | 28.2/11.1/7.8 | 42.1/21.6/15.3 |

cc-pVTZ | 0.9/0.3/0.2 | 1.9/0.9/0.7 | 3.8/1.6/1.2 | 5.9/3.0/2.6 | 11.6/7.8/5.8 | 23.2/15.4/11.0 |

^{a}

Number of CCSDTQ iterations required (cc-pVDZ/cc-pVTZ): H_{2}O (11/11); C_{2}H_{4} (14/15); O_{3} (21/20).

In summary, we conclude that only those non-iterative perturbation models that make an exclusive quadruples correction to the converged CCSDT energy work in practice for any chemical purpose that necessitates high accuracy. Thus, despite being the most accurate among the models that make a correction to the CCSD energy for the combined effect of triple and quadruple excitations, the higher-order models of the CCSD(TQ–*n*) series still perform as poorly as those models that add the [Q] correction to the CCSDT energy. For high-accuracy applications, we find the following hierarchy of models:

In particular, the CCSDT(Q–3) model is found to perform well. While on average the model recovers total CCSDTQ correlation energies less systematically than the more advanced (and expensive) CCSDT(Q–4) model, it benefits from a consistent cancellation of errors which is observed across our test set. The maximum error with respect to atomization energies calculated at the frozen-core CCSDTQ/cc-pVTZ level of theory is thus below 0.3 kJ/mol, and for the most part, results obtained with the CCSDT(Q–3) model are improved when moving from a double-*ζ* to a triple-*ζ* basis set. Although the CCSDT(Q) model offers a reasonable compromise between cost and accuracy, its general applicability is found to be impaired by a few of the molecules of the present test set (e.g., O_{3} and multiply bonded systems such as CO_{2} and N_{2}) — all molecules for which results of CCSDTQ quality are only obtained with the CCSDT(Q–3) and CCSDT(Q–4) models. We also note that its applicability to molecules with heavy atoms has previously been questioned,^{36} but to what extent similar caveats exist for the models of the CCSDT(Q–*n*) series is outside the scope of the present communication. In conclusion, we advocate the use of the CCSDT(Q–3) and CCSDT(Q–4) models as alternatives to the often used CCSDT(Q) model for various branches of thermochemistry, while the performance of the models in determining bond lengths and vibrational frequencies needs further testing.

J.J.E. and P.J. acknowledge support from the European Research Council under the European Union’s Seventh Framework Programme (No. 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 (No. DFG GA 370/5-1). D.A.M. would like to acknowledge support from the US National Science Foundation (NSF) under Grant No. ACI-1148125/1340293.

## REFERENCES

We will not present any explicit numerical comparison against the alternative CCSD(TQ) and CCSD+TQ models, as tests performed by us on the systems of the present test set have shown the differences between the CCSD(TQ_{f}) and CCSD(TQ) models and the CCSD+TQ^{*} and CCSD+TQ models to be negligible.