(Paper 04S-25). Multi-Channel Chemically Activated Reactions: Comparison of Troe’s Modified Strong Collision Model and Exact Solution of the Master Equation by Monte Carlo Method* Ameya V. Joshi,1 Scott G. Davis,2 and Hai Wang1 1 Department of Mechanical Engineering, University of Delaware, Newark, DE 19716, 2 Exponent, Inc., 21 Strathmore Road, Natick, MA 01760
Introduction The temperature and pressure dependence of chemically reacting systems is essential in the modeling of combustion phenomena. Apart from the direct application to simulation of these complex systems, a fundamental understanding of the dependence of the rate coefficients and product yields on pressure and temperature is important for extrapolation of the available experimental data and development of detailed reaction mechanisms. At the high temperatures prevalent in combustion systems, many gas phase reactions proceed through the formation of rovibrationally excited adducts. These reactions are called as chemically activated reactions. The fate of the excited species is determined by several competing processes. Energy is transferred to and from the molecule as a consequence of inelastic collisions with the surrounding bath gas molecules. Between collisions, the molecule can undergo internal rearrangement (isomerization), decompose to form new products or simply dissociate back to the reactants. Much progress has been made in our understanding of collisional energy transfer and its influence on the temperature and pressure dependence of unimolecular and chemically activated reactions. Recognizing the inadequacies of the ‘strong collision’ assumption, Troe1 introduced a collision efficiency factor β c , such that the rate constant at low pressure was related to the strong collision rate constant through k0 ≡ β c k0sc (1) The modified strong collision approach thus allowed for weak collisions, with the efficiency of energy transfer dependent on the nature of the bath gas. While this constitutes a major advance in the treatment of energy transfer, a major drawback of the modified strong collision model is that it provides an oversimplified picture of weak collisions, by treating energization/deenergization as essentially single-step processes. As mentioned before, the adduct formed in a chemically activated reaction often has sufficient internal energy to undergo further reactions, and it takes several collisions to remove that excess energy before it becomes trapped in a potential well. Though the modified strong collision model was intended for use only for single-channel
Presented at the Spring Technical Meeting of the Western State Section Meeting of the Combustion Institute, University of California, Davis, March 29-30, 2004.
unimolecular reactions,1 over the years this model has been widely used for multi-channel chemically activated reactions. The fundamental validity of this approach is certainly questionable. However, the ease with which this model can be applied to chemically activated reactions makes it an attractive approach of rate constant estimation and extrapolation. To this end the accuracy of the modified strong collision treatment of multi-channel chemically activated reactions remains unknown. The purpose of the present study is to compare Troe’s modified strong collision model and the exact solution of master equation for multi-channel chemically activated reactions. In the following section, we describe the details of a program developed to calculate rate coefficients for thermally and chemically activated systems by solving the Master Equation using the Monte Carlo approach. Such a program is certainly not the first of its kind and there exist other codes2-8 with varying approaches and capabilities developed to solve the Master Equation. The thrust of this work is a comparison of the solutions obtained using the Monte Carlo program with those obtained with the modified strong collision model for a couple of chemically activated systems. The Master Equation Rovibrationally excited molecules of a reacting system are subjected to competitive processes of collisional activation/deactivation and unimolecular isomerization and dissociation. The time evolution of such a system is then adequately described by the master equation (see Ref. 9,10 for a detailed review). In the discrete form the master equation is given by11-13
In the above equation, A ( E j ) is a species at a specific energy state Ei,  denotes the concentration, M is a bath-gas molecule, kij is the second-order rate constant for the collisional energy transfer of a molecule from state Ej to state Ei and km ( Ei ) denotes the first-order microcanonical rate constant for a particular isomerization or dissociation channel m. An energy transition probability Pji may be defined as k Pji = ji , (3) Z
d [ A ( Ei )] = ∑ kij [ M ] ⎡ A ( E j ) ⎤ − ∑ k ji [ M ] ⎡ A ( Ei ) ⎤ − ∑ km ( Ei ) ⎡ A ( Ei ) ⎤ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ j dt j m
where Z is the collision rate constant. Obviously, we have be re-written in terms of Pij as
d ⎡ A ( Ei ) ⎤ ⎣ ⎦ dt
= 1 . The master equation may
⎧ ⎫ = Z [ M ] ⎨∑ Pij ⎡ A ( E j ) ⎤ − ⎡ A ( Ei ) ⎤ ⎬ − ∑ km ( Ei ) ⎡ A ( Ei ) ⎤ . ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎩ j ⎭ m
Microcanonical rate constants
The unimolecular rate constant k(E) is given by the RRKM expression11-13
k (E) = L
† Q † W ( Evr ) Q hρ ( E∗ )
† In the above equation, W ( Evr ) denotes the sum of states of the transition state with energy ≤ E † , † while ρ ( E ∗ ) is the density of states of the energized molecule at energy E ∗ = Evr + E0 , and E0 is the energy barrier. To compute both the density and sum of states we have used the WhittenRabinovitch approximation.14 To account for the effect of adiabatic rotations, the above expression contains the factor Q † Q , where Q † and Q are partition functions for the adiabatic rotations of the transition state and the energized molecule respectively. Finally, the factor L is introduced to account for the reaction path degeneracy. Equation (5) is easily implemented for reactions with a pronounced barrier, where the critical geometry is independent of the internal energy as opposed to barrierless reactions, for which the variational transition state theory is imperative. The rate constant at the high-pressure limit is calculated using the expression
k∞ = L
kT Q † exp(− E0 / kT ) h Q
where k is the Boltzmann constant and h the Plank constant. The RRKM program used here15 has been widely tested for a number of earlier studies.
The transition probability matrix
Various models have been proposed to calculate the transition probabilities, such as the exponential down model, the stepladder model, the Poisson model, the Gaussian model and the biased random walk model. All of these must satisfy the normalization condition ∑ Pji = 1 and j detailed balancing Pji ρi exp ( − Ei kT ) = Pij ρ j exp ( − E j kT ) (7)
which ensures that the distribution of a non-reacting system at long times is the Boltzmann distribution. Here, ρi denotes the density of states of the ith energy grain. In this work, we adopt an exponential down model. The transition probability for this model is given by
Pij = C j exp ⎡ −α ( E j − Ei ) ⎤ ⎣ ⎦
Pij = Ci
ρi exp ⎡ − (α + β ) ( Ei − E j ) ⎤ ⎣ ⎦ ρj
where Ci s are the normalization factors, α is the reciprocal of the average energy transferred in the ‘down’ transitions, α = 1 ∆E d , and β = 1 kT . In the Monte Carlo treatment of collision energy transfer, the normalization factors are first calculated for the transition probability. The algorithm employed here is the recursive algorithm outlined in Ref. 11. Since the probabilities depend on the density of states and temperature, the matrix has to be separately calculated for each isomer separately and at each temperature of interest. Figure 1 presents the variation of Pij on the energy of the final state for an excited HOCO molecule. The Boltzmann distribution is approximately reached after multiple collisions as seen in Figure 2. Computationally the challenge of implementing the collision probability matrix lies in the balance of computational speed and accuracy as well as the memory requirements, as this matrix is often excessively large. We have devised a method to reduce the number of elements in this matrix by considering that the energy transfer probability vanishes exponentially with the energy step size and thus the likelihood of large energy transfer or the non-diagonal elements of the transition probability matrix are negligibly small. The approximation employed in the present work is then to use only the elements close to the diagonal, as illustrated in Fig. 3. For this example case, the number of elements required are reduced from 4200 × 4200 to 117 × 117. To ensure that the accuracy is maintained, we examine the normalization condition is met. While in the full matrix ∑ Pji = 1 , in the reduced matrix, j ∑ Pji = 0.991. The above procedure enables us to calculate and store the transition probabilities j keeping the energy grain size uniform over the entire energy range. The current approach is quite different from that of the Multiwell6-8 collection of programs, where double arrays are used to reduce the memory requirements.
Selection of initial energies
The first step of the stochastic process in the Monte Carlo algorithm for any given molecule is the determination of its initial energy. This is easily done if the energy distribution of the adduct formed by the chemical activation process is known. Consider the system A+B AB* +M AB Rate of formation of AB is
ν r ( E ) dE = kr ,∞ [ A][ B ] g ( E ) dE ,
where g ( E ) dE is the fraction of collisions leading to the formation of AB* in the energy range E to E + dE and kr ,∞ is the combination rate constant at the high-pressure limit. It can be shown that11
g ( E ) dE = ka ( E ) ρ ( E ) exp ( − E kT ) dE
∫ k ( E ) ρ ( E ) exp ( − E kT ) dE
Knowing the fractional distribution of the collisions, initiating the stochastic process is straightforward. Specifically, a random number R1 (0 ≤ R1 ≤ 1) is chosen from the unit interval uniform distribution and the initial energy E is selected such that
R1 = ∫ g ( E ′)dE ′ .
Figure 4 shows the histogram for the initial energies selected for 105 molecules for a particular case considered in this study. The narrow energy distribution is typical for a chemically activated system.
Implementation of stochastic method
Once the initial energy of the chemically activated molecule is determined, the isomerization or collisional energization/de-energization of this hot molecule is determined by a stochastic algorithm, till the molecule stabilizes or decomposes to form products. The stochastic approach is similar to Multiwell,6-8 based on the algorithm proposed by Gillespie.17,18 Specifically events along the path of reaction are treated as Markovian. The time scale τ at which the system of state k changes state is given by
− ln( R2 ) , ∑ ank
where R2 (0 ≤ R2 ≤ 1) is a random number drawn from unit interval uniform distribution and ank is the probability rate for the transition to the accessible state n. For a unimolecular or chemically activated process, ank in the above equation are replaced by the microcanonical rate constants and the collision frequency. To determine the outcome of the state change or whether the molecule undergoes isomerization or activation/deactivation, another random number R3 (0 ≤ R3 ≤ 1) is chosen from unit interval uniform distribution and the state n is chosen such that
∑ an′k < R3 ∑ ank ≤ ∑ an′k .
n′ =1 n n′=1
If the molecule undergoes isomerization the species is accordingly updated and the above process is repeated at the same energy. All reactions except for the decomposition channels are treated reversibly and rates are calculated in either direction for any given channel. If the collisional energy transfer process is chosen, activation or deactivation is decided as follows. For a molecule in a given state (or energy grain) i, the total probability for activation is given by
Pact ,i = ∑ Pji
Similarly, the total probability for deactivation is
Pdeact ,i = ∑ Pji
We proceed to choose yet another random number R4 and select activation if R4 ≤ Pact , and deactivation otherwise. The amount of energy transferred is then calculated using yet another random number and the transition probabilities previously calculated. The final energy state n upon collision is given by
R5 = ∑ Pji for activating collisions and
j ≥i n n
R5 = ∑ Pji for deactivating collisions
The final energy is found by interpolation of the discrete probabilities. While the above procedure can ideally be applied to a system composed of many molecules, the memory requirements will be prohibitive and instead we resort to simulating the random walk for one molecule at a time and then averaging the results for an ensemble of stochastic trials. A difficulty faced in modeling the stabilizing process is the criteria for declaring an isomer to be ‘stable.’ While the collisional energy transfer can continue ad-infinitum, a molecule trapped at the bottom of a potential well rarely gathers enough energy to climb the potential barrier and it may be termed ‘stable’. One of the ways to recognize that the molecule has stabilized may simply be to record enough statistics and compare this with the Boltzmann distribution. This is however computationally expensive. We follow here a simpler approach by marking the energy grain i for which the total probability for activation is equal to that for deactivation i.e. Pact ,i = Pdeact ,i . Heuristically, once a molecule is deactivated to this energy level, the probability of it being re-energized is equal to that of its being de-energized and it is highly unlikely that it will climb up the energy ladder significantly, and as such the molecule can be considered to be stabilized. Figure 2 shows the location of the sink for a typical case considered for this study, calculated using the above criterion. It should be recognized that the criterion
discussed above is more rigorous than placing the sink at an arbitrary level below the reaction threshold. Counters are initiated for each isomer in the reaction network as also for the reactants and products, to which the hot molecule may decompose. The stochastic process is repeated for several thousand molecules, with the accuracy of results increasing with the increase in the number of trials. The rate constant for the jth channel is given by
k j = F j k r ,∞ = F j
kd ,∞ K eq
where Fj denotes the fraction of molecules and K eq is the equilibrium constant for the entrance channel.
Application to representative systems
This section discusses the application of the above code to several chemically activated reactions. Comparison will be made with the solution obtained by using Troe’s modified strong collision approximation.
CO2 + H → CO + OH
This is the simplest of the test cases considered here, with only one isomer and one exit channel,
CO2 + H
CO + OH
HOCO The energy diagram for this reaction system can be seen in Figure 5. The reaction is known to proceed through the formation of the intermediate, HOCO, and the stabilization of this complex upon collision with the bath gas is responsible for the complex temperature and pressure dependence of this reaction. The reaction has been extensively studied experimentally and computationally (see, for example Ref. 19). Our intention here is simply to observe the difference, if any, between the rates obtained upon solving the Master Equation and those with the modified strong collision model, for a single well scenario. For this purpose, the energies and RRKM parameters have been taken from the study of Larson et.al.20 and are listed in Table 1. Figure 6 presents a comparison of the rate constants, obtained by the Monte Carlo method with that using Troe’s modified strong collision model over a temperature range of 300-2000 K at
atmospheric pressure. The calculations use a <Edown> value of 157.1 cm-1 for N2 as the bath gas. For the Monte Carlo simulations, 106 molecules are used. As can be seen in Fig. 6, the modified strong collision model systematically over-predicts the rate constant of the combination reaction, CO2 + H → HOCO over the entire range of temperature. The chemically activated reaction channel CO2 + H → CO + OH is affected by the modified strong collision model for temperature < 800 K, where the stabilization channel becomes increasingly fast.
Reactions on the C3H5 potential energy surface
The ability of the code to handle multiple well systems efficiently and accurately is demonstrated for reactions on the C3H5 potential energy surface. The C2H2 + CH3 has been extensively studied, both experimentally21-24 and computationally.25-27 The reactions considered here include the chemically activated reaction of C2H2 + CH3, the mutual isomerization, stabilization and dissociation of the C3H5 isomers: CH3CHCH, CH3CHCH2 and aC3H5 (allyl).
pC3 H4 + H CH3CCH2* C2H2 + CH3 CH3CHCH* aC3H5 * aC3H4 + H
The potential energy diagram is presented in Figure 7, and the RRKM parameters, taken from Davis et. al.27, are listed in Table 2. Here we compare the rates calculated for this reaction using the Monte Carlo code with the previously published results of Davis et. al,27 obtained using Troe’s modified strong collision model. As in the previous case, the Monte Carlo simulation was run for 106 molecules, for a pressure of 2 bar and with Ar as the bath gas. We have chosen an <Edown> value of 260 cm-1. Figure 8 shows the computed and experimental results for the C2H2 + CH3 → reactions. The experimental data of Hidaka et.al.23 and Kislitsyn et.al.24 are in the high-temperature regime,
where the stabilization reactions are negligible. On the other hand the stabilization channel is expected to be dominant in the low-temperature experiments of Holt and Kerr22 and GarciaDominguez et.al.21 As seen in the figure, the Monte Carlo predicts the rate constants in both the regimes within the uncertainty of the experimental values. Of interest here, is the prediction of the rates for stabilization channels at higher temperature. Conclusions similar to the HOCO case can be drawn. Troe’s modified strong collision model predicts the rate constants well for lower temperatures for all channels, but the overprediction of the stabilization channels is evident at higher temperatures for reaction channels that require collision stabilization. The number of molecules required for the simulation depends on the disparity of the rates of the competing processes. Thus, for a system characterized by both shallow and deep wells, a large number of stochastic trials will be necessary to capture a few molecules in the shallow well and thereby predict the rate coefficients to this channel accurately. For C2H2 + CH3 → reactions, 106 molecules were found to give an adequate description of the system for the purpose of this study, as seen in Fig. 9. It is seen that the relative change in the rate constant values on adding more molecules is small, and convergence is achieved for > 2×105 molecules
A code that employs the Monte Carlo method to solve solving the Master Equation of collision energy transfer and calculate the rate constants of chemically activated and unimolecular reactions of arbitrary complexity is developed and described in detail. The code has been tested for a couple of systems, and the solutions compared with those obtained using Troe’s modified strong collision model. The results demonstrate the inadequacy of the latter approach for collision stabilization channels.
Acknowledgement The work was supported by the National Science Foundation (CTS9874768) and by the Air Force Office of Scientific Research ( FA9550-04-1-0008).
1. Troe J.; J. Chem. Phys. 1977, 66, 4745; 1977, 66, 4758. 2. Mokrushin, V., Bedanov, V., Tsang W., Zachariah, M.R., and Knyazev, V.D.; ChemRate,
Version1.16, National Institute of Standards and Technology: Gaithersburg, MD 20899.
3. Vereecken, L., Huyberechts, G., and Peeters, J.; J. Chem. Phys. 1997, 106, 6564. 4. Klippenstein, S.J., Wagner, A.F., Robertson, S.H., Dunbar, R., and Wardlaw, D.M.;
Variflex Software, Version 1.0, 1999.
5. Smith, S.C., Diau, E.W.-G., and Schranz, H.W.; UNIRATE, 2001. 6. Barker, J.R.; Int. J. Chem. Kinet. 2001, 33, 232. 7. Barker, J.R., and Ortiz, N.F.; Int. J. Chem. Kinet 2001, 33, 246 8. Barker, J.R.; Multiwell Program Suite, version 1.3.2, Ann Arbor, MI 2003. 9. Pilling M.J., and Robertson, S.H.; Annu. Rev. Phys. Chem. 2003, 54, 245. 10. Barker, J.R., and Golden, D.M.; Chem. Rev. 2003, 103, 4577. 11. Holbrook, K.A., Pilling, M.J., and Robertson, S.H.; Unimolecular Reactions, 2nd ed.; Wiley: Chichester, 1996. 12. Gilbert, R.G., and Smith, S.C.; Theory of Unimolecular and Recombination Reactions, Blackwell-Scientific, Oxford, 1990. 13. Forst W.; Theory of Unimolecular Reactions, Academic, New York, 1973. 14. Whitten, G.Z., and Rabinovitch, B.S.; J. Chem. Phys. 1963, 38, 2466. 15. Wang, H.; Detailed Kinetic Modeling of Soot Particle Formation in Laminar Premixed Hydrocarbon Flames. Ph.D. Thesis, 1992, The Pennsylvania State University. 16. Oref, I., and Tardy, D.C; Chem. Rev., 1990, 90, 1407. 17. Gillespie, D.T.; J. Comp. Phys. 1976, 22, 403. 18. Gillespie, D.T.; J. Comp. Phys. 1978, 28, 395. 19. Baulch, D.L., Cobos, C.J., Cox, R.A., Frank, P., Hayman, G., Just, Th., Kerr, J.A., Murrells, T., Pilling, M.J., Troe, J., Walker, R.W., and Warnatz, J.; J. Phys. Chem. Ref.
Data 1992, 21, 737.
20. Larson, C.W., Stewart, P.H., and Golden, D.M.; Int. J. Chem. Kinet. 1988, 20, 27. 21. Garcia-Dominguez, J.A., and Trotman-Dickenson, A.F.; J. Chem. Soc. 1962, 940. 22. Holt, P.M., and Kerr, J.A.; Int. J. Chem. Kinet. 1977, 9, 185.
23. Hidaka, Y., Nakamura, T., Miyauchi, A., Shiraishi, T., and Kawano, H.; Int. J. Chem.
Kinet. 1989, 21, 643.
24. Kislitsyn, M.N., Slagle, I.R., and Knyazev, V.D.; Proc. Combust. Inst. 2002, 29, 1237. 25. Diau, E.W., Lin, M.C., and Melius, C.F.; J. Chem. Phys. 1994, 101, 3923. 26. Wang, B., Hou, H., and Gu, Y.; J. Chem. Phys. 2000, 112, 8458. 27. Davis, S.G., Law, C.K., and Wang H.; J. Phys. Chem. A 1999, 103, 5889. 28. Kee, R.J., Warnatz, J., and Miller, J.A.; Sandia Report SAND 83-8209; Sandia National Laboratories: Albuquerque, NM, 1983.
Table 1: RRKM parameters of the HOCO system
υ , cm-1
3456, 1833, 1260, 1071, 615, 615 0.348 (1,2) external inactive; (1,1) external active
L-J param σ = 4 Å; ε / k = 200 K
υ , cm-1
2400, 2100, 1200, 870, 560 0.339 (1,2) external inactive; 10.8 (1,1) external active 3456, 1900, 900, 425, 175 0.173 (1,2) external inactive; 4.71 (1,1) external active
υ , cm-1
The numbers in parentheses are the symmetry number and the dimension of the rotor, in that order.
Table 2: RRKM parameters of the C3H5 system
υ , cm-1
189 395, 597, 773, 798, 911, 1027, 1079, 1232, 1357, 1437, 1437, 1647, 2943, 2950, 3004, 3049, 3164 0.303 (1,2) external inactive; 1.929 (1,1) external active
L-J paramb CH3CCH2
σ = 4 Å; ε / k = 200 K
175, 305, 468, 856, 874, 913, 1015, 1067, 1344, 1372, 1410, 1425, 1711, 2899, 2960, 2983, 3009, 3075 0.280 (1,2) external inactive; 2.520 (1,1) external active 412, 514, 536, 763, 788, 905, 983, 1008, 1176, 1235, 1376, 1469, 1470, 3055, 3058, 3065, 3158, 3160 0.317 (1,2) external inactive; 1.835 (1,1) external active 448i, 47c, 231, 436, 470, 530, 632, 723, 743, 803, 1365, 1374, 1873, 3014, 3175, 3187, 3299, 3387 0.212 (1,2) external inactive; 1.307 (1,1) external active 665i, 153, 357, 410, 454, 632, 655, 920, 1012, 1027, 1364, 1428, 1432, 2091, 2961, 3033, 3034, 3380 0.278 (1,2) external inactive; 2.229 (1,1) external active 2038i, 147, 240, 352, 590, 805, 882, 1000, 1001, 1341, 1419, 1426, 1815, 2325, 2922, 2974, 3006, 3027 0.280 (1,2) external inactive; 2.556 (1,1) external active 2092i, 444, 616, 670, 869, 893, 914, 1018, 1024, 1096, 1194, 1378, 1600, 1805, 2988,2 3056, 3083, 3107 0.381 (1,2) external inactive; 1.286 (1,1) external active 503i, 59c, 225, 338, 400, 640, 732,5 930, 1011, 1020, 1364, 1426, 1428, 2125, 2951, 3016, 3024, 3374 0.264 (1,2) external inactive; 2.355 (1,1) external active 1864i, 296, 394, 676, 808, 810, 945, 990, 1066, 1109, 1377, 1405, 1644, 2145, 2980, 2995, 3109, 3113 0.289 (1/2d,2) external inactive; 2.575 (1,1) external active 425i, 220, 302, 363, 389, 840, 855, 862, 969, 985, 1067, 1369, 1424, 1967, 3038, 3049, 3114, 3127 0.271 (1,2) external inactive; 2.248 (1,1) external active 722i, 140, 375, 431, 486, 821, 825, 856, 977, 984, 1061, 1375, 1425, 1938, 3047, 3052, 3128, 3141 0.287 (1,2) external inactive; 2.123 (1,1) external active
υ , cm-1
a a -1
CH2CHCH2 υ , cm-1 B0, cm
TS1 υ , cm-1 B0,a cm-1 TS2 υ , cm TS3 υ , cm
TS4 υ , cm-1 B0,a cm-1 TS5 υ , cm-1 B0,a cm-1 TS6 υ , cm-1 B0,a cm-1 TS7 υ , cm-1 B0,a cm-1 TS8 υ , cm
a -1 -1
The numbers in parentheses are the symmetry number and the dimension of the rotor, in that order.
The vibrational motion is
replaced by an internal free rotor, and the frequency is not used in the RRKM calculation.
The symmetry number of ½ was implemented to
account for the asymmetric transition structure which accounts for the proper path degeneracy.
Figure 1. The dependence of Pij on the energy of the final state for a typical excited molecule initially in the 1000th energy grain. The above calculation was done for the HOCO system at 300 K, using the exponential down model for energy transfer.
Energy (cm )
0 100 0 5 10-8 1 10-7 1.5 10-7 2 10-7 2.5 10-7 3 10-7 3.5 10 -7 Time (sec)
Figure 2. Trajectory of a stabilizing molecule: The hot adduct gets trapped in the well in the absence of reactions. Above calculation was done for the HOCO system at 300K, in the absence of exit channels. The solid (horizontal) line shows the ‘sink’ as calculated by the code. Inset shows that the Boltzmann distribution is achieved at long enough times, when the isomer is considered to have stabilized.
No. of elements Full Matrix : 4200 Reduced : 117
0 100 850
Figure 3. Approximation employed to reduce the size of the transition probability matrix. The elements far away from the diagonal are negligibly small, compared to the ones near the diagonal and may be neglected without loss in accuracy. This reflects the notion that large energy transfer is highly unlikely. Solid curve: the exact form of Pij, dotted curve: the approximate, truncated form. The matrix size for this case was reduced from 4200 × 4200 to 117 × 117 elements.
E = 32.5 kcal/mol
Number of molecules
0 0 10 20 30 40 50 60 70
Figure 4. Histogram showing the initial energy distribution of HOCO produced from H addition to CO2 at 300 K. The energy barrier for CO2 + H is 32.5 kcal/mol.
Relative Energy (kcal / mol )
TS 1 32.5
TS 2 33.1 CO + OH
CO2 + H HOCO
Figure 5. Potential energy diagram for CO2 + H → products .
CO2 + H = CO + OH
Bath gas = N2 P = 1 atm
CO2 + H = HOCO
0.5 1.0 1.5 2.0 2.5 3.0 3.5
Figure 6. Comparison of the rate coefficients for CO2 + H → products. Open symbols: Troe’s modified strong collision model , Solid symbols: Values obtained from solution of the Master Equation using the Monte Carlo approach. <Edown> = 157.1 cm-1 N2.
TS3 67.4 TS2 58.8 pC3H4 + H 53.3
TS6 63.8 TS7 57.5
Relative Energy (kcal / mol )
TS5 TS1 56.5 56.0
TS8 60.0 aC3H4 + H 54.4
C2H2 + CH3 46.2
Figure 7. Potential energy diagram of the C3H5 system (from ).
Bath gas = Ar <Edown> = 260 cm-1
k (cm3/mol s)
k (cm3/mol s)
102 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
Figure 8. Comparison of the rate coefficients for reaction C2H2 + CH3 → products. Dashed lines: modified strong collision model, from , Solid lines and filled symbols: solution of the Master Equation. Open symbols: Experimental data (squares, Hidaka et.al.,23 circles, Kislitsyn et.al.,24 diamonds, Diau et.al.,25 remodeling of GarciaDominguez, et.al.,21 triangles, Holt et.al.22).
C2H2 + CH3 -> products (2 atm, 1500 K)
1011 Propyne + H
Allene + H
106 2 4 6 8 10 1.2 106
No. of molecules x 105
Figure 9. Effect of varying the number of molecules on Monte Carlo rate calculations for reactions C2H2 + CH3 → products (M = Argon).