summary: | The fundamental importance of alkanes as organic chemistry building blocks and in industrial chemistry (particularly petrochemistry) is self-evident to any chemist. Linear and branched lower alkanes are the principal components of gaseous and liquid fossil fuels. Accurate knowledge of their thermodynamic properties is essential for reliable computational modeling of combustion processes. (We note that one of us, J.M.L.M., is a member of a IUPAC task group working in this area.1) Aside from their practical relevance, alkanes present some intriguing methodological issues. The importance of accurate zero-point vibrational energies and diagonal Bo -Oppenheimer corrections has been discussed previously,2 and this applies to both wave function ab initio and density functional methods. While post-CCSD(T) computational thermochemistry methods like W4 theory3,4 or HEAT5-7 have no trouble dealing with systems that, from an electronic structure point of view, are much more taxing than alkanes, their steep cost scaling makes application to higher alkanes (or higher hydrocarbons in general) impractical at present. Density functional theory seems to be the obvious alte ative. However, in recent years a number of authors8-15 have pointed to a disturbing phenomenon;16 the error in computed atomization energies of n-alkanes grows in direct proportion to the chain length. In addition, these same authors found that popular DFT methods have significant problems with hydrocarbon isomerization energies in general and alkane isomerization energies in particular. This latter problem appears to be related to the poor description of dispersion by most DFT functionals and can be remedied to a large extent by empirical dispersion corrections.11 For molecules as chemically systematic as alkanes, a computationally more cost-effective approach than brute force atomization energy calculations is the use of bond separation reactions, such as isodesmic17 and homodesmotic18 reactions. Recently, Schleyer and co-workers19,20 discussed the concepts of “protobranching” and of “hypohomodesmotic reactions”, that is, reactions which, in addition to being isodesmic (that is, conserving numbers of each formal bond type), conserve the number of C atoms in each hybridization state and the hapticity (primary, secondary, tertiary, quate ary). The latter is a refinement of the earlier “homodesmotic reaction” concept.18 They established a consistent hierarchy of hydrocarbon reaction types that successively conserve larger molecular fragments, atomization g isogyric g isodesmic g hypohomodesmotic g homodesmotic g hyperhomodesmotic, which provides a converging sequence in the sense that the energetic components of the reaction cancel to a larger extent between reactants and products as the reaction hierarchy is traversed. In the present work, we obtain “quasi-W4” atomization energies for C4-C8 alkanes through the use of isodesmic and hypohomodesmotic reaction cycles that involve only methane, ethane, and propane in addition to one larger alkane. The reaction energies are calculated at the W3.2lite or W1h levels, while for methane, ethane, and propane, W4 benchmark values are used. We shall show that the reaction energies of hypo- * To whom correspondence should be addressed. E-mail: gershom@ weizmann.ac.il. J. Phys. Chem. A XXXX, xxx, 000 A 10.1021/jp904369h CCC: $40.75 XXXX American Chemical Society Downloaded by AUSTRIA CONSORTIA on July 6, 2009 Published on July 1, 2009 on http://pubs.acs.org | doi: 10.1021/jp904369h homodesmotic reactions and judiciously selected isodesmic reactions are well converged even at the W1h level. The use of hypohomodesmotic reactions leads to near-perfect cancellation of valence correlation effects, and the use of judiciously selected isodesmic reactions leads to near-perfect cancellation of post- CCSD(T) correlation effects. We will then proceed to evaluate a number of DFT functionals and composite ab initio thermochemistry methods against the reference values obtained, both for the atomization and for the isomerization energies. II. Computational Methods All calculations were carried out on the Linux cluster of the Martin group at Weizmann. DFT geometry optimizations were carried out using Gaussian 03, revision E.01.21 The B3LYP22-24 DFT hybrid exchange-correlation (XC) functional was used in conjunction with the pc-225 polarization consistent basis set of Jensen. All large-scale self-consistent field (SCF), CCSD, and CCSD(T) calculations26,27 were carried out with the correlation- consistent family of Dunning and co-workers28-32 using version 2006.1 of the MOLPRO33 program system. All singlepoint post-CCSD(T) calculations were carried out using an OpenMP-parallel version of M. Ka´llay’s general coupled cluster code MRCC34 interfaced to the Austin-Mainz-Budapest version of the ACES II program system.35 The diagonal Bo - Oppenheimer correction (DBOC) calculations were carried out using its successor CFOUR.36 The computational protocols of Wn theories W1,37,38 W3.2lite,39 and W43 used in the present study have been specified and rationalized in great detail elsewhere.3,37-39 (Throughout, W3.2lite refers to variant W3.2lite(c) as described in ref 39.) The use of the Wnh variants of the Wn methods, in which the diffuse functions are omitted from carbon and less electronegative elements, is of no thermochemical consequence for neutral alkanes,38 but computer resource requirements are substantially reduced. For the sake of making the paper self-contained, we will briefly outline the various steps in W3.2lite theory and in W4h theory for first-row elements. • The reference geometry and ZPVE correction are obtained at the B3LYP/pc-2 level of theory for W3.2lite and at the CCSD(T)/cc-pVQZ level for W4h. • The ROHF-SCF contribution is extrapolated using the Karton-Martin modification40 of Jensen’s extrapolation formula41 For W3.2lite and W4h, the extrapolations are done from the cc-pV{Q,5}Z and cc-pV{5,6}Z basis set pairs, respectively. • The RCCSD valence correlation energy is extrapolated from these same basis sets. Following the suggestion of Klopper,42 Ecorr,RCCSD is partitioned in singlet-coupled pair energies, tripletcoupled pair energies, andTˆ 1 terms. TheTˆ 1 term (which exhibits very weak basis set dependence) is simply set equal to that in the largest basis set, while the singlet-coupled and triplet-coupled pair energies are extrapolated using A + B/LR with RS ) 3 and RT ) 5. • The (T) valence correlation energy is extrapolated from the cc-pV{T,Q}Z basis set pair for W3.2lite and from cc-pV{Q,5}Z for W4h. For open-shell systems, the We er-Knowles-Hampel (aka, MOLPRO) definition43 of the restricted open-shell CCSD( T) energy is employed throughout, rather than the original Watts-Gauss-Bartlett27 (aka ACES II) definition. • The CCSDT - CCSD(T) difference,Tˆ 3 - (T), in W3.2lite is obtained from the empirical expression 2.6 × cc-pVTZ(no f 1d)(no p on H) - 1.6 × cc-pVDZ(no p on H), where the CCSDT energy is calculated using ACES II. In W4h, it is instead extrapolated using A + B/L3 from cc-pV{D,T}Z basis sets. • The difference between ACES II and MOLPRO definitions of the valence RCCSD(T) definition is extrapolated from ccpVDZ and cc-pVTZ basis sets. One-half of this contribution is added to the final result, as discussed in the appendix of ref 3. • Post-CCSDT contributions in W3.2lite are estimated from UCCSDT(Q)/cc-pVDZ(no p on H) - UCCSDT/cc-pVDZ(no p on H) scaled by 1.1. In W4h, the connected quadruples are obtained as 1.1 × [UCCSDT(Q)/cc-pVTZ - UCCSDT/ccpVTZ + UCCSDTQ/cc-pVDZ - UCCSDT(Q)/cc-pVDZ], while the contribution of connected quintuple excitations is evaluated at the CCSDTQ5/cc-pVDZ(no d) level. • The inner-shell correlation contribution, in both cases, is extrapolated from RCCSD(T)/cc-pwCVTZ and RCCSD(T)/ccpwCVQZ calculations. • The scalar relativistic contribution, again in both cases, is obtained from the difference between nonrelativistic RCCSD(T)/ cc-pVQZ and second-order Douglas-Kroll RCCSD(T)/DK-ccpVQZ calculations. • Atomic spin-orbit coupling terms are taken from the experimental fine structure. • Finally, a diagonal Bo -Oppenheimer correction (DBOC) is obtained at the ROHF/cc-pVTZ level. The main changes in W1h relative to W3.2lite are that (a) the SCF component is extrapolated from the cc-pV{T,Q}Z basis sets, using the formula A + B/L5; (b) the valence RCCSD component is extrapolated from the same basis sets, using A + B/L3.22; (c) the valence parenthetical triples, (T), component is extrapolated from cc-pV{D,T}Z basis sets, using A + B/L3.22; (d) inner-shell correlation contributions are evaluated at the CCSD(T)/MTsmall level; and (e) post-CCSD(T) correlation effects as well as the DBOC are completely neglected. The CCSDTQ5/cc-pVDZ(no d) calculation for propane proved to be too taxing even for our strongest machine (8 core, Intel Cloverton 2.66 GHz, with 32 GB of RAM). For the alkanes for which we do have this term, CH4 and C2H6, it is practically zero (0.00 and 0.01 kcal/mol, respectively); therefore, for propane it was safely neglected. The anharmonic zero-point vibrational energy (ZPVE) of propane, propene, propyne, and allene was calculated using the following equation44 where the cubic, quartic, and kinetic energy terms were computed at the MP2/cc-pVTZ level of theory and the harmonic term was partitioned into valence and core-valence contributions, which were calculated at the CCSD(T)/cc-pVQZ and CCSD(T)/MTsmall levels of theory, respectively. (For propane, we resorted to a CCSD(T)/cc-pVTZ calculation since the CCSD(T)/cc-pVQZ proved too daunting; on the basis of the results for the other systems, this is expected to have little effect; for example, the differences between the harmonic ZPVE calculated with the two basis sets are 0.02, 0.02, 0.02, and 0.01 kcal/mol for methane, ethane, propene, and allene, respectively.) EHF,L ) EHF,∞ + A(L + 1) exp(-9√L) (1) ZPVE ) 1 2Σi ωi - 1 32Σ ijk φiikφkjj ωk - 1 48Σ ijk φijk 2 ωi + ωj + ωk + 1 32Σ ij φiijj + Zkinetic (2) B J. Phys. Chem. A, Vol. xxx, No. xx, XXXX Karton et al. Downloaded by AUSTRIA CONSORTIA on July 6, 2009 Published on July 1, 2009 on http://pubs.acs.org | doi: 10.1021/jp904369h Unless noted otherwise, experimental data for the heats of formation at 0 K were taken from the NIST Computational Chemistry Comparison and Benchmark Database (CCCBDB).45 The atomization energies quoted in CCCBDB assume CODATA46 values for the atomic heats of formation at 0 K; however, particularly for the carbon atom, the ATcT value47 (170.055 ( 0.026 kcal/mol) is significantly higher than the CODATA value (169.98 ( 0.11 kcal/mol). Consequently, using the ATcT value in converting ΔH°f,0K to atomization energy raises the atomization energy over the CCCBDB value by m × 0.075 kcal/mol for a system with m carbon atoms. Thus, throughout the paper, the experimental TAE0 were obtained from the heats of formation at 0 K using ATcT atomic heats of formation at 0 K (C, 170.055 ( 0.026; H, 51.633 ( 0.000 kcal/ mol).47 In cases where only experimental heats of formation at 298 K are available from CCCBDB (n-heptane, n-octane, isoheptane, and isooctane), they were first converted to 0 K using the H298 - H0 for H2(g) of 2.024 ( 0.000 and that for C(cr,graphite) of 0.251 ( 0.005 kcal/mol from CODATA46 and the molecular heat content functions from the TRC (Thermodynamic Research Center) tables,48 which are the source of virtually all CCCBDB enthalpy functions for the species considered in the present paper. In order to facilitate direct comparison with experiment, we have also converted our calculated atomization energies to heats of formation at room temperature, ΔH°f,298K. Rather than mix our calculated atomization energies with the TRC enthalpy functions, we have calculated our own H298 - H0 for the alkanes. The translational, rotational, and vibrational contributions were obtained within the RRHO (rigid rotor-harmonic oscillator) approximation from the B3LYP/pc-2 calculated geometry and harmonic frequencies. Inte al rotation corrections were obtained using the Ayala-Schlegel method,49 again on the B3LYP/pc-2 potential surface. This left us with the issue of correcting for the ensemble of low-lying conformers of the alkanes; by way of illustration, n-butane through n-octane has 2, 4, 12, 30, and 96 unique conformers, respectively. The relative energies of these conformers (which are surprisingly sensitive to the level of theory as they are strongly driven by dispersion) were the subject of a recent benchmark study by our group.50 While large basis set CCSD(T) calculations for all conformers of the heptanes and octanes proved too costly (primarily for those without any symmetry), it was found in ref 50 that the B2K-PLYP-D double hybrid functional51 with an empirical dispersion correction, in conjunction with a sufficiently large basis set, tracks the CCSD(T) reference data50 for n-butane, n-pentane, and n-hexane exceedingly closely, and this is the approach we have followed for all systems with more than one conformer in this work. Dispersion corrections for the DFT energies (denoted by the suffix “-D”) were applied using our implementation of Grimme’s expression52,53 where the damping function is taken as and C6 ij ≈ (C6 i C6 j )1/2; Rr ) RvdW,i + RvdW,j is the sum of the van der Waals radii of the two atoms in question, and the specific numerical values for the atomic Lennard-Jones constants C6 i and the van der Waals radii have been taken from ref 52, whereas the length scaling sR ) 1.0 and hysteresis exponent R ) 20.0 are as per ref 53. Equation 3 has a single functional-dependent parameter, namely, the prefactor s6. This was taken from refs 52 and 53 for BLYP, B3LYP, and PBE and from ref 54 for the double hybrids and was optimized in the present work for the remaining functionals. These were, for the most part, optimized against the S22 benchmark set of weakly interacting systems |