Reply to “Comment on 'When Rate Constants Are Not Enough'” - The

Dec 23, 2015 - Department of Climate and Space Science and Engineering, University of Michigan, Ann Arbor, Michigan 48109-2143, United States...
0 downloads 0 Views 223KB Size
Subscriber access provided by ORTA DOGU TEKNIK UNIVERSITESI KUTUPHANESI

Comment

Reply to: Comment on “When Rate Constants Are Not Enough” John R. Barker, Michael Frenklach, and David M Golden J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.5b06652 • Publication Date (Web): 23 Dec 2015 Downloaded from http://pubs.acs.org on December 23, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 11

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Reply to: Comment by Miller et al. on “When Rate Constants Are Not Enough” by John R. Barker, Michael Frenklach, and David M. Golden John R. Barker,a* Michael Frenklach,b* and David M. Goldenc* a Department

of Climate and Space Science and Engineering, University of Michigan, Ann Arbor, MI 48109-2143 b Department of Mechanical Engineering, University of California, Berkeley, CA 94720-1740 c Department of Mechanical Engineering, Stanford University, Stanford, CA 94305-3032

Introduction The main goals of our paper1 (BFG) were to call attention to the ubiquitous presence of non-steady energy distributions (NSEDs) and to identify diagnostics for their presence. We presented several typical examples illustrating the performance of the diagnostics. We were concerned that proper time-independent elementary rate constants could not be obtained from stochastic simulations if the energy distributions were still evolving. This issue is also a concern in experiments and in practical reaction environments, such as combustion of a fuel, the chemical changes in the atmosphere, etc. We hypothesized that the presence of NSEDs could affect the phenomenological rate constants required to describe the reaction system and, in response to a referee, we attempted to provide estimates of the possible effects. The presence of NSEDs might also confound chemical reaction models that do not account for their presence. How accurate are phenomenological models (of any kind) when the system is disturbed (by a shock wave, for example) and NSEDs exist? Thus we called for testing phenomenological chemical reaction models under conditions when NSEDs are present in order to evaluate their performance and inform users. To the best of our knowledge, the phenomenological models derived using rate constants from the Bartis-Widom2 (BW) procedure, although apparently quite accurate, had never been examined in this way. Miller et al. have responded with a Comment3 on BFG, to which we are replying. (Hereinafter we will refer to the Comment and its authors as Miller et al. without a reference number.)

Discussion of Central Issues We begin by pointing out that Miller et al. make the following sensational claim: "BFG argue that rate constants for the macroscopic reactions modeled cannot be time independent if the energy distribution of reactants and intermediates in the wells is not in a steady state, thus rejecting one of the foundational concepts of *Correspondence: [email protected], [email protected], [email protected] 1

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 11

chemical kinetics… The purpose of this Comment is to demonstrate that it is not necessary to abandon the concept of time-independent rate coefficients and that no revision of this foundational concept is required to describe these systems accurately." This strange statement is a false characterization of our work, since BFG did not advocate adopting time-dependent rate constants. The whole thrust of BFG was to identify conditions under which time-independent elementary rate constants can be obtained from stochastic simulations. Before we address the substance of Miller et al., it is convenient to define some terms. We use the word "phenomenological" generically. For example, we define a phenomenological reaction mechanism as one that is postulated to describe a physical-chemical phenomenon or process. We use the term "phenomenological rate constants" to mean the rate constants for reactions in a phenomenological reaction mechanism. In BFG, we used the term "elementary rate constant" to describe a phenomenological rate constant corresponding to when a steady energy distribution (SED) exists. We use the expression "phenomenological model" to describe the mathematical expression of the phenomenological reaction mechanism and associated phenomenological rate constants; typically this is the usual set of ordinary differential equations. Although one can imagine phenomenological rate constants that depend on time, we do not recommend them and we confine our attention to phenomenological rate constants that are independent of time. We think that these definitions are accepted widely in the chemical kinetics community. The master equation, with specific molecular parameters, may be termed a "mechanistic" model, since it is based on physical-chemical mechanisms . In the same sense, the Lindemann mechanism may also be regarded as "mechanistic": A + M ⇄ A* + M A* → Products

(1, -1) (2)

(Here, A is a thermalized chemical species and A* is vibrationally excited; M is a nonreactive collider gas.) Our generic definition of "phenomenological" is similar to that of Miller et al. The rate constants obtained using the Bartis-Widom (BW) technique and the CSE/IERE technology described by Miller et al. are phenomenological rate constants for a special phenomenological reaction mechanism. In its usual implementations, the BW technique replaces the thousands of eigenvalues of the full transition matrix with just a handful of eigenvalues corresponding exclusively to chemical transformations. Thus the BW phenomenological mechanism does not include excited chemical species like A* and reactions like Reaction 1. It cannot account for energy relaxation, incubation times, and other related phenomena. Because it is a special case, we will refer to BW rate constants, the BW reaction mechanism, and the BW model. Our definition of the BW model includes the IERE/CSE implementation described by Miller et al.

2

ACS Paragon Plus Environment

Page 3 of 11

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

We regard the BW model as approximate. Because the BW model results from discarding most of the eigenvalues, it cannot reproduce all of the solutions of the master equation. In particular, it fails to accurately represent the early time period during which most of the energy relaxation takes place in the master equation simulation (and in real reaction systems). In their Comment, Miller et al. identify a mistake that we made in our calculation of how much rate constants might be affected as a result of NSED conditions. We agree that the method we used is wrong, because it is based on just one reaction flux coefficient, while the BW phenomenological rate constants depend on more than one. As a result, readers should disregard our numerical estimates of the effects of NSEDs. Fortunately, this error in calculating rate constants has little effect on our main conclusions. Although it does greatly constrict the time over which master equation / rate constant treatments are not valid. Miller et al. go on to give a lengthy exposition of their methods and argue that NSEDs are never important following the initial energy relaxation period. However, they recognize that NSEDs significantly affect the rates of reaction during the initial energy relaxation, but argue that this is not a failure of the BW model, because the model is not "valid" during that period. They admonish users not to apply the BW model when it is not valid. NSEDs are a focus of our paper, BFG. In every one of the isothermal and shock-initiated simulations presented by Miller et al., the chemical reaction rates are affected by NSEDs during the initial relaxation. Miller et al. do not present any phenomenological models that can account for this period. The BW approach, which produces “approaching-to-equilibrium” rate constants, is able to identify the time period when the BW method is valid by comparing phenomenological model results with full master equation simulations, as illustrated in Fig 1b of Miller et al. This analysis requires the computation of all eigenmodes, which is numerically challenging and often expensive, given the large size of the matrices and their large condition numbers.4 From the discussion of Bartis and Widom,2 one can estimate that the transient time needed for relaxation is approximately equal to the time constant for what Miller et al. call the "slowest IERE". Miller et al. argue that the initial energy relaxation is essentially complete in the 2-butene reaction system within a time period ~3-5 times as long. Stochastic simulations do not rely on the eigenvalue-eigenvector decomposition and hence it would be an extra computational expense (and not a small one) just to get the eigenmodes and perform the analysis. Thus this means of estimating the duration of the initial energy relaxation is not practical. Instead, when employing stochastic simulations (and we comment on advantages of the latter below), we suggested using the flux coefficients to identify when NSEDs are important. The examples we presented clearly demonstrate the usefulness of flux coefficients for flagging the initial energy relaxation period, which Miller et al. refer to as the “induction period”. Often, this time period is brief and can be neglected, but not

3

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

always, as is apparent, for example, in Fig VI-3 of the Supporting Information supplied by Miller et al. Miller et al. show that the BW model is quite accurate for describing the cis-butene2 system as it approaches equilibrium, even when NSEDs are present. This result is promising, but of course it does not prove that the BW model will be sufficiently accurate for all conditions and for all multi-well reaction systems. Despite this reservation, we feel that Miller et al. have made a good case that the NSEDs that persist in multi-well systems at times much longer than the expected duration of the initial energy relaxation are probably due to the chemical isomerization reactions. In single-well systems, we do not observe NSEDs that persist longer than the initial energy relaxation. The C2H5+O2 single-well example that we presented in BFG shows that the branching ratio is affected by NSEDs, but the NSED conditions did not persist beyond the initial energy relaxation. Thus our revised hypothesis is that the NSEDs observed in multi-well systems long after the expected initial energy relaxation are due to chemical isomerization reactions, and the BW model may accurately describe the chemical transformations over that period. Of course, this revised hypothesis needs to be tested. The system and conditions used by Miller et al. for their demonstrations of the BW model are for a Knudsen flow reactor,5 which cannot be considered as representative of most combustion and atmospheric conditions. It would be good practice to test phenomenological models (of any type and however derived) for accuracy under conditions representative of where they will be used and to report the results in a form that can be easily utilized. Testing can be carried by simulating model conditions that have features similar to actual experiments (like the shock experiments in BFG), and comparing time histories from the phenomenological model to those obtained using the original master equation. This notion motivated the four examples we presented in BFG. Testing could take a form similar to the calculations presented by Miller et al. and should cover a very wide range of conditions. Of particular interest are systems with transition matrices that are not well-conditioned for BW methods. Such testing is also needed to inform users, who need to decide if the computed phenomenological model is of sufficient accuracy for an intended purpose. An anonymous reviewer observed that if bimolecular reactions take place involving a species subject to NSED conditions, the reaction rates may be affected. This observation is correct. Several recent studies have addressed how master equations can be used to account for bimolecular reactions that compete with vibrational relaxation during and after chemical activation events (for examples, see 6-12). Not surprisingly, the results depend on the particular reaction systems and physical environmental conditions. Such effects can be assessed quantitatively by modeling the individual systems under realistic conditions, but the effects are usually quite small, since the rate constants for collisional thermalization are much larger than those for typical bimolecular reactions, as discussed in the cited studies. 4

ACS Paragon Plus Environment

Page 4 of 11

Page 5 of 11

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Much of the remainder of Miller et al. consists of demonstrations of the methods they have developed. They also make many comments about the performance of their methods, compared to stochastic methods, and discuss other issues as well. Because their comments range far beyond the scope of BFG, we select only a few for replies. For example, we briefly present our opinions about the pros and cons of the various methods for solving the master equation and for arriving at rate constants, observations that can guide the systematic inter-comparisons that will be required for authoritative evaluation of the various techniques.

Discussion of Peripheral Issues Rate Constants from Global Fitting This subject was not discussed in BFG. The methods described by Miller et al. (in their Supporting Information, Section III) can be added to the list of other quantitative methods for analyzing complicated reaction systems. Some example global fitting methods can be found elsewhere13-17 and references therein. Note that the curve fitting methods described by Miller et al. fail to describe reaction during the initial energy relaxation, when NSEDs exist. Thus the diagnostics discussed by BFG are needed to identify the initial energy relaxation period. Rate Constants from Flux Coefficients The treatment described by Miller et al. (in their Supporting Information, Section IV) is new to us and may be useful. We look forward to seeing it published following the normal peer-review process. However, it appears that this method also fails to describe reaction during the initial energy relaxation, when NSEDs exist. Chemical Activation Miller et al. raise several questions about our treatment of Chemical Activation. Their discussion was motivated by our section entitled "Ethyl + O2: A single-Well System. Chemical Activation". The heading is somewhat misleading, because this section describes how one can obtain thermal rate constants from stochastic simulations with an initial energy distribution that is not thermal. (The italicized title “Chemical Activation” should have appeared on the following line, similar to the next such heading “Shock Initiated”.) It really doesn't matter how the simulations are initiated, because, after the initial relaxation (i.e. after ~2 µs), the system has lost all "memory" of the initial distribution. One can determine the total phenomenological rate constant from the least squares fit of the C2H5O2 decay (BFG Fig. 1) and the ratio of rate constants from the branching ratio, after taking into account the products produced prior to ~2 µs (BFG Fig. 2). The values thus obtained are essentially the same as those in the Supporting Information provided by Miller et al. In addition, this same calculation shows that under atmospheric conditions the reaction of C2H5 with O2 produces largely C2H5O2 and even under combustion conditions the C2H5O2 produced is not negligible, demonstrating the inadequacy of representing the reaction as simply C2H5+O2C2H4+HO2.

5

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Stochastic simulation of Chemical Activation is not relevant to BFG and the descriptions presented elsewhere18-20 will not be repeated here. However, we note that the Discussion provided by Miller et al. (Section VI. in their Supporting Information) is not correct in some details. We briefly address those details here. We assume that the nascent chemically activated species exists at t=0 because we base our calculations on statistical rate theory, which assumes instantaneous intramolecular vibrational redistribution (IVR) and which does not consider the dynamics of individual collisions. A dynamical theory might be more accurate, but this model produces results in extremely good agreement with experimental data21 and thus seems to be sufficient for most purposes. Vibrationally excited species are produced by chemical activation. Thus NSEDs are intrinsic to the chemical activation process. It is correct that the probability of branching in the C2H5O2 example is not time dependent, but if the simulations are terminated too early, the computed branching ratio will be incorrect. (This is in contrast to our mistaken statement in BFG.) The stochastic simulations of chemical activation are special, since they are used to acquire the average time history of an ensemble of bimolecular reaction events that are all initiated at precisely the same instant. This differs from experiments (and from the usual simulations), in which bimolecular reactions are initiated at random times. In order to determine accurate probabilities for calculating the rate constants, the time history of the chemical activation special ensemble must be followed until thermalization has been completed; if terminated too early, the probabilities will be incorrect. Thus diagnostics for NSEDs are useful to determining when to terminate the simulations. In the stochastic simulations, the time scales are determined entirely by the energy distribution (which depends on time) and the k(E)s (microcanonical RRKM rate constants) utilized in the master equation for unimolecular reaction of the excited molecules. The strict upper limit to the rate of production of product from the chemical activation special ensemble mentioned above is equal to the largest k(E) in the numerical calculation. In contrast, the bimolecular capture rate constant is an average quantity connected via the equilibrium constant to the high pressure limit of the corresponding unimolecular dissociation reaction. The latter is the thermal average of the same k(E)s used in the stochastic simulations. Some k(E)s are higher than the average, others are lower. The chemical activation rate constants (averaged quantities) are equal to the capture rate constant (averaged quantity) multiplied by the probabilities extracted from the stochastic simulations. Since each probability is ≤ unity, the capture rate constant is a strict upper limit to each chemical activation rate constant (averaged quantities), as it should be. Eigenvalue vs. Stochastic Methods Inevitably, some methods perform better under some conditions than others. Knowledge of these characteristics is useful to both software developers and users. The best way to compare different methods is to carry out systematic inter6

ACS Paragon Plus Environment

Page 6 of 11

Page 7 of 11

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

comparisons using test cases designed to reveal the strengths and weaknesses. This is not easy, because performance measures are multi-dimensional and often strongly coupled. Note that the examples in BFG were not selected as test cases. Appropriate test cases would probe the limits of validity of the methods. As far as we know, no comprehensive inter-comparisons have been carried out to compare the relative merits of the stochastic and eigenvalue approaches (the latter include the BW method as implemented by the CSE/IERE technology described by Miller et al.). In the absence of inter-comparisons and because Miller et al. have been so free in expressing their opinions, we decided to do the same, but with the following cautionary note: All three of us have had a good deal of experience using stochastic approaches, but little or no experience using eigenvalue approaches. Two of us (JRB and DMG) have been involved in some informal inter-comparisons with the eigenvalue approach, but the inter-comparisons were not at all comprehensive. Thus our experience is fragmentary and somewhat anecdotal (hence our call for systematic inter-comparisons). It is our opinion that both approaches to solving the master equation can produce accurate time histories of concentrations (or probabilities), internal energy, and other properties. The time histories are valid starting from t=0. Rate constants can be extracted from the time histories. To extract rate constants, the time histories must be analyzed in ways analogous to analyzing experimental data, including the use of curve fitting and other methods. Typically, total rate constants and branching ratios can be determined relatively easily from the time histories, but determining the rate constants of linked isomerization reactions is more challenging. We think that any rate constant that can be measured in an experiment can also be extracted from an appropriate simulation. The BW method (and associated CSE/IERE technology), which does not rely on simulating experiments and generating time histories, can generate all of the rate constants for large complex systems (including linked isomerizations). The BW method is particularly suitable for mass production of rate constants intended for populating huge reaction mechanisms. This is an excellent method when highly detailed chemical models are required and the initial energy relaxation is not important. We do not address questions regarding the utility and accuracy of such large mechanisms, when few of the rate constants have been (or ever will be) measured. The stochastic approach does not have this direct capability, since it was developed to simulate experiments and is tied to time histories. With regard to energy resolution, the stochastic approach scales much better than the eigenvalue approach. This is of particular interest for problems requiring high energy resolution (small energy grain sizes). Such applications include almost all simulations at low temperatures, and chemical activation and photoexcitation at low to moderate temperatures. At high temperatures, the energy distributions are so broad that energy resolution is not as critical. The execution time of a stochastic simulation formally scales as ~O(1/G), where G is the energy grain size, while the 7

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

matrix inversions in the eigenvalue approach scale as ~O(1/G3).22 Thus, replacing G with G/2 formally increases execution time of the stochastic method by a factor of 2, while the execution time of the matrix inversions increases by a factor of 8. The codes implementing the two classes of methods may perform differently from this. For example, because MultiWell19,23,24 employs a hybrid method for solving the ME, its execution time scales by ~O(1/G0.6) for the cis-butene-2 case. We do not know the corresponding scaling factor for the eigenvalue implementation. On the basis of a limited number of intercomparisons, we are under the impression that, for common reaction conditions, codes based on the two approaches perform roughly about the same, on the basis of CPU time. The ratio of CPU times is not constant, but varies with conditions and with the systems being modeled. For example, MultiWell is not well suited to simulating slow reactions at high pressures, and eigenvalue methods encounter numerical challenges at low temperatures and where high energy resolution is required.22 MESMER has been parallelized (as has VariFlex, we think), but MultiWell has not, although it is expected to be "embarrassingly parallel" and at least two groups are working on it. Miller et al. report an impressive "clock time" performance, but it is not clear whether this refers to total CPU time, or whether it should be multiplied by the number of CPUs utilized in a large parallel execution job. In either case, it is impressive. In our opinion, both approaches to solving the master equation, the methods for extracting rate constants, and all of the various implementations have both strengths and weaknesses. We feel that the two general approaches are somewhat complementary and there is value in having more than one method for solving the same problem. If two very different methods obtain the same solution, confidence in the solution is increased accordingly. For certain applications, there are clear distinctions between approaches and codes, and users should choose accordingly. For other applications, any of the approaches will deliver satisfactory results, and other factors will influence a selection. These factors may include unique features offered by a particular code, the match to available computer resources, ease of use (which depends on user experience), personal preferences, etc.

Summary In this Response, we reply to the Comment of Miller et al., which ranges far beyond the scope of BFG. We acknowledge that we made a mistake in calculating the possible effects of NSEDs on rate constants, but this error did not affect our main conclusions. In light of the new results presented by Miller et al., we agree that the persistent NSEDs observed in multi-well systems accompany the chemical isomerization reactions and that the BW method (and associated CSE/IERE methodology) is able to provide a good description of the chemical system. However, the BW model fails to describe the period when the initial energy relaxation is taking place. Accurate treatment of this initial period will require alternative phenomenological models or direct master equation simulations.

8

ACS Paragon Plus Environment

Page 8 of 11

Page 9 of 11

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

In this Reply, we also address several other issues, including the need to further test the hypothesis (just mentioned) by conducting evaluations based on realistic simulations for phenomenological models. This is just good practice. Other issues we have addressed briefly are our poor labeling and discussion in the section on chemical activation, some technical details regarding chemical activation simulations, and the merits of different approaches to solving the master equation and arriving at rate constants. In the end, after the extensive discussion that appears in Miller et al. and in this Reply, the question still remains: How accurate are phenomenological models (of any kind) when the system is disturbed (by a shock wave, for example) and NSEDs exist? This question is far too broad to be answered by a single example. It can be addressed by comprehensive testing of phenomenological models in ways similar to those described by BFG and Miller et al.

Acknowledgements JRB acknowledges financial support from the National Science Foundation (Atmospheric and Geospace Sciences, Contract No. 1231842) and NASA (Upper Atmospheric Research Program, Contract No. NNX12AE03G). MF acknowledges financial support by the US Army Corps of Engineers, Humphreys Engineering Center Support Activity (Contract No. W912HQ-11-C-0035). DMG acknowledges support from the Department of Mechanical Engineering, Stanford University. We also thank the anonymous reviewers for useful comments.

References (1) Barker, J. R.; Frenklach, M.; Golden, D. M. When Rate Constants Are Not Enough. J. Phys. Chem. A 2015, 119, 7451−7461, DOI: 10.1021/acs.jpca.5b00640. (2) Bartis, J. T.; Widom, B. Stochastic models of the interconversion of three or more chemical species. J. Chem. Phys. 1974, 60, 3474-3482, doi: 10.1063/1.1681562. (3) Miller, J. A.; Klippenstein, S. J.; Robertson, S. H.; Pilling, M. J.; Shannon, R.; Zádor, J.; Jasper, A. W.; Goldsmith, C. F.; Burke, M. P. Comment on “When Rate Constants Are Not Enough, by John R. Barker, Michael Frenklach, and David M. Golden". J. Phys. Chem. A, Submitted. (4) Templates for the Solution of Eigenvalue Problems: A Practical Guide; Bai, Z.; Demmel, J.; Dongarra, J.; Ruhe, A.; van der Vorst, H., Eds.; Society for Industrial and Applied Mathematics, 2000. (5) Bedanov, V. M.; Tsang, W.; Zachariah, M. R. Master Equation Analysis of Thermal Activation Reactions: Reversible Isomerization and Decomposition. J. Phys. Chem. 1995, 99, 11452- 11457. (6) Glowacki, D. R.; Lockhart, J.; Blitz, M. A.; Klippenstein, S. J.; Pilling, M. J.; Robertson, S. H.; Seakins, P. W. Interception of Excited Vibrational

9

ACS Paragon Plus Environment

The Journal of Physical Chemistry

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(7)

(8)

(9)

(10) (11)

(12)

(13)

(14) (15)

(16)

(17)

(18)

(19)

Quantum States by O2 in Atmospheric Association Reactions. Science 2012, 337, 1066-1069, DOI:10.1126/science.1224106. Gonzalez-Garcıa, N.; Olzmann, M. Kinetics of the chemically activated HSO5 radical under atmospheric conditions – a master-equation study. Phys. Chem. Chem. Phys. 2010, 12, 12290–12298, DOI: 10.1039/C0CP00284D. Green, N. J. B.; Robertson, S. H. General master equation formulation of a reversible dissociation/association reaction. Chem. Phys. Letters 2014, 605, 44-46, 10.1016/j.cplett.2014.05.012. Maranzana, A.; Barker, J. R.; Tonachini, G. Master Equation Simulations of Competing Unimolecular and Bimolecular Reactions: Application to OH Production in the Reaction of Acetyl Radical with O2. Phys. Chem. Chem. Phys. 2007, 9, 4129 - 4141, DOI:10.1039/b705116f. Olzmann, M. On the role of bimolecular reactions in chemical activation systems. Phys. Chem. Chem. Phys. 2002, 4, 3614–3618. Peeters, J.; Müller, J.-F.; Stavrakou, T.; Nguyen, V. S. Hydroxyl Radical Recycling in Isoprene Oxidation Driven by Hydrogen Bonding and Hydrogen Tunneling: The Upgraded LIM1 Mechanism. J. Phys. Chem. A 2014, 118, 8625–8643, DOI: 10.1021/jp5033146. Pfeifle, M.; Olzmann, M. Consecutive Chemical Activation Steps in the OHInitiated Atmospheric Degradation of Isoprene: An Analysis with Coupled Master Equations. Int. J. Chem. Kinet. 2014, 46, 231-244, DOI: 10.1002/kin.20849. Andraos, J. A Streamlined Approach to Solving Simple and Complex Kinetic Systems Analytically. J. Chem. Educ. 1999, 76, 1578-1583, DOI: 10.1021/ed076p1578. Boyd, S.; Vandenberghe, L. Convex Optimization; Cambridge University Press: Cambridge, UK, 2004. Frenklach, M.; Packard, A.; Feeley, R. Optimization of Reaction Models with Solution Mapping. In Modeling of Chemical Reactions; Carr, R. W., Ed.; Elsevier: Amsterdam, 2007; Vol. 42; pp 243–291. Lewis, E. S.; Johnson, M. D. The Reactions of p-Phenylene-bis-diazonium Ion with Water. J. Am. Chem. Soc. 1960, 82, 5399–5407, DOI: 10.1021/ja01505a027. Mucientes, A. E.; Peña, M. A. d. l. Kinetic Analysis of Parallel-Consecutive First-Order Reactions with a Reversible Step: Concentration–Time Integrals Method. J. Chem. Educ. 2009, 86, 390-392, DOI: 10.1021/ed086p390. Barker, J. R. Monte-Carlo Calculations on Unimolecular Reactions, EnergyTransfer, and IR-Multiphoton Decomposition. Chem. Phys. 1983, 77, 301318. Barker, J. R. Multiple-well, multiple-reaction-path unimolecular reaction systems. I. MultiWell computer program suite. Int. J. Chem. Kinet. 2001, 33, 232-245.

10

ACS Paragon Plus Environment

Page 10 of 11

Page 11 of 11

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

(20) Golden, D. M. What Do We Know About the Iconic System CH3+CH3+M↔C2H6+M? Z. Phys. Chem. 2011, 225, 969–982, DOI 10.1524/zpch.2011.0139. (21) Weston, R. E., Jr.; Nguyen, T. L.; Stanton, J. F.; Barker, J. R. HO + CO Reaction Rates and H/D Kinetic Isotope Effects: Master Equation Models with ab Initio SCTST Rate Constants. J. Phys. Chem. A 2013, 117, 821-835, http://dx.doi.org/10.1021/jp311928w,. (22) Glowacki, D. R.; Liang, C. H.; Morley, C.; Pilling, M. J.; Robertson, S. H. MESMER: An Open-Source Master Equation Solver for Multi-Energy Well Reactions. Journal of Physical Chemistry A 2012, 116, 9545-9560, 10.1021/jp3051033. (23) Barker, J. R. Energy Transfer in Master Equation Simulations: A New Approach. Int. J. Chem. Kinet. 2009, 41, 748-763. (24) MultiWell-2014.1b Software. Barker, J. R.; Ortiz, N. F.; Preses, J. M.; Lohr, L. L.; Maranzana, A.; Stimac, P. J.; Nguyen, T. L.; Kumar, T. J. D., 2014, http://clasp-research.engin.umich.edu/multiwell/.

11

ACS Paragon Plus Environment