Connectivity-Based Parallel Replica Dynamics for Chemically

Oct 23, 2013 - Reactive force field methods such as AIREBO, ReaxFF and COMB, are extremely useful for studying physical and chemical interactions betw...
1 downloads 16 Views 852KB Size
Download PDF
Letter pubs.acs.org/JPCL

Connectivity-Based Parallel Replica Dynamics for Chemically Reactive Systems: From Femtoseconds to Microseconds Kaushik L. Joshi,† Sumathy Raman,*,‡ and Adri C. T. van Duin† †

Department of Mechanical and Nuclear Engineering, The Pennsylvania State University, University Park, Pennsylvania 16802, United States ‡ ExxonMobil Research and Engineering, Annandale, New Jersey, 08801, United States ABSTRACT: Reactive force field methods such as AIREBO, ReaxFF and COMB, are extremely useful for studying physical and chemical interactions between molecules and materials. However, many chemical reactions have relatively high activation energies, putting them beyond the times-scale of conventional molecular dynamics (MD) at modest temperatures. To capture the low-temperature long-lived radical chemistry in atomistic simulations, we have developed a new transition detection scheme for performing Reactive Parallel Replica Dynamics (RPRD) simulation enabling an extended MD time-scales, essentially up to a microsecond using ReaxFF. In the newly implemented event detection scheme, the transition events are identified whenever there is a change in connectivity of any atom. 1-Heptene pyrolysis is chosen as a model system, and RPRD simulations are performed at temperatures as low as 1350K for up to 1 μs for a system consisting of 40 heptene molecules. The chemical mechanism and the product distribution that were obtained from the RPRD simulation are in reasonable agreement with shock tube experiments. SECTION: Molecular Structure, Quantum Chemistry, and General Theory

A

calculations are still mostly restricted for smaller systems containing less than 200 atoms. In addition, if one wants to study the dynamical evolution of a reactive system, then the ab initio MD simulations can simulate only up to picoseconds on smaller systems. Molecular dynamics (MD) method using empirical and analytical force fields, however, can simulate much bigger systems consisting of millions of atoms and processes up to nanosecond time-scale. To capture chemical events during the simulation, one can use reactive force field methods such as AIREBO,7 ReaxFF,8 and COMB9 that can model bond formation and bond breaking during a MD simulation. The ReaxFF reactive force field method has been applied for a wide range of materials such as metallic,10,11 covalent,12,13 and ionic metal oxide/hydride/carbide systems.14,15 The reactive molecular dynamics simulations (RMD) can be used very effectively to study the chemical evolution of the system at the atomistic level. However, to resolve atomic vibrations accurately, one has to take a very small time step (of the order of 0.25 or even 0.10 fs at high temperatures) and as a result, these simulations have a severe time-scale restriction and one can go only up to a couple of nanoseconds. Many chemical processes, especially that

ccurate understanding, both in terms of kinetics and reaction mechanisms, of various catalytic processes such as hydrocarbon cracking inside HZSM-5 zeolite, methanol to olefin conversion, alkylation of aromatics, selective oxidation of alkenes, etc., and the thermal pyrolysis and combustion of linear and branched chain hydrocarbons, is of great importance either to improve or to design new materials and operating conditions for many industrial applications. Different experimental and computational techniques are commonly used synergistically to gain a better insight into such complex chemical systems. Although experimental techniques are very useful for understanding the overall chemistry, it is still challenging and cumbersome to provide detailed chemical information of elementary reaction steps and is largely due to the existence of wide range of chemical events involving a large number of intermediate and metastable species. The majority of the classical kinetic models that solve the system of coupled differential equations of complex reaction networks require an a priori knowledge of detailed reaction mechanism, activation energies, and frequency factors. For such information, one has to refer to either experiments1−3 or to quantum calculations.4−6 Electronic structure calculations can provide useful information at atomic scale by providing the optimized structure and energies of reactants and transition states, which can then be used further with statistical machinery to extract kinetic parameters. However, even on modern computers, these © XXXX American Chemical Society

Received: September 6, 2013 Accepted: October 23, 2013

3792

dx.doi.org/10.1021/jz4019223 | J. Phys. Chem. Lett. 2013, 4, 3792−3797

The Journal of Physical Chemistry Letters

Letter

thermal decomposition of 1-heptene at temperatures as low as 1350 K. The new connectivity-based event detector is implemented in LAMMPS28 and is currently coupled to the ReaxFF reactive force field method.29 It is known that reactive molecular dynamics are roughly 2 orders of magnitude slower compared to a simple Lennard-Jones (LJ) potential (http://lammps. sandia.gov/bench.html) due to the construction of the bond connectivity table after every time step. Consequently, the simple approach of comparing the entire bond connectivity matrix to identify the occurrence of a chemical reaction can be a computationally intensive task especially for bigger systems. In reactive systems, bond breaking or bond formation is accompanied by a change in coordination of the atoms that are involved in chemical reaction. Hence, in this implementation, the change in number of neighbors surrounding any atom is used as a main criterion for driving the system from one energy basin to the next basin. This newly implemented event detector is capable of exploiting different functionalities of LAMMPS and is very flexible. For example, in order to neglect isomerization reactions, the event detector can also check for the change in total number of molecules present in the system. This criterion is very useful to filter out low barrier chemical events such as intramolecular hydrogen transfer in C2H5, C3H5, etc. Such low barrier conformational events, if allowed, occur very frequently and subsequently reduce the overall efficiency of RPRD algorithm significantly. Furthermore, to neglect other low barrier events such as the formation and breaking of hydrogen bonds, one can also create a group of atoms that does not contribute to hydrogen bonds and can restrict an event detector only for that group of atoms. Thus, one can customize this event detector to model chemical systems with heterogeneous energy landscape without compromising the efficiency of PRD algorithm. In the first part of this paper, we evaluate the assumption that molecules wander randomly in a particular basin for a very long time, and hence they forget how they reached a particular state in a basin before making an escape to the next basin. This is followed by the PRD results obtained on heptene pyrolysis, the definition of PRD efficiency, and a discussion around how we choose our dephasing and correlation times. For studying the thermal decomposition of heptene, we used the following two criteria for event detection, namely, a change in the number of surrounding neighbors on any atom in the system and a change in the number of molecules present in the system. Before performing actual PRD simulations, it is essential to ensure that the pyrolysis process obeys first-order kinetics. Hence, we performed modified PRD simulations to calculate probability distribution of escape time for each of the two steps. In the modified PRD simulations, PRD runs were started on all replicas in a normal way. However, whenever any replica found any reactive event, instead of broadcasting that event to all the remaining replicas and searching for any correlated event, it recorded the event time and the event coordinates and then waited until all the remaining replicas found their first reactive event. The events were identified using a newly implemented connectivity-based event detector. The simulation was stopped when all the replicas registered their respective first event. From these recorded times, we calculated the probability distribution function for the first reactive event. The thermal decomposition of hydrocarbons consists of initiation, propagation, and termination steps. For the initiation escape time distribution, the simulations were performed on a system containing 39

require thermal activation, often require significantly longer time scales. To circumvent this time-scale limitation, one can perform RMD simulation at elevated temperatures. At elevated temperatures, RMD simulations can provide some useful information and can capture chemical reactions in a short simulation time-scale by increasing the rate of all reactions. However, increasing the simulation temperature often increases not only the rate constants but also changes the relative ratio of rate constants between the low barrier and high barrier reactions. Consequently, at elevated temperatures, entropically favored reactions dominate the low temperature chemistry, thereby affecting the prediction of the overall product distribution. In order to capture the low-temperature chemistry in RMD simulations on systems containing thermally activated events, it is imperative to advance and use accelerated molecular dynamics methods with reactive force fields. Accelerated molecular dynamics methods are designed to help simulate much longer time scales, up to the order of microseconds, while retaining full atomistic details. Most of these methods are derived from transition state theory and are suitable for system with infrequent events. Except parallel replica dynamics (PRD) method, the other remaining methods involve some form of modified dynamics, e.g., in “on-the-fly KMC” and “Temperature Accelerated Dynamics” method,16,17 one has to calculate the activation energy of the chemical event on the fly to advance the system from one state to another while in Hyperdynamics;18 state-to-state transition is accelerated by biasing the original potential to reduce the height of original barrier between the two states and in Metadynamics;19 and a repulsive Gaussian is added at a specific rate to a selected collective variable to increase the occurrence of a rare event. By contrast, PRD achieves temporal parallelization of state-to-state dynamics by simulating number of independent trajectories simultaneously. If the probability distribution of transition events follows the first-order process, then it can be shown that both the sequence of transition events and the transition times of this parallel simulation are indistinguishable from a simulation on a single processor. In that aspect, PRD is the most powerful and the most accurate accelerated molecular dynamics method provided that system under consideration follows first order escape kinetics. The details of the PRD algorithm can be found in the original article written by Art Voter.20 The method has been applied to different types of applications such as diffusion of Li through polymer matrix,21 fracture process of metals,22 stick−slip friction,23 stretching of carbon nanotubes,24 diffusion of defects in plutonium,25 and graphene formation from small nanotube fragments.26 However, all the previous studies were done using nonreactive force fields. The only effort to perform reactive PRD simulation using reactive potential was made by Kum et al.27 in which authors investigated the pyrolysis of n-hexadecane using AIREBO potential. In this study, state-to-state transition was detected whenever there was a change in bond connectivity matrix. The authors reported that using PRD, they were able to see initiation events at temperatures that are 1000 K lower than the direct molecular dynamics simulations. However, their study was done only on single particle system and as a result, the authors observed only unimolecular initiation reactions of n-hexadecane. In this Letter, we present a connectivity-based RPRD implementation that can model multiparticle systems at temperatures approaching experimental conditions using ReaxFF. Using this new event detector, we simulated the 3793

dx.doi.org/10.1021/jz4019223 | J. Phys. Chem. Lett. 2013, 4, 3792−3797

The Journal of Physical Chemistry Letters

Letter

Figure 1. Probability distribution function for the first transition event for two stages of pyrolysis process.

molecules of heptene at 1600 K on 192 replicas. For the propagation step, we chose the system consisting of 38 heptene molecules, one C4H9 radical, and one C3H5 radical. The probability distribution was evaluated at 1400 K based on 192 replica simulation. Figure 1 shows the obtained probability distribution function for both the steps. It can be seen that for both the steps, the probability distribution function is exponential for the first escape event and the obtained distribution can be described by P(t) = 1 − e−kt indicating a first order escape kinetics. The actual heptene cracking process involves more stages and in principle, we can obtain exponential probability distribution for each stage by repeating such calculations for different molecular compositions. However, we believe that investigating these two stages is enough to justify the use of RPRD for pyrolysis of paraffins and olefins. After verifying that pyrolysis process obeys first order kinetics, we performed PRD simulations to study thermal decomposition of heptene. The PRD system consisted of 40 heptene molecules placed in a periodic box with a density of 0.05gm/cm3. The C−H parameters of ReaxFF were taken from Chenworth et al.12 The equations of motion were integrated using velocity-verlet algorithm. The simulations were done in NVT ensemble with time step of 0.1 fs. Nosé−Hoover30 thermostat with damping constant of 100 fs was employed for controlling the system temperature. For choosing an appropriate measure of the dephasing time, we calculated the velocity-autocorrelation function of the same 40 molecules system at 1300 K, 1500 K, and 1800 K. At all the temperatures, the velocity autocorrelation function was decaying toward zero in less than 1 ps. This indicates that the system loses its memory about the initial momentum state in less than 1 ps, and hence we set the dephasing time to 5 ps for all the PRD simulations. The correlated event search time was set to 10 ps. With these settings, we conducted PRD simulations at different temperatures. Table 1 shows the number of replicas used and the total time simulated at each temperature to achieve 47.5% fuel conversion. The reported lowest MD temperature at which pyrolysis of C7 hydrocarbon was studied is above 2000 K.31 However using PRD, on 180 replicas, we were able to observe considerable pyrolysis reactions at 1350 K in 1 μs simulation time. The direct consequence of simulating such a long time scale at a low T is that we are now able to observe pyrolysis reactions at temperatures that are inaccessible

Table 1. Number of Replicas Used and Total Simulation Time Achieved from PRD Simulations of Heptene Pyrolysis at Different Temperatures temperature (K)

no. of replicas

simulated time (ns)

% efficiency

1800 1500 1450 1400 1350

10 80 80 120 180

7.39 174.72 386.18 650.07 992.3 ≈ 1 μs

59.7 64.44 75.69 89.94 93.55

in conventional MD simulations. Table 1 also shows efficiency of PRD, which is calculated using the following formula: % Efficiency =

Actual simulated time Simulation time on each replica No. of replicas × without correlated event

In a reactive atomic scale simulation, the time between consecutive reactive events depends on the temperature and system size (species concentration). In many reactive systems, this time interval between consecutive events varies inversely with temperature, and the events become more infrequent at lower temperatures. The efficiency values reported here are system dependent and can vary from system to system. However, the overall trend that efficiency increases by reducing the system temperature should be valid for most of the reactive systems. As a result, our calculations indicate that, at lower temperatures, one can use higher number of replicas to exploit the time acceleration of PRD to a maximum extent and can simulate longer time scales. From the chemistry point of view, the simulation temperature has significant influence on the overall product distribution and observed reaction mechanism. Figure 2 shows the % mass fraction of (C1 + C2) products as a function of temperature. It can be seen that by lowering the simulation temperature, the amount of methane and ethylene formed reduces progressively. This trend is qualitatively in good agreement with shock-tube experiments.2 Such an enhancement in the product distribution is due to many reasons. At the higher temperatures that are generally employed in conventional MD, entropically favored reactions (ΔS > 0) dominate the overall chemistry. To investigate this entropy effect, we calculated the change in entropy of each reaction that was 3794

dx.doi.org/10.1021/jz4019223 | J. Phys. Chem. Lett. 2013, 4, 3792−3797

The Journal of Physical Chemistry Letters

Letter

Table 2. RPRD Results on the Number and Percentage of Reactions with ΔN < 0 as a Function of Temperature at a Fixed Conversion (47.5%) of 1-Heptenea temperature (K)

simulation time (ns)

no. of reactions with ΔN < 0

% of reactions with ΔN < 0

1800 1500 1400 1350

7.4 174.8 650.0 992.3

13 20 21 23

24.0 27.4 28.0 30.7

a

ΔN is the change in stoichiometry of the reaction.

Figure 2. Effect of temperature on product distribution and entropy change of reaction events.

observed in RPRD simulations. Entropies of products and reactants at each temperature were obtained from NASA’s thermodynamic data.32 It can be seen that (Figure 2) the contribution of entropically less favored (ΔS < 0) reactions increases progressively by reducing the simulation temperature. This gradual shift in reactive events, from entropically more favored to entropically less favored, is essential in obtaining a better product distribution that agrees well with the experiments aimed in understanding the low T cracking/pyrolysis chemistry of hydrocarbons. Different classes of reactions such as initiation reactions, radical propagation reactions, and termination reactions contribute to this shift in chemical events; e.g., in the presence of ethyl radical, heptene molecule can be initiated by H-transfer and dehydrogenation reactions as shown below.

Figure 3. Effect of temperature on average life of ethyl radical.

Besides the entropy effect, temperature also affects the lifetime of intermediate radicals, and Figure 3 shows the effect of temperature on average lifetime of C2H5 radicals that are formed during the simulation. It can be seen that the ethyl radicals have 3 times longer lifetime at 1350 K than at 1500 K. This trend indicates that at lower temperatures, the intermediate radicals are able to survive longer. If the intermediate radicals survive longer, then they can diffuse through the system and can explore different chemical pathways other than just undergoing fragmentation/β-scission reactions. This increase in average radical lifetime is responsible for lowering the amount of C1 and C2 products during the pyrolysis simulation. Figure 4 shows the reaction network from RPRD simulations for ethyl and propyl radicals at 1350 K and 1800 K. At 1350 K, in addition to heptene pyrolysis, we also observed initiation of butene and propene by ethyl radical. However, such events were not seen at 1800 K. It can be seen that the reaction pathways that lead to formation of heavier fragments like C5H10, C6H14, C7H12, etc. are more likely to occur at lower temperatures than at higher temperatures. This trend further highlights the entropy effect and shows that longer simulation time can assist atomistic simulations in capturing low temperature chemistry. We also tried to evaluate the apparent rate constants and the corresponding Arrhenius parameters for the pyrolysis of heptene. The rate constant at each temperature was determined from the dynamic trajectory of the heptene being fitted to the first order decay expression, ln No − ln Nt = kt. Subsequent plot of ln k versus 1/T enabled the determination of Arrhenius parameters, namely the activation energy (Ea) and preexponential factor A in the expression k(T) = A exp(−Ea/ RT). Table 3 shows the comparison between the calculated parameters and the experimental results. It can be seen that the activation barriers calculated from PRD simulations are in reasonably good agreement with experimental values.33 We have successfully implemented a connectivity-based event detector in LAMMPS for performing reactive PRD

C7H14 + C2H5 → C7H15 + C2H4

C7H14 + C2H5 → C7H13 + C2H4 + H 2

Our reaction analysis shows that contribution of dehydrogenation pathway (ΔS > 0) to thermal decomposition decreases progressively with simulation temperature; e.g., at 1800 K, 26% of pyrolyzed fuel was initiated by dehydrogenation. However, that contribution went down to 21% at 1450 K and to 10% at 1350 K. Efforts were also made to examine the results from the reactive MD simulations for the increased probability for reactions with ΔN, the change in the number of molecules between products and reactants, being less than zero with decreasing reaction T. The trend in ΔN is taken as a rough measure of entropy effects in the reactive system. We first identified the time “tT” needed for a fixed extent of conversion of reactants as a function of reaction temperature, T. We then counted the total number of reactive events (Ntot) at each temperature and further categorize the reactive events into reactions with positive (or) negative change in number of molecules (ΔN as +ve (or) −ve). The results tabulated in Table 2 reveal the increased percentage of reactions with negative ΔN with decreasing T, suggesting that an atomistic simulation performed at reasonably low enough T has a potential to predict a more accurate product distribution. 3795

dx.doi.org/10.1021/jz4019223 | J. Phys. Chem. Lett. 2013, 4, 3792−3797

The Journal of Physical Chemistry Letters

Letter

Figure 4. Reaction pathways of ethyl and propyl radicals obtained from RPRD simulations. The thickness of the arrow is proportional to the frequency of that pathway observed in RPRD.

Table 3. Comparison of Arrhenius Parameters between PRD Calculations and Experiments Arrhenius parameters

PRD

experiment

Ea (kcal/mol) pre-exponential factor, A

50.01 7.17 × 1013

54 6.3 × 1011

ACKNOWLEDGMENTS

■

REFERENCES

A.C.T.v.D. acknowledges funding from NSF/SI2 grant#1047857, Air Force Office of Scientific Research, under AFOSR-MURI Grants FA9550-10-1-0563 and AFOSR Grant FA9550-11-1-0158.

(1) Depeyre, D.; Flicoteaux, C. Modeling of Thermal Steam Cracking on n-Hexadecane. Ind. Eng. Chem. Res. 1991, 30, 1116−1130. (2) Garner, S.; Sivaramakrishnan, R.; Brezinsky, K. The HighPressure Pyrolysis of Saturated and Unsaturated C7 Hydrocarbons. Proc. Combust. Inst. 2009, 32, 461−467. (3) Nageswara Rao, P.; Kunzru, D. Thermal Cracking of JP-10: Kinetics and Product Distribution. J. Anal. Appl. Pyrolysis 2006, 76, 154−160. (4) Sharma, S.; Raman, S.; Green, W. H. Intramolecular Hydrogen Migration in Alkylperoxy and Hydroperoxyalkylperoxy Radicals: Accurate Treatment of Hindered Rotors. J. Phys. Chem. A 2010, 114, 5689−5701. (5) Swisher, J. A.; Hansen, N.; Maesen, T.; Keil, F. J.; Smit, B.; Bell, A. T. Theoretical Simulation of n-Alkane Cracking on Zeolites. J. Phys. Chem. C 2010, 114, 10229−10239. (6) Viskolcz, B.; Lendvay, G.; Seres, L. Ab Initio Barrier Heights and Branching Ratios of Isomerization Reactions of a Branched Alkyl Radical. J. Phys. Chem. A 1997, 101, 7119−7127. (7) Stuart, S. J.; Tutein, A. B.; Harrison, J. A Reactive Potential for Hydrocarbons with Intermolecular Interactions. J. Chem. Phys. 2000, 112, 6472−6486. (8) Van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396−9409. (9) Shan, T.-R.; Devine, B. D.; Kemper, T. W.; Sinnott, S. B.; Phillpot, S. R. Charge-Optimized Many-Body Potential for the Hafnium/Hafnium Oxide System. Phys. Rev. B 2010, 81, 125328. (10) Ojwang, J. G. O.; van Santen, A., R.; Kramer, G. J.; van Duin, A. C. T.; Goddard, W. A., III. Predictions of Melting, Crystallization and Local Atomic Arrangements of Aluminum Clusters Using a Reactive Force Field. J. Chem. Phys. 2008, 129, 244506. (11) Jarvi, T. . T.; Kuronen, A.; Hakala, M.; Nordlund, K.; van Duin, A. C. T.; Goddard, W. A., III; Jacob, T. Development of a ReaxFF Description for Gold. Eur. J. Phys. B 2008, 66, 75. (12) Chenworth, K.; van Duin, A. C. T.; Goddard, W. A. ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A 2008, 112, 1040−1053.

computations using ReaxFF. With this new implementation, we can perform multiparticle reactive atomic scale simulations up to microseconds. The implemented event detector for identifying barrier crossing instances is flexible and can be applied to different chemical systems that obey first-order kinetics. We have shown that for most of the reactive systems, the efficiency of the PRD algorithm will increase by reducing the simulation temperature. As shown from heptene pyrolysis, extending the time scale of simulation has a significant effect on the simulation temperature and the reaction chemistry. By lowering the simulation temperature, we observed significant improvement in the product distributions and overall chemistry description that can be compared qualitatively with experiments. In addition, we can estimate kinetic parameters at temperatures that are much closer to experiments. RPRD efficiency will be very high for systems characterized by infrequent events, and it can be improved by carefully designing the algorithm for event detection. Caution must be exercised to ensure that the probability distribution follows a first order decay while discarding the low barrier events. Efforts are underway to explore the applicability of RPRD to capture the catalytic chemistry of methanol to hydrocarbon (MTH) conversion over ZSM-5 and SAPO catalyst while additional event detection schemes are being considered for investigation of porous zeolites before releasing the fixes to LAMMPS user community for extensive exploration of this accelerated algorithm.

■

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest. 3796

dx.doi.org/10.1021/jz4019223 | J. Phys. Chem. Lett. 2013, 4, 3792−3797

The Journal of Physical Chemistry Letters

Letter

(13) Van Duin, A. C. T.; Strachen, A.; Stewman, S.; Zhang, Q.; Xu, X. ReaxFF Reactive Force Field for Silicon and Silicon Oxide Systems. J. Phys. Chem. A 2003, 107, 3803−3811. (14) Chenworth, K.; van Duin, A. C. T.; Persson, P.; Cheng, M. J.; Oxgaard, J.; Goddard, W. A. Development and Application of a ReaxFF Reactive Force Field for Oxidation and Dehydrogenation on Vanadium Oxide Catalysts. J. Phys. Chem. C 2008, 112, 14645−14654. (15) Joshi, K.; van Duin, A. C. T.; Jacob, T. Development of a ReaxFF Description of Gold Oxides and Initial Application to Cold Welding of Partially Oxidezed Gold Surfaces. J. Mater. Chem. 2010, 20, 10431−10437. (16) Xu, L.; Henkelman, G. Adaptive Monte-Carlo for First Principles Accelerated Dynamics. J. Chem. Phys. 2008, 129, 114104. (17) Voter, A. F.; Montalenti, F.; Germann, T. C. Extending the Time Scale in a Atomistic Simulation of Materials. Annu. Rev. Mater. Res. 2002, 32, 321−346. (18) Perez, D.; Uberuaga, B. P.; Shim, Y.; Amar, J. G.; Voter, A. F. Accelerated Molecular Dynamics Methods: Introduction and Recent Developments. Annu. Rep. Comput. Chem. 2009, 5, 79−98. (19) Laio, A.; Gervasio, F. L. Metadynamics: A Method to Simulate Rare Events and Reconstruct the Free Energy in Biophysics, Chemistry and Material Science. Rep. Prog. Phys. 2008, 71, 126601. (20) Voter, A. Parallel Replica Method for Dynamics of Infrequent Events. Phys. Rev. B 1998, 57, R13985−R13988. (21) Duan, Y.; Halley, J. W.; Curtiss, L.; Redfern, P. Mechanisms of Lithium Transport in Amorphous Polyethylene Oxide. J. Chem. Phys. 2005, 122, 54702-1−8. (22) Warner, D. H.; Curtin, W. A.; Qu, S. Rate Dependence of Crack-Tip Processes Predicts Twinning Trends in FCC Metals. Nat. Mater. 2005, 6, 876−81. (23) Martini, A.; Dong, Y.; Perez, D.; Voter, A. F. Low-Speed Atomistic Simulation of Stick-Slip Friction Using Parallel Replica Dynamics. Tribol. Lett. 2009, 36, 63−68. (24) Uberuaga, B. P.; Stuart, S. J.; Voter, A. F. Parallel Replica Dynamics for Driven Systems: Derivation and Application to Strained Nanotubes. Phys. Rev. B 2007, 75, 014301-1−9. (25) Uberuaga, B. P.; Valone, S. M.; Baskes, M. I.; Voter, A. F. Accelerated Molecular Dynamics Study of Vacancies in Pu. AIP Conf. Proc. 2003, 673, 213−5. (26) Uberuaga, B. P.; Stuart, S. J.; Windl, W.; Masquelier, M. P.; Voter, A. F. Fullerene and Graphene Formation from Carbon Nanotube Fragments. Comput. Theor. Chem. 2012, 987, 115−121. (27) Kum, O.; Dickson, B. M.; Stuart, S. J.; Uberuaga, B. P.; Voter, A. F. Parallel Replica Dynamics with a Heterogeneous Distribution of Barriers: Application to n-Hexadecane Pyrolysis. J. Chem. Phys. 2004, 121, 9808−19. (28) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (29) Aktulga, H. M.; Fogarty, J. C.; Pandit, S. A.; Grama, A. Y. Parallel Reactive Molecular Dynamics: Numerical Methods and Algorithmic Techniques. Parallel Comput. 2012, 38, 245−259. (30) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (31) Ding, J.; Zhang, L.; Zhang, Y.; Han, K.-L. A Reactive Molecular Dynamics Study of n-Heptane Pyrolysis at High Temperature. J. Phys. Chem. A 2013, 117, 3266−3278. (32) Zehe, M. J.; Gordon, S.; McBride, B. J. CAP: A Computer Code for Generation of Tabular Thermodynamic Functions from NASA Lewis Coefficients; NASA: Hanover, MD, 2001; NASA TP-2001-210959. (33) Blades, A. T.; Sandhu, H. S. The Arrhenius Factors for Some Six-Center Unimolecular Reactions. Int. J. Kinet. 1971, III, 187−193.

3797

dx.doi.org/10.1021/jz4019223 | J. Phys. Chem. Lett. 2013, 4, 3792−3797