Direct observation of realistic-temperature fuel combustion mechanisms in atomistic simulations

Abstract

Atomistic simulations can in principle provide an unbiased description of all mechanisms, intermediates, and products of complex chemical processes. However, due to the severe time scale limitation of conventional simulation techniques, unrealistically high simulation temperatures are usually applied, which are a poor approximation of most practically relevant low-temperature applications. In this work, we demonstrate the direct observation at the atomic scale of the pyrolysis and oxidation of n-dodecane at temperatures as low as 700 K through the use of a novel simulation technique, collective variable-driven hyperdynamics (CVHD). A simulated timescale of up to 39 seconds is reached. Product compositions and dominant mechanisms are found to be strongly temperature-dependent, and are consistent with experiments and kinetic models. These simulations provide a first atomic-level look at the full dynamics of the complicated fuel combustion process at industrially relevant temperatures and time scales, unattainable by conventional molecular dynamics simulations.

---

Introduction

A detailed understanding of pyrolysis and combustion is of great technological and industrial importance. A fundamental insight in (bio)fuel decomposition chemistry is essential to improve the selectivity of cracking and reforming processes and increase the efficiency of combustion engines and minimize their production of pollutants. For example, low temperature combustion (LTC) strategies can significantly decrease production of particulate matter (PM) and nitrogen oxides (NO_x) in engines, but additional insights into their operation are required for further optimization.^1,2 To screen and improve possible operating conditions, kinetic modeling can be used to explain and guide experimental investigations.^3–6 It is, however, extremely challenging to create sufficiently complete and accurate kinetic models due to the wealth of possible intermediates and pathways that can all contribute significantly to the overall process, which generally limits their predictive power.

Atomistic simulation techniques can be used to bridge the gap between experimental results and kinetic models. Molecular dynamics (MD) simulations do not require any a priori knowledge of all possible reaction mechanisms and intermediates but generate the natural system evolution by explicitly integrating the equations of motions of all atoms. Therefore, MD simulations can be used to predict product compositions and to discover new unexpected pathways and intermediates, without any bias introduced by an incomplete reaction set. Not only can MD simulations be used to predict the outcome of a complex chemical process, but the thus obtained fundamental knowledge can also be used to extend and improve existing kinetic models. Crucial to the success of a MD simulation is the accuracy of the interatomic potential; in particular, the ReaxFF potential⁷ has been successfully applied to various pyrolysis and combustion reactions.^8–16 Nevertheless, a significant limitation of MD simulations is the short (up to nanosecond) time scale they are able to reach; previous MD studies therefore invariably used very high (>2000 K) simulated temperatures to be able to observe appreciable pyrolysis or combustion within the short MD time scale. The main drawback of this approach, however, is that it is difficult to correlate insights from high-temperature simulation with industrially relevant processes at lower temperatures, such as alkane cracking at ∼1000 K or low-temperature diesel engines. In order to reach these lower operating temperatures, the simulation time scale must be drastically extended.

Applying accelerated simulation methods to fuel decomposition is extremely challenging. Methods that require saddle-point searching, such as temperature-accelerated dynamics (TAD)¹⁷ or on-the-fly kinetic Monte Carlo,¹⁸ have difficulties handling liquid- or gas-like systems, whereas force-bias Monte Carlo simulations have been successful primarily in relaxing amorphous solids.¹⁹ The parallel replica (ParRep) method^20,21 imposes almost no constraints on the simulations and has been applied to the thermal decomposition of n-hexadecane²² and 1-hexene.²³ In the latter case, pyrolysis could be simulated at 1350 K over a simulated time of ∼1 μs by using up to 180 replicas. A further extension of the ParRep time scale to the millisecond-to-second range necessary for capturing processes at temperatures of 1000 K or lower is, however, impractical. Indeed, because the acceleration by ParRep is proportional to the number of processors, simulating this kind of process would put unrealistic demands on available computational resources.

In principle, longer time scales can be reached with hyperdynamics, at a much smaller cost.^24,25 This method operates by applying a bias potential ΔV to the potential energy surface, “filling” energy minima and consequently lowering the reaction activation energy. Designing a suitably general and efficient expression for ΔV is also the most challenging aspect of hyperdynamics. A practical problem of fuel pyrolysis and combustion simulations is the large separation of reaction barriers (and associated reaction time scales) that can be encountered during the process, ranging from ∼30 kcal mol⁻¹ for alkyl radical β-scissions to ∼80 kcal mol⁻¹ for initiation reactions of alkane pyrolysis. This has a major impact on the applicability of hyperdynamics, since a simple “static” bias potential can only be designed to work well for a small range of possible barriers; a bias that achieves a good acceleration or boost factor of β-scissions will still fall short in bringing the initiation reaction within reach. In some specific cases, a conventional hyperdynamics scheme can be sufficient: Cheng et al. exploited the very fast radical chemistry in hydrogen combustion, only applying a predefined bias potential to radical initiaton.²⁶ However, the much longer lifetimes of hydrocarbon radicals²³ and the employed ReaxFF-specific concepts render this approach not generally applicable.

In this work, we apply our recently proposed self-learning variant of the hyperdynamics algorithm, collective variable-driven hyperdynamics (CVHD) method²⁷ to the initial phase of n-dodecane pyrolysis and combustion to, for the first time, uncover detailed atomic-level fuel decomposition pathways under realistic conditions. These simulations are the first direct atomistic simulations of fuel pyrolysis and combustion chemistry under realistic conditions and provide an additional validation of contemporary mechanistic insights.

Computational methodology

The CVHD method

In the CVHD method,²⁷ which combines hyperdynamics with aspects of metadynamics,²⁸ a suitable bias potential can be slowly “grown” during the simulation until a transition is observed. A detailed discussion of the CVHD method and a comparison with other adaptive accelerated MD methods is available in ref. 27, but here we briefly summarize its main aspects.

Crucial to the success of a CVHD simulation is the choice of an appropriate collective variable (CV) that includes the relevant degrees of freedom s, and their distortions from equilibrium χ(s), associated with the to-be-boosted process. CVHD uses a general functional form that is inspired by the work of Tiwary and van de Walle,²⁹ and a generalization of the bond boost method.³⁰ Bias and system-specific dynamics are therefore cleanly separated: any complicated dynamics is projected on a single CV η as a value between 0 (no distortion) and 1 (maximal distortion receiving a bias), which is the only variable on which the bias explicitly depends. Furthermore, in contrast to the bond boost method, degrees of freedom other than bond elongations can be biased. For example, the folding of a model polymer was studied with CVHD by calculating η from dihedral angles rather than bond lengths.²⁷

As in metadynamics, a history-dependent bias potential is constructed by adding Gaussian-shaped “hills” w exp((η − η(t_i))²/2δ²) at intervals t_i. New hills are continuously added during the simulation to strengthen the bias, until a transition is observed. The criterion to detect a transition is time-based: if η remains 1 during a predefined waiting time t_w, the system is assumed to have undergone a transition. Then, the bias deposition procedure is reinitiated from scratch in the new state. No bias potential is added to the system when η = 1, so that the correct sequence of state-to-state transitions is preserved by construction.²⁴ CVHD is partially inspired by infrequent metadynamics, in which conventional metadynamics CVs are used but the bias deposition is made very slow in order to keep transition states relatively bias-free (instead of explicitly enforcing this, as is the case with CVHD's η).^31,32 Due to the different choice of CVs and biasing parameters, the two methods do not share the same focus: infrequent metadynamics can be used to calculate highly accurate rate estimates of a specified (slow) reaction (by repeated sampling of this transition), whereas CVHD is meant to capture the natural long time scale state-to-state evolution of the full system, discovering new reaction channels on-the-fly (and not necessarily sampling any encountered reaction beyond the first pass).

The combination of a bond length-based CV and an adaptive bias potential allows CVHD to handle complicated reactive processes with a wide distribution of barriers. As a first successful application to a chemical process, the CVHD method has already been used to simulate nickel-catalyzed methane decomposition, a process of which individual steps have barriers ranging from 8 to 32 kcal mol⁻¹, and time scales of several ps to ms at 800 K.²⁷

Simulation parameters

All simulations were carried out with LAMMPS³³ and the colvars module,³⁴ using the ReaxFF potential⁷ with the Chenoweth et al. parameter set⁸ and QEq charge equilibration,³⁵ as implemented in LAMMPS.³⁶ The equations of motion were integrated with a time step of 0.1 fs, and the system was initially equilibrated at the target temperature with a Langevin-type thermostat.³⁷ Further sampling in the NVT ensemble was achieved through application of a Nosé–Hoover chain³⁸ with a relaxation time of 0.1 ps, whereas for isotropic NPT simulations, the Martyna–Tobias–Klein (MTK) equations of motion³⁹ were integrated through the scheme of Tuckerman et al.,⁴⁰ using a relaxation time of 1 ps.

For pyrolysis simulations, the system consisted of 24 alkane molecules in a 50 × 50 × 50 ų periodic box, corresponding to a density of about 0.05 g cm⁻³. As local degrees of freedom we used C–C and C–H bond lengths, with as distortion function the strain χ_i = (r_i − r^min_i)/(r^max_i − r^min_i) for every bond, as described in ref. 27. The r^min_i and r^max_i parameters were respectively 1.55 and 2.20 Å for C–C bonds, and 1.05 and 1.65 Å for C–H bonds. The r^max_i values were specifically chosen to be smaller than the lengths of the breaking C–C and C–H bonds in the transition states of radical β-scissions and intramolecular hydrogen atom transfers, respectively, to ensure these states remain unbiased. This choice of CV means that only events involving bond breaking are accelerated, and conformational changes are unbiased; following a similar reasoning as in previous work, low-barrier conformational dynamics can be considered to have reached equilibrium well within the time spent while waiting for a reaction.²² Gaussian hills of width δ = 0.025 and height w = 0.25 kcal mol⁻¹ were added every t_i = 0.2 ps; the waiting time to detect events was t_w = 1 ps. CVHD simulations were carried out between 1000 and 1800 K; for comparison, unbiased MD simulations were conducted at a temperature of 2500 K.

Constant density combustion simulations were carried out for a 40 × 40 × 40 ų box containing 5 n-dodecane and 100 oxygen molecules, corresponding to a fuel-lean mixture with a density of about 0.1 g cm⁻³. The CVHD parameters are the same as those of the pyrolysis simulations, with all interactions involving oxygen atoms being described by the corresponding values for carbon. Biased simulations were carried out between 700 and 1800 K, and conventional MD was again performed at 2500 K. The average pressures in these simulations range from ∼200 bar at 700 K to almost 500 bar at 1800 K.

In order to capture the pressure dependence of the oxidation process over the range of pressures relevant to practical combustion applications, we also carried out a set of constant pressure simulations at 1000 K and pressures between 10 and 500 bar. A particularly important complication of CVHD simulation of gas-phase systems is that lowering the pressure also lowers the collision frequency in the system. Therefore, to prevent excessive buildup of bias between possible reactive collisions, and an overestimation of the time scale, the Gaussian deposition stride must be lowered accordingly. While 0.2 ps suffices for the high-density NVT simulations, we found that the 10 bar simulation requires a deposition interval of 0.5 ps, a value we used in all NPT simulations. When applying CVHD to other gas-phase systems, care must again be taken to choose an appropriate deposition stride.

Unless noted otherwise, all comparisons between simulations at different temperatures, such as of product compositions and time scales, are made at a fixed conversion level. For pyrolysis, analysis was performed at 50% fuel conversion. Combustion simulations were carried out until 20% of the O₂ molecules were consumed. For every condition, two independent trajectories were calculated to obtain reliable statistics. Error intervals, if reported, reflect the 90% confidence level.

Results and discussion

Accessible time scale

The dynamic self-learning nature of the CVHD method is illustrated in Fig. 1, which shows the evolution of the applied bias potential in the first stages of a pyrolysis simulation: the bias strength is slowly increased until an event is detected, and the biasing procedure is restarted. It can also be seen that the initiation reaction, which is a C–C bond fission, is the slowest event that requires the largest bias potential, whereas subsequent radical isomerizations and β-scissions have lower barriers. Thus, the bias strength is automatically tuned to be optimal for the current stage of the simulation. As summarized in Table 1, application of CVHD allows us to observe alkane pyrolysis and combustion at temperatures as low as 1000 and 700 K, respectively; the largest boost factor in our simulations is 8 × 10⁶ larger than that of the longest pyrolysis ParRep simulation.²³ The longest simulated physical time is therefore almost 40 s.

Fig. 1 Applied maximal bias potential during the initial steps of a 1000 K CVHD pyrolysis simulation. The time scales of the two displayed distinct regimes are also shown.

Table 1 Lowest temperatures achieved in the CVHD simulations of n-dodecane pyrolysis and combustion, and corresponding physical times and boost factors

	Pyrolysis	Combustion
Lowest temperature	1000 K	700 K
Longest simulated time	57 ms	39 s
Largest boost	6.3 × 10^6	1.3 × 10^9

---

Pyrolysis

In general, the alkane decomposition chemistry observed in the CVHD simulations is similar to previous high-temperature (>2000 K) MD simulations of alkane pyrolysis.^10,11 Most reactions of large alkyl radicals are either isomerization by intramolecular H-transfer, or decomposition to 1-alkenes through β-scission (the Rice–Kossiakoff mechanism). At high temperatures, the entropically favored decomposition reactions are the dominant process: ethylene is by far the dominant reaction product, in agreement with previous high-temperature MD simulations.^10,11 In contrast, low-barrier isomerization occurs much more frequently at low temperatures, forming more stable secondary radicals which give rise to the formation of larger 1-alkenes after eventually undergoing β-scission. Therefore, lower pyrolysis temperatures yield larger product molecules, as shown in Fig. 2. In contrast to the 2500 K simulation, where the C2 fraction is dominant and higher fractions are negligible, heavier molecules (C3 and higher) comprise about 50% of the products at 1000 K. Similarly, we observe that low-temperature propagation reactions involving H-abstraction by small radicals such as H, CH₃ and C₂H₅ constitute the main consumption channel of unreacted alkanes, but at high temperatures unimolecular initiation through bond fission gains importance.

Fig. 2 Products of the CVHD n-dodecane pyrolysis simulations at different temperatures.

The relative stability of C–C and C–H bonds is also found to be temperature-dependent. Because a C–H bond is about 25 kcal mol⁻¹ stronger than a C–C bond, unimolecular initiation at low temperatures only occurs through C–C dissociation; at high temperatures, considerable C–H dissociation is also observed, resulting in highly reactive free H atoms, in agreement with earlier high-temperature simulations of n-heptane pyrolysis.¹¹ A constant supply of free H radicals has a large impact on the overall reactivity of the system and the propagation rate, again illustrating the temperature-dependence of the pyrolysis mechanism. At low temperatures, C–H dissociation is only observed in radicals: ReaxFF predicts that C–H bonds vicinal to a radical site are about 50 kcal mol⁻¹ weaker than those in alkanes (dissociation energies of ∼50 and ∼100 kcal mol⁻¹, respectively), significantly facilitating their dissociation. Especially the ethyl radical, which has a C–H bond dissociation energy of 45 kcal mol⁻¹, frequently decomposes into C₂H₄ + H.

Oxidation

More complicated mechanisms are observed in the oxidation simulations, of which there are two distinct limiting cases, summarized in Fig. 3a. In the low temperature mechanism, the oxidation process is always initiated by hydrogen abstraction by an oxygen molecule and the subsequently formed alkyl radical combines with another oxygen molecule to form a peroxy radical ROO˙. Further isomerization leads to a hydroperoxyalkyl radical ˙QOOH, which can react further through a variety of pathways. Additionally, further H-abstractions by O₂ or reactive oxygen species from alkenes, radicals and carbonyl-containing compounds lead to the formation of compounds such as (conjugated) alkenes, ketenes and keto-hydroperoxides. At high temperatures, on the other hand, initial steps are essentially a pyrolysis process initiated by unimolecular C–C bond fission and subsequent β-scissions, forming primarily C₂H₄ which is further oxidized in a later stage.

Fig. 3 (a) Limiting temperature-dependent initial oxidation mechanisms, and (b) products of the constant-density CVHD n-dodecane oxidation simulations at different temperatures. Unox species do not contain oxygen, whereas Ox do; large products are C3 or heavier. The mass fraction is that of carbon only, and reflects how carbon is distributed over the various species.

At intermediate oxidation temperatures, both mechanisms are at play: below 1500 K, alkanes are initiated by H-abstraction but then easily break down into olefins, whereas from 1000 K and lower, C–C bond fission only rarely occurs in the initial oxidation stages. These temperature-dependent mechanisms are reflected by the product distributions of Fig. 3b. High temperatures primarily produce C₂H₄ and its oxidation products, whereas lowering the temperature suppresses dissociation events. In agreement with the findings of the pyrolysis simulations, alkyl radical β-scissions become less likely at lower temperatures, but the formed 1-alkenes are larger due to the relatively increased isomerization rate so that the mass fraction of produced hydrocarbons remains almost constant.

The temperature also has an impact on the formation of hydrogen peroxide and water. A first hydrogen atom transfer to O₂ forms a hydroperoxyl radical, HO₂, which can subsequently either abstract another hydrogen atom and form H₂O₂, or transfer its hydrogen atom to another radical. The further reactivity of H₂O₂ is strongly temperature-dependent, as it is found to be stable at low temperatures, whereas at high temperature, dissociation in two OH radicals occurs within a short time. These highly reactive OH radicals can then carry out an additional hydrogen abstraction to form H₂O. Therefore, at low temperatures, the kinetically stable H₂O₂ tends to accumulate whereas high temperatures favor the formation of OH and water. Indeed, at 700 K, the H₂O₂ fraction accounts for ∼17% of the non-O₂ oxygen atoms to be compared with ∼14% in the H₂O fraction. At 1000 K, this ratio is already 20/10 and from 1200 K onwards, the H₂O₂ fraction is negligible while the H₂O fraction contains about 30% of all reacted O₂. These observations are in agreement with conceptual models of low-temperature diesel engines.⁴¹

Combustion chemistry is also affected by pressure, and CVHD simulations can be used to investigate this effect. If the pressure-dependent reaction rate is proportional to pⁿ and, at constant temperature and assuming ideal gas behavior, the average reaction time 〈t〉 ∼ p/rate, the overall reaction order n can be determined by fitting ln〈t〉 = m ln p + ln〈t〉_p=1, in which m = 1 − n. This way, we obtained n = 2.07 ± 0.07, indicating that the rate-determining step of the oxidation is of second order, most likely involving hydrogen abstraction. Average oxidation time scales ranged from 0.6 ms at 500 bar, to 45 ms at 10 bar. The pressure effect on the relative importance of uni- and bimolecular processes is also reflected by the product distribution, as depicted in Fig. 4. Although this effect is less pronounced than the influence temperature has on the oxidation process, it can be seen that pyrolytic mechanisms are favored at low pressures, but suppressed in denser systems.

Fig. 4 Products of CVHD n-dodecane oxidation simulations at different pressures, at 1000 K. Presentation of the data is the same as in Fig. 3.

Comparison with experiments and unbiased MD

Our CVHD simulations also compare well with experimental results and existing kinetic models. The product distribution of the 1000 K pyrolysis process can be compared with a recent experimental study at the same temperature, in which the product distribution 0.08/0.44/0.23/0.25 of C1 through >C3 was obtained, in good agreement with our results.⁴² Moreover, the half-life of n-dodecane was found to be in the order of 20–40 ms, which compares well with the results in Table 1. The temperature-dependent oxidation mechanisms observed in CVHD simulations are also consistent with generally accepted models^3,4 and experiments.⁴²

There exists some discrepancy between our simulations and oxidation experiments. While experimentally, an early pyrolytic stage is already observed at 1050 K, our simulations suggest that this requires higher temperatures above 1200 K. This can be attributed to the high pressures in our constant density simulations, which will favor bimolecular over unimolecular reactions and thus a relative decrease of β-scissions over alkyl radical reactions with oxygen-containing species. Indeed, as shown earlier, lowering the pressure in our 1000 K CVHD simulation suppresses bimolecular reactions and gives rise to an early pyrolytic stage at lower temperatures than suggested by high-pressure simulations.

Finally, apparent first order rate constants for pyrolysis and combustion were computed from the C₁₂H₂₆ and O₂ consumption rates, respectively. By fitting the Arrhenius equation, prefactors A and activation energies E_A were obtained, and activation enthalpies Δ^‡H and entropies Δ^‡S were calculated from the Eyring equation, which are collected in Table 2. The pyrolysis parameters are consistent with other ReaxFF pyrolysis studies of n-dodecane, in which values of E_A between 56 and 66 kcal mol⁻¹ and A from 10¹⁵ to 10¹⁶ s⁻¹ are found,¹⁰ and with a unimolecular C–C dissociation as rate-determining step, as the positive entropy of activation indicates. For combustion, the activation energy matches that of a hydrogen abstraction by O₂. Indeed, experimental barriers of hydrogen atom transfers from alkanes to O₂ lie between 44 and 51 kcal mol⁻¹ (ref. 43) and ReaxFF predicts a barrier of ∼50 kcal mol⁻¹ for O₂-mediated hydrogen abstraction from methane.⁸ The negative Δ^‡S value for combustion is also in line with a bimolecular mechanism. This means that hydrogen atom transfers to O₂ are rate-determining at all temperatures, regardless of the different temperature-dependent initial reaction steps. Furthermore, as can be seen from the Arrhenius plots in Fig. 5, the CVHD values are also consistent with unbiased MD simulations: extrapolation of the CVHD results to higher temperatures agree with the MD results, therefore further validating the application of CVHD to pyrolysis and combustion.

Table 2 Kinetic parameters of n-dodecane pyrolysis and combustion as obtained from fitting apparent first order Arrhenius and Eyring equations

	Pyrolysis	Combustion
Temperature range (K)	1000–1800	700–1800
E A (kcal mol^−1)	70 ± 5	46 ± 1
A (s^−1)	5 × 10^15 to 2 × 10^17	7 × 10^11 to 3 × 10^12
Δ^‡H (kcal mol^−1)	68 ± 5	44 ± 1
Δ^‡S (cal mol^−1 K^−1)	12 ± 4	−8 ± 1

---

Fig. 5 Arrhenius plots of the apparent first order rate constants of n-dodecane pyrolysis and combustion as obtained from CVHD simulations. Filled symbols at 2500 K are unbiased MD simulations that were not included in the fit.

Conclusions

We have applied a recently developed self-learning hyperdynamics implementation, the CVHD method, to pyrolysis and combustion of the n-dodecane model fuel. Owing to the unprecedented long time scale of our simulations, we were able to conduct the first explicit verification of temperature- and pressure-dependent pyrolysis and combustion mechanisms through direct atomistic simulations. Reaction pathways uncovered by CVHD simulations agree well with experiments and kinetic models and suggests CVHD's ability to extend and supplement chemical kinetic models. Moreover, these results show that a flexible accelerated molecular dynamics method such as CVHD can give access to the long timescale dynamics of complex chemical processes, and how it can further extend the interpretive and predictive power of atomistic simulations by bridging the gap between theory and experiment.

Acknowledgements

K. M. B. is funded as PhD fellow (aspirant) of the FWO-Flanders (Fund for Scientific Research-Flanders), Grant 11V8915N. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center) and the HPC infrastructure of the University of Antwerp (CalcUA), funded by the Hercules Foundation and the Flemish Government – department EWI. The authors would also like to thank S. Banerjee for assisting with the interpretation of the experimental results.

References

1. M. Yao, Z. Zheng and H. Liu, Prog. Energy Combust. Sci., 2009, 35, 398–437.
2. R. D. Reitz and G. Duraisamy, Prog. Energy Combust. Sci., 2015, 46, 12–71.
3. F. Battin-Leclerc, E. Blurock, R. Bounaceur, R. Fournet, P.-A. Glaude, O. Herbinet, B. Sirjean and V. Warth, Chem. Soc. Rev., 2011, 40, 4762–4782.
4. C. K. Westbrook, W. J. Pitz, O. Herbinet, H. J. Curran and E. J. Silke, Combust. Flame, 2009, 156, 181–199.
5. K. Kohse-Höinghaus, P. Oßwald, T. Cool, T. Kasper, N. Hansen, F. Qi, C. Westbrook and P. Westmoreland, Angew. Chem., Int. Ed., 2010, 49, 3572–3597.
6. L. S. Tran, B. Sirjean, P.-A. Glaude, R. Fournet and F. Battin-Leclerc, Energy, 2012, 43, 4–18.
7. A. C. T. van Duin, S. Dasgupta, F. Lorant and W. A. Goddard III, J. Phys. Chem. A, 2001, 105, 9396–9409.
8. K. Chenoweth, A. C. T. van Duin and W. A. Goddard III, J. Phys. Chem. A, 2008, 112, 1040–1053.
9. K. Chenoweth, A. C. T. van Duin, S. Dasgupta and W. A. Goddard III, J. Phys. Chem. A, 2009, 113, 1740–1746.
10. Q.-D. Wang, J.-B. Wang, J.-Q. Li, N.-X. Tan and X.-Y. Li, Combust. Flame, 2011, 158, 217–226.
11. J. Ding, L. Zhang, Y. Zhang and K.-L. Han, J. Phys. Chem. A, 2013, 117, 3266–3278.
12. T. Qi, C. W. Bauschlicher, J. W. Lawson, T. G. Desai and E. J. Reed, J. Phys. Chem. A, 2013, 117, 11115–11125.
13. C. Zou, S. Raman and A. C. T. van Duin, J. Phys. Chem. B, 2014, 118, 6302–6315.
14. E. Salmon, A. C. T. van Duin, F. Lorant, P.-M. Marquaire and W. A. Goddard III, Org. Geochem., 2009, 40, 416–427.
15. F. Castro-Marcano, A. M. Kamat, M. F. Russo, A. C. T. van Duin and J. P. Mathews, Combust. Flame, 2012, 159, 1272–1285.
16. M. Zheng, X. Li, J. Liu and L. Guo, Energy Fuels, 2013, 27, 2942–2951.
17. M. R. Sørensen and A. F. Voter, J. Chem. Phys., 2000, 112, 9599–9606.
18. G. Henkelman and H. Jónsson, J. Chem. Phys., 2001, 115, 9657–9666.
19. K. M. Bal and E. C. Neyts, J. Chem. Phys., 2014, 141, 204104.
20. A. F. Voter, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, R13985–R13988.
21. D. Perez, B. P. Uberuaga and A. F. Voter, Comput. Mater. Sci., 2015, 100, 90–103.
22. O. Kum, B. M. Dickson, S. J. Stuart, B. P. Uberuaga and A. F. Voter, J. Chem. Phys., 2004, 121, 9808–9819.
23. K. L. Joshi, S. Raman and A. C. T. van Duin, J. Phys. Chem. Lett., 2013, 4, 3792–3797.
24. A. F. Voter, J. Chem. Phys., 1997, 106, 4665–4677.
25. A. F. Voter, Phys. Rev. Lett., 1997, 78, 3908–3911.
26. T. Cheng, A. Jaramillo-Botero, W. A. Goddard III and H. Sun, J. Am. Chem. Soc., 2014, 136, 9434–9442.
27. K. M. Bal and E. C. Neyts, J. Chem. Theory Comput., 2015, 11, 4545–4554.
28. A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 12562–12566.
29. P. Tiwary and A. van de Walle, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 87, 094304.
30. R. A. Miron and K. A. Fichthorn, J. Chem. Phys., 2003, 119, 6210–6216.
31. P. Tiwary and M. Parrinello, Phys. Rev. Lett., 2013, 111, 230602.
32. K. L. Fleming, P. Tiwary and J. Pfaendtner, J. Phys. Chem. A, 2016, 120, 299–305.
33. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
34. G. Fiorin, M. L. Klein and J. Hénin, Mol. Phys., 2013, 111, 3345–3362.
35. A. K. Rappé and W. A. Goddard III, J. Phys. Chem., 1991, 95, 3358–3363.
36. H. M. Aktulga, J. C. Fogarty, S. A. Pandit and A. Y. Grama, Parallel Computing, 2012, 38, 245–259.
37. G. Bussi and M. Parrinello, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 056707.
38. G. J. Martyna, M. L. Klein and M. Tuckerman, J. Chem. Phys., 1992, 97, 2635–2643.
39. G. J. Martyna, D. J. Tobias and M. L. Klein, J. Chem. Phys., 1994, 101, 4177–4189.
40. M. E. Tuckerman, J. Alejandre, R. López-Rendón, A. L. Jochim and G. J. Martyna, J. Phys. A: Math. Gen., 2006, 39, 5629.
41. M. P. B. Musculus, P. C. Miles and L. M. Pickett, Prog. Energy Combust. Sci., 2013, 39, 246–283.
42. S. Banerjee, R. Tangko, D. A. Sheen, H. Wang and C. T. Bowman, Combust. Flame, 2016, 163, 12–30.
43. W. Tsang, J. Phys. Chem. Ref. Data, 1990, 19, 1–68.

---