Automatic estimation of pressure-dependent rate coefficients w

View Article Online / Journal Homepage / Table of Contents for this issue PCCP Dynamic Article Links Cite this: Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 PAPER www.rsc.org/pccp Automatic estimation of pressure-dependent rate coefficientsw Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C Joshua W. Allen,a C. Franklin Goldsmithab and William H. Green*a Received 30th August 2011, Accepted 27th October 2011 DOI: 10.1039/c1cp22765c A general framework is presented for accurately and efficiently estimating the phenomenological pressure-dependent rate coefficients for reaction networks of arbitrary size and complexity using only high-pressure-limit information. Two aspects of this framework are discussed in detail. First, two methods of estimating the density of states of the species in the network are presented, including a new method based on characteristic functional group frequencies. Second, three methods of simplifying the full master equation model of the network to a single set of phenomenological rates are discussed, including a new method based on the reservoir state and pseudo-steady state approximations. Both sets of methods are evaluated in the context of the chemically-activated reaction of acetyl with oxygen. All three simplifications of the master equation are usually accurate, but each fails in certain situations, which are discussed. The new methods usually provide good accuracy at a computational cost appropriate for automated reaction mechanism generation. 1 Introduction Thermal unimolecular reactions—the set of reactions involving a single reactant molecule—are unusual in that they require an inert third body to proceed. The third body provides or removes the energy needed for reaction via collisions with the reactant molecule. The rate of these collisions depends on the concentration of the inert, which in turn is related to the pressure of the system; therefore, under conditions where this collision rate is rate-limiting, the observed phenomenological rate coefficient k(T,P) is a function of both temperature T and pressure P. An activated species (isomer) can be formed either as the product of an association reaction (chemical activation) or via collisional excitation (thermal activation). Once activated, multiple isomerization and dissociation reactions may become competitive with one another and with collisional stabilization; these combine to form a network of unimolecular reactions. Often the reactive events occur more rapidly than collision events, and an excited molecule will traverse multiple reactive events—appearing to ‘‘skip’’ over intermediate species— before being collisionally stabilized. This suggests that there is a net reaction rate from each isomer and product to every other isomer and product set in the network, not just those directly adjacent. a Dept. of Chemical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Ave, Cambridge, MA 02139, USA. E-mail: [email protected]; Fax: 617-324-0066; Tel: 617-253-4580 b Dept. of Inorganic Chemistry, Fritz Haber Institute, Berlin, 14195, Germany w Electronic supplementary information (ESI) available. See DOI: 10.1039/c1cp22765c This journal is c the Owner Societies 2012 Our interest in unimolecular reaction networks is colored by our intended application: automatic reaction mechanism generation. Our automatic mechanism generation code, RMG,1 utilizes a rate-based approach to model enlargement, and it is therefore important to account for the pressure dependence of chemically-activated reaction networks in order for RMG to generate a valid mechanism. RMG often considers thousands of these reaction networks in the course of constructing a reaction mechanism, so we require a method that is fully automatable, successful over a wide range of conditions, reasonably accurate, and not too computationally intensive. At the time it is making an estimate of k(T,P), RMG has only limited information available about the molecules and reactions, as discussed below. This makes it more challenging to accurately estimate all the k(T,P). The importance of bimolecular collisions in unimolecular reactions was first proposed by Lindemann in 1922.2 It was soon recognized by Hinshelwood and others that a rigorous treatment of these processes required consideration of molecular energy levels.3 The RRKM expression for the microcanonical rate coefficient k(E) was derived in the early 1950s.4–6 In the late 1950s master equation models of chemical systems began appearing,7–11 including an early linear integraldifferential equation formulation by Widom.12 Analytical solutions for a variety of simple models soon followed,13–15 as did the first numerical approaches.16 Numerical methods— which are required for complex unimolecular reaction networks— became much more attractive in the 1970s with the appearance of new algorithms, including Gear’s method for solving stiff systems of ordinary differential equations17 and efficient algorithms for calculating the density of states.18–20 In the 1990s computing power had increased to the point where it Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 1131 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online was practical to solve them numerically by discretizing the integrals over energy. Methods for directly solving the discretized master equations include integration using a stiff ODE solver and stochastic simulation using the Gillespie method. Both of these methods have the benefit of not introducing additional approximations (in the latter only if sufficient sampling of events is taken). However, both methods suffer when wide ranges of time scales are encountered, as they often are in chemical kinetics. Furthermore, with either approach it is not clear how to extract the time-independent rate constants k(T,P) from the full solution. Barker et al. have developed the stochastic approach into the software package MultiWell.21 Over the last ten to fifteen years a number of numerical techniques have been developed with the goal of simplifying a full master equation formulation into a single set of phenomenological rate coefficients k(T,P). The methods vary in the simplifying assumptions made and the numerical procedures used, which in turn affects the computational expense required and the accuracy and robustness of the result. A sampling of methods and associated software packages is provided in the following paragraphs. Several approaches are based on computing the eigenmodes of the full master equation matrix. When the conventional phenomenological-kinetics description of the chemical system is applicable, certain ‘‘chemically-significant’’ eigenmodes become separable from the large number of internal energy relaxation eigenmodes. The chemically-significant eigenvalues match the eigenvalues of the matrix of phenomenological rates, and therefore can be used to extract values for each k(T,P), although this may not be straightforward.22 Eigenvalue determination is computationally expensive and numerically challenging for large, stiff systems; e.g. at low temperatures and high pressures, solvers often return nonphysical or complex eigenvalues. There are multiple software packages that utilize eigenvalue approaches, including UNIMOL,23 ChemRate,24 Variflex,25 and the recently-released MESMER.26 On the opposite end of the scale are methods based around the strong collision approximation, wherein all collisions are assumed to irreversibly stabilize the energized species. This decoupling of the energy grains in the system is a huge computational savings, at the cost of neglecting the effects of collisions which cause small energy changes. Despite this gross simplification, the modified strong collision method appears to provide reasonably accurate results, typically within an order of magnitude of the true rate coefficients.27 This may be partially due to the modification of the collision frequency by a temperature-dependent efficiency factor.28–30 For multiwell systems, this method was codified by Chang, Bozzelli, and Dean into the software package CHEMDIS.31 The modified strong collision method is implemented in the automatic mechanism generation software package XMG, and in early and current versions of RMG. More recently, Green and Bhatti have proposed a new approximate method for computing k(T,P) based around the pseudo steady-state approximation at high energies coupled with the approximation of a Boltzmann population of low energy levels.32 This method maintains the collision model of the full analysis for the important high-energy states, so a 1132 Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 more accurate result is expected when compared to the modified strong collision method. However, it is not clear yet just how much more expensive and more accurate this method is. Also, the multiple-well formulation of the reservoir state method as described by Green and Bhatti is cumbersome, requiring the method to be repeated using each isomer well and reactant set as a starting point. The present work will develop a formulation that does not require this repetition. A key input to all k(T,P) calculations is the density of states for each species. In the context of automated reaction mechanism generation it is not yet possible to determine this information by performing quantum chemistry calculations on every species, as these are still too slow for large molecules, and nontrivial to automate. Bozzelli and coworkers developed a method for estimating the density of states from heat capacity data (obtained from group additivity estimates) by fitting three pseudofrequency-degeneracy pairs.33 The fitted frequencies do not necessarily reflect the real vibrational degrees of freedom of the molecule, and the harmonic-oscillator form is not a good representation for hindered internal rotors. The primary objective of the current work was to develop a general framework for automatic estimation of pressuredependent rate coefficients in reaction networks of arbitrary size and complexity. Within this framework, various methods of estimating the density of states and extracting phenomenological rate constants from the master equation matrix could be developed and evaluated. In the former, a new method based on characteristic functional-group frequencies was created to more accurately reflect the real molecular degrees of freedom, including hindered rotors. In the latter, the chemically-significant eigenvalues, modified strong collision, and reservoir state methods were directly compared using identical input parameters. Section 2 provides a brief review of the formulation of the master equation model. Section 3 describes our framework for estimating the phenomenological rate coefficients for unimolecular reaction networks. Section 4 describes the methods used to estimate the density of states, focusing on a new method based on characteristic functional group frequencies. Section 5 describes the methods used to estimate the phenomenological rate coefficients from the full master equation formulation. Section 6 evaluates the methods of the previous sections using a case study: the reaction of acetyl radical with oxygen, a process important in atmospheric chemistry. 2 Background: the master equation Descriptions of the full master equation formulation are quite ubiquitous,34–38 and so only a brief summary will be given here. The system of interest is a set of chemically reactive molecular configurations—local minima on a potential energy surface—divided into unimolecular isomers and bimolecular reactants or products. In our vernacular, reactants can associate to form an isomer, while such association is neglected for products. These configurations are connected by chemical reactions to form a network. The system also consists of an excess of inert gas M, representing a thermal bath; this allows This journal is c the Owner Societies 2012 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online for neglecting all collisions other than those between an isomer and the bath gas. An isomer molecule at sufficiently high internal energy E can be transformed by a number of possible events:  The isomer molecule can collide with any other molecule, resulting in a change in energy from E to E 0  The isomer molecule can isomerize to an adjacent isomer at constant E  The isomer molecule can dissociate into any directly connected bimolecular reactant or product channel The dependent variables describing the isomers are pi(E,t), the population distribution of isomer i over energy E at time t. From statistical mechanics we expect a Boltzmann distribution bi(E,T) to result at long times: lim pi ðE; tÞ ¼ xi1 t!1 ri ðEÞebE ¼ xi1 bi ðE; TÞ Qi ðbÞ ð1Þ Above, xiN is the total population of isomer i at equilibrium, ri(E) and Qi(b) are the rovibrational density of states and corresponding partition function of isomer i, and b  (kBT)1. Assuming that inelastic collisions are much more common than reactive collisions, we can treat the bimolecular reactants as thermalized, and represent the concentrations of the nth reactant configuration by ynA(t) and ynB(t), where A and B distinguish between the two reactants. Neglecting the dependencies of the microcanonical rates and other quantities on the angular momentum quantum number J or any other quantum numbers besides E,39 and neglecting reactions which do not proceed through one of the isomers, the master equation can be written as Z 1 d pi ðE; tÞ ¼ oi ðT; PÞ Pi ðE; E 0 ; TÞpi ðE 0 ; tÞdE 0 dt 0  oi ðT; PÞpi ðE; tÞ þ N isom X þ N reac X jai n¼1 kij ðEÞpj ðE; tÞ  jai kji ðEÞpi ðE; tÞ ynA ðtÞynB ðtÞfin ðEÞbn ðE; TÞ Nreac þNprod X  N isom X n¼1 gni ðEÞpi ðE; tÞ ð2Þ d d ynA ðtÞ ¼ ynB ðtÞ dt dt Z 1 N isom X gni ðEÞpi ðE; tÞdE ¼ i¼1  0 N isom X i¼1 ynA ðtÞynB ðtÞ Z 0 1 ð3Þ fin ðEÞbn ðE; TÞdE where oi(T,P) is the collision frequency; P(E,E 0 ,T) is the probability of collisional transfer from energy E 0 to E; kij(E), fin(E), and gni(E) are the microcanonical rate coefficients for isomerization, association, and dissociation, respectively; This journal is c the Owner Societies 2012 bn(E,T) is the Boltzmann distribution for bimolecular reactant channel n; and Nisom, Nreac, and Nprod are the numbers of isomers, bimolecular reactant channels, and bimolecular product channels, respectively. In eqn (2), the first pair of terms correspond to collision, the second pair to isomerization, and the final pair to association/dissociation. Eqn (2) applies to isomers, while eqn (3) applies to bimolecular reactants. Note that there are almost always additional reactions creating and destroying all the species, so eqn (2) and (3) are for an idealized situation. There is some subtlety to the formulation of the master equation, in that some treatments are based on the total energy, while others use the active internal energy. This has been discussed in the literature previously.40,41 This work utilizes the one-dimensional total energy model as described by Miller and Klippenstein. In particular, the densities of states contain all of the vibrational and rotational modes of the system, which implies that we are treating all rotational modes as active; this approximation appears to be increasingly accurate as temperature increases.41 Eqn (2) and (3) are nonlinear, both due to the presence of the bimolecular reactant terms and because both oi and P(E,E 0 ) depend on the composition, which is changing with time. The rate coefficients can be derived from considering the pseudo-first-order situation where ynA(t) { ynB(t), and all y(t) are negligible compared to the bath gas M. From these assumptions the changes in o,P(E,E 0 ), and all ynB can be neglected, which yields a linear equation system. To extract the phenomenological rate coefficients numerically, it is helpful to discretize the energy E into Ngrains grains {Er}. This converts the linear integro-differential equations into a system of first-order ordinary differential equations with the form d P ¼ MP dt 2 M1 6 6 K21 6 6 6 . 6 . 6 . 6 M¼6 6 ðg ÞT 6 11 6 6 6 ðg21 ÞT 6 4 .. . 3 2 p1 7 6 6 p 7 6 2 7 7 6 6 . 7 6 . 7 6 . 7 7 6 P¼6 7 6y 7 6 1A 7 7 6 7 6 6 y2A 7 7 6 5 4 .. . K12 ... F11 b1 y1B F12 b2 y2B ... M2 ... F21 b1 y1B F22 b2 y2B .. . .. . .. . .. . ðg12 ÞT ... h1 0 ðg22 ÞT ... 0 h2 .. . .. .. . .. . 7 ...7 7 7 .. 7 7 .7 7 7 ...7 7 7 7 ...7 7 5 .. . . Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 3 ð4Þ 1133 View Article Online where the elements of the vectors pi are such that (pi)r = pi(Er). The diagonal matrices Kij and Fin and the vectors gni contain the microcanonical rate coefficients for isomerization, association, and dissociation, respectively:  R Er þDEr =2 1 ðKij Þrs ¼ DEr Er DEr =2 kij ðEÞdE r ¼ s ð5aÞ 0 ras Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C ðFin Þrs ¼  R Er þDEr =2 1 DEr Er DEr =2 fin ðEÞdE 0 ðgni Þr ¼ 1 DEr Z Er þDEr =2 Er DEr =2 r¼s ras ð5bÞ gni ðEÞdE ð5cÞ Above, r is the index for the discretized energies of isomer i, and s is the index for the discretized energies of isomer j or bimolecular channel n. The matrices Mi represent the collisional energy transfer probabilities of isomer i minus the rates of reactive loss to other isomers and to reactants and products: 8 isom >   NX > > oi ðT; PÞ ðPi Þ  1  ðKij Þrs > rr > > > jai > < r¼s ðMi Þrs ¼ Nreac þNprod X > > > >  ðgni Þr > > > n¼1 > : oi ðT; PÞðPi Þrs ras ð6Þ The collisional energy transfer probabilities (Pi)rs for isomer i are given by Z Er þDEr =2 Z Es þDEs =2 1 ðPi Þrs ¼ Pi ðE; E 0 ; TÞdE 0 dE DEr DEs Er DEr =2 Es DEs =2 ð7Þ The scalars hn are simply the total rate coefficient for loss of reactant channel n due to chemical reactions: hn ¼  N grains isom NX X i¼1 r¼1 ynB ðFin Þrr bn ðEr Þ ð8Þ 3 A general framework for processing reaction networks The procedure, outlined in Fig. 1, is used to estimate phenomenological rate coefficients k(T,P) for a reaction network of arbitrary connectivity and complexity based exclusively on high-pressure-limit parameters. The parameters needed are:  The species in the network. For each species the thermodynamic parameters (DHof (298 K), DSof (298 K), and a model for Cp(T)) are needed. For those species that are isomers in the network, the chemical structure (connectivity only), molecular weight, and Lennard-Jones parameters are also needed.  The reactions in the network. For each reaction the highpressure limit rate coefficient kN(T) is needed.  The bath gas. The molecular weight and the LennardJones parameters are required.  The collisional transfer probabilities model. Throughout this paper the single exponential down model is used due to its dependence on a single parameter; more complex models can 1134 Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 be used given sufficient information, but often this is not available. Other inputs include the set of temperatures and pressures at which to estimate the phenomenological rate coefficients, and an interpolation model to be fitted to the computed k(Tn,Pm) values to give smooth functions k(T,P). All of the numerical methods used to estimate the phenomenological rate coefficients require discretization of the energy domain into a vector E. Such a discretization is characterized by a minimum energy Emin, maximum energy Emax, and either a number of grains Ngrains or a energy spacing DE. The minimum energy Emin is straightforward to determine: it is simply the ground-state energy of the lowest molecular configuration on the potential energy surface. (The ground-state energies on the potential energy surface are estimated as the enthalpy at 0 K using the thermodynamic parameters for each species.) The number of grains Ngrains or energy spacing DE depend on the user’s desired level of accuracy, and are left as parameters to be input. However, the choice of Emax is not straightforward. A choice of Emax that is too low may cause the accuracy of the resulting rates to suffer due to neglect of important energies. A choice of Emax that is too high may cause numerical effort to be wasted on energies that do not contribute significantly; worse, the master equation matrix becomes increasingly stiff due the faster reaction rates observed at higher energies. Nonetheless, an automatic procedure for determining Emax is required for our desired application of automatic mechanism generation. Our chosen procedure is as follows: 1. Guess an initial value Emax,0 that is a certain number of kBT above the highest energy on the potential energy surface, usually a transition state or a bimolecular reactant or product channel. 2. Estimate the density of states and the equilibrium distribution for the isomer with the highest ground-state energy. For a calculation performed at multiple temperatures, determine the equilibrium distribution at the maximum temperature, as this gives the widest distribution. 3. Determine the energy at which the tail of the equilibrium distribution is some desired fraction of the maximum of the distribution. 4. If a suitable energy is found, add to it the difference between the highest energy on the potential energy surface and the isomer’s ground-state energy. This is the maximum energy Emax. 5. If no energy in the current range qualifies, repeat this procedure starting from a higher value of Emax,0 until a suitable Emax is found. In the interest of speed it is preferred to minimize the number of times the density of states is estimated, so our first choice of Emax,0 is deliberately chosen to be significantly higher than the maximum energy in most cases. Once a suitable set of energy grains has been selected, the density of states r(E) is calculated for all isomers in the network. For multiple temperature and pressure calculations, this step need only be performed once. Two methods for estimating the molecular degrees of freedom needed to compute r(E), the three-frequency and functional-group frequency methods, are described in Sections 4.1 and 4.2, respectively. The input includes a set of temperatures and pressures at which to estimate the phenomenological rate coefficients. This journal is c the Owner Societies 2012 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online Fig. 1 A general procedure for estimating pressure-dependent rate coefficients for reaction networks of arbitrary size and complexity. There are certain steps that are only dependent on the current temperature and not on pressure, so we place the iteration over temperature as the outer loop and the iteration over pressure as the inner loop. For example, the equilibrium distributions are functions of temperature only and are calculated using the formula bi ðEr ; TÞ ¼ PN ri ðEr ÞebEr grains s¼1 This journal is c ri ðEs ÞebEs the Owner Societies 2012 ð9Þ for isomer i at energy Er for energy grain r, where ri(E) is the vibrational-rotational density of states for the isomer. The microcanonical rate coefficients k(E) can be computed using a variety of methods. If reactant and transition state energies, geometries, and vibrational frequencies are available (e.g. from a quantum chemistry calculation), then the RiceRamsperger-Kassel-Marcus (RRKM) method gives excellent results. However, such information is often expensive to obtain and nontrivial to automate, since it requires the sum Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 1135 View Article Online Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C of states for the transition state. Furthermore, anharmonicity effects—such as hindered internal rotations—become increasingly important for larger polyatomic molecules, and anharmonic corrections are often necessary to achieve rate coefficients at accuracies better than an order of magnitude. An alternative approach for when such information is not available is to utilize the inverse Laplace transform to transform the canonical rate k(T) in the high-pressure limit into the microcanonical rate coefficient k(E). The relationship between these, recognized by Slater,42 is k(E)r(E) = L1[Q(b)kN(b)] (10) 37 Exact formulas exist for simple Arrhenius kinetics and modified Arrhenius kinetics kN(T) = ATn exp(Ea/RT) for n 4 1/2.31 The inverse Laplace transform method assumes that the given k(T) expression is valid over the temperature range from zero to infinity, and that the activation energy Ea is physically equivalent to the reaction barrier height E0, both of which are nontrivial approximations. Note that, for simple Arrhenius kinetics (n = 0) the quantum Rice-RamspergerKassel (QRRK) method is essentially the same as the inverse Laplace transform method. The RRKM or inverse Laplace transform method is used to determine the microcanonical rate coefficient for the forward reaction only. The reverse rate coefficient is determined from detailed balance. For an isomerization reaction A - B the detailed balance equation is kf(E)rA(E) = kb(E)rB(E) (11)  Constants of proportionality in the density of states become unimportant, as they cancel when taking the ratio r(E)/Q(b), e.g. when computing the equilibrium distribution b(E,T). In Section 4 we will use this to include an arbitrary active K-rotor in the density of states expression. Thus, our implementation of the master equation uses eqn (13) and (14). The collision model is made up of the collision frequencies oi(T,P) and the collisional transfer probabilities model Pi(E,E 0 , T). As is common within the master equation literature, here oi is calculated using the Lennard-Jones collision model. This should be a reasonable choice for collisions between polyatomic molecules and weak colliders such as noble gases and diatomics, but less accurate for collisions between two large polyatomics or two very polar molecules.43 The collisional transfer probabilities function Pi(E,E 0 , T) comes in a variety of forms, all of which share the general trend that a collision is more likely to cause a small transfer of energy than a large one. The most common form of Pi(E,E 0 , T) is the single exponential down model Pi(E,E 0 , T) = Ci(E 0 , T)exp[ai(T)(E 0  E)] E o E 0 (15) where the activating equivalent is determined from detailed balance. C is determined by normalization of the probability distribution. The lone parameter aI(T) is the inverse of hDEdi. To be consistent with thermodynamics, detailed balance must be satisfied:   ri ðE 0 Þ E0  E 0 Pi ðE; E 0 ; TÞ EoE 0 exp  Pi ðE ; E; TÞ ¼ ri ðEÞ kB T ð16Þ and for a dissociation reaction A - B + C the equation is kf(E)rA(E) = kb(E)rBC(E) (12) where rBC(E) is the convolved vibrational-rotational density of states for species B and C, and also includes relative translational motion. An alternative formulation incorporates the macroscopic equilibrium coefficient Keq(T) and equilibrium distributions bi(E,T) at each temperature. For isomerization and dissociation, respectively, the alternative formulation is kf(E)bA(E,T) = Keq(T)kb(E)bB(E,T) (13) kf(E)bA(E,T) = Keq(T)kb(E)bBC(E,T) (14) where bBC(E,T) is the combined equilibrium distribution for species B and C. These two formulations are equivalent; however, there are multiple reasons to use eqn (13) and (14) instead of eqn (11) and (12):  Only the density of states of the unimolecular isomers need be computed. This is a result of the assumption of thermalized bimolecular channels, which means that we only need to compute the product kb(E) bBC(E,T), and not the individual values of kb(E) and bBC(E,T). (In eqn (4) the required products are the Fimbm terms).  Only the reactive vibrational-rotational modes need be included in the density of states. Missing modes will not affect the observed equilibrium because we are imposing the macroscopic equilibrium via Keq(T). This is particularly important in automatic mechanism generation, where we do not yet have an efficient, accurate method of estimating the external rotational constants. 1136 Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 At this point we have all of the information needed to construct the full master equation matrix, and we are ready to extract the phenomenological rate coefficients. Three methods to do this, the modified strong collision (MSC), reservoir state (RS), and chemically-significant eigenvalue (CSE) methods, are described in detail in Sections 5.1, 5.2, and 5.3, respectively. 4 Methods for estimating the density of states As described in Sections 1 and 2, formulation of the master equation requires the density of states rvk for each isomer. The relevant rovibrational degrees of freedom are those degrees of freedom in which the energy of the excited species can be randomized, i.e. the constant J and M modes are not included. In RMG, these degrees of freedom are taken to be the harmonic-oscillator (HO) vibrational frequencies, hindered internal rotors, and a one-dimensional external rotation (i.e. the K-rotor for a symmetric top approximation, which is expected to strongly couple with vibrations via Coriolis terms). Furthermore, it is assumed that the degrees of freedom are independent. RMG assumes that each molecule may be approximated as a symmetric top and that the density of states for the one-dimensional K-rotor is:  1=2 1 1 rr ðEÞ ¼ ð17Þ  s jA  BjE where s is the rotational symmetry number for this rotor, A is the unique rotational constant, and B% is the average of the two This journal is c the Owner Societies 2012 View Article Online Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C most similar rotational constants. As stated in Section 3, since % is a proportionality constant in the above expression, s|A  B| we need not be concerned with its actual value, as it will cancel when we evaluate equilibrium distributions. A discussion of the accuracy of the K-rotor model is available from Zhu et al.44 The density of states for hindered internal rotors is given by the Laplace transform of the hindered-rotor partition function. RMG assumes a semi-classical approximation for the hindered-rotor partition function: QHRsemi ¼ QHRclassical ¼  V0 p 2s2 Br QHOquantum QHOclassical 1=2 eau 2a e I0 ðuÞ 1  e2au 1=2 u u ð18Þ where Br is the rotational constant corresponding to the reduced moment of inertia, V0 is the rotational barrier height assuming a single cosine potential, s is the symmetry number of the internal rotor, n is the normal mode frequency corresponding to the torsion, I0 is the modified Bessel function of the first kind, u  V0/2kBT, and a  n/V0. Like with the K-rotor, the rotational constant Br appears only as a proportionality constant, so its value will not be important for our master equation calculation. If detailed information about the species or transition state is available, we compute the reduced moment of inertia using the method of Pitzer45 which is described as I(2,3) by East and Radom.46 The density of states corresponding to the classical hindered rotor partition function was described by Knyazev.47 Unfortunately there is no analytical expression for the density of states corresponding to the semi-classical hindered rotor partition function, eqn (18), since there is no analytical inverse Laplace transform for the denominator in the quantum harmonic oscillator partition function. However, the BeyerSwinehart method gives the exact density of states of the quantum harmonic oscillator. By replacing this denominator with a truncated Taylor expansion, the result from Knyazev can be extended to include the density of states corresponding to a semi-classical hindered rotor partition function: 8 pffiffiffiffiffiffiffiffiffi ~ 0 2K E=V > < pffiffiffiffiffiffiffiffi ps Br V0  ffi rHRsemi ðEÞ ¼ 2K pffiffiffiffiffiffiffiffi V0 =E~ > : pffiffiffiffiffiffi ps Br E~ for Eo max ðV0 ; hn=2Þ otherwise ð19Þ Above, Ẽ  E  hn/2, and K is the complete elliptical integral of the first kind. The density of states eqn (19) is calculated for each hindered rotor. The combined density states for multiple hindered internal rotors is calculated numerically by repeated application of the convolution integral. Similarly, the density of states for internal rotors and active external rotation is calculated by convoluting the density of states for the active K-rotor with the density of states for all the hindered rotors. Finally, following the method of Astholz,20,34 the BeyerSwinehart algorithm18 is initialized with the K-rotor/ hindered-rotor density of states, which convolutes this density of states with the vibrational modes. The result is the complete rovibrational density of states, r(E). This journal is c the Owner Societies 2012 Accurate determination of the vibrational frequencies and the potentials for internal rotation requires a threedimensional geometry and a quantum calculation, which is not routinely available in the RMG framework. Instead, these parameters are estimated from the heat capacity. The heat capacity may be separated into contributions from external translation, external rotation, and the internal rovibrational degrees of freedom: Ctotal = Ctranslation + Crotation + Crovibration V V V V In this presentation, the K-rotor is included in the external rotation. The vibrational frequencies and hindered rotor parameters can be calculated by fitting a function to the vibrational contribution to the heat capacity. In most cases, however, the number of vibrational modes exceeds number of known values of the heat capacity as a function of temperature, and therefore these frequencies cannot be determined uniquely. Two methods for approximating the vibrational frequencies are described below. 4.1 The three frequency model The three frequency model for estimating r(E) from heat capacity data was proposed by Bozzelli, Chang, and Dean.33 This approach assumes that the molecule can be described by the HO model. The heat capacity of a single mode as a function of its HO frequency is given by the following formula from statistical mechanics:  2 CVHO ðT; n i Þ n i =kB T ¼ eni =kB T n =k T ð20Þ R e i B 1 where ni is the vibrational frequency for the ith mode. Bozzelli et al. fit three vibrational frequencies to the heat capacity at seven temperatures (300, 400, 500, 600, 800, 1000, 1500 K), using the following formula: CVvibration ðTÞ ¼ 3 X i¼1 gi CVHO ðT; n i Þ ð21Þ g3 ¼ s  g1  g2 where s is the number of vibrational modes (e.g. 3N  6 for a nonlinear molecular with no internal rotors). In the original implementation, the three degeneracies gi were allowed to be non-integers. Later implementations of this approach forced the degeneracies to be integers, so r(E) could be computed efficiently using the Beyer-Swinehart algorithm, and that is what we have done here.27 Thus, five parameters are fit to the heat capacity data: three frequencies and two (integer) degeneracies. This method is simple and computationally efficient. However, there are two problems with this approach. First, the resulting pseudo-frequencies may not be physically realistic. Second, it assumes that the heat capacity necessarily increases monotonically with temperature. For molecules with hindered internal rotors, the true heat capacity may actually have a local maximum with respect to temperature. 4.2 The functional-group frequency model An alternative to the three-frequency model, inspired by infrared spectroscopy, is to assign frequencies to functional groups. Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 1137 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online For each functional group, the number of characteristic frequencies is limited to 3Natoms  Nrotons  6. For example, a heavy atom attached to a methyl group contains 5 atoms and one single bond between heavy atoms, and thus has 9 characteristic frequencies: 3 C–H stretches, 1 R–C stretch, 2 R–C–H bending modes, 2 R–C–H rocking modes, and 1 umbrella mode. Additionally, one of the six subtracted modes is the overall rotation of the methyl group, and this torsional mode is treated as a hindered internal rotor. RMG automatically recognizes that single bonds between heavy atoms that are not within a ring structure may be hindered rotors. Each vibrational mode is associated with a range of frequencies, and RMG automatically selects the appropriate number of frequencies, taking them to be evenly spaced between the two bounds. Thus, for a molecule with two methyl groups, RMG would assume 6 C–H stretching frequencies that are evenly spaced between 2750–2850 cm1. A list of the functional groups, the corresponding range of frequencies, and the number of frequencies per group are available in the Supporting Materials.w By summing over the functional groups, the model predicts the majority of the HO vibrational frequencies. The heat capacity for the remaining degrees of freedom (e.g. vibrational frequencies not explicitly accounted for, as well as hindered rotors) is obtained by subtracting the heat capacity of the functional group frequencies from the rovibrational heat capacity: CVremaining degrees of freedom ðTÞ ¼ CVrovibration ðTÞ Nfunctional group frequencies X  i¼1 CVHO ðT; n i Þ ð22Þ The remaining vibrational frequencies and the parameters for hindered internal rotors are determined from this heat capacity data. Unlike the three frequency model, this method includes a function for hindered internal rotors and therefore can accommodate heat capacities with a local maximum. Nremaining CVremaining degrees of freedom X ¼ i¼1 þ NX rotors j¼1 CVHO ðT; n i Þ CVHR ðT; V0;j ; n j Þ where Nremaining  3N  6  Nrotors  Nfunctional group frequencies. The heat capacity of a hindered internal rotor is derived from the semi-classical model for the hindered rotor partition function, eqn (18). The heat capacity for this partition function is:   d ln QHRsemi CVHR ðT; V0 ; nÞ d ¼ T2 dT R dT d2 ln QHRsemi du2 "   # 1 I1 ðuÞ I1 ðuÞ 2 u   þu u 2 I0 ðuÞ I0 ðuÞ ¼ u2 þ 1138 2uaeau 1  e2au  2 Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 ð23Þ Although the functional group method will estimate most of the frequencies, it may be the case that not all of the remaining parameters can be determined uniquely from the heat capacity data. Typically, the CV(T) data from a group additivity calculation has enough information content to determine about six parameters. If there are more than six unknown parameters, then some of the rigid-rotor harmonic-oscillator frequencies or the hindered rotor parameters (or in some cases both) may be lumped into a single frequency or hindered rotor, respectively, so that no more than six parameters are fitted: N X i¼1 M X j¼1 CVHO ðT; n i Þ ! NCVHO ðT; nÞ CVHR ðT; V0;j ; n j Þ ! MCVHR ðT; V0 ; nÞ Alternatively, one can sometimes obtain a better fit and/or fit additional parameters directly by first fitting an interpolation model to the known CV(T) points, and choosing a higher density of points to use when performing the frequency fitting. We use this approach with the acetyl + oxygen case study discussed in Section 6, using the Wilhoit polynomial model48,49 to directly fit up to ten parameters. In order to make the resulting parameters as realistic as possible, the fitted parameters are given physically meaningful constraints: the harmonic oscillator frequencies are bound between 180 and 4000 cm1, the hindered rotor torsional frequencies are bound between 40 and 600 cm1, and the hindered rotor barrier heights are bound between 10 and 10 000 cm1. The FORTRAN 90 code DQED50,51 was used for the bounded, constrained, non-linear least-squares fitting procedure. 5 Methods for estimating phenomenological rate coefficients The objective of each of the methods described in the following sections is to transform the large master equation matrix from eqn (4) into a small number of phenomenological rate coefficients k(T,P). All of the methods share a common formalism in that they seek to express the population distribution vector pi for each unimolecular isomer i as a linear combination of the total populations xj(t) and ymA(t) ymB of unimolecular isomers Cj and bimolecular reactant channels Am + Bm: pi ðtÞ ¼ N isom X j¼1 xj ðtÞuij þ N reac X m¼1 ymA ðtÞymB vim ð24Þ The vectors uij represent the portion of the population distribution of energy states of unimolecular isomer i that tracks the population of isomer j. In the modified strong collision and reservoir state methods, this is because the energy levels of isomer i are in pseudo-steady-state relationships with isomer j. The interpretation is a bit different in the chemically-significant eigenvalues approach, but the form of the equations is the same. Similarly, the vectors This journal is c the Owner Societies 2012 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online vim represent the population distribution of energy states of unimolecular isomer i that tracks the population of bimolecular reactant channel m. Both uij and vim are independent of time. Notice that, even for the best uij and vim, eqn (24) is an approximation to the true pi(t) given by the solution of the full master equation. This approximation is a form of model reduction, and is done so that we can use the conventional phenomenological kinetics equations, and do not have to keep track of the evolving energy distributions of each species. The phenomenological rate coefficient matrix can be constructed from the full master equation matrix M shown in eqn (4) and the population distribution vectors uij and vim: Ngrains kij ðT; PÞ ¼ X s¼1 ðMi uij Þs þ N grains isom NX X ‘ai Ngrains kim ðT; PÞ ¼ X s¼1 ðMi vim Þs þ s¼1 ðKi‘ u‘j Þs N grains isom NX X ‘ai s¼1 ð26aÞ ðKi‘ v‘m Þs ð25bÞ Ngrains þ X s¼1 ð25aÞ knm ðT; PÞ ¼ þ N isom X ‘¼1 N isom X ‘¼1 gn‘  u‘j ð25cÞ gn‘  v‘m n¼1 ð26bÞ ynA ðtÞynB wn ðEÞ 31 As described by Chang, Bozzelli, and Dean, the modified strong collision method utilizes a greatly simplified collision model that allows for a decoupling of the energy grains. In the simplified collision model, collisional stabilization of a reactive isomer is treated as a single step process, ignoring the effects of collisional energy redistribution within the reactive energy space. An attempt to correct for the effect of collisional energy redistribution is made by modifying the collision frequency oi(T,P) with a collision efficiency bi(T), following early work by Troe on the single-isomer ‘‘fall-off’’ case.30,52 By approximating the reactive populations as existing in pseudo-steady state, a matrix equation is formed at each energy through which uij and vim vectors, and subsequently the k(T,P) values, can be determined. the Owner Societies 2012 mi ðEÞ ¼  oi ðT; PÞbi ðTÞ  ð25dÞ The modified strong collision method c N reac X where xi(t) reflects the total population of the nonreactive energies of isomer i, dij is the Kronecker delta, and Above, the indices i and j represent unimolecular isomers of the initial adduct, m represents bimolecular reactants, n represents bimolecular reactants and products, and s represents an energy grain. Thus, the rate coefficients above are for isomerization, association, dissociation, and bimolecular reactions Am + Bm - An + Bn, respectively. The various methods provide different values of uij and vim. This journal is N isom X d pðE; tÞ ¼ LðEÞpðE; tÞ þ xi ðtÞzi ðEÞ dt i¼1 ðFim bm Þs knj ðT; PÞ ¼ 5.1 After applying the modified strong collision approximation, a population balance performed at a reactive energy E gives 2 2 3 32 3 p1 ðE; tÞ m1 ðEÞ k12 ðEÞ . . . p1 ðE; tÞ 6 6 7 76 7 6 7 76 7 d6 6 p2 ðE; tÞ 7 ¼ 6 k21 ðEÞ m2 ðEÞ . . . 76 p2 ðE; tÞ 7 6 6 6 7 7 7 dt 4 4 5 54 5 .. .. .. .. .. . . . . . 2 3 oi ðT; PÞbi ðTÞbi ðE; TÞdi1 6 7 N isom X 6 o ðT; PÞb ðTÞb ðE; TÞd 7 i i i2 7 i xi ðtÞ6 þ 6 7 4 5 i¼1 .. . 2 3 f1n ðEÞbn ðE; TÞ 6 7 N reac X 6 f ðEÞb ðE; TÞ 7 7 2n n ynA ðtÞynB 6 þ 6 7 4 5 n¼1 .. .  Nreac þNprod X n¼1 N isom X jai kji ðEÞ ð27Þ gni ðEÞ is the total rate of loss from isomer i. Applying the pseudosteady state approximation to the above gives LðEÞpðEÞ ¼  N isom X i¼1 xi ðtÞzi ðEÞ  N reac X n¼1 ynA ðtÞynB wn ðEÞ ð28Þ Solving these gives p(i)(E) = [L(E)]1zi(E) p(n)(E) = [L(E)]1wn(E) (29) where p(i)(E) is the pseudo-steady population of each isomer at reactive energy E resulting from the thermal activation of isomer i, and p(n)(E) is the pseudo-steady population of each isomer at reactive energy E resulting from the chemical activation of reactant channel n. These can then be used to compute the k(T,P) values: Ngrains kij ðT; PÞ ¼ oi ðT; PÞbi ðTÞ knj ðT; PÞ ¼ N grains isom NX X i¼1 r¼r0;i X r¼r0;i ðjÞ pi ðEr Þ ðjÞ gni ðEÞpi ðEr Þ ð30aÞ ð30bÞ Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 1139 View Article Online Ngrains kim ðT; PÞ ¼ oi ðT; PÞbi ðTÞ Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C knm ðT; PÞ ¼ N grains isom NX X j¼1 r¼r0;j X r¼r0;i ðmÞ pi ðEr Þ ð30cÞ ðmÞ ð30dÞ gnj ðEr Þpj ðEr Þ The indices above have the same meaning as in eqn (25): i and j represent unimolecular isomers, m represents bimolecular reactants, n represents bimolecular reactants and products, and r represents energy grains. The modified strong collision method can be expressed in terms of the notation presented in the introduction to this section, eqn (24):  bi ðEr Þdij Er oEj;crit uij ðEr Þ ¼ ð31aÞ ðjÞ pi ðEr Þ Er  Ej;crit vim ðEr Þ ¼  0 Er oEm;crit ðmÞ pi ðEr Þ Er  Em;crit ð31bÞ where dij is the Kronecker delta. The critical energy that divides the low-energy Boltzmann grains from the high-energy steady-state grains is determined as the first reactive energy grain, i.e. the lowest energy where any microcanonical rate coefficient k(E) is nonzero. Inserting the above equations into eqn (25) gives identical k(T,P) values as those in eqn (30c), (30d), (30a) and (30b). 5.2 The reservoir state method An alternative to approximating collisional stabilization as a single-step process, proposed by Green and Bhatti,32 is to assume that, except for a depleted region near the transition state energy, the low energy grains are Boltzmann distributed. Therefore, we partition the energy grains for each isomer; the low-energy grains form the reservoir, while the high-energy grains form the active-state. Green and Bhatti suggest the cutoff be placed a few kBT below the lowest transition-state energy adjacent to that isomer; our work suggests that a better choice is simply to place the cutoff at the lowest adjacent transition-state energy, as the intervening grains are better modeled as Boltzmann distributed than as pseudo-steady. (This choice only matters at temperatures where the equilibrium population of these grains is significant). Here we neglect all reactions (including tunneling) from these reservoir grains; if tunneling below the transition state energy is important, one should set the reservoir cutoff a little lower. As the development below implies, the mathematics of the reservoir state method parallels that of the modified strong collision method. The primary difference is that, for the reservoir state method, we will be working with all high-energy grains for all isomers at once, rather than being able to treat each energy grain independently. A more detailed treatment, including a demonstration that the resulting k(T,P) values are in fact independent of the initial condition, is provided in the Supporting Materials.w First, a Boltzmann distribution bri is imposed on the reservoir grains pri for each isomer i, i.e. pri 1140 = xi(t)bri Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 (32) leaving a single time-dependent constant of proportionality xi(t) which is related to the total population of isomer i. Note that a significant fraction of the equilibrium population is in PNres;i r the activated grains at high temperatures, i.e. s¼1 ðbi Þs o1. r The quantity bi is not renormalized, however, in order to ensure that the computed k(T,P) values satisfy macroscopic equilibrium. With the reservoir approximation, the master equation for the active-state grains pai for each isomer is 2 aa 2 a3 32 a 3 p1 p1 Ka12 . . . M1 6 6 7 76 7 6 a 76 7 7 aa d6 6 pa2 7 ¼ 6 K21 M2 . . . 76 pa2 7 6 76 7 7 dt 6 4 54 5 4 5 .. .. .. .. .. . . . . . 2 ar r 3 M1 b1 d1j 6 7 N isom X 6 Mar br d 7 2 2 2j 7 xj ðtÞ6 ð33aÞ þ 6 7 4 5 j¼1 .. . 2 a a 3 F1m bm 6 7 N reac X 6 Fa ba 7 6 7 2m m ymA ðtÞymB 6 þ 7 4 5 m¼1 .. . N N isom reac X X d a xj ðtÞzj þ ymA ðtÞymB wm p ¼ Lpa þ dt j¼1 m¼1 ð33bÞ As with the modified strong collision method, the pseudosteady state approximation is applied to the active-state grains to give a matrix equation Lpa ¼  N isom X j¼1 xj ðtÞzj  N reac X m¼1 ymA ðtÞymB wm ð34Þ for which the solution is pa(j) = L1zj pa(m) = L1wm (35) a(i) where p is the pseudo-steady active-state population of each isomer resulting from the thermal activation of isomer i, and pa(n) is the pseudo-steady active-state population of each isomer resulting from the chemical activation of reactant channel n. These can then be used to compute the k(T,P) values via X aðjÞ kij ðT; PÞ ¼ ðMra ð36aÞ i pi Þr r knj ðT; PÞ ¼ kim ðT; PÞ ¼ N isom X i¼1 X knm ðT; PÞ ¼ r aðjÞ ð36bÞ aðmÞ ð36cÞ gni  pi ðMra i pi N isom X i¼1 This journal is Þr aðmÞ gni  pi c ð36dÞ the Owner Societies 2012 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online where pa(j) and pia(m) contain the components of pa(j) and pa(m) i relating to isomer i. Once again, indices i and j represent unimolecular isomers, index m represents bimolecular reactants, index n represents bimolecular reactants and products, and index r represents energy grains. Note that we have not needed the common pseudo-first-order assumption on bimolecular reactants (i.e. ymA { ymB) to arrive at the above result. However, the steady-state approximations used in both the modified strong collision and reservoir state methods implicitly assume that both ymA and ymB are small enough that ymA ymBFimbm is small. The terms of the pseudo-steady state vector pai and the reservoir populations bri can be used to construct the uij and vim vectors for use in eqn (24): uij ¼ dij brj aðjÞ pi vim ¼ 0 aðmÞ pi ð37Þ Again, using the above expressions with eqn (25) gives identical k(T,P) values as those from eqn (36). Finally, a note on efficient solving of eqn (34). Collisional energy transfer models favor small transfers of energy over large ones. For the single exponential down model, P(E,E 0 ) 0 exponentially as |E 0  E| - N. This suggests that the probabilities are negligible for sufficiently large energy transfers, which in turn results in a matrix Pi that is strongly banded. We can take advantage of this bandedness by using an indexing scheme that iterates over energies as the outer loop and isomers as the inner loop. The equations in this manuscript utilize the opposite indexing scheme for ease of understanding, but our software implementation of these equations uses the more efficient indexing. A typical sparsity pattern for an active-state matrix L using the two indexing schemes is shown in Fig. 2 for a three-isomer network. The grainmajor indexing scheme, Fig. 2b, yields a matrix with a much narrower bandwidth than the isomer-major indexing scheme, Fig. 2a. 5.3 The chemically-significant eigenvalues method The following description of the chemically-significant eigenvalues method is based on the works of Pilling and Robertson53,54 and Miller and Klippenstein.22,55,56 First, the eigenvalues and eigenvectors of the master equation matrix in eqn (4), denoted here as M, are computed: M = QLQ1 (38) The eigenvalue decomposition is usually performed after symmetrizing M, which ensures that the calculated eigenvalues K are real. Since this is a physical system, the eigenvalues are nonpositive, with one eigenvalue of exactly zero corresponding to the equilibrium distribution. (If irreversible product channels are included, however, the zero eigenvalue will not be present). If the system contains Nchem isomers and reactant channels, then we expect the master equation matrix to consist of Nchem chemically-significant eigenvalues (including the zero eigenvalue, if present). Under conditions of high pressure and low temperature, these eigenvalues should be distinct and lower in magnitude (i.e. less negative) than the other eigenvalues, which This journal is c the Owner Societies 2012 correspond to internal energy relaxation. The distinctness criterion is especially important: if the largest magnitude chemical eigenvalue is too close to the smallest magnitude internal energy eigenvalue, the calculated k(T,P) values will be nonsensical. We only consider the case where Nchem chemically-significant eigenmodes can be identified. The separation of the chemical and internal energy eigenvalues means that we can select a time after the initial energy modes have relaxed but before significant chemical transitions have occurred. After this short-time period, the trajectory of the solution will be contained within the manifold defined by the chemically-significant eigenvectors. Therefore we ought to be able to construct the k(T,P) values using only the chemically-significant eigenpairs, independent of the initial condition. First we construct an Nchem  Nchem matrix Z from the chemically-significant eigenvectors by summing all terms in each eigenvector for each isomer: X Qs‘ i ¼ 1; 2; . . . ; Nisom ; ‘ ¼ 1; 2; . . . ; Nchem ð39Þ Zi‘ ¼ s2i For reactant channels there is only one element in each eigenvector for that channel, so it is transferred to Z as-is: Znl = Qnl n = 1,2,. . .,Nreac, l = 1,2,. . .,Nchem (40) Physically the matrix Z corresponds to the eigenvector matrix obtained from diagonalizing the matrix of phenomenological rate coefficients K, i.e. K = ZL 0 Z1. For systems satisfying the assumed separation of timescales, this allows us to directly express the phenomenological rate coefficients for all reactions between isomers and reactant channels: kij ðT; PÞ ¼ NX chem ‘¼1 l‘ Zi‘ Z‘j1 i; j 2 1; 2; . . . ; Nisom þ Nreac ð41Þ To obtain k(T,P) values for dissociation and bimolecular reactions to a product channel (infinite sink), we will need an additional Nprod  Nchem matrix Y with elements Ym‘ ¼ N isom X i¼1 ðiÞ gmi  q‘ m ¼ 1; 2; . . . ; Nprod ð42Þ Above, ql(i) contains only the elements of chemically-significant eigenvector ql involving isomer i. The k(T,P) values for reactions to the product channels are then kmj ðT; PÞ ¼ Ncse X ‘¼1 Ym‘ Z‘j1 j ¼ 1; 2; . . . ; Nisom þ Nreac ð43Þ The time-independent population distribution vector sets uij and vim can be constructed by combining the full and reduced chemically-significant eigenvectors via uij(Er) = (QchemZ1)(i,r),j vim(Er) = (QchemZ1)(i,r),m (44) where Qchem is a rectangular matrix containing only the chemically-significant eigenvalues, and the subscripts indicate the row of QchemZ1 corresponding to isomer i at energy grain r, and the column corresponding to isomer j or reactant channel m. Note that, since the original master equation Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 1141 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online Fig. 2 A typical sparsity pattern for an active-state matrix of a three-isomer network using (a) an isomer-major indexing scheme and (b) an energy grain-major indexing scheme. The latter gives a significantly more banded matrix, which can be used to accelerate the reservoir state linear solve. matrix M depends on ymB, in the eigenvalue method the k(T,P), uij, and vim depend on ymB. However, in most cases only small ymB are of interest, for which these dependencies are very weak. The only new assumption utilized by the chemicallysignificant eigenvalues method is that the chemical and internalenergy relaxation timescales are distinct. When this is the case, as it often is at low and moderate temperatures and high and moderate pressures, the k(T,P) values will accurately reflect the system behavior at chemical timescales. For high temperatures and/or low pressures, one or more chemical timescales become indistinguishable from the internal-energy relaxation timescales. To apply the chemically-significant eigenvalues method in this case, it is necessary to lump the configurations that participate in the fast chemical relaxations together into a single pseudospecies. To our knowledge, no one has yet automated this process; although it certainly would be possible to do so, it is beyond the scope of this paper. Finally, Miller and Klippenstein discuss two ways of extracting the k(T,P) values from the chemically-significant eigenpairs. The above discussion presents the so-called ‘‘long-time’’ method, which has two advantages: it numerically performs better as a chemical eigenvalue approaches the continuum of internal energy eigenvalues, and it makes clear that the k(T,P) values are not dependent on the initial populations. The other method, called the ‘‘initial-rate’’ method, gives identical results at almost all conditions, and may be slightly faster to execute since it avoids the small matrix inversion required by the long-time method. A brief discussion of the initial-rate method is provided in the Supporting Materials.w 6 Case study: the acetyl + oxygen system As an application of the various methods for estimating the density of states and the phenomenological rate coefficients, we turn our attention to the chemically-activated reaction of acetyl radical with oxygen. This system is important to atmospheric chemistry as a step in the conversion of acetaldehyde to peroxyacetylnitrate (PAN), the latter of which is a 1142 Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 secondary air pollutant.57 It is also potentially important in the ignition chemistry of ethanol. The potential energy surface for this reaction was computed originally by Lee, Chen, and Bozzelli.58 The values in Fig. 3 are updated results, described below. Our updated values are consistent with those recently published by Villano et al.59 Michael, Keil, and Klemm demonstrated that the acetyl radical is a major product of the reaction of acetaldehyde with OH radical, which preferentially abstracts the weakly-bonded aldehydic hydrogen. They also observed pressure-dependent regeneration of OH radical when acetyl radical is reacted with oxygen, ranging from nearly complete regeneration at low pressures to minimal regeneration at high pressures.60 Later experimental efforts generally agree on the relevant pressure range being about 0.001 to 1 bar at 300 K.61–65 Recent experiments have confirmed that the a-lactone is the carboncontaining product associated with OH regeneration.66 On the theoretical side, this system was previously studied by Lee, Chen, and Bozzelli using a three-frequency model to estimate the density of states and the modified strong collision method to estimate the phenomenological rate coefficients. Their calculations significantly underpredicted the regeneration of the OH radical.58 More recently Maranzana, Barker, and Tonachini studied the system using quantum chemistry calculated modes as input to the density of states and the MultiWell stochastic master equation solver. They allowed the collisional energy transfer parameters to vary so as to match the experimental data, and thus were able to match the regeneration of OH.67 The calculations of Lee et al. predicted that the dominant bimolecular product channel was a biradical + OH. Hou et al. highlight the unlikeliness of this product after recomputing the potential energy surface using the G3MP2 compound method.63 Based on their results, we do not consider the biradical as a stable configuration, and do not include it in our potential energy surface or subsequent calculations. Additionally, Hou et al. found a saddle point geometry for the association reaction and used a rigid-rotor harmonic oscillator model to compute the association rate coefficient. Depending upon which method they chose—MP2, CCSD, G3MP2, CASPT2—the barrier height This journal is c the Owner Societies 2012 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online Fig. 3 The RQCISD(T)/CBS//B3LYP/6-311++G(d,p) potential energy surface for the acetyl oxidation network. varied between +13 kcal mol1 and 1 kcal mol1, relative to the reactants. Given that the C–O bond distance in the saddle point was 2.5 Å, the CASPT2(7,7) DE0 of 1.7 kcal mol1 is the most reliable. At these distances the modes between the two reacting moieties will be highly anharmonic, and a rigid-rotor harmonic oscillator model is ill-advised. Using DFT, we were able to locate a saddle point for this reaction as well, but there were four vibrational modes with frequencies less than 20 cm1. These coupled, anharmonic modes correlate to the relative orientation and separation of the two fragments. A better way to treat this entrance channel would be to use variable reaction coordinate transition state theory.68–70 Ideally the interaction potential would be computed at the CASPT2 level, with an active space of at least 7 electrons and five orbitals—six electrons in four orbitals for O2, and one electron in one orbital for the unpaired electron in acetyl, with the possibility of adding additional electrons and orbitals to account for the CQO p-bond and/or any oxygen lone pair electrons in acetyl. However, these calculations are beyond the scope of this work, and we have adopted the high-pressure limit for acetyl + O2 recommended by Lee et al. for simplicity. We performed our own RQCISD(T)/CBS//B3LYP/ 6-311++G(d,p) calculations to generate the PES in Fig. 3 Table 1 Calculated high-pressure limit rate coefficients for acetyl + O2a Reaction A b Acetyl + O2 - acetylperoxy Acetylperoxy - hydroperoxyl-vinoxy Acetylperoxy - ketene + HO2 Hydroperoxyl-vinoxy - ketene + HO2 Hydroperoxyl-vinoxy - lactone + OH 4.4 2.3 2.6 5.3 1.9      12 10 109 109 1016 1017 n Ea 0.0 0.75 1.2 1.0 1.1 0.0 23.2 34.1 29.5 27.2 The units for A are s1 for unimolecular reactions and in cm3 molecule1 s1 for bimolecular reactions, with Ea in kcal mol1. The rate coefficient is k = A(T/1[K])nexp(Ea/RT). Computed from RQCISD(T) variational TST calculations; see text for details. b Taken from ref. 53. a This journal is c the Owner Societies 2012 using the methodology of Goldsmith et al.71 All DFT and CBS-QB3 calculations were done using Gaussian03.72 All MP2 and QCISD(T) calculations were done using MOLPRO.73 The high-pressure limit for each reaction rate coefficient except acetyl + O2 was computed using variational transition state theory, as implemented in Variflex.25 These results are provided in Table 1. Although these are high-level calculations, it is doubtful that these numbers are the final word on this important reaction system. 6.1 Density of states comparison Fig. 4 shows the density of states, partition function, and heat capacity for acetylperoxy and hydroperoxylvinoxy as calculated from vibrational and hindered rotor data estimated using the three-frequency and group frequency approaches and calculated from quantum chemistry, with the goal of judging how the two approximate models compare to the detailed quantum chemistry model. The parameters for the internal modes as determined for each model are given in tables in the Supporting Materials.w The quantum chemistry parameters were obtained from geometry optimization and frequency calculations using density functional theory with the B3LYP functional and 6-311++G(d,p) basis set, as described in the preceding section. The three-frequency and functional groupfrequency methods estimate the internal (vibrational and hindered rotor) modes, and a K-rotor for the coupled external rotational mode. For the purposes of this comparison, we do the same with the quantum chemistry calculations. Both acetylperoxy and hydroperoxylvinoxy have 18 internal degrees of freedom to be represented by the three pseudofrequency-degeneracy pairs. Acetylperoxy contains 16 vibrational modes and two hindered rotor modes, while hydroperoxylvinoxy contains 15 vibrations and three hindered rotors. Nonetheless, the three frequency method will represent all of these as three pseudofrequency-degeneracy pairs. As shown in Fig. 4, we see that this method performs somewhat Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 1143 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online Fig. 4 Comparison of densities of states, partition functions, and heat capacities for internal degrees of freedom for the (a) acetylperoxy and (b) hydroperoxylvinoxy species using a three-frequency model (TF), group frequency model (GF), and quantum chemistry calculations (QC). Both the TF and GF models reproduce the QC heat capacity data, but GF more consistently predicts the density of states than TF. better for hydroperoxylvinoxy than for acetylperoxy. Our quantum chemistry-calculated barrier heights for the acetylperoxy hindered rotors are 1.22 and 6.12 kcal mol1, while those for hydroperoxylvinoxy are 4.78, 6.26, and 14.7 kcal mol1. The high barriers to internal rotation make the harmonic oscillator approximation fairly accurate for hydroperoxylvinoxy. For acetylperoxy, the group frequency approach provides twelve of the eighteen internal modes as characteristic vibrational frequencies: nine from the CH3 group and three from the COO group. For hydroperoxylvinoxy, eleven of the eighteen internal modes are set via characteristic frequencies: six from the CH2 group and five from the COOH group. In both cases the remaining modes—four vibrations and two hindered rotors for acetylperoxy, four vibrations and three 1144 Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 hindered rotors for hydroperoxylvinoxy—are few enough that their parameters can be fit directly to the remaining heat capacity data. The results in Fig. 4 show an improved estimate of the density of states for both isomers when compared to the three frequency method. In particular, the fitted hindered rotor barrier heights for acetylperoxy using the group frequency method are 0.99 and 4.68 kcal mol1, correctly predicting a low-barrier and a moderate-barrier hindered rotation. 6.2 Comparison of the phenomenological rate coefficient methods There are several means by which the three methods of extracting phenomenological rate coefficients from the master This journal is c the Owner Societies 2012 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online equation may be evaluated. Most pertinent is direct comparison of the k(T,P) values, as these are what we wish to use in our kinetic models. The downside of this approach is that there is no ‘‘true’’ set of k(T,P) values available to use as a basis for comparison, since the solution of the master equation does not exactly match the conventional reduced-order k(T,P) model. To enable comparison with the true solution, we must turn to comparison of concentration profiles, which we can obtain by direct integration of the master equation. However, we can go one step further by using eqn (24) to reconstruct the approximate population distributions pi(E,t) for each isomer, and compare those with one another and with the full master equation solution. In this section we apply these methods of comparison to the acetyl + oxygen reaction network. 6.2.1 Rate coefficients. The phenomenological rate coefficients for the acetyl + oxygen system were calculated by the modified strong collision, reservoir state, and chemicallysignificant eigenvalue methods using the same r(E) and k(E) from quantum chemical calculations, the same energy discretization (constant 0.5 kcal mol1 grain spacing), and hDEdowni = 200(T/300)0.85 cm1.74 Fig. 5 shows the k(T,P) values for CH3CO + O2 - products at a variety of temperatures and pressures. (The equivalent values for acetylperoxy products and hydroperoxylvinoxy - products are available in the Supporting Materials.w) As expected, at high pressures and low temperatures, the dominant pathway is collisional stabilization to the acetylperoxy isomer. At low pressures, the lactone + OH product channel is dominant at low and moderate temperatures, while the ketene + HO2 channel is dominant at high temperature. All of this is consistent with our physical understanding of the system. The modified strong collision, reservoir state, and chemicallysignificant eigenvalue methods give very similar values of the rate coefficients at high and moderate pressures and at low and moderate temperatures. However, there are several regions of disagreement. The effects of temperature and pressure on the agreement varies with the type of product. For reaction to the isomers acetylperoxy and hydroperoxylvinoxy, there is good agreement between the methods at low temperature and high pressure; increasing temperature and/or decreasing pressure causes the methods to diverge somewhat. For reactions to the bimolecular product channels, the reverse effect is observed: all methods disagree somewhat at low temperature and high pressure, while increasing temperature and/or decreasing pressure causes the methods to converge. Note that the methods always seem to show consensus on the rate coefficients for the major pathways, and only disagree on the minor channels. However, depending on what aspect of the system is of physical interest, having a good value for the minor rate constants can be important. 6.2.2 Concentration profiles. The most direct method of evaluating the performance of the three methods is to use the k(T,P) values to generate predicted concentration profiles for each approximate method and compare them to those produced by solving the full time-dependent master equation with a stiff ODE solver. The objective here is to evaluate how well the xi(t) and ymA(t) reflect the total populations of isomer i This journal is c the Owner Societies 2012 and reactants/products m. For these calculations, the initial concentrations were set to a small amount of acetyl in an atmospheric ratio of oxygen and nitrogen such that the total was consistent with the temperature and pressure chosen. The lactone + OH and ketene + HO2 product channels were treated as irreversible sinks. Both low-temperature (500 K) and high-temperature (1500 K) conditions were studied; in each case the pressure was chosen to be as low as possible without obtaining unphysical results with the chemicallysignificant eigenvalue calculation. Fig. 6 compares the concentration profiles from each method with those from the full solution at the same sets of conditions. Under all conditions shown, all three methods exhibit a period at short times before their profiles match those from the full master equation solution. This time period reflects the time scale of collisions, and is generally of the order of 10o1, where o is the mean collision frequency at the conditions being studied. Once a suitable number of collisions have occurred, all of the methods predict the total concentrations of reactants, products, and isomers very well. For the modified strong collision and reservoir state methods, the poor performance at short times comes from use of the pseudosteady state approximation, while for the chemically-significant eigenvalues method, this comes from treating the internal energy eigenmodes as completely relaxed. The first set of conditions, 500 K and 0.1 bar, represent a low temperature, moderately low pressure system. Under these conditions acetylperoxy is the primary product, although the pressure is low enough that a significant amount of lactone is also produced. All of the methods approximate the total concentrations of reactants, products, and isomers very well after the short-time period of about 108 s. The modified strong collision method does deviate from the others slightly (most notably for hydroperoxylvinoxy), but remains within an order of magnitude of the true profiles. Also interesting is the effect of high temperature, demonstrated using conditions of 1500 K and 1 bar. This time the primary product is the ketene channel, with relatively little of each isomer produced along the way. All of the reactant and product channels are well-approximated by all of the methods, but show some disagreement for the total concentrations of the isomers. Of the three methods, the chemically-significant eigenvalues method is clearly the most accurate in predicting the isomer profiles, while the modified strong collision and reservoir state methods are only slightly worse. 6.2.3 Population distributions. The concentration profile comparison is valuable because it allows for comparison of the three methods with the ‘‘true’’ solution obtained from the full master equation. However, the ‘‘true’’ solution provides the full population distributions pi(E,t) for each isomer, while the approximate solutions only provide the total populations xi(t). In comparing concentration profiles, we simply summed the population distributions to obtain xi(t), but we can also do the reverse: reconstruct the approximate pi(E,t) from xi(t) and ymA(t) ymB using eqn (24). Thus, we can also examine how the approximate population distributions evolve in time to see in more detail each method’s strengths and weaknesses. Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 1145 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online Fig. 5 Comparison of rate coefficients versus (a) pressure and (b) temperature for CH3CO + O2 - products estimated using the modified strong collision (MSC), reservoir state (RS), and chemically-significant eigenvalue (CSE) methods. In the plots, A = acetylperoxy, B = hydroperoxylvinoxy, C = ketene + HO2, and D = lactone + OH. The error between the three methods is generally within an order of magnitude. However, over the range where the CSE method can be applied, the RS method gives rate constants that better match the CSE values than the MSC method. The CSE method runs into numerical difficulties at high temperature and low pressure, and at low temperature; see text for details. Fig. 7 shows the approximate population distributions pi(E,t) versus energy for the acetylperoxy isomer from each method and the full solution at several times. Under both sets of conditions, the high-energy grains of acetylperoxy are populated before collisional stabilization begins to be important. The collisions tend to cause only small changes in energy, so even after ten collision times there have not been enough collisions to significantly populate the lowenergy grains. The result is that the low-energy grains are not yet Boltzmann-distributed, and so none of the approximate methods match the true populations at low energies. As time progresses to 100 and 1000 collsion times, the low-energy grains become increasingly thermalized, and so the approximate methods much more accurately reflect the full solution. 1146 Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 Of the three methods, the chemically-significant eigenvalue method is the first to accurately match the entirety of the ‘‘true’’ population distribution. As the eigenvalue method makes no assumption about the nature of the energies, it is able to capture the non-Boltzmann, non-steady state behavior of the high reservoir and near threshold grains (about 15 to 0 kcal mol1). (A side effect of this is that the eigenvalue method predicts that the lowest-energy grains have slightly negative populations at short times; this is why the distribution drops sharply off at certain times.) The reservoir state method distributions have (by assumption) a Boltzmann distribution at low energies, which in reality is only established after many collision times. However, only the most populated low-energy This journal is c the Owner Societies 2012 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online Fig. 6 Concentration profiles at low-temperature and high-temperature conditions, comparing the modified strong collision (MSC), reservoir state (RS), and chemically-significant eigenvalue (CSE) methods with the full master equation solution, for the acetyl + oxygen system as shown in Fig. 3. In the plots, A = acetylperoxy, B = hydroperoxylvinoxy, C = ketene + HO2, D = lactone + OH, and E = acetyl + O2. The reservoir state and chemically-significant eigenvalue methods both perform very well at low temperatures. All three methods are comparably accurate at high temperatures even though they disagree with the full solution at short times. The vertical line indicates the average time between collisions with the bath gas. grains need to be Boltzmann-distributed before the total population is predicted accurately. The modified strong collision method population distributions are quite similar to the reservoir state method distributions despite the significantly different collision models used. The main discrepancy is in the low reactive energies (around 6 to 0 kcal mol1 below the reactant energy, but above the lowest isomerization barrier), which the reservoir state method predicts more accurately since it allows for collisional energy transfer between grains. From the plots in Fig. 6 and 7, we see that the reservoir state method performs most consistently at low temperatures, and is still competitive at moderately high temperatures. The kinetics of the major pathways are consistently predicted across all methods; the disagreement appears only in the minor pathways. The chemically-significant eigenvalue method also performs very well under the conditions shown—often equal to or better than the reservoir state method—but would not be successful at lower pressures due to the blending of the chemical eigenmodes into the internal energy eigenmodes, as shown in Fig. 8. 6.2.4 Eigenvalue analysis. In Fig. 5b the k(T,P) curves given by the chemically-significant eigenvalues method appear to stop as temperature decreases at all pressures. To examine why, let us look in more detail at the eigenvalues themselves. Treating the two bimolecular product channels as irreversible sinks, there are three chemically-significant eigenvalues, and no zero eigenvalue. The computed eigenvalues are plotted in Fig. 8a and b as a function of temperature and pressure, respectively. The figures show the four lowest-magnitude eigenvalues, with l4 representing the lowest internal-energy eigenmode. In Fig. 8b the magnitudes of the chemical and internalenergy eigenvalues become increasingly disparate as temperature This journal is c the Owner Societies 2012 decreases, indicating that the master equation matrix is becoming increasingly stiff. Once the temperature becomes too low, the large separation in time scales makes the lowmagnitude eigenvalues numerically difficult to calculate when using double-precision arithmetic—as was used throughout this work—resulting in one or more unphysical positive eigenvalues.75,76 A quadruple-precision eigenvalue algorithm would allow for successful determination of the eigenvalues at low temperature, but we believe this to be unnecessary due to the success of the reservoir state method at low temperatures. The reservoir state method is not immune to numerical issues due to stiffness, but this is generally only problematic for much larger reaction networks than acetyl + oxygen. 6.3 Perturbations to the acetyl + oxygen potential energy surface Although there are a few regions in Fig. 5 where the three approximate methods for estimating k(T,P) values disagree, the magnitude of the disagreement is typically within an order of magnitude. In order to exaggerate the circumstances under which one or more methods differ significantly, we introduce one or more artificial perturbations to the potential energy surface. These perturbations involve raising the ground-state energies or lowering the transition-state barrier heights for one or more configurations on the surface, in order to vary the degree of chemical activation and the relative rates of isomerization. Several such perturbations will be discussed in what follows. For all perturbations we have provided plots of k(T,P) versus pressure for rate coefficients corresponding to acetyl + oxygen as the reactants. Many other plots of k(T,P) values and some concentration profiles are available in the Supporting Materials.w 6.3.1 Acetyl + 20 kcal mol1. We begin with the case where the acetyl + oxygen ground-state energy is artifically Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 1147 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online Fig. 7 Population distributions for acetylperoxy versus energy after 10, 100, and 1000 collision times, comparing the modified strong collision (MSC), reservoir state (RS), and chemically-significant eigenvalue (CSE) methods with the full master equation solution. At the conditions used above, o1 is equal to 3.4  109 and 6.9  1011 s, respectively. The short-lived, high-energy reactive states immediately reach a steady state which is accurately replicated by all three methods. After 100–1000 collision times, each method provides a reasonable approximation of the population distribution of the full solution. The CSE method is the most accurate, while the RS and particularly MSC methods underpredict the population of the adduct states with energies just below the reactant energy (6 to 0 kcal mol1) at lower temperatures. The zero of energy is acetyl + O2. The arrows indicate the energy of the lowest transition state, which is the top of the ‘‘reservoir’’ in the RS method. increased by 20 kcal mol1. This perturbation results in a system that is much more chemically activated than before by creating a significant range of energies below the ground state of acetyl + oxygen but well above all of the transition barriers. The potential energy surface now reflects a very chemically-activated system, as the entire acetyl + oxygen configuration is many hDEdowni above all other reaction barriers in the system. This means that an excited acetylperoxy adduct can survive several collisional stabilization events and remain reactive. Thus we are not surprised to see that very high pressure is required to trap the adducts, as seen in Fig. 9a. However, only the reservoir state and chemically-significant eigenvalue methods capture this effect, while the modified strong collision method does not. The modified strong collision method replaces the collision model from the full master equation with a greatly simplified single-step activation/ deactivation process. This is not a great model of this perturbed 1148 Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 system, as the excited adduct can be stabilized by many kcal mol1 without losing reactivity. The collision efficiency b(T) is intended as a correction to improve the single-step approximation; this perturbation demonstrates a situation where b(T) does not provide an adequate correction. Concentration profile comparisons at both low and high temperature conditions clearly show the modified strong collision method provides a significantly less accurate prediction for the isomer profiles, while the other methods perform noticeably better, especially at low temperatures. 6.3.2 Acetylperoxy + 20 kcal mol1. Next we perturb the system by artificially raising the ground-state energy of the initial adduct acetylperoxy by 20 kcal mol1, creating a much shallower well on the potential energy surface with a much smaller density of states, and so faster k(E). This has the effect of pushing more of the Boltzmann distribution of acetylperoxy into the higher-energy reactive grains, and also encourages This journal is c the Owner Societies 2012 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online Fig. 8 Eigenvalues of the master equation matrix at (a) several pressures at 500 and 1000 K and (b) several temperatures at 0.01 and 1 bar. As the pressure is decreased or the temperature is increased, one or more of the chemical eigenvalues (l1 to l3) become indistinguishable from the internal energy eigenvalues (l4). At low temperatures the eigensolver gives unphysical positive eigenvalues; the last correctly computed low-magnitude eigenvalues are circled. chemically-activated well skipping from acetyl + oxygen to hydroperoxylvinoxy. In Fig. 9b we see that the k(T,P) values from acetyl + oxygen to hydroperoxylvinoxy have increased as a result of this perturbation. At 500 K, this increase is so significant that hydroperoxylvinoxy is predicted to be the dominant product over the range of 1 to 100 bar. All three methods predict this behavior, although a chemical eigenvalue (equilibration of the isomers acetylperoxy and hydroperoxylvinoxy) is often blended into the collisional energy relaxation eigenspectrum, so the chemically-significant eigenvalue method does not succeed with three chemical eigenvalues. For situations like this, Miller and Klippenstein recommend treating the two equilibrated isomers as a single pseudo-species, which would correspond to two chemical eigenvalues.22 When concentration profiles are compared under low pressure conditions—where the k(T,P) value for acetyl + oxygen - acetylperoxy has fallen off—there is some disagreement between the methods on the acetylperoxy profile. This is shown at both low-temperature and high-temperature conditions in Fig. 10. At low temperature, the modified strong collision method consistently overpredicts the acetylperoxy concentration, while This journal is c the Owner Societies 2012 the reservoir state method consistently underpredicts it. The chemically-significant eigenvalue method is significantly more accurate than both of the other two approximate methods in this instance. However, the modified strong collision and reservoir state methods provide very good profiles for the other species as well, with the reservoir state profiles slightly more accurate. Furthermore, under these conditions acetylperoxy is a minor species, so more error in the observed concentration may be acceptable depending on the application. At high temperature the primary disagreement between the modified strong collision and reservoir state profiles is again for acetylperoxy. It is unsurprising that the reservoir state differs from the full master equation solution, especially at short times, as a significant fraction of the Boltzmann population of acetylperoxy now resides in the reactive grains. What is more surprising is how accurate the modified strong collision method is for the same conditions. The single-step collision model used in the modified strong collision method is perhaps more accurate when there are fewer nonreactive energies, especially when coupled with an accurate collision efficiency b(T). For this perturbed potential energy surface, the highpressure-limit rate coefficients kN(T) for reactions involving Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 1149 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online Fig. 9 Comparison of rate coefficients versus pressure for CH3CO + O2 - products estimated using the modified strong collision (MSC), reservoir state (RS), and chemically-significant eigenvalue (CSE) methods. In the plots, A = acetylperoxy, B = hydroperoxylvinoxy, C = ketene + HO2, and D = lactone + OH. (a) Increasing the ground-state energy of acetyl + oxygen by 20 kcal mol1 causes the MSC method to disagree more with the RS and CSE methods. (b) Increasing the ground-state energy of acetylperoxy by 20 kcal mol1 leads to significant discrepancies between the methods, particularly at high temperature and high pressure. The CSE method fails to distinguish three chemical eigenvalues at many of these conditions. acetylperoxy as the reactant are the largest. At 500 K the kN(T) values for acetylperoxy - hydroperoxylvinoxy, acetylperoxy - acetyl + O2, and acetylperoxy - ketene + HO2 are 9.7  109 s1, 3.4  108 s1 and 2.0  106 s1, respectively. At the pressure of 10 bar used for the profiles in Fig. 10, these are all smaller than the average collision frequency of o = 7.1  1010 s1. This implies that the collision and reaction processes do indeed occur at different timescales, and the chemically-significant eigenvalue method is able to resolve three chemical eigenvalues. However, at 1500 K the kN(T) values for the same reactions are 1.9  1011 s1, 1.5  1012 s1 and 2.0  1011 s1, respectively, which are all well above the 1150 Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 average collision frequency of o = 3.2  1010 s1 at 10 bar. This suggests that the reaction time scale blends into the collision time scale, and thus we are not surprised that the chemically-significant eigenvalues method is unable to resolve three chemical eigenvalues at these conditions. 6.3.3 Hydroperoxylvinoxy + 20 kcal mol1. A related perturbation is to artificially raise the ground-state energy of the other isomer, hydroperoxylvinoxy, by 20 kcal mol1, again creating a much shallower well on the potential energy surface. Similar to the last perturbation, this reduces the density of states r(E), increases k(E), and causes more Boltzmann This journal is c the Owner Societies 2012 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online Fig. 10 Concentration profiles at low-temperature and high-temperature conditions, comparing the modified strong collision (MSC), reservoir state (RS), and chemically-significant eigenvalue (CSE) methods with the full master equation solution for the perturbed acetyl + oxygen system where the ground-state energy of acetylperoxy is artificially increased by 20 kcal mol1. In the plots, A = acetylperoxy, B = hydroperoxylvinoxy, C = ketene + HO2, D = lactone + OH, and E = acetyl + O2. At low temperature, all three methods agree with each other and the full master equation solution except for acetylperoxy at short times, for which the CSE method provides the best match of the full solution. At high temperature, the MSC profiles are noticeably more consistent with the full master equation solution than the RS profiles, particularly for acetylperoxy. The CSE method is unable to find three chemical eigenvalues at the high temperature conditions, so its profiles are not shown). population of hydroperoxylvinoxy to exist in the reactive grains. Overall this encourages well skipping from acetylperoxy to the product channels, and reduces the lifetime of the hydroperoxylvinoxy. The result is the k(T,P) values for reactions that produce hydroperoxylvinoxy are depressed, as seen in Fig. 11a. The only major disagreement in k(T,P) values between the modified strong collision and reservoir state methods are for those involving hydroperoxylvinoxy as a reactant or product. (Again, the chemically-significant eigenvalue method for Nchem = 3 has no results because a chemical eigenvalue has merged into the internal energy relaxation eigenspectrum). Concentration profile comparisons to full solutions of the master equation again show that the reservoir state method performs slightly better under most conditions for all configurations except the very unstable species hydroperoxylvinoxy, for which the modified strong collision method is slightly more accurate. 6.3.4 Isomerization  20 kcal mol1. As a final perturbation, we lower the reaction barrier for the isomerization reaction by 20 kcal mol1. The results of this perturbation are shown in Fig. 11b. The low barrier causes the isomerization process to be very fast, to the point where a chemical eigenvalue blends into the internal energy relaxation eigenvalues at nearly all conditions shown, causing the chemicallysignificant eigenvalues method to have no result with Nchem = 3. The modified strong collision and reservoir state methods are competitive, although at high temperatures the former performs noticeably better than the latter in predicting both isomer profiles (figure available in the Supporting Materialsw). 6.3.5 Discussion. A few common themes emerge from the above analysis. First, all three methods are fairly accurate for This journal is c the Owner Societies 2012 the major channels across a wide range of conditions. Second, the modified strong collision method performs consistently in all systems, nearly always giving reasonable results, although sometimes noticeably less accurate than the other methods. The modified strong collision method is also numerically inexpensive and robust. For initial explorations of a system, the modified strong collision method is a reasonable choice. Third, there is a wide range of conditions, especially at low and moderate temperatures, where the reservoir state method provides better estimates for the phenomenological rate coefficients than the modified strong collision method. Under these conditions, the equilibrium population distributions for each isomer are narrow enough that only a small fraction of the Boltzmann population exists in reactive grains. When the ground-state energies of either isomer were raised, the reservoir state method provided less accurate estimates for the k(T,P) values involving the unstable isomer, under conditions where the unstable isomer is a minor channel. When the isomerization barrier was lowered, reservoir state method estimates involving both isomers became less accurate. Thus, we conclude that the reservoir state method is most accurate when the isomer wells are ‘‘deep’’, i.e. the fraction of the Boltzmann distribution that exists in reactive grains is insignificant. Fourth, at conditions where the chemically-significant eigenvalues method was able to resolve three chemical eigenvalues, it was consistently the most accurate in its k(T,P) estimates. This is the best method to use for high-accuracy k(T,P) calculations. However, there are a wide range of conditions where this method was unsuccessful due to one or more chemical eigenvalues being merged with the internalenergy relaxation eigenvalues, or because the numerical eigensolver returned an unphysical positive eigenvalue. The eigenvalue-merging issue, i.e. the breakdown of the assumed Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 1151 Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C View Article Online Fig. 11 Comparison of rate coefficients versus pressure for CH3CO + O2 - products estimated using the modified strong collision (MSC), reservoir state (RS), and chemically-significant eigenvalue (CSE) methods. In the plots, A = acetylperoxy, B = hydroperoxylvinoxy, C = ketene + HO2, and D = lactone + OH. (a) Increasing the ground-state energy of hydroperoxylvinoxy by 20 kcal mol1 encourages well skipping between acetylperoxy and the two product channels. The CSE method is unable to determine three chemical eigenvalues at all conditions due to the very fast equilibrium between the isomers. (b) Decreasing the isomerization barrier height by 20 kcal mol1 creates a situation of a fast equilibrium between the two isomers except at extreme pressure and low temperature, causing the CSE method to be unable to resolve three chemical eigenvalues. The discrepancies in predicted isomer yields are usually not important, since the thermal equilibrium of the isomers is fast (B1010 s1). The discrepancies in predicted yields of the product channels (C and D) are more important; the RS method is more accurate at low temperature, while the MSC method is more accurate at high temperature. separation of collision and reaction timescales, has been discussed by Miller and Klippenstein; they suggest that under these conditions, the concept of a phenomenological rate coefficient breaks down.22 Unfortunately, both issues arise frequently even for networks of modest size. 6.4 Computational cost considerations Our results show that the reservoir state and chemicallysignificant eigenvalues methods are somewhat more accurate 1152 Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 than the modified strong collision method in predicting rate coefficients for minor channels. The chemically-significant eigenvalues method has a stronger theoretical basis, but the reservoir state method is nearly as accurate. An estimate of the relative execution time of the methods, shown in Table 2, comes from considering the time needed by each method to generate the data in Fig. 5. The reservoir state method has about 2–3x the computational cost of the modified strong collision method, while the cost of chemically-significant This journal is c the Owner Societies 2012 View Article Online Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C Table 2 CPU time to generate Fig. 5 Method Time (s) Modified strong collision Reservoir state Chemically-significant eigenvalues 4.8 9.8 310 eigenvalues method is nearly 100x the modified strong collision method. These factors will only get larger as more isomer wells are considered, or if finer energy discretization is employed. The high computational expense and difficulties inherent in the chemically-significant eigenvalues method at low and high temperatures make it less useful for automatic estimation of large numbers of k(T,P) parameters. Given the approximate nature of several of the inputs to the master equation—especially the collisional transfer probabilities model and parameters—there is at this point some question as to whether the reservoir state method provides enough of an improvement in accuracy compared with the modified strong collision method to be worth the extra computational effort. This is especially true when we apply these methods using automatically-generated estimated parameters, which themselves can have a significant error. Furthermore, the reservoir state method, while generally faster and more robust than the chemically-significant eigenvalues method, is not immune to numerical difficulties. The reservoir state method involves a linear solve for the pseudo-steady activestate grain populations via Gaussian elimination. For very large networks this elimination can be unsuccessful due to the large matrices involved. We generally observe failures of the LAPACK Gaussian elimination routines when there are more than around thirty unimolecular isomers in the reaction network (typically we allocate about 200 active-state grains per isomer). For these reasons, our automatic mechanism generation code, RMG, provides the user the option of either the modified strong collision or reservoir state methods, and automatically falls back to the modified strong collision method if the reservoir state method is unable to complete the Gaussian elimination step. 7 Conclusion A general framework for estimating the phenomenological rate coefficients for unimolecular reaction networks of arbitrary size and complexity using only macroscopic parameters has been presented. Two important steps in the calculation have been explored in detail: the modeling of the molecular degrees of freedom used to determine the density of states, and the methods for estimating the phenomenological rates from the full master equation model. In the former, the functional group frequency method has been shown to be generally superior to the three frequency method, giving a density of states that is more consistently close to that of the full quantum chemistry-based model. In the latter, the reservoir state method has been shown to be superior to the modified strong collision method except at very high temperatures, giving concentration profiles that more closely follow those of the full master equation solution. These methods are compared with the chemically-significant eigenvalues method This journal is c the Owner Societies 2012 of Klippenstein and Miller and with the full master equation solution in a way that we hope sheds some light on the relationships between these methods. All three approximations to the master equation perform well at estimating the k(T,P) values of the major pathways of a potential energy surface. However, for strongly chemicallyactivated systems, significant disagreement between the methods is observed for minor pathways. The chemicallysignificant eigenvalues method is noticeably more computationally expensive than the other methods, but also gives the most accurate estimates of the k(T,P) values. However, this method has numerical difficulty at low temperatures due to large separation of chemical and internal energy relaxation timescales, and currently requires human adjustment to handle fewer chemical eigenvalues for systems with a fast equilibration. This could perhaps be automated, but to our knowledge this automation is not currently available. The reservoir state method is much less expensive than the chemicallysignificant eigenvalues method, and works very well at low and moderate temperatures. However, when the reservoir approximation is not valid—i.e. for high temperatures and shallow wells—the reservoir state method gives worse estimates for the k(T,P) values than even the modified strong collision method. This could perhaps be addressed by manually removing the unimolecular wells which no longer correspond to stable species, but such a procedure has not yet been made automated. The modified strong collision method is not as accurate as the other methods in the ranges where those methods work well—particularly for strongly chemicallyactivated potential energy surfaces with deep unimolecular wells—but is computationally inexpensive and much more robust over the full temperature range. In practical applications of automatic mechanism generation, most pressure-dependent reaction networks considered are very complex, with a high likelihood of shallow wells being present. For this reason, we recommend the modified strong collision method for general use with automatic mechanism generation. There is room for improvement by using the reservoir state or chemically-significant eigenvalues methods in the ranges where they are reliable, e.g. for reactions which the model is sensitive to. In particular, the more rigorous methods could become viable if one could automate removal of the very short-lived wells at high temperature. However, one would need good error control in the software to catch and recover from cases where the more rigorous methods are not reliable. The result of this work is a pair of modules for estimating molecular degrees of freedom and phenomenological rates that are utilized by and distributed with the open-source automatic mechanism generation code RMG. RMG uses the functionalgroup frequency method exclusively and offers users a choice between the modified strong collision and reservoir state methods. Preliminary versions of both of these codes, written in Fortran 90, are currently available in the RMG package.1 Acknowledgements This work was funded in part by the US Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences, under Contract Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 1153 View Article Online DEFG02-98ER14914.; and by Award No. KUS-I1-010-01, made by King Abdullah University of Science and Technology (KAUST). Their financial support is gratefully acknowledged. Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C References 1 J. W. Allen, R. W. Ashcraft, G. J. Beran, C. F. Goldsmith, M. R. Harper, A. Jalan, G. R. Magoon, D. M. Matheu, S. S. Merchant, J. D. Mo, S. Petway, S. Ruman, S. Sharma, K. M. Van Geem, J. Song, J. Wen, R. H. West, A. Wong, H.-W. Wong, P. E. Yelvington, J. Yu and W. H. Green, RMG (Reaction Mechanism Generator) version 3.2, 2010, http://rmg. sourceforge.net/. 2 F. A. Lindemann, Trans. Faraday Soc., 1922, 17, 598–606. 3 C. N. Hinshelwood, Proc. Royal Soc. A, 1926, 17, 230–233. 4 O. K. Rice and H. C. Ramsperger, J. Am. Chem. Soc., 1927, 49, 1617–1629. 5 L. S. Kassel, J. Phys. Chem., 1928, 32, 1065–1079. 6 R. A. Marcus and O. K. Rice, J. Phys. Chem., 1951, 55, 894–908. 7 A. J. F. Siegert, Phys. Rev., 1949, 76, 1708–1714. 8 A. F. Bartholomay, Bull. Math. Biophys., 1958, 20, 175–190. 9 E. W. Montroll and K. E. Shuler, Adv. Chem. Phys., 1958, 1, 361–399. 10 I. M. Krieger and P. J. Gans, J. Chem. Phys., 1960, 32, 247–250. 11 P. J. Gans, J. Chem. Phys., 1960, 33, 691–694. 12 B. Widom, J. Chem. Phys., 1959, 31, 1387–1394. 13 J. Keck and G. Carrier, J. Chem. Phys., 1965, 43, 2284–2298. 14 J. Troe and H. G. Wagner, Ber. Bunsenges. Phys. Chem., 1967, 71, 937. 15 J. Troe, Ber. Bunsenges. Phys. Chem., 1973, 77, 665. 16 D. C. Tardy and B. S. Rabinovitch, Ber. Bunsenges. Phys. Chem., 1973, 45, 3720–3730. 17 C. W. Gear, Commun. ACM, 1973, 14, 176–179. 18 T. Beyer and D. F. Swinehart, Commun. ACM, 1973, 16, 379. 19 S. E. Stein and B. S. Rabinovitch, J. Chem. Phys., 1973, 58, 2438–2444. 20 D. C. Astholz, J. Troe and W. Wieters, J. Chem. Phys., 1979, 70, 5107–5116. 21 J. R. Barker, Int. J. Chem. Kinet., 2001, 33, 232–245. 22 J. A. Miller and S. J. Klippenstein, J. Phys. Chem. A, 2006, 110, 10528–10544. 23 R. G. Gilbert, S. C. Smith and M. J. T. Jordan, UNIMOL program suite (calculation of fall-off curves for unimolecular and recombination reactions), available from: School of Chemistry, Sydney University, New South Wales, Australia, 1993. 24 V. Mokrushin, V. Bedanov, W. Tsang, M. R. Zachariah and V. D. Knyazev, ChemRate Computer Program, version 1.10, 1999. 25 S. J. Klippenstein, A. F. Wagner, R. C. Dunbar, D. M. Wardlaw and S. H. Robertson, Variflex Version 1.00, 1999, http://ftp.tcg. anl.gov/pub/variflex/Summary.vrfx. 26 S. H. Robertson, D. R. Glowacki, C.-H. Liang, C. Morley and M. J. Pilling, MESMER (Master Equation Solver for Multi-Energy Well Reactions), 2008; an object oriented C++ program for carrying out ME calculations and eigenvalue-eigenvector analysis on arbitrary multiple well systems, 2008, http://sourceforge.net/ projects/mesmer. 27 D. M. Matheu, PhD thesis, Massachusetts Institute of Technology, 2003. 28 J. Troe, J. Phys.Chem., 1979, 83, 114–128. 29 J. Troe, Ber. Bunsenges. Phys. Chem., 1983, 87, 161–169. 30 R. G. Gilbert, K. Luther and J. Troe, Ber. Bunsenges. Phys. Chem., 1983, 87, 169–177. 31 A. Y. Chang, J. W. Bozzelli and A. M. Dean, Z. Phys. Chem., 2000, 214, 1533–1568. 32 N. J. B. Green and Z. A. Bhatti, Phys. Chem. Chem. Phys., 2007, 9, 4275–4290. 33 J. W. Bozzelli, A. Y. Chang and A. M. Dean, Int. J. Chem. Kinet., 1999, 29, 161–170. 34 R. G. Gilbert and S. C. Smith, Theory of Unimolecular and Recombination Reactions, Blackwell Scientific, Oxford, England, 1st edn, 1990. 1154 Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 35 T. Baer and W. L. Hase, Unimolecular Reaction Dynamics, Oxford University Press, 1996. 36 K. A. Holbrook, M. J. Pilling and S. H. Robertson, Unimolecular Reactions, John Wiley and Sons, 2nd edn, 1996. 37 W. Forst, Unimolecular Reactions: A Concise Introduction, Cambridge University Press, 2003. 38 M. J. Pilling and S. H. Robertson, Annu. Rev. Phys. Chem., 2003, 54, 245–275. 39 S. H. Robertson, A. I. Shushin and D. M. Wardlaw, J. Chem. Phys., 1993, 98, 8673–8679. 40 S. C. Smith and R. G. Gilbert, Int. J. Chem. Kinet., 1988, 20, 307–329. 41 J. A. Miller and S. J. Klippenstein, J. Phys. Chem. A, 2002, 106, 4904–4913. 42 N. B. Slater, Proc. Leeds Philos. Lit. Soc., 1960, 6, 691–694. 43 A. Fernández-Ramos, J. A. Miller, S. J. Klippenstein and D. G. Truhlar, Chem. Rev., 2006, 106, 4518–4584. 44 L. Zhu, W. Chen, W. L. Hase and E. W. Kaiser, J. Phys. Chem., 1993, 97, 311–322. 45 K. S. Pitzer, J. Chem. Phys., 1946, 14, 239–243. 46 A. L. L. East and L. Radom, J. Chem. Phys., 1997, 106, 6655–6674. 47 V. D. Knyazev, I. A. Dubinsky, I. R. Slagle and D. Gutman, J. Phys. Chem., 1994, 98, 5279–5289. 48 R. C. Wilhoit, Thermodynamics Research Center Current Data News, NIST, 1975, vol. 3, ch. 2. 49 W. C. Gardiner and A. Burcat, in Combustion Chemistry, ed. W. C. Gardiner, Springer-Verlag, 1984, ch. 8. 50 R. Hanson and F. Krogh, DQED, http://www.netlib.org/opt/dqed.f. 51 J. Burkardt, DQED Fortran90, http://people.sc.fsu.edu/burkardt/ f_src/dqed/dqed.html. 52 J. Troe, J. Chem. Phys., 1977, 66, 4745–4757. 53 M. A. Blitz, K. J. Hughes, M. J. Pilling and S. H. Robertson, J. Phys. Chem. A, 2006, 110, 2996–3009. 54 S. H. Robertson, M. J. Pilling, L. C. Jitariu and I. H. Hillier, Phys. Chem. Chem. Phys., 2007, 9, 4085–4097. 55 S. J. Klippenstein and J. A. Miller, J. Phys. Chem. A, 2002, 106, 9267–9277. 56 J. A. Miller and S. J. Klippenstein, J. Phys. Chem. A, 2003, 107, 2680–2692. 57 J. H. Seinfeld and S. N. Pandis, Atmospheric Chemistry and Physics–From Air Pollution to Climate Change, John Wiley and Sons, 2nd edn, 2006. 58 J. Lee, C. J. Chen and J. W. Bozzelli, J. Phys. Chem. A, 2002, 106, 7155–7170. 59 S. M. Villano, N. Eyet, S. W. Wren, G. B. Ellison, V. M. Bierbaum and W. C. Lineberger, J. Phys. Chem. A, 2010, 114, 191–200. 60 J. V. Michael, D. G. Keil and R. B. Klemm, J. Chem. Phys., 1985, 83, 1630–1636. 61 G. S. Tyndall, J. J. Orlando, T. J. Wallington and M. D. Hurley, Int. J. Chem. Kinet., 1997, 29, 655–663. 62 M. A. Blitz, D. E. Heard and M. J. Pilling, Chem. Phys. Lett., 2002, 365, 374–379. 63 H. Hou, A. Li, H. Y. Hu, Y. Li, H. Li and B. S. Wang, J. Chem. Phys., 2005, 122, 224304. 64 R. K. Talukdar, M. E. Davis, L. Zhu and A. R. Ravishankara, 19th International Symposium on Gas Kinetics, 2006. 65 G. Kovács, J. Zádor, E. Farkas, R. Nádasdi, I. Szilágyi, S. Dóbé, T. Bérces, F. Márta and G. Lendvay, Phys. Chem. Chem. Phys., 2007, 9, 4142–4154. 66 S.-Y. Chen and Y.-P. Lee, J. Chem. Phys., 2010, 132, 114303. 67 A. Maranzana, J. R. Barker and G. Tonachini, Phys. Chem. Chem. Phys., 2007, 9, 4129–4141. 68 Y. Georgievskii and S. J. Klippenstein, J. Phys. Chem. A, 2003, 107, 9776–9781. 69 S. J. Klippenstein, J. Chem. Phys., 1992, 96, 367–371. 70 S. J. Klippenstein, Y. Georgievskii and L. B. Harding, Phys. Chem. Chem. Phys., 2006, 8, 1133–1147. 71 C. F. Goldsmith, S. J. Klippenstein and W. H. Green, Proc. Combust. Inst., 2011, 33, 273–282. 72 M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. J. A. Montgomery, T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Scalmani, G. Cossi, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, This journal is c the Owner Societies 2012 View Article Online 73 74 75 76 W. Chen, M. W. Wong, C. Gonzalez and J. A. Pople, Gaussian 03, 2004. H.-J. Werner, P. J. Knowles, R. Lindh, F. R. Manby and M. Schuetz, et al., MOLPRO, version 2006.1, a package of ab initio programs, 2006. P. Zhang, S. J. Klippenstein, H. Sun and C. K. Law, Proc. Combust. Inst., 2011, 33, 425–432. T. J. Frankcombe and S. C. Smith, Comput. Phys. Commun., 2001, 141, 39–54. T. J. Frankcombe and S. C. Smith, Faraday Discuss., 2001, 119, 159–171. Downloaded by Fritz Haber Institut der Max Planck Gesellschaft on 21 November 2012 Published on 07 December 2011 on http://pubs.rsc.org | doi:10.1039/C1CP22765C O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, This journal is c the Owner Societies 2012 Phys. Chem. Chem. Phys., 2012, 14, 1131–1155 1155