Variational calculation of the rate of dissociation of ... - ACS Publications

Nov 1, 1991 - Elisabeth A. Wade, Axel Mellinger, Melonie A. Hall, and C. Bradley Moore ... Donald G. Truhlar , Bruce C. Garrett , Stephen J. Klippenst...
0 downloads 0 Views 1017KB Size
Download PDF
9882

J. Phys. Chem. 1991, 95, 9882-9889

Variational Calculation of the Rate of Dissociation of CH,CO into 'CH, and CO on an ab Initio Determined Potential Energy Surface Jenwei Yu and Stephen J. Klippenstein*

Chemistry Department, Case Western Reserve University, Cleveland, Ohio 44106 (Received: June 3, 1991; In Final Form: July 22, 1991) An ab initio potential energy surface is developed and used in microcanonical variational RRKM theory calculations of the rate of dissociation of CHzCO into 'CH, and CO. The CC bond length is first demonstrated to provide a better description of the reaction coordinate than the distance separating the centers of mass of the two fragments via calculations on a model potential energy surface. The potential energy is then evaluated at more than 200 different configurations for each of four different bond lengths, with these configurations spanning the ranges of the interfragment orientational coordinates. The four bond lengths are chosen to span the expected range of transition-state location for the energies at which experimentally determined rate constants are available. In the ab initio calculations a 6-31G* basis set is used and electron correlation effects are included through the application of second-order Maller-Plesset perturbation theory. Comparison is made between theoretical results employing analytical fits to the ab initio determined energy points of varying accuracy and experimental results. Reasonable agreement with experiment is obtained, particularly when approximate correction factors for the reactant anharmonicities and the usual ab initio overestimate of force constants are included.

I. Introduction The dissociation of CHzCO has been the focus of considerable recent experimental interest.'-7 There are two low-energy dissociation channels corresponding to the lowest singlet and triplet electronic states of CHz with the triplet state lower in energy by 3 147 cm-IS8 The dependence of the total dissociation rate constant on excess energy in the range from approximately 500 to 6000 cm-I above the threshold for the singlet channel has been determined in the time-resolved experimental studies of Zewail and co-workers.' Meanwhile, the studies of Moore and co-workers have yielded detailed product rovibrational state distributions for 1CH224and C0.596 The studies of Moore and co-workers have also yielded dissociation rate constantsS at energies ranging from below the singlet threshold to a few hundred cm-' excess and singlet/triplet branching ratios6 at energies ranging from the singlet threshold to 2500 cm-l excess energy. The combination of the singlet/triplet branching ratios with the measured total dissociation rate constants yields individual singlet and triplet contributions to the dissociation rate.6 These combined experimental studies'-7 provide one of the most detailed sets of information regarding a unimolecular dissociation. As such, they provide an excellent test for statistical theories of unimolecular reactions,"2 probing both the "inner" transition-state region and the dynamics from this inner transition-state region on to products. The inner transition-state region is expected on the basis of previous calculations to be at CC separation distances on the order of 2-3 A.I3J4 This previous theoretical work13,14 (1) Potter, E. D.; Gruebele, M.; Khundkar, L. R.; Zewail, A. H. Chem. Phys. Leu. 1989, 164, 463. (2) Chen, I.-C.; Green, W. H., Jr.; Moore, C. B. J . Chem. Phys. 1988,89, 314. (3) Green, W.H.,Jr.: Chen, 1.4.; Moore, C. 9. Eer. Eumenges. Phys. Chem. 1988. 92. 389. (4) Green, W. H., Jr.; Mahoney, A. J.; Cheng, Q.-K.; Moore, C. B. J. Chem. Phys. 1991, 94, 1961. ( 5 ) Chen, 1 . C ; Moore, C. 9. J. Phys. Chem. 1990, 94, 263, 269. (6) Kim, S. K.; Choi, Y. S.; Pibel, C. D.; Zheng, Q.-K.; Moore, C. B. J. Chem. Phys. 1991, 94, 1954. (7) Lovejoy, E. R.; Kim, S. K.; Alvarez, R. A,; Moore, C. B. J . Chem. Phys. 1991, 95, 4081. (8) Bunker, P. R.; Jensen, P.; Kraemer, W. P.; Beardsworth, R. J . Chem. Phys. 1986, 85, 3124. (9) Marcus, R. A. J . Chem. Phys. 1952, 20. 359; 1%5,43, 2658. (IO) Pechukas, P.: Light, J. C. J . Chem. Phys. 1%5,42,3281. Pechukas, P.; Rankin, R.; Light, J. C. J . Chem. Phys. 1966, 44, 794. Klots, C. E. J . Phys. Chem. 1971, 75, 1526. ( I I ) Quack, M.; Troe, J. Eer. Eunrenges. Phys. Chem. 1974,78,240,1975, 79, 170, 469. (12) Truhlar, D. G.; Isaacson, A. D.; Garrett, B. C. In Theory of Chemical Reaction Dynamics; Baer, M.. Ed.; CRC: Boca Raton, FL, 1985; Vol. 4, p 65. Hase, W. L. Acc. Chem. Res. 1983, 16, 258. Truhlar, D. G.; Hase, W. L.;Hynes, J. T. J . Phys. Chem. 1983,87, 2664. (13) Klippenstein, S. J.; Marcus, R. A. J. Chem. Phys. 1989, 91, 2280.

has indicated that the various aspects of the experimental results for the 'CHz channel can be understood in terms of a recently developed implementation of variational RRKM theory.lsJ6 However, the potential energy surface used in these theoretical calculations was based on a simple model involving a combination of a bonding potential for the CC interaction and nonbonding van der Waals potentials for the remaining interactions. Additionally, more recent theoretical workI7J8has indicated that a bond length reaction coordinate may in general provide a better description of the reaction coordinate than the center-of-mass separation distance reaction coordinate used in the previous theoretical studies of this dissociation process.13J4 Thus, in this paper a further analysis of this dissociation process utilizing ab initio determined potential energy surfaces is provided. Consideration is also given to the effect of the variation of the assumed reaction coordinate between the two above- mentioned choices. Previous a b initio studies for the dissociation of C H 2 C 0 have been reported by Yamabe and Morokuma19 and by Allen and Schaefer.20.2' The studies by Allen and Schaefer20.21have considered the equilibrium parameters for various excited electronic statesZoand the transition-state parameters for the triplet-state dissociatiomZ1The focus of the present paper will be with regard to the dissociation on the lowest singlet surface, and thus the latter calculations20*z'will not be considered further here. In their study Yamabe and Morokuma first indicated that a C, dissociation is symmetry forbidden and thus of high energy.19 Correspondingly, the lowest energy reaction path must involve a bending motion. They then proceeded to examine the reaction energetics for both bent in plane and bent out of plane reaction paths.Ig Using a STO-3G minimal basis set and the restricted Hartree-Fock method, they demonstrated that the bent out of plane reaction path provided the lowest energy route to products for the ground singlet state. They also suggested that there should be little or no barrier for the reverse recombination of 'CHz and CO to form ketene.19 This absence of a reverse barrier is in good agreement with the experimental results of Moore and co-workers.2" In the present work the energetics for the dissociation of CHzCO into 'CHz and CO is examined in greater detail using a larger (14) Klippenstein, S. J.; Marcus, R. A. J . Chem. Phys. 1990, 93, 2418. (15) Wardlaw, D. M.; Marcus, R. A. Chem. Phys. Leu.1984, 110, 230. Wardlaw, D. M.; Marcus, R. A. J. Chem. Phys. 1985,83, 3462. Wardlaw, D. M.; Marcus, R. A. J . Phys. Chem. 1986, 90, 5383. Wardlaw, D. M.; Marcus, R. A. Adu. Chem. Phys. 1988, 70, 231. (16) Klippenstein, S. J.; Marcus, R. A. J . Phys. Chem. 1988, 92, 3105, 5412. (17) Klippenstein, S. J. Chem. Phys. Leu. 1990, 170. 71. (18) Klippenstein, S. J. J . Chem. Phys. 1991, 94, 6469. (19) Yamabe, S.; Morokuma, K. J. Am. Chem. Soc. 1978, 100, 7551. (20) Allen, W. D.; Schaefer, H. F., 111. J. Chem. Phys. 1986, 84, 2212. (21) Allen, W. D.; Schaefer, H. F., 111. J . Chem. Phys. 1988, 89, 329.

0022-365419 112095-9882%02.50/0 0 1991 American Chemical Society

Dissociation of C H 2 C 0 into 'CH2 and CO basis set and including correlation effects through the application of second-order Maller-Plesset perturbation theory. Detailed consideration is given to the potential energy surface in the transition-state region corresponding to C-C bond lengths ranging from 2.2 to 3.1 A to allow for an accurate determination of the rate constant for the singlet dissociation channel within variational RRKM theory. Attention is focused on four specific values of the reaction coordinate with a separate potential energy surface developed for the motions "orthogonal" to the reaction coordinate for each value. Within this approach the reaction coordinate must of course be specified before potential energy surfaces may be developed for the motions orthogonal to it. Thus, to begin, in section 11 a comparison between the results of variational calculations based on bond length and center-of-mass separation distance reaction coordinates are compared using the model surface described in ref 13. This comparison leads to a choice of the bond length as the reaction coordinate in the subsequent a b initio calculations. The methodology for the current a b initio calculations and potential fitting for the transition-state region are described in section 111. In section IV the results of the present a b initio and variational RRKM theory calculations are presented and discussed. A comparison of the calculated rate constants with the corresponding experimental results is also provided there. The present calculational method may also be used to consider the product state distributions, although this extension has not been performed here. Some concluding remarks are made in section V. 11. Comparison of Reaction Coordinates Dissociation reactions for which there is little or no reverse barrier present an interesting challenge for statistical theory treatments. The lack of a reverse barrier makes it difficult to know the location of the transition state a priori. In fact, the transition-state location is often found to depend quite strongly on the excess energy and total angular momentum. In addition, the location of this transition state may occur at separation distances where some of these modes are not properly described as either harmonic oscillators or free rotors. These facts suggest that a variational treatment which includes the hindered rotational nature of these modes may be required for accurate results to be obtained. A few years ago Wardlaw and Marcus developed a variational RRKM method that provides an accurate treatment within classical statistical mechanics of the coupled hindered rotor nature of the transition state.I5 Previous statistical calculations of the rate constant and product state distributions for the dissociation of C H 2 C 0 were described in refs 13 and 14 and provided good agreement with experiment. These calculations were based on a variation of the original method of Wardlaw and Marcus in which the phase space integral for the classical density of states for the transitional modes is treated in terms of standard Euler-angle coordinates and their conjugate momenta in place of action-angle coordinates.I6 One of the inherent assumptions in these specific calculations (and in the method of Wardlaw and Marcus in general) is that the reaction coordinate is given by the distance separating the centers of mass of the two departing fragments. A recent study provided a method for implementing a more arbitrary reaction coordinate into variational RRKM theory calculation^.'^^^^ In particular, a detailed scheme for determining the transition-state number of available states for an assumed reaction coordinate given by the distance between any two fixed points in the separate fragments was described. The bond length and center-of-mass separation distance reaction coordinates correspond to two specific limits of this more general reaction coordinate. A consideration of the calculated rate constant for the dissociation of N C N O into N C and NO indicated that the bond length reaction coordinate is to be preferred over the center-of-mass separation distance reaction coordinate for that r e a ~ t i o n . ~ ~ . ~ * Here, consideration is given to the effect of the choice of reaction coordinate on the statistically calculated rate constant for the dissociation of C H 2 C 0 into 'CHI and CO. The model surface described in ref 13 is used since the primary purpose of these initial

The Journal of Physical Chemistry, Vol. 95, No. 24, 1991 9883

TABLE I: Comparison of Mm Calculations Using a Model Potential Enerzy Surface

energy, cm-'

IO 20 40 100 200 500 1000 2000 3500 5000 6000

N(=)" N&, R&? A NLdr RLd,bA N& / N L d 6.4 (2), 3.65 8 5 , 3.13 7.6 14 1.3 (2). 3.1 1 5.6 78 7.3 (2), 3.57 2.1(2),3.05 5.1 4.3(2) l.l(3),3.54 6.8 (2). 2.96 3.4 3.1 (3) 2.3 (3), 3.42 2.1 (3), 2.89 2.7 1.4 (4) 5.5 (3), 3.39 1 . 1 (4), 2.62 2.0 9.0 (4) 2.1 (4), 3.17 3.7 (4), 2.53 2.0 3.7 ( 5 ) 7.2 (4), 3.00 1.6 ( 5 ) . 2.49 1.9 1.7 (6) 2.9 ( 5 ) , 2.94 6.8 ( 5 ) , 2.37 1.9 7.8 (6) 1.3 (6), 2.84 2.2 (6), 2.32 1.6 2.5 (7) 3.6 (6), 2.75 4.2 (6), 2.29 1.5 5.0 (7) 6.2 (6), 2.71

"A'(-) denotes the number of available states at infinite separation of the two fragments. Numbers in parentheses denote the power of IO. h C M and bond subscripts denote calculations assuming a reaction coordinate given by the distance separating the two fragment centers-ofmass and the bond length, respectively. N' and R' denote the minimum number of available states and the value of the reaction coordinate at the inner transition state, respectively.

calculations is to provide an indication of the preferred choice of reaction coordinate to be used in the determination of the ab initio potential energy surfaces. In Table I NE,, the minimum in the fixed energy E and total angular momentum J number of available states as a function of the value of the reaction coordinate, is given for both of the reaction coordinates in the inner transition-state region. Also given there is R*, the value of the reaction coordinate at the inner minimum in NEj(R), for each of these two reaction coordinates. A final quantity presented in Table I is NEA-), the number of available states when the two product fragments are infinitely separated. All of these calculations were performed for a Jvalue of 3, and the classical treatment of the transitional modes described in ref 18 was used in every instance. The results presented in Table I indicate that the bond length once again provides an improved description of the reaction coordinate in the inner transition-state region. The difference between the two reaction coordinates is particularly pronounced at the lowest energies. At energies below -30 cm-I the correct NLJ for this model potential energy surface is the NEj(-) value, and as a result the maximum effect of the choice of reaction coordinate on the rate constant for this dissociation process is approximatel a factor of 5 . The R' values for the bond length are 0.4-0.5 lower than those for the center-of-mass separations, as might be expected on the basis of the distance from the CO center of mass to its C atom. The location of the inner transition state is seen to decrease by approximately 1 A from the singlet threshold up to 6000 cm-I of excess energy. Statistical estimates of kEJ,the E and J conserved singlet dissociation rate constant, may be obtained from the RRKM expression

K

A harmonic approximation is made here for the vibrational contribution to the reactant density of states pEJ using the frequencies listed in ref 13. In keeping with the classical treatment of the transitional modes in the determinations of NiJ the rotational contribution to the reactant density of states was also treated classically here. The variational minimum of NEAR)will be taken as the minimum of the inner and outer transition-state minima.I6 All other details of the calculations are as described in ref 13. In Figure 1 a plot of the logarithm of the rate constant vs energy is given for calculations assuming alternatively a bond length and a center-of-mass separation distance reaction coordinate. The experimental results of Zewail and co-workers1*22 and Moore and c o - ~ o r k e r sfor ~ - the ~ overall and singlet dissociation rate constants (22) In all the comparisons with the experimental results of Zewail and co-workers (ref 1) presented here we consider only those based on the laserinduced fluorescence from the ground vibrational state of CH2 since there is some uncertainty as to the effect vibrational excitation might have on the measured dissociation rates.

Yu and Klippenstein

9884 The Journal of Physical Chemistry, Vol. 95, No. 24, 1991 12

I

I

I

I

I

A

A

11

+

A

0

+ 10

S

Q

+ +

it

Y

A

.

9

+h A

X X Q

'BOND'

0

'EXPT'

+

'CM'

0

' EXPT2

X

'INF'

A

a

I 0

I

I

1

I

I

1000

2000

3000 E (CM-1)

4000

5000

6000

Figure 1. Plot of log of the rate constant k vs excess energy E for the dissociation of CH2C0 into 'CHI and CO. The pluses and crosses denote the experimental results of refs 1 and 6, respectively. The diamonds and squares denote the RRKM calculationsassuming bond length and center-of-mass separation distance reaction coordinates,respectively. These calculations used the model potential energy surface given in ref 13. The triangles denote calculations using the number of states at infinite separation as the transition-state number of states.

of CH2C0, respectively, and the results of an NEJ(m)-based calculation are also presented in Figure 1. The difference between the two experimental results in the excess energy region from 450 to 2500 cm-I is simply an indication of the magnitude of the singlet/triplet branching ratios estimated by Moore and coworkers.6 The results presented in Figure 1 indicate that the bond length calculated rate constants are in reasonable agreement with experiment at the higher energies. However, a t the lower energies the bond length calculated results are considerably lower than the experimental results of Moore and co-workers.6 The fact that the bond length calculated NIJ becomes lower than the N E l ( m ) calculation by ==30cm-' also suggests that the calculations of the product state distributions presented in ref 14 would no longer be in agreement with the corresponding experimental results of Moore and co-workers! In particular, the bond length calculations reported here suggest that the change in slope observed in the vibrationally excited photofragment excitation spectra should occur at ==30 cm-' instead of at the experimentally observed energy of =IO0 cm-Ie4 The failure of the present better optimized variational RRKM theory calculations to provide agreement with experiment suggests the importance of reconsidering these previous calculations using an ab initio determined potential energy surface. Studies on the association of CH, and H to form CHI provide the only previous variational treatments of these "barrierless" associations employing similar completely a b initio based potential energy

surface^.^^-^^ 111. Methodology for ab Initio Calculations and Potential Energy Surface Fitting In the ab initio computations to be reported here a 6-31GS basis set was used and electron correlation effects were included through the application of second-order MaIIer-Plesset perturbation theory (23) Hase, W. L.; Mondro, S.L.; Duchovic, R. J.; Hirst, D. M. J . Am. Chem. Sor. 1987, 109, 2916. (24) LeBlanc, J. F.; Pacey, P.D. J . Chem. Phys. 1985,83, 431 1. (25) Aubanel, E. E.; Wardlaw, D. M.J . Phys. Chem. 1989, 93, 3117.

(MP2).26 These computations were carried out using the Gaussian 86 package of programs.n Recent calculations on the dissociation of CH4 have suggested that a more appropriate method for including correlation effects is through the use of configuration interaction calculations based on multiconfiguration SCF determined wave functions.28 One of the key evidences for the need to include multiple configurations in the S C F wave functions for the dissociation of CH4 was provided by the failure of a single configuration unrestricted Hartree-Fock (UHF) wave function to properly describe the spin state and energies of the system at extended bond length^.^*-^^ This failure remained even when fourth-order Mdler-Plesset perturbation theory (MP4) corrections were included. Interestingly, when the spin contaminant in the single configuration U H F wave function was projected out, the MP4 corrected energies were then found to be in reasonable agreement with the multiconfiguration S C F plus CI calculated energies.29 In the present case of C H 2 C 0 there is no evidence for any spin contamination in the U H F wave function, and the projected and unprojected energies are essentially identical for the range of bond lengths of interest here. As a result one may expect that an MP2 corrected single configuration wave function should provide an accurate estimate of the energetics, and the present calculations are based on this approach. A general difficulty that arises in the a b initio based study of the rate constants for these unimolecular dissociations with no reverse barrier is the need for a potential energy surface for all degrees of freedom, with a semiglobal description required for many of them. Here this requirement is reduced in two ways. (26) Szabo, A.; Ostlund, N.S . In Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory; McGraw-Hill: New York, 1982. (27) Frisch, M. J.; Binkley, J.

S.;Schlegel, H.B.; Raghavachari, K.; Melius, C. F.; Martin, R. L.; Stewart, J. J. P.; Bobrowicz, F. W.; Rohlfing, C. M.; Kahn, L. R.;Defrees, D. J.; Setger, R.; Whiteside, R. A,; Fox, D. J.; Fleuder, E. M.;Pople, J. A. Gaussion 86; Carnegie-Mellon Quantum Chemistry Publishing Unit: Pittsburgh, PA, 1984. (28) Brown, F. B.; Truhlar, D. G . Chem. Phys. Letr. 1985, 113, 441. (29) Schlegel, H. B. J . Chem. Phys. 1986. 84,4530.

The Journal of Physical Chemistry, Vol. 95, No. 24, 1991 9885

Dissociation of CHzCO into 'CH, and CO First, the previous calculations on a model potential energy surface are used as a guide to suggest for which reaction coordinate values to perform our calculations. These considerations lead to the development of a separate potential energy surface for each of the four bond lengths of 2.2, 2.5, 2.8, and 3.1 A. These latter values were chosen on the basis of the results presented in Table I, which suggest that this range should span the transition-state region at the energies of interest. Next, an assumption of a separability of modes into the "conserved" modes and the 'transitional" modes is made. The conserved modes correspond to the vibrational modes which change very little during the reaction process, Le., the modes which correlate with the CH2 and C O vibrations. The remaining modes, excluding the overall center-of-mass translation and the reaction coordinate, correlate with the free rotations of the products and of the line connecting them and correspond to the transitional modes. This assumed separability allows us to restrict our attention to only the transitional modes in the development of the potential energy surfaces. In the transition state the transitional modes are the overall free rotation and the relative rotation/vibration of the two fragments. These relative motions are best described in terms of coupled hindered rotors, and it is for only these modes that a detailed globally correct analytic potential energy surface is required. As a result, in general, at most five modes need to be considered for each transition-state potential energy surface with each surface corresponding to the interaction energy between two rigid fragments with their structures fixed a t the configurations corresponding to the reaction coordinate constrained minimum potential energy (or alternatively some other average configuration). Meanwhile, the conserved modes may be treated either in terms of an assumed correlation between reactants and products values or else in terms of a harmonic vibrational analysis for each reaction coordinate value. Here the latter approach has been applied, although the former would lead to essentially identical results. In the current case of the dissociation of CH$O there are four conserved modes-the two C H stretches, the HCH bend, and the C O stretch. The remaining modes are then the three overall translational modes, the reaction coordinate, and the seven transitional modes. Three of the latter seven modes correspond to overall rotation and four to relative motions of the two fragments. A four-dimensional potential energy surface is then required to describe the potential energy for these relative motions. The specific coordinates used in describing these transitional modes were the CCO bending angle Octo, the two HCC bending angles OHCC,, i = I , 2, and a torsion angle T describing the angle between the H H midpoint and the 0 atom viewed along the CC bond. In the development of the potential energy surfaces an optimized structure was first obtained at each of the four bond length values. A quadratic force field and corresponding vibrational frequencies were next obtained via numerical differentiation of energies obtained at nearby configurations. In the calculations of the frequencies at intermediate bond lengths internal coordinates force constants were used and the CC bond length contribution was projected out. Finally, the conserved mode coordinates were held fixed while electronic energies were calculated for four different sets totaling just over 200 different transitional mode configurations for each bond length. In these four sets alternatively one, two, three, or four of the transitional mode coordinates were varied from their optimized values over a grid of points. The majority of the calculations involved two and three coordinate variations. The ab initio determined electronic energies were then fit to trigonometric Fourier series expansions in the four transitional mode coordinates V~CCO, ~HCC,, 7 ) = Z [aic0 cos (kA0cco) + k,l,m.n

+

6fC0 sin (kAOcco)][arCClcos (IAO,cc,) b/"'l Sin ( ~ A ~ H C[Ut;fcc2 C ~ ) COS ] (mA8Hcc2) 6t;fcc1Sin ( ~ A ~ H c c , )COS ] [ ~(nAs) : 6: Sin AT)] ( 2 ) where A0 and AT denote the displacement from the respective bond length dependent optimized angles and the a and 6 coefficients

+

are parameters to be obtained in a least-squares fitting. This expansion was truncated by forcing the sum k 1 m n to be less than either 2 or 3 providing two analytical potentials with 28 and 8 2 free parameters, respectively, after the relevant symmetries were considered. A least-squares fit of the ab initio determined energies to these two Fourier series expansions was performed for each of the four bond length values. The most important region for determining the number of available states for a given bond length is of course the region near its optimized geometry. Thus, in the fitting procedure some additional constraints were implemented to force an accurate description of this region. In particular, the analytical potential was constrained to reproduce the quadratic force field for the transitional modes at the optimized geometry as determined in the ab initio calculations. This constraint reduced the number of free parameters to 17 and 71 for the two truncations of the k,l,m,n sum described above. Additionally, a weighting procedure described by Harding'O that decreases the weight for higher energy points by up to a given factor S was implemented with one minor change. In this procedure Harding introduced a weighting factor given by FT/, where

+ + +

(3) Ei is the ab initio determined energy of the ith configuration, and Eminand E,,, are the energy at the optimized geometry and the maximum sampled a b initio energy for the given bond length, respectively. The net result of the introduction of this weighting factor is the reduction in the weight of the higher energy points by a factor ranging from 1 to S. When this weighting procedure was implemented in the current calculations, an additional difficulty was encountered in that the higher energy points then obtained such a low weighting that the fitted energies were on occasion lower than that for the optimized geometry. Replacement of Ei in eq 3 with the minimum of the ab initio determined E, and the fitted Ei was found to remove this difficulty by forcing the high-energy points to at least correctly maintain a qualitatively correct high energy without forcing them to be accurately fit. The value of S used in the actual fits was then chosen by trial and error on the basis of a consideration of the accuracy of the fit in different energy ranges and the importance of these energy ranges in NE, evaluations. The results to be reported here were based on S values ranging from IO'O to loso. Typical accuracies for the root mean square errors in these fits at the energies of greatest importance were then on the order of 10% for the 17-parameter fits and 5% for the 71-parameter fits. Greater accuracies could of course be obtained through the inclusion of more parameters in the fitting procedure. However, a 71-parameter fit likely provides a reasonable upper limit to the number of parameters when fitting on the order of 200 ab initio points. Two additional potential energy surfaces with no free parameters are also considered in the statistical calculations in section IV. These latter two potential energy surfaces are simply the quadratic force field about the optimized geometry, V,d, and a hindered rotor expression that reproduces this quadratic force field, V,, (4) vhr

=

Vopt(R) 4-

ni/(1 - cos i

ABi)

+ ij,i#j G j / 2 ) sin ABi sin AOj ( 5 )

where ABi = Oi - 6,Oi ranges over 8 ~ ~O 0H C,C ~ ~and , T , and and VOptcorrespond to the bond length dependent optimized values for the angles and minimum potential energy, respectively. IV. Results and Discussion a. Structure Evaluations. The results from the geometry op timization along the reaction coordinate are summarized in Table 11. Energies and structures are reported therein for optimizations (30)Harding, L. 8. J . Phys. Chem. 1989, 93, 8004.

Yu and Klippenstein

9886 The Journal of Physical Chemistry, Vol. 95, No. 24, 1991 TABLE 11: ODtimized Structure as a Function of Bond Length’

R 1.318 1.35 1.40 1.50 1.70 1.90 2.2 2.5 2.8 3.1 m

1.316 m

~ C H

rco

~HCH

ecco

~HCC

E(m)-E(R)

1.081 1.080 1.079 1.077 1.093 1.103 1.108 1.108 1.108 1.108 1.109

1.180 1.179 1.178 1.176 1.164 1.153 1.148 1.148 1.148 1.149 1.150

Calculated 120.77 180.00 120.77 180.00 120.77 180.00 120.77 180.00 113.58 159.87 106.74 156.07 103.06 156.63 102.21 159.62 102.01 159.63 101.97 163.23 102.07

119.61 119.10 119.10 119.10 103.63 96.98 91.77 88.37 86.92 85.36

34653 34416 33227 28791 18470 10532 4844 2710 1683 1066 0.0

1.077 1.11

1.161 1.128

Experimentalb 122.0 180.00 102.0

119.0

32246 0.0

‘All distances are in units of angstroms, angles in degrees, and energies in cm-’. *Experimental values from refs 2 and 31-35.

of the reactants, products, the four bond lengths of relevance to the transition-state region, and some additional bond lengths intermediate between the reactants and the transition state. The symmetry of the optimized structures was C, for all intermediate bond lengths 1.50 8, or larger and C , for the others. The results reported in Table I differ considerably from those reported by Yamabe and Morokuma, as might be expected with the minimal basis set used there.19 However, the most important findings of no reverse barrier and a bent out of plane reaction path providing the lowest energy dissociation route were reproduced in the present calculations. The structural parameters of the reactants and products are in reasonable agreement with the corresponding experimental measurement^.^^-^^ It is interesting to consider how the bond length dependence of the optimized energies compares with some commonly used analytical representations such as the Morse, Varshni, Lippincott, and Rydberg f ~ n c t i o n s : ~ ~ . ~ ’ VM,,(R)

= D(1 - exp[-a(R -

A plot of the optimized energies and these analytical functions vs bond length is provided in Figure 2. The parameters D and Re for each of these potentials are provided by the calculated dissociation energy and equilibrium CC bond length of CH2C0. The single remaining parameter in each of these analytical functions is chosen to provide the calculated energy at a bond length of 1.5 A. In their study of the dissociation of CH4 Brown and Truhlar indicated that the Lippincott potential provides the best agreement with the calculated results in the intermediate range region and the Varshni potential provided the best overall agreement.2* A consideration of the results presented in Figure 2 indicates that none of the potential forms in eqs 6-9 provides a correct description in the 2-3-A region and the Lippincott potential provides the worst agreement. The failure of these (31) Moore, C. B.; Pimentel, G. C. J. Chem. Phys. 1963, 38, 2816. (32) Mallinson, P. D.; Nemes, L. J. Mol. Spectrosc. 1976, 59, 470. (33) Huber, K. P. Constants of Diatomic Molecules; Van Nostrand Reinhold: New York, 1979. (34) Feldman, D.; Meier, K.; Schmiedl, R.; Welge, K. H.Chem. Phys. Lett. 1978, 60, 30. (35) Petek, H.; Nesbitt, D. .I.Moore, ; C. B.; Birss, F. W.; Ramsay, D. A. J . Chem. Phys. 1987, 86, 1189. Petek, H.;Nesbitt, D. J.; Ogilby, P. R.; Moore, C. B. J. Phys. Chem. 1983.87, 5367. (36) Varshni, Y. P. Rev. Mod. Phys. 1957, 29, 664. (37) Stele, D.;Lippincott, E.R.; Vanderslice, J. T. Reu. Mod.Phys. 1962, 34, 239.

TABLE III: Frequencies as a Function of Bond Length conserved modes” asym C H sym C H Rb stretch stretch C O stretch H C H bend 1.318 3367 (3166) 3263 (3070) 2243 (2152) 1454 (1388) 1475 3007 2139 2.2 3103 2141 1492 3094 3002 2.5 2138 1499 2.8 3096 3006 2132 1500 3090 3006 3.1 m 3086 (2864) 3002 (2806) 2125 (2143) 1499 (1352) transitional modes’ Rb CHI rock C H 2 wag C C O bend C C O bend 443 (438) 588 (588) 528 (528) 1.318 1022 (977) 190 185 1001 812 2.2 589 136 117 2.5 79 1 2.8 58 1 40 1 106 85.8 82.3 67.6 415 272 3.1 reaction coordinate‘ C C stretch 1176 (1118) .

I

.

I

.

I

.

I

“All frequencies are in units of cm-I. Numbers in parentheses denote experimental fundamental frequency values from refs 31-35. b C C bond length in angstroms. TABLE IV: Force Constants as a Function of Bond Length

R“ 2.2 2.5 2.8 3.1

fllb

f12

h2

As

h4

2.1 0.94 0.54 0.35

1.5 0.87 0.49 0.28

8.0 4.7 2.4 1.2

1.2 1.0 0.72 0.42

0.29 0.13 0.052 0.026

fu 0.24 0.12 0.057 0.028

“ C C bond length in angstroms. bh,corresponds to the i, j force constant in units of cm-’/deg2. Mode 1 is the C C O bend, mode 2 is one HCC bend, mode 3 is the other HCC bend, and mode 4 is the OCCX torsion. All force constants not reported above are either related by symmetry to those reported or else zero.

analytical forms to provide an accurate description of the minimum energy dissociation path may be a result of the fact that a double bond is being broken and/or that the dissociation path involves a large measure of a bending coordinate. Alternatively, this failure may be an indication of an overestimate of the interaction strengths in this region in the a b initio calculations. In Tables 111 and IV the ab initio determined frequencies and transitional mode force constants are presented for the four bond length values of interest here. Also reported in Table 111 are the corresponding experimental anharmonic fundamental frequenc i e ~ . ~ IThe - ~ ~frequencies and force constants for the transitional modes roughly follow an exponential decay with decay constants ranging from 0.8 to 1.2 A-‘ for the frequencies and from 1.2 to 2.7 A-’ for the force constants. The off-diagonal force constants deviate the most from the average decay constant of 2.0 A-’ for the force constants. Meanwhile, the conserved mode frequencies which are expected to interpolate between their reactant and product values are already very close to the product values at a bond length of 2.2 A, suggesting a large decay constant of approximately 3 A-’ for these conserved mode frequencies. Only the HCH bend frequency at a bond length of 2.2 A is seen to differ significantly from the product value. This large decay constant is in agreement with the results of the previous theoretical studiesl3~l4 using a model potential energy surface which suggested that a decay constant of at least 1.8 A-1 is required for these conserved modes. The fluctuations arising in the calculated conserved mode frequencies are most likely a result of inaccuracies in the numerical differentiations and have no discernible effect on the present statistical calculations. b. Statistical Calculations. The analytic potential energy surfaces described in section 111 were implemented into the bond length based method for calculating NE, described in refs 17 and 18. In these calculations the bond length dependence of the conserved mode frequencies and geometrical parameters were

The Journal of Physical Chemistry, Vol. 95, No. 24, 1991 9881

Dissociation of CH2C0 into ICH2 and CO 0

-5000

-10000

'MEP' 'MORSE' 'VARSHNI ' 'LIPPIN' VRYD#

-15000 rl I

0

---

-....

-5 w

-20000

-25000

-30000

-35000

1

1.5

2

2.5

3

3.5

R (A)

Figure 2. Plot of minimum potential energy vs bond length. The diamonds denote the ab initio calculated energies. The short-, medium-, and long-dashed lines denote the Lippincott, Varshni, and Morse analytical potential fits, respectively. The short-long-short-dashedline denotes the Rydbcrg analytical potential fit.

TABLE V: Comparison of p,,Calculations Using Various Analytical Fits to the ab Initio Determined Energies

energy, cm-' 50

IO0 200 350 500 1000 2000 3500 5000 6000

Ni,,

71 param 1.7 (4), 3.06 2.0 (4), 3.03 2.8 (4), 2.98 4.2 (4), 2.90 5.8 (4), 2.83 1.4 ( 5 ) , 2.69 5.1 (5), 2.59 2.2 (6), 2.41 6.2 (6), 2.20 1.0 (7), 2.02

17 Param 1.6 (4), 3.07 1.9 (4), 3.06 2.6 (4), 3.02 3.9 (4), 2.96 5.5 (4), 2.90 1.4 (9,2.74 4.8 (5), 2.60 1.9 (6), 2.47 5.4 (6), 2.39 9.8 (6), 2.36

R" hr 1.2 (4), 3.07 1.4 (4), 3.07 1.9 (4), 3.06 2.8 (4), 3.05 4.0 (4), 3.02 1.1 (5), 2.90 4.1 (5), 2.64 1.6 (6), 2.57 5.2 (6), 2.49 1 . 1 (7), 2.42

Quad 1.0 (4), 1.2 (4), 1.7 (4). 2.6 (4), 3.8 (4). 9.8 (4). 3.2 (S), 1.4 (a), 4.6 (6), 9.2 (6),

3.13 3.13 3.10 3.06 3.00 2.81 2.64 2.56 2.48 2.42

frea 3.5 (3), 3.19 4.0 (3), 3.19 5.7 (3), 3.17 9.1 (3), 3.14 1.4 (4). 3.06 3.8 (4), 2.83 1.3 (5), 2.62 5.7 (5), 2.48 1.7 ( 6 ) . 2.37 3.2 (6),2.32

"Numbers in parentheses denote power of IO, and R* denotes the bond length value in angstroms at the minimum in NE,. obtained as an interpolation between the reactant and product values via

Xi = X,(Re)CaPt+

[Xi(R)"lc - Xi(Re)Ca'c] [Xi(m)exP' - x i(Re)CXP'l [Xi( ")talc - Xi(Re)Ca'c]

where i ranges over the four conserved modes, the calc and expt superscripts denote the present a b initio calculated and experimental parameters, respectively, and Re, m, and R denote the reactants, products, and variable bond length configurations, respectively. A J value of 3 was used in the present calculations for the reasons described in ref 13. Values of NiJ and R* were determined from the values of NEJ computed at the four bond lengths via a consideration of the minimum in a quadratic fit of the values for the three bond lengths nearest to the expected R*. The results of these calculations of N i , are presented in Table V. All values of R* reported in Table V are a result of this quadratic fitting procedure, and any values that are either greater than 3.1 A or less than 2.2 A should be viewed as lower and upper limits, respectively. The estimated Monte Carlo error bars are on the order of 5% for all the phase space integral based results reported here.

An additional calculation, labeled "freq" in Table V,corresponds to a standard tight transition-state variational RRKM calculation in which a rigid rotor harmonic oscillator description of the transition state is used.38 In this calculation we have combined a classical treatment of the transitional modes motions with a quantum treatment of the conserved modes. The transitional motions were treated as four harmonic vibrations and three overall rotations with the rotational constants corresponding to those for the bond length dependent optimized structure while the conserved modes were treated as harmonic vibrations. This mixed quantum classical treatment was used here to keep this treatment on the same footing as the other phase space integral based calculations. An analogous freq calculation in which the transitional mode vibrations were treated as quantum harmonic oscillators of the same frequency yielded results within 5% of these classical freq results. The transition-state frequencies used in these freq calculations are those reported in Table 111. A consideration of the results reported in Table V indicates that the 17-parameter and 71-parameter calculations are nearly identical, differing by at most 13%. This finding suggests that (38) Hasc, W. L.;Duchovic, R. J. J . Chem. Phys. 1W5,83, 3448.

9888 The Journal of Physical Chemistry, Vol. 95, No. 24, 1991

Yu and Klippenstein

0

11

0 0

+

0

+

Q

10

0

+

0

0 Y

+

+ x

*

++

t

+

2

x

X

'k-71' 0 'EXPT' + ' k-quad' 0 ' EXPTZ ' X

9

8

7

3000 4000 5000 6000 E (CM-1) Figure 3. Plot of log k vs E for the dissociation of CH2C0 into CH2 and CO. The pluses and crosses once again denote the experimental results of

0

1000

2000

refs 1 and 6 . The diamonds and squares denote the bond length based variational RRKM calculations using the 71-parameter and quadratic potential fits to the ab initio determined energies, respectively. a satisfactory fit to the potential has been obtained at the 17parameter level (at least for the purposes of the present statistical calculations). The hindered rotor and quadratic results are also in reasonably good agreement with the 7 I-parameter results, differing by at most 33 and 42%, respectively. Meanwhile, the freq results are as much as a factor of 5 below the 71-parameter results and a factor of 3 below the analogous quad calculations. The primary difference between the freq calculation and the quadratic expansion calculation is that the latter calculation correctly amounts for angular momentum couplings in both the kinetic energy and the angular momentum; Le., the calculations labeled "quad" are based on the phase space integral formalism described in refs 15-1 8, where the kinetic and angular momentum couplings within the different transitional motions and to the reaction coordinate are explicitly treated in a classically accurate manner. Meanwhile, the standard variational RRKM calculations (see, e&, ref 38) and the freq calculations here consider the angular momentum solely in terms of a rigid rotor. This treatment neglects a variety of possibly important effects such as the dependence of the rotational constants on the relative orientation of the two fragments and also the effect of the vibrational angular An additional difference between the two calculations is that the quad calculation correctly treats the restrictions on the range of the orientational coordinates. However, this last difference is expected to decrease rather than increase the quad calculated N i , relative to the freq calculation and thus cannot provide an explanation for the results found in Table IV. Additionally, simple calculations indicate that this effect is relatively minor in the present case. These findings suggest the importance of using the phase space integral based approach described in refs 15-18 instead of a simple rigid rotor treatment,' 5-18.39-41 (39) Wilson, E. B., Jr.; Decius, J. C.; Cross, P. C . Molecular Vibrations; Dover: New York, 1955. (40) Miller, W. H.; Handy, N. C.; Adams, J. E.J. Chem. Phys. 1980, 72,

Aubanel and Wardlaw have presented a similar comparison between the results of freq-like calculations with phase space integral based calculations for the dissociation of CHI into CH, and HSz5Interestingly, they found very little difference between the two calculations. There are a number of differences between the conditions in their calculations and the present calculations that might provide an explanation for this fact. In particular, one of the products is a single atom which may have a considerable effect on how the angular momentum couples. In addition, this single atom is a hydrogen atom for which the orbital angular momentum energies are anomalously large. Interestingly, hydrogen atom loss channels have previously been observed to have anomalously high kinetic energies.42 Finally, their calculations have been carried out for a canonical ensemble in which case the most important total angular momentum is generally quite large. The effect of the vibrational angular momentum on the conservation of J may then be relatively minor, leading to little difference between the calculated number of states. However, sample calculations carried out for the CHICO dissociation at a J value of 20 indicated roughly the same magnitude difference between the quad and freq calculations. This result indicates that the similarity between the two calculations in the case of CH, is not solely a result of the thermal averaging. It would be useful to consider these factors in greater detail. At this point it is interesting to consider a comparison with the experimental results of Zewail and co-workers.'qz2 and Moore and c o - w o r k e r ~ .Rate ~ ~ ~ constants were calculated according to eq 1 using a harmonic approximation for the reactant density of states and a classical approximation for the rotational contribution to this density. The transition-state number of states is assumed to be given by the absolute minimum of the inner transition-state number of states and the infinite separation number of states. A plot of the log of the rate constants calculated using the 71-parameter fit and the quadratic potential energy surfaces is given in Figure 3 along with the experimental data. Consideration of

99.

(41) Brocks, G.; Van der Avoird, A.; Sutcliffe, B. T.; Tennyson, J. T. Mol. Phys. 1983, 50, 1025.

(42) Chesnavich, W. J.; Bass, L.; Su, T.; Bowers, M.T. J . Chem. fhys. 1981,74,2228. Park, J.; Bersohn, R.; Oref, I. J . Chem. Phys. 1990,93,5700.

Dissociation of C H 2 C 0 into 'CH2 and C O the results plotted there indicates that the 71-parameter calculation, which is expected to be the most accurate of the present calculations, yields a result which is too high by a factor of about 2-3 for all excess energies at which the inner transition state is the dominant transition state (Le., above about 400 cm-I). At the very lowest energies the calculated result corresponds to the infinite separation number of states and is slightly too low at the lowest energy. This last result may be due to the strong dependence of the calculated rate constant on total angular momentum at these low energies. Some possible reasons for the discrepancy in the calculations at the higher energies are now considered. One likely explanation is simply that the difference arises as a result of the inaccuracies of the ab initio calculations; Le., an overestimate of the attractive interaction between the two departing fragments would push the transition state out to larger separations where the interfragment bending frequencies are relatively smaller, thereby leading to an increased number of states for the transition state. Another likely failure of the a b initio calculations is in overestimating the vibrational frequencies, a rough estimate of this factor for calculations at the present MP2/6-31G* level of theory being 5%.43 This fact has been explicitly taken into account for the conserved mode frequencies by correlating with the experimental fundamental vibrational frequencies as described in eq 10. However, for the transitional modes the present treatment involves the use of a highly anharmonic potential energy surface rather than simple force constants and thus a corresponding treatment could not be used. A simple classical harmonic oscillator estimate of this effect suggests that the rate constant should be increased by a factor of 1 .054= 1.22, resulting in a somewhat greater difference between theory and experiment. Another factor that has not been considered here is the effect of anharmonicities on the reactant density of states which may reduce the rate constant by a factor of approximately 1.5.4 Also, the reaction coordinate has not been fully optimized. Recent calculations for the dissociation of N C N O in which a three-dimensional optimization of Ni, is carried out have yielded a 10% decrease.4s The considerably larger difference between the results of the bond length and center-of-mass separation distance calculations for CHICO as compared to N C N O suggests that a corresponding optimization for C H 2 C 0 may lead to an even greater reduction. A combination of these effects may well explain the observed discrepancy, and thus there is no reason to invoke other explanations such as the possible presence of nonstatistical effects. These correction factors should also cause the transition from a dominant outer to a dominant inner transition state to occur (43) Bock, C. W.; Panchenko, Y. N.; Pupyshev, V. 1. J . Compur. Chem. 1990, 1 1 , 623. (44) See, e&., the calculations of the effect of the anharmonicities on the density of states for H2CO in: Barker, J. R. J . Phys. Chem. 1987, 91, 3849. Also in ref 5 it was indicated that a factor of 2 increase in the density of states was required to explain the observed triplet-state dissociation rates in terms of ab initio based RRKM theory calculations. (45) Klippenstein, S.J. J . Chem. Phys., in press.

The Journal of Physical Chemistry, Vol. 95, No. 24, 1991 9889 at about 100 cm-', suggesting that calculation^^^ based on the statistical theory of Marcus4 for the product state distributions might again be in satisfactory agreement with the experimental

result^.^

V. Concluding Remarks The dissociation of CHzCO presents an excellent test case for theoretical studies of unimolecular dissociations with detailed experimental results available for both the rate constants and the product state distributions. Here, a b initio electronic structure calculations have been used as the basis for a development of analytical potential energy surfaces in the transition-state region corresponding to bond lengths of 2.2-3.1 A. These analytical potential energy surfaces were subsequently utilized in statistical calculations of the rate constant for the singlet dissociation process. A bond length reaction coordinate was demonstrated to provide a more optimum reaction coordinate than the more commonly used center-of-mass separation distance and was thus implemented in the present ab initio based calculations. The specific approach used here should be generally applicable to the study of dissociation processes for which there is no reverse barrier. The results calculated on this basis were found to be a factor of 2-3 larger than the experimental results. A variety of correction factors were described whose inclusion may be expected to considerably reduce this discrepancy. Calculations based on a hindered rotor fit to the quadratic force field potential were in reasonably good agreement with the 7 l-parameter potential results, suggesting that a useful simple approach in future calculations might involve a determination of only the R-dependent minimum energy and force field for the transitional modes. Further ab initio calculations of analogous potential energy surfaces a t a higher level of theory would certainly be of interest. Of particular importance would be a consideration of the difference between multireference and single-reference calculations and of the effect of the basis set size on the strength of the interfragment interactions for bond lengths in the 2-3-A region. Theoretical calculations of the actual effect of the anharmonicities on the reactant density of states would also be of interest. Finally, these calculations have also demonstrated the importance of a correct treatment of the transitional mode angular momentum couplings as indicated by a difference of a factor of 3 between a rigid rotor harmonic oscillator treatment and the phase space integral treatment employing the equivalent quadratic potential energy surface and the classically correct kinetic energy and angular momentum terms.

Acknowledgment. Acknowledgment is made to the donors of the Petroleum Research Fund, administered by the American Chemical Society, for partial support of this research. Additional support was provided by the Ohio Board of Regents through a Research Initiation Grant and is gratefully acknowledged. Registry No. CH2C0, 463-51-4. (46) Marcus, R. A. Chem. Phys. Lett. 1988, 144, 208.