Ab Initio Computations and Active Thermochemical Tables Hand in

Article pubs.acs.org/JPCA

Ab Initio Computations and Active Thermochemical Tables Hand in Hand: Heats of Formation of Core Combustion Species Stephen J. Klippenstein,* Lawrence B. Harding,* and Branko Ruscic* Chemical Sciences and Engineering Division, Argonne National Laboratory, Argonne, Illinois 60439, United States S Supporting Information *

ABSTRACT: The fidelity of combustion simulations is strongly dependent on the accuracy of the underlying thermochemical properties for the core combustion species that arise as intermediates and products in the chemical conversion of most fuels. High level theoretical evaluations are coupled with a wide-ranging implementation of the Active Thermochemical Tables (ATcT) approach to obtain well-validated high fidelity predictions for the 0 K heat of formation for a large set of core combustion species. In particular, high level ab initio electronic structure based predictions are obtained for a set of 348 C, N, O, and H containing species, which corresponds to essentially all core combustion species with 34 or fewer electrons. The theoretical analyses incorporate various high level corrections to base CCSD(T)/cc-pVnZ analyses (n = T or Q) using H2, CH4, H2O, and NH3 as references. Corrections for the complete-basis-set limit, higher-order excitations, anharmonic zeropoint energy, core−valence, relativistic, and diagonal Born−Oppenheimer effects are ordered in decreasing importance. Independent ATcT values are presented for a subset of 150 species. The accuracy of the theoretical predictions is explored through (i) examination of the magnitude of the various corrections, (ii) comparisons with other high level calculations, and (iii) through comparison with the ATcT values. The estimated 2σ uncertainties of the three methods devised here, ANL0, ANL0-F12, and ANL1, are in the range of ±1.0−1.5 kJ/mol for single-reference and moderately multireference species, for which the calculated higher order excitations are 5 kJ/mol or less. In addition to providing valuable references for combustion simulations, the subsequent inclusion of the current theoretical results into the ATcT thermochemical network is expected to significantly improve the thermochemical knowledge base for less-well studied species.

1. INTRODUCTION Recent years have seen enormous progress toward the predictive simulation of internal combustion engines (ICE), which is one of the ultimate goals of the combustion modeling community. Such simulations require the coupling of chemical models for the conversion of the fuel into combustion products with numerical treatments of the fluid dynamics of reacting flows. Until recently, global simulations of the combustion process required enormous simplifications in one or the other aspect of the problem. Consequently, these two components (chemistry and fluid dynamics) of combustion modeling have evolved as independent research efforts with little communication between them. Advances in computational algorithms (see, e.g., refs 1−6) and hardware now allow for simulations that employ meaningfully accurate treatments of both aspects of the problem.7−11 There is now considerable interest in the use of ICE simulations as engine design tools,12,13 but to be truly effective, the predictive accuracy of such simulations needs to continue to improve. On the chemical side, the models should describe not only the conversion of the fuel into oxidation products but also the formation of various pollutants and the heat released at various stages of the conversion. Thus, the global chemical models typically consist of thermochemical and transport properties for hundreds of species together with rate coefficients for the thousands of reactions that connect these species within the combustion environment. For consistency purposes, rate © 2017 American Chemical Society

coefficients are generally entered for only the forward direction, with thermochemistry and reversibility arguments used to determine rate coefficients for the reverse direction. The thermochemical properties are also important in describing the heat release arising from the combustion process and thus the temperature fields for the reacting flows. The fidelity of the simulations naturally depends on the accuracy of the parameters that make up the chemical model. Although mechanistic studies generally focus on improving the description of the rate coefficients, combustion observables also show strong sensitivity to the thermochemical parameters.14−18 Modeling efforts conclude that combustion mechanisms are hierarchical in nature,19 with the mechanisms for larger fuels built up from, and strongly dependent on, component mechanisms for smaller fuels. The mechanism for the H2/O2 system is generally taken as a base mechanism, with essentially all combustion processes showing strong sensitivity to its parameters. Similarly, mechanisms for CH2O, CH3OH, and CH4 are central to our understanding of hydrocarbon combustion. The next step in the hierarchy involves models for the oxidation of the C2 hydrocarbons C2H6, C2H4, and C2H2.20 Finally, adding in models for C3 hydrocarbons (C3H8, C3H6, C3H4) and oxygenated C2 hydrocarbons (C2H5OH, Received: June 16, 2017 Revised: July 27, 2017 Published: July 31, 2017 6580

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A

other theoretical results and with ATcT values, the latter being available for about 40% of the species studied here. These comparisons provide valuable data for estimating the accuracy of our theoretical methods. They also highlight cases where discrepancies indicate the need for further study. Importantly, the current ATcT values already include substantial theoretical information in addition to the experimental data. Indeed, for more than half of the species of relevance to the present work, more than half of the provenance comes from theoretical studies, with many of these having more than 90% of their provenance from theoretical studies. This prominent role of theory creates some difficulties for the definitive benchmarking of the present methods. For example, many of the discordant cases in the current study correlate with cases where prior, generally more limited, theoretical analyses have played a significant, if not dominant role, in the ATcT predictions. Nevertheless, appropriate restrictions in the comparisons can be used to provide fairly definitive indications of the accuracy of the present theoretical schemes. These indications match well with the expectations from the convergence of the various corrections and from the comparisons with prior high level theory. An important outcome of the present analysis will be the future inclusion of the present theoretical predictions, with their benchmarked uncertainties, in subsequent versions of the Active Thermochemical Tables. The schemes employed here are closely related to other literature schemes.35−39 We make only minor modifications based on our specific computational resources, our expectations and observations of the relative importance of certain factors, our desire for simplicity, and our desire to treat a large number of species. Notably, we are also generally interested in making kinetic predictions for the reactions of such core chemical species. Quantitative kinetic predictions generally require accurate evaluations of transition state properties. The particular schemes we implement here are also meant to be effective for such evaluations. The present ab initio schemes are described in section 2. This description is followed in section 3 with a brief review of the ATcT approach. The predictions for the heats of formations are then presented and discussed in section 4, while some concluding remarks are made in section 5.

CH3CHO, CH3OCH3) yields what might be called a core mechanism for combustion. The combustion of a larger fuel (e.g., isooctane, dodecane, butylbenzene, pentanol, methylinoleate, etc.) relies heavily on the components of this core mechanism. Correspondingly, it is particularly important to have a high accuracy description of the parameters in the core mechanism. The Active Thermochemical Tables (ATcT) approach21,22 provides a particularly effective approach for accurately determining the requisite thermochemical properties. This approach involves the formation of a thermochemical network that includes all of the available experimental and theoretical determinations, each with its inherent uncertainties. The ATcT thermochemical values are that set of values which is most consistent with the knowledge content of the thermochemical network. The ATcT approach has been implemented for a number of important combustion species.23−27 However, there are still many species for which the limitations in the available data result in thermochemical estimates with large uncertainties. Ab initio electronic structure theory provides an alternative method for predicting thermochemical properties. For many years, a goal of chemical accuracy, which was specified as an uncertainty of 1 kcal/mol (in terms of 95% confidence limits28), provided a common rallying point for such predictions. In recent years, there has been remarkable further progress in the accuracy of such ab initio predictions, and a number of schemes, such as the W4,29 HEAT,30 focal point,31 and Feller−Peterson−Dixon32 methods, now yield uncertainties of about 1 kJ/mol or less. These high level schemes share many elements in common: they generally consist of geometry optimizations and harmonic vibrational analyses at the CCSD(T) level, complete basis set extrapolations employing Dunning’s correlation consistent basis sets,33 and corrections accounting for the effects of higher order electronic excitations, vibrational anharmonicity, core−valence interactions, relativity, and electron−nuclear coupling. Unfortunately, such high level ab initio based schemes are very computationally intensive. Nevertheless, they are readily applied to many of the species of importance in the core combustion mechanisms for C0−C3 fuels. In this work, we apply high-level schemes to the prediction of the 0 K heats of formation for the ground electronic states of essentially all combustion relevant species (or more precisely C/N/O/H species) with 34 or fewer electrons. Our inclusion of nitrogen containing species should contribute to the accurate treatment of NOx formation, a common goal in applied combustion research. The extension of these predictions to the thermodynamic properties at temperatures of relevance to combustion is left for future work. Extensions of the present effort to other moderately larger combustion species, such as those containing up to 50 electrons, will also be considered in future work. The range of species considered here covers most molecules with four or fewer heavy atoms, as well as a few with five or six heavy atoms. The set of species studied correlates closely with those studied in the pioneering bond-additivity corrected thermochemical study of Goldsmith et al.34 Previous applications and benchmarking of subchemical accuracy ab initio thermochemical schemes,29,35−39 have generally been considerably more limited in scope. As part of this effort, we examine in detail the magnitudes of the various corrections and make extensive comparisons with

2. THEORETICAL METHODOLOGY The coupled cluster method provides a very effective means for obtaining high accuracy electronic structure predictions, at least for systems where multireference effects are not too great, as is the case for most minima on potential energy surfaces. Indeed, essentially all high level ab initio schemes include an estimate for the complete basis set (CBS) limit of the CCSD(T) method40 (coupled cluster theory with singles, doubles, and perturbative inclusion of triples contributions) as a central element of their analysis. The present schemes similarly focus on the quantitative evaluation of this CCSD(T)/CBS limit. Uncertainties in the vibrational zero-point energy (ZPE) have generally been a major component of the uncertainty in composite schemes for obtaining thermochemical predictions. Early composite approaches employed Hartree−Fock-based schemes.41,42 Subsequent inclusion of B3LYP density functional theory based geometries and harmonic vibrational frequencies43,44 (or alternatively MP2 ones45) yielded significant reductions in the error estimates. However, even with such methods, the errors in the predicted ZPE are still significant. The underlying errors in the vibrational frequencies, even for 6581

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A more recent functionals such as the M06-2X functional,46 also map into significant errors in the canonical partition functions, which affects the temperature dependence of the thermochemical properties and kinetic predictions. Recent efforts to obtain subchemical accuracy29−32 have turned to CCSD(T) based evaluations of the geometries and harmonic frequencies, which generally compare favorably with spectroscopy based determinations. Although such CCSD(T) based frequency evaluations are much more computationally intensive, distributed parallel evaluations make them feasible, at least for the modest sized systems of interest here. Alternatively, double-hybrid density functional methods, such as the B2PLYP-D3 method of Grimme,47 which is available in the Gaussian software,48 show promise for reducing the geometric and ZPE uncertainties at reduced cost.49 For example, we find in a preliminary exploration for a large subset of the species considered here that CCSD(T)-F12/cc-pVQZ-F12 energies are on average lower at B2PLYP-D3/cc-pVTZ geometries than at CCSD(T)/cc-pVTZ geometries. Although such density functional methods are not explored here, they may prove useful in extending the present schemes to larger system sizes. Vibrational anharmonicities also have an important effect on the ZPE and the uncertainties in such anharmonic effects are perhaps the largest remaining contribution to the overall uncertainty for high level methods.50−52 Spectroscopic perturbation theory (VPT2) and vibrational configuration interaction (CI) provide effective means for estimating the anharmonic contribution to the ZPE, particularly when implemented with accurate electronic structure methods.53 The CCSD(T) method is sometimes inaccurate when the wave function contains substantial multireference character. For many of these cases, higher order coupled cluster methods, such as the CCSDT(Q) method,54 which includes a perturbative treatment of quadruple excitations of the electrons from the reference configuration, provide an effective means for treating such shortcomings, at least for moderate multireference effects. The inclusion of such corrections, facilitated by the valuable MRCC code of Kállay,55,56 has played a leading role in recent efforts to achieve subchemical accuracy. Importantly, such higher-order corrections might be expected to be particularly important in transition state evaluations, where significant multireference effects are more commonplace. Other corrections for core−valence interactions, relativistic effects, and Born−Oppenheimer assumptions are also now regularly included in the highest accuracy calculations.29−32 Core−valence corrections, which are readily obtained from CCSD(T) evaluations where all electrons are treated as active, are the most important of these. Relativistic and Born− Oppenheimer corrections, which are also readily evaluated with current electronic structures codes, are generally of relatively minor importance for C/N/O/H combustion relevant species. The evaluation of the diagonal Born−Oppenheimer correction has a secondary usefulness in that it indicates where the effects of conical intersections may be significant. 2.1. Specific Schemes. We consider two primary schemes here, which, for ease of reference, we label as ANL0 and ANL1. The ANL1 scheme is designed to be of higher accuracy, employing larger basis sets in various aspects of the calculation. However, these larger basis sets limit the applicability of the method. Thus, the ANL1 scheme is applied to only a limited subset of the species, consisting of most species with 22 or fewer electrons and a few high symmetry tetratomic species

with more electrons. Meanwhile, the ANL0 scheme is readily applied to the full set of 348 species considered here. The full ANL0 and ANL1 schemes, which are described in detailed below, may be denoted as EANL0 = ECCSD(T)/CBS(a ′QZ,a ′5Z)//CCSD(T)/TZ + E 0,har CCSD(T)/TZ + (E 0,anh B3LYP/TZ − E 0,har B3LYP/TZ) + (ECCSDT(Q)/DZ − ECCSD(T)/DZ) + (ECCSD(T,Full)/CBS(cTZ,cQZ) − ECCSD(T)/CBS(cTZ,cQZ)) + ΔE DKH/CCSD(T)/acTZ + ΔE DBOC/HF/TZ + ESO

(1)

EANL1 = ECCSD(T)/CBS(a ′5Z,a ′6Z)//CCSD(T)/QZ + E 0,har CCSD(T)/CBS(TZ,QZ) + (E 0,anh B3LYP/TZ − E 0,har B3LYP/TZ) + (ECCSDT(Q)/TZ − ECCSD(T)/TZ) + (ECCSDTQ(P)/DZ − ECCSDT(Q)/DZ) + (ECCSD(T,Full)/CBS(cTZ,cQZ) − ECCSD(T)/CBS(cTZ,cQZ)) + ΔE DKH/CCSD(T)/acTZ + ΔE DBOC/HF/TZ + ESO

(2)

EANL0 − F12 = EANL0 + ECCSD(T) − F12/CBS(TF,QF)//CCSD(T)/TZ − ECCSD(T)/CBS(a ′QZ,a ′5Z)//CCSD(T)/TZ

(3)

The geometries and harmonic vibrational frequencies are first determined with CCSD(T) evaluations employing the ccpVTZ (TZ) basis33 for ANL0 and the cc-pVQZ (QZ) basis33 for ANL1. For ANL1, the harmonic zero-point energy (ZPE), E0,har, is obtained from a CBS limit of the TZ and QZ values. Reference CCSD(T)/CBS limit energies are then determined from simple extrapolations of explicit evaluations for the aug′cc-pVnZ basis set series (where the prime indicates that diffuse functions are included for only the s and p orbitals in C, N, and O and the s function in H). For the ANL0 and ANL1 approaches, the extrapolations employ explicit calculations for the aug′-cc-pVQZ (a′QZ) and aug′-cc-pV5Z (a′5Z) or the aug′-cc-pV5Z and aug′-cc-pV6Z (a′6Z) basis sets, respectively. A minor modification of the ANL0 scheme, labeled ANL0-F12, employs an alternative CCSD(T)/CBS limit obtained from CCSD(T)-F12b57 calculations for the cc-pVTZ-F12 (TZ-F) and cc-pVQZ-F12 (QZ-F) bases.58 For these F12 calculations we employ MOLPRO’s default auxiliary and complementary auxiliary basis sets. Many procedures have been presented for extrapolating to the complete basis set limit (see, e.g., refs 59−62). Here we take a particularly simple approach of extrapolating the total CCSD(T) energies based on two point formulas: ECBS = Ea′nZ + α (Ea′nZ − Ea′n‑1Z). Various procedures for choosing α have been presented in the literature, but these usually focus on separate extrapolations for various components of the energy such as the SCF, CCSD, and (T) contributions, and even then there does not appear to be a well-accepted standard procedure (or parameters). Here, we simply note that the root-meansquare deviation (RMSD) between the ANL1 predictions and 6582

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A the ATcT values is minimized for α near 1.1 for n = 6. Meanwhile, the RMSD between the ATcT values and ANL0F12 approach that replaces the a′5Z,a′6Z CCSD(T) extrapolation with a TZ-F,QZ-F CCSD(T)-F12 extrapolation is minimized for α near 0.5. Similarly, an extrapolation based on a′QZ and a′5Z CCSD(T) values yields a minimal RMSD for α near 0.75. These values, which are similar to what would be expected from prior studies, are reasonably reproduced by an 1/l3.7 extrapolation (where l is the maximum angular momentum quantum number in the basis set) and that extrapolation is used throughout this work. This expression correlates with α = (0.29, 0.53, 0.78, 1.04) for n = 3, 4, 5, and 6, respectively. Notably, the mean error in the a′5Z, a′6Z based results is only very weakly dependent on the extrapolation coefficient, while that for the TZ-F,QZ-F based results is quite strongly dependent, and the a′QZ, a′5Z based results are somewhere in-between. For example, varying the extrapolation coefficient (in the l based extrapolations) from 3 to 5 varies the mean errors for a restricted 79 point comparison set (cf. discussion below) with ATcT from −0.33 to −0.34 kJ/mol for ANL1, from −0.18 to −0.33 kJ/mol for ANL0, and from 0.07 to −0.37 for ANL0-F12. Meanwhile, the RMSD from ATcT are all more weakly dependent on the extrapolation coefficient, with all values in the range from 0.62 to 0.78 kJ/mol for this range of coefficients. Corrections for higher order excitations (HOE) were obtained from CCSDT(Q) calculations63 with the cc-pVDZ basis for ANL0 and the cc-pVTZ basis for ANL1. CCSDTQ(P)/cc-pVDZ corrections were also included for ANL1. HOE corrections are occasionally large (e.g., > 5 kJ/mol), and it is not clear that they are well converged with respect to basis set even for the ANL1 approach. Note that other high level schemes such as HEAT, W4, and FPD include separate corrections for intermediate CCSDT, CCSDTQ, ..., calculations. We prefer to focus solely on the perturbative corrections. Ideally, one would employ larger basis sets for these corrections, but even for the present basis sets such evaluations were generally the most CPU intensive part of the calculations. The ready availability of an effective parallel version of the MRCC code would be a great benefit to future high-level thermochemical kinetics work. It may also be worth exploring alternative approaches such as multireference coupled cluster approaches64 for specific problem cases. The remaining corrections were treated equivalently in the ANL0 and ANL1 approaches. The core electrons were generally treated as frozen in the coupled cluster evaluations. The correction for core−valence interactions was obtained from CCSD(T,full)/CBS calculations (where the core electrons are included as active) based on extrapolation of results for the cc-pcVTZ (cTZ) and cc-pcVQZ (cQZ) basis sets.65 Restriction to the cTZ basis set would have yielded significantly different results, with numerous variations outside the norm of the uncertainties. A relativistic correction (ΔEDKH) was obtained from the difference in the CCSD(T) energy with and without the Douglas−Kroll one-electron integrals for calculations with the aug-cc-pcVTZ-DK basis (acTZ-DK).66 Diagonal Born−Oppenheimer corrections (ΔEDBOC)67 were obtained at the HF/cc-pVTZ level. The implementation of VPT2 to obtain anharmonic ZPE corrections, E0,anh − E0,har, requires the evaluation of the cubic and quartic force field. Even for the moderate sized species of interest here, the evaluation of such force fields through numerical displacements can be prohibitively expensive due to

the large number of single point evaluations that are required. Thus, density functional theory based evaluations of these force fields are appealing, particularly since the presence of analytic second derivatives greatly decreases the required number of evaluations. Unfortunately, for many of the recent/more accurate density functionals (e.g., M06-2X, B2PLYP-D3) such DFT-based evaluations of the anharmonic effects appear to be numerically unstable. Thus, the anharmonic vibrational zeropoint energy (ZPE) corrections were calculated at the B3LYP/ cc-pVTZ level via VPT2,68 where such instabilities were largely absent. A recent paper explores the effectiveness of such VPT2 based ZPE evaluations via comparison with CCSD(T)-based diffusion Monte Carlo evaluations.52 The spin−orbit corrections, ESO, are obtained from experiment. For atoms, they correspond to the straightforward spin− orbit correction involving the weighted center of gravity of the experimental spin−orbit terms. For diatomics they are what is sometimes called “rotational ZPE”, since in these the spin− orbit term is intimately coupled to the rotational angular momentum, and shifts are obtained by considering the appropriate rovibronic Hamiltonian rather than the spin− orbit constant Ae alone. (This usage explains why, for example, CH has such a small and positive correction even though Ae is 28.05 cm−1;69 the positive shift of 0.13 cm−1 is the difference between the vibrational ZPE corresponding to N = 0 and no spin−orbit and the lowest actual rotational level, N = 0, J = 1/2, of 2Π1/2 using Hougen’s or Van Vleck’s Hamiltonian; similarly, OH has Ae = −139.21 cm−1,70 but the lowest level N = 1, J = 3 /2 of 2Π3/2 is 38.1 cm−1 below the vibrational non-spin−orbit ZPE;30 O2 is particularly cute26 since technically the lowest level of 3Σg− would be N = 0, J = 1 at −1.78 cm−1, but it is wiped out by nuclear spin statistics in 16O2, making N = 1, J = 0 the lowest level at −1.09 cm−1, a situation not dissimilar to the state of affairs related to the thermodynamically correct adiabatic ionization energy of methyl.71) The present RUCCSD(T) calculations generally employed RHF wave functions within the UCCSD(T) formalism implemented in MOLPRO72 while the UUCCSDT(Q) calculations employed UHF wave functions as required by MOLPRO’s implementation of Kállay’s MRCC code.56 The vibrational frequencies were determined with MOLPRO’s numerical differentiation procedures, but with tighter integration convergence parameters to improve the accuracy for low frequency modes. The DBOC were obtained with the CFOUR code of Stanton and Gauss.73 The density functional theory calculations were performed with G09.48 2.2. Exceptions. Near degeneracies in the electronic states, which often correlate with the need for multireference wave functions, occasionally cause problems for the implementation of the base ANL0 and ANL1 schemes. In particular, the RUCCSD(T) vibrational analyses, the B3LYP spectroscopic perturbation theory analyses and the DBOC analyses each occasionally produce values that fall far outside our expectations based on correlations from similar species. To resolve these problems, we occasionally employ UUCCSD(T) harmonic vibrational analyses, MP2/TZ spectroscopic perturbation theory or simple correlations for the anharmonic vibrational correction, and a bound for the maximum DBOC. The resolutions to these problems, as outlined below, provide a further prescription of the ANL methods. These exceptions, which are limited here to ∼30 species, become very commonplace (and are even more severe) for electronically excited species. Thus, the only electronic excited state 6583

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A considered in this study is the 1A1 state of CH2, which is not plagued by degeneracy issues and is included due to its central importance to combustion modeling. Geometries and Harmonic Vibrational Frequencies. Various difficulties were encountered in the harmonic vibrational analyses. For some of the resonantly stabilized molecules, the RUCCSD(T) vibrational analysis fails due to singularities in the energies arising from symmetry breaking or related wave function convergence issues. For these species (i.e., 3CHCCH, 4 CHCHCH, CH2CHCH2, CH2CHCHCH3, CH2C(OH)CH2, CH2CHCHOH, NCN, CHCNH, and c-CHN2) we instead employ a UUCCSD(T) analysis. For CH2CCH, where both RU and UU approaches appear appropriate, sample CCSD(T)/cc-pVDZ(d/s) calculations (where d/s denotes neglect of polarization functions for H atoms) indicate a difference of 0.4 kJ/mol between the RU and UU calculated ZPEs of 105.3 and 105.7 kJ/mol. Meanwhile, the UUCCSDT(Q) calculated ZPE is even lower at 105.0 kJ/mol. For CH2CHCH2 the UUCSDT(Q) calculated ZPE is 0.4 kJ/ mol below the UUCCSD(T) value, again for the cc-pVDZ(d/ s) basis. Thus, it appears that the use of the UU rather than RU formalism should introduce differences on the order of 0.5 kJ/ mol. For the CH 2 CHCHCH 3 , CH 2 C(OH)CH 2 , and CH2CHCHOH species the UU vibrational analysis was restricted to the cc-pVDZ basis, which likely introduces further shortcomings on the order of 1 kJ/mol. Similarly, for the 3 CHCCH, 4CHCHCH, and CHCNH species the UUCCSD(T) analysis was restricted to the cc-pVTZ basis, which has a more modest effect on their ANL1 results. For C2nH species, the Σ and Π electronic states are nearly degenerate, with the Σ state being lowest for n < 2 and the Π state being lowest for n > 2. For C4H, the Σ and Π states are very nearly degenerate, with the predicted ground state varying with method. We presume a Σ ground state as observed experimentally,74,75 and use the UUCCSD(T) method to perform the vibrational analysis. Apparently the near degeneracy of the Σ and Π states creates difficulties for the RUCCSD(T) vibrational analysis, leading to a degenerate pair of imaginary frequencies. A UUCCSD(T) analysis instead produces all real frequencies, and the corresponding ZPE is used for the ANL analyses. Note that we calculate the Π state in linear symmetry with ANL0 to be 0.9 kJ/mol lower than the Σ state. An ANL0 analysis for the Π state at a nonlinear geometry obtained by RUCCSD(T)/TZ optimizations suggests a further 13.3 kJ/mol lower energy for the nonlinear geometry. However, we are not able to treat the coupled competing effects of Renner−Teller and spin−orbit splitting and this inadequacy could correlate with significant errors in the ZPE. Thus, we simply note that our ANL0 prediction for the C4H heat of formation is subject to significant uncertainty. For C6H, we find that the Π state is lower, as expected, but we have difficulty calculating the CCSDT(Q) correction due to UHF wave function convergence problems arising from orbital degeneracies. Assuming this correction is the same as for the linear Π state of C4H, yields a C6H heat of formation of 1010.9 kJ/mol. Similar series behavior is expected for the C2n‑1N series. However, for this series the Σ state is still clearly the ground state for C5N, which is the last member of the series that we consider here. A few species have an imaginary frequency for high symmetry with the cc-pVTZ basis set whereas they do not in the CBS limit (i.e., C2v symmetry for OCHO, Cs symmetry for NH2CO, NH2CHO, and linear for C5N). Thus, for OCHO,

NH2CO, NH2CHO, we employ cc-pVQZ values, while for C5N we extrapolate the cc-pVDZ and cc-pVTZ values. For cCH2CH2CH2O in C2v symmetry, there is a 22 cm−1 imaginary frequency for both cc-pVDZ and cc-pVTZ bases, suggesting some minor lowering from ring puckering. We ignore that lowering here. For OC(O)O we obtain 2 imaginary frequencies at the CCSD(T) level. For this species, we employ Davidson corrected multireference configuration interaction frequencies obtained for a six electron−six orbital active space. Vibrational Anharmonicity. The anharmonic corrections were generally determined with second order spectroscopic perturbation theory using the B3LYP density functional and the cc-pVTZ basis set. However, in many instances this approach did not yield realistic anharmonic corrections, and alternative corrections were then incorporated. These “failures” are associated with the presence of either a negative fundamental (most commonly) or a large positive (>600 cm−1) anharmonic correction for at least one mode. For a number of species [i.e., CCH, CH 3 C(CH 2 ) 2 , OCH 2 OH, ON(O)O, CHNO, CH3OON, CHCCCCN] simply replacing the B3LYP/ccpVTZ evaluations with MP2/cc-pVTZ ones yielded reasonable anharmonic corrections.76 For the reference species CH4, H2O, and NH3, higher accuracy anharmonic values were obtained in the theoretical studies of Scwhenke and co-workers.77−79 Notably, these higher accuracy values differ from their B3LYP/cc-pVTZ counterparts by only 0.01, −0.03, and −0.07 kJ/mol, respectively. The corresponding differences from MP2/cc-pVTZ values are similarly small, i.e., −0.10, −0.05, and −0.06 kJ/mol, respectively. Thus, the present B3LYP and MP2 evaluations of the anharmonic corrections should generally be quite satisfactory. For a number of other species [i.e., 3CHCCH, CH2CHC, C 4 H, CH 3 CCCH 2 , CH 3 C(CH 2 ) 2 , CH 3 O, CH 3 CH 2 O, CCCHO, C(OH) 3 , CH 2 CHCH 2 O, c-CH 2 CHCH 2 O, CH3CH2CH2O, OCHO, CH2OCHO, CCCCCN, NCNCH, NH2O, NH2NO, NCO, CNO, and NH2C(O)NH2] the B3LYP/cc-pVTZ and MP2/cc-pVTZ methods both appear to yield unsatisfactory/incorrect anharmonic corrections. These “failures” are again associated with negative fundamentals or a large positive anharmonic correction. These problem cases generally correlate with O radicals, linear radicals, or resonantly stabilized radicals. As illustrated in Figure 1, the calculated anharmonic corrections for the other molecules are reasonably well reproduced (2σ deviation of 0.7 kJ/mol) by the expression −0.481nH − 0.146 (kJ/mol), where nH is the number of hydrogen atoms in the molecule. Only the linear species CHCCCCCH and CHCCCO fall more than 3σ above the line, while only the two O-centered radicals, CH3CH(O)CH3 and OHCH2CH2O fall more than 3σ below the line. This expression is used to predict the anharmonic corrections for the problem species where both B3LYP and MP2 VPT2 calculations fail. Diagonal Born−Oppenheimer Correction. The magnitude of the DBOC contribution to the heat of formation (calculated with respect to the values for the reference species) is generally less than 0.8 kJ/mol. For a few species, the calculated DBOC contribution appears to be much larger, with these species tending to be ones with degenerate or nearly degenerate electronic states. Further inspection for a few such cases indicates that their calculated DBOC tends to be very sensitive to the electronic structure method used, with CCSD based calculations often yielding very different, much lower values. 6584

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A Δf H °0,ANL (CaOb NcHd)

= [aΔf H °0 (CH4) + bΔf H °0 (H 2O) + c Δf H °0 (NH3) + (d /2 − 2a − b − 3c /2)Δf H °0 (H 2)] − [EANL(CaOb NcHd) − aEANL(CH4) − bEANL(H 2O) − cEANL(NH3) + (2a + b + 3c /2 − d /2)EANL(H 2)] (5)

H2 is, of course, the thermodynamic reference for hydrogen with a defined zero enthalpy of formation, and the enthalpies of formation for the three hydrides are well-known from ATcT analyses,21−27 with uncertainties of only ±0.056, ± 0.027, and ±0.030 kJ/mol, respectively. Furthermore, this set is particularly appropriate as a computational reference due to the dominantly single reference nature of the ground electronic states for these molecules. Limited considerations of other possible references suggested somewhat increased errors for the overall comparison. For example, the alternative choice of atomic references yields roughly a factor of 2 increase in the root-mean-square error, in spite of the fact that from ATcT analyses the enthalpies of formation of H, C, O, and N are currently known24,26,27 to ±0.00001, ± 0.050, ± 0.0021, and ±0.023 kJ/mol, respectively. We also report here predicted heats of formation from CBSQB3, G4, BAC-QCISD(T)/CBS, and W4 calculations. For the first three of these methods, the optimization process inherent to the methods presumes atomic references. Meanwhile, for the W4 method, the original presentations focus on atomization energies. Thus, for each of these cases, we choose to use atomic reference species.

Figure 1. Plot of the correlation between the VPT2 anharmonicity and the number of H atoms in the species. The points denote the VPT2 data, while the dashed line denotes the best fit to the data.

The species with an apparent contribution greater than 1.3 kJ/ mol are c-CHCHCH (1.84), CH 2 CHCH 2 O (2.26), CH3CH2CH2O (1.63), CHOO (1.80), OCCCO (1.46), OC(O)O (2.01), CH2CN (1.30), and NO (1.34), where the numbers in parentheses denote the calculated contribution to the heat of formation in kJ/mol. For these species the DBOC contribution was reset to 1.3 kcal/mol, i.e., the upper end of the smooth distribution of values. For some other species, [C3H in linear and ring forms, CH 2 CCCH, C 5 H, 3 CHCHO, CH 3 CH 2 O, CCCHO, CHCHCHO, CH 3 C(OH)CH 3 , CH 2 CH 2 CH 2 OH, 3 CH 2 OO, CHCHOOH, CH 2 CCN, NC HC CH, CH CN CH, CH 2 CH CH N, C HC HN O, CH3CHNO] SCF convergence problems were encountered. For these species the DBOC contribution to the heat of formation was set to zero. For CH3CN an aug-cc-pVDZ basis was used since SCF convergence problems were encountered with the cc-pVTZ basis. Finally, for the C4H9 species, the disk requirements for the HF/cc-pVTZ DBOC calculations exceeded our resources and the calculations were instead performed with a cc-pVDZ basis. 2.3. Reference Species. The theoretical prediction of heats of formation requires the specification of reference values for a set of species containing the elements of interests. Most theoretical methods employ atomic species as references, due, at least in part, to historical reasons dating back to the introduction of composite electronic structure methods.80 However, the character of the wave function for atomic species is quite different from that for a typical bound species, which can magnify any shortcomings in the theoretical methods. Thus, the use of molecular species as a reference can yield a considerable reduction in the average error of the predictions. Here we consider H2, CH4, H2O, and NH3 as reference species for the H, C, O, and N elements, respectively. In particular, we consider a “working reaction” for each CaObNcHd species

3. ACTIVE THERMOCHEMICAL TABLES The development of new electronic structure methods generally benefits from the availability of a set of accurate and reliable benchmarks. So far, the most accurate set of thermochemical quantities has been obtained by the Active Thermochemical Tables (ATcT) approach, which has essentially set a new standard of accuracy in thermochemistry by systematically lowering thermochemical uncertainties, sometimes by one or more orders of magnitude. The ATcT approach is a new paradigm for obtaining accurate, reliable, and internally consistent thermochemical quantities.21,22 In order to derive thermochemical quantities such as enthalpies of formation, traditional thermochemistry relies on a sequential approach (A begets B, B begets C, etc.). Such an approach focuses on one species at time, examines the available data, and selects the “best” available determination that leads to the definition of the thermochemistry of the target species, taking the enthalpies of formation of the species entertained in previous steps as already known and fixed. While relatively straightforward, the sequential approach suffers from a number of problems, such as incomplete utilization of the existent knowledge, the inherent cumulative error entailed by a method that uses a sequence of discrete steps where the value determined during one step is fixed and kept constant in subsequent steps, the inability to a posteriori propagate the consequences of new determinations because of the hidden progenitor-progeny relationships within the derived set of thermochemical values, etc. These disadvantages are sidestepped by the ATcT approach, which relies on constructing and solving a thermochemical network (TN) that contains all

aCH4 + bH 2O + c NH3 + d /2H 2 → CaOb NcHd + (2a + b + 3c /2)H 2

(4)

with the computed enthalpy of formation given by 6585

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A Table 1. Comparison of Calculated ΔfHo(0 K) (kJ/mol) with ATcT Valuesa species

ATcTb

unc.c

ΔANL1d

H H2 C CH 3 CH2 1 CH2 CH3 CH4 C2 CCH CHCH CH2C CH2CH CH2CH2 3 CH3CH CH3CH2 CH3CH3 c-CCCH c-CHCCH 3 CHCCH CH2CC CH2CCH c-CHCHCH CH3CC c-CH2CHC CH3CCH CH2CCH2 c-CHCH2CH CH2CHCH2 CH3CCH2 CH3CHCH c-CH2CHCH2 CH3CHCH2 c-CH2CH2CH2 CH3CHCH3 CH3CH2CH2 CH3CH2CH3 CHCCCH CH2CHCHCH2 CH3CCCH3 c-CH2CH2CHCH CH3CH2CCH (CH3)2CCH2 CH3CHCHCH3 CH2CHCH2CH3 c-CH2CH2CH2CH2 CH3CH2CH2CH2 (CH3)2CHCH3 CH3CH2CH2CH3

216.0 0 711.4 592.8 391.0 428.6 149.8 −66.6 820.2 563.9 228.9 411.3 301.1 61.0 361.5 130.9 −68.3 710.7 497.1 543.4 554.4 354.0 491.0 529.0 528.1 192.8 197.3 292.5 179.6 262.9 278.0 302.8 35.0 70.8 105.3 117.9 −82.7 458.6 125.1 159.0 176.8 179.4 3.9 9.3 20.9 52.2 101.5 −106.2 −98.6

0.0 0 0.1 0.1 0.1 0.1 0.1 0.1 0.3 0.1 0.1 0.3 0.3 0.1 0.8 0.3 0.1 1.1 0.5 0.6 0.4 0.4 0.8 0.8 1.3 0.2 0.3 0.5 0.5 0.7 0.8 1.6 0.2 0.5 0.7 0.6 0.2 0.9 0.4 0.6 0.9 0.7 0.4 0.4 0.4 0.5 1.0 0.3 0.3

0.1 ref −0.6 −0.1 −0.3 0.1 −0.0 ref −0.5 −0.8 −0.3 −0.5 −0.6 −0.8 −1.1 0.1 −0.7 0.4 −0.9 −3.9 −1.5 0.1 −0.4 −0.8 −1.2 −0.4 −0.7 −0.3

O OH H2O CO HCO COH CH2O CHOH CH2OH

246.8 37.3 −238.9 −113.8 41.4 217.3 −105.3 112.4 −10.3

0.0 0.0 0.0 0.0 0.1 0.7 0.1 0.3 0.3

−0.8 −0.1 ref −0.2 −0.5 0.4 −0.4 0.1 −0.6

ΔANL0e

ΔANL0-F12f

ΔW4g

hydrocarbons: CxHy 0.1 0.0 ref ref ref −0.1 0.0 0.5 ref 0.3 0.7 −0.2 0.1 0.4 −0.2 0.5 0.8 −0.1 0.1 0.3 −0.2 ref ref −0.2 −2.2 −1.4 0.7 −1.1 −0.8 −0.7 −0.5 −0.5 −0.4 0.2 0.4 −0.3 −0.2 −0.1 −0.3 −0.8 −0.7 −0.4 −0.4 −0.2 0.4 0.5 −0.6 −0.6 −0.4 0.9 1.3 −0.5 −0.3 −2.0 −2.1 −0.8 −0.8 0.0 −0.3 −0.1 0.0 −0.8 −0.7 −1.2 −0.8 −0.5 −0.7 0.3 −0.4 −0.5 0.2 −0.2 −0.1 0.0 0.0 0.1 0.2 0.4 0.5 1.6 1.8 −0.5 −0.5 −0.2 0.7 0.7 −0.2 −0.2 0.3 0.3 −0.5 −0.5 −0.3 −0.6 −1.0 −0.1 −0.2 0.2 −0.3 1.0 1.0 0.8 0.5 0.5 0.4 0.3 0.2 0.4 0.4 −0.6 −0.7 1.7 1.7 0.3 0.2 0.2 0.1 oxygenated hydrocarbons: CxOyHz 0.3 1.2 ref 0.5 0.9 −0.4 ref ref −0.1 −0.1 −0.1 0.2 −0.7 −0.7 0.4 0.6 0.8 −0.3 −0.5 −0.2 0 0.5 0.7 −0.3 −0.2 6586

ΔW3-XLh

ΔBACi

ΔAvej

ΔG4k

ΔCBS-QB3l

ref.

ref

ref. −0.9 −0.3 −0.4 −0.6 −0.2

ref −0.2 ref −0.7 1.0 −0.3 0.6 0.7

ref −2.4 −0.2 −1.3 −1.9 −1.4

ref −1.6 ref −4.5 −1.5 −3.2 −2.0 −0.1

ref −4.9 ref −0.0 4.7 −0.1 2.1 0.2

−0.2 1.6 0.3 −0.4 −0.2 −1.9 −0.8 −1.4

−0.7 0.5 −0.1 0.7 −0.3 0.8 1.1 0

2.8 3.3 2.3 0.1 0.3 0.6 0.1 −1.2

−5.0 0.4 −1.3 −2.6 −0.1 −0.8 −0.5 1.0 −2.9 −1.3 −2.3 −2.9 −2.4 −1.5 −5.8 −3.7 0.0 −0.2 1.1 −0.7 −3.2 −2.9 −0.6 1.1 1.2 −0.4 −0.8 2.0 0.1 1.6 1.7 3.6 2.1 1.7 2.2 2.3 1.8 −0.5 2.3 3.1

4.5 6.0 5.7 3.4 3.2 6.0 4.6 1.9 7.2 6.9 1.3 5.6 4.8 7.8 4.4 6.1 5.9 4.9 6.4 4.2 4.7 5.8 8.5 5.4 6.0 5.8 6.5 3.8 9.9 8.9 7.2 10.1 9.1 7.4 7.5 7.7 6.4 9.3 5.7 6.1

ref −1.2 1.8 −3.8 −2.7 −3.0 −2.7 −1.5 −0.2

ref −0.4 −1.7 −1.8 −1.2 −0.6 −4.9 −1.3 0.2

−2.7 −2.5 1.0 −1.6 −1.0 0.0 −0.2 −0.1

−0.8 −0.7

−1.9 1.2

−0.3 −1.2 −1.0 2.0 0.5 0.6 2.2 −0.8 0.8 0.9 1.2 −0.6 −0.1 0.6

−0.4 2.0 3.5 1.8 1.4 0.8 2.6

0.8 1.2

−0.8 4.2

0.0 0.0 0.7 −0.8 −0.9 −0.5

0.6 −0.3 4.5 0.0

1.7 4.6 3.0 1.2 1.3 1.4

−1.4 −1.8

2.2 0.2 0.1

−0.5 −0.1

−0.5 0.3 −0.4 −0.9 −0.6 −0.8 −0.4

ref −0.1 −0.9 0.4 0.7 −0.3 0.6 1.3

ref −0.8 −1.6 −1.6 −0.4 −2.8 −0.8 −0.1

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A Table 1. continued species CH3O CH3OH 3 CCO HCCO CH2CO CHCOH c-CHOCH CH3CO CH2CHO c-CHOCH2 CH3CHO CH2CHOH c-CH2CH2O CH3CHOH CH2CH2OH CH3CH2O CH3CH2OH CH3OCH3 CH3C(O)CH3 CH3CH2CHO (CH3)2CH(OH) CH3CH2CH2OH O2 HO2 OHOH OCO OCOH OCHO OCHOH c-CH2OO CH3OO OHCH2OH CH3OOH OCHCHO CH3OCO CH2OCHO CH3C(O)OH CH3OCHO O3 OOOH OHOOH OHOOOH OC(OH)2 OCHOOH N 3 NH NH2 NH3 N2 NNH NHNH NH2N NHNH2 NH2NH2 CN HCN CNH CH2N

ATcTb

unc.c

ΔANL1d

28.9 −189.8 376.3 176.8 −45.5 94.9 272.8 −3.3 22.5 172.8 −155.0 −112.6 −40.0 −42.2 −12.8 1.1 −216.9 −166.5 −199.5 −170.0 −248.8 −231.4 0 15.2 −129.5 −393.1 −181.0 −124.8 −371.1 9.4 22.3 −379.1 −114.9 −207.0 −149.8 −148.7 −418.4 −344.8 144.4 25.3 −81.4 −33.5 −602.8 −279.8

0.3 0.2 1.0 0.6 0.1 1.3 1.0 0.4 0.8 1.3 0.3 0.8 0.4 0.6 0.6 0.5 0.2 0.4 0.4 0.9 0.4 0.3 0 0.2 0.1 0.0 0.5 0.5 0.2 0.5 0.5 0.9 0.7 0.5 1.9 1.5 0.5 0.6 0.0 0.1 0.7 2.3 0.8 1.1

−1.0 −0.3 0.5 −0.2 −0.2 −0.3 −1.5

470.6 358.7 188.9 −38.6 0.00 252.2 207.3 308.0 235.3 111.8 436.7 129.7 192.0 242.2

0.0 0.2 0.1 0.0 0.0 0.5 0.4 0.7 0.8 0.5 0.1 0.1 0.4 0.8

−0.5 −0.7 0.4 0.0 −0.4 −4.0

0.1 1.8 −0.5

−0.9 −0.2 0.0 ref −0.1 −0.5 −0.2 −0.5 −1.0 0.0 −0.7 −0.6 −0.3 −0.0

ΔANL0e

ΔANL0-F12f

ΔW4g

oxygenated hydrocarbons: CxOyHz 0.6 0.7 −0.4 −0.4 −0.1 0.6 0.7 −0.2 −0.3 −0.1 −0.4 0.4 −0.5 −0.7 −2.4 −2.1 −0.7 −1.2 −1.2 0.8 0.6 0.4 0.6 0.1 −0.2 0.2 −1.9 −1.9 0.1 0.1 0.5 −0.9 −0.9 0.7 0.7 0.2 0.4 −0.5 −0.6 0.1 0.3 0.1 −0.7 −1.0 −0.5 −0.8 0.9 0.8 0.9 0.7 −0.5 0.7 0.5 −0.3 0.3 0.1 0.2 0.5 0.6 −0.3 −0.8 0.2 −0.6 −0.8 −4.8 −5.0 0.0 −0.3 0.8 −0.6 −0.3 0.4 0.7 0.4 0.3 0.2 0.3 −0.9 −1.2 −0.7 −1.6 −1.9 −0.7 −1.2 0.8 0.4 −0.6 −1.1 0.9 −0.8 0.6 1.0 2.0 2.8 6.0 −0.9 −0.3 −1.3 −0.4 0.6 0.1 0.6 0.5 nitrogen compounds: CxNyOzHw −0.2 0.5 ref 0.3 0.7 0.0 0.2 0.3 −0.0 ref ref −0.1 −0.3 −0.2 −0.3 −0.5 −0.4 −0.7 −0.5 −0.3 −0.1 −0.8 −0.9 −0.9 −0.9 −0.7 −0.6 −0.1 −1.3 −1.2 −0.5 −0.7 −0.7 −0.0 −0.2 −0.4 −0.1 0.6 0.6 −0.3

6587

ΔW3-XLh

ΔBACi

ΔAvej

ΔG4k

ΔCBS-QB3l

−2.7 −1.9

−0.8 −1.1 0.6 −0.5 −0.2 −0.3 1.6 0.5 3.4 2.8 −0.1 −2.0 0.8 0.4 0.0 −0.8 −1.6 0.1 −1.2 −0.8 −1.3 −1.3 0.0 −0.8 1.2 −0.3 1.0 0.2 −0.2

−2.6 −1.8

−3.5 −0.0 −1.5 −2.7 −0.0 1.0 −2.5 −3.3 −2.2 −2.9 −0.5 0.7 −1.3 −1.7 −1.8 −2.4 1.0 −0.3 −0.1 −0.2 0.8 1.1 1.0 0.3 3.2 −3.6 −2.9

−3.2 −3.3 2.1 0.8 −0.1 1.0 −3.1 −0.4 −2.0 −0.1 −1.9 −0.0 −3.6 0.0 1.1 −2.4 −1.6 −4.9 −1.1 −1.2 2.6 3.1 −3.4 −3.9 −5.2 −7.2 −4.9

−0.2 −2.6 −1.9 1.0 0.8 −3.6 −6.7 −6.5 1.6 −2.6 1.6 7.4 2.8 2.8 1.2 2.1

−5.9 −9.0 −5.5 −6.1 −8.1 −7.2 −7.0 −7.7 −4.3 −9.4 0.1 −1.2 −10.5 −15.1 −7.9 −8.8

ref −3.4 −0.2 2.8 −1.9 −0.0 −0.5 −2.0 1.1 3.0 2.3 −1.0 −1.5 −2.8

ref 2.0 2.2 1.4 2.9 −0.7 0.9 −2.1 1.2 2.0 5.8 3.4 −0.1 0.4

−1.4 −0.4 0.0 −1.7

−1.1 −2.3 −0.7

−3.9 −1.1 −2.7

−2.2 0.9 −1.0 −0.6 −1.3 −2.4 −0.8 −1.8

−0.9 −2.4

1.3 −1.0 −0.3 0.6 0.8 0.4 −0.5

−11.4

−1.5 −0.9 1.3 −0.3 −2.4

−1.8 −0.5 −2.2 −1.6 −1.0 −2.8 −3.4 −2.6 −2.4

−2.8 6.0 −2.0 −3.7 −3.0 −1.8 −2.0 −4.0

−1.8 −5.7 −11.1

−1.3 ref

ref

0.0 0.4

ref −1.5 1.0 1.8

1.0 1.5 −0.3 −0.4 0.2 −0.1 1.8 0.7

2.0 2.6 1.0 2.3 3.5 4.8 2.2 0.2

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A Table 1. continued species

ATcTb

unc.c

ΔANL1d

CHNH CH2NH 3 CH3N CH2NH2 CH3NH CH3NH2 CH3NHCH3 3 NCN NCNH CHNN NCCN CNCN NO HNO 3 NOH NH2O NHOH NH2OH NNO ONO ONOH ONHO ON(O)O ON(OH)O ONOOH NCO CNO NHCO NCOH CHNO CNOH CH2N(O)O CH3N(O)O CH3ONO

275.9 96.5 323.8 159.4 187.9 −6.0 5.9 450.4 324.7 472.8 308.2 410.7 90.6 109.9 217.5 71.3 101.5 −33.1 86.0 36.9 −73.0 −37.2 79.4 −124.5 −6.0 127.0 389.4 −115.9 −12.2 170.9 235.4 138.8 −60.9 −55.5

0.8 0.6 1.6 0.6 0.6 0.5 0.7 0.7 1.7 1.9 0.4 1.6 0.1 0.1 0.9 0.8 1.0 0.5 0.1 0.1 0.1 1.4 0.2 0.2 0.4 0.3 1.4 0.3 0.5 0.5 0.5 1.2 0.5 0.5

−0.2 −0.2 −0.9 0.1 1.1 −0.0 −0.3 −3.3 −7.1 −1.3 0.5 −0.1 0.1 −0.8 −0.6 0.1 0.6 −0.6 0.6 1.3

−0.4 0.1 0.2 −0.2 −0.7 −0.2

ΔANL0e

ΔANL0-F12f

ΔW4g

nitrogen compounds: CxNyOzHw −0.4 −0.4 −0.3 −0.3 −0.4 −0.3 −0.1 0.2 0.4 0.4 −0.4 1.4 1.4 −1.2 −0.3 −0.3 −0.5 −0.4 −0.5 0.1 0.2 −3.2 −3.4 −7.3 −7.4 −1.6 −1.9 −0.7 1.8 1.3 0.6 0.9 0.0 −0.2 0.1 0.2 0.4 0.7 −0.2 −0.1 −0.7 −0.5 −0.2 −0.2 −0.6 −0.8 −0.2 −0.2 0.2 0.3 −0.3 0.1 −0.8 0.7 1.0 −2.1 −1.0 0.4 0.8 −1.3 −0.6 −0.6 −0.8 0.4 0.4 −0.4 −0.8 0.4 −0.5 −0.7 0.3 −1.1 −1.3 0.1 −0.3 −0.4 0.2 2.6 2.8 0.3 0.4 0.9 1.0

ΔW3-XLh

ΔBACi

ΔAvej

ΔG4k

ΔCBS-QB3l

0.1 0.4

0.6 0.6

0.4 0.2 −0.6 −1.4

2.6 1.8 0.5 −0.8

1.8 3.8

−0.4 −1.1

0.7

−0.3

−0.6 0.1

1.1 −0.5 1.1

0.3

−0.2

−1.6 −0.7 −0.4

−4.6 −3.5

−1.1 0.2 0.7 0.2

−1.6 −0.4 −3.1 −2.9

−2.6 −1.2 −3.3 1.6 −0.4 1.4 0.3 −1.1 −5.7 −0.59 −5.6 −5.1 −1.8 −2.1 −0.6 −0.5 −0.3 2.0 −3.4 −4.8 −0.4 −2.5 −4.8 −3.2 2.7 −4.1 −4.7 −0.8 −1.5 −3.3 −4.7 −5.1 −4.3 −0.0

1.8 1.0 2.7 4.4 3.7 1.1 0.5 −2.1 −3.8 −8.1 2.4 0.3 −4.1 −4.3 −1.5 −5.6 −2.4 −3.3 −6.2 −9.2 −5.3 −7.6 −16.6 −13.3 −8.1 −4.1 −4.4 −2.9 −2.3 −4.4 −5.2 −8.1 −9.9 −6.9

a Zero K heats of formation in kJ/mol. bATcT values for ΔfHo(0 K). cUncertainty in the ATcT values. dDifference between ANL1 and ATcT values for ΔfHo(0 K). eDifference between ANL0 and ATcT values for ΔfHo(0 K). fDifference between ANL0-F12 and ATcT values for ΔfHo(0 K). g Difference between W4 and ATcT values from Karton et al.85 for ΔfHo(0 K). hDifference between BAC-QCISD(T)/CBS from Goldsmith et al.34 and ATcT values for ΔfHo(0 K). iDifference between averaged values from Simmie and co-workers89−92 and ATcT values for ΔfHo(0 K). jDifference between G4 and ATcT values for ΔfHo(0 K). kDifference between CBS-QB3 and ATcT values for ΔfHo(0 K).

very distributed provenances,26,27 in contrast to sequentially derived thermochemical values, whose provenances most frequently involve the thermochemistry of only one chemical reaction per target species, combined with a set of assumed (or previously determined) thermochemical values for the other species in such reaction. The ATcT results used in this study are based on TN ver. 1.122b, which contains ∼1200 chemical species interlinked by more than 20,000 thermochemical determinations. Its immediate predecessor, TN ver. 1.122, was used to determine all the sequential bond dissociation enthalpies of methane, ethane, and methanol.27 TN ver. 1.122b is nearly the same, and has been obtained by slightly expanding TN ver. 1.122 in order to derive the definitive value for the energy difference between HCN and its isomer HNC.83 The full set of 150 ATcT 0 K enthalpies of formation used in the present study is given in Table 1 in the second (the values) and third (the accompanying uncertainties) columns. The accompanying uncertainties correspond to 95% confidence limits, in keeping with the generally accepted standard for

available thermochemical determinations (both from experiment and theory) relevant for defining the thermochemistry of the included chemical species, such as reaction enthalpies or free energies, bond dissociation energies, adiabatic ionization energies, dissociative photoionization thresholds, adiabatic electron affinities, relative gas-phase acidities, etc.81,82 These determinations, together with their uncertainties, effectively present a set of conditional constraints that the resulting thermochemistry should simultaneously satisfy. Since inconsistent determinations (essentially determinations with too optimiztic uncertainties) will tend to skew the final results, prior to attempting a solution ATcT statistically analyzes the TN in order to locate such determinations and iteratively augments their uncertainties until the epistemic content of the TN is self-consistent. Once self-consistency is achieved, the TN is solved simultaneously for all chemical species. 23−27 Compared to traditional thermochemistry, the ATcT solutions have superior accuracy because they are based on the entire knowledge content of the TN. The ATcT solutions also exhibit a much larger degree of robustness, since they typically have 6588

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A expressing uncertainties in thermochemistry.28 The selection of species was guided by the desire to present here the ATcT values for all the species in common between the set of core combustion species computed during this study and the set of species existing in ATcT TN ver. 1.122b. This approach to selection is rather different than that used in previous ATcT benchmarkings of highly accurate electronic structure methods (such as W4,29 HEAT,30 FPD,84 etc.), where the process of crafting a benchmark set took also into account the accuracy of the selected ATcT enthalpies of formation. In contrast, the full set used here includes species for which there is very little in terms of prior experimental and/or theoretical determinations. Consequently, the individual uncertainties ui of the enthalpies of formation in the complete ATcT set presented here cover a rather broad span, ranging from less than ±0.001 kJ/mol all the way to ±1.9 kJ/mol, with a median of ±0.46 kJ/mol. Using the whole set for benchmarking (as opposed to, for example, using only those species that have an uncertainty lower than a particular cutoff, say ±0.60 kJ/mol) presents some interesting challenges, as discussed below, as well as later during the analysis of the results. Namely, as pointed out previously,28 the apparent accuracy of a method obtained by benchmarking (i.e., type A uncertainty), uA, is a composite of the inherent accuracy of the method, umethod, and the overall accuracy of the benchmark set, ubench: uA ≈

umethod 2 +

ubench 2

Figure 2. Dependence of the estimated average uncertainty of the ATcT benchmark set, ubench, on the maximum uncertainty umax of the included benchmarks.

predictions via examinations of (i) the size and convergence of the corrections, (ii) the distribution of energy differences between W4 and ANL results, and (iii) the distribution of energy differences between ANL and ATcT results. The ANLtheory and the ATcT databases are also used to explore the accuracy of other more limited theoretical methodologies. Each of these explorations provides some perspective on the accuracy of the present and related theoretical calculations. 4.1. Corrections and Basis Set Convergence. The cumulative distribution functions for the various corrections to the computed ΔfHo(0 K) in ANL0 are illustrated in Figure 3. Recall that these corrections are computed relative to the corresponding corrections for the reference species (cf. section 2.3). The HOE [T(Q)/cc-pVDZ], core−valence, and vibrational anharmonicity corrections are of similar magnitude, with mean corrections of −1.8, 2.1, and 2.3 kJ/mol and RMSD from this mean of 2.4, 1.2, and 1.2 kJ/mol. Although not shown in the plot, for the HOE correction there is a tail to the distribution that significantly affects the average and RSMD values. In particular, the OC(O)O, CH3OON, O3, ON(O)O, CH3OOO, and OOOH species have HOE corrections of −24.1, −15.7, −13.6, 12.9, −12.4, and −9.9 kJ/mol, respectively. Simply removing those 6 species reduces the mean value and RMSD from the mean for the T(Q)/cc-pVDZ correction to −1.5 and 1.6 kJ/mol, respectively. The relativistic and diagonal Born−Oppenheimer corrections are considerably smaller in magnitude with mean corrections of −0.63, and 0.23 kJ/mol and RMSD relative to these means of 0.51 and 0.37 kJ/ mol, respectively. The large magnitude of the HOE, core− valence, and vibrational anharmonicity corrections is indicative of the difficulty of obtaining overall uncertainties much below 4 kJ/mol using methods that ignore any of these correction terms. Distribution functions for various combinations of these corrections are illustrated in Figure 4. The mean values for the total, HOE + core−valence + anharmonicity, core−valence + anharmonicity, and relativistic + DBOC corrections are 2.4, 3.4, 4.4, and −0.41 kJ/mol, respectively, while the corresponding RMSD from the means are 2.1, 2.1, 1.7, and 0.47 kJ/mol. The large negative HOE corrections for a few molecules again results in large RMSD values for the total and HOE + core− valence + anharmonicity corrections, while 95% of the values actually fall within 2.4 kJ/mol of the mean. The relatively low

(6)

where ubench can be estimated by averaging in quadrature the uncertainties of the individual benchmarks in the set ubench 2 =

1 n

n

∑ ui 2 i=1

(7)

Clearly, uA ≈ umethod only if ubench is sufficiently small. In other words, if umethod and ubench happen to be comparable, then the apparent u A will be larger than the intrinsic u method, approximately by a factor of 1.4. In practice, considering the fact that generic method uncertainties obtained by benchmarking are seldom more accurate than ±10%, in most cases it is sufficient that ubench < 0.5 umethod, in which case uA ≈ 1.1 umethod. The dependence of the average uncertainty of the benchmark set, ubench, on the maximum allowed uncertainty umax within the current ATcT benchmark set is depicted in Figure 2. For the complete set of 150 species, ubench = 0.64 kJ/mol; for a set pruned in such a way that no uncertainty is larger than umax = ± 1.0 kJ/mol (137 species), ubench = 0.50 kJ/mol, if pruned so that umax = ± 0.6 kJ/mol (99 species), ubench = 0.33 kJ/mol, and for a further reduced set with a maximum uncertainty umax = ± 0.5 kJ/mol (87 species), ubench = 0.28 kJ/mol.

4. RESULTS AND DISCUSSION The present ANL0, ANL0-F12, and ANL1 predictions for the 0 K heat of formation, ΔfHo(0 K), are compared with the corresponding ATcT determinations in Table 1. Also reported therein are values from W4,29,85 W3X-L,86 BAC-QCISD(T)/ CBS,34 G4,87 and CBS-QB388 theoretical methods, and an average of CBS-QB3, CBS-APNO, G3, and G4 (as promulgated by Simmie89−92). For molecules not in the ATcT data set, the ANL0, ANL0-F12, and ANL1 predictions for ΔfHo(0 K) are reported in Table 2. In the following discussion, we explore the expected accuracy of these 6589

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A Table 2. ΔfHo(0 K) (kJ/mol) for Species not in the ATcT Data Seta species CH2CHC 4 CHCHCH CH3CCH3 3 CH3CH2CH CCCCH CH2CCCH CHCHCCH c-CHCHCHC CH2CHCCH CH2CCCH2 c-C(CH2)CHCH c-CHCHCHCH CH3CCCH2 CH2CHCCH2 CH3CHCCH c-CH2CHCHCH CH2CHCHCH CH2CCHCH3 CH2CHCHCH3 CH3C(CH2)2 CH2CHCH2CH2 CH3CCHCH3 c-CH2CH2CH2CH (CH3)2CCH CH3CH2CHCH (CH3)3C CH3CH2CHCH3 (CH3)2CHCH2 3 CHCCCCH CH2CCCCH CHCCHCCH CH3CCCCH CH2CCHCCH CH2CCCCH2 CHCCH2CCH 3 c-CHCHCHCHC CHCCCCCH 3 CHCHO CH2COH CHCHOH CH3COH CH3OCH CH3OCH2 CHCCO CCCHO CH2CCO CHCCHO CH2CHCO CH2CCHO CHCHCHO CH2CHCHO CH3CHCO c-CH2C(O)CH2 c-CH2OCHCH CH3C(O)CH2 CH3CH2CO CH3CHCHO CH2CHCHOH

ANL0 556.4 656.6 315.7 349.5 802.2 495.7 547.4 655.5 295.7 327.7 395.7 434.9 318.5 327.3 329.3 342.1 371.6 176.2 156.0 158.4 226.3 241.0 247.4 252.2 264.2 75.5 90.9 97.8 727.1 568.5 577.3 417.9 441.7 451.9 459.5 525.4 682.4 261.3 121.5 137.5 69.5 129.9 13.7 290.0 468.3 131.9 135.2 102.8 178.1 190.9 −55.4 −53.4 31.8 65.3 −17.9 −17.3 −14.2 2.7

ΔANL0-F12 hydrocarbons: CxHy 0.2 0.3 0.1 0.2 −0.3 −0.1 0.4 −0.3 −0.5 −0.1 0.3 −0.5 −0.2 −0.3 0.1 0.0 −0.2 −0.2 −0.2 0.0 −0.1 0.0 0.0 0.0 −0.1 0.0 0.0 −0.7 −0.9 −0.5 −0.7 −0.4 −0.8 −0.5 0.0 −0.9 oxygenated hydrocarbons: CxOyHz 0.0 0.1 0.0 0.1 0.1 −0.1 −0.5 0.0 −0.3 −0.4 −0.3 −0.3 −0.2 −0.3 −0.3 −0.2 −0.1 −0.2 −0.1 −0.3 −0.1 6590

ΔANL1

ΔBAC

ΔSimmie

−0.7

1.4 −0.2 −0.1 −0.6

1.5 1.4 0.3 −1.5 −1.0 −0.6 −1.9 0.1

0.6 1.5 1.6 0.4

−0.3 1.1 0.6

2.1 −3.7 0.2 0.0 0.5 1.0 0.6 −0.5 0.0 2.3 −0.3 2.6 1.0 DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A Table 2. continued species CH2C(OH)CH2 CH2CH2CHO CH2OCHCH2 CH3C(OH)CH CH3CCHOH CH2CHCH2O c-CH2CH2CHO c−CH(CH2)OCH2 CH2CCH2OH c-CH2CHCH2O CH3C(OH)CH2 CH3CHCHOH CH2CHCH2OH c-CH2CH2CH(OH) c−CH(CH3)CH2O c-CH2CH2CH2O CH3C(OH)CH3 CH3CH2CHOH CH3CH(OH)CH2 CH3CHCH2OH CH2CH2CH2OH CH3CH(O)CH3 CH3CHOCH3 CH3CH2OCH2 CH3CH2CH2O CH2CH2OCH3 CH3CH2OCH3 3 CCCCO CHCCCO CCCCOH OCCHCCH CH2CCCO CHCCCOH c-CHCHCHCHO CH2OO OCH2OH 3 OCCO OCCHO OCCOH CCOOH CHCOO c−C(O)OCH2 OHCHCO OHCCOH CC(OH)2 CH2C(O)OH CH3C(O)O OCHCH2O c−CH(O)CH2O c−CH(CH2)OO c-CHOOCH2 CH2CHOO CHCHOOH CH2C(OH)2 OCHCH2OH OHCHCHOH CH3OCOH CH2CHOOH CH3CHOO

ANL0 3.8 33.1 101.0 101.4 103.4 112.4 116.1 120.3 124.5 142.7 −152.0 −133.8 −107.8 −82.1 −75.1 −60.0 −78.3 −55.8 −43.2 −39.5 −30.9 −25.4 −19.4 −12.2 −11.7 14.2 −194.2 606.5 381.0 641.3 207.5 225.3 323.2 −21.7 108.8 −158.5 14.1 −63.2 21.3 179.5 363.5 −164.5 −145.3 −19.3 141.6 −224.8 −182.2 −69.0 −33.1 161.1 203.5 118.5 230.5 −307.4 −303.6 −273.5 −171.1 −26.7 49.1

ΔANL0-F12 oxygenated hydrocarbons: CxOyHz −0.2 −0.2 −0.2 0.0 −0.1 0.0 −0.1 0.0 0.0 −0.1 −0.2 −0.2 −0.1 −0.1 −0.2 −0.2 −0.1 0.0 0.0 −0.1 0.0 0.1 −0.2 −0.1 0.0 −0.2 −0.3 −0.8 −0.6 −0.6 −0.5 −0.7 −0.6 −0.4 0.2 0.0 −0.5 −0.3 −0.2 0.0 0.3 −0.3 −0.3 −0.4 0.1 −0.4 −0.2 −0.1 −0.1 0.3 0.3 0.3 0.3 −0.2 −0.3 −0.2 −0.1 0.1 0.0

6591

ΔANL1

ΔBAC

ΔSimmie

3.0 5.0 3.2

−1.1 1.7 2.3 −0.4 −1.0 −0.4 −3.2 −0.7 0.3 0.0 1.4 0.8 −0.3 −0.3 −0.6 0.7 1.8 −0.3 −2.3 −1.8 −0.3

2.3

−3.2

0.9 −1.4 3.6 2.4 2.7

0.3 −0.9 1.3 0.8 4.5 2.0 6.1

5.7 −0.4 −2.0 −1.2 −0.7

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A Table 2. continued species OHCHCH2OH CH3OCHOH OHCH2OCH2 OHCH2CH2O CH3OCH2O CH3CH2OO CH2CH2OOH CH3CH2OOH CH3OOCH3 OCCCO OC(O)O c−C(O)OO OHC(O)O OCHOO OCOOH c−CH(OH)OO C(OH)3 OHCH2OO OCH2OOH CH3OOO CH(OH)3 OHCH2OOH CH3OOOH CH2CN CHCNH CH3CN CH3NC CH2CNH CHCNH2 c-CH2CHN 3 CH2CHN CH3CHN CH2CHNH CH3CNH CH2NCH2 CH3NCH CH2CNH2 CH3CHNH CH2CHNH2 CH3NCH2 CH3CHNH2 CH3CH2NH CH3NHCH2 CH2CH2NH2 CH3NCH3 CH3CH2NH2 C3 N CHCCN CH2CCN CHCHCN CHCCHN CHCNCH CH2CHCN CH2CHCHN c-CHCHCHCHNH CCCCCN CHCCCCN c-CHN2 NCCHN

ANL0 −196.1 −171.0 −166.7 −147.2 −130.9 −6.1 65.2 −141.1 −101.3 −92.9 −147.0 −161.5 −361.3 −100.3 −73.3 −216.7 −401.1 −154.2 −77.2 27.6 −588.0 −299.6 −71.0 264.0 396.3 80.1 183.8 193.5 253.2 280.8 369.2 210.7 219.7 233.7 252.2 278.6 294.8 55.9 69.0 95.1 133.4 167.4 169.5 175.2 176.7 −28.1 717.3 372.3 416.3 447.2 498.6 564.4 193.7 312.5 124.9 943.8 599.8 489.3 412.6

ΔANL0-F12 oxygenated hydrocarbons: CxOyHz 0.0 −0.2 −0.2 0.0 −0.1 0.2 0.2 0.1 −0.1 −1.1 0.0 0.0 −0.3 0.2 0.1 0.1 −0.1 0.3 0.2 0.4 −0.3 0.1 0.4 nitrogen compounds: CxNyOzHw −0.2 −0.2 −0.2 −0.4 −0.2 −0.3 0.1 0.1 0.0 −0.1 0.0 −0.2 −0.2 0.0 −0.1 −0.1 −0.2 0.0 0.0 −0.1 0.1 −0.1 0.0 0.0 −0.4 −0.2 −0.0 −0.1 −0.3 −0.2 0.0 −0.4 −0.5 −0.8 0.3 0.0 6592

ΔANL1

ΔBAC

ΔSimmie

0.5 2.0 1.4 −1.3 −0.2 1.1 −0.9 −1.7 −1.9

0.9 0.3

0.1 −2.2

−2.0

−0.3 −0.0 0.1 −0.1 0.1 0.1 0.4 −0.4

−0.5

0.2

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A Table 2. continued species NCNCH CHCNN NHNO NHNOH NH2NO NH2N(O)O NHN(O)OH NH2ONO CH2NO CHNOH NH2CO NH2CHO CH3NO CHCHNO CH2CHNO CH3CHNO CH3CH2NO NCOO CH2N(O)OH CH3OON CH3N(O)OH NH2C(O)NH2 a

ΔANL0-F12

ANL0 515.3 565.1 208.7 87.7 88.7 17.1 54.0 95.4 160.4 257.1 −8.7 −184.7 81.4 423.9 175.2 125.8 59.7 282.3 −3.0 81.1 −13.4 −214.3

ΔANL1

nitrogen compounds: CxNyOzHw −0.4 −0.1 0.2 0.2 0.1 0.2 0.3 0.4 0.0 0.1 −0.3 −0.5 0.2 0.3 0.1 −0.1 0.2 0.4 0.0 0.3 0.1 −0.6

ΔBAC

ΔSimmie

1.6

The 0 K heats of formation in kJ/mol. bT(Q) Correction.

Figure 3. Plot of the cumulative distribution function for the ANL0 corrections for higher order excitations [CCSDT(Q)/DZ], core− valence interactions, vibrational anharmonicity, relativistic, and Born− Oppenheimer effects.

Figure 4. Plot of the cumulative distribution function for various combinations of the corrections for higher order excitations [T(Q) denotes CCSDT(Q)/DZ-CCSD(T)/DZ], core−valence interactions (CV), vibrational anharmonicity (Anh.), relativistic (rel.), and Born− Oppenheimer (DBOC) effects.

mean correction and narrow distribution for the sum of the DBOC and relativistic corrections suggests that simply neglecting those two terms may not strongly degrade the accuracy of the predictions. The ANL1 predictions provide an opportunity to explore the convergence of the HOE correction with respect to degree of excitation and basis set size. As illustrated in Figure 5, both the CCSDT(Q) DZ to TZ basis set correction and the CCSDTQ(P) correction for higher excitations are generally greatly reduced in magnitude from the CCSDT(Q)/DZ correction of the ANL0 scheme. The mean values of the total, CCSDT(Q)/ DZ, CCSDT(Q) basis set, CCSDTQ(P), and the sum of the CCSDT(Q) basis set and CCSDTQ(P) corrections (which corresponds to the difference between the ANL1 and ANL0 HOE corrections) are −1.3, −1.7, 0.04, 0.29, and 0.34 kJ/mol,

while the corresponding RMSD from those mean values are 1.8, 2.1, 0.36, 0.51, and 0.64 kJ/mol. Overall, the distribution function for HOE corrections within the ANL1 scheme, which includes the CCSDT(Q)/DZ correction of the ANL0 scheme together with a CCSDTQ(P)/DZ correction and the DZ → TZ CCSDT(Q)/basis set correction, is rather similar to the distribution for the CCSDT(Q)/DZ HOE correction of the ANL0 scheme. The scatter plot provided in Figure 6 illustrates the correlation between the ANL0 HOE correction and the extra contribution to the HOE correction in the ANL1 scheme. Clearly, as the magnitude of the HOE increases, there is more uncertainty in the ANL0 prediction for it. In effect, this corresponds to saying that the more multireference character 6593

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A

Figure 5. Plot of the cumulative distribution function for various components of the correction for higher order excitation within the ANL1 scheme. The total term denotes the sum of (i) the CCSDT(Q)/DZ − CCSD(T)/DZ correction denoted as [CCSDT(Q)/DZ], (ii) the CCSDT(Q)/TZ-CCSDT(Q)/DZ correction denoted as [T(Q)/TZ-T(Q)/DZ], and (iii) the CCSDTQ(P)/DZCCSDT(Q)/DZ correction denoted as [TQ(P)/DZ-T(Q)/DZ].

Figure 7. Plot of the cumulative distribution function for the CBS limit of the core−valence correction as well as for the difference between the corrections obtained for DZ, TZ, and QZ basis sets.

values for the difference between the QZ and CBS limit suggests that this correction should be converged with respect to basis set to about 0.1 kJ/mol. Similarly, the low mean and narrow distribution for the difference between the CCSD(T) calculations for the a′5Z and CBS limits (cf. Figure 8) suggests that the CBS limit in ANL0 is

Figure 6. Scatter plot of the difference between the HOE correction for the ANL1 and ANL0 schemes versus that for the ANL0 scheme.

there is, the more difficult it is to quantitatively evaluate its effect with coupled cluster theory. Thus, one might rationally consider the uncertainties in our final ANL0 scheme to be a function of the HOE. Here we simply note that removing the species with an ANL0 HOE below −4 kJ/mol reduces the mean and RMSD from the mean for the ANL1 extra HOE contribution to 0.24 and 0.54 kJ/mol, respectively. Further restricting the set to include only the species with ANL0 HOE values ranging from −3.0 to 0.6 kJ/mol yields mean and RMSD from the mean values of 0.20 and 0.43 kJ/mol. For reference, we note that the species with differences in the HOE values (ANL1-ANL0) of magnitude greater than 1.0 kJ/mol are CH2OO (3.2), CHNO (2.2). C2 (2.1), NNO (1.8), OCHO (1.5), NHNO (1.4), CO2 (1.2), NHCO (1.2), O (−1.2), O3 (1.1), ONO (1.1), and OHNO (1.0). The basis set convergence of the distribution of core− valence corrections is illustrated in Figure 7. The mean values are 2.1, −0.05, −0.36, and −0.12 kJ/mol for the CBS limit and the DZ, TZ, and QZ differences from the CBS limit, respectively. The corresponding RMSD from the mean are 2.5, 1.2, 0.42, and 0.15 kJ/mol. The small mean and RMSD

Figure 8. Plot of the cumulative distribution function for the difference between the CBS limit and values for the a′TZ, a′QZ, and a′5Z basis sets.

reasonably well converged. The mean differences for the a′DZ, a′TZ, a′QZ, and a′5Z basis sets relative to the estimated CBS limit are −10.7, −0.94, 1.5, and 0.66 kJ/mol, respectively. The RMSD from these mean values are 12.3, 4.3, 1.8, and 0.77 kJ/ mol. The basis set convergence for the DZ-F12, TZ-F12, and QZ-F12 series is similar (cf. Figure 9) with the QZ-F12 convergence roughly correlating with the a′5Z convergence. The mean values for the differences from the F12-CBS limit are −4.31, −1.63, and −0.56 kJ/mol and the RMSD from the mean are 2.7, 0.96, and 0.33 kJ/mol. However, unfortunately, the CBS-F12 and CBS-a′nZ limits are not identical with mean differences and RMSD from the mean of −0.06, and 0.31 kJ/ mol, respectively. The calculations with the a′6Z basis set in the ANL1 scheme provide an opportunity to further examine the basis set convergence. The mean deviation between the {a′5Z,a′6Z} and 6594

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A

the CCSD(T) analysis, and the anharmonicity analysis. We estimate a slightly larger 2σ uncertainty of 1.5 kJ/mol in the ANL0-F12 method due to the increased deviation from the {a′5Z,a′6Z} CBS limit. For the ANL1 method we estimate a slightly smaller 2σ uncertainty of 0.9 kJ/mol due to the more well-converged HOE correction and CBS limits. 4.2. Comparison with W4 Theory. The Weizmann-4 (W4) method is similar in spirit to the present methods, but with a somewhat more extensive consideration of the corrections for higher order excitations and the use of the best available vibrational ZPE, which mostly, but not always means a ZPE computed at a higher level than in the present work. Some time ago, W4 theory29 was shown to reproduce an ATcT set of atomization energies for 35 molecules with a 3σ uncertainty of under 1 kJ/mol.85 With this high level of accuracy, these predictions provide a useful database for testing the present ANL schemes. The full W4 analysis of Karton et al.29,85 included the prediction of the atomization energies for a set of 140 small first-and second-row species, with 68 of those species overlapping with the species of interest here. These atomization energies were converted to the W4 heats of formation that are reported in Table 1 via the ATcT heats of formation for the reference atoms. The corresponding cumulative distribution functions for the differences between the W4 and ANL predictions for ΔfHo(0 K) are illustrated in Figure 11. The mean differences relative to

Figure 9. Plot of the cumulative distribution function for the difference between the F12-CBS limit and values for the DZ-F12, TZ-F12, QZF12 basis sets and the CBS limit from {a′QZ,a′5Z} calculations.

{a′QZ,a′5Z} CBS limits is 0.031 kJ/mol, while that from the (TZF,QZF) F12-CBS limit is 0.066 kJ/mol. Apparently, the {a′QZ,a′5Z} CBS estimate is slightly better converged than the (TZF,QZF) F12-CBS estimate, although the difference may also be indicative of a systematic shortcoming of the a′nZ series. The difference between the B3LYP- and MP2-based VPT2 vibrational anharmonicity corrections provides a qualitative indication of the uncertainties in the predicted anharmonicity correction. The cumulative distribution function for this difference (for the cases where both methods appear to yield reasonable results) is plotted in Figure 10. The mean and

Figure 11. Plot of the cumulative distribution function for the differences between the W4, ANL1, and ANL0 predictions for the ΔfHo(0 K).

ANL1 and ANL0 are 0.24 and 0.28 kJ/mol, respectively, while the RMSD are 0.80 and 0.88 kJ/mol. These RMSDs are biased by the presence of particularly large deviations for a handful of species, e.g., 94% of the ANL1 values lie within 1.0 kJ/mol of the W4 values. For the other 4 species (OOOH, CH3NH, ONOH, and C2) the W4 deviations from ANL1 are 4.2, −2.6, −1.4, and 1.2 kJ/mol, while the corresponding deviations from ANL0 are 4.0, −2.6, −0.5, and 2.9 kJ/mol. For OOOH, the majority of the difference comes from the W4 focus on the cis isomer, while the trans isomer is predicted here to be 2.8 kJ/mol lower in energy after ZPE correction. Nevertheless, it is worth noting that the ANL1 HOE correction is extraordinarily large (−9.7 kJ/mol) indicating extensive multireference character and significant larger uncertainty than normal for the ANL values. For CH3NH, the W4 calculations employ a scaling of the B3LYP/cc-pVTZ calculated ZPE for the

Figure 10. Plot of the distribution of differences between VPT2 anharmonic corrections to the zero-point energy obtained with the MP2/cc-pVTZ and B3LYP/cc-pVTZ methods.

RMSD from the mean are 0.06 and 0.42 kJ/mol, but a few outliers again strongly affect this RMSD. 68% of the values fall within 0.22 kJ/mol of the mean, while 95% fall within 0.61 kJ/ mol. Crude estimates of the overall uncertainty in the ANL schemes can be obtained from these various observations. In particular, we estimate the 2σ uncertainty in the ANL0 scheme as 1.4 kJ/mol from crudely estimated 1σ uncertainties of 0.5, 0.08, 0.3, and 0.3 for the uncertainties in the HOE correction, the CBS limit of the core−valence correction, the CBS limit of 6595

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A

CH2OO, NHNO, and C2) from the test set increases the magnitude of the mean deviation slightly from −0.01 to −0.08 kJ/mol, while the RMSD decreases from 0.56 to 0.46 kJ/mol. 4.3. Comparison with ATcT Values. The cumulative distribution function of the apparent errors for the ANL0 predictions (relative to the ATcT values) is plotted in Figure 12. The mean deviations for the ANL0, ANL0-F12, and ANL1

anharmonic ZPE. The 1.3 kJ/mol difference between this W4 ZPE and the present ANL1 ZPE, suggests that the W4 value may not be particularly reliable in this case. C2 is well-known to involve a high degree of multireference character. A recent high level theoretical analysis84 yields a ΔfHo(0 K) for C2 of 819.5 ± 0.8 kJ/mol, which is in good agreement with the present ANL1 value (819.7 kJ/mol) and reasonably satisfactory agreement with the ANL0 value (818.0 kJ/mol). In contrast, the W4 value of 820.9 kJ/mol lies 1.4 kJ/mol higher, while the higher level W4.4 value is 820.0 kJ/mol. The high level calculations of ref 84 generally employ more extended basis sets and also consider separability and other issues more extensively, and so can be considered to be more reliable. Nevertheless, the current ATcT value of 820.2 lies closer to the W4.4 value, due to the statistical reinforcement between additional thermochemical cycles available in the TN and a number of such calculations and their presumed lower error bars. For HONO the difference between the W4 and ANL1 predictions again appear to be due to a difference of −1.8 kJ/mol in the predicted ZPE. However, it is difficult to rationalize this difference, as the primary difference in the two ZPE treatments for this case appears to involve the inclusion of core−valence effects in the W4 treatment. However, our own all-electron CCSD(T)/cc-pcVTZ ZPE calculations yields a negligible difference of 0.02 kJ/mol. Removing these 4 outliers reduces the RMSD differences in the W4-ANL1 numbers to 0.48 kJ/mol and W4-ANL0 to 0.59 kJ/ mol, while the mean differences are only slightly reduced to 0.23 and 0.24 kJ/mol, respectively. For the W4 to ANL0 comparison there are 3 additional apparent outliers (with differences greater than 1.0 kJ/mol); cCHOCH, O3, and CHNO with deviations of 1.7, 1.8, and 1.2 kJ/mol, respectively. O3 is again a case with substantial multireference character as indicated by the ANL0 HOE correction of −13.6 kJ/mol. For CHNO, the primary difference is related to a difference in the ZPE values of 1.5 or 0.7 kJ/mol for W4 relative to ANL1 and ANL0, respectively. Notably, in this case, there is considerable ambiguity in the anharmonic ZPE correction, with a B3LYP/cc-pVTZ VPT2 analysis yielding a positive anharmonicity correction of 0.61 kJ/mol, while an MP2/cc-pVTZ VPT2 analysis yields a negative anharmonicity correction of −0.95 kJ/mol. We employ the MP2/cc-pVTZ analysis. For c-CHOCH, the W4 ZPE is 0.2 and 0.8 kJ/mol higher than the ANL1 and ANL0 values, respectively. Overall, it appears that the different ZPE treatments contribute most strongly to the differences between the W4 and ANL treatments. Adding O3 and CHNO to the list of removed species further reduces the W4-ANL1 and W4ANL0 RMSD only slightly to 0.46 and 0.53 kJ/mol, while the mean deviations are also slightly further reduced to 0.21 and 0.20 kJ/mol. For the ANL1 to ANL0 comparison there are 8 species with absolute differences greater than 0.9 kJ/mol: CH2OO (2.3), 3 CHCCH (−1.9), C2 (1.7), CH3O (−1.6), NHNO (1.6), NNO (1.2), O (−1.1), and O3 (0.9), where the numbers in parentheses denote the ANL1-ANL0 energies in kJ/mol. A number of these species are again clearly multireference in nature (O3, CH2OO, NHNO, and C2, with CCSDT(Q)/DZ corrections of −13.6, −7.4, −6.8, −4.9, and −3.6 kJ/mol, respectively). The differences for NHNO and O are also primarily due to differences between the ANL1 and ANL0 HOE corrections, while for 3CHCCH and CH3O changes in the ZPE are the largest contributors to the difference. Removing the four clearly multireference molecules (O3,

Figure 12. Plot of the cumulative distribution function (solid) and Gaussian fit (dashed) to the peak for the errors in the ANL0 predicted ΔfHo(0 K) relative to ATcT values.

methods are −0.22, −0.15, and −0.44 kJ/mol, respectively. The corresponding RMSD are 1.08, 1.12, and 1.16 kJ/mol. The distribution is seen to be nearly Gaussian in shape for small deviations (up to 1 or 2 kJ/mol), but contains a significant tail with larger deviations. Some understanding of these apparent outliers is central to properly evaluating the accuracy of the present ANL methods. A list of 15 “problem cases”, i.e., the set with errors for any one of the three ANL methods greater than 1.4 kJ/mol in magnitude, is provided in Table 3 together with their ATcT uncertainty, T1 and %TAE[(T)] diagnostics,29,93 HOE corrections within ANL0 and ANL1, and a label noting whether the ATcT determination is primarily based on theoretical calculations. The two diagnostics are commonly taken as a measure of multireference effects, while the corrections provide an alternative direct indication of such effects. These problem cases are seen to correlate with either large HOE corrections (the average HOE correction for the full set of 348 ANL0 species is −1.8 kJ/mol) or large uncertainty in the ATcT values (as discussed above, the average ATcT uncertainty for the full set of 150 species is 0.64 kJ/mol). For OOOH and ON(O)O, the problems appear to be related to large multireference character as noted by the extraordinarily large magnitude for the higher order corrections to the enthalpy. The largest magnitude ANL0 HOE corrections for the full set of studied species were for OC(O)O (−24.1), CH3OON (−15.7), O3 (−13.6), ON(O)O (−12.9), CH3OOO (−12.4), OOOH (−9.9), CH2OO (−7.5), NHNO (−7.1), ONOOH (−6.4), 3CCCCO (−6.3), CH3C(O)O (−6.3), CH3CHOO (−6.3), CH3CHNO (−6.0), CHNN (−6.0), OCHO (−5.7), CHCHNO (−5.6), CCCHO (−5.5), ONO (−5.3), OHC(O)O (−5.3), CHCCCO (−5.2), CCCCOH (−5.1), CN (−5.0), and C2 (−5.0), (with ANL0 HOE corrections in parentheses in kJ/mol). As one might expect, the species with large HOE corrections correlate with species 6596

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A Table 3. ΔfHo(0 K) (kJ/mol) and Diagnostics for Apparent Problem Speciesa

Key: (a) ATcT value for ΔfHo(0 K) in kJ/mol. (b) Uncertainty in the ATcT values in kJ/mol. Species with uncertainties >1.0 kJ/mol are highlighted in bold blue. (c) Difference between ANL0 and ATcT values for ΔfHo(0 K) in kJ/mol. (d) Difference between ANL0-F12 and ATcT values for ΔfHo(0 K) in kJ/mol. (e) Difference between ANL1 and ATcT values for ΔfHo(0 K) in kJ/mol. (f) T1 Diagnostic. Species with T1 diagnostics >0.035 are highlighted in bold green. (g) %TAE[(T)] diagnostic. Species with %TAE[(T)] > 4.5 are highlighted in bold purple. (h) CCSDT(Q)/cc-pVDZ − CCSD(T)/cc-pVDZ correction to the enthalpy in kJ/mol. Species with T(Q) correction < −5.0 kJ/mol are highlighted in bold red. (i) CCSDTQ(P)/cc-pVDZ − CCSDT(Q)/cc-pVDZ + CCSDT(Q)/cc-pVTZ − CCSD(T)/cc-pVTZ correction to the enthalpy in kJ/ mol. (j) A yes indicates that greater than 50% of the provenance for this species in ver. 1.122b of the ATcT thermochemical network comes from prior theoretical calculations. a

schemes such as CBS-QB3, G4, etc.) for a number of separate data entries, though such correlations may be allayed by the independent empirical correctors incorporated in these methods. As discussed elsewhere28 and reiterated in section 3, the presence of uncertainties in the benchmark set that are comparable to the expected accuracy of the theoretical predictions complicates the error analysis, because the accuracy obtained by benchmarking is a composite of the inherent accuracy of the benchmarked method and the accuracy of the benchmarks. The plots in Figure 14 illustrate the dependence of the calculated RMSD on the particular set of species included in the comparison between ANL0 and ATcT.

that have large T1 and/or %TAE[(T)] diagnostics (cf., the scatter plot in Figure 13).

Figure 13. Scatter plot of the correlation between the ANL0 HOE correction and the T1 and %TAE[(T)] diagnostics.

Only the 3CHCCH, CH2CHOH, and NCCN species do not have either an ATcT uncertainty that exceeds 1.0 kJ/mol or an ANL0 HOE correction that exceeds 5.0 kJ/mol. Interestingly, these three problem cases are all estimated to have a theory provenance of greater than 50%. Perhaps not surprisingly, the set of species with relatively large ATcT uncertainties correlates with the set of species whose ATcT provenance is primarily theoretical. In principle, the presence of correlations in the contributing theoretical data could result in artificially low predicted error bars. Method correlations might arise, for example, from the neglect of HOE corrections (as in many

Figure 14. Plot of the RMSD for the ANL0 − ATcT difference as a function of the maximum ATcT uncertainty in the comparison set and the maximum absolute value for the HOE ANL0 correction. 6597

DOI: 10.1021/acs.jpca.7b05945 J. Phys. Chem. A 2017, 121, 6580−6602

Article

The Journal of Physical Chemistry A Table 4. Error Properties for ΔfHo(0 K) (kJ/mol) from the ANL Methodsa data setb

no. of species

ATcT unc. mean 2σ

method

mean error

RMSD error

ATcT ∩ ANL0

150

0.68

ATcT ∩ ANL0 ∩ HOE < 5

141

0.68

ATcT unc. < 1.0 ∩ ANL0 ∩ HOE < 5

126

0.50

ATcT unc. < 0.6 ∩ ANL0 ∩ HOE < 5

91

0.33

ATcT ∩ ANL1

93

0.61

ATcT ∩ ANL1 ∩ HOE < 5

86

0.59

ATcT unc. < 1.0 ∩ ANL1 ∩ HOE < 5

79

0.46

ATcT unc. < 0.6 ∩ ANL1 ∩ HOE < 5

61

0.29

ANL0 ANL0-F12 ANL0 ANL0-F12 ANL0 ANL0-F12 ANL0 ANL0-F12 ANL1 ANL0 ANL0-F12 ANL1 ANL0 ANL0-F12 ANL1 ANL0 ANL0-F12 ANL1 ANL0 ANL0-F12

−0.22 −0.15 −0.10 −0.07 −0.12 −0.09 −0.12 −0.09 −0.44 −0.39 −0.27 −0.34 −0.26 −0.15 −0.33 −0.25 −0.15 −0.28 −0.22 −0.12

1.08 1.12 0.77 0.82 0.63 0.70 0.50 0.60 1.16 1.19 1.25 0.80 0.72 0.80 0.70 0.62 0.69 0.50 0.48 0.60

a

ANL0 and ANL1 methods uncertainties and errors relative to ATcT values in kJ/mol. bData sets correspond to intersection of ATcT set (either full set or set with uncertainties less than 1.0 or 0.6 kJ/mol) with either ANL0 or ANL1 sets (either full set or with the magnitude of the HOE less than 5.0 kJ/mol).

Notably, the RMSD decreases dramatically when the ANL0 HOE corrections >5 kJ/mol are excluded from the set, but then stays more or less constant for further reductions. The RMSD also gradually decreases as the set is restricted to ATcT values with smaller and smaller maximum values of the uncertainties, but with a plateau near 0.6 kJ/mol when the limiting maximum ATcT uncertainty is in the range from 0.7 to 1.6 kJ/mol. The mean and RMSD for the ANL1, ANL0, and ANL0-F12 predicted ΔfHo(0 K) are tabulated in Table 4 for various data sets. The cumulative distribution functions for the errors in the ANL0, ANL0-F12, and ANL1 calculated ΔfHo(0 K) are provided in Figure 15 for a data set with 84 species (consisting of the intersection between the ATcT data set and the ANL1 data set, but excluding the problem species listed in Table 3). Restricting the data sets to those for which the present methods

are expected to be accurate (i.e., with HOE < 5 kJ/mol) and for which the ATcT provides a proper validation set (i.e., with ATcT uncertainties