G3X-K theory: A composite theoretical method for thermochemical kinetics

Chemical Physics Letters 558 (2013) 109–113

Contents lists available at SciVerse ScienceDirect

Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett

G3X-K theory: A composite theoretical method for thermochemical kinetics Gabriel da Silva Department of Chemical and Biomolecular Engineering, The University of Melbourne, Victoria 3010, Australia

a r t i c l e

i n f o

Article history: Received 2 October 2012 In final form 21 December 2012 Available online 2 January 2013

a b s t r a c t A composite theoretical method for accurate thermochemical kinetics, G3X-K, is described. This method is accurate to around 0.5 kcal mol1 for barrier heights and 0.8 kcal mol1 for enthalpies of formation. G3X-K is a modification of G3SX theory using the M06-2X density functional for structures and zeropoint energies and parameterized for a test set of 223 heats of formation and 23 barrier heights. A reduced perturbation-order variant, G3X(MP3)-K, is also developed, providing around 0.7 kcal mol1 accuracy for barrier heights and 0.9 kcal mol1 accuracy for enthalpies, at reduced computational cost. Some opportunities to further improve Gn composite methods are identified and briefly discussed. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Composite electronic structure theory methods have gained widespread use for the efficient calculation of chemical reaction energies. At one end, methods from the Gn [1–4] and CBS [5,6] classes can provide thermochemical properties with mean unsigned errors (MUEs) at around ‘chemical accuracy’ of 1 kcal mol1 (4 kJ mol1) for even moderately large compounds (tens of heavy atoms). At the other end, methods such as those in the Wn [7,8] and HEAT [9,10] classes, based on extrapolating coupled cluster energies to the complete basis set limit, can provide ‘benchmark accuracy’ of 1 kJ mol1 (0.25 kcal mol1) but at such an increased cost that they can only be routinely applied to molecules with less than say 6 heavy atoms. The ‘chemically accurate’ methods in particular have been widely used to calculate barrier heights for theoretical kinetic modeling, although relatively little attention has been paid to developing [11] (or even benchmarking [12]) composite theoretical methods for thermochemical kinetics. In contrast, the last decade has seen a number of exchange–correlation functionals developed with kinetics applications in mind [13–15]. Zheng et al. [12] have recently developed a representative set of barrier heights, DBH24/08, and used it to benchmark a wide selection of wavefunction theory, density functional theory (DFT) and composite (or multilevel) methods. Of the composite methods examined, G3SX and G3SX(MP3) [16] were the best-performed, both achieving a MUE of 0.57 kcal mol1. The accuracy of these methods is comparable to CCSD(T) theory with a correlation consistent polarized valence triple-zeta basis set, calculations that are well over an order of magnitude more computationally expensive. The G3SX and G3SX(MP3) methods also offered a E-mail address: [email protected] 0009-2614/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cplett.2012.12.045

slight improvement over the more-expensive G4 theory (MUE of 0.58 kcal mol1) and significantly outperformed the popular CBS-QB3 method (MUE of 1.62 kcal mol1). It is thus apparent that G3SX theory offers an excellent balance between cost and accuracy for barrier heights of moderately sized systems, and for this reason it has been used extensively in our recent work [17–24]. The results of Zheng et al. come with a caveat, however; they are for optimized geometries obtained with QCISD theory and a triple-zeta quality Gaussian basis set (MG3). The original G3SX method [16] uses B3LYP/6-31G(2df,p) optimized structures, and it is well-known that the B3LYP functional performs poorly for the non-bonded interactions necessary to accurately describe the energy, and to a lesser extent structure [25], of transition states. Here we show that B3LYP structures seriously degrade the performance of G3SX theory for the DBH24/08 barrier heights. To rectify this we have implemented a modification of G3SX theory that uses the M06-2X density functional [26] (which does perform well for barrier heights) in place of B3LYP. The resultant method, referred to as G3X-K theory, provides a MUE of less than 0.5 kcal mol1 for the DBH24/08 barrier heights and around 0.8 kcal mol1 for G3/99 test set [27] atomization enthalpies of formation. A cost-effective reduced perturbation order variant, G3X(MP3)-K, is also described. Finally, the implications of this work for the development of future composite methods are discussed. 2. Method description The G3X-K and G3X(MP3)-K methods represent relatively minor modifications of the G3SX and G3SX(MP3) methods, which are described in detail in Ref. [16]. There are, indeed, three essential differences between the original and the ‘kinetics’ versions of these methods:

110

G. da Silva / Chemical Physics Letters 558 (2013) 109–113

Table 1 Scale factors in the G3X-K and G3X(MP3)-K composite methods.a

SHF SE234 SCC SE2 SE3 SE4

G3X-K

G3X(MP3)-K

1.0945 1.0712 1.2347 1.1467 1.2249 0.4135

1.1014 1.0642 1.2559 1.1643 1.0879 –

152.7 kcal mol1) are in error by several kilocalories per mole. Enthalpies of formation are calculated using the experimental 0 K enthalpies of formation and 298 K enthalpy corrections for atoms employed with the other Gn methods [32]. Note that the empirical scaling factors (vs. higher level corrections, HLCs) are employed in this Letter; because the HLC largely cancels for barrier height determinations (and completely cancels for isomerizations) relative to atomization reactions it is less suitable for kinetics applications.

a Scale factors are optimized to minimize the MUE in the 23 barrier heights and 223 enthalpies of formation in the test set.

3. Results and discussion (i) The use of M06-2X/6-31G(2df,p) optimized structures and scaled zero point energies. (ii) The use of coupled cluster theory energies in place of quadratic configuration interaction theory [28]. (iii) Determination of new scaling factors that minimize the MUE over the DBH24/08 barrier heights and G3/99 atomization enthalpies. In G3X-K theory M06-2X optimized structures are used in the series of calculations HF/G3XLarge, MP2(full)/G3Large, MP4/631G(2df,p), MP4/6-31+G(d), CCSD(T)/6-31G(d), with the final energy given by: E0(G3X-K) = HF/6-31G(d) + SHF[HF/G3XL  HF/d] + SE234 [E2/d + E3/d + E4/d] + SCC[ECC/d] + SE2[E2(FU)/G3L  E2/d] + SE3 [E3/+  E3/d] + SE3[E3/2dfp  E3/d] + SE4[E4/+  E4/d] + SE4 [E4/2dfp  E4/d] + E(SO) + E(ZPE) For the G3X(MP3)-K method the energy is determined as follows, where calculations with the 6-31+G(d) basis set are no longer required, and the 6-31G(2df,p) calculations are evaluated only up to the MP3 level: E0(G3X(MP3)-K) = HF/6-31G(d) + SHF[HF/G3XL  HF/ d] + SE234[E2/d + E3/d + E4/d] + SCC[ECC/d] + SE2[E2(FU)/G3L  E2/ d] + SE3[E3/2dfp  E3/d] + E(SO) + E(ZPE) In both methods the energy increment terms are defined as: E2/d = MP2/6-31G(d)  HF/6-31G(d) E3/d = MP3/6-31G(d)  MP2/6-31G(d) E4/d = MP4/6-31G(d)  MP3/6-31G(d) ECC/d = CCSD(T)/6-31G(d)  MP4/6-31G(d) E2(FU)/G3L = MP2(full)/G3Large  HF/G3large E3/+ = MP3/6-31+G(d) – MP2/6-31+G(d) E3/2dfp = MP3/6-31G(2df,p)  MP2/6-31G(2df,p) E4/+ = MP4/6-31+G(d) – MP3/6-31+G(d) E4/2dfp = MP4/6-31G(2df,p) – MP3/6-31G(2df,p). The zero point energy, E(ZPE), is scaled by 0.972 [26] and the spin–orbit correction, E(SO), is applied to atoms and diatomic radicals as per G4 theory [3]. Post-SCF calculations exclude core electrons (frozencore), except for the MP2(full)/G3Large job step in which all electrons are included in the correlation calculation (full). All energy calculations were performed using GAUSSIAN 09 [29], and minima optimized at the M06-2X level were confirmed to have zero imaginary vibrational frequencies. In both methods the S terms are scaling factors (Table 1), optimized to minimize the MUE for 246 energies in the test set. This test constitutes 223 enthalpies of formation from the G3/99 database and 23 barrier heights from the DBH24/08 database. The DBH24/08 barrier heights are from benchmark quality ab initio calculations. The experimental G3/99 enthalpies of formation are largely used as is, even though several of these enthalpies have been since refined by small amounts. The only changes are for vinyl chloride (5.4 kcal mol1 [30]) and CF2O (145.3 kcal mol1 [31]) where the values used in the G3/99 test set (8.4 and

Barrier heights in the DBH24/08 database have been calculated using the G3X-K and G3X(MP3)-K methods, with MUEs reported in Table 2. Included in this table are MUEs obtained using the original G3SX method as well as the results of Zheng et al. using QCISD/ MG3 optimized structures, which we denote G3SX//QCI. As discussed above, it is apparent from the work of Zheng et al. that G3SX theory is capable of providing an accurate description of barrier heights, with an MUE of only 0.57 kcal mol1 for G3SX//QCI. However, when B3LYP optimized structures are used, as in the original G3SX method, the MUE increases significantly to 1.11 kcal mol1. Note also that this excludes barrier heights for the CH3 + FCl M CH3F + Cl reaction, where B3LYP theory fails to locate the transition state altogether, as well as the F  CH3Cl ? Cl  CH3F barrier, where the reactant complex is unstable. Curtiss et al. [33] have demonstrated that Gn methods employing B3LYP/6-31G(2df,p) optimized structures are unsuitable for evaluating barrier heights for reactions involving fluorine; when these reactions are excluded from the database the G3SX MUE drops to 0.89 kcal mol1, which represents a significant improvement yet is still well short of the G3SX//QCI results. The above illustrates that B3LYP DFT structures present a serious limitation to accurate barrier height determinations within the Gn ansatz. To overcome this we have replaced the B3LYP functional with M06-2X (retaining the 6-31G(2df,p) basis set), a method that does perform well for barrier heights, in order to arrive at the G3X-K and G3X(MP3)-K methods. As an illustration of the accuracy of M06-2X vs. B3LYP, barrier heights in DBH24/08 are reproduced with a MUE of 0.98 kcal mol1 at the M06-2X/MG3S level of theory, but only 4.15 kcal mol1 at the B3LYP/MG3S level of theory (again using QCISD/MG3 structures). Note however that the weak F  CH3Cl complex still could not be located at the M06-2X/ 6-31G(2df,p) level, and this barrier is again excluded (reducing the total set to 23 barrier heights). Comparing the DFT structures with the high-level QCISD/MG3 ones, we find that for a set of 25 critical interatomic distances (see Supplementary Data) the M06-2X calculations achieve a mean unsigned deviation of 0.04 Å, compared to 0.06 Å for B3LYP. Moreover, the B3LYP calculations achieve a relatively large mean unsigned deviation of 0.04 Å (vs. 0.001 Å for M06-2X), indicating a systematic over-prediction of interatomic distances. This is consistent with the known poor treatment of dispersion forces in B3LYP, which is considerably improved in M062X. Although the conventional thinking may be that the B3LYP functional performs adequately for transition state structures (if Table 2 Mean unsigned errors (MUEs) for barrier heights in the DBH24/08 database using four variants of G3SX theory. MUE (kcal mol1) G3SX G3SX//QCIa G3X-K G3X(MP3)-K a

1.11 0.57 0.47 0.71

Denotes G3SX theory on QCISD/MG3 optimized structures.

G. da Silva / Chemical Physics Letters 558 (2013) 109–113

not energies), it is apparent from the present work that the errors in B3LYP structures are significant when chemical accuracy is desired. The results presented in Table 2 demonstrate that G3X-K theory offers a significant improvement over G3SX for barrier heights, yielding a MUE of 0.47 kcal mol1 for the DBH24/08 barrier heights. Note also that G3X-K theory surpasses the G3SX//QCI variant, which provides testimony to the general accuracy of the M062X optimized geometries. The reduced-order G3X(MP3)-K model chemistry also performs well, with MUE of 0.71 kcal mol1, although now this method does not outperform the G3SX(MP3)// QCI method (where MUE = 0.57 kcal mol1). Figure 1 provides a histogram of the barrier height errors for the two new methods as well as for G3SX theory. We can see that the G3X-K errors are grouped tightly, with all but five barrier heights falling within ±0.6 kcal mol1 of their benchmark values. With the G3X(MP3)-K variant now only about half of the barriers are within ±0.6 kcal mol1 although the errors are still centered around zero. For G3SX theory the error distribution is broader yet and there also appears to be a small systematic underprediction of the barrier heights, with a mean signed error (MSE) of 0.14 kcal mol1. Results for the G3X-K and G3X(MP3)-K methods against the complete test set of 23 barrier heights and 223 enthalpies of formation are summarized in Table 3. For all energies the MUE is 0.82 kcal mol1 with G3X-K. Note that with the original G3SX scale factors the MUE increases to 0.94 kcal mol1, illustrating the significance of re-optimizing these terms. The G3X-K method achieves a MUE of 0.85 kcal mol1 across the 223 enthalpies of formation, with a small mean signed error, in addition to the MUE of 0.47 kcal mol1 for the barrier heights discussed above. Not unexpectedly, the G3X-K method provides similar performance to G3SX

111

theory for enthalpies of formation, whereas it achieves considerably improved performance for barrier height determinations. This is consistent with B3LYP performing well for the predominantly covalently-bound compounds in the enthalpy test set. It is of relevance to note also that when geometry effects are excluded atomization enthalpies of formation present a considerably more challenging calculation than barrier heights, with the large change in energy upon atomization resulting in little cancelation of errors across the reaction. As such it is not surprising that G3X-K performs better for barrier heights than atomization enthalpies. Furthermore, for most chemical reaction energies (where many bonds are conserved) the barrier height uncertainty of around 0.5 kcal mol1 is probably more representative. Examining the different sub-categories in the test set it is clear that G3X-K does not perform particularly well for non-hydrogens, with a MUE of 1.55 kcal mol1. It is known that G4 theory significantly improves over previous Gn iterations for nonhydrogens, providing a better alternative for these systems, although at substantially increased computational cost (and without a scaling factor variant for kinetics applications). Note also that although G3X-K theory performs reasonably well for the radical subset it does demonstrate a systematic over-estimation, with MSE of +0.57 kcal mol1. This could likely be improved by using restricted-open-shell wavefunctions, as in the Gn-RAD methods [34]. For barrier heights, performance is relatively even across the various sub-categories, with the worst-performing category (heavy atom transfer, MUE = 0.70 kcal mol1) biased by a single outlier. Across the entire test set the largest outlier is the P4 enthalpy of formation, predicted to be too-large by 10.56 kcal mol1. This molecule is known to present a challenging case for electronic structure theory methods [35,36]. Several other problem cases in the

Figure 1. Histograms of the DBH24/08 barrier height errors with the G3X-K, G3X(MP3)-K, and G3SX methods.

112

G. da Silva / Chemical Physics Letters 558 (2013) 109–113

Table 3 Mean unsigned errors (MUEs), mean signed errors (MSEs) and maximum unsigned errors (MaxUEs) for DBH24/08 barrier heights and G3/99 enthalpies of formation with the G3X-K and G3X(MP3)-K methods. G3X-K

G3X(MP3)-K

MUE

MSE

MaxUE

Barrier heights (23) Heavy-atom transfer (6) Nucleophilic substitution (5) Unimolecular and association (6) Hydrogen transfer (6)

0.47 0.70 0.27 0.41 0.48

MUE

MSE

0.05 0.29 0.10 0.29 0.46

2.16 2.16 0.41 1.06 1.31

0.71 0.62 1.10 0.61 0.59

0.02 0.22 0.05 0.33 0.59

Enthalpies of formation (223) Hydrocarbons (38) Substituted hydrocarbons (91) Radicals (31) Inorganic hydrides (15) Non-hydrogens (48) All (246) All except non-hydrogens (198)

2.23 1.60 2.23 1.28 1.21

0.85 0.62 0.61 0.78 0.57 1.64

0.11 0.04 0.07 0.57 0.19 0.08

10.56 2.61 2.92 2.18 2.42 10.56

0.89 0.72 0.59 0.90 0.62 1.68

0.12 0.15 0.13 0.54 0.01 0.10

11.72 2.38 2.98 2.89 2.33 11.72

0.82 0.62

0.10 0.14

10.56 2.92

0.87 0.68

0.11 0.12

11.72 2.98

non-hydrogens include C2F4 (5.07 kcal mol1) and SF6 (3.96 kcal mol1). Note though that the experimental enthalpy of formation for C2F4 has been brought into question, possibly being too-large by up to 2.5 kcal mol1 [37]. Outside of the nonhydrogens the largest enthalpy of formation outlier is pyrazine (+2.92 kcal mol1), in the substituted hydrocarbons subset. Since the submission of this Letter, however, a revised experimental enthalpy for pyrazine has been published (48.6 kcal mol1) [38], in better agreement with the G3X-K value (49.8 kcal mol1). For barrier heights the most significant outlier is for OH + N2 ? N2O + H (error of 2.16 kcal mol1), which is not surprising considering that this also corresponds to the largest barrier in the database (82.47 kcal mol1). Much of the above analysis for G3X-K theory applies to G3X(MP3)-K as well, although with respective MUEs in barrier heights and enthalpies of formation typically 0.1 and 0.2 kcal mol1 larger, respectively. The MUE for enthalpies of formation is 0.88, or 0.68 kcal mol1 with non-hydrogens excluded, whereas the value for barrier heights is 0.68 kcal mol1. The largest outlier, overall, is again P4. Within the barrier heights nucleophilic substitution of Cl for F now provides the largest error, followed by the OH + N2 reaction (1.82 kcal mol1). Both forward and reverse barriers in the CH3F + OH M CH3OH + F reaction have errors larger than 1 kcal mol1 as well, whereas they were accurately characterized by G3X-K theory. The present work provides some insight into means of improving Gn-like composite theoretical methods, other than to just avoid B3LYP optimized structures. Examining the scaling factors presented in Table 1 we see that the smallest term is SE4 (0.4135), which scales the contributions from MP4 theory (i.e., the [MP4  MP3] components). The relatively small contribution of these [MP4  MP3] energy increments can explain the excellent performance of the MP3-based reduced-order variant G3X(MP3)K, which is similar to G3X-K with SE4 = 0. More importantly, however, this reveals that relative to the unscaled energy increments the most computationally expensive step in G3X-K theory, the MP4/6-31G(2df,p) energy, actually provides the smallest contribution to the total energy! Much larger contributions relative to their unscaled energies arise from the scaled MP2, MP3 and QCISD(T) energy terms, which are all less computationally expensive to determine. Future methods may be able to exploit this by increasing the number of basis functions used in these job steps and perhaps even removing the MP4/6-31G(2df,p) calculation altogether; indeed, the correlation consistent Composite Approach (ccCA) [39,40] features a complete basis set extrapolation on HF and MP2 energies and can achieve 1 kcal mol1 accuracy without an MP4 energy or the need to incorporate empirical corrections.

MaxUE

Within G3X-K theory the small contribution of the MP4 energy and importance of the MP3 energy can perhaps be attributed to convergence issues with perturbation theory as well as to the nature of the MP4 single (S), double (D), triple (T) and quadruple (Q) excitations. We know that MPn energies do not necessarily converge to the full configuration interaction limit as n ? 1 [41–43] and for any given molecule it is possible that perturbation theory is better accounting for electron correlation at lower orders. Within MP4 theory the S, D, and Q terms only describe correlation effects up to pairs of electrons, as do double excitations in MP2 and MP3, whereas the T excitation term introduces (and can overestimate) three-electron correlation effects [42,43]. The positive Q energy term in MP4 is responsible for reducing electron pair correlation treatment that is over-estimated by MP2 theory although it cancels to some extent with the negative (stabilizing) triples excitation energy that first introduces (and again, probably overestimates) three-electron correlation. Thus, with Gn methods, scaled-up E3 = [MP3  MP2] energies may be able to provide a better description of electron pair correlation than E4 = [MP4(SDTQ)  MP3] energy terms, which are left to provide the only account of three electron correlation outside of coupled cluster theory, and thus need to be scaled-down heavily. This highlights why it may be more efficient for composite methods to treat electron correlation using a lower level of Møller–Plesset perturbation theory than MP4 (i.e., MP2 and/or MP3) but with a more complete basis set. Acknowledgments This work was supported by the Australian Research Council through Discovery Project DP110103889. Computational resources provided in part by the Victorian Partnership for Advanced Computing. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.cplett.2012.12.045. References [1] L.A. Curtiss, K. Raghavachari, G.W. Trucks, J.A. Pople, J. Chem. Phys. 94 (1991) 7221. [2] L.A. Curtiss, P.C. Redfern, K. Raghavachari, V. Rassolov, J.A. Pople, J. Chem. Phys. 110 (1999) 4703. [3] L.A. Curtiss, P.C. Redfern, K. Raghavachari, J. Chem. Phys. 126 (2007) 084108. [4] B. Chan, M.L. Coote, L. Radom, J. Chem. Theory Comput. 6 (2010) 2647. [5] J.W. Ochterski, G.A. Petersson, J.A. Montgomery, J. Chem. Phys. 104 (1996) 2598.

G. da Silva / Chemical Physics Letters 558 (2013) 109–113 [6] J.A. Montgomery, M.J. Frisch, J.W. Ochterski, G.A. Petersson, J. Chem. Phys. 110 (1999) 2822. [7] J.M.L. Martin, G. de Oliveira, J. Chem. Phys. 111 (1999) 1843. [8] A.D. Boese, M. Oren, O. Atasoylu, J.M.L. Martin, M. Kallay, J. Gauss, J. Chem. Phys. 120 (2004) 4129. [9] A. Tajti, P.G. Szalay, A.G. Csaszar, M. Kallay, J. Gauss, J. Chem. Phys. 121 (2004) 11599. [10] Y.J. Bomble, J. Vazquez, M. Kallay, C. Michauk, J. Chem. Phys. 125 (2006) 064108. [11] B. Chan, J. Deng, L. Radom, J. Chem. Theory Comput. 7 (2011) 112. [12] J. Zheng, Y. Zhao, D.G. Truhlar, J. Chem. Theory Comput. 5 (2009) 808. [13] A.D. Boese, J.M.L. Martin, J. Chem. Phys. 121 (2004) 3405. [14] Y. Zhao, B.J. Lynch, D.G. Truhlar, J. Phys. Chem. A 108 (2004) 2715. [15] Y. Zhao, N.E. Schultz, D.G. Truhlar, J. Chem. Theory Comput. 2 (2006) 364. [16] L.A. Curtiss, P.C. Redfern, K. Raghavachari, J.A. Pople, J. Chem. Phys. 114 (2001) 108. [17] G. da Silva, C. Graham, Z.-F. Wang, Environ. Sci. Technol. 44 (2010) 250. [18] G. da Silva, J.A. Cole, J.W. Bozzelli, J. Phys. Chem. A 114 (2010) 2275. [19] G. da Silva, Phys. Chem. Chem. Phys. 12 (2010) 6698. [20] G. da Silva, J. Phys. Chem. A 114 (2010) 6861. [21] G. da Silva, Angew. Chem. Int. Ed. 49 (2010) 7523. [22] G. da Silva, J. Phys. Chem. A 115 (2011) 291. [23] G. da Silva, A.J. Trevitt, Phys. Chem. Chem. Phys. 13 (2011) 8940. [24] G. da Silva, J. Phys. Chem. A 116 (2012) 5317. [25] L. Simon, J.M. Goodman, Org. Biomol. Chem. 9 (2011) 689. [26] Y. Zhao, D.G. Truhlar, Theor. Chem. Acc. 120 (2008) 215.

113

[27] L.A. Curtiss, K. Raghavachari, P.C. Redfern, J.A. Pople, J. Chem. Phys. 112 (2000) 7374. [28] L.A. Curtiss, P.C. Redfern, K. Raghavachari, J.A. Pople, J. Chem. Phys. 359 (2002) 390. [29] M.J. Frisch et al., Gaussian 09, revision B.01, Gaussian Inc., Wallingford, CT, 2010. [30] N.S. Shuman, M.A. Ochieng, B. Sztaray, T. Baer, J. Phys. Chem. A 112 (2008) 5647. [31] W.F. Schneider, T.J. Wallington, J. Phys. Chem. 98 (1994) 7448. [32] L.A. Curtiss, K. Raghavachari, P.C. Redfern, J.A. Pople, J. Chem. Phys. 106 (1997) 1063. [33] L.A. Curtiss, P.C. Redfern, K. Raghavachari, Chem. Phys. Lett. 499 (2010) 168. [34] D.J. Henry, M.B. Sullivan, L. Radom, J. Chem. Phys. 118 (2003) 4849. [35] A. Karton, J.M.L. Martin, Mol. Phys. 105 (2007) 2499. [36] N.L. Haworth, G.B. Bacskay, J. Chem. Phys. 117 (2002) 11175. [37] B. Ruscic, J.V. Michael, P.C. Redfern, L.A. Curtiss, K. Raghavachari, J. Phys. Chem. A 102 (1998) 10889. [38] S.P. Verevkin, V.N. Emel’yanenko, R. Notario, M.V. Roux, J.S. Chickos, J.F. Liebman, J. Phys. Chem. Lett. 3 (2012) 3454. [39] N.J. DeYonker, T. Grimes, S. Yockel, A. Dinescu, B. Mintz, T.R. Cundari, A.K. Wilson, J. Chem. Phys. 125 (2006) 10411. [40] N.J. DeYonker, B.R. Wilson, A.W. Pierpont, T.R. Cundari, A.K. Wilson, Mol. Phys. 107 (2009) 1107. [41] J. Olsen, O. Christiansen, H. Koch, P. Jorgensen, J. Chem. Phys. 105 (1996) 5082. [42] D. Cremer, Z. He, J. Phys. Chem. 100 (1996) 6173. [43] Y. He, D. Cremer, Theor. Chem. Acc. 105 (2000) 110.