Accurate Anharmonic Zero-Point Energies for Some Combustion

May 17, 2017 - Full-dimensional analytic potential energy surfaces based on CCSD(T)/cc-pVTZ calculations have been determined for 48 small combustion-...
1 downloads 9 Views 3MB Size
Download PDF
Article pubs.acs.org/JPCA

Accurate Anharmonic Zero-Point Energies for Some CombustionRelated Species from Diffusion Monte Carlo Lawrence B. Harding,* Yuri Georgievskii, and Stephen J. Klippenstein Chemical Sciences and Engineering Division, Argonne National Laboratory, Argonne, Illinois 60439, United States ABSTRACT: Full-dimensional analytic potential energy surfaces based on CCSD(T)/ cc-pVTZ calculations have been determined for 48 small combustion-related molecules. The analytic surfaces have been used in Diffusion Monte Carlo calculations of the anharmonic zero-point energies. The resulting anharmonicity corrections are compared to vibrational perturbation theory results based both on the same level of electronic structure theory and on lower-level electronic structure methods (B3LYP and MP2).

coordinate force field, as is the case with most implementations of vibrational self-consistent field and VCI. Probably the most straightforward, general approach to obtaining high-accuracy ZPEs is with Diffusion Monte Carlo (DMC) simulations.9 In principle, DMC will give the exact Born−Oppenheimer ZPE for a given potential energy surface (PES) regardless of the complexity of the potential surface. It does however require a full-dimensional (but not necessarily global) potential surface that accurately describes the region around the minimum. For many molecules, a quartic, internal coordinate force field may suffice. However, for fluxional molecules with multiple equivalent or near-equivalent minima, surface fitting becomes more challenging. For example, the ethyl radical has six equivalent minima separated by barriers of less than 50 cm−1. Recently Braams et al.10 developed a new technique for fitting PESs using permutationally invariant polynomials. The key idea is to use as fitting functions the complete set of internuclear distances. This allows one to straightforwardly impose permutational invariance of like atoms. The Bowman group11 has used this approach to fit a large number of global (reactive) surfaces with expansions in Morse variable functions of the internuclear distances. In the present approach, because we are only interested in fitting the PES in the vicinity of the minimum, we replace the Morse variable expansion with a simple Taylor series expansion and only impose permutational invariance on those atoms that can permute via internal rotations and/or inversions. For the ethyl radical mentioned above, even a quadratic, Taylor series expansion in the 21 internuclear distances, imposing permutational invariance among the three CH3 hydrogens and between the two CH2 hydrogens (62 symmetry-adapted polynomials),

1. INTRODUCTION High-accuracy, ab initio thermochemistry is an active area of research with a great deal of recent emphasis on increasing the accuracy of the underlying electronic structure calculations.1 There has been less emphasis on increasing the accuracy of the calculated zero point energies (ZPEs). This is in spite of the fact that almost 10 years ago Karton et al.2 noted that anharmonic corrections to ZPEs are one of the main factors limiting the accuracy of calculations on small organic molecules. One notable exception to this lack of emphasis on ZPEs is the recent study by Pfeiffer et al.3 comparing results of vibrational perturbation theory (VPT2) and vibrational configuration interaction (VCI) calculations on a set of 26 small molecules. Many early, high-accuracy methods simply scale harmonic ZPEs from relatively low-level calculations. For example, early Gn, Wn, and ccCA approaches all used scaled B3LYP harmonic frequencies. The W4 approach of Karton et al.4 and the focal point analysis (FPA) of Allen and co-workers5 do not have precise prescriptions for the ZPEs, typically getting them from the best-available literature sources. For example, in the FPA study of the quasi-linear molecule fulminic acid (HCNO), Schuurman et al.6 relied on a variational ZPE calculation reported by Pinnavaia et al.7 based on an MP2 anharmonic force field. Both the HEAT approach of Tajti et al.8 and the Feller−Peterson−Dixon1 (FPD) approach employ secondorder, rovibrational, perturbation theory (VPT2). In the HEAT approach, the perturbation theory is based on CCSD(T)/cc-pVQZ anharmonic force fields, while in the FPD approach, it is based on a combination of CCSD(T) harmonic ZPEs and MP2-based anharmonic corrections. For many species, these approaches are expected to yield very accurate ZPEs. However, for some species, particularly those in which the ground vibrational wave function is spread over multiple minima, perturbation theory can fail catastrophically. The same will be true of any approach that relies on a normal © 2017 American Chemical Society

Received: March 31, 2017 Revised: May 12, 2017 Published: May 17, 2017 4334

DOI: 10.1021/acs.jpca.7b03082 J. Phys. Chem. A 2017, 121, 4334−4340

Article

The Journal of Physical Chemistry A

where ΔXi is the maximal allowed displacement of the ith configurational coordinate (we use Z-matrix coordinates) and ξi is a random number between −1 and 1. The new configuration is accepted if its energy Enew < Elast. If Enew > Elast, the new configuration is accepted with the Boltzmann probability last new e(E −E )/T. The parameters ΔXi and T are adjusted for optimal performance. To avoid undesirable configurations (for example, undesirable conformational changes), the configurational coordinates Xi could be constrained within certain limits Xmin i and Xmax i . 2.3. Surface Fitting. The sampled points are least-squared fit to fourth-order Taylor series expansions of the complete set of symmetry-adapted internuclear distances. Typically, we include only points with energies less than 40−50 kcal/mol above the global minimum. The points are weighted in the least-squares fit such that the highest energy points have weights 105−106 times smaller than the lowest-energy point. In preliminary test calculations on hydrogen peroxide, we fit 30 000 points having an energy range of 3000 cm−1 using fourth-order expansions based on Taylor series, Simons−Parr− Finlan19 functions, and Morse variables. Each of these expansions involves 210 terms. The root-mean-square (RMS) errors in these three fits are 7, 45, and 52 cm−1, respectively. On the basis of these results, we have chosen to use Taylor series expansions for this study. We take methanol as a typical example of our fitting procedure. The quartic Taylor series expansion for methanol involves 905 symmetric polynomials. The surface is sampled and fit as described above using 50 000 points. An additional 75 000 points are then evaluated and used to test the accuracy of the fitted PES. The test points cover an energy range up to 15 000 cm−1 above the minimum. In Figure 1, we show a

will yield a PES having the desired six equivalent minima. Of course, this small fit will not be adequate for predicting anharmonicity corrections as it neglects the cubic and quartic terms, but even a full quartic expansion in these coordinates involves only 1774 terms. The primary focus of this paper then will be the determination of accurate, anharmonic corrections for some fluxional molecules relevant to combustion chemistry. We also report DMC ZPEs for some nonfluxional molecules for comparison to perturbation theory results. In the next section, we describe the electronic structure, surface fitting, and DMC approaches used. The results are then presented and discussed in sections 3 and 4 and briefly summarized in section 5.

2. THEORETICAL METHODS The general procedure followed in this work was to (i) sample the PES using high-level ab initio electronic structure methods, (ii) perform a least-squares fit of the ab initio energies using fourth-order Taylor series expansions in permutationally invariant combinations of internuclear distances, and (iii) carry out DMC calculations on the fitted PESs. In the next three subsections, we describe the electronic structure methods, the surface sampling and fitting procedures, and finally the details of the DMC calculations. 2.1. Electronic Structure Methods. The PESs were all fit to frozen core CCSD(T) energies. The CCSD(T) calculations employed Dunning12 correlation-consistent triple-ζ basis sets (cc-pVTZ). We note that although the CCSD(T)/cc-pVTZ level of theory is generally recognized as adequate for highaccuracy anharmonic corrections to ZPEs, it is less adequate for the harmonic component of the ZPE. High-accuracy harmonic ZPEs are often calculated with larger basis sets and/or explicitly correlated CCSD(T) methods and can often include corrections for core correlations and other effects.1,2 The focus of this paper is on high-accuracy anharmonic corrections and for this, we feel that CCSD(T)/cc-pVTZ represents a reasonable compromise between accuracy and computational cost. Except for some pathological cases discussed below, we expect that limitations in the electronic structure calculations contribute ∼10 cm−1 or less to the uncertainty in the anharmonic corrections,3 with similar or smaller errors inherent in the surface fitting and DMC calculations. All of the CCSD(T) calculations were carried out with the Molpro program package.13−16 In addition to the DMC calculations based on CCSD(T) PESs, we also report, for comparison, ZPEs from four lowerlevel calculations, G4 and three levels of second-order vibrational perturbation theory (VPT2). In G4 calculations, the ZPE is simply taken to be the B3LYP/6-31G(2df,p) harmonic frequency scaled by 0.9854. Our VPT2 calculations are based on CCSD(T)/cc-pVTZ, B3LYP/cc-pVTZ, and MP2/cc-pVTZ calculations. The B3LYP- and MP2-based VPT2 calculations were carried out using the GAUSSIAN09 program17 (revision D.01). The int = ultrafine option was used in the B3LYP analyses. The CCSD(T)/cc-pVTZ-based VPT2 calculations were done using Molpro, which is limited to asymmetric top molecules. 2.2. Surface Sampling. To sample statistically important regions of the PES, the Metropolis algorithm is used.18 The new configuration is obtained from the last one according to the following equation

Figure 1. Scatter plot of the ab initio energies vs fitted energies of 75 000 test points for the methanol potential surface.

scatter plot comparing the ab initio and fitted energies for the test points. The points within 1000 cm−1 of the minimum are fit with a RMS error of 1 cm−1. The RMS error for all 75 000 test points is 20 cm−1. 2.4. Diffusion Monte Carlo Calculations. The DMC approach to solving the Schroedinger equation was first proposed by Anderson.20 This approach is now widely used and has been recently reviewed by McCoy.9 In our calculations,

Xinew = Xilast + ΔXiξi 4335

DOI: 10.1021/acs.jpca.7b03082 J. Phys. Chem. A 2017, 121, 4334−4340

Article

The Journal of Physical Chemistry A

Table 1. ZPEs and Anharmonic Corrections (cm−1) for Nonfluxional, Closed-Shell Molecules Based on CCSD(T)/cc-pVTZ PESs differences in anharmonic corrections relative to DMC molecule H2O NH3 CH4 HCCH HCN HNC H2CCH2 trans-HNNH dioxirane H2CO CH2OO HOCN HONC HNCO CH3CN RMS difference a

harmonic ZPE

DMC total ZPE

DMC anharmonic corrections

B3LYP/VPT2

MP2/VPT2

CCSD(T)/VPT2

G4

4727 7575 9833 5784 3493 3407 11149 6214 7154 5857 6732 4692 4439 4655 9909

4653 7446 9690 5731 3462 3371 11006 6106 7063 5777 6643 4647 4391 4611 9789

−74 −129 −143 −53 −31 −36 −143 −108 −91 −80 −89 −45 −48 −44 −120

1 8 −3 23 −3 −4 −3 2 1 2 6 −5 −7 −5 −59 16.9(4.4)a

2 7 7 71 0 2 1 5 4 3 3 −5 −4 −2 −3 18.7(4.0)a

0

−30 −15 −10 76 129 88 21 0 −40 −24 107 43 38 32 5 57.6

−5 0 −4 −2 −1 −20 −7 −5 7.6

Numbers in parentheses do not include the two outliers HCCH and CH3CN.

Table 2. ZPEs and Anharmonic Corrections (cm−1) for Nonfluxional, Open-Shell Molecules Based on CCSD(T)/cc-pVTZ PESs differences in anharmonic corrections relative to DMC molecule 3

CH2 CH3 HCO HO2 trans-HOCO CH2CH CH2N 3 CH3N CH2CCH RMS difference

harmonic ZPE

DMC total ZPE

DMC anharmonic corrections

B3LYP/VPT2

MP2/VPT2

CCSD(T)/VPT2

G4

3805 6534 2851 3121 4619 8011 5521 8165 8911

3745 6465 2792 3070 4581 7902 5440 8031 8829

−60 −69 −59 −51 −38 −109 −81 −134 −82

−6 −15 3 2 −30 −13 −4 −28 0 15.4

−1 24 29 7 −24 −12 −26 13 32 21.2

−7

−4 −52 17 −11 −37 −5 6 −75 3 33.7

we start 10 000 “walkers” at the equilibrium geometry and propagate for 50 000 steps at 4 atu per step. The energy of each DMC “trajectory” will fluctuate around the exact ZPE as a function of the time steps. We then average the energies of each time step, ignoring the first 200 time steps to allow for initialization of the trajectory. We typically run nine such calculations for each molecule. For all of the molecules, this is sufficient to give a standard deviation of less than 10 cm−1 (for most, less than 5 cm−1). For most molecules, we repeat this procedure, substantially increasing the number of points, to verify that the results are converged with respect to the number of points in the fit. To give the reader an idea as to the computer resources required, for the methanol problem described above, a single DMC trajectory takes about the same amount of time as 15 000 ab initio energy evaluations. Thus, for this particular problem, the computational requirements for the nine DMC trajectories are similar to the requirements for all of the electronic structure calculations.

1 0 −27 −15 −7 22 14.8

pVTZ calculations, the total anharmonic ZPE from the DMC calculations, the anharmonic correction from the DMC calculations (the difference between the harmonic and anharmonic ZPEs), and the differences between the anharmonic corrections from VPT2 and G4 and those from the DMC calculations. We note that the anharmonic correction for G4 is not well-defined as the scale factor used in G4 calculations is designed to make up for deficiencies in both the harmonic and anharmonic components of the ZPE and, therefore, the appropriate harmonic reference does not exist. Here we define the G4 anharmonic correction to be the difference between the G4 ZPE and the harmonic CCSD(T)/ cc-pVTZ ZPE. This is equivalent to taking the difference between the total anharmonic ZPE from the DMC and G4 calculations. Alternatively, one could assume that the G4 scale factor accounts for only anharmonicity,21 in which case the G4 anharmonicity correction would be just −0.0146 times the unscaled B3LYP ZPE. For most closed-shell, nonfluxional molecules (Table 1), the DMC anharmonic corrections agree with the B3LYP-, MP2-, and CCSD(T)-based VPT2 results within the estimated ±10 cm−1 uncertainty in the DMC calculations. For the B3LYP and MP2/VPT2 results, the two primary exceptions are acetylene and CH3CN, which will be discussed in more detail in the next

3. RESULTS The results of the calculations are shown in Tables 1−4. In each case, we show the harmonic ZPE from CCSD(T)/cc4336

DOI: 10.1021/acs.jpca.7b03082 J. Phys. Chem. A 2017, 121, 4334−4340

Article

The Journal of Physical Chemistry A

Table 3. ZPEs and Anharmonic Corrections (cm−1) for Fluxional, Closed-Shell Molecules Based on CCSD(T)/cc-pVTZ PESs differences in anharmonic corrections relative to DMC molecule CH3CH3 NH2NH2 HOOH CH3OH NH2OH HCOOH HCNO CH3OCH3 NH2NO NH2CHO CH3NO CH3OOH RMS difference

harmonic ZPE

DMC total ZPE

DMC anharmonic corrections

B3LYP/VPT2

MP2/VPT2

16401 11773 5829 11311 8886 7461 4194 17554 7169 9978 9477 12045

16185 11569 5724 11132 8733 7363 4168 17328 7004 9847 9347 11862

−216 −204 −105 −179 −153 −98 −26 −226 −165 −131 −130 −183

−12 −8 12 13 −4 −3 77 −21 209 68 −4 5 67.9

−15 −1 3 8 −1 −2 −53 −16 −319 569 1 11 189.1

CCSD(T)/VPT2 −6 0 0 −6 −4 −385 −24 −1 −209 0 138.8

G4 −61 −35 15 −56 −32 −33 −8 −122 −30 −52 −24 −35 50.6

Table 4. ZPEs and Anharmonic Corrections (cm−1) for Fluxional, Open-Shell Molecules Based on CCSD(T)/cc-pVTZ PESs differences in anharmonic corrections relative to DMC molecule CH3CH2 3 CH3CH CH2OH CH3NH CH3CO NH2O NH2CO CH2NO HCCO CH3CHCH3 CH3OO cyclopropyl RMS difference

harmonic ZPE

DMC total ZPE

DMC anharmonic corrections

B3LYP/VPT2

MP2/VPT2

13052 10415 8244 10760 9469 5896 7149 6946 4094 19354 9506 14731

12826 10211 8010 10590 9335 5732 7111 6867 4056 19050 9332 14505

−226 −204 −234 −170 −134 −164 −38 −79 −38 −304 −174 −226

27 11 29 0 −5 −104 −38 −6 −13 −19 41 15 37.0

18 33 54 13 6 376 116 −48 24 −18 48 16 117.3

CCSD(T)/VPT2 22 23 43 0 0 −35 −11 −14 40 25.8

G4 −35 −40 81 −52 −28 32 −7 27 40 −58 −59 96 51.9

the small tunneling splitting associated with the torsional motion. One exception is hydrogen peroxide, where the trans torsional barrier is predicted to be 365 cm−1 and the torsional harmonic frequency is 372 cm−1. Pfeiffer et al.3 do not report VCI calculations on hydrogen peroxide. As a check of our DMC calculations, we carried out a VCI calculation on hydrogen peroxide at the CCSD(T)/cc-pVTZ level, obtaining a ZPE of 5729 cm−1, 5 cm−1 higher than the DMC result. While this difference is within the expected uncertainty in the DMC calculations, it is also interesting to note that it is ∼1/2 the known tunneling splitting in hydrogen peroxide.22 In the next section, we examine in some detail those cases where there are large differences between the VPT2 and DMC results.

section. For the CCSD(T) VPT2 results, a significant discrepancy exists for HOCN. Note that we typically find no significant difference between the B3LYP, MP2, and CCSD(T) VPT2 corrections, while the differences for the G4 corrections are on average more than 10 times larger. For open-shell, nonfluxional molecules (Table 2), the differences between the DMC corrections and VPT2 are 2−4 times larger. The RMS differences, relative to DMC, for B3LYP/VPT2, MP2/VPT2, CCSD(T)/VPT2, and G4 are 15, 21, 15, and 34 cm−1, respectively. For fluxional molecules, we find, as expected, much larger differences between the DMC and VPT2 results. Here we loosely define fluxional molecules to include molecules, such as ethane and hydrogen peroxide, having relatively low-barrier internal rotations. For the 12 closed-shell fluxional molecules in Table 3, the RMS differences relative to DMC for B3LYP/ VPT2, MP2/VPT2, CCSD(T)/VPT2, and G4 are 68, 189, 139, and 51 cm−1, respectively. Similar differences are found for the open-shell, fluxional molecules in Table 4. Our results are in good agreement with the VPT2 and VCI results of Pfeiffer et al.3 for the 10 molecules that we have in common. One might question why VPT2 (and VCI) calculations do so well for molecules such as CH3CH3, NH2NH2, HOOH, CH3OH, and NH2OH, with multiple torsional minima. For most of these molecules, the torsional barriers are significantly larger than the torsional frequencies, and therefore, the VPT2 and VCI calculations are just missing

4. DISCUSSION 4.1. HCCH. The DMC and VPT2 anharmonicities for acetylene differ widely, 5−10 times the RMS deviation for the set of 15 closed-shell, nonfluxional molecules in Table 1. It has been noted 3 that the VPT2 implementation in the GAUSSIAN09 program is problematic for certain highsymmetry molecules including methane and all linear molecules such as acetylene. There have been two previous high-level determinations of the acetylene zero point by Zou et al.23 and Xu et al.24 In both cases, the vibrational problem was solved using full-dimensional, converged, quantum mechanical calculations. Zou et al. fit a full-dimensional PES using CCSD(T)/ 4337

DOI: 10.1021/acs.jpca.7b03082 J. Phys. Chem. A 2017, 121, 4334−4340

Article

The Journal of Physical Chemistry A

that their predicted ΔH°f,0(HCNO) should be increased from +40.9 to +41.7 kcal/mol. 4.4. NH2NO. The CCSD(T)/cc-pVTZ calculations predict the equilibrium geometry of NH2NO to be nonplanar with HNNO and HHNN dihedral angles of 14 and 150°, respectively, and a barrier to planarity of only 40 cm−1. The barriers to rotation about the NN bond are predicted to be much larger, 5593 and 6241 cm−1, depending on the direction of rotation. B3LYP calculations predict the molecule to be planar, while MP2 calculations predict a barrier to planarity of only 19 cm−1. The B3LYP and MP2/VPT2 calculations fail catastrophically, in one case giving the wrong sign for the anharmonicity correction and the other in error by a factor of 3. 4.5. NH2CHO. High-level calculations predict formamide to be planar.26 However, at the level of theory used here, formamide is predicted to have a slightly nonplanar equilibrium geometry, a very small barrier to planarity, and larger barriers to internal rotation about the CN bond. At the CCSD(T)/ccpVTZ level, the barrier to planarity and the two barriers to internal rotation are 5, 5372, and 5926 cm−1, respectively. The B3LYP calculations predict a planar, equilibrium geometry, as in NH2NO. All three VPT2 calculations fail catastrophically. The VPT2/B3LYP and VPT2/CCSD(T) results have the correct sign, but the magnitudes are more than a factor of 2 too small or too large, respectively. The VPT2/MP2 calculations have the wrong sign, and the magnitude is off by more than a factor of 3. 4.6. CH3CH2. The internal rotation−CH2 inversion potential of the ethyl radical exhibits six equivalent minima separated by six equivalent saddle points, as shown in Figure 2. In this plot,

aug-cc-pVTZ calculations, while Xu et al. scaled the Zou PES to match certain experimentally observed vibrational levels. These calculations yield ZPEs of 5716 and 5749 cm−1, respectively, in good agreement with our own result of 5731 cm−1. It is worth noting that the bend frequencies in acetylene are unusually sensitive to the basis set. Upon going from the cc-pVTZ to the cc-pVQZ basis set, the harmonic ZPE increases by 28 cm−1 and the anharmonic ZPE increases by 15 cm−1. Our best estimate for the anharmonic ZPE of acetylene (DMC based on a CCSD(T)/cc-pVQZ PES) is 5746 cm−1, within 3 cm−1 of the experimentally derived result noted above. 4.2. CH3CN. Among the closed-shell, nonfluxional molecules, CH3CN is also an outlier, with the DMC and VPT2/ B3LYP anharmonicities differing by 59 cm−1. It is not clear why this is the case, although we note that the DMC and VPT2/ MP2 results are in very good agreement. 4.3. HCNO. Fulminic acid, HCNO, is a quasi-linear molecule thought to be an important intermediate in NOreburning mechanisms to control NOx emissions.25 At the CCSD(T)/cc-pVTZ level, the HCN angle is 157° and the barrier to linearity is only 42 cm−1. Both of these quantities are quite sensitive to the level of theory. For example, increasing the basis set to cc-pVQZ increases the equilibrium angle to 165° and decreases the barrier to just 7 cm−1. The errors in the VPT2/B3LYP and MP2 anharmonicities relative to DMC are similar in magnitude but opposite in sign, while the VPT2/ CCSD(T) calculation fails catastrophically. The B3LYP calculations incorrectly predict the molecule to be linear, and therefore, the GAUSSIAN09 VPT2 calculations may be open to question, as noted above for acetylene. HCNO is a good example of the potential pitfalls of combining the harmonic zero point determined with a large basis set calculation and an anharmonic correction from a smaller basis set calculation. The results of DMC calculations with cc-pVDZ, cc-pVTZ and cc-pVQZ basis sets are shown in Table 5. The total anharmonic ZPE from the DMC calculations Table 5. Basis Set Dependence of Harmonic and Anharmonic ZPEs (cm−1) for HCNO basis set

harmonic ZPE

anharmonic ZPE

anharmonic correction

cc-pVDZ cc-pVTZ cc-pVQZ

4244 4194 4159

4169 4168 4183

−75 −26 +24

is found to vary only slightly (15 cm−1) as a function of the basis set. The anharmonic corrections however vary by ∼100 cm−1 from −75 cm−1 with the smallest basis set to +24 cm−1 with the largest basis set. The most thorough previous estimate for the ZPE of HCNO is the variational calculation of Pinnavaia et al.7 Their calculations, based on four different internal coordinate force fields all fit to MP2/DZP energies, yielded ZPEs ranging from 3904 to 4115 cm−1, 50−250 cm−1 below the current estimate. Some of this difference might be attributable to the differences in the electronic structure theory employed. However, their harmonic ZPE is significantly higher than the present calculations, 4330 vs 4194 cm−1. This means that their anharmonic correction is 10−20 times that from the DMC calculations. As noted above, the FPA of Schuurman et al.6 employed a ZPE of 3901 cm−1 (based on Pinnavaia et al.7). This is 267 cm−1 below the present DMC result and implies

Figure 2. Two-dimensional plot of the PES for the ethyl radical as a function of the internal rotation about the CC bond and inversion of the CH2 moiety. See the text for the definition of the coordinates.

the inversion angle is defined relative to a pseudo-C3 axis, P. This pseudo-C3 axis is defined in such a way that if Ha, Hb, and Ca are the atoms in the CH2 moiety, then the Ha−Ca−P, Hb− Ca−P, and Cb−Ca−P angles are all equal. The inversion angle is then 90° minus this angle. Thus, an inversion angle of 0° corresponds to a planar CH2 moiety. The C−C rotation angle is then defined to be the dihedral angle Hm−Cb−Ca−P, where Hm is one of the methyl hydrogens. The saddle points separating the minima are predicted to lie 32 cm−1 above the minima. Despite these very small barriers and the highly anharmonic nature of the PES, the VPT2 4338

DOI: 10.1021/acs.jpca.7b03082 J. Phys. Chem. A 2017, 121, 4334−4340

Article

The Journal of Physical Chemistry A

set calculation and an anharmonic correction from a smaller basis set calculation. The results of DMC calculations with ccpVDZ, cc-pVTZ, and cc-pVQZ basis sets are shown in Table 6.

calculations do not fail catastrophically. The errors in the VPT2 ZPEs are of the same order of magnitude as that found for nonfluxional, open-shell molecules. 4.7. CH2OH. The internal rotation−CH2 inversion potential of the hydroxy methyl radical exhibits two pairs of equivalent minima, as shown in Figure 3. The coordinates are analogous to

Table 6. Basis Set Dependence of Harmonic and Anharmonic ZPEs (cm−1) for NH2CO basis set

harmonic ZPE

anharmonic ZPE

anharmonic correction

cc-pVDZ cc-pVTZ cc-pVQZ

7235 7149 7189

7071 7111 7128

−164 −38 −61

As with HCNO, the results show that although the anharmonic zero point varies smoothly with basis set the harmonic zero points (and thus the anharmonic corrections) do not. The problem here is that with the cc-pVDZ and cc-pVTZ basis sets the molecule is nonplanar, having a double-minimum potential with a small barrier to inversion (61.5 and 0.03 cm−1, respectively). With the cc-pVQZ basis set, however, the molecule is planar. As a result, if one were to combine the harmonic zero point from the cc-pVQZ calculation with an anharmonic correction from either the cc-PVDZ or cc-PVTZ calculation, the resulting estimated anharmonic zero point would be less accurate than the anharmonic zero point from either the cc-pVDZ or cc-pVTZ calculations alone.

Figure 3. Two-dimensional plot of the PES for the hydoxymethyl radical as a function of the internal rotation about the CO bond and inversion of the CH2 moiety. See the text for the definition of the coordinates.

5. SUMMARY One important conclusion from this study for the nonfluxional molecules is that there appears to be no advantage to using the more expensive MP2-based VPT2 over the much cheaper B3LYP-based VPT2. In fact, the B3LYP results are somewhat more accurate than the MP2 results. Even the much more expensive CCSD(T)-based VPT2 method shows only very modest improvement over B3LYP VPT2. Molecules with simple hindered rotations, even those with low barriers, are adequately treated with VPT2. Presumably, the VPT2 calculations neglect the effect of the multiple minima torsional potentials; however, because the torsions tend to be low-frequency modes, small errors in these modes do not greatly impact the ZPE. Quasi-linear and quasi-planar molecules are more problematic for VPT2. In several cases, the perturbation calculations are found to fail catastrophically. For these molecules, simply scaling harmonic ZPEs, as done in G4, tends to be more accurate than VPT2 calculations.

those described for the ethyl radical above. Within each pair, the minima are connected by inversion of the CH2 moiety. The two pairs are in turn connected by rotation about the CO bond. The inversion barrier is predicted to be 139 cm−1, while the barrier to internal rotation is predicted to be 1670 cm−1. The most extensive previous calculation on the hydroxyl methyl radical was by Marenich and Boggs.27,28 In these calculations, a partial, quartic, internal coordinate force field was fit to ∼1800 ab initio points at the CCSD(T)/cc-pVTZ level, with the quadratic part of the force field scaled to account for limitations of the cc-pVTZ basis set. They further assumed a separation between in-plane and out-of-plane modes. They then performed a variational calculation to obtain a ZPE of 7900 ± 7 cm−1 (the uncertainty comes from analysis of the limitations in the size of the basis set used in the electronic structure calculations). This is significantly different than the DMC result of 8010 cm−1. The most likely explanations for the difference are the use of an incomplete quartic force field and/ or the assumed separation of in-plane and out-of-plane modes. 4.8. NH2O. NH2O is another quasi-planar molecule having a barrier to inversion of 40 cm−1. It is thought to be one of the dominant products of the reaction of OH with HCNO,25 a key step in NO-reburning. VPT2/B3LYP correctly predicts a negative anharmonicity correction, but the magnitude of the correction is almost a factor of 2 too large. VPT2/MP2 incorrectly predicts a positive anharmonicity correction. 4.9. NH2CO. Similar to NH2NO and NH2CHO, NH2CO is slightly nonplanar, having a predicted barrier to planarity of less than 1 cm−1 (at the CCSD(T)/cc-pVTZ level) and a much larger barrier to rotation about the NC bond, 5623 cm−1. As with NH2NO and NH2CHO, the B3LYP and MP2/VPT2 calculations have large errors with opposite signs. NH2CO is another good example of the potential pitfalls of combining a harmonic zero point determined with a large basis

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Lawrence B. Harding: 0000-0002-5867-4716 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Anne McCoy for providing us with a DMC code. This work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences, under Contract Number DEAC02-06CH11357. 4339

DOI: 10.1021/acs.jpca.7b03082 J. Phys. Chem. A 2017, 121, 4334−4340

Article

The Journal of Physical Chemistry A

■

(23) Zou, S.; Bowman, J. M. A New Ab Initio Potential Energy Surface Describing Acetylene/Vinylidene Isomerization. Chem. Phys. Lett. 2003, 368, 421−424. (24) Xu, D.; Guo, H.; Zou, S.; Bowman, J. M. A scaled Ab Initio Potential Surface for Acetylene and Vinylidene. Chem. Phys. Lett. 2003, 377, 582−588. (25) Feng, W.; Meyer, J. P.; Hershberger, J. F. Kinetics of the OH Plus HCNO Reaction. J. Phys. Chem. A 2006, 110, 4458−4464. (26) Alessandrini, S.; Puzzarini, C. Structural and Energetic Characterization of Prebiotic Molecules: The Case Study of Formamide and Its Dimer. J. Phys. Chem. A 2016, 120, 5257−5283. (27) Marenich, A. V.; Boggs, J. E. A Variational Study of Nuclear Dynamics and Structural Flexibility of the CH2OH Radical. J. Chem. Phys. 2003, 119, 3098−3105. (28) Marenich, A. V.; Boggs, J. E. Structural and Thermochemical Properties of the hydroxymethyl (CH2OH) Radical: A High Precision Ab Initio Study. J. Chem. Phys. 2003, 119, 10105−10114.

REFERENCES

(1) Peterson, K. A.; Feller, D.; Dixon, D. A. Chemical Accuracy in Ab Initio Thermochemistry and Spectroscopy: Current Strategies and Future Challenges. Theor. Chem. Acc. 2012, 131, 1079−1099. (2) Karton, A.; Ruscic, B.; Martin, J. M. L. Benchmark Atomization Energy of Ethane: Importance of Accurate Zero-Point Vibrational Energies and Diagonal Born-Oppenheimer Corrections for a ‘Simple’ Organic Molecule. J. Mol. Struct.: THEOCHEM 2007, 811, 345−353. (3) Pfeiffer, F.; Rauhut, G.; Feller, D.; Peterson, K. A. Anharmonic Zero Point Vibrational Energies: Tipping the Scales in Accurate Thermochemistry Calculations. J. Chem. Phys. 2013, 138, 044311. (4) Karton, A.; Rabinovich, E.; Martin, J. M. L.; Ruscic, B. W4 Theory for Compuational Thermochemistry: In Pursuit of Confident Sub-kJ/mol Predictions. J. Chem. Phys. 2006, 125, 144108. (5) East, A. L. L.; Allen, W. D. The Heat of Formation of NCO. J. Chem. Phys. 1993, 99, 4638. (6) Schuurman, M. S.; Muir, S. R.; Allen, W. D.; Schaefer, H. F. Toward Subchemical Accuracy in Computational Thermochemistry: Focal Point Analysis of the Heat of Formation of NCO and [H,N,C,O] Isomers. J. Chem. Phys. 2004, 120, 11586−11599. (7) Pinnavaia, N.; Bramley, M. J.; Su, M.-D.; Green, W. H.; Handy, N. C. A Study of the Ground Electronic State of the Isomers of CHNO. Mol. Phys. 1993, 78, 319−343. (8) Tajti, A.; Szalay, P. G.; Csaszar, A. G.; Kallay, M.; Gauss, J.; Valeev, E. F.; Flowers, B. A.; Vazquez, J.; Stanton, J. F. HEAT: High Accuracy Extrapolated Ab Initio Thermochemistry. J. Chem. Phys. 2004, 121, 11599−11613. (9) McCoy, A. B. Diffusion Monte Carlo Approaches to Investigating the Structure and Vibrational Spectra of Fluxional Systems. Int. Rev. Phys. Chem. 2006, 25, 77−107. (10) Braams, B. J.; Bowman, J. M. Permutational Invariant Potential Energy Surfaces in High Dimensionality. Int. Rev. Phys. Chem. 2009, 28, 577−606. (11) Bowman, J. M.; Czako, G.; Fu, B. High-Dimensional Ab Initio Potential Energy Surfaces for Reaction Dynamics Calculations. Phys. Chem. Chem. Phys. 2011, 13, 8094−8111. (12) Dunning, T. H., Jr. Gaussian-Basis Sets for use in Correlated Molecular Calculations. 1. The Atoms Boron Through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (13) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manbt, F. R.; Schutz, M.; et al. MOLPRO, version 2010.1, a package of ab initio programs. http://www.molpro.net (2015). (14) Hampel, C.; Peterson, K.; Werner, H.-J. A Comparison of the Efficiency and Accuracy of Quadratic Configuration-Interaction (QCISD), Coupled Cluster (CCSD), and Brueckner Coupled Cluster (BCCD) Methods. Chem. Phys. Lett. 1992, 190, 1−12. (15) Deegan, M. J. O.; Knowles, P. J. Perturbative Corrections to Account for Triple Excitations in Closed and Open-Shell CoupledCluster Theories. Chem. Phys. Lett. 1994, 227, 321−326. (16) Knowles, P. J.; Hampel, C.; Werner, H.-J. Coupled-Cluster Theory for High-Spin, Open-Shell Reference Wave-Functions. J. Chem. Phys. 1993, 99, 5219−5227. (17) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; et al. GAUSSIAN 09; Gaussian, Inc.: Wallingford, CT, 2010. (18) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. Equation of State Calculations by Fast Computing Machines. J. Chem. Phys. 1953, 21, 1087−1092. (19) Simons, G.; Parr, R. G.; Finlan, J. M. New Alternative to the Dunham Potential for Diatomic Molecules. J. Chem. Phys. 1973, 59, 3229−3234. (20) Anderson, J. B. A Random Walk Solution of the Schroedinger Equation: H+3. J. Chem. Phys. 1975, 63, 1499−1503. (21) Kesharwani, M. K.; Brauer, B.; Martin, J. M. L. Frequency and Zero-Point Vibrational Energy Scale Factors for Double-Hybrid Density Functionals (and Other Selected Methods): Can Anharmonic Force Fields be Avoided? J. Phys. Chem. A 2015, 119, 1701−1714. (22) Harding, L. B. Theoretical Studies of the Hydrogen Peroxide Potential Surface 1. An Ab Initio Anharmonic Force Field. J. Phys. Chem. 1989, 93, 8004−8013. 4340

DOI: 10.1021/acs.jpca.7b03082 J. Phys. Chem. A 2017, 121, 4334−4340