¹

²

³

1 Introduction↩︎

Total atomization energies (TAEs) are the thermochemical cognates of heats of formation. Computing TAEs at the sub-kJ\(\cdot\)mol\(^{-1}\)accuracy level required for such projects as Active Thermochemical Tables (ATcT)[1]–[3] requires accounting for post-CCSD(T) correlation effects. (For reaction energies that involve less drastic chemical transformations, judicious use of hierarchies of (homo/iso)desmotic reactions — see, e.g., Refs. — may obviate the need of most higher-order corrections.)

Valence post-CCSD(T) corrections are routinely evaluated in such protocols as the HEAT (high-accuracy extrapolated ab initio thermochemistry) approach of the late lamented Stanton and coworkers[4]–[7] and W4 (Weizmann-4) of Karton et al.[8], [9], as well as its explicitly correlated[10] and \(\Lambda\)-coupled cluster[11]–[14] variants.[15] Thanks to the NCC (New Coupled Cluster[16], [17]) module, built on the TBLIS tensor contraction library,[18] in the CFOUR[19] electronic structure program system, CCSDT(Q),[20] CCSDT(Q)\(_\Lambda\),[21], [22] and full CCSDTQ[23] calculations can now be carried out with unprecedented efficiency.

However, for still higher connected excitations, the only practically available implementation remains the fully general coupled cluster code in MRCC.[24]

Significant contributions of connected quintuple excitations for some molecular properties were first reported by Musiał et al.[25] (see also Refs.). Karton et al.[8] in 2006 reported that the total atomization energy (TAE, i.e., the sum of all bond energies) of ozone has a connected quintuples component of 0.3–0.5 kcal\(\cdot\)mol\(^{-1}\). Both the HEAT protocol[4]–[7] and the W4 (Weizmann-4) approach[8], [9] include connected quintuples contributions.

The extremely steep CPU time scaling \(\propto O^5V^7\) for CCSDTQ5 (where \(O\) and \(V\) represent numbers of occupied and virtual orbitals, respectively), and \(\propto O^5V^6\) for CCSDTQ(5) and CCSDTQ(5)\(_\Lambda\), makes the evaluation of \(T_5\) contributions a very costly affair. In the course of attempting to assess the importance of quintuples for a well-known multireference organic chemistry problem (vide infra) we ran into a computational scaling wall. As a field expedient, we attempted to carry out the same calculation in a reduced set of virtuals — which ultimately led us to the present study.

The idea of frozen natural orbitals (FNO) ultimately can be traced back to the iterative natural orbitals approach of Bender and Davidson.[26] Following the realization that iteratively refining natural orbitals has very poor convergence properties in practice, Meyer[27] proposed a single-shot set of NOs followed by configuration interaction (CI) in that reduced basis set; Sosa et al.[28] considered the same for many-body methods.

FNO-CC in the modern sense of the word was first proposed by Taube and Bartlett[29], [30] for ground-state coupled cluster and by Krylov[31], [32] for excited-state EOM-CC (equations-of-motion coupled cluster). The method was also implemented[33] in the MRCC[24] program system.

2 Computational Details↩︎

All electronic structure calculations in this work were carried out on the CHEMFARM HPC system at the Weizmann Institute of Science. The FNO-CCSDTQ(5)\(_\Lambda\) calculations have been carried out using the arbitrary-order coupled cluster code[21], [22], [34] in a development version of MRCC[24] kindly provided by Prof. Mihály Kállay.

The only basis sets considered in the present work are the Dunning correlation consistent[35], [36] polarized double, triple, quadruple, and quintuple-zeta basis sets (cc-pVDZ, cc-pVTZ, cc-pVQZ, and cc-pV5Z), as well as the cc-pVDZ and cc-pVTZ basis sets with the highest angular momentum removed, which we denote cc-pVDZ(p,s) and cc-pVTZ(d,p), respectively.

Reference geometries for the W4-08[37] thermochemical benchmark were taken verbatim from the ESI of the respective references. However, in order that we could reuse the CCSDTQ(5)\(_\Lambda\)/cc-pVDZ data from the W5preview paper[38], the species in Table 2 were evaluated at the CCSD(T)/cc-pwCVQZ optimized geometries from the ESI of the said paper. For subsets of W4-08, CCSDTQ(5)\(_\Lambda\)/cc-pVTZ and cc-pVTZ(d,p) calculations were carried out at the same geometries in this work.

Some CCSDTQ/cc-pVDZ and CCSDT(Q)\(_\Lambda\)/cc-pVTZ calculations on the benzyne isomers were carried out using a development version of the CFOUR[19] electronic structure program system.

In principle, one could evaluate (5)\(_\Lambda\) contributions to the total atomization energy; however, with FNO this prompts the question as to size consistency, since the sum of the numbers of FNOs for the separated atoms often does not match the number of molecular FNOs retained. However, for quintuple excitation contributions this is a nonissue for any compound in the {H,B,C,Al,Si} system as the atomic quintuples contribution will be identically zero. For the remaining elements, as seen in Table 1, quintuples contributions to the atomic energies are quite negligible, with the possible exception of chlorine atom, and even the latter would drown in remaining possible error sources (not least of which are anharmonic zero-point corrections). We have hence elected to carry out the comparison entirely in terms of the molecular correlation energy increments.

For B, C, Al, Si corrections are 0 by definition if cores are frozen.

3 Results and discussion↩︎

3.1 W4-08 thermochemical benchmark↩︎

4pt

(a) RMS value is 0.113 ; (b) values taken from ESI of Ref.; (c) Evaluated at CCSD(T)/cc-pV(Q+d)Z geometries from ESI of Ref.

.

Complete results are available in the ESI. Table 2 presents FNO-CCSDTQ(5)\(_\Lambda\) data with the cc-pVDZ basis set for a representative subsample, compared with untruncated CCSDTQ(5)\(_\Lambda\)/cc-pVDZ. The final columns also display what percentage of the full number of virtuals is included at a given FNO threshold.

A cutoff of 0.01 proved to yield too erratic results, with just a single virtual included in some cases, hence we abandoned it. FNO(0.005) includes a geometric mean of 37% of the virtual orbital space, but with a very wide spread — from just 25% for and to a whopping 65% for . Saturated hydrides aside, the second-row species clock in at markedly higher percentages than their first-row counterparts. This is probably best illustrated by Figure 1.

FNO(0.005) only recovers about half the quintuples contribution for cases like , , . FNO(0.0025) performs markedly better for such species, but seems to add little for , , . A clear performance leap occurs for FNO(0.001), with the FNO result often being very close to the full-basis number. It should be noted that, while 71% is the geometric mean percentage of the full orbital space which is included in the FNO calculations, the spread runs from 42% for methane to 96% for , with the second-row species again tending to higher percentages than their first-row counterparts. We attempted pushing further to FNO(0.0005), but clearly this only avails for (mostly first-row) cases where 0.001 retained just 60-70% of the full space. Comparison of NO histograms between cc-pVDZ and cc-pVTZ for the same species (see Figure 2 for the example of ozone) reveals that the extra basis functions in cc-pVTZ make an entry at or just below FNO(0.001), so we conclude that cutoffs below 0.001 are not meaningful unless one also wants to upgrade the underlying basis set to cc-pVTZ. (In Figure 2, a similar phenomenon is apparent for cc-pVQZ and cutoffs below 0.0001, as well as for cc-pV5Z and 0.00001, respectively.) Indeed, as seen in Table 3, already FNO(0.001) in a cc-pVTZ basis set can retain more FNOs than there would be virtuals available in a cc-pVDZ basis set.

Energy-wise, FNO(0.001) is actually a quite good approximation of the CCSDTQ(5)\(_\Lambda\)/cc-pVDZ results, with RMSD=0.013 kcal\(\cdot\)mol\(^{-1}\)out of an RMS quintuples contribution of 0.113 kcal\(\cdot\)mol\(^{-1}\).

So let us now switch gears to the cc-pVTZ basis set (Table 3). FNO(0.001) here achieves a comparable accuracy, RMSD=0.016 out of 0.099 kcal\(\cdot\)mol\(^{-1}\). Tightening to FNO(0.0005) reduces the RMSD to 0.011, and further to 0.007 for FNO(0.00025). For the nonhydrogen diatomics, we pushed one step beyond to FNO(0.0001): only for and did this change results even just at the 0.02 kcal\(\cdot\)mol\(^{-1}\)level.

The largest quintuples contributions in our sample exceed 0.4 kcal\(\cdot\)mol\(^{-1}\)for ozone and tetrasulfur, and 0.3 kcal\(\cdot\)mol\(^{-1}\)for dicarbon. Several more species have contributions of 0.2 kcal\(\cdot\)mol\(^{-1}\)  such as , , , and the like. These are also the species for which FNO(0.0025) still falls a bit short, while this appears to be adequate for the species with smaller contributions.

Typically, the quintuples contribution to the TAE rises as the FNO threshold is tightened, then plateaus close to the full-basis values, with sometimes a mild oscillation remaining. An exception is the notoriously multireference BN diatomic, where convergence is actually from above — coarse cutoffs lead to an overestimate of its importance.

A reviewer inquired about error estimates. One possibility would, of course, be to take twice the RMS error as an approximate 2-sigma error estimate. Another, more tailored to the specific system at hand, would be to consider the change between two successive FNO thresholds; this is somewhat akin philosophically to the focal point approach (FPA).[39]

Can somewhat useful results be obtained through extrapolation from fairly inexpensive FNO(0.0025)/cc-pVDZ and FNO(0.005)/cc-pVDZ values? We experimented with a Petersson-style pair natural orbital extrapolation,[40], [41] to no avail. Then we attempted an empirical two-point extrapolation of the form E[0.0025]+\(a\times\)(E[0.0025]-E[0.005]): surprisingly, simply setting \(a=1\) (in sort-of a naive Richardson extrapolation) turned out to yield the optimum RMSD=0.024 kcal\(\cdot\)mol\(^{-1}\), which may prove useful for somewhat larger systems such as those discussed anon.

A reviewer requested timings. In response, Figure 3 present a double-log plot of wall clock times on 16 Intel Ice Lake cores for three representative molecules. As can be seen there, the effective CPU time scaling is fairly close to a power law in the number of FNOs retained, with effective exponents hovering around 6.

3pt

(a) RMS value is 0.099

A reviewer wondered whether the FNO approach would also work for connected sextuples, E[\((6)_\Lambda\)]=E[CCSDTQ5\((6)_\Lambda\)]–E[CCSDTQ5]. Canonical results with the unpolarized cc-pVDZ(p,s) basis set can be found in the ESI of Ref.[15]: only a handful of systems in W4-08 have nontrivial E[\((6)_\Lambda\)] to begin with, and even those do not exceed 0.05 kcal/mol (for and ). Furthermore, since we cannot take advantage here of a split between faster and slower codes the way we could for E[CCSDTQ\((5)_\Lambda\)]–E[CCSDT\((Q)_\Lambda\)], the more meaningful quantity is really E[CCSDTQ5\((6)_\Lambda\)]–E[CCSDTQ\((5)_\Lambda\)]. Some results for , , , , and are presented in Table 4. If anything, it appears that convergence with respect to the FNO threshold is still faster than for quintuples; on the other hand, the values in question are so small, particularly for E[CCSDTQ5\((6)_\Lambda\)]–E[CCSDTQ\((5)_\Lambda\)], that their calculation is only desirable if all other contributions to the total atomization energy have been truly tied down. If however one is curious about the importance of sextuples in a species but cannot afford even a cc-pVDZ(p,s) calculation, FNO(0.001) or even FNO(0.0025) may be alternatives to consider.

3.2 Proof of principle: the benzyne isomers↩︎

\(T_3-(T)\) and (Q) extrapolated from cc-pV{D,T}Z according to Table 5 in Karton[42]

The isomers of benzyne (didehydrobenzene) are a classic multireference problem: for a detailed review of theoretical and experimental work on these species, see Vu et al.[43]

All three species have singlet ground states. However,[44] while the singlet-triplet splitting in o-benzyne and m-benzyne is quite substantial, it drops to just a few kcal\(\cdot\)mol\(^{-1}\)for p-benzyne, which in organic chemistry classes is typically treated as a biradical.

For these calculations, we used the CCSD(T)/cc-pVQZ optimized geometries from Ref. Results are presented in Table 5, while static correlation diagnostics can be found in Table 6.

The \(\mathcal{T}_1\) and %TAE[(T)] diagnostics[8], [45] for p-benzyne are actually deceptively low, just 0.019 and 4.2%, respectively. These suggest some static correlation but nothing like the alarming largest doubles amplitude max|\(T_2\)|, which rises from 0.35 at the CCSD level to a startling 0.55 at the CCSDT level (with two triples amplitudes of 0.012), rising still further to 0.60 at the CCSDTQ level (the two largest triples amplitudes now being 0.014 and 0.013). No such drastic ‘jump’ from CCSD to CCSDT is seen in o- and m-benzyne; since the newly proposed[46] \(\Delta \overline{I_\textrm{ND}}\)[(T)] and \(r_I\)[(T)] are both essentially proxies for the change in the density upon introduction of connected triples, these diagnostics do give fair warning about the p-isomer’s much more pronounced multireference character. Incidentally, so the doubles-based \(\mathcal{D}_2\) diagnostic[47], Matito’s[48] \(I_\mathrm{ND}^\mathrm{max}\) and especially \(\Delta I_\mathrm{ND}^\mathrm{max}\)[(T)] introduced in Ref., Truhlar’s M diagnostic[49] from CCSD natural orbital occupations, and the change in M upon introducing (T). (\(\mathcal{D}_2\) was evaluated using MOLPRO[50].)

Full CCSDTQ/cc-pVDZ using CFOUR on 16 Intel Ice Lake cores took about 2.5 hours per iteration for the \(D_{2h}\) p-benzyne, and about 4 hours per iteration for the lower-symmetry (\(C_{2v}\)) o-benzyne and m-benzyne isomers, compared about 10 hours per iteration for the p-benzyne triplet. (Geometries were optimized at the CCSD(T)/cc-pVQZ level using MOLPRO.)

The differences between CCSDTQ and CCSDT(Q)\(_\Lambda\) are all negative, as expected[15], becoming more so from -0.38 for ortho to -1.04 kcal\(\cdot\)mol\(^{-1}\)for para. It is known from experience[15], [51]that a degree of error compensation exists between higher-order quadruples and connected quintuples; this is borne out here too.

A full CCSDTQ(5)\(_\Lambda\)/cc-pVDZ(p,s) calculation proved elusive, but we were able to complete FNO(0.0025)/cc-VDZ for all three species. A simple linear extrapolation from FNO(0.005) and FNO(0.0025) to FNO(0) suggests quintuples contributions of {0.3,0.5,0.8} kcal\(\cdot\)mol\(^{-1}\)for {ortho,meta,para}-benzyne: these cancel the lion’s share of the \(T_4-\)(Q)\(_\Lambda\) contribution but not all of it.

FNO(0.001) proved possible only for the para isomer (thanks to its \(D_{2h}\) computational symmetry). We obtained 0.701 kcal\(\cdot\)mol\(^{-1}\), which is quite close to the value extrapolated at much smaller expense from the {0.005,0.0025} data points, 0.78 kcal\(\cdot\)mol\(^{-1}\).

It is noteworthy that across all three isomers — as indeed seen more broadly for general thermochemistry[38] and spectroscopic properties[51] — connected quintuples cancel to a large extent with \(T_4-(Q)_\Lambda\). This once more shows ‘the power of Lambda’ (John F. Stanton, 1961-2025).[51]

Overall, the sum of all the post-CCSD(T) correlation terms works out to: ortho -0.21, meta +0.01, para +1.34 kcal\(\cdot\)mol\(^{-1}\). The corresponding totals if we neglected all post-CCSDT(Q)\(_\Lambda\) terms would have been -0.11, +0.10, +1.60 kcal\(\cdot\)mol\(^{-1}\).

4 Conclusions↩︎

Connected quintuple excitations can be computationally more important than generally appreciated. For molecules with significant static correlation such as ozone, their importance in the total atomization energy can exceed 0.4 kcal\(\cdot\)mol\(^{-1}\). Unfortunately, the steep CPU time scaling of CCSDTQ5 as \(O^5V^7\), or of CCSDTQ(5)\(_\Lambda\) as \(O^5V^6\), soon puts this out of reach for all but the smallest molecules.

However, FNO-CCSDTQ(5)\(_\Lambda\) turns out to be a viable approach to surmount this daunting scaling through truncation of a frozen natural orbital space. Since the contribution of connected quintuples to atomic total energies is either zero (for B, C, Al, Si) or negligible (N–Ne, P–Ar), this truncation will not cause (meaningful) deviations from size-consistency.

FNO(0.0025)/cc-pVDZ appears to be the lowest-cost alternative that provides a meaningful estimate. FNO(0.001)/cc-pVDZ or FNO(0.001)/cc-pVTZ will yield results quite close to the full orbital space limit. With cc-pVDZ and cc-pVTZ basis sets, cutoffs below 0.001 and 0.0001, respectively, are contraindicated as they enter NO occupation regions where basis functions from the next larger basis set are truly required.

A naive extrapolation to zero cutoff from {0.005,0.0025} works surprisingly well, and (as demonstrated for the p-benzyne example) may be an attractive low-cost option for not-so-small molecules.

Aside from saturated hydrides, our results clearly indicate that FNO(Threshold) calculations on second-row species require a significantly larger percentage of the full orbital space than for the isovalent first-row molecules. This is consistent with the often-stated observation (see, e.g., Ref.) that basis set convergence is slower for second-row species than for first-row species.

We illustrate the utility of the FNO-CCSDTQ(5)\(_\Lambda\) approach through the isomers of benzyne. Connected quintuple excitations are found to reach nearly 1 kcal\(\cdot\)mol\(^{-1}\)for the para isomer, and to largely, but not perfectly, cancel with higher-order connected quadruples contributions.

CRediT authorship contribution statement↩︎

GHJ: conceptualization (equal), methodology (supporting), data interpretation (equal); writing - review and editing (equal); ADB: investigation (equal), data curation (equal), conceptualization (supporting), methodology (supporting), writing - review and editing (supporting); MS: funding acquisition (supporting), writing - review and editing (equal); JMLM: funding acquisition (lead), conceptualization (equal), supervision, methodology (lead), investigation (equal), writing - original draft, writing - review and editing (equal)

Acknowledgments↩︎

Work at Weizmann was supported by the Minerva Foundation, Munich, Germany. GHJ was funded by the National Science Foundation (grant CHE-2430408, “Advances in Coupled Cluster Theory”; deceased PI: John F. Stanton; current PI: Alberto Perez). The authors would like to thank Prof. Mihály Kállay (Technical University of Budapest) for kindly providing a prerelease version of MRCC, and Dr. James H. Thorpe (Argonne National Laboratory) and Peter R. Franke (QTP) for helpful discussions. JMLM thanks the Quantum Theory Project at the University of Florida for their hospitality. Computational work was carried out on the Faculty of Chemistry’s high-performance computing facility CHEMFARM, which is supported in part by the Ben May Center for Chemical Theory and Computation.

Conflict of interest↩︎

The authors declare no potential conflicts of interest.

Supporting information↩︎

Microsoft Excel workbook with the relevant total and interaction energies.

Additional raw data may be obtained from the corresponding author upon reasonable request.

Table of Contents Graphic (at magnification 100%)↩︎

References↩︎