Elaborated thermochemical treatment of HF, CO, N₂, and H₂O: Insight into HEAT and its extensions

In the last few decades, a number of quantum chemical protocols have been developed to compute molecular bond energies and enthalpies of formation to within 1 kJ/mol. This level of accuracy has been called “sub-chemical”¹ and is roughly a fourfold improvement over the 1 kcal/mol associated with the term “chemical accuracy.”² While these protocols are typically used for straightforward thermochemical applications, they are also used to compute quantities such as ionization potentials,^3,4 electron affinities,^5–7 barrier heights for chemical reactions,^8–11 and molecular structures.^11–13 Among electronic structure specialists, it is appreciated that these methods necessarily benefit from some degree of error cancellation; direct (i.e., non-extrapolated and non-composite) determination of the electron correlation energy contribution to a bond energy or enthalpy of formation to within 1 kJ mol⁻¹ (about 80 cm⁻¹) remains beyond the reach of computation for most molecules due to practical limitations. One purpose of this paper is to explore some aspects of the systematics of this error cancellation for the title molecules. This work targets such issues within the context of the HEAT (High-Accuracy Extrapolated Ab Initio Thermochemistry) family of methods; for additional studies similar in spirit to the present study, the reader is referred to the work of Martin and Karton on the quite similar Wn family of methods^14–16 and the analyses of Feller et al.^17–19 

In addition to HEAT, which was developed in collaboration involving groups in the US, Germany, and Hungary,^20–23 the Wn family of methods by Martin’s group at the Weizmann Institute,^6,24–27 the focal point analysis (FPA) by Allen and co-workers,^1,28–30 FPD by Feller, Peterson, and Dixon,^31–33 and ANL-n by Klippenstein, Harding, and Ruscic³⁴ can achieve sub-chemical accuracy on a reasonably routine basis. All of these approaches attempt to approximate the full-CI energy at the complete basis-set limit. While HEAT and Wn are most closely related, all of these “model chemistries” employ a broadly similar strategy: the total energy is partitioned into a number of contributions that are assumed to be additive, and each is computed to some necessarily approximate level of theory. Estimates and extrapolations are then carried out to refine the values of each contribution. The FPD and FPA methods differ from HEAT, Wn, and ANL-n in that the former foregoes a fixed-recipe prescription in favor of flexible, molecule-specific energy components. All such approaches, fixed-recipe or not, are sometimes termed “composite methods.” The electronic energy [consisting of Hartree–Fock Self-Consistent-Field (HF-SCF) and electron correlation contributions] is the most critical and generally uncertain aspect of this procedure. Beyond the vibrational zero-point energy (ZPE), which is significant in magnitude and important at both the “sub-chemical” and “chemical” accuracy levels, there are smaller contributions from the effects of special relativity, spin–orbit coupling (for degenerate states of radicals, the difference between the lowest spin–orbit level and the averaged state typically computed in quantum chemical calculations), and the so-called diagonal Born–Oppenheimer correction.³⁵ It is tacitly assumed, as will be done here, that additional and smaller effects (such as quantum electrodynamics and non-adiabatic corrections) can be neglected.

Among the protocols discussed above, the HEAT family of methods are the most computationally demanding, largely because they do not separate core and valence correlation contributions at the coupled-cluster singles, doubles, and perturbative triples, or CCSD(T), level. HEAT generally provides sub-chemical accuracy for “well-behaved” species (those that exhibit only mild multireference character and are well represented by the Born–Oppenheimer approximation) that comprise first- and second-row atoms (Z = 1 through 10). Some recent and, as yet, unpublished work in our group on applying HEAT to molecules containing sulfur and chlorine, however, has provided convincing evidence that a logical extension of the protocol for third-row atoms (Z = 11 through 18) is flawed, serving as a reminder that the protocol as first designed is “lucky” in the sense of benefiting from fortuitous error cancellation.

This work, although originally motivated by these findings, remains focused on molecules comprised of first- and second-row atoms and looks at these cancellation effects in some detail. Chosen for the study is a set of small molecules: HF, CO, N₂, and H₂O. While these are all closed-shell species containing light atoms, the set spans a range of chemical properties: fluorine-containing compounds are notoriously difficult problems for quantum chemistry, N₂ and CO contain triple bonds, and H₂O contains two highly polar bonds. The virtue of the species chosen here is that they are small enough to be studied at extremely high levels of theory with extensive basis sets. In addition, the level of precision with which the exact thermochemical quantities are known, which became possible only with the advent of Active Thermochemical Tables (ATcT),³⁶ is better than 0.05 kJ mol⁻¹, so it is possible to know the “right” answer at an accuracy that makes the present study meaningful.

The procedure is as follows: By appealing either to experiment or to improved computational methods for comparison, we study how close the standard HEAT treatments of the zero-point energies, diagonal Born–Oppenheimer, and relativistic corrections are to convergence. While this information is interesting by itself, precise knowledge of these parameters together with the overall ATcT thermochemistry permits one to determine the contribution to the bond energies/enthalpies of formation that strictly result from the Born–Oppenheimer electronic energy. These “exact” electronic total atomization energy (eTAE) contributions are then used as benchmarks to guide our investigation of electron correlation and basis-set effects (particularly, the coupling between these), beyond those included in the HEAT protocol, and to explore the systematics that underscore the standard HEAT protocol for the title molecules.

All quantum chemical calculations reported here were performed with the cfour^37,38 and mrcc^39,40 packages. The former was used for all calculations apart from contributions beyond CCSDTQ for closed-shell molecules and post-CCSDT calculations for the (open-shell) atoms, both of which used the general coupled-cluster implementation of mrcc. The largest all-electron CCSD(T) and CCSDT(Q)_Λ calculations reported here, with up to 1138 and 362 basis functions, respectively, were made possible by the new NCC module of cfour.^41–43 All atomic open-shell energies used an unrestricted Hartree–Fock (UHF) reference function. Symmetry equivalencing of the atomic orbitals was not performed.

Basis sets used in the evaluation of the correlation energy comprised of the standard aug-cc-pCVXZ (through X = 6),^44–46 as well as extensions for X = 7 and 8 that were constructed for the purpose of this study. Details of the septuple- and octuple-zeta basis sets for hydrogen, carbon, nitrogen, oxygen, and fluorine are given in the supplementary material. The use of such large basis sets in polyatomic calculations carries with it an inherent tendency to suffer from linear dependency issues, in addition to the rapidly prohibitive computational cost.

While zero-point energies (ZPEs) for the diatomic species (HF, CO, and N₂) can be determined via coefficients of the Dunham expansion,⁴⁷ a highly accurate ZPE for water must be calculated. To this end, the “reference” ZPE for water was determined from a well converged, variational, discrete variable representation vibrational calculation using the nitrogen program^48,49 and spectroscopic quality potential energy surfaces available in the literature.^50,51

Reference scalar-relativistic contributions to the atomic and molecular energies were taken as the energy difference between a spin-free Dirac–Coulomb (SFDC)^52,53 and equivalent non-relativistic calculation for a given species. Uncontracted aug-cc-pCVTZ and aug-cc-pCVQZ basis sets were used to determine basis-set dependence, and calculations were performed at the SCF, coupled cluster singles and doubles (CCSD), and CCSD(T) levels of theory to gauge the convergence of this contribution with respect to the treatment of electron correlation.

The full details of the HEAT family of methods are discussed at length elsewhere,^20–23 but are summarized here. The HEAT treatment is based on the following partitioning of the energy, each term of which is determined at the equilibrium geometry obtained for the molecule in question at the all-electron CCSD(T)/cc-pVQZ level of theory:

For the present discussion, these contributions are conveniently grouped as below.

In HEAT, the SCF energy is extrapolated using the three-point exponential formula previously advocated by Feller,⁵⁴ 

using three successive members of the aug-cc-pCVXZ sequence of basis sets (X = 3 for aug-cc-pCVTZ, X = 4 for aug-cc-pCVQZ, etc.). In HEAT, two choices are used; the extrapolation in HEAT-345 uses the basis sets with X = T, Q, 5, while HEAT-456 uses X = Q, 5, 6.

The SCF energy is then augmented by an all-electron CCSD(T) correlation correction, extrapolated by the two-point scheme of Helgaker et al.,⁵⁵ viz.,

where X has the same meaning as above and $ΔEcorr.∞$ is the estimated basis-set limit of the all-electron CCSD(T) correlation energy. In HEAT-345, this involves the aug-cc-pCVXZ, X = Q, 5 basis sets, while HEAT-456 uses X = 5, 6.

The remaining contributions to the correlation energy, those beyond CCSD(T), are resolved into two parts. The first is the difference between the extrapolated, frozen-core CCSDT and CCSD(T) energies, using the cc-pVTZ and cc-pVQZ basis sets, ΔE_CCSDT. Finally, higher-level corrections (HLC) are approximated by unextrapolated, frozen-core energy differences of higher-order cluster excitations such as CCSDT(Q)–CCSDT [HEAT-ABC(Q), where ABC = 345 or 456], CCSDTQ–CCSDT (HEAT-ABCQ), or CCSDTQP–CCSDT (HEAT-ABDQP), which aims at accounting for differences between CCSDT and the full-CI limit. These latter calculations are performed using the smallest correlation consistent basis set, cc-pVDZ. This choice of basis was essential at the time of the original HEAT protocol given the considerable cost of higher-order correlation calculations. Both of these post-CCSD(T) corrections are done in the frozen-core approximation, in contrast to the CCSD(T) correlation energy contribution that correlates all electrons.

The nuclear motion terms in HEAT come from an all-electron, CCSD(T)/cc-pVQZ level VPT2 calculation of the zero-point energy and an all-electron, CCSD/aug-cc-pCVQZ (or frozen-core, aug-cc-pVTZ) diagonal Born–Oppenheimer correction (the so-called “DBOC”).^56–58 

The spin–orbit correction in HEAT, which represents the difference between the lowest spin–orbit level of an atom or molecule and the averaged spin–orbit state that is computed by conventional quantum chemical methods, is obtained from experiment, when available. The scalar-relativistic correction is calculated with one- and two-electron mass–velocity and Darwin corrections, using first-order perturbation theory. These calculations are carried out at the CCSD(T) level with the aug-cc-pCVTZ basis.

In the original HEAT formulation,²⁰ and for historical reasons based on program limitations, the (perhaps sub-optimal) combination of all-electron CCSD(T) and the cc-pVQZ basis set were used to determine the molecular geometries at which thermochemical properties were calculated. At the time, this struck a working balance between accuracy and computational expense and did not appear to introduce significant errors; the molecular enthalpies of formation matched those of ATcT to within 1 kJ mol⁻¹ for the most part. Here, however, where accuracy on the order of 10–20 cm⁻¹ (0.12–0.24 kJ mol⁻¹) is desired, the choice of molecular geometry becomes important. To this end, the calculations presented in the rest of this paper were performed at the established experimental or highly converged computational equilibrium geometries.^59,60

The goal of the present study is to investigate systematic deficiencies with various elements of the HEAT energy partitioning [Eq. (1)], particularly those associated with the treatment of electron correlation. Accordingly, we have adopted the following strategy. For the small systems studied here, it is possible to obtain essentially exact (within a few cm⁻¹) results for the zero-point energy, relativistic correction, and DBOC contributions to the total atomization energies (TAEs). The zero-point energies can be taken from experimental studies of diatomics or well converged variational calculations on spectroscopically accurate potential energy surfaces of water. The accuracy of the HEAT-based relativistic corrections is assessed in Sec. III A of this report, as is that of the relatively insignificant DBOC correction.

With the essentially exact values for the contributions mentioned in the preceding paragraph, the very precise TAEs for these species that are available from ATcT (utilizing ATcT Thermochemical Network version 1.124⁶¹) allow one to accurately determine the contribution of the non-relativistic electronic energy to the TAE (eTAE for short). We then present a series of calculations that reproduce these eTAE values to within 12 cm⁻¹. This set of data is then used to inspect the performance of the standard HEAT-345(Q) procedure and, by extension, other members of the HEAT family. Finally, some aspects of the electronic energy calculations are discussed, and suggestions for future, highly accurate model chemistries are proposed.

The relativistic contributions to the total atomization energies of HF, CO, N₂, and H₂O are given in Table I. While the first-order atomic spin–orbit corrections are obtained from experiment, the scalar-relativistic corrections for both the atoms and molecules must be calculated. Reference SCF, CCSD, and CCSD(T) SFDC relativistic corrections were obtained as described in Sec. II. This treatment of the scalar-relativistic contribution to the TAE appears to be converged to, roughly, <3 cm⁻¹ (<0.04 kJ mol⁻¹), given the slight basis-set dependence of the CCSD correlation contributions (<1 cm⁻¹ for all species) and the small contribution of perturbative triples to the CCSD correlation correction (<4 cm⁻¹ for all species). Previous studies⁶² have shown that basis-set convergence of relativistic contributions to molecular energies can be relatively slow. However, the scalar-relativistic contributions to the barrier height of water (see above reference) and to the total atomization energies studied here appear to converge far faster than the corrections to the individual atoms/molecules, especially with uncontracted aug-cc-pCVXZ basis sets. It should also be noted here that second-order spin–orbit coupling has been neglected. However, it is anticipated that this contribution is negligible for these species.

It appears that the relativistic corrections in standard HEAT approaches closely approximate these reference values: the first-order spin–orbit correction is “exact,” and the MVD2 (mass–velocity as well as one- and two-electron Darwin corrections), all-electron CCSD(T)/aug-pCVTZ treatment of the scalar-relativistic term only differs from the converged spin-free Dirac–Coulomb approach by at most 3 cm⁻¹ (for N₂). These molecules consist of relatively “light” atoms; however, beyond the second row, it is likely that the HEAT treatment of relativity will be degraded and ultimately prove inadequate.

The nuclear motion contributions to the total atomization energies are given in Table II. The standard HEAT treatment of the DBOC appears to be converged to within a few cm⁻¹, with a larger and systematic error for the zero-point energy correction. The DBOC changes by at most −10 cm⁻¹ (for water) between SCF and CCSD wavefunctions with aug-cc-pCVTZ and aug-cc-pCVQZ basis sets; correlation effects beyond CCSD also appear to be very small (less than 2 cm⁻¹ for all species), and we estimate that the converged results collected in Table II are accurate to 1 cm⁻¹.

As mentioned in previous publications, the zero-point energy is a significant source of systematic error in HEAT. The all-electron, CCSD(T)/cc-pVQZ ZPE is consistently larger than the “exact” values contained in Table II, leading to an underbinding of the molecule (lowered TAE). The worst offender here is water, with an error of over 28 cm⁻¹, or 0.34 kJ mol⁻¹ compared to the reference value. Using frozen-core CCSD(T)/ANO2, a combination of theory and basis that has been shown to accurately reproduce fundamental frequencies,⁶³ improves these errors (the largest is then 8 cm⁻¹ for H₂O) and distributes them around zero (ZPEs for HF and H₂O are too large, while those for CO and N₂ are too small). It is likely that this error will grow with system size; the numbers given here clearly scale roughly linearly with the magnitude of the ZPE.

Variational calculations for these species were also performed, where the Watson Hamiltonian was expanded in a direct product basis of 11 harmonic oscillator functions for each normal mode, and the necessary potential energy integrals were evaluated by Gauss–Hermite quadrature at the frozen-core CCSD(T)/ANO2 level of theory. The Coriolis and pseudopotential integrals were calculated using the VPT4 expressions for the Hamiltonian.⁶⁴ For the title molecules, the zero-point energy provided by VPT2 agrees with this variational calculation to within 0.1 cm⁻¹, implying that the error in the ZPE for these molecules arises from the quality of the potential energy surface or underlying rovibrational Hamiltonian rather than deficiencies of second-order vibrational perturbation theory. That said, these are relatively well-behaved species from the point of view of VPT2; the anharmonicity is dominated by stretching anharmonicity (as VPT2 is exact for a Morse potential, it tends to treat stretching anharmonicity with great fidelity), they do not exhibit large amplitude motion that can plague the theory,⁶⁵ and there are no couplings between modes beyond three-body terms. Vibrational Diffusion Monte Carlo (DMC) calculations by Harding et al. indicate that such good agreement between variational calculations and VPT2 should not be expected in general.⁶⁶ 

Given these accurate zero-point energies, experimental spin–orbit values, and the essentially converged DBOC and scalar-relativistic corrections, we estimate that the errors in the nuclear motion and relativistic contributions to the total atomization energies of the species studied here are on the order of 5 cm⁻¹ (0.06 kJ mol⁻¹). This is a fairly conservative estimate, which qualitatively can be regarded as being consistent with the 95% confidence interval (two standard deviations) of the standard ATcT uncertainties. Thus, reckoning the empirical TAEs from the ATcT TAEs by subtracting off the nuclear and relativistic contributions and using standard error propagation with the 0.06 kJ mol⁻¹ estimate above and the 95% confidence interval from ATcT should give overall eTAE uncertainties roughly consistent with the total TAE uncertainties from ATcT. The largest error estimate for the electronic TAE is 0.08 kJ mol⁻¹ or 6 cm⁻¹, for CO and N₂, as seen in Table III. We now turn to the task of developing a computational strategy capable of recreating these non-relativistic, electronic TAE values.

In keeping with the HEAT procedure, electronic energy contributions are partitioned as follows: the SCF energy, a CCSD(T) correlation correction, the T–(T) correction, and a two-part higher-level correlation correction. The guiding principle is the same as that of HEAT: as far as possible, all correlation corrections will be obtained without the frozen-core approximation and will be based on extrapolations using large aug-cc-pCVXZ basis sets. The resulting values are presented in Table IV.

In Table IV, the SCF contribution to the eTAE of our test species is expanded in basis sets up to aug-cc-pCV8Z. This contribution seems to be converged to within 2 cm⁻¹ by aug-cc-pCV5Z and to below 1 cm⁻¹ by aug-cc-pCV6Z. As has been commented on before,²² the three-point exponential SCF extrapolation using aug-cc-pCV{T,Q,5}Z serves to hinder, rather than help, the SCF TAE contribution. The larger extrapolations make less than 1 cm⁻¹ difference to this TAE contribution. We have elected to take the unextrapolated aug-cc-pCV8Z value as the best $EHF∞$ contribution to the TAE. We estimate that these values are converged to within 1 cm⁻¹.

To determine the best CCSD(T) correlation contribution to the eTAE, this term was evaluated up to the very large aug-cc-pCV8Z basis set. Following the analysis of Feller et al.,^18,19 we have elected to employ the two-point, $(X+12)−4$ extrapolation formula of Martin,⁶⁷ rather than the two-point, X⁻³ formula of Helgaker⁵⁵ used in a typical HEAT calculation. This basis-set extrapolation is critical even when using basis sets up to aug-cc-pCV8Z: the extrapolation length (distance from the largest raw calculation and extrapolated limit) for all four molecules is above 25 cm⁻¹ and is nearly 100 cm⁻¹ (or 1.18 kJ mol⁻¹) in the case of N₂. We have selected the aug-cc-pCV{7,8}Z extrapolated CCSD(T) correlation energy as the best $ΔECCSD(T)∞$ contribution to the TAE and estimated it to be converged to at most roughly 12 cm⁻¹, based on the differences between the aug-cc-pCV{6,7}Z and the aug-cc-pCV{7,8}Z extrapolated values.

The so-called “full-T” TAE contribution used here differs significantly from that of HEAT. For the desired high-accuracy TAEs, the all-electron (rather than valence-only) CCSDT and CCSD(T) energies are calculated with aug-cc-pCVXZ basis sets up through aug-cc-pCV5Z and extrapolated with the two-point formula of Martin discussed above. Perhaps unsurprisingly, the aug-cc-pCVDZ basis produces wildly inaccurate T–(T) contributions to the TAE and serves to contaminate the aug-cc-pCV{D,T}Z extrapolated values. However, the basis sets beyond double-ζ produce smoothly convergent T–(T) corrections, and the aug-cc-pCV{Q,5}Z extrapolation length for this quantity is less than 6 cm⁻¹ for all four species. We have selected the aug-cc-pCV{Q,5}Z extrapolated T–(T) as the best choice for this contribution to the TAE and estimate it to be converged to within 6 cm⁻¹.

The adopted higher-level correlation contributions to the TAE are obtained by a two-step process. First, an all-electron CCSDT(Q)_Λ–CCSDT correction is applied to account for $T̂4$ contributions of the coupled-cluster expansion. Then, following the recommendation of Karton et al.,¹⁴ effects beyond $T̂4$ are estimated via a valence-only CCSDTQ(P)_Λ correction. It should be pointed out that the choice of CCSDT(Q) vs CCSDTQ vs CCSDT(Q)_Λ for the $T̂4$ correction does not actually matter, provided that the subsequent $T̂5$ correction sufficiently accounts for the shortcomings of the selected $T̂4$ approach. We elected to use CCSDT(Q)_Λ because testing proved that it resulted in very small $T̂5$–$T̂4$ TAE contributions, a desirable property when basis-set convergence is too expensive to obtain.

As with the T–(T) corrections, the (Q)_Λ–T TAE contributions are obtained in all-electron calculations using basis sets up to aug-cc-pCV5Z, which is probably the most computationally demanding step in the procedure that has been adopted here. While not to the same extent as the T–(T) correction, the aug-cc-pCVDZ basis set still performs poorly, overestimates the (Q)_Λ–T TAE contributions by up to 26 cm⁻¹ in the case of N₂, and once again contaminates the aug-cc-pCV{D,T}Z extrapolated values. However, all values obtained with basis sets beyond aug-cc-pCVDZ agree to within 10 cm⁻¹, and the aug-cc-pCV{Q,5}Z extrapolation length is at most 3 cm⁻¹ for CO. Thus, we have selected the extrapolated aug-cc-pCV{Q,5}Z value for the all-electron (Q)_Λ–T portion of the ΔE_HLC contribution and estimate it to be converged within 2 cm⁻¹.

Even with modern hardware, the valence-only (P)_Λ–(Q)_Λ TAE contributions proved exceedingly expensive to evaluate beyond cc-pVTZ (for CO). It is not clear if the cc-pVDZ basis set can be trusted for extrapolating this contribution, but it is comforting that the cc-pVTZ and cc-pVDZ (P)_Λ–(Q)_Λ values agree to 2 cm⁻¹ for these four species. We have selected the cc-pVTZ value for the valence-only (P)_Λ–(Q)_Λ portion of ΔE_HLC contribution and estimated it to be converged within 2 cm⁻¹.

The results of our best current ab initio, HEAT, and empirical non-relativistic, electronic total atomization energies are presented in Table V. The HEAT electronic TAE differs from the nearly exact empirical eTAE by at least a half kJ mol⁻¹ for three of the four species. Inclusion of the nuclear and relativistic contributions results in a cancellation of errors that brings the HEAT values more in line with the ATcT TAE: the largest HEAT error is 0.57 kJ mol⁻¹ for N₂ (0.32, 0.09, and 0.31 kJ mol⁻¹ for HF, CO, and H₂O, respectively). While the HEAT-345(Q) TAEs shown in Table V are still quite accurate, the set of calculations chosen for the purposes of this study far exceeds HEAT-345(Q) in performance (and expense). Our current best electronic TAE estimates match the empirical values (and, therefore, the full TAE and enthalpies of formation) to a remarkable degree: 0.14 kJ mol⁻¹ (12 cm⁻¹) or better for all four species (12, −11, −2, and −1 cm⁻¹ for HF, CO, N₂, and H₂O, respectively).

Propagating the errors of the individual components of the best ab initio results (see Table IV) yields estimated 95% confidence intervals of ±12, 8, 13, and 13 cm⁻¹ for HF, CO, N₂, and H₂O, respectively. These uncertainties are roughly consistent with the agreement of said values and those of ATcT. A larger test suite of molecules would be required to draw any conclusions about the general performance of the method presented above.

Using the results obtained with the approach used here to calculate very accurate TAE contributions, we now take a look at the systematics of what makes HEAT “work” in the first place.

For the moment, we will focus solely on HEAT-345(Q). As demonstrated in Table VI, HEAT-345(Q) and other HEAT variants owe some of their typically quite excellent performance to a fortuitous cancellation of errors between various components. We will now go through these one by one.

As mentioned above, the ZPE is one source of systematic error in the HEAT procedure. For all four molecules here, the HEAT ZPE results in TAEs that are too negative compared to experimental (or theoretically nearly exact) values. The order of magnitude varies, but this effect is strongest for HF and H₂O (−0.17 and −0.34 kJ mol⁻¹, respectively, or −14 and −28 cm⁻¹). This error is due to the sub-optimal combination of the valence cc-pVQZ basis set used in combination with an all-electron CCSD(T) wavefunction, which results in geometries that are too “tight” and consequently overestimated ZPE values, a problem that likely scales with the magnitude of the ZPE.

The next source of error in HEAT-345(Q) is the SCF energy extrapolation to the CBS (complete basis set) limit. As commented upon in a previous publication,²² and demonstrated in Table IV, the exponential extrapolation formula employed in HEAT actually hinders the convergence of the SCF energy. The resulting TAEs are too negative by more than 0.13 kJ mol⁻¹ for all the molecules here, and a surprising 0.27 kJ mol⁻¹ for CO. The SCF TAE contribution has very nearly converged by the aug-cc-pCV5Z basis set, and inclusion of the less accurate aug-cc-pCVTZ basis set in the 345 extrapolation only serves to corrupt this value. This error can be entirely avoided by simply dropping the extrapolation and going with the energy obtained in the largest basis set used in the calculation.

The CCSD(T) CBS correlation energy, usually the first or second most significant contribution to the TAE, is also the largest source of error in HEAT-345(Q). This contribution differs from our best values by more than 0.70 kJ mol⁻¹ for three of the four species studied here, and there is an alarming 1.54 kJ mol⁻¹ discrepancy for N₂. This arises from a combination of two factors: the inherently slow convergence of the CCSD component of the CCSD(T) correlation energy when expanded in a single-particle basis, and a tendency for the Helgaker, X⁻³ extrapolation employed in HEAT to overshoot the actual basis-set limit. The former can be ameliorated by employing a sufficiently large basis: up to 7-, 8-ζ in this case. Other authors have commented that the $(X+12)−4$ extrapolation proposed by Martin (and adopted here) has generally good performance.^18,19,67 In brief, the X⁻³ extrapolation systematically “over-extrapolates” the correlation energy of the molecules relative to the atoms, causing the TAE to be too positive. The $(X+12)−4$ extrapolation is apparently better-balanced and leads to TAEs that are too small, but by a lesser magnitude than the X⁻³ values are too large. It is probably best to exploit explicitly correlated methods, which are ideally suited for two-particle correlation methods such as MP2 and CCSD, to obtain the latter, and then deal with the triples separately, an approach that has been advocated by Martin, but such an investigation in the present context lies outside the scope of this work. Interested readers are directed to a wealth of prior literature^{18,19,27,67–72} for discussion of strategies for the CCSD and CCSD(T) energy contribution in composite methods.

One of the more expensive components of the HEAT procedure is the valence, cc-pV{T,Q}Z extrapolated T–(T) correction, which attempts to account for potential deficiencies of the treatment of triple excitations in CCSD(T). Implicit here are two assumptions: that the T–(T) correction is sufficiently converged when extrapolating from the cc-pVTZ and cc-pVQZ basis sets and that core-correlation beyond CCSD(T) can be ignored. Both assumptions are investigated in Table VI. It seems that while the HEAT treatment of the T–(T) TAE contribution does capture the CBS extrapolation of the frozen-core energies for these molecules, the neglected core-correlation for both CO and N₂ contributes −0.14 and −0.24 kJ mol⁻¹ to the electronic TAE, respectively. These findings are consistent with those found in studies based on W4, where core-correlation beyond CCSD(T) is accounted for explicitly.²⁵ Here, where the goal is accuracy on the order of 10–20 cm⁻¹, the difference between CBS all-electron and frozen-core T–(T) TAE contributions is significant. It should be noted that the X⁻³ and the $(X+12)−4$ extrapolations produce only minimal differences for correlation contributions beyond CCSD(T) when using sufficiently large basis sets, as is consistent with the recent findings of Karton.¹⁶ 

The last component of the HEAT treatment of the electronic TAE is the so-called “higher-level correlation” correction, which attempts to account for the remaining valence correlation beyond CCSDT via an energy difference between CCSDT(Q), CCSDTQ, or CCSDTQP (depending on the variant of HEAT employed) and CCSDT using the cc-pVDZ basis set. In Sec. III B, this contribution was calculated using a two-step process: a converged CCSDT(Q)_Λ–CCSDT component and a subsequent, valence-only, CCSDTQ(P)_Λ–CCSDT(Q)_Λ correction. In order to analyze the performance of HEAT, we have also performed HLC calculations using a (Q)–T and (P)_Λ–(Q) sequence with the same basis sets, which is more directly comparable.

It turns out that the use of the practical, but quite small, cc-pVDZ basis set in the HLC correction is responsible for the second largest error in the HEAT procedure for these molecules. As shown in Table VI, the (Q)–T TAE correction computed using cc-pVDZ is far from the frozen-core, CBS limit calculated using aug-cc-pCVXZ, X = {Q,5} basis sets. The sign of this error is not consistent: the cc-pVDZ (Q)–T TAE correction for both HF and H₂O is too positive by 0.25 and 0.14 kJ mol⁻¹, respectively, whereas for CO and N₂ (comprising triple bonds and, therefore, more demanding of electron correlation), it is 0.38 and 0.69 kJ mol⁻¹ smaller than the frozen-core CBS results. Additionally, there remains a small degree of core-correlation in the (Q)–T correction. This is most important for CO where it accounts for −0.08 of the −0.45 kJ mol⁻¹ total error for this term.

The final, remaining difference between our best values and HEAT-345(Q) is the “missing” (P)_Λ–(Q) correction, which accounts for a non-negligible 0.22 and 0.17 kJ mol⁻¹ of the CO and N₂ TAE, respectively. However, this “missing” correction serves to offset the −0.45 and −0.74 kJ mol⁻¹ errors accrued in the (Q)–T correction of CO and N₂. Fortunately, this can be lessened significantly by employing CCSDT(Q)_Λ rather than CCSDT(Q) (or CCSDTQ) in the first part of the HLC corrections for these molecules. It would be desirable to also test the use of CCSDTQP in the HLC corrections, but getting anywhere close to basis-set convergence for this calculation remains beyond the reach of computational resources, even for these small molecules.

The calculations reported here shed some light on aspects of the systematic error cancellation from which the HEAT-345(Q) method benefits. Specifically, for the molecules studied here, the analysis of Sec. III D reveals the following features of HEAT calculations of total atomization energies:

• ZPE contributions are systematically too small (up to −0.34 kJ mol⁻¹ for H₂O).

• SCF contributions are systematically underestimated (up to −0.27 kJ mol⁻¹ for CO).

• CCSD(T) contributions systematically overshoot (up to 1.54 kJ mol⁻¹ for N₂).

• T–(T) contributions miss core-correlation of up to −0.24 kJ mol⁻¹ for N₂.

• (Q)–T contributions suffer from basis-set incompleteness of up to −0.69 kJ mol⁻¹ for N₂.

• (P)_Λ–(Q) corrections are unaccounted for up to 0.22 kJ mol⁻¹ for CO.

The most conspicuous aspect of this behavior is that the most significant error among the contributions to the composite energy, that of the CCSD(T) correlation energy contribution, has one sign, and the bulk of the remaining deficiencies in the methodology have errors in the opposite direction. As a result, the overall accuracy of HEAT-345(Q) atomization energies is significantly better than what might have been expected were this cancellation not operative. Indeed, long ago when the corresponding author of this manuscript first compared HEAT-based enthalpies of formation to the (then new) ATcT values, he was stunned by a level of accuracy that seemed unwarranted by the approximations made. Looking at the numbers in the manner outlined in this work is certainly illuminating; the method as originally devised, indeed, had shortcomings in its contributing terms, but an underlying tendency toward cancellation that was certainly not obvious. Therefore, while one could definitely improve certain contributions–for example, the overestimation of the vibrational zero-point energy due to reliance on a molecular geometry that tends to be “too tight” would clearly be improved by using a method that does a better job of getting vibrational frequencies correct (such as CCSD(T)/ANO1⁶³)–it is not at all clear that this would improve the overall performance of the method. For molecules comprised of first- and second-row atoms, HEAT-345(Q) works reasonably well, and it should probably be left alone. To obtain a systematic improvement over this method requires that all contributions be improved simultaneously, in a manner accomplished in this work. The ZPE needs to be improved, the intrinsic shortcomings of extrapolation need to be muted through use of very large basis sets, one must go well beyond a double-zeta frozen core calculation for high-level correlation effects, core-correlation beyond CCSD(T) is required, one should consider dealing with CCSD and (T) separately (the former through explicitly correlated calculations to minimize the extrapolation error), and so on.

In the context of the above, it is now possible to perform significantly more quantitative calculations of thermochemical parameters than was possible a decade ago. The revolutionary mrcc program of Kállay brought general high-level coupled-cluster calculations to practitioners of quantum chemistry, and continued improvements in software and hardware now enable such calculations to be done with very large basis sets. The treatment of post-CCSD(T) correlation energy (in terms of basis-set completeness) done in this work is well beyond that of the standard model chemistries, and it will not be long before such calculations can be done routinely on somewhat larger molecules. At present, we are in the process of carrying out calculations similar to those presented here for all members of the original HEAT atomic and molecular dataset.

See the supplementary material for the atomic masses used in the above calculations, along with the details of the aug-cc-pCV7Z and aug-cc-pCV8Z basis sets.

This project would not have been possible if not for the highly efficient ncc module of Devin Matthews (SMU) in the cfour program package,^37,41–43 and the mrcc program, of Mihály Kállay (Budapest) and co-workers.⁴⁰ Additionally, we would like to thank Filippo Lipparini (Pisa) for discussions about large basis-set calculations and Lan Cheng (Johns Hopkins) for advice on relativistic treatments. J.H.T and J.F.S. were supported by the U.S. Department of Energy, Basic Energy Science, under Contract No. DE-SC0018164, and by the NSF under Award No. CHE-1664325. P.B.C. was supported by the NSF under Award No. AST-1908576. The work at Argonne National Laboratory was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences and Biosciences, under Contract No. DE-AC02-06CH11357, through the Computational Chemical Sciences Program (D.H.B.) and the Gas-Phase Chemical Physics Program (B.R.). Most of the authors are members of Active Thermochemical Tables Task Force One.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.