J. Phys. Chem. 1995,99, 15881-15889
15881
Further Studies of the Classical Dynamics of the Unimolecular Dissociation of RDX Candee C. Chamberst and Donald L. Thompson* Department of Chemistry, Oklahoma State University, Stillwater, Oklahoma 74078 Received: February 13, 1995; In Final Form: August 17, 1995@
Classical trajectories have been used to compute the unimolecular dissociation rates of RDX (hexahydro1,3,5-trinitro-1,3,5-triazine) for total energies over the range 200-450 kcavmol. The purpose of the study was to further refine the potential-energy surface. Thus, calculations were carried out on various surfaces, and results are reported for three of them. The three surfaces are similar except for the barrier height for the triple bond-fission reaction which gives 3CH2NN02 as products; rates are reported for barriers of 37, 56, and 71 kcaYmol for this channel. The energy required for the other primary dissociation channel, simple N-N bond fission, is about 48 k c d m o l for the three surfaces. The two surfaces with the higher barriers to ring fission yield branching ratios (ring fission to bond rupture) that appear to be inconsistent with the value determined from the infrared multiphoton dissociation experiments of Zhao et al. (J. Chem. Phys. 1988, 88,
801).
Introduction We have carried out several studies on the classical dynamics of RDX (hexahydro-1,3,5-trinitro-1,3,5-triazine):
These studies have focused on the unimolecular reaction dynamics,' conformational change^,^,^ and the basic behavior of the vibrational dynamic^.^-^ The challenge in modeling the chemical dynamics of RDX is the formulation of an accurate potential-energy surface (PES). Ab initio methods are of limited use in determining the global potential that includes the reaction channels of the size of the molecule. Thus, we have developed a PES that is based mainly on experimental data. Some ab initio results for the equilibrium reactant and product molecules are available in the literature and were used in adjusting the PES parameters. In an earlier study' of the unimolecular reaction dynamics of RDX, we computed the rates for the ring fission reaction that gives 3H2CNN02 molecules and the simple N-N bond rupture reaction giving the radicals NO2 and C3H&04, and the product energy distributions, over the energy range 250350 kcal/mol. These reaction channels were proposed by Zhao et al.' They used infrared multiphoton dissociation (IRMPD) to study RDX in a molecular beam. Their results indicate that the ratio of the rate of the molecular elimination reaction to that of the N-N bond rupture reaction is about 2 at a total energy (about 80 kcal/mol zero-point plus excitation) estimated to be in the range 150-170 kcal/mol. That is, they found that the two channels are competitive. In our earlier study,' we used potential-energy surfaces for which the barrier to N-N fission is 47.8 kcal/mol (the reaction
' Present address: Department of Chemistry, University of Minnesota, Minneapolis, MN 55455 Abstract published in Aduance ACS Abstracts, October 1, 1995. @
0022-3654/95/2099-15881$09.00/0
endothemicity) and barrier heights of 35.7,38.3, and 41.0 kcaV mol for the ring fission reaction (most of the calculations were done for the 38.3 kcaumol barrier). The calculated branching ratios are in the range 1.0-2.4. These results and the computed product energy distributions are in qualitative agreement with the molecular beam results of Zhao et al.' The calculations are for-energies well in excess of those in the experiments, and thus it is necessary to extrapolate the computed values to lower energies for the comparison. Because of the statistical error in the calculated values of the branching ratio, it is not clear how quantitative the agreement is. Nevertheless, the comparison with experiment suggests that the PES is qualitatively correct. However, some ab initio calculations have given values on the order of 70-80 kcdmol for the ring fission reaction.* This seems high in light of the Zhao et aL7 conclusion that this channel is competitive with the N-N bond fission reaction, which has an energy requirement of about 46-50 kcal/mol. In the present study we have carried out extensive classical trajectory calculations on similar but new PESs. The main reason for further study of the classical reaction dynamics of RDX is to help resolve the question of the barrier height for the ring fission reaction. That is, can a reasonable PES with a ring fission energy barrier on the order of 70 kcal/mol (as predicted by ab initio calculations) yield a branching ratio of 2:l for ring fission to N-N bond fission (which has a barrier of less than 50 kcal/mol)? Most of the calculations reported here were done for surfaces with ring fission barrier heights of 37, 56, and 71 kcaVmo1 (which we will denote as PES I, 11, and 111, respectively). The energy required for N-N bond rupture is 47.8 kcal/mol for all three PES. Trajectories were calculated over the energy range 200-450 kcal/mol. We also address problems in reconciling reported estimates of the RDX bond energies with the thermodynamics of the reaction to give the secondary products and the bond energies of the secondary product molecules. Also, we have performed a small number of calculations to investigate changes in the reaction rates if the RDX bond energies are adjusted to be in reasonable agreement with the thermodynamics of observing the secondary reaction products (as reported by Zhao et aL7). RDX has been the subject of considerable interest and some controversy lately. For example, there has been a problem of reconciling the results of experiments on condensed-phaseRDX with the Zhao et al.' results obtained in a molecular beam. The 0 1995 American Chemical Society
Chambers and Thompson
15882 J. Phys. Chem., Vol. 99, No. 43, 1995 latter experiments indicate that the triple bond fission reaction is competitive with the simple N-N bond rupture channel; however, the results for liquid and solid RDX do not indicate reaction via the concerted molecular elimination channel. This difference may be due to the importance of the volume of activation, as recently pointed out by Wight and Botcher? which may prohibit .the ring fission reaction in the confines of a liquid or solid. (It should be noted that there have been difficulties with regard to interpreting different experiments for condensedphase RDX reactions because of the observation of different sets of products.) As already discussed, another controversy, which we address in this study, concerns the energy barrier to the ring fission reaction. Some quite definitive experiments on the thermal decomposition of RDX (and HMX) have recently been reported by Behrens and Bulusu.'O They used combined thermogravimetric analysis and time-of-flight mass spectrometry to measure the products resulting from heated RDX samples. They determined the primary reaction pathways for the decomposition of RDX in the solid and liquid phases. They find that the initial reaction in solid RDX is the breaking of a N-N bond to form NO;?. Their data indicate that only about 10% of the decomposition could occur by the concerted molecular elimination mechanism, which has been observed in the gas phase;' however, their data are consistent with RDX H2CN NO2 2H2C=NN02 followed by the dissociation of the H2C=NN02 molecular fragments. The Wight and Botcher9 data show no evidence for the concerted molecular elimination of H;,C=NN02. The Behrens and Bulusu'O data show that following the elimination of N02, l-nitroso-3,5-dinitro-hexahydro-s-triazine (ONDNTA) is formed by reaction of the ring fragment with NO, which is postulated to come from either CH20 NO2 H20 CO NO or the decomposition of HONO which is formed by a concerted mechanism involving H-atom migration (from a C to an 0 atom) followed by N-N bond fission in RDX. The CH20 product disperses in the solid, leading to a weakening of the lattice structure and the formation of molten RDX, for which the decomposition reactions are much faster than in the solid. Recently, Wight and Botcher9 have reported a study of the photolysis (by a pulsed CO;, laser) of thin films of RDX deposited on an infrared-transparent window. They also indicate that the N-N simple bond fission reaction to give NO2 is the initial step in the thermal decomposition of RDX in the solid. Several other groups have reported relevant studies of RDX. Brill and co-workers'' have studied the surface region of energetic materials (RDX and HMX) heated by a T-jump method. They analyzed the product gases as functions of time by FTIR. Fifer and co-workersI2 pyrolyzed RDX and used chemiionization to produce ions which were then mass analyzed. A summary statement of the experimental results for the decomposition of RDX in the solid and molten phases is that the mechanisms are relatively complicated. However, these experiments have provided a relatively clear picture of the mechanisms. Adams and Shawl3 have recently reviewed this subject. There are now enough experimental data with which to make comparisons and thus provide guidance for the development of models for realistic simulations of RDX.
-
+
+
+
+
+
-
Computational Methods We have constructed three anharmonic potential-energy surfaces that include the two reaction pathways, NN bond rupture, and concerted ring fission. The surfaces were constructed by using available experimental and theoretical data. The three surfaces have barrier heights of 37,56, and 71 kcaV
mol for the ring fission reaction and 47.8 kcdmol for the N-N bond fission reaction. Trajectories were calculated on all three surfaces over the energy range 200-450 kcaymol, including zero-point energy (80 kcal/mol). Unimolecular dissociation rate coefficients were calculated for the two reactions. Our general classical trajectory code GenDynI4 was used to calculate the trajectories. Hamilton's equations of motion were integrated by using a fourth-order Runge-Kutta-Gill routine in a space-fixed Cartesian coordinate system. The integration step size was 1.0 x s and the trajectories were calculated for 50 ps or until the molecule dissociated. Initial conditions for the trajectories were selected by Monte Carlo quasiclassical and classical microcanonical methods; in the former, the energy is projected onto the normal modesI5 and in the latter classical initial conditions were generated by using a Metropolis procedure developed by Nordholm and co-workers.I6 The total angular momentum is zero in all of the calculations. The fiist-order rate coefficients, k, for the overall decay of RDX were computed by fitting the lifetimes to
(1)
ln(NJN,) = -kf
where Nt is the number of unreacted molecules at time t and NO is the number of trajectories in the ensemble. Since the molecule decays by two competing reactions: RDX
f
'
3H2C=NN02
C ~ H ~ N S NO2 O~+
k(ring) k(NN)
(R2)
the rate coefficients for the two channels are obtained from
k = k(ring) k(NN) = k/[1
+ k(NN)
+ N(ring)/N(NN)]
(3)
where k(ring) and k(NN) are the rates of reactions R1 and R2, respectively, and N(ring) and N(NN) are the numbers of trajectories which led to reactions R1 and R2, respectively. The branching ratio of concerted ring fission to N-N bond rupture is branching ratio = k(ring)/k(NN)
(4)
Potential-Energy Surface Thermochemistry of Reactants, Primary, and Secondary Products. The dissociation of RDX into its final products is believed to involve several reaction pathways. Zhao et al.' have shown that the primary pathways are simple NN bond rupture and concerted ring fission in the gas-phase. The products of the initial reactions further decompose. One of the products of the decomposition of CH2NNO2 is HCN HON0.'7.'s Determination of the values of the bond energies in RDX based upon primary products leads to dramatically different well depths' than does a similar determination based upon the secondary prod~cts.'','~-~'In this study, most of our calculations were done using PESs in which the well depths for the reactants and the primary products are based on the observed thermochemistry of the primary products. A few calculations were performed using much weaker bond strengths derived from the known thermochemistry of the secondary product^.'^^'^ The calculated branching ratios are approximately the same for the two different "parametrizations" of the PES. Sewell and Thompson' reviewed the thermochemistry of the reactant and products. The sum of the RDX bond energies calculated using the Sewell-Thompson PES is 1811 kcal/mol,
+
Unimolecular Dissociation of RDX RDX
3[CH$"NzI
J. Phys. Chem., Vol. 99, No. 43, 1995 15883 3[HCN+HONO]
Units are kcalhol
Figure 1. Range of possible thermochemistry (available in the literature) of the RDX molecule and one of its primary decomposition pathways (concerted ring fission) followed by further decomposition to secondary products. Maximum and minimum values for the endothermicity are obtained from different combinations of published values of the heat of formation of solid RDX, the heat of sublimation of RDX, the heat of formation of gaseous RDX, and the heat of formation of CH~NNOZ.The maximum and minimum exothermicities of the secondary reaction shown illustrate the range of values available in the literature. The asterisk marks the endothermicity suggested by Zhao et al.'
and the sum of the bond energies of the secondary products is 483 kcal/mo1.20.21The bond energies of the secondary products (three HCN HONO can be produced from one RDX molecule) gives product bond energies of 1450 kcdmol. This leaves a resulting reaction endothermicity from reactants to secondary products of 362 kcal/mol. The exothermicity of the reaction CHZNNOZ HCN HONO has been reported by Beyer and MorganI9 to be 10.3 kcallmol and by Mowrey et al.I7 to be 23.8 kcal/mol. These values, multiplied by three, give exothermicities of 30.9 and 71.4 kcallmol, respectively. These reported bond energies of RDX and the secondary products are not reconcilable with the reported thermochemical values. Figure 1 contains an illustration of the ranges of values of thermochemical quantities for the decomposition of RDX via the concerted ring fission pathway to primary and secondary products that are obtained using the various values given in the literature. We take the classical ground-state of RDX to be the zero of energy. The products 3CHzNN02 have been estimated to be as low as 24.4.' Since it has been suggested that the barrier height for concerted ring fission may be as high as 72-75 kcall mol,8 supposing that there is no barrier to the back reaction, we might assume an upper limit of 75 kcal/mol. Then, the secondary products are 30.9-71.4 kcaYmol(3 times the values reported by Beyer and MorganIg and Mowrey et al.,I7 respectively) lower in energy than the primary products. The maximum endothermicity of the reaction from the reactant (RDX) to the products (HCN HONO) is only 44.1 kcallmol and the maximum exothermicity is 47.0 kcal/mol; this is a broad range. The asterisk in Figure 1 marks the suggested endothermicity of Zhao et aL7 of 32 kcal/mol (see reaction VI1 in ref 7). If the range of exothermicities of the reaction of primary products to secondary products is correct, the value given by Zhao et aL7 for the entire reaction endothermicity is not possible unless the endothermicity of the primary reaction is at least 62.9 kcal/mol. To reconcile these differences, we adjusted the reactant and product (CH2NN02) bond strengths to be consistent with the available experimental and theoretical thermochemical data but not consistent with the available information on the bond
strengths. Calculations of the branching ratio gave values that are similar to those for the PES for which the bond strengths were derived from only thermochemical data for RDX and primary products as well as values suggested elsewhere in the literature.'qZ2 There are two points that should be noted. This comparison is for calculations performed on PESs with a barrier to concerted ring fission of 71 kcal/mol. The other is that while at this point we cannot conclusively say which set of bond strengths for RDX and CH2NNOz is correct, we have elected to use those suggested in the literature. We note that they may not be correct. However, the branching ratios calculated on the two surfaces are comparable, indicating that using the other values of the well depths would yield similar results. In other words, the branching ratio does not appear to be very sensitive to these parameters. Analytical Form of the Potential-Energy Surface. We have used potential-energy surfaces based on anharmonic stretches (Morse functions), harmonic bends, and truncated cosine series for torsions. We have previously shown with a series of s t u d i e ~that ~ ~such , ~ ~simple forms of potentials suffice to realistically describe the dynamics of various molecules. The analytical form is
bends
waes
+
-
+
+
It is the sum of Morse functions for the 21 bonds, harmonic oscillators for the 36 angles and 3 wags, and six-term cosine series for 12 torsional angles. We constructed three reactive surfaces by modeling the reactant molecule and the products based on information available in the literature. Geometries and frequencies of the products are available for only one set of the products, those due to concerted ring fission. Such information was not available for the fragment C3H&04, which results from simple NN bond rupture, thus we attenuated the values of the adjacent " 0 , CNN, and C2N-N angle force constants to zero as the NN bond distance becomes large; the other potential parameters (e.g., bond parameters) of the fragment were not attenuated. Information is available to construct a CHzNNOz force field.I7sz4 The reactant and product force fields were then used to construct a global potential by using switching functions (given in Table 1). The equilibrium geometries of the reactantsz5and productsI7 are shown in Table 2. All of the surfaces used in this study have the same values of these geometric parameters. The bond strengths in RDX and CHzNNOz are the same on all but one of the surfaces. (As discussed above, we calculated a single ensemble with bond strengths different from those used in the remainder of the study and found essentially the same results.) The bond energies used here (for PESs I, 11, and 111) are those that have been suggested elsewhere for RDX'~22~26 and CH2NN02'~22and are given in Table 3 along with the force constants. The frequencies which were obtained from a normal mode analysis (eq 5) and those obtained from experiment27and ab initio calculationsz8 for RDX are given in Table 4. The eigenvalues listed in Table 4 are not assigned to particular experimental frequencies, since the experimental normal mode assignments are not available. Yet, an attempt was made to fit the eigenvectors (and eigenvalues) to assignments suggested by the ab initio calculations.28 Table 5 compares the normal-mode
15884 J. Phys. Chem., Vol. 99, No. 43, 1995
Chambers and Thompson
TABLE 3: Potential Parameters
TABLE 1: Attenuation on the RDX PES coordinate
parameter based on switching function
Morse Parameters RDX De (kcaYmo1) a (,&-I)
bond
85.0" 47.8" 98.0b 95.0"
rCN TNN TNO rCH
CHzNzOz De (kcaumol) a (k')
1.954 422 2.651 632 1.951 689 1.941 979
153.07' 30.60" 98.0 108.0
1.734 45 1 2.675 695 2.576 902 1.861 508
Harmonic Force Constants
ke (kcaVmol rad?) angle OCNC ONCN eCNN ONNO
@ON0 ONCH OHCH
YCNC-N
RDX
CHzN202
189.593 228 145.334 935 231.652 616 241.113 503 205.484 71 1 58.274 924 56.180 756 95.001 720
357.1 18 978 319.754 009 245.089 04 63.025 996 52.472 672 0.0
Cosine Series Coefficients [kcaV(mol rad2)] dihedral angle CH2NN02 RDX" bonds (A) CN NN NO CH angles (deg) CNC NCN CNN NNO ON0 NCH HCH NCH-N CNCN, NCNC CNNO
1.464 1.413 1.213 1.089 123.7 109.4 116.3 117.25 125.5 110.571 577 105.1 20.091 558 36.340 12 19.1
this study 1.271 1.441 1.213 1.089
115.7 117.25 125.5 120.05 120.0
ab initiob 1.271 1.441 1.240, 1.303' 1.073, 1.07 1'
115.7 114.5, 121.1' 124.4 124.1, 116.0' 120.0
a Shishkov et aLZ5 Mowrey et al." Two bonds or angles reflect differences.between similar coordinates.
SCNCN ~CNNO
ao
al
0.425221 1.17345
0.0 0.0
a2 -0.430135 -1.65016
a3
a4
a5
0.0 0.0
0.361213 0.524958
0.0 0.0
a Melius and Binkley.22 Sumpter and Thompson.26 Sewell and Thompson.'
TABLE 4: RDX Normal-Mode Frequenciee calcb 36 (2) 39 46 49 (2) 104 (2) 138 280 (2) 302 328 (2) 360 508 525 (2) 588 (2) 598 641 (2) 688 698(2) 800(2) 803
expt'
ab initiod
calcb
expt'
ab initiod
22(2) 52
889(2) 93 1 972 1023 (2) 1018 1098 1137 (2) 1257 1251 (2) 1323 (2) 1339 1444 1477 (2) 1575 1606 (2) 2993 (2) 2993 3051 (2) 3054
910 885 935 1015 1045
995 (2) 912 1053 1070(2) 1117 1215 1121 (2) 1238 1309 (2) 1348 (2) 1274 1342 1364(2) 1481 1470 (2) 2954 (2) 2961 3066 (2) 3068
64
610 595
87(2) 198(2) 29 1 346 (2) 415 380 (2) 426 535 537(2) 606(2) 665 666 (2) 706 758(2) 912(2) 833
1230 1270 1320 1392 1435 1460 1550 1585 2980
frequencies calculated from the normal-mode analysis for CH2"02 and the values from an ab initio calc~lation.'~ The forms of the switching functions used to connect the reactant and product wells are the same for all the surfaces studied. The parameters of the switching functions were adjusted so as to give the desired barrier height to concerted ring fission. The switching functions are those given in eqs 18 and 19 of ref 1; they are
a Units are cm-'. The normal-mode frequencies were obtained by using second derivatives of the potential. The zero-point energy is 80 kcaYmo1. Iqbal et aL2' Coffin et al.** The ab initio frequencies have been scaled by a factor of 0.9.
Equations 6 and 7 are used when r > fl; when r < fl, Sl(r) and S2(r) are unity. (See ref 1 for a discussion of the calculation of the coefficients a(r,qR,qP).Table 6 contains the values of the parameters used to define the switching functions for the various surfaces. The rate at which a parameter such as a well depth, force constant, equilibrium bond length, or equilibrium angle changes during reaction is controlled by the shape of the switching functions that smoothly connect reactant and product parameters. We used 198 switching functions to incorporate the reaction
channels. Different sets of values of the switching function parameters were used to give barriers of 37, 56, and 71 kcaY mol (PES I, 11, and 111,respectively) for the ring fission reaction (see Table 6). Barriers to Reaction. Figure 2 shows plots of the minimumenergy paths (MEP) for concerted ring fission using the three different sets of values of the switching function parameters (see Table 6). We calculated the MEP as described in ref 1. The endothermicity of the reaction, 24.4 kcdmol, is the same on all three surfaces. The MEPs were calculated by extending the altemate CN bonds by 0.01 8, and then allowing the remaining coordinates to relax to the local minimum corresponding to the particular extension of the breaking CN bonds. Starting with the equilibrium geometry of RDX the altemate
740 794 750
3080
Unimolecular Dissociation of RDX
J. Phys. Chem., Vol. 99, No. 43, 1995 15885
TABLE 5: CHzNN02 Normal-Mode Freauencies ~
~
calc" (cm-')
~
~~
~~
ab initiob (cm-I)
calc" (cm-])
ab initiob (cm-I)
90.0 351.9 504.5 558.7 598.6 784.7 806.1 1098.4
1038.7 1140.5 1371.4 1479.9 1618.2 3022.3 3144.8
1135.1 1179.5 1422.0 1450.7 1655.1 3024.5 3134.3
C
370.1 C
553.0 589.7 779.7 780.5 928.3
Barrier: 37 kcaUmol
-1810
13
The normal-mode frequencies were obtained by using second derivatives of the potential. Mowrey et al." The ab initio frequencies have been scaled by a factor of 0.9.
-1740
2.0
25
1
P Sl(TCN) SWCN)
n
P n P
surface 1"
surface IIb
surface 111'
2 1.o 4 5.5 2 3.5
2 1.o 4 5.5 2 2.8
4 1.o 4 4.5
3.0
35
4.0
45
Barrier: 56 kcaUmol1
TABLE 6: Switching Function Parameters n
1
4 8.5
U
Concerted ring dissociation barrier height of 37 kcaumol. Concerted ring dissociation barrier height of 56 kcal/mol. Concerted ring dissociation barrier height of 71 k c d m o l .
2.0
25
3.0
35
4.0
45
a
CN bonds, bonds 1, 3, and 5 (if the ring bonds are numbered 1-6) are extended while bonds 2,4, and 6 are allowed to relax to their equilibrium values. The MEP is obtained by plotting the minimum energy as bonds 1, 3, and 5 are extended. There are several features of the MEP that should be noted. When the switching function parameters were varied to obtain the desired barrier heights, other features of the MEP changed as well. The initial regions of the MEP differ. The MEP with a barrier of 37 kcdmol (PES I) increases more gradually with increasing reaction coordinate distance than do the MEPs with the higher barriers. Another distinguishing feature is the position of the barrier. The MEP with heights of 37 and 56 kcdmol have barriers that occur at greater reaction-coordinate distances than the barrier of the surface with a barrier height of 71 kcaYmol. Data are currently not available to fit these features. Thus, we cannot judge their accuracy except by comparing the trajectory results with experiments. Note that with the same switching function parameters as used in the previous study,' but with a different analytical form of the PES, the profile of the low barrier, 37 kcdmol (see Figure 2) has a different shape at low values of the reaction coordinate from the previous surface (see ref 1, Figure 4). Comparison to the Sewell-Thompson Potential-Energy Surface. In the previous study of the unimolecular dissociation of RDX' a similar PES was used, although there are a number of important differences from those used here. The new PESs presented here were constructed with data that were not available when the previous surface was constructed. There are differences from the Sewell-Thompson surface' in the reactants, products, and regions between the reactants and products (due to the fact that the attenuation of the PESs differ from that used in ref 1). Additionally, the new PESs include more of the intemal motions; thus they are more realistic. The analytical form of the Sewell-Thompson' PES is also based on a valence force field; however, nonbonded interactions were used instead of cosine series. The PES in ref 1 is a sum of Morse functions for the 21 bonds, harmonic oscillators for the 36 angles and 3 wags, and 9 Lennard-Jones (6-12) potentials for the nonbonded interactions. The nonbonded interactions, not included in the present PESs (see eq l), are interactions between the C and N atoms across the ring, C and
-1740
c -1750
8
A
Barrier: 71 kcaUmol
-1770
O'-'"tJ -1810 I 15
I
2.0
25
3.0
35
4.0
45
Ructioo Coordinate (Anlptroarr)
Figure 2. Plots of energy profiles along the minimum-energy paths (MEPs) for ring fission in RDX for PES I, 11, and 111. The endothermicity of the reaction, 24.4 kcaymol, is the same on all three surfaces. The MEPs were calculated by extending the alternating CN bonds by 0.01 A and then allowing the remaining coordinates to relax to the local minimum. Starting with the equilibrium geometry of RDX we extended alternate bonds 1, 3, and 5 (where we number the ring bonds 1-6) while allowing bonds 2, 4,and 6 to relax. The MEP is obtained by plotting the minimum energies as bonds 1, 3, and 5 are extended. Branching ratios (knng:k") calculated from ensembles of trajectories at a series of energies are shown in Figure 4.
the N atoms of the -NO2 groups across the ring, and the N atoms of the -NO2 groups interacting with the N atoms of the other -NO:! groups. The PES used in this study includes 12 terms to represent the torsional forces in the ring ( ~ C N C Nand ZNCNC) and the forces for the -NO2 group twisting motions (ZCNNO). These terms were not used in the Sewell-Thompson surface. Different equilibrium geometries for RDX and CH2NN02 were used for the new PESs. The Sewell-Thompson PES' was based on the crystal-phase structure reported by Choi and Prince.29 A tetrahedral ring geometry was assumed as well as C2N-N wag angles of 0" with -NO2 groups which were neither equatorial nor axial. The new PES$ are based on the gas-phase geometry determined by electron diffraction by Shishkov et al.25 Two of the most significant differences are that the nitro groups are in axial positions and that the C2N-N wag angle is 19.9" for the new PESs. This difference in the positioning and orientation of the nitro groups probably affects the dynamics.30 The only vibrational frequencies available at the time of the construction of the Sewell-Thompson surface' were those from the experimental work by Iqbal et aL2' There are 13 normal-
'
Chambers and Thompson
15886 J. Phys. Chem., Vol. 99, No. 43, 1995
mode frequencies’ for which experimental values are still not available. For the PESs used in the present studies, we fit the eigenvalues to the frequencies reported by Iqbal et al.27and to the ab initio frequencies (which we scaled by a factor of 0.9) reported by Coffin et a1.28 Additionally, we attempted to fit the eigenvectors using the reported contributions of intemal coordinates to the normal modes from the ab initio calculations.28 Comparisons of the eigenvectors for the new PESs with those calculated from the Sewell-Thompson PES reveal that they are quite different except for the high-frequency CH stretch modes. The attenuations of the potential parameters as the system goes from reactants to products are similar for the two surfaces. The forms of the switching functions (see eqs 6 and 7) are identical. For PES I, which has the lowest barrier (37 kcaY mol) to concerted ring fission, even the values of the switching function parameters are the same (although some of the asymptotic limits differ). The only difference is that four parameters (UNO, koNo, ~ H C H ,and YCNC-N) are attenuated on PES I which are not attenuated on the Sewell-Thompson PES.3’ Finally, there are some differences in the CH2NN02 product molecule portions of the surfaces as well. The differences are similar to those described for the reactant portion of the surfaces. Well depths for the bonds are the same in all the PESs. The geometry used for them was that calculated by Mowrey et aLi7 The Sewell-Thompson PES’ used CH2NN02 geometries and frequencies from the calculations at the stationary points for NN bond rupture (Table 111, ref 17). Whereas the current PES uses CH2NN02 geometries and frequencies from the calculations at the stationary points for concerted dissociation (Table I, ref 17). These variations in the construction of the PESs lead to different values of the force constants.
TABLE 7: Reaction Rates and Branching Ratio@ energy (kcaVmol) 200
225
250
300
350
400
450
RMSf
IIId
Ib
&ring) k(NN) branching ratio k(ring) k(NN) branching ratio k(ring) k(NN) branching ratio &ring) k(NN) branching ratio k(ring) k(NN) branching ratio &ring) k(NN) branching ratio k(ring) k(NN) branching ratio k(ring) k(NN)
I‘e
0.0029 0.0031 0.93 (112) 0.0086 0.018 0.47 (159) 0.0 0.035 0.055 0.051 0.69 (1 16) O.O(l9)
0.0 0.24 O.O(l21)
0.021 0.032 0.66 (100)
0.060 0.0052 0.00069 0.082 0.19 0.11 0.039 0.23 0.36 (300) 0.028 (339) 0.0063 (320) 1.5 (300) 0.011 0.0072 0.11 0.33 0.3 1 0.47 0.25 (200) 0.033 (186) 0.023 (200)
0.16 0.072 2.3 (200)
0.022 0.023 0.21 0.56 0.65 0.67 0.31 (150) 0.039 (325) 0.035 (300) 0.044 0.78 0.056 (300)
0.40 0.91 0.44 (300)
0.02 0.16
0.00 0.02
0.01 0.13
Rates are in ps-’ and branching ratios are k(ring):k(NN). Numbers in parentheses are ensemble sizes. Surface with a barrier to concerted ring fission of 37 kcal/mol (PES I). Surface with a barrier to concerted ring fission of 56 kcal/mol (PES 11). Surface with a barrier to concerted ring fission of 71 kcaUmol (PES 111). e PES I. A ring mode was initially excited. f See text for discussion.
Results and Discussion The rate coefficients and branching ratios were calculated over the energy ranges 200-450 (PES I), 250-400 (PES II), and 250-450 kcaYmol (PES 111), including about 80 kcaYmol of zero-point energy. The results are summarized in Table 7. The number of trajectories in each ensemble is given in parantheses in Table 7. The rate coefficients were obtained by fitting the trajectory results to eq 1. Plots of ln(N,/iV,,) vs time show that the decomposition occurs by single first-order decay, that is, the decay curves are linear. The root-mean-square deviations of the values of the rate coefficients calculated from eqs 1-3 from the values of rate coefficients calculated from the fits to eq 8 are given in Table 7. We did not observe the ring fission reaction at 250 kcaVmol for ensembles of 50 trajectories integrated for 50 ps on PESs LI and 111; therefore, we did not calculate trajectories at lower energies on these PESs. We calculated an ensemble of 29 trajectories at 160 kcdmol (the estimated energy in the Zhao et ale7experiments) on PES I, but found no reaction in either channel in the 50 ps trajectory calculation time. On the basis of extrapolation of the higher energy results to this energy, this is not surprising, since the decay time of the concerted ring fission would be on the order of (2-3) x lo3 ps. Figure 3 contains RRK plots for the two reaction channels, NN bond rupture (panels a, c, and e), and concerted ring fission (panels b, d, and f). The data were fit to ln[k(E)I = ln[Al
+ (s - 1) ln[l - (EdE)]
(8)
where EOwas taken to be 47.8 kcdmol for the NN bond fission reaction for all three fits and & was set equal to 37,56, and 71 kcaVmol for the concerted ring reaction channel for PES I, 11,
and 111, respectively. Fits of the data to eq 8, give values of the frequency factor, A, of 2220, 597, and 112 cm-’ for the NN reaction on PESs I, 11, and 111, respectively. Values of the frequency factor of 373, 34.3, and 2760 cm-’ were calculated from fits to the rates of the concerted ring fission reaction as a function of energy for surfaces I, 11, and 111, respectively. The calculated effective degrees of freedom, s, are 36, 28, and 16 for the NN reaction channels and 40, 27, and 44 for the concerted ring fission channels for surfaces I, LI, and 111, respectively. In both cases, the computed s values are less than the theoretical value of 57. The values for the concerted reaction are consistently larger than those for NN bond rupture. This is as expected since the former reaction is more “global” and is expected to involve more molecular modes than the latter, that is, there are fewer “inactive modes” in the latter. The calculated values of the effective degrees of freedom from the SewellThompson surface’ (43.4 and 38.4 for R2 and R1, respectively) are similar to those obtained here for the PES with a similar barrier height to concerted ring fission (PES I). The values of the frequency factor, on the other hand, differ significantly. Values of the frequency factor from the Sewell-Thompson surface were 9673 and 1433 cm-’ for reactions R2 and R1, respectively. Figure 4 contains plots of the branching ratios as a function of energy for the three PESs (the experimentally determined’ branching ratio at approximately 160 kcaVmol is shown as a solid triangle). One can see that PESs I1 and 111, with barriers to concerted ring fission of 56 and 71 kcaYmo1, respectively, give values of the branching ratio that are quite small and, for the ranges of energies studied, are essentially independent of the energy. We did not observe the ring reaction at 250 kcall mol on either of these surfaces. It does not appear, based on
Unimolecular Dissociation of RDX
J. Phys. Chem., Vol. 99, No. 43, 1995 15887
A
TABLE 8 Comparison of Actual Ensemble Sizes and the Initial Populations Predicted by the Intercepts of the Decay Curves for PES-I initial population total energy (kcaymol) actual predicted, eq 1 200 112 1 13.6 225 159 158.9 250 116 123.3 300 300 227.1 350 200 173.3 400 150 118.7 450 300 250.3
*
? A
0 0 PES1 PESII A PES 111
+ -6 -0.3
1
0 I
-0.1
-0.2
I n[l -Eo/E] 0
-2
1
0 0
-
0
0
-4 -
+
A
A
0
0 PES1 PESII A PES 111
+
1 ’
-8 -I -0.30
I
I
-0.20 -0.10 In[l-Eo/E] Figure 3. RRK plots calculated for the reaction R2 (panel a) and reaction R1 (panel b) for the PESs I (circles), I1 (crosses), and I11 (triangles).
r,l I
0 Barrier:37kcaVmol
0
m
A
0
Barrier: 56 kcaVmol
A A
Barrier: 71 kcaVmol Experiment
I
S
. I
“ 1
O l 100
0
O -
m
. p
200
0 a
m I
I
300
400
. II
A
500
Energy (kcal/mol) Figure 4. Calculated branching ratio (kin&”) as a function of energy for PES I, 11, and 111. Note that the ring reaction on PES I1 and I11 goes to zero at 250 kcal/mol. These results, along with the kingand k” values, are listed in Table 7. these results, that they would give a branching ratio of 2:l at 160 kcal/mol. The computed branching ratios for PES I are shown as solid circles in Figure 4. At energies above 300 kcal/mol the ratio is essentially independent of energy. As the energy decreases from 300 to 250 kcdmol the ratio increases (the statistical error causes some nonmonotonic behavior) such that a reasonable extrapolation to 160 kcal/mol would give a value in accord with the experimental value (2), shown as a solid triangle in Figure 4.
We can use RRK fits of the calculated rates at high energy to make extrapolations to the energy range (-160 kcal/mol) of the IRMPD experiment^.^ The accuracy of the prediction of the rate at the lower energy depends on the validity of the statistical assumption of the RRK theory, that is, that dynamical bottlenecks are not present over the range of the fit and extrapolation. The simplest way to determine this is by comparing the actual initial populations (i.e., trajectory ensemble sizes) to the “predictions” of the linear portion of the decay curves (eq 1). If the intercept of the decay curves equals the actual ensemble size, there are no significant bottlenecks and thus the RRK extrapolation to lower energy should be valid. Table 8 shows this comparison for PES-I for the energies for which ensembles of trajectories were computed. At the lower energies, there is good agreement between the intercept of the semilog plot of the calculated populations vs time and the actual initial populations of the ensembles; however, there are some differences at the higher energies. If we use only the first two points (200 and 225 kcal/mol) of the RRK plots to predict the branching ratio at 160 kcal/mol, we obtain a branching ratio of 5:l (ring fission to bond fission). If we use all seven values, we obtain a branching ration of 1.4: 1. Making predictions using all the combinations of the points, excluding those for 400 and 450 kcal/mol, give values of the branching between 1.4 and 1.8. Thus, our calculations for PES-I predicts the branching ratio to be between 1.4 and 5. Since the decay time of RDX at 160 kcal/mol (the energy in the experiments) is longer than the times for which we can integrate trajectories, it is difficult to obtain reaction rates at that energy. Despite this, we attempted to observe reaction at this energy by calculating a small (29 trajectories) ensemble for 50 ps; we observed no reactions. On the basis of the RRK fits, we would expect the rates of ring and NN reaction at 160 kcaVmol (on surface I) to be in the ranges (0.5-4) x and (1-3) x ps-I, respectively. Thus, it is not surprising that we did not observe the concerted ring fission reaction. Although we were unable to perform calculations of the branching ratio at the energy of the experiment,’ reasonable extrapolations of the computed results for PES I to that energy suggest that the potential with a 37 kcal/mol barrier to ring fission is in accord with the measurement. On the other hand, the other two surfaces give results that indicate that the barriers significantly higher than that of PES I are not. The results of all of our calculations, including those published earlier,’ indicate that the height of the barrier to ring fission is lower than that for NN bond fission. We did a few calculations for initial conditions with modespecific excitations. The rates and branching ratios resulting from calculations on PES I with the ring mode which has a frequency of 972 cm-I (see Table 4) initially excited are given in the last column of Table 7. The branching ratios as a function of energy for both the mode-specific and random initial condition ensembles are shown in Figure 5. At higher energies,
15888 J. Phys. Chem., Vol. 99, No. 43, 1995
0
Random
A
Mode Specific
Chambers and Thompson
A
a
a 2
m J
A
C
I
E 1
V
’
100
200
300
400
500
Energy (kcal/mol) Figure 5. Calculated branching ratio as a function of energy for PES I (which has a barrier to concerted ring fission of 37 kcaymol), with classical microcanonical (circles) and mode-selective (triangles) initial conditions. In the mode-selective calculations, the molecule was assigned zero-point energy, and then the excitation energy was placed in the ring mode with a frequency of 972 cm-’ (see Table 4).
initial excitation of the ring mode enhances rather dramatically the knng:k” branching ratio. The energy dependencies of the ring fission, rate for the randomly and nonrandomly distributed energy initial conditions (see Table 7) are roughly the same, and the values are comparable. However, there is a large difference in the energy dependencies of the two sets of rates for the NN bond fission reaction. It is interesting that for the mode-selective initial conditions that the rate of NN bond fission is essentially unaffected by increasing the energy. Thus, k(ring) increases with energy much as it does for nonselective initial conditions, while k(NN) remains almost constant. This results in the rapid increase in the branching ratio with energy (see Figure 5).
Conclusions We have presented the results of classical trajectory calculations of the unimolecular decomposition of RDX. We have computed the rates for the two reactions, the molecular elimination of 3H2CNN02 due to triple ring-bond rupture and the simple N-N bond rupture giving the radicals NO2 and CzhN404, proposed by Zhao et aL7 to explain the results of their infrared multiphoton dissociation (IRMPD) molecular beam experiments. Zhao et al. reported that the ratio of the triple bond-rupture reaction to simple N-N bond rupture is about 2 at a total energy estimated to be around 160 kcdmol. It is this branching ratio that we focused on in this study. The challenge is to develop a PES that accurately represents the forces in RDX by using simple analytical forms with the values of the parameters determined from the sparse information (e.g., spectroscopic, thermochemical, and kinetic data, and some ab initio results) that is available in the literature. In a previous study’ we developed such a potential and carried out classical trajectory calculations of the rates and product energy distributions. Comparison of the calculated results with the experimental data7 indicated that the PES is qualitatively correct. The present study was performed because new information became available that could be used to improve the PES. Also, it has been suggested that the barrier to the ring fission reaction may be twice as large as that proposed by us in our first study. The barrier height to N-N bond rupture in RDX seems to be well-established as 47.8 kcaVmol. However, the barrier height for the ring fission reaction is not as well-known. We concluded in our earlier study‘ that it is on the order of 40 kcaV
mol or less. However, some recent ab initio results,8 as well as those of Zhao et al.,7 suggest values ranging from 72 to 80 kcdmol. Thus, we have carried out further calculations of RDX in order to determine what the value of the ring fission barrier must be in order to calculate a value of the branching ratio that is in accord with experiment. We carried out calculations on three PESs which are quite similar except that they have barriers of 37,56, and 71 kcaVmol for ring fission reaction. The results of these calculations support the conclusions of our first study.’ That is, for the various PESs of the type we have considered, those with the higher barriers (40 kcdmol) for the ring fission reaction yield results that appear to be inconsistent with the experimentally determined branching ratio. Our studies indicate that the barrier height is on the order of 37 kcaVmo1 and that the barrier heights that are significantly greater than this do not give results in accord with the e~periment.~ A more general conclusion that might be drawn from this study is that the reaction dynamics of large molecules can be simulated with relatively simple PESs. Even our earlier study,’ that used a PES that describes the molecular vibrations less well than does the present one gave results in qualitative agreement with experiment. In both cases, the PESs were parametrized by using the limited information that is available in the literature, which is less than that available for many other systems.
Acknowledgment. This work was supported by the U.S. Army Research Office (Grants DAAH04-93-G-0450 and DAAL03-92-G-0358). We also thank the Oklahoma State University Computer Center for providing some of the computer time. C.C.C. would also like to thank Dr. Thomas D. Sewell for many helpful discussions. References and Notes (1) Sewell, T. D.; Thompson, D. L. J. Phys. Chem. 1991, 95, 6228. (2) Wallis, E. P.; Thompson, D. L. Chem. Phys. Lett. 1992, 189, 363. (3) Wallis, E. P.; Thompson, D. L. J. Chem. Phys. 1993, 99, 2661. (4) Sewell, T. D.; Chambers, C. C.; Thompson, D. L.; Levine, R. D. Chem. Phys. Lett. 1993, 208, 125. ( 5 ) Chambers, C. C.; Sewell, T. D.; Thompson, D. L., manuscript in preparation. (6) Chambers, C. C. Ph.D. Thesis, Oklahoma State University, 1994. (7) Zhao, X.; Hintsa, E. J.; Lee, Y. T. J . Chem. Phys. 1988, 88, 801. (8) Habibollahzadeh, D.; Grodzicki, M.; Seminario, J. M.; Politzer, P. J. Phys. Chem. 1991, 95, 7699. (9) Wight, C. A.; Botcher, T. R. J . Am. Chem. SOC. 1992, 114, 8303; J . Phys. Chem. 1993, 97, 9149. (10) Behrens, R., Jr.; Bulusu, S. J. Phys. Chem. 1992, 96, 8877; 1992, 96, 8891. (1 1) Oyumi, Y.; Brill, T. B.; Rheingold, A. L.; Haller, T. M. J . Phys. Chem. 1985,89,4317. Mesaros, D. V.; Oyumi, Y.; Brill, T. B.; Dybowski, C. J . Phys. Chem. 1986, 90, 1970. Oyumi, Y.; Brill, T. B. J . Phys. Chem. 1987, 91, 3657. (12) Snyder, A. P.; Kremer, J. H.; Liebman, S. A,; Schroeder, M. A.; Fifer, R. A. Org. Muss. Spectrom. 1984, 24, 15. Synder, A. P.; Liebman, S. A.; Schroeder, M. A.; Fifer, R. A. fbid. 1990,25, 15. (13) Adams, G. F.; Shaw, R. W., Jr. Ann. Rev. Phys. Chem. 1992,43, 311. (14) GenDyn is a general classical dynamics program developed at Oklahoma State University; for a description of the basic components of the code see: Bintz, K. L. M.S. Thesis, Oklahoma State University, 1986. (15) Raff, L. M.; Thompson, D. L. In: Theory of Chemicul Reaction Dynamics Vol. 3, Baer, M., Ed.;CRC Press: Boca Raton, FL, 1985; Vol. 3, p 1. (16) Schranz, H. W.; Nordholm, S.; Nyman, G. J. Chem. Phys. 1991, 94, 1487. Nyman, G.; Nordholm, S.; Schranz, H. W. J . Chem. Phys. 1990, 93, 6767. (17) Mowrey, R. C.; Page, M.; Adams, G. F.; Lengsfield, B. H., I11 J . Chem. Phys. 1990, 93, 1857. (18) Rice, B. M.; Adams, G. F.; Page, M.; Thompson, D. L. J . Phys. Chem., submitted. (19) Beyer, R. A.; Morgan, C. U. Sixteenth JANNAF Combustion Meeting Proceedings 51, 1979. (20) Guan, Y.; Lynch, G. C.; Thompson, D. L. J . Chem. Phys. 1987, 87, 6957. Guan, Y.; Thompson, D. L. Chem. Phys. 1989, 139, 147.
Unimolecular Dissociation of RDX (21) Murrell, J. N.; Carter, S.;Halonen, L. 0. J. Mol. Spectrosc. 1982, 93, 307.
(22) Melius, C. F.; Binkley, J. S.Twenfy-$rst Symposium (International) on Combustion 1953, 1986. (23) (a) Sumpter, B. G.; Thompson, D. L. J . Chem. Phys. 1985, 82, 4557; 1985, 86, 2805; 1987, 86, 3301; 1987, 87, 5809. (b) Bintz, K. L.; Thompson, D. L. J . Chem. Phys. 1986, 85, 1848; Chem. Phys. Lett. 1986, 131, 398; J. Chem. Phys. 1987.86.4411. (c) Guan, Y.;Thompson, D. L. J . Chem. Phys. 1988, 88, 2355; Chem. Phys. 1990.92, 313. (d) Uzer, T.; MacDonald, B. D.; Guan, Y.; Thompson, D. L. Chem. Phys. Lett. 1988, 152,405. Ahances in Molecular Vibrations and Collision Dynamics; Vol. lB, Bowman, J. M., Ratner, M., Eds.; JAI Press: Greenwich, CT, 1991; p 81. (e) Gai, H.; Thompson, D. L.; Fisk, G. A. J . Chem. Phys. 1989, 90, 7055. (0 Gai, H.; Thompson, D. L. Chem. Phys. Lett. 1990,168, 119. (8) Schranz, H. W.; Raff, L. M.; Thompson, D. L. J. Chem. Phys. 1991, 94, 106. (h) Qin, Y . ; Thompson, D. L. J. Phys. Chem. 1992, 96, 1992. (i) Chang, X.; Thompson, D. L.; Raff, L. M. Chem. Phys. Lett. 1993, 206, 137.
J. Phys. Chem., Vol. 99, No. 43, 1995 15889 (24) Rice, B. M.; Chabalowski, C. F.; Adams, G. F.; Mowrey, R. C.; Page, M. Chem. Phys. Lett. 1991, 184, 335. (25) Shishkov. I. F.: Vilkov. L. V.: Kolonits, M.: Rozsondai. B. Srruct. Chem.'1990, 2 , 57. (26) SumDter. B. G.: Thomoson. D. L. J . Chem. Phvs. 1987.86. 3301. (27) Iqbal, Z.; Suryananraianan, K.; Bulusu, S.; Autera, J. R.; Army Technical Report 4401, 1972. (28) Coffin, J. M.; Newton, S.Q.; Ewbank, J. D.; Schafer, L.; Alsenoy, C. V.; Siam, K. J . Mol. Struct. 1991, 251, 219. (29) Choi, C. S.; Prince, E. Acta Ctystallogr. 1972, 828, 2857. (30) Burov, Yu. M.; Nazin, G. M. Kinet. Katal. 1982, 23, 12. Maksimov. Yu. Ya.: Apal'kova. V. N.: Braveman, 0. V.; Solov'ev, A. I. Zh. Fiz. Khim. 1985, f9,342. (31) There is a twJograDhica1 error in Table 7 of ref 1. The functional dependence of the force constant for wag-angle displacements is given incorrectly. The correct fonn, which was used in the calculations, is ky' = I
,
1
__
~
~~OSII(~NN)SII(~'CN)SII(~'+'CN). JP950395M
