Assessment of Density Functional Theory for Thermochemical

Article pubs.acs.org/JPCA

Assessment of Density Functional Theory for Thermochemical Approaches Based on Bond Separation Reactions Dirk Bakowies* Laboratory of Physical Chemistry, ETH Zürich, CH 8093 Zürich, Switzerland S Supporting Information *

ABSTRACT: The recently proposed ATOMIC protocol is a fully ab initio thermochemical protocol that rests upon the concept of bond separation reactions (BSRs) to correct for systematic errors of composite wave function approaches. It achieves high accuracy for atomization energies and derived heats of formation if basis set requirements for all contributing components are balanced carefully. The present work explores the potential of density functionals as possible replacements of composite wave function approaches in the ATOMIC protocol. Twenty density functionals are examined for their accuracy in thermochemical predictions based on calculated bond-separation energies and precomputed highlevel data for the small parent molecules entering BSRs. The best density functionals outperform CCSD (coupled cluster with singles and doubles excitations), but none reaches the accuracy of well-balanced composite wave function approaches that consider quasiperturbational connected triples excitations at least with small basis sets. Some functionals show unexpected problems with bond separation reactions and are analyzed further with a model of empirically calibrated bond additivity corrections. Finally, the benefit of adding empirical dispersion terms to common density functionals is analyzed in the context of BSR-corrected thermochemistry. particular, the Gaussian (Gn)31−34 and CBS families of methods,35−38 and the more recent correlation-consistent composite approach (ccCA).39−41 Some, but not all, of these approaches use empirically calibrated correction terms to remove residual systematic errors observed in applications to atomization processes. The recently introduced ATOMIC protocol42,43 is comparable in computational expense to the popular G3 and G4 methods but follows a quite different strategy to achieve high accuracy without resorting to excessively expensive calculations. It is based on the concept of bond-separation reactions (BSRs) and treats the process of bond formation with high-level methods embedded in a lower-level treatment of the perturbation experienced through the chemical environment. BSRs were introduced by Pople and co-workers more than 40 years ago,44,45 they essentially describe the dissection of a molecule into its bonds, represented by molecules containing only two non-hydrogen atoms each. Stoichiometric balance is restored through the addition of simple hydrides on the reactant side. The BSR of oxirane into methanol and ethane may serve as an example:

1. INTRODUCTION Ab initio quantum chemistry nowadays offers a rich tool set for accurate thermochemistry. Several methods,1,2 including the Weizmann3−6 and HEAT protocols,7−9 focal point analysis,10−16 and the procedures of Klopper and various collaborators17−22 and of Dixon and Feller,2,23−26 allow one to calculate very precise atomization energies and derive heats (enthalpies) of formation to an accuracy that often challenges direct experimental determination. The earliest among these methods, focal point analysis, is designed to provide uncertainties of computed results as well, as it systemtically probes for incremental convergence of both the one- and Nparticle descriptions toward the “focal point” of a theoretically exact treatment. The currently most reliable thermochemical data and associated uncertainties are obtained from thermochemical networks established from either experimental (ATcT)27 or theoretical (NEAT)28 reaction enthalpies and related data. All of the mentioned theoretical methods approach in some way the complete basis set limit of coupled cluster theory with single and double excitations and quasiperturbational treatment of triples excitations (CCSD(T)),29,30 typically using either basis set extrapolation or explicitly correlated wave functions. Beyond this, details differ between various methods, but relativistic corrections are often added and, depending on the level of accuracy needed, also higher-order excitation corrections and diagonal Born−Oppenheimer terms. The significant computational expense of most methods unfortunately limits their use to rather small molecules. Less elaborate approaches have thus found more widespread use, including, in © 2012 American Chemical Society

C2H4O + 2CH4 + H 2O → C2H6 + 2CH3OH

Simple hydrides (here: methane, water) and product molecules (here: ethane, methanol) will be collectively referred to as BSR parent molecules or BSR prototypes. The ATOMIC Received: October 30, 2012 Revised: December 6, 2012 Published: December 6, 2012 228

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

Table 1. Atomization Energies (AE): RMS Errors (kcal/mol)a,b with BSR corrections DZ

TZ

QZ

5Z

13.6

10.8

9.9

9.4

BP86 BLYP PBE

8.6 7.3 9.3

6.7 5.7 7.3

6.0 5.1 6.5

VSXC TPSS M06-L

10.3 7.0 9.0

9.5 5.4 8.1

8.7 4.9 7.7

B3LYP B3PW91 PBE0 TPSSh M05 M06

5.4 6.2 6.5 5.9 8.6 6.7

3.5 4.4 4.4 4.2 7.5 5.7

2.8 3.7 3.6 3.8 6.9 5.0

BMK M05-2X M06-2X M06-HF

5.2 4.4 4.3 7.8

2.4 2.3 2.8 8.2

2.0 1.9 2.5 7.6

B2-PLYP B2-PLYP-D B2GP-PLYP

5.9 5.9 5.6

3.9 3.8 3.6

3.1 3.0 2.7

HF MP2c CCSDc,d CCSD(T)c,d model B3e

8.8 7.4 3.8 4.1

8.5 5.3 2.4 1.6

8.7 4.6 2.5 0.6

SVWN3

regular A

B

DZ

TZ

QZ

LDA Functionals 9.4 112.7 119.8 120.7 GGA Functionals 5.7 5.8 5.6 19.5 25.3 26.1 4.8 5.0 4.8 16.2 15.6 16.1 6.1 6.2 5.9 26.5 31.4 32.1 Meta-GGA Functionals 8.3 8.6 8.3 8.0 4.5 4.5 4.7 4.8 4.7 7.7 8.0 8.6 7.3 7.7 7.7 8.8 6.4 6.6 Hybrid Functionals (Less than 40% Exact Exchange) 2.6 2.7 2.6 11.1 3.9 3.8 3.4 3.6 3.4 8.0 3.1 3.7 3.3 3.4 3.2 6.9 4.3 5.0 3.6 3.7 3.5 11.3 4.7 4.3 6.3 6.6 6.5 7.1 5.0 5.0 4.5 5.0 4.9 9.0 4.9 4.5 Hybrid Functionals (More than 40% Exact Exchange) 1.9 1.8 1.7 11.2 3.6 3.6 1.9 1.9 1.9 12.6 5.6 4.6 2.4 2.6 2.5 9.7 4.0 3.4 7.6 8.5 8.7 16.4 12.3 9.3 Double Hybrid Functionals 2.7 2.9 2.9 21.0 3.5 3.6 2.6 2.8 2.8 20.3 3.1 3.6 2.3 2.6 2.6 26.9 5.4 3.1 Wave Function Theory 8.8 201.4 188.6 186.9 4.3 35.5 7.7 12.8 2.7 62.5 33.8 24.5 0.3 51.9 19.1 9.7 0.3 9.6

emp 5Z

A

B

5Z

120.3

121.4

120.9

4.7

25.7 15.6 31.3

26.3 15.9 32.1

25.4 15.2 30.8

3.3 3.3 3.4

4.2 8.4 6.9

4.6 8.9 7.3

4.3 8.6 7.0

3.5 3.1 4.1

3.9 3.7 4.8 4.5 5.0 5.1

3.8 4.2 5.4 3.7 6.1 6.2

4.0 3.6 4.7 4.3 4.7 4.6

1.7 1.9 1.5 2.3 4.3 3.0

4.0 4.9 3.6 9.0

4.0 4.7 3.7 9.3

3.6 4.8 3.4 9.2

1.2 1.6 2.2 4.3

4.2 4.4 2.8

4.0 3.6 4.9

3.5 3.5 3.1

1.5 1.5 1.3

186.6 16.0 20.8 5.9 3.3

6.2 2.9 1.6 0.3 0.3

a

Reference data are taken from highly accurate calculations at the CCSD(T)(full) level, extrapolated to the complete basis set limit and with all electrons correlated. The three column blocks refer to results that include, respectively: “with BSR corrections”, BSR corrections toward the CBS limit of CCSD(T)(full); “regular”, no corrections at all; “emp”, empirical bond terms. bCorrelation-consistent basis sets cc-pVXZ are designated DZ, TZ, QZ, 5Z. Basis set A refers to 6-311+G(3df,2p); basis set B refers to aug-pc-2. cWave function results (MP2, CCSD, CCSD(T)) consider electron correlation only at the frozen-core level. dWave function results reported in italics refer to truncated basis set treatments as detailed in the text. eModel B3 results have arbitrarily been reported in the “5Z” column; the ATOMIC/B3 protocol is defined as BSR-corrected model B3 (left).

protocol implements BSRs in an ab initio fashion, it uses highlevel atomization energies precomputed for the parent molecules and combines them with lower-level calculations of the bond separation energies to obtain accurate estimates of atomization energies for larger molecules (here: for oxirane). High-level atomization energies for prototypes are based on accurate complete-basis set extrapolations of CCSD(T) energies with all electrons correlated and include further corrections to account for higher-order correlation effects (up to CCSDTQ), for scalar relativistic effects, and for diagonal Born−Oppenheimer terms. So far the protocol has been implemented for neutral, closed-shell molecules containing the atoms H, C, N, O, and F, and application to any of these molecules only requires that a (unique) Lewis valence structure can be drawn. Various “low-level” composite wave function methods have been suggested for use with ATOMIC, which reproduce reference atomization energies to high accuracy if BSR corrections are applied. The currently preferred composite model, ATOMIC/B3, requires MP2 calculations with up to quadruple-ζ, CCSD

calculations with mixed triple/double-ζ, and (T) evaluations with double-ζ basis sets. Although standardly applicable to molecules with up to 10−20 non-hydrogen atoms, it is still fairly expensive for larger molecules. Here we wish to explore the possibility of replacing composite wave function approaches with computationally less-demanding but still reasonably accurate approaches. Density-functional theory appears to be an attractive choice, in particular because many of the functionals proposed in the past decade have shown exceptional promise for thermochemical applications. Only a few systematic assessments of density functional theory 4 6 −4 9 (and of composite wave function approaches)47,50,51 for BSRs have been reported, although various other types of isodesmic reactions have regularly been used in computational thermochemistry.52 The most extensive study has been published by Raghavachari et al.,46 who assessed seven density functionals for a total of 40 molecules containing CC single, double, and triple bonds, CO single and double bonds, CN single bonds, as well as CH, OH, and NH bonds, concluding, in particular, that B3LYP53 shows promising 229

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

(4.2) application of BSR corrections, focusing both on differences between functionals and basis sets, and comparing to wave function theory. Regular BSRs will then be contrasted with empirically calibrated systems of bond corrections to analyze the poor performance observed for some density functionals (4.3). Section 4.4 highlights the benefit of empirical dispersion terms in the context of BSR correction schemes.

potential when used with BSRs. Other studies have been restricted to hydrocarbons or to fewer density functionals. The limitation to few bond types arises from the necessity of having accurate experimental reference data available for all relevant BSR parent molecules. The use of high-level ab initio reference data as implemented in the ATOMIC protocol lifts this restriction and allows for a much more comprehensive assessment of density functionals, which is the purpose of the present work. Density functionals are usually classified according to the level of detail with which the unknown exchange−correlation functional is approximated. Using the biblical metaphor of Jacob’s ladder,54 Perdew and Schmidt term functionals of the local (spin) density as first-rung (local density approximation, LDA) and those depending also on density gradients as secondrung (generalized gradient approximation, GGA). Third-rung functionals further include a dependency on either Laplacians of the spin density or on kinetic energy densities (meta-GGA), whereas those of the fourth and fifth rung are nonlocal and include also exact exchange (hybrid functionals) or both exact exchange and exact partial correlation, respectively. We have selected 20 density functionals from all five rungs for our study, focusing on commonly used ones and those reported to be accurate in thermochemistry, among them firstprinciple functionals as well as those with a certain number of empirically calibrated parameters. Specifically we considered one standard first-rung functional (SVWN355,56) and three popular second-rung functionals (BP86,56−58 BLYP,58,59 PBE60), the latter also augmented with exact exchange (PBE061). One early semiempirical (VSXC62) and one nonempirical (TPSS63) meta-GGA have been chosen along with the hybrid version of the latter (TPSSh64) for their reported accuracy64−66 in reproducing atomization energies and heats of formation. The extremely popular three-parameter hybrid functional B3LYP53 has been selected, along with its sibling B3PW91,67 because it has become the de facto standard in much of applied quantum chemistry. We further include some of the more recent highly parametric hybrid functionals that have been extensively calibrated for applications in thermochemistry and kinetics (BMK,68 M05,69 M0670), sometimes specifically for main group elements (M05-2X71 and M06-2X,70). The latter four functionals belong to the Minnesota family, which also encompass the meta-GGA M06L66 calibrated with a similar focus, and the hybrid functional M06-HF72 developed for spectroscopic applications. We included these for comparison. Finally, we considered two highly accurate semiempirical double hybrid functionals (B2PLYP,73 B2GP-PLYP74), one (B2-PLYP-D75) with additional empirical dispersion corrections, which include both exact exchange and second-order correlation based on Kohn−Sham orbitals and may thus be considered special cases of fifth-rung functionals. All considered density functionals are listed in Table 1. The present work examines the suitability of these 20 density functionals for atomization energies and derived heats of formation based on bond separation reactions as implemented in the ATOMIC protocol. Two reference data sets are used that encompass a large variety of neutral, closed shell molecules containing the atoms H, C, N, O, and F, and a total of 34 different bond types.76 Section 2 outlines the theoretical concept and introduces the two reference data sets used in this study. Computational details are given in section 3. Section 4 discusses results for density functionals before (4.1) and after

2. THEORETICAL CONCEPT AND REFERENCE DATA SETS The ATOMIC protocol42 uses geometries optimized at the RIMP2(fc)/cc-pVTZ level of theory,77,78 and it is based on a decomposition of reference “bottom-of-the-well” atomization energies, Eref A,e, obtained at the complete-basis set limit of CCSD(T)(full), into eight components. Two of these components are Hartree−Fock terms, three valence-shell electron-correlation terms and three terms represent core− core and core−valence electron-correlation effects. An approximation or composite model is obtained by reducing basis set requirements for one, several, or all eight components. The ATOMIC protocol then uses bond separation reactions to obtain atomization energies, as demonstrated for the simple example of oxirane (see Introduction), ATOMIC/model EA,e (C2H4O) model ref model = EA,e (C2H4O) + EA,e (C2H6) − EA,e (C2H6) ref model ref + 2EA,e (CH3OH) − 2EA,e (CH3OH) − 2EA,e (CH4) model ref moded + 2EA,e (CH4) − EA,e (H 2O) + EA,e (H 2O)

(1)

and chooses basis set requirements of the model such that errors do not exceed a certain limit for each of the eight components, and for the total atomization energy. A hierarchy of three levels was suggested, which, after BSR correction, reproduce reference atomization energies to within 0.1 kcal/ mol (A), 0.3 kcal/mol (B1,B2,B3), and 1.0 kcal/mol (C) for a set of 73 neutral closed-shell molecules with at least 3 heavy atoms (C, N, O, F) and hydrogen. This reference data set covers a wide range of different bond types and will henceforth be referred to as “AE set” (atomization energies). In a second step these atomization energies are further refined to include first-order spin−orbit corrections, which for closed-shell molecules only arise from the atoms, and corrections “beyond CCSD(T)”, i.e., terms representing relativistic effects (MVD, one-electron mass−velocity−Darwin terms, EMVD A,e (M)), the departure from the Born−Oppenheimer approximation (DBOC, diagonal Born−Oppenheimer correcDBOC tion, EA,e (M)), and higher order electron correlation corrections (CCSDTQ-CCSD(T), ECC A,e (M)). All these corrections are derived from a trivially simple model that assumes thermoneutral bond separation reactions and, thus, according to theory analogous to eq 1, does not require any calculation for the molecule under consideration. Instead corrections are composed entirely of precalculated values for the BSR parent molecules. This model works very well for MVD terms, which are local in nature, is acceptable for higher order electron correlation correction terms and typically accounts for at least some of the generally much smaller DBOC effects. Further details are given in ref 42. Standard heats of formation (T = 298.15 K) are then obtained from the negative CCSD(T) atomization energies through the addition of (a) minus the above-mentioned 230

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

that applications are not limited to cases for which accurate experimental data exist. Second, ATOMIC applies the concept of bond separation reactions to atomization energies rather than to heats of formation. This implies that no BSR corrections to zero-point energies or thermal enthalpies are taken into account.

corrections, (b) zero-point vibrational energies, ZPE(M), and differences of thermal enthalpies between 0 K and T, ΔH0(M,T), both computed from scaled harmonic frequencies, (c) atomic heats of formation, ΔH0f,k(M,0) for 0 K, and (d) negative elemental thermal enthalpy differences, −ΔH0k(M,T); see ref 43 for details. ΔHf0,ATOMIC/model(M,T ) =

ATOMIC/model (M) −EA,e

−

SO EA,e (M)

−

3. COMPUTATIONAL DETAILS Single point energy calculations were performed for all molecules contained in the AE and HOF test sets and all relevant parent molecules and simple hydrides appearing in the respective bond separation reactions. Six different basis sets81−83 have been considered, including correlation-consistent basis sets cc-pVXZ (X = D, T, Q, 5),78 which were used mainly to assess basis set convergence, 6-311+G(3df,2p),84,85 an augmented triple-ζ quality basis set often used in density functional calculations,46,79,86 and aug-pc-2,87−89 another one specifically designed for density functional theory. Tables report basis sets as DZ, TZ, QZ, 5Z (cc-pVXZ), A (6-311+G(3df,2p)), and B (aug-pc2) for brevity of notation. Availability of particular density functionals varies between different quantum chemistry packages, we have used the Gaussian90 and Turbomole91,92 (version 5.10) suites of programs, and where possible ran calculations with both to check implementation. Some functionals come in different flavors;93,94 we have used, in particular, VWN-III56 for local correlation in both the SVWN functional (denoted as SVWN3) and the B3LYP functional (Gaussian implementation), and VWN-V56 in the BP86 functional, which is implemented in Turbomole and available as “BVP86” in Gaussian. The unrestricted Kohn−Sham formalism was used for atoms. All energy calculations were performed for geometries previously optimized at the RI-MP277,95/cc-pVTZ level of theory,42 and all further data necessary for the present evaluation were taken from our earlier studies.42,43 This includes reference-level energies for BSR prototypes, and in the case of heats of formation, additionally zero-point energies and thermal enthalpies obtained from RI-MP2/cc-pVTZ calculations, as well as experimental thermochemical data for the atoms, and estimates of thermal enthalpies of n-alkanes arising from conformational averaging. Following the specification of the ATOMIC protocol, all fully BSR-corrected heats of formation additionally include spin−orbit coupling terms for the atoms as well as BSR terms computed for higher-order electron correlation corrections, scalar relativistic, and diagonal Born−Oppenheimer correction terms, all taken from our previous study. For comparison we also report composite model B3,42 Hartree−Fock, and correlated wave function (MP2(fc), CCSD(fc), CCSD(T)(fc), fc = frozen core approximation) results taken form our earlier work.42 Unfortunately some of the molecules in the HOF set had been too large to be treated at the CCSD/cc-pV5Z, CCSD(T)/cc-pVQZ, or CCSD(T)/ccpV5Z levels. Hence basis sets were consistently truncated in the spirit of composite methods, replacing individual contributions to the atomization energy as follows: X = T for both X = Q and X = 5 in (CCSD(T)/ccpVXZ − CCSD/cc-pVXZ), and X = Q for X = 5 in (CCSD/cc-pVXZ − MP2/cc-pVXZ). Using the nomenclature introduced earlier,42 CCSD/cc-pV5Z, CCSD(T)/cc-pVQZ, and CCSD(T)/cc-pV5Z are thus replaced by [5 5 4 - | - - - -], [4 4 4 3 | - - - -], [5 5 4 3 | - - - -], respectively. These approximations have fairly little effect on results as none of them ever changed any of the important error statistics

MVD EA,e (M)

DBOC CC (M) − EA,e (M) + ZPE + ΔH 0(M,T ) − EA,e atoms

+

∑ k∈M

0 (ΔHf,k (M,0) − ΔHk0(M,T ))

(2)

ATOMIC heats of formation have been compared42 to experimental data for the subset of neutral, closed-shell species of the G3/99 benchmark,79 excluding molecules with two nonhydrogen atoms or less and discarding also tetrafluoroethylene and carbonyl fluoride, for which experimental reference data are likely to be in error.42 This data set covers 103 mostly organic molecules and will henceforth be referred to as HOF (heat of formation) set. Here we simply replace the composite wave function model encoded in the superscript “model” by a particular combination of density functional and basis set. The assessment of BSRcorrected density functional calculations for the AE and HOF sets is thus performed in complete analogy to the previous one for wave function models.42 Results will be contrasted with those obtained from regular density functional calculations, which do not contain any of the BSR corrections toward CCSD(T) (eq 1) or beyond CCSD(T) (eq 2, −ESO A,e (M) − DBOC CC EMYD (M) − E (M) − E (M)). Omission of the latter is A,e A,e A,e standard practice in density functional studies, although spin− orbit coupling terms are occasionally considered80 and one might argue in favor of both MVD and DBOC,74 whereas the term −ECC A,e (M) reflects only the particular partitioning of correlation energies in ATOMIC. To remain consistent, none of the “regular” results reported herebe they density functional or wave function for comparisoncontain any of these terms. The AE set had been assembled with a focus on diversity of bond types42 but was necessarily limited to fairly small molecules with up to 6 non-hydrogen atoms (C, N, O, F) to be able to generate highly accurate theoretical reference data. It includes many examples with several heteroatoms (N, O, F) and with multiple bonds. The HOF set as subset of the standard G3/99 benchmark had been chosen to validate the ATOMIC protocol for heats of formation; it contains accurate experimental data and is thus somewhat biased toward typical organic molecules. Molecules are twice as large on average as in the AE set.42 Although the two sets share a total of 31 common molecules, they are still different enough to be considered independent benchmarks. Any differences observed in statistics for the two data sets point at differences in performance for larger vs smaller molecules and those for typically organic molecules vs those including many heteroatoms and multiple bonds. To ease comparison with earlier work,42 we have not eliminated any common molecules from either test set. Finally we note that the implementation of bond separation reactions in ATOMIC differs from their traditional use in two aspects. First, ATOMIC uses high level ab initio rather than experimental data for reference. This has the obvious advantage 231

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

(mean signed, mean unsigned, rms, maximum) by more than 9% for the AE set, for which untruncated calculations had been possible throughout. Percentile deviations may be somewhat larger after BSR correction but are negligible in absolute terms (up to 0.2 kcal/mol for maximum errors, up to 0.1 kcal/mol for all other statistical measures). To ease notation in text, tables, and figures, we still refer to all CCSD and CCSD(T) calculations as if they had been obtained without basis set truncation.

4. RESULTS A total of 73 atomization energies (AE set) and 103 heats of formation (HOF set) have been calculated both with and without BSR corrections applied, each for a total of 120 combinations of density functionals and basis sets. Considering the large number of data generated, the analysis must focus on a few important statistical parameters. Figures 1 and 2 present

Figure 2. Normal distributions of the errors (kcal/mol) for the HOF set. See the caption of Figure 1 for details. BSR-corrected results also include spin−orbit coupling terms for the atoms, and BSR corrections for terms beyond CCSD(T), i.e., MVD, DBOC, and CCSDTQCCSD(T).

⎡ 1 sn = ⎢ ⎢⎣ n − 1

ρ(e) =

normal distributions ρ(e) of errors for all density functionals using the largest of the considered basis sets (cc-pV5Z), thus essentially permitting assessments of intrinsic accuracies, as residual basis set errors are expected to be minor at the DFT level of theory. Normal distributions are calculated from mean signed errors e ̅ and corresponding standard deviations sn 1 n

n

∑ ei i=1

∑ (ei − e ̅ ) i=1

⎤1/2

2⎥

⎥⎦

⎡ 1 1⎛e − exp⎢ − ⎜ ⎢⎣ 2 ⎝ sn sn 2π

(4) 2⎤ e̅ ⎞ ⎥ ⎟ ⎠ ⎥⎦

(5)

Tables 1 and 2 focus on one statistical parameter instead (rms errors), but they show results for all basis sets analyzed. A complete set of statistical data is available as Supporting Information (Table S1). Note that Figures 1 and 2 and Tables 1 and 2 also list wave function results for comparison. 4.1. Density Functional Theory without BSR Corrections. Assessments of density functional theory for atomization energies and heats of formation abound in the literature. Here they serve the purpose of setting a reference against which BSR-corrected results may be judged, so we keep discussion relatively short. We shall initially focus on differences between functionals, normally using the largest basis set (cc-pV5Z) as reference, and then comment on basis set effects. When discussing signed errors (individual or mean signed), recall that the AE and HOF sets collect atomization energies and heats of formation, respectively. Hence negative errors indicate underbinding if reported for the AE set and overbinding if reported for the HOF set. 4.1.1. Comparison of Functionals. Results collected in Tables 1, 2, S1 (Supporting Information), and Figures 1 and 2 confirm the well-known excessive overbinding tendency of the local-density approximation, on a level comparable to the underbinding of Hartree−Fock. GGA functionals fare much

Figure 1. Normal distributions of the errors (kcal/mol) for the AE set. Dashed curves represent uncorrected results; solid curves represent BSR-corrected results. The cc-pV5Z basis has been used for all density functionals, HF, and MP2, whereas partially truncated basis sets have been used for CCSD and CCSD(T).

e̅ =

n

(3) 232

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

Table 2. Heats of Formation (HOF): RMS Errors (kcal/mol)a,b fully BSR corrected DZ

TZ

QZ

10.1

8.1

7.4

BP86 BLYP PBE

5.5 4.6 5.8

4.3 3.8 4.5

3.9 3.6 4.0

VSXC TPSS M06-L

9.9 4.6 6.0

10.8 3.9 5.6

10.3 3.8 5.3

B3LYP B3PW91 PBE0 TPSSh M05 M06

4.0 4.6 4.6 4.2 5.2 4.3

3.0 3.6 3.4 3.5 4.5 3.3

2.9 3.4 2.9 3.4 4.0 2.8

BMK M05-2X M06-2X M06-HF

4.0 3.1 3.0 5.3

2.3 1.8 2.0 6.4

2.2 1.7 2.0 5.8

B2-PLYP B2-PLYP-D B2GP-PLYP

4.0 3.9 4.0

2.9 2.6 2.9

2.4 2.0 2.4

HF MP2c CCSDc,d CCSD(T)c,d model B3e

7.4 5.8 4.2 3.2

7.4 5.4 3.0 1.4

7.6 5.2 2.9 1.0

SVWN3

5Z

regular A

B

DZ

TZ

QZ

LDA Functionals 7.2 187.0 200.5 202.2 GGA Functionals 3.8 3.9 3.8 17.2 29.8 31.3 3.7 3.7 3.8 29.0 15.7 14.6 3.7 3.9 3.6 23.9 33.1 34.3 Meta-GGA Functionals 10.3 10.5 10.0 13.5 4.9 4.6 3.9 3.8 3.9 10.9 5.9 6.6 4.7 5.4 5.4 14.2 6.6 6.2 Hybrid Functionals (Less than 40% Exact Exchange) 3.0 2.8 3.0 23.0 8.0 6.5 3.3 3.3 3.3 13.8 4.5 5.2 2.8 2.8 2.7 9.1 9.4 10.6 3.5 3.4 3.5 15.6 4.5 4.3 3.7 3.9 3.8 9.2 5.9 4.5 2.4 3.0 2.8 14.6 4.6 3.4 Hybrid Functionals (More than 40% Exact Exchange) 2.3 2.1 2.2 18.2 4.1 3.8 1.8 1.7 1.8 17.4 4.2 4.7 2.0 2.0 2.0 15.7 3.8 3.4 6.1 6.2 7.2 23.4 7.9 12.8 Double Hybrid Functionals 2.3 2.5 2.7 38.4 5.5 3.3 1.9 2.3 2.6 35.1 2.6 3.3 2.2 2.6 3.1 44.0 5.6 2.7 Wave Function Theory 7.8 326.1 308.7 306.6 4.9 67.1 12.9 14.3 3.1 102.6 50.1 34.5 1.1 88.4 30.5 14.8 0.9 7.1

7.4

emp 5Z

A

B

5Z

201.5

202.4

202.8

2.9

30.7 15.1 33.2

30.7 16.0 33.7

29.7 16.5 32.0

2.2 2.5 2.1

4.9 6.3 6.4

5.1 6.3 6.8

4.6 6.6 7.2

2.5 2.5 2.7

6.8 5.2 10.4 4.5 4.0 4.1

7.3 5.3 10.6 3.7 5.1 5.1

7.9 5.0 10.0 4.4 6.1 3.8

1.8 1.9 1.5 2.2 3.2 2.1

4.9 4.4 3.4 10.3

4.2 5.3 3.7 6.7

3.8 4.4 3.5 6.5

1.7 1.6 1.9 3.2

3.5 5.1 3.6

8.1 4.9 8.8

3.6 3.8 3.3

1.3 1.1 1.2

306.1 18.9 28.6 9.0 2.6

4.7 2.2 1.6 0.8 0.8

Reference data are taken from experiment. The three column blocks refer to results that include, respectively: “fully BSR corrected”, BSR corrections toward the CBS limit of CCSD(T)(full); “regular”, no corrections at all; “emp”, empirical bond terms. Results reported under “fully BSR corrected” and “emp” also include spin−orbit coupling terms for the atoms, and BSR corrections for terms beyond CCSD(T), i.e., MVD, DBOC, CCSDTQ-CCSD(T). bCorrelation-consistent basis sets cc-pVXZ are designated DZ, TZ, QZ, 5Z. Basis set A refers to 6-311+G(3df,2p); basis set B refers to aug-pc-2. cWave function results (MP2, CCSD, CCSD(T)) consider electron correlation only at the frozen-core level. dWave function results reported in italics refer to truncated basis set treatments as detailed in the text. eModel B3 results have arbitrarily been reported in the “5Z” column; the ATOMIC/B3 protocol is defined as BSR-corrected model B3 (left). a

hybrid functionals B3LYP and B3PW91 for small molecules (AE set) that certainly contributed to the popularity of B3LYP,97 (c) the comparably poor performance of B3LYP for larger molecules (HOF set) that has more recently been criticized,47,79,98,99 (d) the good statistics of TPSS, and to a lesser degree of TPSSh, for large molecules (HOF set) relative to that for small molecules (AE),64,65 (e) the generally good performance of some highly parametrized modern hybrid-metaGGA functionals such as M05-2X, M06-2X, and BMK,49,70,100 and (f) the potential superiority of double-hybrid functionals (B2-PLYP, B2GP-PLYP) compared to the related “simple” hybrid B3LYP,74,100 provided that basis sets are chosen sensibly (see below). One result at variance with the literature concerns Zhao’s and Truhlar’s M06-HF functional, which shows rather large errors. It contains full Hartree−Fock exchange and has been developed with spectroscopic applications in mind; so its fairly poor thermochemical performance would come at no surprise. However, the authors report that M06-HF is more accurate for atomization energies than B3LYP,72 which is

better although they still show prohibitively large residual errors of around 15−30 kcal/mol. Both BP86 and PBE overbind consistently, whereas BLYP tends to overbind molecules with many heteroatoms and underbinds many typical organic molecules (individual results not shown). Both types of molecules are represented in either test set (AE and HOF), but to different extents. This leads not only to overall smaller rms errors for BLYP compared to BP86 and PBE but also to the observed difference in behavior for the two test sets (overbinding for AE, underbinding for HOF, Table S1 (Supporting Information) and Figures 1 and 2). Third to fifth-rung, i.e., meta-GGA to double hybrid, functionals all perform significantly better than simple GGA functionals. Typically, but not consistently, there is an increase in accuracy with every rung up on Jacob’s ladder. The analysis of Tables 1 and 2 and Figures 1 and 2 confirms several observations previously reported in the literature, including (a) the outstanding performance of VSXC among the nonhybrid functionals,66,86,96 (b) the reliability of the early semiempirical 233

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

4.2.1. Comparison of Functionals. Data collected in Tables 1 and 2 and Figures 1 and 2 show promising results for the two simplest approximations of density functional theory, LDA and GGA. Application of BSR corrections cuts root-mean-square errors by factors of 3−10 for the AE set and by even more for the HOF set. Error distributions are much narrower and centered at values significantly closer to zero, now almost fully contained within the −15 to +15 kcal/mol range displayed in Figures 1 and 2. The improvement observed for LDA (SVWN3) compares to that observed for HF, and GGA functionals show an improvement similar to that observed for MP2. Like for wave function theory, we further notice a significantly reduced method dependence. LDA is now significantly closer in performance to GGA overall, and all three GGA functionals now show the same modest tendency of overbinding for the AE set and an essentially balanced performance for the HOF set with mean signed errors between only −1 and +1 kcal/mol. Recall that BLYP behaved quite differently from BP86 and PBE in the absence of BSR corrections (see above). Of course, residual rms errors of 3−4 kcal/mol (HOF set) and 5−6 kcal/mol (AE set) preclude the use of BSR-corrected GGA calculations for accurate thermochemistry, and functionals of higher than second rung are significantly more interesting. Unfortunately, though, statistical analysis draws a quite discouraging picture here. Improvement through BSR corrections is not nearly as systematic for higher-rung functionals, and in fact, several of them, including all metaGGAs, M05, and M06-HF, do not show any advantage over GGA if BSR corrections are applied. Furthermore, statistics for the AE set seems to indicate an increased rather than decreased method dependence for BSR-corrected calculations at the meta-GGA or higher level. There are only six density functionals of third or higher rung that benefit consistently and strongly from the application of BSR corrections. One is the nonempirical meta-GGA TPSS, for which residual rms errors unfortunately remain on a fairly high level. PBE0 and B3LYP both show particularly encouraging results for the HOF set, but using larger basis sets, they also benefit from BSR corrections for the AE set. The remaining three functionals are the hybrid-meta GGAs BMK, M05-2X, M06-2X, which all attain a respectable accuracy with rms errors of around 2 kcal/mol after BSR correction. M05-2X benefits from a consistently 2-fold or higher cut in rms errors for both benchmarks, irrespective of basis set. Four functionals (B3PW91, TPSSh, M06, M06-HF) show little to no improvement through BSR corrections for the AE set and at best moderate improvements for the HOF set. The unsatisfactory performance of B3PW91 compared to B3LYP may surprise, but it has in fact been observed in an earlier study using bond separation energies based on experimental reference data for a small set of typical organic molecules.46 The double hybrid functionals take a special position in that BSR corrections mainly seem to reduce basis set effects. Hence significant improvements are seen for some basis sets, none or even slight deteriorations for others. However, BSR-corrected large-basis set results are always better than the best (not necessarily large basis set) uncorrected results. Two of the remaining three functionals (M06-L, M05) show a small but consistent increase in rms error for the AE set and only very minor improvements for the HOF set. The biggest surprise is the excessively poor performance of the VSXC functional, especially because it is the best among the

clearly not confirmed here. Closer inspection of individual results indicates that the discrepancy is at least partly due to some very large errors for individual species not included in the benchmark of Zhao and Truhlar, such as O3, N2O3, and F2O2 (all >20 kcal/mol for cc-pV5Z basis set). 4.1.2. Basis Set Effects. Error statistics are expectedly much less basis set dependent for density functionals than they are for correlated wave function theory. Still we note that the cc-pVDZ basis set is generally too small for any density functional to provide results of useful accuracy. Even beyond double-ζ, basis set considerations may still be important, as has been reported repeatedly in the literature. Curtiss et al. have noticed a surprising dependence of density functional heats of formation on the particular set of polarization functions added to a triple-ζ (6-311+G) basis set,98 whereas Wang and Wilson have concluded that atomization energies computed with various density functionals and correlation-consistent basis sets are often not converged until the quadruple-ζ level.101 The abovenoted M06-HF functional, not considered by these authors, seems to be a fairly extreme example demonstrating such behavior, as binding increases by 7 kcal/mol (AE set) and 12.3 kcal/mol (HOF set) on average, respectively, when going from cc-pVTZ to cc-pVQZ. This increase in binding leads to a noticeable improvement of results for the AE set (mean signed error: −1.6 instead of −8.6 kcal/mol) and to an appreciable deterioration of results for the HOF set (mean signed error: −9.0 instead of +3.3 kcal/mol), clearly visible also in the rms errors reported in Tables 1 and 2. Finally, we note that doublehybrid functionals are the only functionals studied here for which rms errors below 3 kcal/mol are within reach. Unfortunately, however, best results are not obtained at the basis set limit, and basis set sensitivity is somewhat more pronounced than for other functionals, reflecting the MP2-like admixture of electron correlation. Hence judicious basis set choices need to be made and a fair degree of error cancellation be accepted if one wishes to profit from the better error statistics (Tables 1 and 2). 4.1.3. Comparison to Wave Function Theory. Perusing the data collected in Tables 1 and 2 quickly shows that, for a given basis set, any density functional beyond GGA outperforms MP2 and CCSD and typically even CCSD(T). Still the observed accuracy is hardly sufficient for reliable calculations of atomization energies and heats of formation, especially because maximum errors almost always exceed 10 kcal/mol for one or both test sets, with only two marginal exceptions: B2-PLYP-D/ cc-pVTZ (9.4 kcal/mol, dicyanoacetylene, AE) and B2GPPLYP/cc-pVQZ (9.2 kcal/mol, F2O2, AE; Table S1, Supporting Information). What makes wave function theory still usable for treating atomization processes is the observation that residual errors largely reflect the very slow but very regular convergence of electron correlation energies toward the complete basis set limit, which can be overcome by suitable extrapolation procedures.102 4.2. Density Functional Theory with BSR Corrections. In the development of the ATOMIC protocol we noted not only that BSR corrections reduced residual errors of wave function type approaches significantly but also that results became much less method and basis set dependent upon introduction of BSR corrections.42 The reduced basis set dependence was in fact the key observation that allowed us to lower the basis set requirements of more expensive components to the electron correlation energy without sacrificing overall accuracy. 234

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

nonhybrid functionals for thermochemistry not based on bond separation schemes. However, BSR corrections lead to significant overbinding (Figures 1 and 2) and thus consistently double the rms errors for all basis sets beyond cc-pVDZ. In fact, VSXC performs worse in BSR schemes than any of the GGA functionals considered and, for the HOF set, even worse than SVWN3. This issue will be further analyzed later on. 4.2.2. Basis Set Effects. The only general improvement observed for all functionals is a systematically reduced basis set dependence. In particular we observe (a) a significant reduction of any excessive error found for the small double-ζ basis set, (b) a notable reduction of the larger basis set dependence of double hybrid functionals, especially for the HOF set (see basis set A, in particular), and (c) a correction of the initially puzzling large basis set dependence of M06-HF (s.a.). Apart from these individual observations, rms error statistics seem less conclusive. Still BSR corrections lead to a notable overall reduction in basis set effects as careful analysis of all individual data (73 AE data and 103 HOF data each for 20 functionals, total of 3520 data) reveals. Discarding the small and unreliable cc-pVDZ basis set, the variation among the remaining five basis sets is less than 1.0, 2.0, and 5.0 kcal/mol, respectively, for only 5, 42, and 88% of all AE data if no corrections are applied, but for 53, 79, and 98%, if BSR corrections are applied. Values for the HOF data set confirm this picture (uncorrected 3, 28, 73%; with BSR corrections 64, 88, 98%). Finally it is noteworthy that the small remaining basis set effects generally mean reduction of rms errors with basis set extension. This tendency is fairly modest, however, and never seems to justify the computational expense of quintuple−ζ basis sets. In general, the Pople-style basis set A (6-311+G(3df,2p)) seems to be the optimal choice; it is similar in size to cc-pVTZ but yields results in quality closer to cc-pVQZ. Basis set B (aug-pc-2) is the largest of all triple-ζ quality basis sets considered here, it generally performs slightly better than A for the AE set, but not consistently so for the HOF set. 4.2.3. Comparison to Wave Function Theory. Comparison to wave function theory is favorable for the simple LDA and GGA functionals where the level of improvement achieved through BSR corrections is similar to that observed for HF and MP2 theory, respectively. Higher rungs of density functional theory often outperform even CCSD(T) in finite basis set calculations; however, residual errors are much less systematic and thus less correctable in bond separation schemes. The best density functionals still achieve an accuracy comparable to or slightly better than CCSD after BSR correction but are far behind CCSD(T). Inclusion of connected triples excitations, however, was found to be important for accurate ab initio thermochemistry. Note that the computationally most attractive BSR-corrected wave function model, ATOMIC/B3, uses fairly small basis sets for electron correlation contributions beyond MP2, and only a double-ζ basis sets for the triples correction,42 but still shows rms errors of only 0.3 and 0.9 kcal/ mol for the AE and HOF sets, respectively. BSR corrections essentially annihilate the mean signed error and significantly reduce the width of the error distribution as well. This level of accuracy seems out of reach for any current density functional, as inspection of Figures 1 and 2 readily shows. 4.3. Empirical Bond Corrections. The poor performance of VSXC and some other functionals for bond separation reactions is of some concern and deserves further analysis. Obviously errors are fairly unbalanced between larger molecules and smaller prototypes, but lack of consistency may also be

more general. Empirical calibration of bond increments, often referred to as bond additivity corrections (BAC),103,104 should provide some insight. Note that regular BSR corrections to atomization energies can be mapped exactly onto a system of corrective bond terms.42 The mapping only requires that bond terms are derived from energy differences (reference minus production level) for the BSR parent molecules; the term for CH bonds would thus be one-quarter of the energy difference for methane, and the term for CC single bonds would be the energy difference for ethane minus 6 times the CH bond term. Alternatively, one may calibrate these corrective bond terms so that they provide the best possible corrections for the entire data set of larger molecules. Then they do not represent bond separation reactions anymore, but they will lead to substantially better results if errors are imbalanced between small BSR parent molecules and larger molecules but otherwise still fairly characteristic of bond type. Empirical calibration of bond terms unfortunately meets the practical problem that a number of bond types occur for only one or a few molecules in the data set. This not only introduces undue bias, as calculations for some of these molecules would necessarily reproduce reference data exactly or almost so, but also may even lead to linear dependencies among parameters. We have thus limited calibration, somewhat arbitrarily, to those 14 bond types that occur for at least 5 molecules of either of the two data sets (AE, HOF), and constrained bond terms to be equal to their BSR value for the remaining 20 bond types. 4.3.1. Technical Aspects. The empirical model estimates the atomization energy EA,e(M) of a particular molecule M from the computed DFT value Emodel A,e (M) and the sum of all relevant bond terms, some of which are constrained to their BSR value BSR (ΔEA,e [j]) and others that are free to be optimized emp (ΔEA,e [k]). The reference function to be minimized is then written as Z=

ref model BSR (M) − (EA,e (M) + ∑j njBSR ∑ (EA,e ,M ΔEA,e [j] M

+

2

emp ∑k nkemp ,M ΔEA,e [k ]))

(6)

Eref A,e(M),

nBSR j,m ,

nemp k,m

where and denote the reference atomization energy for M, and the number of bonds of type j and k occurring in M, respectively, with index j running over all bond types with fixed BSR increments, and index k running over all bond types with optimizable empirical parameters. The optimization problem is of linear least-squares type, and solutions are easily obtained by matrix inversion (see, e.g., ref 105), emp ΔEA,e [k] = (A−1B)k

(7)

where matrix A and vector B are defined by their elements and components

(A)lk =

emp ∑ nlemp ,M nk ,M M

(B)l =

(8)

ref model ∑ nlemp ,M (EA,e(M) − EA,e (M) M

−

BSR ∑ njBSR ,M ΔEA,e [j]) j

(9)

The reference function Z encompasses all 176 data from the AE (73) and HOF (103) sets, the latter of course enter with 235

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

Figure 3. Empirical bond parameters and BSR increments. The graphs illustrate the difference between bond parameters ΔEemp A,e , obtained by linear least-squares fit, and BSR increments ΔEBSR A,e computed from atomization energies of the relevant BSR parent molecules. Least squares fits have been performed as documented in the text. The panel on the left-hand side indicates the bond types106 and the number of molecules in the (AE/HOF) sets for which they occur as well as their total occurrence in the (AE/HOF) sets. 0,ATOMIC/model negative sign (−ΔH0,ref (M) + f (M) and −[ΔHf ATOMIC/model model ref (M) − EA,e (M)] instead of EA,e (M) and EA,e Emodel A,e (M), respectively) because bond increments are defined as corrections to atomization energies. As written, density functional data for the HOF set already include the corrections “beyond the CCSD(T) level”, i.e., BSR terms for MVD, DBOC, and CCSDTQ-CCSD(T), as well as atomic spin−orbit terms, compare to eqs 1 and 2. We note that some molecules occur in both the AE and HOF sets and hence are included with two (different) reference data in Z; however, only this ensures that the optimization does not carry an artificial bias toward one of the two data sets. It should be stressed that the calibration efforts described here serve the sole purpose of exploring the possible benefits of a BAC model; they are not intended to provide a set of recommended parameters, which would require a significantly larger and more balanced reference data set. 4.3.2. Observations. Optimizations were performed for all 20 density functional and 4 wave function models using the largest basis set studied, as well as for composite model B3. Resulting rms errors for the AE and HOF sets are included in the rightmost columns of Tables 1 and 2. Inspection of the data shows that many density functionals benefit substantially from empirically calibrated bond corrections. Two functionals stand out in particular, SVWN3 and PBE0, as they already performed significantly better with BSR corrections than without and experience a further halving of rms errors through empirical calibration. Using these bond terms, the simple LDA functional is competitive with the best uncorrected hybrid functionals, whereas PBE0 is better than any of the tested functionals with or without BSR corrections. Note that both functionals are intrinsically parameter-free; hence they were not fitted to reproduce experimental data. The only empirical aspect in their use here is that of a posthoc application of calibrated bond corrections. Reassuringly empirical calibration of bond terms helps in the problematic case of VSXC. Root mean square errors are now

halved for the HOF set instead of doubled as was the case with proper BSR corrections. Improvements seems less satisfactory for the AE set (from 4.2 to 3.5 kcal/mol), but this is almost entirely due to two molecules in the AE set, which contain only bonds not considered in the calibration process due to scarcity of reference data (F2O2 and O3, residual errors 15.5 and 11.2 kcal/mol). Skipping these two molecules in the statistics reduces the rms error from 3.5 to 2.7 kcal/mol. At least for VSXC we can thus conclude that errors overall are still reasonably characteristic of bond type and that the poor performance in bond separation reactions chiefly reflects an inconsistency of error between the very smallest molecules (BSR prototypes) and larger molecules. However, there are also functionals that do not permit such a conclusion. The performance of M05, for example, seems fairly resistant to empirical bond corrections. Empirical calibration certainly improves statistics, as it should, but not much is gained compared to uncorrected results. Compared to VSXC, the AE statistics is not as much flawed by the inclusion of F2O2 and O3 (residual errors 9.0, 6.6 kcal/mol), but it rather reflects a number of highly strained molecules for which M05 performs fairly poorly (see below). In absolute terms double hybrid functionals are the strongest competitors of wave function theory, if empirical bond corrections are accounted for. For the HOF set, they exhibit residual rms errors only about 50% higher than those of CCSD(T) or carefully designed composite approaches (B3). Results seem less promising for the AE set. At this level of accuracy one needs to be careful with interpretation of differences, however, as reference data here are based on CCSD(T)(full) calculations at the complete basis set limit. Though generally very accurate, this reference level of theory has problems with typical multireference cases, and with ozone in particular.42 Nevertheless, rms errors calculated for the remaining 72 molecules (1.1−1.3 kcal/mol) still indicate that double hybrids are some way behind the best wave function 236

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

Table 3. Heats of Formation of Hydrocarbons (kcal/mol)a,b deviation (computed − reference) VSXC/5Z molecule C3H4 C3H4 C3H4 C3H6 C3H6 C3H8 C4H6 C4H6 C4H6 C4H6 C4H6 C4H6 C4H8 C4H8 C4H10 C4H10 C5H8 C5H8 C5H10 C5H12 C5H12 C6H6 C6H8 C6H8 C6H12 C6H14 C6H14 C7H8 C7H16 C8H8 C8H18 C10H8 C10H8 RMS

allene cyclopropene propyne cyclopropane propene propane 1,2-butadiene 1,3-butadiene 2-butyne bicyclo[1.1.0]butane cyclobutene methylenecyclopropane cyclobutane isobutene isobutane n-butane isoprene spiropentane cyclopentane n-pentane neopentane benzene 1,3-cyclohexadiene 1,4-cyclohexadiene cyclohexane 3-methylpentane n-hexane toluene n-heptane 1,3,5,7-cyclooctatetraene n-octane azulene naphthalene (33 data)

M05/5Z

M05-2X/5Z

reference

BSR

reg

emp

BSR

reg

emp

BSR

reg

emp

45.5 66.2 44.2 12.7 4.8 −25.0 38.8 26.3 34.8 51.9 37.5 47.9 6.8 −4.0 −32.1 −30.0 18.0 44.3 −18.3 −35.1 −40.2 19.7 25.4 25.9 −29.5 −41.1 −39.9 12.0 −44.9 70.7 −49.9 69.1 35.9

−4.4 −6.8 −2.1 −7.4 −2.1 −1.7 −4.8 −4.7 −3.3 −15.0 −6.3 −11.8 −7.5 −4.8 −5.2 −3.1 −10.0 −15.6 −6.8 −4.3 −11.9 −13.5 −10.0 −8.1 −11.5 −13.0 −5.9 −15.5 −7.3 −15.2 −8.7 −27.9 −28.0 11.1

−5.5 −0.8 −0.9 1.9 0.1 3.7 −2.8 −2.7 1.1 1.1 2.8 −2.7 4.8 0.4 3.3 5.4 −5.0 3.6 8.6 7.2 −0.4 −4.8 2.0 3.9 6.9 1.6 8.7 −3.7 10.4 −3.5 12.1 −9.5 −9.7 5.3

−1.2 −0.2 0.6 −1.4 0.4 0.1 0.4 0.5 1.5 −3.0 2.3 −3.2 0.4 −0.3 −1.4 0.7 −2.9 −1.6 3.1 1.4 −6.2 0.5 3.3 5.2 0.3 −5.3 1.8 0.5 2.4 3.5 3.0 −0.5 −0.6 2.4

−3.7 −7.7 −1.5 −5.6 −0.2 0.4 −3.3 −0.2 −2.2 −11.2 −2.2 −9.2 −1.9 0.6 2.1 1.3 1.3 −11.8 2.9 2.4 4.8 −1.5 2.3 1.5 5.1 4.8 3.1 −0.6 4.0 2.9 4.9 −1.5 −1.3 4.4

−1.1 −7.0 −0.6 −7.4 −0.1 −1.9 −1.3 1.8 −1.9 −13.1 −2.2 −9.1 −4.4 0.2 −0.8 −1.6 2.7 −14.3 −0.1 −1.1 1.3 2.3 3.6 2.8 1.4 0.7 −1.0 2.6 −0.7 7.9 −0.5 5.3 5.5 4.8

−4.3 −7.2 0.1 −6.2 −1.9 −2.4 −4.1 −1.0 −0.9 −9.9 −2.0 −8.9 −2.8 −1.3 −0.9 −1.7 0.2 −10.8 1.8 −0.9 1.6 0.5 3.2 2.4 3.8 1.4 −0.4 1.2 0.3 5.6 1.0 4.0 4.2 4.1

−1.4 1.5 −0.6 −0.6 −0.3 −0.2 −1.0 0.1 −0.6 1.5 3.3 −2.5 1.4 −0.2 0.1 0.1 0.1 −1.0 1.3 0.5 0.4 −2.2 0.9 0.3 0.8 0.7 0.5 −1.8 0.8 2.4 1.0 3.4 −2.8 1.4

−1.5 −0.1 0.8 −1.9 0.0 0.3 −1.5 −0.4 0.4 −1.9 1.3 −4.4 −0.3 −0.4 0.2 0.2 −0.8 −4.8 −0.8 0.1 0.1 −5.5 −2.1 −2.8 −1.7 0.0 −0.2 −5.5 −0.4 −2.0 −0.6 −4.0 −10.1 2.8

−0.7 1.5 −0.7 −1.0 0.1 −0.2 −0.4 0.7 −0.8 0.8 3.2 −2.6 0.9 0.0 0.0 −0.1 0.6 −1.9 0.7 0.2 0.1 −1.8 1.0 0.4 0.1 0.3 0.1 −1.4 0.3 3.0 0.4 3.8 −2.3 1.4

Reference data are taken from experiment. The three column blocks refer to results that include, respectively: “BSR”, BSR corrections toward the CBS limit of CCSD(T)(full); “reg”, no corrections at all; “emp”, empirical bond terms. Results reported under “BSR” and “emp” also include spin− orbit coupling terms for the atoms, and BSR corrections for terms beyond CCSD(T), i.e., MVD, DBOC, CCSDTQ-CCSD(T). bThe correlation consistent basis set cc-pV5Z has been used. a

Parameter shifts are consistently below 1 kcal/mol, and improvements beyond proper BSR corrections are moderate at best, indicating that BSR corrections tap the full potential of removing systematic error. M06-2X behaves similarly except for two types of bonds between N and O, whereas BMK benefits noticeably from somewhat larger, but overall still very modest, parameter shifts. 4.3.3. Illustrative Examples. It is probably instructive to analyze some individual results in more detail. We have chosen VSXC as an example, where BSR corrections fail and empirical corrections work well, M05 as a case where neither choice produces really satisfactory results, and M05-2X as an example, where already BSR corrections work well and further empiricism shows no advantage. Table 3 shows results for the 33 hydrocarbons taken from the HOF set. M05-2X already performs fairly well without any corrections, and the major remaining weakness of overbinding polycyclic and multiply unsaturated molecules, i.e., species with few hydrogens, is corrected well in BSR schemes. This good

models listed in Table 1, even if empirical bond corrections are applied. The values of all 14 optimized bond parameters ΔEemp A,e [k] are displayed in Figure 3, as deviation from regular BSR values ΔEBSR A,e [k]. It is obvious that parameter shifts are often quite substantial, exceeding a tolerable amount of perhaps 1 kcal/mol regularly. Note that the displayed range of −6.5 to +6.5 kcal/ mol is even exceeded in several cases and that many methods show considerable shifts across the entire range of bonds. This simply means that knowledge on bond errors obtained from small BSR prototypes is often not transferable to larger molecules. Accurate wave function theory (CCSD(T) and model B3) behaves quite differently in this respect. Here all sources of remaining systematic error, mainly basis set incompleteness and lack of core−core and core−valence correlation, are quite local in nature and thus easily correctable through BSR corrections obtained from small model systems. Among all density functionals considered here, M05-2X is the closest in this respect to accurate wave function theory. 237

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

Table 4. Atomization Energies of Compounds with NO Bonds (kcal/mol)a,b deviation (computed − reference) VSXC/5Z molecule CHNO CHNO CH3NO2 CH3NO2 FNO HNO2 HNO2 HNO2 HNO3 H2N2O N2O N2O3 RMS

formonitrile oxide, HCNO isofulminic acid, HONC methylnitrite nitromethane nitrosyl fluoride nitrous acid (cis) nitrous acid (trans) nitrous acid (H-NO2) nitric acid nitrosamine nitrous oxide dinitrogen trioxide (12 data)

M05/5Z

M05-2X/5Z

reference

BSR

reg

emp

BSR

reg

emp

BSR

reg

emp

363.9 349.9 598.9 601.5 214.6 311.6 312.0 304.5 388.0 386.5 269.2 390.0

11.3 2.7 11.6 10.3 9.8 8.2 8.1 8.5 13.8 8.6 12.3 24.5 11.9

5.6 −3.0 −0.2 −0.8 4.8 −0.8 −1.0 0.9 1.7 −2.0 3.4 9.2 3.7

−0.5 2.7 1.0 −0.5 −2.0 −1.0 −1.1 −4.0 1.3 −0.2 −0.2 2.9 1.8

9.3 4.9 4.8 6.2 5.2 6.2 6.2 7.5 11.7 6.0 12.6 14.1 8.5

7.1 2.4 0.2 2.1 −0.4 −0.5 −0.6 2.5 4.9 1.5 7.7 6.0 4.0

1.1 3.7 1.0 −0.8 −1.9 −0.3 −0.3 −1.3 1.7 1.0 3.8 0.0 1.8

−1.2 −0.7 0.9 1.8 −0.1 1.4 1.0 1.1 2.7 3.2 −0.4 −2.0 1.6

−4.6 −2.0 −5.2 −4.2 −6.5 −5.7 −6.2 −5.8 −4.5 −2.3 −6.2 −13.2 6.2

−1.3 −1.4 0.9 1.6 −0.5 0.3 −0.2 1.1 1.9 1.6 −0.4 −2.4 1.3

a

Reference data are taken from highly accurate calculations at the CCSD(T)(full) level, extrapolated to the complete basis set limit and with all electrons correlated. The three column blocks refer to results that include, respectively: “BSR”, BSR corrections toward the CBS limit of CCSD(T)(full); “reg”, no corrections at all; “emp”, empirical bond terms. bThe correlation consistent basis set cc-pV5Z has been used.

of the 25 methods analyzed here exhibit rms errors in excess of 2 kcal/mol after application of empirical bond terms (M06-HF, 3.2 kcal/mol and HF, 4.0 kcal/mol). This good performance is not a trivial result reflecting overparametrization, as the 12 listed compounds involve five different types of NO bond (Table 3 of ref 42), of which only two (Figure 3) were included in the empirical calibration. Five molecules contain one of these two bond types106 (3N−O: HCNO, CH3NO(O), H NO(=O), HONO(O), NNO), five the other (1NO: CH3ONO, FNO, cis- and trans-HONO, H2NNO), and one both (ONNO(O)), obviously in quite different environments. Of course, any substantial improvement for molecules in the test set comes at the expense of significantly worse results for the BSR parent molecules. The extreme parameter shifts observed for VSXC, e.g., translate into errors of −12 and −9 kcal/mol for the atomization energies of H3 NO and HNO, respectively. In summary, the analysis shows that many cases of substantial error at the BSR level still reflect to a large extent systematic error characteristic of bond type, which can be reduced significantly through empirical calibration. In these cases small BSR parent or prototype molecules are just not representative enough to correct for systematic error. This is not a particularly satisfactory conclusion theoretically, as it precludes the use of many density functionals in the strictly ab initio framework of the ATOMIC approach. However, it opens possible avenues of improvement if empirical calibration is accepted. Still some of the remaining problems are necessarily resistant to any type of bond correction, among them problems with strained or branched molecules, and, as discussed next, with large molecules. 4.4. Empirical Dispersion Terms. We have included one double-hybrid functional with (B2-PLYP) and without empirical dispersion terms (B2-PLYP-D) in our study, allowing us to assess the effect of dispersion terms in the framework of BSR-corrected thermochemistry. Grimme’s second version of the empirical dispersion model 75,107 has been chosen deliberately for this purpose, as it not only is very simple but also is arguably the most popular one currently in use and, furthermore, seems to have an advantage over the very recent

performance of M05-2X for bond separation energies of hydrocarbons was already demonstrated by Wodrich et al.49 Empirical calibration of bond corrections shifts individual values somewhat but it does not affect overall quality. Uncorrected VSXC, on the other hand, underbinds acyclic alkanes and often overbinds hydrogen-poor hydrocarbons. The BSR scheme grossly overcorrects the first problem and aggravates the second. Empirical calibration of bond terms, however, results in major improvement for almost all hydrocarbons. Still this functional does not reach the quality of M05-2X, reflecting to a large extent its failure to account for alkane branching (see, e.g., 3-methylpentane vs n-hexane), which is of course not correctable in any bond-based scheme. Finally, M05 shows a reasonable performance for most hydrocarbons but significantly overbinds highly strained small rings. This tendency was also observed for examples not among those listed in Table 3, including molecules such as aziridine (6.8 kcal/mol overbinding, AE set), 3H-diazirine (7.2, AE), oxirane (5.5, AE), and tetrahedrane (23.2, AE). Of course, also such a weakness is resistant to any bond correction scheme, thus rendering both BSR and empirical corrections fairly ineffective. Figure 3 indicates that significant parameter shifts through empirical calibration are not limited to CC and CH bonds. Particularly extreme shifts are generally found for NO bonds, where only M05-2X, CCSD(T), and model B3 appear immune to calibration. Table 4 lists results for all 12 molecules of the AE set which contain NO bonds. Overall, bond corrections are surprisingly effective. Again, M05-2X is very accurate already at the BSR level, obviating the need for empirical calibration. The reliable removal of error through BSR corrections here is quite remarkable in view of the complex valence structures that give rise to quite unusual bond separation reactions, such as N2O3 + NH3 + 2NH4 + → HNO + H3NO + H3NNH 2+ + H 2NO+

For both VSXC and M05, we observe that BSR corrections deteriorate moderately useful accuracy, but that empirical calibration of bond parameters improves the situation significantly. Additional analysis further shows that only two 238

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

Table 5. Heats of Formation of Alkanes: Empirical Dispersion Terms and BSR Corrections (kcal/mol)a deviation (computed − reference) B3LYP molecule

reference

BSR

C3H8 C4H10 C4H10 C5H12 C5H12 C6H14 C6H14 C7H16 C8H18

propane isobutane n-butane n-pentane neopentane 3-methylpentane n-hexane n-heptane n-octane

−25.0 −32.1 −30.0 −35.1 −40.2 −41.1 −39.9 −44.9 −49.9

0.6 2.8 1.8 3.1 6.2 6.0 4.1 5.3 6.5

C3H8 C4H10 C4H10 C5H12 C5H12 C6H14 C6H14 C7H16 C8H18

propane isobutane n-butane n-pentane neopentane 3-methylpentane n-hexane n-heptane n-octane

−25.0 −32.1 −30.0 −35.1 −40.2 −41.1 −39.9 −44.9 −49.9

0.6 2.8 1.8 3.1 6.0 6.0 4.0 5.2 6.3

B3LYP-D reg

BSR

cc-pVDZ Basis Set 17.4 −0.3 24.5 −0.1 23.5 −0.2 29.7 0.0 32.8 0.2 37.6 −0.7 35.7 −0.2 41.8 −0.1 47.9 −0.1 cc-pV5Z Basis Set 2.6 −0.3 6.3 −0.2 5.3 −0.3 8.1 −0.1 10.9 −0.1 12.4 −0.8 10.5 −0.3 13.1 −0.3 15.8 −0.3

B2PLYP

B2PLYP-D

reg

BSR

reg

BSR

reg

12.7 16.3 16.2 19.8 20.0 22.5 23.1 26.5 30.0

0.3 1.6 1.0 1.9 3.4 3.6 2.4 3.2 3.9

30.2 39.9 39.4 48.7 50.2 58.7 57.6 66.8 75.9

−0.2 0.0 −0.1 0.2 0.3 0.0 0.2 0.4 0.5

27.8 35.7 35.6 43.5 43.5 50.9 51.0 58.8 66.5

−2.0 −1.9 −2.0 −1.9 −1.9 −2.6 −2.1 −2.1 −2.1

0.1 1.3 0.8 1.5 2.9 3.0 1.9 2.5 3.1

0.3 1.6 1.1 2.1 3.4 3.7 2.7 3.5 4.3

−0.4 −0.3 −0.3 −0.1 −0.3 −0.6 −0.3 −0.3 −0.3

−2.1 −2.7 −2.7 −3.1 −3.3 −4.2 −3.9 −4.5 −5.1

Column designations (BSR for “fully BSR corrected” and regular) are defined in footnote a of Table 2. B3LYP-D data are not reported elsewhere in this paper; they have been generated from B3LYP calculations through adding empirical dispersion terms using the recommended scaling factor for B3LYP of s6 = 1.05 (see text). a

but also shift them toward the “right” origin. And, indeed, we observe that, after BSR correction, B2-PLYP-D always performs slightly better than B2-PLYP, no matter which basis set is used. This comparison clearly shows that BSR corrections largely remove basis set effects, which otherwise obscure the true benefit of including empirical dispersion terms. Being local in nature, BSR corrections themselves cannot properly account for dispersion even if the employed reference level of theory does. However, they do modify empirical dispersion terms just as they affect any other contribution to the total energy of a molecule. These modifications can become fairly large because empirical dispersion terms (Edisp) in the model of Schwabe and Grimme75 are non-negligible for the small parent molecules that enter the bond separation reaction and they have no counterpart at the reference level of theory. For example, we find BSR dispersion contributions of +5.9 kcal/mol for n-octane (=6Edisp(CH4) − 7Edisp(C2H6)) and of +5.5 kcal/mol for naphthalene (=12Edisp(CH4) − 5Edisp(C2H4) − 6Edisp(C2H6), results not shown in detail), which reduce the dispersion term of n-octane from −9.4 to −3.4 kcal/mol and essentially annihilate that of naphthalene (from −5.4 to +0.2 kcal/mol). The latter case shows that dispersion terms may actually lead to increased molecular enthalpies if considered in the framework of BSR-corrected models. An examination of mean signed errors (Table S1, Supporting Information) indeed shows that this is the rule rather than the exception for the rather small molecules contained in the AE set. In summary one should not attribute too much physical meaning to BSRcorrected dispersion terms and simply note that they serve the purpose of improving results overall. Staying with the two examples mentioned above and looking at BSR-corrected B2PLYP employing the largest basis set cc-pV5Z, we note that empirical dispersion terms correct the overestimated heat of formation of n-octane by just the right amount (B2-PLYP,

third version for thermochemical problems involving carbon atoms.108 4.4.1. AE and HOF Sets. Although calibrated with reference data for noncovalent binding energies,75 the dispersion model also showed remarkable success in improving heats of formation in the G3/99 benchmark, with mean absolute deviations decreasing from 2.4 to 1.7 kcal/mol.75 Tables 2 and S1 (Supporting Information) confirm this observation for the HOF data set, which is a subset of G3/99 including only neutral closed-shell (HCNOF) molecules and, expectedly to a lower extent, for the AE set, which focuses on smaller molecules. It is also obvious, however, that the observed improvement is limited to small basis sets of up to triple-ζ quality for which B2-PLYP shows a consistent tendency of underbinding (Table S1, Supporting Information). Previously, Grimme had reported good thermochemical performance for regular B2-PLYP without dispersion terms, using a basis set of quadruple-ζ quality.73 Mean signed errors collected in Table S1 (Supporting Information) demonstrate that B2-PLYP underbinds for cc-pVTZ (AE, by 1.4 kcal/mol; HOF, by 4.1 kcal/ mol), is about balanced for cc-pVQZ and overbinds for ccpV5Z (by 2.4 and 1.5 kcal/mol; see also Figures 1 and 2). Because empirical dispersion terms will always lower the molecular energy, results are less promising at the cc-pV5Z level, where rms errors are larger for B2-PLYP-D than for B2PLYP (Tables 4 and 5). Root mean square errors only tell half the story, however, as they convey no information on error distributions. A closer look at Figures 1 and 2 and Table S1 (Supporting Information) shows that, irrespective of basis set, error distributions always get narrower when empirical dispersion terms are included. BSR corrections cancel or at least reduce systematic deficiencies in the description of bond breaking, and so they typically not only narrow down error distributions even further 239

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

−46.8; B2-PLYP-D, −50.2 kcal/mol; exp, −49.9 kcal/mol), whereas they do not deteriorate the already underestimated heat of formation of naphthalene (B2-PLYP, 28.2; B2-PLYP-D, 28.4 kcal/mol; exp, 35.9 kcal/mol). 4.4.2. Alkanes. It is instructive to analyze results obtained for alkanes in some more detail (Table 5). Because the popular B3LYP functional has often been criticized for its poor performance for alkanes,47,109,110 we have generated B3LYPD data, using damping function parameters calibrated for this particular functional.107 Near the basis set limit (cc-pV5Z), regular B3LYP shows the repeatedly reported increase of error with molecular size. B2-PLYP performs much better in this respect, because the perturbation term partially accounts for medium-range correlation, but it is still far from being perfect. BSR corrections reduce each individual error, still they cannot remove the size dependence completely. On the other hand, addition of empirical dispersion terms always lowers computed heats of formation, but for the cc-pV5Z basis set it obviously carries the “correction” to excess. Consideration of both, empirical dispersion terms and BSR corrections, yields nearperfect agreement with experiment for either functional. In the particular case of n-alkanes BSR corrections appear to eliminate basis set errors extremely well so that the same conclusion would be reached for all other basis sets studied, including the smallest one (cc-pVDZ, Table 5), which cannot normally be recommended for applications. The data further show that relative enthalpies between linear and branched isomers (isobutane vs n-butane, neopentane vs n-pentane, and 3methylpentane vs n-hexane) are much improved upon consideration of empirical dispersion terms. With either of the two basis sets and BSR corrections applied, residual errors do not exceed 0.3 and 0.5 kcal/mol for B2-PLYP-D and B3LYP-D, respectively, whereas they reach up to 1.5 and 3.1 kcal/mol, respectively, for B2-PLYP and B3LYP. Both types of error may be related to the common failure of many currently used density functionals to properly account for “protobranching” stabilization,111 i.e., for the favorable interaction between (CHn) groups in 1,3-position. There is some debate about the origin of this failure, with explanations rooting in either medium-range correlation,110 intramolecular dispersion,109 or long-range exchange effects.112,113 Here it suffices to say that the inclusion of empirical dispersion terms as defined in the model of Schwabe and Grimme75 provides for a very accurate description of linear and branched alkanes, provided that BSR corrections derived from accurate reference data are applied. The exceptional agreement with experiment seems somewhat fortuitous, however, and should certainly not be expected in general. First, a recent statistical assessment of dispersion-corrected density functionals for alkanes114 shows that B3LYP-D and B2-PLYP-D are indeed among the very best to reproduce high-level bond separation energies for alkanes, and that several other functionals perform worse. VSXC, not examined in that work, would likely show severe problems for branched alkanes, because it generally predicts them to be far too stable relative to unbranched ones (section 4.3.3 and Table 3), and dispersion terms would only aggravate that trend. Second, our statistical analysis for the AE and HOF test sets has indicated that gains in accuracy through BSR-corrected dispersion terms are much more moderate overall than for the specific case of alkanes (Tables 1 and 2). Finally, we note that some of the Minnesota functionals, and M05-2X in particular, are known to account for medium-range correlation effects through parametrization.115,116 Our analysis indeed

confirms that BSR-corrected M05-2X performs exceptionally well for linear and branched alkanes (Table 3), showing no need for the inclusion of dispersion terms to treat these molecules.

5. DISCUSSION AND CONCLUSIONS BSR corrections as implemented in the ATOMIC model consistently reduce the large errors that LDA and GGA functionals show for atomization energies and derived heats of formation. Another consistent improvement concerns the basis set incompleteness error. Although basis set dependence is generally lower for DFT than for correlated wave function theory, BSR corrections help to reduce it further. This is particularly noticeable for small basis sets, where excessive error is reduced significantly but becomes important also when one adds empirical dispersion terms to a density functional. Dispersion terms are meant to improve on the description of nonlocal correlation effects which favor compact over extended structures and large over small molecules. However, they always increase atomization energies, and their discriminatory power may be obscured by the use of a density functional that tends to overbind at the basis set limit as we have observed for double hybrid functionals. BSR corrections cancel such basis set effects to a large extent and thus restore the benefit of adding dispersion terms. The analysis of third- and higher-rung functionals, generally more reliable for thermochemistry, showed two unexpected observations: First, not all functionals benefit from the application of BSR corrections, and at least one (VSXC) performs significantly worse. Second, some of the most highly parametrized hybrid functionals (BMK, M05-2X, M06-2X), already performing fairly well for atomization reactions, show the most significant improvement through BSR corrections. After correction, they are generally more accurate than even double-hybrid functionals (B2-PLYP, B2-PLYP-D, B2GPPLYP) or correlated wave function theory at the MP2 and CCSD levels. Still, CCSD(T) outperforms them by a significant margin, and the inclusion of connected triples, even evaluated only with small basis sets, was found to be critical in designing ATOMIC composite models of acceptable accuracy. The midrange model ATOMIC/B3, previously suggested43 as the best compromise between accuracy and computational expense, evaluates connected triples corrections from double-ζ basis sets and it carries a 6−8 times lower rms error than the best BSRcorrected hybrid functionals for the AE benchmark, and a 2− 2.5 times lower rms error for the HOF benchmark. It is obvious that none of the considered density functionals can replace the composite wave function model implemented in the ATOMIC protocol. However, some of the most successful functionals may serve as low-cost alternatives if demands are lower. Balancing cost and accuracy, M05-2X/6-311+G(3df,2p) appears to be the overall best recommendation for such purposes. The suggested basis set is about as large as cc-pVTZ, but it produces BSR-corrected results of cc-pVQZ-quality for this functional (Tables 1 and 2: basis set A). When this suggestion is followed, it may be advantageous to perform BMK/6-311+G(3df,2p) calculations as well. The two functionals perform almost equally well overall but show little correlation between individual errors,117 so that the comparison of results may give some hint on expected accuracies. The alternative M06-2X is less recommended as it shows larger maximum errors than either M05-2X and BMK. (ΔMax − ΔMin in kcal/mol for AE/HOF sets, from Table S1 (Supporting 240

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

functionals show problems, but they are not a panacea for all types of problems observed. One particular example is the poor performance of M05 for highly strained cyclic systems where hyperhomodesmotic reactions are unlikely to show a significant advantage, as all elemental products and reactants are acyclic. And finally, it seems to be a major challenge to generalize the concept of hyperhomodesmotic reactions to molecules other than hydrocarbons. The very recent development of a generalized connectivitybased hierarchy model122,123 avoids the complexity of homodesmotic52,124 and related reactions, and strictly bases formulation of reactions on atom connectivity. It certainly holds promise for relief with respect to some of the challenges mentioned above, and allows, in particular, for a seamless generalization beyond hydrocarbons. Of course, also this approach cannot help with density functional errors encountered for topologies (e.g., small rings) not preserved in the reactions, and, if implemented with high-level ab initio reference data, it still faces the challenge of generating these data, because reactions involve significantly more and larger parent molecules than bond separation reactions do. The concept of bond separation reactions is very simple to implement and it has thus been used for ATOMIC. Application to a particular molecule only requires that there exists a unique valence structure. Slightly more than 30 bond types and about as many BSR parent molecules representing them have been found sufficient to cover most closed-shell neutral molecules containing H, C, N, O, and F, excluding a few special cases such as the elements (e.g., H2).42 This simplicity comes along with high demands placed on both the accuracy of the reference level of theory and the systematic nature of residual errors for the production level. Composite wave function approaches meet the latter requirement to a very high degree, density functionals to varying degrees. The most attractive option among density functionals is M05-2X, not only because BSR corrections lead to significant improvement but also because these ab initio derived bond corrections are close to optimal across the entire range of bonds considered. This robustness is a significant advantage as it obviates the need for calibration.

Information): 12.1/11.2 (BMK), 10.6/13.4 (M05-2X), 19.9/ 17.7 (M06-2X)). One particularly pleasing characteristic of M05-2X (and of BMK, not shown) is the observation that BSR corrections prove to be reliable not only for the allegedly, though in practice rarely, simple case of hydrocarbons but also for more complex valence structures such as those of compounds with NO bonds. Empirically calibrated bond “additivity” corrections (“BAC”) generally help with density functionals for which the application of bond separation reactions shows too little improvement or even deterioration of results. A tentative analysis of BACs has indeed attributed the poor performance of BSR-corrected VSXC to an incompatibility in error between larger molecules and the corresponding small BSR prototypes representing their valence structure. Still, a number of problems relating to nonlocal features (branching, ring strain, etc.) have been identified for several functionals, which are all resistant to any kind of bond-type based correction. M05 is the example of a functional that for this reason benefits fairly little from either proper BSR or empirical bond corrections. On the other hand, all double hybrid functionals attain remarkable accuracy with empirical bond corrections and appear to be the preferred choice in a BAC context. Even though they are still not fully competitive with composite wave function approaches, they produce results of very useful accuracy and are of course much less demanding computationally, at least if used with basis sets smaller than employed for this analysis (cc-pV5Z). Additional tests have shown that such basis set reduction has little effect on statistical results (rms errors for B2-PLYP/B2-PLYP-D/B2GPPLYP with basis set A; AE 1.6, 1.6, 1.4; HOF 1.2, 1.1, 1.2 kcal/ mol; compare to Tables 1 and 2). The calibration of BAC corrections must be considered tentative. Serious parametrization efforts would require a significantly larger reference data set that represents all types of bonds to similar extent and in as many combinations as possible. Such efforts are beyond the scope of this work, which uses calibration chiefly to shed some light on the puzzling failure of several density functionals for bond separation reactions. The concept of bond separation reactions, however, is central to the strictly ab initio design of the ATOMIC protocol, and any method showing poor performance in this regard must be considered inappropriate for use with this protocol. Better results on average may be achieved with isodesmic reaction schemes that keep larger molecular fragments untouched than just bonds between atoms. One particular choice is the set of elemental hyperhomodesmotic reactions suggested for hydrocarbons by Wheeler et al.118 Their use faces a number of challenges, however. First, definition of a unique set is required for the setup of an unambiguous and automatable procedure, and this seems to be far from trivial as evidenced by two revisions following Wheeler’s original paper.119−121 Second, the number of elemental reactants and products grows substantially, as bond-separation reactions require only four for hydrocarbons (methane, ethane, ethylene, acetylene), whereas hyperhomodesmotic reactions require 35.121 Third, these elemental reactants and products are generally larger in size than for bond separation reactions, posing significant problems to the use of very high level ab initio procedures to produce reference data, in particular for electron-correlation corrections beyond CCSD(T). Fourth, hyperhomodesmotic reactions certainly help to deal with phenomena such as protobranching, for which many density

■

ASSOCIATED CONTENT

S Supporting Information *

Table S1 lists the complete set of statistical data (rms errors, mean signed errors with standard deviations, mean unsigned errors, maximum positive and negative errors) for the AE and HOF sets. This information is available free of charge via the Internet at http://pubs.acs.org

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The author declares no competing financial interest.

■

ACKNOWLEDGMENTS The author thanks the Swiss National Supercomputing Center (CSCS, Manno TI, now Lugano TI) for generous grants of computer time.

■

REFERENCES

(1) Helgaker, T.; Klopper, W.; Tew, D. P. Mol. Phys. 2008, 106, 2107−2143. 241

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

(2) Peterson, K. A.; Feller, D.; Dixon, D. A. Theor. Chem. Acc. 2012, 131, 1079. (3) Martin, J. M. L.; de Oliveira, G. J. Chem. Phys. 1999, 111, 1843− 1856. (4) Boese, A. D.; Oren, M.; Atasoylu, O.; Martin, J. M. L.; Kállay, M.; Gauss, J. J. Chem. Phys. 2004, 120, 4129−4141. (5) Karton, A.; Rabinovich, E.; Martin, J. M. L.; Ruscic, B. J. Chem. Phys. 2006, 125, 144108. (6) Karton, A.; Martin, J. M. L. J. Chem. Phys. 2012, 136, 124114. (7) Tajti, A.; Szalay, P. G.; Császár, A. G.; Kállay, M.; Gauss, J.; Valeev, E. F.; Flowers, B. A.; Vázquez, J.; Stanton, J. F. J. Chem. Phys. 2004, 121, 11599−11613. (8) Bomble, Y. J.; Vázquez, J.; Kállay, M.; Michauk, C.; Szalay, P. G.; Császár, A. G.; Gauss, J.; Stanton, J. F. J. Chem. Phys. 2006, 125, 064108. (9) Harding, M. E.; Vázquez, J.; Ruscic, B.; Wilson, A. K.; Gauss, J.; Stanton, J. F. J. Chem. Phys. 2008, 128, 114111. (10) Allen, W. D.; East, A. L. L.; Császár, A. G. In Structures and conformations of non-rigid molecules; Laane, J., Dakkouri, M., van der Veken, B., Oberhammer, H., Eds.; Kluwer: Dordrecht, The Netherlands, 1993; pp 343−373. (11) East, A. L. L.; Allen, W. D. J. Chem. Phys. 1993, 99, 4638−4650. (12) Klippenstein, S. J.; East, A. L. L.; Allen, W. D. J. Chem. Phys. 1996, 105, 118−140. (13) Császár, A. G.; Allen, W. D.; Schaefer, H. F., III. J. Chem. Phys. 1998, 108, 9751−9764. (14) Császár, A. G.; Szalay, P. G.; Leininger, M. L. Mol. Phys. 2002, 100, 3879−3883. (15) Császár, A. G.; Leininger, M. L.; Burcat, A. J. Phys. Chem. A 2003, 107, 2061−2065. (16) Császár, A. G.; Leininger, M. L.; Szalay, V. J. Chem. Phys. 2003, 118, 10631−10642. (17) Bak, K. L.; Jørgensen, P.; Olsen, J.; Helgaker, T.; Klopper, W. J. Chem. Phys. 2000, 112, 9229−9242. (18) Helgaker, T.; Klopper, W.; Halkier, A.; Bak, K. L.; Jorgensen, P.; Olsen, J. In Quantum-mechanical prediction of thermochemical data; Cioslowski, J., Ed., Kluwer: Dordrecht, The Netherlands, 2001; pp 1− 30. (19) Helgaker, T.; Ruden, T. A.; Jørgensen, P.; Olsen, J.; Klopper, W. J. Phys. Org. Chem. 2004, 17, 913−933. (20) Klopper, W.; Ruscic, B.; Tew, D. P.; Bischoff, F. A.; Wolfsegger, S. Chem. Phys. 2009, 356, 14−24. (21) Klopper, W.; Bachorz, R. A.; Haettig, C.; Tew, D. P. Theor. Chem. Acc. 2010, 126, 289−304. (22) Haunschild, R.; Klopper, W. Theor. Chem. Acc. 2012, 131, 1112. (23) Feller, D.; Dixon, D. A.; Peterson, K. A. J. Phys. Chem. A 1998, 102, 7053−7059. (24) Dixon, D. A.; Feller, D.; Sandrone, G. J. Phys. Chem. A 1999, 103, 4744−4751. (25) Feller, D.; Dixon, D. A. J. Phys. Chem. A 1999, 103, 6413−6419. (26) Feller, D.; Dixon, D. A. J. Chem. Phys. 2001, 115, 3484−3496. (27) Ruscic, B.; Pinzon, R. E.; Morton, M. L.; von Laszevski, G.; Bittner, S. J.; Nijsure, S. G.; Amin, K. A.; Minkoff, M.; Wagner, A. F. J. Phys. Chem. A 2004, 108, 9979−9997. (28) Császár, A. G.; Furtenbacher, T. Chem.Eur. J. 2010, 16, 4826−4835. (29) Purvis, G. D.; Bartlett, R. J. J. Chem. Phys. 1982, 76, 1910−1918. (30) Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. Chem. Phys. Lett. 1989, 157, 479−483. (31) Pople, J. A.; Head-Gordon, M.; Fox, D. J.; Raghavachari, K.; Curtiss, L. A. J. Chem. Phys. 1989, 90, 5622−5629. (32) Curtiss, L. A.; Raghavachari, K.; Trucks, G. W.; Pople, J. A. J. Chem. Phys. 1991, 94, 7221−7230. (33) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Rassolov, V.; Pople, J. A. J. Chem. Phys. 1998, 109, 7764−7776. (34) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. J. Chem. Phys. 2007, 126, 084108. (35) Petersson, G. A.; Tensfeldt, T. G.; Montgomery, J. A., Jr. J. Chem. Phys. 1991, 94, 6091−6101.

(36) Montgomery, J. A., Jr.; Ochterski, J. W.; Petersson, G. A. J. Chem. Phys. 1994, 101, 5900−5909. (37) Ochterski, J. W.; Petersson, G. A.; Montgomery, J. A., Jr. J. Chem. Phys. 1996, 104, 2598−2619. (38) Petersson, G. A. In Quantum-mechanical prediction of thermochemical data; Cioslowski, J., Ed.; Kluwer: Dordrecht, The Netherlands, 2001; pp 99−130. (39) DeYonker, N. J.; Cundari, T. R.; Wilson, A. K. J. Chem. Phys. 2006, 124, 114104. (40) DeYonker, N. J.; Grimes, T.; Yockel, S.; Dinescu, A.; Mintz, B.; Cundari, T. R.; Wilson, A. K. J. Chem. Phys. 2006, 125, 104111. (41) DeYonker, N. J.; Wilson, B. R.; Pierpont, A. W.; Cundari, T. R.; Wilson, A. K. Mol. Phys. 2009, 107, 1107−1121. (42) Bakowies, D. J. Chem. Phys. 2009, 130, 144113. (43) Bakowies, D. J. Phys. Chem. A 2009, 113, 11517−11534. (44) Ditchfield, R.; Hehre, W. J.; Pople, J. A.; Radom, L. Chem. Phys. Lett. 1970, 5, 13−14. (45) Hehre, W. J.; Ditchfield, R.; Radom, L.; Pople, J. A. J. Am. Chem. Soc. 1970, 92, 4796−4801. (46) Raghavachari, K.; Stefanov, B. B.; Curtiss, L. A. Mol. Phys. 1997, 91, 555−559. (47) Redfern, P. C.; Zapol, P.; Curtiss, L. A.; Raghavachari, K. J. Phys. Chem. A 2000, 104, 5850−5854. (48) Sivaramakrishnan, R.; Tranter, R. S.; Brezinsky, K. J. Phys. Chem. A 2005, 109, 1621−1628. (49) Wodrich, M. D.; Corminboeuf, C.; Schreiner, P. R.; Fokin, A. A.; Schleyer, P. v. R. Org. Lett. 2007, 9, 1851−1854. (50) Raghavachari, K.; Stefanov, B. B.; Curtiss, L. A. J. Chem. Phys. 1997, 106, 6764−6767. (51) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Stefanov, B. B. J. Chem. Phys. 1998, 108, 692−697. (52) Wheeler, S. E. WIREs Comput. Mol. Sci 2012, 2, 204−220. (53) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623−11627. (54) Perdew, J. P.; Schmidt, K. In Density functional theory and its application to materials; van Doren, V., van Alsenoy, C., Geerlings, P., Eds.; AIP Conference Proceedings, Vol. 577; 2001; pp 1−20. (55) Slater, J. C. The self-consistent field for molecules and solids: Quantum theory of molecules and solids; McGraw Hill: New York, 1974; Vol. 4. (56) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200− 1211. (57) Perdew, J. P. Phys. Rev. B 1986, 33, 8822−8824. (58) Becke, A. D. Phys. Rev. A 1988, 38, 3098−3100. (59) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785−789. (60) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (61) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158−6170. (62) van Voorhis, T.; Scuseria, G. E. J. Chem. Phys. 1998, 109, 400− 410. (63) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Phys. Rev. Lett. 2003, 91, 146401. (64) Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. J. Chem. Phys. 2003, 119, 12129−12137. (65) Csonka, G. I.; Ruzsinszky, A.; Tao, J.; Perdew, J. P. Int. J. Quantum Chem. 2005, 101, 506−511. (66) Zhao, Y.; Truhlar, D. G. J. Chem. Phys. 2006, 125, 194101. (67) Becke, A. D. J. Chem. Phys. 1993, 98, 5648−5652. (68) Boese, A. D.; Martin, J. M. L. J. Chem. Phys. 2004, 121, 3405− 3416. (69) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. J. Chem. Phys. 2005, 123, 161103. (70) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215−241. (71) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. J. Chem. Theory Comput. 2006, 2, 364−382. (72) Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2006, 110, 13126− 13130. (73) Grimme, S. J. Chem. Phys. 2006, 124, 034108. 242

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243

The Journal of Physical Chemistry A

Article

(74) Karton, A.; Tarnopolsky, A.; Lamère, J.-F.; Schatz, G. C.; Martin, J. M. L. J. Phys. Chem. A 2008, 112, 12868−12886. (75) Schwabe, T.; Grimme, S. Phys. Chem. Chem. Phys. 2007, 9, 3397−3406. (76) A complete listing of all bond types may be found in Table 2 of ref 42, which also includes data for three bond types not present in any of the reference data sets considered here (C-H and C-O with divalent carbon, O-H(+) with trivalent oxygen). (77) Møller, C.; Plesset, M. S. Phys. Rev. 1934, 46, 618−622. (78) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007−1023. (79) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Pople, J. A. J. Chem. Phys. 2000, 112, 7374−7383. (80) Graham, D. C.; Menon, A. S.; Goerigk, L.; Grimme, S.; Radom, L. J. Phys. Chem. A 2009, 113, 9861−9873. (81) Some of the basis sets were obtained from the EMSL basis set library.82,83 (82) Feller, D. J. Comput. Chem. 1996, 17, 1571−1586. (83) Schuchardt, K. L.; Didier, B. T.; Elsethagen, T.; Sun, L.; Gurumoorthi, V.; Chase, J.; Li, J.; Windus, T. L. J. Chem. Inf. Model. 2007, 47, 1045−1052. (84) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650−654. (85) Frisch, M. J.; Pople, J. A.; Binkley, J. S. J. Chem. Phys. 1984, 80, 3265−3269. (86) Rabuck, A. D.; Scuseria, G. E. Chem. Phys. Lett. 1999, 309, 450− 456. (87) Jensen, F. J. Chem. Phys. 2001, 115, 9113−9125. (88) Jensen, F. J. Chem. Phys. 2002, 116, 7372−7379. (89) Jensen, F. J. Chem. Phys. 2002, 117, 9234−9240. (90) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. (91) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Chem. Phys. Lett. 1989, 162, 165−169. (92) Treutler, O.; Ahlrichs, R. J. Chem. Phys. 1995, 102, 346−354. (93) Hertwig, R. H.; Koch, W. Chem. Phys. Lett. 1997, 268, 345−351. (94) Furche, F.; Perdew, J. P. J. Chem. Phys. 2006, 124, 044103. (95) Weigend, F.; Häser, M. Theor. Chem. Acc. 1997, 97, 331−340. (96) Zhao, Y.; Pu, J. Z.; Lynch, B. J.; Truhlar, D. G. Phys. Chem. Chem. Phys. 2004, 6, 673−676. (97) Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Pople, J. A. J. Chem. Phys. 1997, 106, 1063−1079. (98) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. J. Chem. Phys. 2005, 123, 124107. (99) Schreiner, P. R.; Fokin, A. A.; Pascal, Robert, A., J.; de Meijere, A. Org. Lett. 2006, 8, 3635−3638. (100) Korth, M.; Grimme, S. J. Chem. Theory Comput. 2009, 5, 993− 1003. (101) Wang, N. X.; Wilson, A. K. J. Chem. Phys. 2004, 121, 7632− 7646. (102) Bakowies, D. J. Chem. Phys. 2007, 127, 084105/1−23. (103) Berry, R. J.; Burgess, D. R. F., Jr.; Nyden, M. R.; Zachariah, M. R.; Melius, C. F.; Schwartz, M. J. Phys. Chem. 1996, 100, 7405−7410. (104) Petersson, G. A.; Malick, D. K.; Wilson, W. G.; Ochterski, J. W.; Montgomery, J. A.; Frisch, M. J. J. Chem. Phys. 1998, 109, 10570− 10579. (105) Besler, B. H.; Merz, K. M.; Kollman, P. A. J. Comput. Chem. 1990, 11, 431−439. (106) The designation of bond types includes the number of free valencies on either atom as superscript. An appended “(+)” indicates that the formal total charge on the two atoms equals +1. See also ref 42. (107) Grimme, S. J. Comput. Chem. 2006, 27, 1787−1799. (108) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. J. Chem. Phys. 2010, 132, 154104. (109) Wodrich, M. D.; Corminboeuf, C.; Schleyer, P. v. R. Org. Lett. 2006, 8, 3631−3634. (110) Grimme, S. Angew. Chem., Int. Ed. 2006, 45, 4460−4464.

(111) Wodrich, M. D.; Wannere, C. S.; Mo, Y.; Jarowski, P. D.; Houk, K. N.; Schleyer, P. v. R. Chem.Eur. J. 2007, 13, 7731−7744. (112) Brittain, D. R. B.; Lin, C. Y.; Gilbert, A. T. B.; Izgorodina, E. I.; Gill, P. M. W.; Coote, M. L. Phys. Chem. Chem. Phys. 2009, 11, 1138− 1142. (113) Song, J.-W.; Tsuneda, T.; Sato, T.; Hirao, K. Org. Lett. 2010, 12, 1440−1443. (114) Karton, A.; Gruzman, D.; Martin, J. M. L. J. Phys. Chem. A 2009, 113, 8434−8447. (115) Zhao, Y.; Truhlar, D. G. Org. Lett. 2006, 8, 5753−5755. (116) Grimme, S. WIREs Comput. Mol. Sci 2011, 1, 211−228. (117) A linear regression analysis (y = b·x + a) between (y) BMK/6311+G(3df,2p) errors after BSR correction and (x) M05-2X/6311+G(3df,2p) errors after BSR correction yields a = 0.13, b = 0.80, with a correlation coefficient of R2 = 0.42 and a = 0.47, b = 0.48, R2 = 0.28 for the HOF and AE benchmarks, respectively. (118) Wheeler, S. E.; Houk, K. N.; Schleyer, P. V. R.; Allen, W. D. J. Am. Chem. Soc. 2009, 131, 2547−2560. (119) Wodrich, M. D.; Corminboeuf, C.; Wheeler, S. E. J. Phys. Chem. A 2012, 116, 3436−3447. (120) Fishtik, I. J. Phys. Chem. A 2012, 116, 8792−8793. (121) Wodrich, M. D.; Gonthier, J. F.; Corminboeuf, C.; Wheeler, S. E. J. Phys. Chem. A 2012, 116, 8794−8796. (122) Ramabhadran, R. O.; Raghavachari, K. J. Chem. Theory Comput. 2011, 7, 2094−2103. (123) Ramabhadran, R. O.; Raghavachari, K. J. Phys. Chem. A 2012, 116, 7531−7537. (124) George, P.; Trachtman, M.; Bock, C. W.; Brett, A. M. Theor. Chim. Acta 1975, 38, 121−129.

243

dx.doi.org/10.1021/jp310735h | J. Phys. Chem. A 2013, 117, 228−243