Sum over Histories Representation for Chemical Kinetics - The

As with Feynman path integrals, chemical pathway theory is often more .... obtained by integrating the species decay probability computed from the tot...
0 downloads 0 Views 557KB Size
Download PDF
Letter pubs.acs.org/JPCL

Sum over Histories Representation for Chemical Kinetics Shirong Bai,† Dingyu Zhou,† Michael J. Davis,‡ and Rex T. Skodje*,† †

Department of Chemistry and Biochemistry, University of Colorado, Box 215, Boulder, Colorado 80309, United States Chemical Science and Engineering Division, Argonne National Laboratory, 9700 South Cass Avenue, Argonne, Illinois 60439, United States

‡

S Supporting Information *

ABSTRACT: A new representation for chemical kinetics is introduced that is based on a sum over histories formulation that employs chemical pathways defined at a molecular level. The time evolution of a chemically reactive system is described by enumerating the most important pathways followed by a chemical moiety. An explicit formula for the pathway probabilities is derived and takes the form of an integral over a time-ordered product. When evaluating long pathways, the time-ordered product has a simple Monte Carlo representation that is computationally efficient. A small numerical stochastic simulation was used to identify the most important paths to include in the representation. The method was applied to a realistic H2/O2 combustion problem and is shown to yield accurate results.

M

dX i = Fi(X ) dt

ass action chemical kinetics plays an important role in understanding and predicting the behavior of large chemical networks. These networks are built on underlying mechanisms of potentially large numbers of elementary chemical reactions acting in concert to describe an overall chemical change. Mechanisms describing problems such as hydrocarbon combustion, surface catalysis, atmospheric chemistry, and biological chemistry can consist of thousands of elementary steps and hundreds of independent species. Numerical kinetic models based on these mechanisms provide important predictive tools that are in widespread use across the field of chemistry.1,2 Unfortunately, complex kinetics mechanisms can be quite difficult to physically understand and can become computationally prohibitive if they are too large. Here we offer a new approach to the analysis and computation of chemical kinetics based on chemical pathways that yields new insight and computational power. Our approach is based on the familiar notion of a chemical pathway in which a chemical moiety is followed through a sequence of steps where the product of one reaction becomes the reagent for the next. Our method goes well beyond heuristic reaction path diagrams because we develop a quantitative scheme to actually compute kinetic observables such as concentrations using this representation. There have been two standard approaches to the solution of chemical kinetics problems. First, the continuous concentrations (number of molecules per unit volume) of the distinct species Si as functions of time, Xi(t) ≥ 0, can be expressed as solutions of differential rate equations.3 For the case of a spatially homogeneous system of N-species in a closed system, the rate equations are © XXXX American Chemical Society

(1)

The rate expression for species i, Fi, is composed of the source and sink terms from each elementary reaction in the mechanism, that is sources

Fi(X ) =

∑

sinks

νi , jR j(X ) −

j

∑ νi̅ ,jR j(X ) (2)

j

The elementary rate law for reaction j is given by Rj(X). The molecularities, νi,j, enforce elemental mass conservation for a balanced chemical equation for each reaction reactants

∑ i

products

νi̅ , jSi ⇄

∑

νi , jSi

i

reaction j (3)

This rate equation approach should be regarded as a mean-field theory that is valid in the limit of large concentrations. A second approach is based on a discrete atomistic representation in which individual molecules are allowed to react in a stochastic model.4 This viewpoint, pioneered by the early work of McQuarrie5 and later by Gillespie and coworkers,6−8 makes use of Monte Carlo (MC) sampling methods to propagate large ensembles of discrete molecules. The atomistic approach has proven useful in problems, such found in biochemistry, where concentrations are very low and hence number fluctuations can play a role. A convenient framework for this method involves the chemical master Received: October 22, 2014 Accepted: December 18, 2014

183

DOI: 10.1021/jz502239v J. Phys. Chem. Lett. 2015, 6, 183−188

Letter

The Journal of Physical Chemistry Letters equation (CME), where one solves for the probability distribution function P(x,t|x 0 ,t 0 ). The variables x = (x1,x2,...xN) represent the number of species molecules in a given volume and P(x,t|x0,t0) is the probability that the species are populated by x at time t given that they were populated as x0 at time t0. Obviously, the solution to the CME contains vastly more information than the mean field eq 1. It is also much more difficult to solve the CME than the concentrationbased kinetics equations. The most effective numerical techniques for solving the CME involve averaging many stochastic trajectories that simulate the time evolution of a large ensemble of molecules. We propose a new method to simulate the time dependence of a chemical network that draws its motivation from Feynman’s “sum over histories” approach to quantum mechanics.9 The basic idea of the method is to evaluate kinetic observables as a sum over all relevant chemical pathways that connect an initial state to some observable final state. As with Feynman path integrals, chemical pathway theory is often more computationally intensive than traditional approaches but can provide new insight or numerical advantage for select problems.10,11 A chemical pathway is multistep sequence of elementary chemical reactions represented at a molecular level in which the product of one reaction becomes the reactant for the next reaction. For example, in methane oxidation the three-step chemical pathway +OH

+OH

+O

R1

R2

R3

Figure 1. Schematic diagram showing three chemical pathways that deliver a tagged atom from species S3 to S5. The transformations between species occurs at random times based on the survival probability the species.

problem of determining the concentration of some species Si at a time tf, Xi(tf), given the initial conditions X0 at time t0. We follow a specific atom found in Xi(tf) that may potentially originate from one of several initial species, Sk, each of which may contain more than one chemically equivalent followed atom, viz. μk. The pathways are labeled by “mk” where k specifies the species of origin for the path of the followed atom. The probability of the atom following pathway mk from initial species Sk arriving at the target species Si at time tf is given by Pmk(tf). Then, the concentration of Si is given by the sum over paths leading to Si

CH4 ⎯⎯⎯⎯→ CH3 ⎯⎯⎯⎯→ CH 2 ⎯→ ⎯ CH 2OCH 2 O yields CH 2 O.

The pathway is defined through a sequence of species and reactions that connect them. The notion of the chemical pathway has been often employed to provide mechanistic insight into complex networks of reactions.12−17 While there are a variety of schemes that may be used to define a chemical pathway, a method that is simple to apply involves an atom following procedure. Hence, a given atom is “tagged” at time t0 and is followed as it passes from one species to the next until a final time tf, for example, the carbon-atom in the path above. The history of a given atom is thus an ordered sequence of species that contain the atom and reactions that cause the atom to migrate from one species to another. Hence, if an atom hops from species S0 to S1 to S2 (using the variables Si to denote the ith species along the path) due to reactions R1 and R2 (using variables Rj to denote the jth reaction along the path), the pathway is abbreviated as (S0 R1 S1 R2 S2). We can further specify the path by noting the reaction times, t1 and t2, along with the initial and final time, that is, (t0|S0 R1(t1) S1 R2(t2)S2|tf). The chemical histories of an atom that goes from an initial to a final species may be, in principle, of any length “n” passing through n + 1 species via n reaction steps occurring at times tf > tn > tn−1...> t1 > t0. In Figure 1, we show a schematic of several chemical pathways that may deliver a tagged atom from species S3 to S5. It should be noted that an apparent ambiguity can arise in the atom following procedure if a given species contains the followed type of atom at two or more chemically distinct sites. For example, following the bold H atom in the insertion reaction of formaldehyde, we might have either step H+HHCO → HHCOH or H+HHCO → HHCOH or some combination of the two. However, it is important to realize that no experiment will distinguish such identical atoms, and any choice of branching ratio that resolves the ambiguity that is consistent with the kinetics equations defines a valid pathway model. To illustrate the direct relationship between chemical pathways and conventional kinetics, consider the standard

Xi(tf ) ·μi =

∑ X 0,kPm (tf )·μk k

mk

(4)

For simplicity of notation, we define Xi to be concentration of species Si. In this method, each chemically equivalent followed atom in a reagent species is sampled equally, and the sum of all probabilities for a given followed atom is 1, and thus ∑kmkfixedPmk = 1. Hence, the solution to kinetics problems can generally be reduced to enumerating the relevant chemical pathways and determining their probabilities. We note that unlike a conventional kinetic simulation where the full problem must be solved all at once, the contributions from pathways can be added one-at-a-time to obtain convergence of a kinetic attribute as in eq 4. The key to the pathway representation is the development of an efficient scheme to generate the pathway probabilities. We shall compute the probability for any given pathway using arguments from the theory of sequences of random events, that is, Markov chains. We note that some of our methods bear a similarity to stochastic simulation kinetics although the aims and detailed methodologies are different. Furthermore, the present pathway method can be vastly more efficient. We wish to compute the probability for an n-step chemical path (t0| S0R1(t1)S1...Rn(tn)Sn|tf). We can build up this total probability by combining the probabilities for the individual steps along the path. We define the contingent probability for the ith step of 184

DOI: 10.1021/jz502239v J. Phys. Chem. Lett. 2015, 6, 183−188

Letter

The Journal of Physical Chemistry Letters the path wherein species Si−1 reacts via Ri during a narrow interval of time [ti,ti + dti] to produce Si assuming that the Si−1 was created at ti−1 and that no other reaction occurs in the interval [ti−1,ti] to be i th step prob = ρR (Si , ti|Si − 1 , ti − 1) × dti i

Ppath = ( −1)n n i=1

(5)

≡ exp( −

∫t

tb sin ks j

A i (t ) d t )

a

(6)

where we abbreviate the subscript of P as “i” for brevity. The quantity Ai,j(t) is similar to a pseudo-first-order rate coefficient (with dimension 1/time) and is defined by the time-dependent expression A i , j (t ) =

νi̅ , jR j(X (t )) X i (t )

(7)

Ppath =

where i labels the species and j labels the reaction. The concentrations X(t) are evaluated along the solution to eq 1. The input of a reference trajectory, X(t), appears to be an external input to the pathway theory, although an alternative approach is discussed later. The contingent reaction probability density of species Si−1 reacting via Ri to form species Si at time ti is given by total reaction probability density of Si−1, −dPi−1(ti−1,ti)/dti, times the branching ratio for the process Ri

≡−

dPi − 1(ti − 1 , ti) Γi − 1(ti) dti

dPi − 1(ti − 1 , ti) Ai − 1,J(ti) γ dti Ai − 1(ti) i − 1

(8)

For convenience, we use the index J to denote the selected reaction Ri. The number γi, where 0 ≤ γi ≤ 1, is the fractional probability that the followed atom goes to the Si species. Usually, γi = 1, but if there are followed atoms occurring in several product species of Ri, then γi ≠ 1. Any convenient choice for γi may be employed that is consistent with eq 1. Because the occurrence of each reaction constitutes a random event, we must integrate the probability density over all allowed times consistent with the time ordering. We now consider the full pathway consisting of n reactions and n + 1 species t1

t2

tn

R1

R2

Rn

t3 0

dt 2

∫t

t2

dt1

0

⎛ dPi − 1(ti − 1 , ti) ⎞ Γi − 1(ti)⎟Pn(tn , tf ) dti ⎝ ⎠

(9)

1 N

N

n

k−1 q (tk ))) ∑ Pn(tnq , tf ) ∏ (Γk− 1(tkq)(1 − Pmin q=1

k=1

where the index q is used to label each string of n-random numbers. There are a number of tractable schemes that might be used to generate pathways that lead from S0 to Sn. The reaction mechanism itself, eq 3, connects the various species through elementary reaction steps. In terms of graph theory,18,19 a pathway consists of a sequential list of molecules (nodes) and reactions (directed edges) connecting the species in which the followed atom resides. Abstract graph theoretic methods can be used to generate atom-following pathways in a manner consistent with the mechanism, although this approach does not identify the most probable paths. The most important pathways can be inferred by running a small number of test trajectories. A followed atom in the t = t0 reactant species, S0, is allowed to hop freely from species to species using a Markov chain of reactions, and its pathway is recorded. The MC sampling is carried out, as above, by uniformly sampling the survival probability of species Si in the interval (Pi−1 min(ti−1),1) and inverted to obtain ti, the reaction time. However, instead of efficiently averaging the analytical expression for the branching ratio, as in eq 9, a second random number is used to select the branch the atom follows, and thus each trajectory maps out a single path. By using a modest size ensemble of trajectories, a set of the most important pathways can be defined and roughly ordered in relative probability. If an extremely large number of MC trajectories is run, the pathway probabilities may be actually converged in this way. In a previous work on simple surface reactions, we converged the pathway probabilities using this method.20 However, we generally choose to compute the

ti

i

∫t

(10)

Si − 1 → Si , that is

ρR (Si , ti|Si − 1 , ti − 1) = −

0

dtn...

This integral of the time-ordered product is the central formula of the pathway representation. We note the final term, Pn(tn,tf), is included to account for the survival of species Sn from tn to the final time tf. The integral over the time-ordered product in eq 9 may be numerically evaluated using quadrature or, if n is large, by a Monte Carlo integration scheme. If the latter approach is adopted, an importance sampling can be implemented by switching the integration variables to the survival probabilities using (dPk/dtk)dtk = dPk. The overall sign in eq 9 is canceled because the direction of integration in dPk is reversed from the direction of dtk. The sampling would involve the selection of nrandom numbers in the following way. A random number, y1, is chosen in the interval (P0min(t0),1) where P0min(t0) = P0(t0,tf) is the minimum survival probability of species S0 in the full integration interval. Then, y1 = P0(t0,t1) is inverted to obtain t1, which is used to evaluate the branching ratio ((A1,J(t1))/ (A1(t1)))γ1 in eq 9. A second random number, y2, in the interval (P1min(t1),1), with P1min(t1) = P1(t1,tf), is chosen and used to find t2 and the branching ratio ((A2,J(t2))/(A2(t2)))γ2. This process is continued until the last step of the path, where the additional “leftover” survival probability Pn(tn,tf) must also be included in the integrand. We thus have the MC expression

∫t ∑ Ai ,j(t ) dt ) a

tb

tf

∏⎜

The time-resolved reaction probability density ρRi can be constructed from the kinetic rate eqs 1 and 2. It is convenient to define the species survival probability, PSi(ta,tb), which is the probability that species Si undergoes no reaction in the general time interval (ta,tb), where PSi(t,t) = 1. Using the rate equations, the survival probability for any species Si is obtained by integrating the species decay probability computed from the totality of sink reactions, that is Pi(ta , tb) = exp( −

∫t

S0 → S1 → S2...Sn − 1 → Sn where the reactions obey the time

ordering tf > tn−1 >...t1 > t0. We construct the full n-step pathway probability by combining the single-step probabilities 185

DOI: 10.1021/jz502239v J. Phys. Chem. Lett. 2015, 6, 183−188

Letter

The Journal of Physical Chemistry Letters

pathway representation is a convergent method for solution to the chemical kinetics and also provides insight into the mechanism of H2O2 formation. To test the convergence of the method, we compute the concentration of H2O2 versus time using XH2O2(t)·2 = XH2(0)∑Pi·2, where probabilities Pi are for H-following paths starting at H2 and terminating in H2O2 and the factors of two reflect the presence of two identical H atoms in H2 and H2O2. The time is scaled so that the chain branching explosion occurs at t = 1, which in real time occurs at 4.4 × 10−3 s. In Figure 2,

probabilities by the much more efficient scheme of eq 9 (or eq 10). In this way, the branching ratios are exactly determined without the need to be generated by MC. We note that the computation of the pathway probability using eq 9 implicitly requires a solution to the kinetics equations, X(t), which we term a reference trajectory. This X(t) is needed for the time dependence of the pseudo-first-order rate coefficients, Ai,j(t). The reference trajectory need not be computed for a purely linear system or a problem in steady state where the A i,j (t) are time-independent and the probabilities are found analytically, as shown in the Supporting Information. In the general case, pathway representation would be most useful for mechanistic interpretation and analysis as opposed to being a method to directly solve for the concentrations. Indeed, such interpretive analysis is a major motivation behind our work and there are significant advantages over the conventional kinetics. The pathway formalism, for example, is ideal for the calculation of occurrence probability for rare events such as the production of trace species. It is also possible to use the pathway representation as a tool to understand the results of a chemical sensitivity analysis. Thus, a reaction may prove sensitive when it is a rate-limiting step along a key pathway or if a reaction controls the branching between alternative pathways. However, there are strategies for using the pathway representation to solve for the kinetics without recourse to a prior solution to the ODEs. Note that the concentrations are given as functions of pathway probabilities in eq 4 and the probabilities are given as functions of concentrations in eq 9. Thus, these functional equations can be combined to solve for the probabilities. We have verified that iteration is stable and gives the correct probabilities for several very simple test problems. Specifically, eqs 4 and 9 are written M M+1 as XM+1 = F((PM = Gk(XM), where M is the 1 ,...,PK )) and PK iteration number and k = 1,...,K labels the pathways. We have not yet explored the general utility of this approach for realistic mechanisms, which would then replace the need for the explicit solution to the kinetics eqs 1. Our emphasis in this initial work is to demonstrate the logical consistency of the method and its use in providing mechanistic insight. As a numerical example, we consider a model of the hydrogen combustion. Here we choose a stoichiometric mixture composed initially of H2 and O2 at T = 1000 K and P = 5.6 bar. The mechanism (presented in the Supporting Information) consists of 21 reversible reactions involving eight distinct chemical species, H2, O2, O, H, OH, HO2, H2O2, and H2O. The hydrogen peroxide species (H2O2) is unique in that the weak O−O bond ruptures very easily at high temperatures so it exists only in very low concentrations and its formation can be often regarded as a rare event. Nevertheless, H2O2 is known to play an important role in the growth of the radical pool necessary for a chain branching chemical explosion. A simulation was carried out including the top 100 H atom following chemical pathways connecting the H2 reactant species to the H2O2 target species obtained using a reference trajectory. These pathways were identified using a “small sample” (∼106) Markov chain trajectory simulation starting from the initial mixture. The pathway probabilities were then computed >1% accuracy using the MC evaluation of the path integral representation given by eq 10. We note that a calculation of similar accuracy for the pathway probability using the CME would require a factor of 106 to 108 more numerical computations because the pathway probabilities are often in the range 10−8 to 10−6. The calculation demonstrates that the

Figure 2. Concentration versus time of the unstable H2O2 species in a realistic chemical simulation computed using the pathway method. The simulation was carried out for a stoichiometric mixture of H2 and O2 at 1000 K and 5.62 bar. The time is in units of ignition delay time, 4.4 × 10−3 s. The concentrations are rescaled to probabilities so that [H2] = 1/2 at t = 0. In the lower panel, the converged numerical solution to eq 1 using the full mechanism is shown with symbols. The concentrations obtained using 1, 5, 10, and 100 chemical pathways combined cumulatively via eq 4 are shown with the colored lines. In the middle panel, the individual concentrations of the five most important paths, A−E defined in the text, are shown on a logarithmic scale. In the top panel, the same concentrations from the five most important paths are shown on a linear scale.

we show the concentration of H2O2 computed by sum of pathway probabilities along with the exact simulation result obtained by solving the kinetics eqs 1 using an accurate integrator. It is clearly seen that the pathway representation converges to the true result. Furthermore, we note that at t = 0.5 the use of just five paths gives the result within 20% and the use of 10 paths gives the result within 8%. Paths with extremely low probabilities can be computed using eq 9 (or eq 10) with no more difficulty than more probable paths, unlike the stochastic trajectory method. Several explicit pathways obtained with MC are shown in Figure 3. In the left panel, one trajectory for each of the three most important pathways is shown (A, B, and C defined later). In the right panel, trajectories for path C obtained with using random times are shown. These three histories that follow path C are generated using distinct sets of random numbers in eq 10 and convey a sense of when individual species may be created during the ignition process. We can achieve a mechanistic understanding of the system by examining the details of the contributing pathways. The 186

DOI: 10.1021/jz502239v J. Phys. Chem. Lett. 2015, 6, 183−188

Letter

The Journal of Physical Chemistry Letters

Because probabilities are implicit functions of the concentrations through the pseudo-first-order rate coefficients, an iterative scheme was suggested to obtain the concentrations. A numerical example of the use of the representation was provided by the H2/O2 combustion problem where the reference trajectory was explicitly computed. The chemical pathways that create the unstable, but very important, H2O2 species were ordered in importance using the pathway representation. The present pathway representation provides an important tool for the analysis of chemical systems that goes beyond conventional treatments that give only instantaneous “snapshots” of chemical paths.21 In particular, by following individual molecules from reagents to products, much can be learned about the sensitivity of kinetic observables to rate coefficients. A comparison between the conventional reaction path method and the present approach is provided in the Supporting Information. Finally, it is also envisioned that by the enumeration of the most important chemical pathways, large mechanisms may be simplified by eliminating reactions that are not involved in the dominant pathways.22

Figure 3. In the left panel, the three most important H-atom following pathways leading to H2O2 in the H2/O2 combustion system are shown. In the right panel, three random stochastic trajectories are shown for the single pathway H2 → HO2 → H2O2. These three “path C” trajectories show how the species HO2 can be generated at a wide variety of times on route to the H2O2 target species.

■

Details of the H2/O2 combustion mechanism are presented. The exact analytical solution for the pathway probabilities in systems of first-order reactions or for systems in steady state are presented. A comparison of the sum over histories representation and the results of traditional reaction path analysis are provided. This material is available free of charge via the Internet at http://pubs.acs.org.

H2O2 molecule can, in principle, be created by any of six possible reactions, listed in the Supporting Information. It is found, however, that only two of these reactions are part of any pathway with appreciable probability, viz. HO2 + H2→ H2O2 + OH and HO2 + HO2→ H2O2 + OH. All of the 100 top pathways generate H2O2 by one of these steps. The most important paths that deliver an H atom from H2 to H2O2 are +OH

+O2

+HO2

H 2 ⎯⎯⎯⎯→ H( +H 2O) ⎯⎯⎯→ HO2 ⎯⎯⎯⎯⎯→ H 2O2 ( +O2 )

H 2 ⎯⎯⎯⎯⎯→ H 2O2 ( +H) +O2

+HO2

H 2 ⎯→ ⎯ H( +OH) ⎯⎯⎯→ HO2 ⎯⎯⎯⎯⎯→ H 2O2 ( +O2 ) +HO2

+O2

+HO2

H 2 ⎯⎯⎯⎯⎯→ H( +H 2O2 ) ⎯⎯⎯→ HO2 ⎯⎯⎯⎯⎯→ H 2O2 ( +O2 ) +OH

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected].

(A)

Notes

+HO2

+O

ASSOCIATED CONTENT

S Supporting Information *

+O2

+H 2

H 2 ⎯⎯⎯⎯→ H( +H 2O) ⎯⎯⎯→ HO2 ⎯⎯⎯→ H 2O2 ( +H)

(B)

The authors declare no competing financial interest.

(C)

ACKNOWLEDGMENTS The work of D.Y.Z., M.J.D., and S.R.B. was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences, under contract number DE-AC0206CH11357.

■

(D) (E)

■

In the top panel of Figure 2, it is seen that the abstraction pathway, B, is more prominent in the early stages of the reaction, while toward the latter stages the radical recombination pathways, A, C, and D, take over. The growth of the radical pool with time and the large value of the rate coefficient for HO2 + HO2 explain this crossover behavior. In summary, we have introduced the basic formulas underlying a new method to simulate the time evolution of chemically reactive systems. Using the idea of a chemical pathway defined through a followed atom approach, the chemical composition of a reacting mixture is obtained using the sum of probabilities for the pathways. Similar in motivation to Feynmann’s sum over histories approach, the chemistry is described by enumerating the pathways that lead from initial reagents to a target species. An explicit formula for the pathway probabilities is derived and takes the form of the integral of a time-ordered product. For long pathways, this integral has a simple MC representation. The pathways themselves were identified numerically using a small sample of stochastic trajectories that naturally follow the most probable paths.

REFERENCES

(1) Miller, J. A.; Kee, R. J.; Westbrook, C. K. Chemical Kinetics and Combustion Modeling. Annu. Rev. Phys. Chem. 1990, 41, 345−387. (2) Hammer, B.; Norskov, J. K. Theoretical Surface Science and Catalysiscalculations and Concepts. Adv. Catal. 2000, 45, 71−129. (3) Laidler, K. Chemical Kinetics, 3rd. ed.; Harper Collins: New York, 1987. (4) Van Kampen, N. Stochastic Processes in Physics and Chemistry; North-Holland Personal Library: Amsterdam, 2007. (5) McQuarrie, D. A. Stochastic Approach to Chemical Kinetics. J. Appl. Probab. 1967, 4, 413−478. (6) Gillespie, D. T. A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions. J. Comput. Phys. 1976, 22, 403−434. (7) Gillespie, D. T.; Hellander, A.; Petzold, L. R. Perspective: Stochastic Algorithms for Chemical Kinetics. J. Chem. Phys. 2013, 138, 170901. (8) Gibson, M. A.; Bruck, J. Efficient Exact Stochastic Simulation of Chemical Systems with Many Species and Many Channels. J. Phys. Chem. A 2000, 104, 1876−1889.

187

DOI: 10.1021/jz502239v J. Phys. Chem. Lett. 2015, 6, 183−188

Letter

The Journal of Physical Chemistry Letters (9) Feynman, R. P.; Hibbs, A. R.; Styer, D. F. Quantum Mechanics and Path Integrals; Dover Publications: New York, 2010. (10) Makri, N. Time-Dependent Quantum Methods for Large Systems. Annu. Rev. Phys. Chem. 1999, 50, 167−191. (11) Cao, C. S.; Voth, G. A. The Formulation of Quantum-Statistical Mechanics Based on the Feynman Path Centroid Density. 1 Equilibrium Properties. J. Chem. Phys. 1994, 100, 5093−5105. (12) He, K.; Ierapetritou, M. G.; Androulakis, I. P. A Graph-Based Approach to Developing Adaptive Representations of Complex Reaction Mechanisms. Combust. Flame 2008, 155, 585−604. (13) Lehmann, R. An Algorithm for the Determination of All Significant Pathways in Chemical Reaction Systems. J. Atmos. Chem. 2004, 47, 45−78. (14) Feng, H.; Han, B.; Wang, J. Dominant Kinetic Paths of Complex Systems: Gene Networks. J. Phys. Chem. Lett. 2010, 1, 1836−1840. (15) Chern, J. M.; Helfferich, F. G. Effective Kinetic Modeling of Multistep Homogeneous Reactions. AIChE J. 1990, 36, 1201−1208. (16) Bendtsen, A. B.; Glarborg, P.; Dam-Johansen, K. Visualization Methods in Analysis of Detailed Chemical Kinetics Modelling. Comput. Chem. 2001, 25, 161−170. (17) Sun, W. T.; Chen, Z.; Gou, X. G.; Ju, Y. G. A Path Flux Analysis Method for the Reduction of Detailed Chemical Kinetic Mechanisms. Combust. Flame 2010, 157, 1298−1307. (18) Bollobas, B. Modern Graph Theory; Springer Verlag: Berlin, 1998. (19) Temkin, N. O.; Zeigarnik, A. V.; Bonchev, D. G. Chemical Reaction Networks: A Graph-Theoretical Approach; CRC Press: Cleveland, OH, 1996. (20) Kramer, Z. C.; Gu, X. K.; Zhou, D. D. Y.; Li, W. X.; Skodje, R. T. Following Molecules through Reactive Networks: Surface Catalyzed Decomposition of Methanol on Pd(111), Pt(111), and Ni(111). J. Phys. Chem. C 2014, 118, 12364−12383. (21) Kee, R. J.; et al. CHEMKIN-PRO; Reaction Design, Inc.: San Diego, CA, 2008. (22) Lu, T.; Law, C. K. On the Applicability of Directed Relation Graphs to the Reduction of Reaction Mechanisms. Combust. Flame 2006, 146, 472−483.

188

DOI: 10.1021/jz502239v J. Phys. Chem. Lett. 2015, 6, 183−188