Extension of the ReaxFF Combustion Force Field toward Syngas

Article pubs.acs.org/JPCA

Extension of the ReaxFF Combustion Force Field toward Syngas Combustion and Initial Oxidation Kinetics Chowdhury Ashraf and Adri C.T. van Duin* Department of Mechanical and Nuclear Engineering, The Pennsylvania State University, University Park, Pennsylvania 16802, United States S Supporting Information *

ABSTRACT: A detailed insight of key reactive events related to oxidation and pyrolysis of hydrocarbon fuels further enhances our understanding of combustion chemistry. Though comprehensive kinetic models are available for smaller hydrocarbons (typically C3 or lower), developing and validating reaction mechanisms for larger hydrocarbons is a daunting task, due to the complexity of their reaction networks. The ReaxFF method provides an attractive computational method to obtain reaction kinetics for complex fuel and fuel mixtures, providing an accuracy approaching ab-initio-based methods but with a significantly lower computational expense. The development of the first ReaxFF combustion force field by Chenoweth et al. (CHO-2008 parameter set) in 2008 has opened new avenues for researchers to investigate combustion chemistry from the atomistic level. In this article, we seek to address two issues with the CHO-2008 ReaxFF description. While the CHO-2008 description has achieved significant popularity for studying large hydrocarbon combustion, it fails to accurately describe the chemistry of small hydrocarbon oxidation, especially conversion of CO2 from CO, which is highly relevant to syngas combustion. Additionally, the CHO-2008 description was obtained faster than expected H abstraction by O2 from hydrocarbons, thus underestimating the oxidation initiation temperature. In this study, we seek to systemically improve the CHO-2008 description and validate it for these cases. Additionally, our aim was to retain the accuracy of the 2008 description for larger hydrocarbons and provide similar quality results. Thus, we expanded the ReaxFF CHO-2008 DFT-based training set by including reactions and transition state structures relevant to the syngas and oxidation initiation pathways and retrained the parameters. To validate the quality of our force field, we performed high-temperature NVT-MD simulations to study oxidation and pyrolysis of four different hydrocarbon fuels, namely, syngas, methane, JP-10, and n-butylbenzene. Results obtained from syngas and methane oxidation simulation indicated that our redeveloped parameters (named as the CHO-2016 parameter set) has significantly improved the C1 chemistry predicted by ReaxFF and has solved the low-temperature oxidation initiation problem. Moreover, Arrhenius parameters of JP-10 decomposition and initiation mechanism pathways of n-butylbenzene pyrolysis obtained using the CHO-2016 parameter set are also in good agreement with both experimental and CHO-2008 simulation results. This demonstrated the transferability of the CHO-2016 description for a wide range of hydrocarbon chemistry.

1. INTRODUCTION The development of a comprehensive reaction mechanism and kinetic models for gas-phase oxidation of hydrocarbons has been an active area of research1−13 for the last couple of decades. Detailed insight of the key reactive events related pyrolysis and oxidation of hydrocarbon fuels will not only facilitate our understanding of combustion chemistry but also be helpful in developing comprehensive kinetic models. However, there are a couple of areas where the existing experimental and numerical techniques fail to offer important details. First, comprehensive kinetic models are available for smaller hydrocarbons only. In fact, for C3 and lower hydrocarbons we have detailed kinetic models1,7,9,10 with a maximum of 469 reactions and 71 species. With the increasing number of carbons the reaction network becomes more complex, making most of the reactions intractable. Recently, there have been some efforts to develop models for larger hydrocarbons,2,8,14−18 yet kinetic models for complex fuels and © 2017 American Chemical Society

fuel mixtures including jet surrogates are not available due to an increasingly complicated reaction network. Even if a model is developed, comprehensiveness of that model would be difficult to validate. Second, almost all of the available kinetic models are only applicable for low- or intermediate-pressure combustion, ignoring high-pressure and condensed-phase oxidation. This indicates the importance of a computationally feasible atomistic tool that not only can simulate the chemistry and physics involved in complex combustion dynamics but also requires no direct user input of possible chemical reactions. This tool should also be able to accurately predict energies and barriers for important chemical reactions. Additionally, it should be capable of handling relatively large systems for a longer time scale to generate reasonably accurate statistics. Received: December 9, 2016 Published: January 10, 2017 1051

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A

combustion dynamics, the first ReaxFF combustion description (CHO-2008)44 was developed by Chenoweth et al. in 2008. Initially, the CHO-2008 force field was tested to study hightemperature oxidation of methane, propane, o-xylene, and benzene. The correct reactivity trend and good agreement in reaction pathways with QM results were obtained. Furthermore, to evaluate the transferability of this force field, a followup study was performed to investigate JP-10 pyrolysis and initial oxidation pathways.43 Activation energy predicted by ReaxFF closely matched with experimental results; initial oxidation pathways also showed reasonable agreement with QM predictions. As the training set did not have data for any highly strained aliphatic molecule like JP-10, it clearly demonstrated the transferability of the ReaxFF CHO-2008 description in studying different hydrocarbon fuels. Since then it has been implemented in a wide range of applications in the scientific community for studying pyrolysis and oxidation of variety of fuels.42−44,46,61−68 These parameter sets are further extended to other elements including nitrogen/boron69 and various transition metals,70 which makes it a highly versatile atomistic-scale tool to study combustion dynamics. Thus, over the past few years, the ReaxFF CHO-2008 description has been used to study the pyrolysis and combustion of n-dodecane,67 1heptane,39 n-octanol,62 toluene,42 Illinois No. 6 coal,41 1,6dicyclopropane-2,4-hexyne,66 and liginin.61 The recent development of GPU-enabled ReaxFF71 has made it suitable to simulate large complex systems like coal72 and lignin68 pyrolysis. Apart from the hydrocarbon fuel, this description has also been used to investigate a wide range of aspects related to carbon-based materials including the oxidation of graphene,73 structural and chemical properties of graphene oxide,74 chemomechanics of crack propagation in graphene,75 and dynamics of carbon nano-onion formation.76 For all cases, generally a good agreement with experimental results in terms of the initiation mechanism and barriers were observed. Thus, the CHO-2008 description is considered to be an important tool to study highly reactive events during combustion through molecular dynamics (MD) simulation. For a more detailed review of this force field we refer to the recent paper by Dontgen et al.77 Despite its significant applicability in studying combustion chemistry in recent years, we identified three major areas where the CHO-2008 parameter set needs further refinement. First, though this parameter set gave a reasonable description for graphite and diamond in terms of chemical stability, it failed to describe their mechanical properties accurately. This is due to the fact that this description was never parametrized to describe the dynamics of carbon at condensed phase. Second, we will demonstrate here that though the majority of the applications of this description were focused to study larger hydrocarbons, it fails to properly describe the chemistry of smaller hydrocarbons, especially C1 chemistry. In particular, this description was found insufficient to describe the reactions involving CO and conversion of CO2 from CO. Lastly, we noticed that O2 is very reactive in this description, and it abstracts hydrogen from a hydrocarbon molecule at a lower than expected temperature. Since this is the most probable oxidation initiation reaction for most hydrocarbons, this issue deserves attention. The first issue was addressed by Srinivasan et al. (ReaxFF C-2013),78 where the ReaxFF training set was extended to describe carbon condensed phases more accurately. This provides am improved force field for studying mechanical property and defect chemistry for graphene.

Although quantum mechanical (QM) based methods are relatively accurate in predicting reaction energies and barriers, the combination of system size and time scale put a significant burden on those types of methods and are only limited to smaller systems (only couple of hundred atoms) within a shorter time scale. For example, Wang et al.19 simulated polymerization of acetylene and “Urey−Miller”20 experiment in an ab-initio nanoreactor using first-principle-based ab-initio molecular dynamics (AIMD). Their effort provided significant chemical insight; however, in spite of having a relatively small system size (156 atoms for acetylene and 228 atoms for Urey− Miller simulation), total computational cost for these simulations was very high (41 700 and 132 400 CPU and GPU hours for acetylene and Urey−Miller simulation, respectively). On-the-fly quantum chemistry calculations are also used in an automated reaction mechanism generator (RMG)21 to develop kinetic models for fuels such as butanol,22 JP-10,23 ketone biofuels,24 and neopentane;25 however, quantum chemistry calculations are only used to estimate thermodynamic parameters for different species. Alternative to QM methods, empirical methods like tight binding26 or force field (FF) based methods have been developed which are both computationally inexpensive and can provide a reasonable approximation of QM results. Since the FF-based methods can simulate a large system for a longer time scale, potentials such as DREIDING,27 MM3,28−30 MM4,31 COMPASS,32 etc., have been used to perform atomistic-scale simulations of hydrocarbon fuels. However, due to the rigid connectivity requirement, these potentials are nonreactive, i.e., unable to simulate chemical reactions. In addition, a number of reactive potentials have been developed, such as first-33 and second-generation34 reactive empirical bond order (REBO), charge-optimized many-body (COMB) potential,35,36 modified embedded-atom method (MEAM),37 and reactive force field (ReaxFF)38 which allow one to dynamically simulate bond formation and bond breaking. In particular, the ReaxFF has been used extensively39−46 to investigate complex combustion phenomena. Since ReaxFF solves Newton’s equation of motion instead of the more complex Schrodinger equation, it is a number of magnitudes faster than the QM-based methods. Additionally, ReaxFF aims to retain the accuracy of QM results in describing energies and barriers for chemical reactions as its parameters are trained against QM-based calculations. It includes connection-dependent terms in the force field description, where the interatomic potential is developed based on bond orders calculated directly from interatomic distances. Moreover, a polarizable charge description and bond order-dependent 3and 4-body terms are also included in the force field description; the combination of all these makes ReaxFF uniquely applicable for both metallic and covalent systems. Thus, ReaxFF can generate a fully dynamical description of chemical reactions related to combustion and other reactive systems, which has enabled the researchers to study previously inaccessible computational chemistry problems. Apart from studying combustion problems, the ReaxFF method has been extensively used to investigate a wide range of applications in materials,47−53 catalysts,54,55 and other chemical systems.56−60 Getting atomic insight of every minute detail during the combustion reaction is important as even a simple reaction that generates a short-lived species might prove significant in the overall reaction mechanism, which are practically intractable through experiments. To facilitate our understanding on 1052

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A

found in the previous literature.38,44 For recent reviews of the ReaxFF method, its applications to combustion, and related methods for large-scale simulations of chemical reactions, please refer to Senftle et al.81 and Dontgen et al.77 ReaxFF force field parameters are trained against QM-based density functional theory (DFT) or experimental data for reaction energies and transition states. We put all the relevant data in a training set and reparametrize the force field until a desired error level is achieved. The functional form of the error level that we optimize during parametrization is as follows

In this article, we extend the ReaxFF CHO-2008 potential to obtain a force field which can (i) accurately describe C1 chemistry, particularly reactions related to CO, (ii) correctly mimic the oxidation initiation reaction, (iii) properly describe carbon condense phases, and (iv) most importantly, retain the overall quality and accuracy of the results predicted by the CHO-2008 description for complex hydrocarbon pyrolysis and oxidation reactions. The newly developed description is then used to study the pyrolysis and oxidation of four hydrocarbon fuels−syngas, methane, JP-10, and butyl benzene. First, two fuels are chosen to illustrate the improved C1 chemistry captured by the new description, while the remaining two fuels allow us to demonstrate the quality of the new force field for complex hydrocarbon pyrolysis and how it compares with the previous description. Thus, the main objective of this study is to demonstrate the quality and transferability of the latest force field in studying pyrolysis and oxidation of varieties of hydrocarbon fuels. This paper is organized as follows. In section 2, the ReaxFF method is described along with the QM calculations done to train the ReaxFF parameters. Also, the simulation setups of different hydrocarbon systems are described in section 2. Results are discussed in section 3, while conclusions are summarized in section 4. Additional comparisons of energies between QM and ReaxFF are provided in the Supporting Information.

⎛ Xi , t − Xi ,ReaxFF ⎞2 Errori = ⎜ ⎟ σi ⎝ ⎠

where Xi,t and Xi,ReaxFF are the target quantum mechanical/ experimental and ReaxFF values of the ith entry of the training set and σi the inverse weight that determines how accurately the data needs to be fitted with the QM/experimental data. The summation of all errors in the training set provides the overall error, which is one of the measures of the quality of the force field. Details about the training are reported in section 3, while the force field parameters are supplied in the Supporting Information. From hereon, the reoptimized parameter sets will be called as ReaxFF CHO-2016 description. 2.2. Quantum Mechanical Calculation. ReaxFF parameters are optimized against a training set containing data from quantum mechanical (QM) calculation. To begin with, we combined the complete training sets of Chenoweth et al.44 and Srinivasan et al.78 This resultant training set ensures having data for both condensed-phase carbon and pyrolysis and oxidations of larger hydrocarbons. Readers are referred to these publications44,78 for details regarding the QM calculation. However, to summarize, gas-phase combustion data in the CHO-2008 training set were obtained using the B3LYP82,83 hybrid DFT functional and Popel 6-311G** basis set84 calculated in the Jaguar85 package. Equation of states of graphite and diamond were calculated in the VASP86 package using the Perdew−Burke−Ernzerhoff (PBE)87 exchange correlation with DFT-D288 parameters. In addition to the combined training set, since one of the main objectives of this work is to improve C1 chemistry, we include a number of reactions related to CO and HCO. The QM data for HCO + O2 reaction energies and barriers are taken from Hsu et al.89 Furthermore, we also add data for reaction energies and barriers for H2/O2 chemistry.40 Using an in-house python code, we generated all of the ReaxFF transition state structures related to our target reactions. We fed a probable transition state structure of each reaction calculated using low-temperature constrained MD simulation (known as the bond restraint method) in to our python code. It was made sure that our initial guess of the transition state structure has multiple high negative frequencies (vibrational frequency). Next, our python code systematically moves all atoms in a direction given by the eigenvectors of the second largest negative frequency. Doing that repeatedly gives us a structure that has only one high negative vibrational frequency, which is a characteristic feature of the transition state structure. 2.3. Molecular Dynamics Simulations. Once the force field is developed, we perform molecular dynamics (MD) simulations with four different fuels including C1 hydrocarbons (syngas and methane) and large hydrocarbons (JP-10 and nbutyl benzene). To obtain combustion dynamics, we performed NVT-MD simulations for each fuel. In MD simulations, NVT

2. COMPUTATIONAL METHOD 2.1. ReaxFF Reactive Force Field Method. In contrast to traditional, nonreactive force fields, ReaxFF can simulate bond formation and bond breaking during simulation. The connectivity between every pair of atoms is determined at each time step during the MD simulation, which enables ReaxFF to simulate reaction on-the-fly. ReaxFF is a bond order33,79 dependent interatomic potential where bond order is a function of interatomic distance between the pair of atoms. It uses the following expression to derive the forces on each atom Esystem = E bond + Eover + Eunder + E lp + E val + Etor + EvdWaals + Ecoulomb

(2)

(1)

where Ebond, Eover, Eunder, Elp, Eval, Etor, EvdWaals, and Ecoulomb represent bond energy, overcoordination energy penalty, undercoordination stability, lone pair energy, valence angle energy, torsion angle energy, van der Waals energy, and Coulomb energy, respectively. All the connectivity-dependent terms like bond, angle, and torsion are made bond order dependent, and ReaxFF functions are designed in such a way that their contribution diminishes upon bond breaking. ReaxFF also calculates the nonbonded interactions like Coulomb and van der Waals between every pair of atoms irrespective of their connectivity, and any excessive short-range nonbonded interactions are avoided by including a shielding term. ReaxFF calculates atomic charges using a geometry-dependent charge calculation scheme and uses the electronegativity equalization method (EEM)80 for this purpose. Additionally, to eliminate discontinuities in the nonbonded interaction energies and to reduce the range of the Coulomb interactions, a seventh-order Taper function is employed,38 which ensures that all nonbonded terms, together with their first, second, and third derivatives, go to zero at the nonbonded cutoff distance, which is typically picked to be 10 Å. A more detailed description of the ReaxFF energy terms can be 1053

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A

Figure 1. (a) Snapshot of the equilibrated syngas−O2 system; (b) final configuration after NVT-MD simulation is performed for 0.5 ns at 3000 K. Blue, red, and white spheres represent carbon, oxygen, and hydrogen atoms, respectively.

results. To do that, first, we minimized the system, and then the system was equilibrated for 25 ps at temperatures from 2000 to 2600 K with an interval of 100 K. After this equilibration, we chose 10 different starting configurations at each temperature and performed NVT-MD simulations for 50 ps. 2.3.4. n-Butylbenzene. Similar to JP-10, we investigated the initiation reaction related to n-butylbenzene pyrolysis. We put 24 n-butylbenzene molecules in a cubic box of 25 Å side length. NVT-MD simulation was performed at 2200 K for 25 ps after minimizing the system. Once the system is equilibrated, we chose 10 different initial configurations and performed longtime NVT-MD simulation at 2200 K. In all four cases, first a short (10 ps) NVT-MD simulation is performed to equilibrate the system. For oxidation cases, this is performed at low temperature with switched-off C−O and H− O bond parameters to prevent any reaction from occurring. This enables us to use a time step of 0.25 fs, which is typically larger than acceptable time steps of MD simulation. However, since the reactions usually occurred during the long-time NVT simulations, which can be referred to as the “production period”, it is important to choose the time step carefully. To capture the shortest vibrational frequency, usually a time step of 0.5-1 fs is required. However, as ReaxFF allows one to change charges and bond orders at every iteration, it requires even smaller time steps during the production period. Previous ReaxFF literature43,44 shows that a time step of 0.1 fs is small enough to capture all features in ReaxFF; thus, we used 0.1 fs as the time step during the production period. A bond order cutoff of 0.3 is used to identify different species during the simulation. The choice of bond order cutoff only effects the interpretation of the results without effecting the simulation itself. It should be noted carefully that the pressure of the simulations performed ranges from 0.1 to 0.5 GPA, which is very high compared to typical combustion conditions though not uncommon. For instance, high pressure can easily occur during rocket combustion and shock-induced processes. Usually in hydrocarbon combustion, ignition delay takes between microseconds to milliseconds. However, in reactive MD simulations, due to the computational cost we can only simulate up to couple of nanoseconds. Thus, in order to accelerate the process, we use higher pressure and temperature than the typical combustion conditions.

designates the condition at which simulations are performed with a constant number (N) of atoms at a constant volume (V). The temperature (T) of the system is also kept constant by using a thermostat. We typically use a Berendsen90 thermostat for this purpose with a damping constant of 100 fs. Details of the simulation for each system are as follows 2.3.1. Syngas. Syngas is a mixture of carbon monoxide and hydrogen gas. Since our main focus is to improve the C1 chemistry, especially CO oxidation reactions, syngas becomes an obvious choice for us to simulate. We put 40 CO, 40 H2, and 60 O2 in a cubic box of 20 Å side length. The number of molecules of different species made it a stoichiometric mixture. The system was first energy minimized, and then NVT-MD simulation was performed to equilibrate the system at 2000 K for 25 ps. After the equilibration, the system temperature was systematically increased from 2000 to 3000 K at a rate of 0.005 K/iteration. Once the system reached 3000 K, NVT-MD simulation was performed for a 500 ps. As all oxidation reactions took place during this period, it could be referred to as the production period. To get statistical data on the initiation reaction, we ran a series of 10 NVT-MD simulations, each of which had a unique starting configuration. Figure 1 shows a snapshot of such a system. Additionally, we simulated a similar system with the CHO2008 force field to compare the dynamics. 2.3.2. Methane. To simulate methane oxidation, we took the same system as reported in Chenoweth et al.44 We put 1 methane (CH4) molecule with 100 O2 molecules in a 25 Å cubic box. After minimization, we equilibrated the system using NVT-MD simulation at 2000 K for 25 ps; next, the temperature was ramped to 3000 K at a rate of 0.005 K/iteration. After that NVT-MD equilibrium was performed at 3000 K for 900 ps. Similar to syngas, we ran a series of 10 NVT-MD simulations to determine imitation reaction and reaction pathways of methane oxidation. Additionally, to compare the dynamics, we put 20 CH4 and 100 O2 molecules in a cubic box of 25 Å side length and performed ReaxFF NVT-MD simulations at 3000 K using both the CHO-2008 and the CHO-2016 descriptions. 2.3.3. JP-10. To study the pyrolysis of JP-10 fuel, we put 40 JP-10 molecules in a 36 Å × 36 Å × 18 Å box. Since the CHO2008 force field has already shown its capability to predict Arrhenius parameters43 related to the JP-10 initiation reaction, our aim is to obtain similar parameters using the CHO-2016 force field in order to see how they compare with previous 1054

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A

Figure 2. QM (6-311G**/B3LYP) and ReaxFF (CHO-2008 and CHO-2016) energies for (a) H−CO bond distance (the CHO-2008 description gives a barrier of around 15 kcal/mol, while CHO-2016 description gives a barrier of ∼4.5 kcal/mol, which is closer to the QM value (∼1.8 kcal/ mol)) and (b) H−C−O angle (HCO is a linear molecule in the CHO-2008 description, while it has an angle of 124.5° in CHO-2016 description (DFT value is 124°)).

3. RESULTS AND DISCUSSION 3.1. Modification of the ReaxFF Valence Angle Term. One of major objectives of reparameterization of the ReaxFF C/H/O force field is to accurately describe the C1 chemistry as the CHO-2008 description fails to do so. For example, in syngas combustion simulation, it cannot capture all possible pathways of CO2 formation, i.e., conversion of CO2 from CO, especially the pathways involving HCO radicals. HCO is arguably the most important radical in syngas combustion kinetics, which separates it from H2−O2 kinetics.91 The energy barrier of the reaction H + CO = HCO is very high in the CHO-2008 description (Figure 2a), which prohibits HCO formation during the simulation. CHO-2008 also makes it a linear molecule, while the QM calculation89 predicts it as angular with an H−C−O angle of ∼124° (Figure 2b). As this is the primary channel of HCO formation, a high barrier for this reaction causes a complete absence of HCO radicals during the simulation of syngas oxidation using the CHO-2008 description. Thus, an entire branch of CO2 conversion pathways from CO is missing in the simulations conducted using the CHO-2008 description of ReaxFF. During the initial stage of the force field training we notice that the valence angle formulation of ReaxFF makes HCO a linear molecule, which cannot be fixed by training only. For that reason, we modify the valence angle formulation of ReaxFF which was previously described44 as

angles within the ring structure of aromatic hydrocarbons; this leads to formation of a weak C−C bond between the distant C atoms. As this is not expected, we switch off this modification for the C−C−C angle only. With this modification in place and upon reparameterization we manage to get the desired energy profiles of the H−C−O angle and H−CO bond distance. Figure 2 compares the energy profiles from QM,89 CHO-2008, and the recently developed CHO-2016 force field for these two cases. It is evident that CHO-2016 agrees much better with the QM result than CHO2008 in terms of angle and bond distance. Also, it significantly reduces the energy barrier for the H + CO reaction, enhancing the possibility of producing HCO radicals during simulation. Details of the CHO-2016 force field training will be discussed in the next section. 3.2. Reparameterization of the ReaxFF Combustion Force Field. The quality of a molecular dynamics (MD) simulation vastly depends on the accuracy of the force field as the interactions between the atoms in the simulation are defined by the force field. A carefully scrutinized force field against available/calculated quantum mechanical (QM) and/or experimental results in a small representative system is required before simulating practical large-scale systems. In ReaxFF, a successive one-parameter search optimization technique92 is used to parametrize the force field against an extensive training set which contains quantum mechanical/experimental data of various quantities like bond lengths, valence angles, torsion angles, charges, heats of formation, equations of states, binding energies, diffusion pathways, and transition state calculations. The different parameters that require refitting are listed in a “params” file. This file is an essential input file for the force field optimization program, which in addition contains the maximum and minimum allowable values of different parameters and their search intervals. ReaxFF performs four individual optimization steps for each entry of the params file. In each run, the error for each entry of the training set is calculated in eq 2. During the four-step optimization process for each parameter, ReaxFF finds an optimized value of that parameter and visits the next parameter in the params file. As all parameters are strongly coupled in the force field, it is desirable to have multiple entries of the entire list of parameters so that each parameter can be optimized several times taking account the most recent values of other parameters. To derive ReaxFF parameters for combustion chemistry, we started with the parameters for the ReaxFF combustion force field (CHO-2008).44 However, since one of the main objectives

neighbors(j)

∑

SBO =

(BOjnπ + BOjnππ )

n=1 neighbors(j)

+ [1 −

∏

exp( −BOjn8 )]( −Δangle − pval8 . nlp , j) j

n=1

(3)

The equation is modified by removing the angle undercoordination term (Δjangle). The modified equation now has the following functional form neighbors(j)

SBO =

∑

(BOjnπ + BOjnππ )

n=1 neighbors(j)

+ [1 −

∏ n=1

exp( −BOjn8 )]( −pval8 . nlp , j)

(4)

Later, during the simulation of aromatic hydrocarbons we find that this modification allows one to make smaller C−C−C 1055

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A

Figure 3. QM (6-311G**/B3LYP) and ReaxFF CHO-2016 bond dissociation energies for (a) C ≡ O, (b) CO, (c) C−O, (d) O−O, (e) OO, and (f) O−H on representative molecules. Cyan, red, and white spheres represent carbon, oxygen, and hydrogen atoms, respectively.

Figure 4. QM (6-311G**/B3LYP) and ReaxFF CHO-2016 valence angle distortion energies for (a) OC−O angle, (b) O−C−O angle, (c) O CO angle, (d) O−O−O angle, (e) C−O−H angle, (f) C−O−O angle, (g) H−O−H angle, and (h) H−O−O angle. Cyan, red, and white spheres represent carbon, oxygen, and hydrogen atoms, respectively.

so that the resulting force field can describe both the solid and the gaseous phases of carbon simultaneously. As mentioned before, we include the training sets of both CHO-2008 and C2013 in our training. Additionally, to improve the chemistry related to smaller hydrocarbon combustion, we add bond dissociation energies for single, double, and triple bonds, angle bending energies for representative molecules, torsion energies,

of this reparameterization is to describe the mechanical properties of graphene and condensed phases of carbon properly, we took the C atom, C−C bond and nonbonded interactions, C−C−C valence angle, and C−C−C−C torsional angle parameters from ReaxFF C-2013.78 Though ReaxFF C2013 parameters were solely developed for carbon condensed phases, we started our force field training with those parameters 1056

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A

Figure 5. QM (6-311G**/B3LYP) and ReaxFF CHO-2016 dihedral angle energies for (a) O−C−O−O, (b) OC−O−O, (c) H−O−CO, and (d) H−C−O−O dihedrals. Cyan, red, and white spheres represent Carbon, oxygen, and hydrogen atoms, respectively.

Figure 6. Experimental and ReaxFF CHO-2016 (a) reaction energies and (b) energy barriers for H2−O2 subsystem and (c) reaction energies of CO oxidation reaction.

heat of formation, relative energies of a wide variety of representative molecules, transition state energies of different important reactions related to the H2−O2 subsystem, and CO/ HCO oxidation reactions. For bond dissociation, valence angle, and torsion angle, we use quantum mechanical (QM) data. However, we use experimental data of reaction energies and energy barriers for the H2−O2 subsystem and CO oxidation, since proper transition states in QM are difficult to identify and experimental data is readily available. Below we describe our training results related to C1 hydrocarbon combustion and carbon condensed phase only; remaining results can be found in the Supporting Information.

(a) For bond dissociation energies of various bonds containing C, H, and O atoms, we use the same QM data that was used in the CHO-2008 training set. All dissociation energies are calculated in multiple spin states (singlet/triplet or doublet/quartet). Since ReaxFF does not contain the concept of the spin state, it is only trained against the lowest energy spin state. In CHO2008, the authors obtained a C−O single bond dissociation energy curve using methanol as the reference molecule; similarly, they obtained a CO double bond dissociation energy curve from formaldehyde, CO triple bond dissociation energy from CO, and O−O 1057

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A

Figure 7. QM (6-311G**/B3LYP) and ReaxFF reaction energies and energy barriers for HCO + O2 reaction along with the structures identified in ReaxFF.

2016 force field can successfully reproduce the experimental values. (e) It is already mentioned that due to a very high reaction barrier, the ReaxFF CHO-2008 description cannot produce HCO radicals through the H + CO reaction. With valence angle modification and optimizing the required parameters, ReaxFF CHO-2016 now can describe the HCO formation chemistry accurately. Since HCO is an intermediate species to convert CO to CO2, we train our force field to successfully mimic the oxidation reactions of HCO, which produces either CO2 or CO as a product. Though the reactions look simple, QM study89 shows that the HCO + O2 reaction can assume a very complex route, and in order to produce CO or CO2, it can pass through four different transition states. CHO-2008 parameters seem to be good for the reaction energies related to these reactions; however, it cannot reproduce either the transition state structures or the energy barriers as observed in QM calculation. According to the QM calculation,89 when HCO and O2 reacts they make an intermediate “Int1” (Figure 7). Production of this intermediate species is a barrierless process; however, CHO-2008 gives a small barrier. Also, transition states TS3 and TS4 can only be observed if Int1 can pass through another transition state TS2 which has very high in energy in CHO-2008. Thus, these paths are highly unlikely to occur in the CHO-2008 description. Figure 7 shows that CHO-2016 description can reproduce the QM data, and the transition state structures and energies are in very good agreement with QM results. (f) Methane oxidation pathways using CHO-2008 have been reported by Chenoweth et al.,44 which is consistent with kinetic modeling studies. According to Chenoweth et al. the initiation reaction for methane oxidation is as follows

single and OO double bond dissociation energy from H2O2 and O2, respectively. Also, O−H single bond dissociation energy was calculated using the water molecule as the reference molecule. Figure 3 shows the result of our fitting, which indicates that the ReaxFF CHO-2008 description closely reproduces QM results for different bond dissociation energies. (b) We optimize valence angle parameters of different possible combinations of C/H/O against the QM data calculated from representative small molecules. Some cases also included the single- and double-bond combinations, e.g., O−C−O angle distortion energies were calculated for all three possible combinations: O− C−O, OC−O, and OCO. While optimizing the force field, only the angle of interest is constrained while the remaining structure is allowed to relax. Figure 4 shows a comparison between ReaxFF and QM for eight possible combinations of angles relevant to C 1 combustion chemistry. However, in total, 18 possible combinations of C, H, and O angles are included in the training set. (c) In the training set, we include 21 different torsion angle cases sufficient to describe all possible X−C−C−X, X− C−O−X, and X−O−O−X, which in turn can describe C−C, C−O, and O−O rotational barriers. We obtain the rotational barrier by constraining the appropriate dihedral angle and then relaxing the whole structure. Additionally, we also consider all possible combinations of single and double bonds around the central atoms. Figure 5 shows a comparison of ReaxFF and DFT for four dihedral cases related to C1 combustion chemistry. (d) Reaction energies and energy barriers related to hydrocarbon combustion are also added to our training set. We add all reaction energies for the H2−O 2 subsystem and CO oxidation mechanism. In either case, experimental values of reaction energies93 and energy barriers94 are used for training. Transition states are carefully calculated for the relevant cases using an inhouse python code. Figure 6 shows the ReaxFF CHO-

CH4 + O2 → CH3 + HO2

(R1)

Using a low-temperature (25 K) constrained MD simulation to drive the reaction coordinate, which is 1058

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A

Figure 8. (a) Comparison of minimum energy pathways (MEP) of reaction R3 using QM/experimental and ReaxFF descriptions where BR and sim represent bond restraint and simulation, respectively. (B) Transition state (TS) and stable intermediate (Int) structures obtained using different versions of ReaxFF description. Cyan, red, and white spheres represent carbon, oxygen, and hydrogen atoms, respectively.

Figure 9. Equation of state (EOS) data for (a) the isotropic variation of the lattice constant along graphite basal plane directions, (b) the variation of the lattice constant normal to the graphite basal plane directions, (c) the cubic expansion of a 2 × 2 × 2 diamond supercell containing 16 atoms, and (d) the uniaxial expansion of a 2 × 2 × 2 diamond supercell containing 16 atoms in the (001) direction.

the initiation reaction from real-time simulation and followed the trajectory of a CH4 and O2 molecule which took part in reaction at 2500 K MD simulation. We take several snapshots just before and after the reaction and calculate the vibrational frequency of each of the structures. This lead us to a structure which has only one strong negative frequency but much less energy than the transition state structure obtained by simply scanning or restraining the bond of interest (bond restrain (BR) method) using CHO-2008. If the reaction R1 passes through this transition state structure, it will have a barrier of ∼34 kcal/mol, whereas the enthalpy of reaction of reaction R1 is ∼38 kcal/mol.96 This monotonic increase of energy from reactant to product does not make this structure a transition state (TS) structure unless there is a complex stable intermediate between the TS and the products. With the help of an in-house python code which moves the atoms along the reaction coordinate, we manage to find a stable intermediate that has slightly lower energy than the TS structure. Further moving the atoms in the same direction gives us the products of the reaction. This confirms that the both the

typically known as the bond restraint method, they reported that this reaction has a barrier of ∼50 kcal/mol in the CHO-2008 description, which is reasonably close to the theoretical value of ∼52 kcal/mol.95 They performed 40 independent NVT-MD simulations of this system at 2500 K, and the initiation time of this reaction was below 1 ns for all cases they considered. They even observed this reaction at 2250 K, but the initiation time was much higher. However, a reaction with this high-energy barrier is less likely to occur at 2250 K in our molecular dynamics simulation time scale. Thus, we hypothesize that this faster than expected H abstraction by O2 is due to the presence of a favorable low-energy pathway in the CHO-2008 description. Since the low-temperature constrained MD biases the atoms to go through the desired directions to forcefully mimic the reaction, it is highly possible to move away from the reaction coordinate and the minimum energy pathway. In contrast, during the simulation, molecules can approach from any direction and take part in reaction. Thus, in order to prove our hypothesis of a low-energy pathway in CHO-2008, we moved our attention to track 1059

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A

perform a series of NVT-MD simulations of high-temperature (3000 K) oxidation of syngas (CO + H2). As mentioned in the computational details section, CO and H2 molecules in a 1:1 ratio are exposed to oxygen molecules with an equivalence ratio of 1. Additionally, to compare the initiation mechanism and combustion dynamics, we simulate the same system with the CHO-2008 force field. Figure 10 shows the distributions of

TS and the stable intermediate (Int1) exists along the minimum energy pathway (MEP) of reaction R1 when the CHO-2008 description is used. Figure 8a compares the MEP obtained through QM calculations with various ReaxFF cases. It clearly demonstrates that there exists a favorable low-energy pathway in CHO-2008 through which reaction R1 can occur during low-temperature simulation. Figure 8b shows the TS structures for various cases along with the stable intermediate (Int1) captured using CHO-2008. The stable intermediate Int1 has a complex structure where a H atom is highly overcoordinated (has a bond order of 1.76) as it makes a bond with both nearby C and O atoms. This high overcoordination then facilitates C− H bond breaking and results in faster than expected H abstraction by O2 at low temperatures. In contrast, CHO-2016 predicts a reaction barrier of ∼48 kcal/mol for this reaction. The energy barrier and TS structure for CHO-2016 are initially calculated using the bond restraint method, which is further confirmed by the inhouse python code. Additionally, simulations using CHO-2016 confirm the nonexistence of any other lowenergy pathway for reaction R1, which is further discussed in the next section. Since this is a highly relevant initiation reaction in hydrocarbon combustion, fixing this problem in CHO-2016 will lead toward more accurate prediction of the initial chemistry. (g) Additionally, we perform minor training on the carbon parameters obtained by Srinivasan et al. (C-2013) in order to capture the proper chemistry of aromatic hydrocarbons. We keep the C atom, C−C bond, and nonbonded interaction parameters the same as the C2013 description while only modifying the C−C−C angle and C−C−C−-C torsion parameters. Figure 9 shows the equation of state of graphene and diamond using QM, C-2013, and CHO-2016 parameters. The curves obtained from C-2013 and CHO-2016 are very similar, indicating that the newly developed parameters will be as accurate as C-2013 parameters to describe carbon condensed phases and mechanical properties. Furthermore, the CHO-2016 force field has been trained against the heat of formation of varieties of hydrocarbon molecules, reaction energies of relevant combustion reactions, and additional bond, valence angle and torsion angle energies of different combinations of C, H, and O atoms. A comparison of ReaxFF and DFT/experimental results for these cases is shown in the Supporting Information. Additionally, ReaxFF CHO2016 can reproduce the correct trend of the heat of formation of different C4H8O2 isomers (SI Figure 7, DFT data from Verstraelen et al.97). Moreover, the water density at room temperature in this description (1.05 g/cm3) is in reasonable agreement with experiment (1.00 g/cm3), whereas the density of water predicted by the CHO-2008 description was way higher (1.83 g/cm3) (SI Figure 8). Though this force field was never trained for room-temperature water, the density of water predicted by the CHO-2016 description is very promising. Please note, however, that this force field is not developed for liquid water reactivity simulationsfor such simulation targets we recommend using the ReaxFF parameter sets on the “water branch” as defined in Senftle et al.81 3.3. Molecular Dynamics Simulations of Hydrocarbon Pyrolysis/Oxidation. 3.3.1. Oxidation of Syngas. We

Figure 10. Distribution of reactants, products, and intermediate species during temperature ramp (25 ps) and NVT-MD equilibration simulation (500 ps) using (a) CHO-2008 and (b) CHO-2016 description. In the CHO-2008 description, H2−O2 chemistry starts at low temperature (during temperature ramp stage); however, for CHO2016 description, no reaction is observed during the temperature ramp.

reactants, products, and intermediate species during the temperature ramp and NVT-MD simulation at high temperature using both CHO-2008 and CHO-2016 descriptions. One representative case is chosen for each description. Significant differences are observed in both the initiation mechanism and combustion dynamics. Initiation Mechanism. As direct oxidation of CO has a very high energy barrier, the initiation mechanism for syngas combustion always starts with H2−O2 chemistry. For CHO2008, H2−O2 chemistry kicks off at a very low temperature (∼2100 K). Two different initiation reactions are observed at an almost equal ratio

1060

H 2 + O2 + M → H 2O2 + M; H 2O2 → 2OH

(R2)

H 2 + O2 → HO2 + H

(R3) DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A

Figure 11. Reaction network of CO2 conversion from CO in syngas simulation for (a) CHO-2008 and (b) CHO-2016.

Though according to the NIST kinetic database the reaction barrier for (R2) (∼50 kcal/mol) is slightly smaller than that of reaction R3 (∼60 kcal/mol), the trimolecular reaction R2 has a very slow rate.93 Thus, reaction R2 is less likely to initiate the H2−O2 chemistry at low temperatures (∼2100 K) in typical MD time scale. This leaves reaction R3 to be the most probable initiation reaction for syngas combustion as it produces the relatively stable HO2 radical. The CHO-2016 description in consistent with this observation as we observe reaction R3 to be the predominant initiation reaction. Additionally, no reaction is observed during the temperature ramp stage as the initiation reaction R3 has a high-energy barrier, which requires either high temperature or sufficient residence time to take place. Thus, this reaction is expected to take place at high temperature, and our simulation results with CHO-2016 observe this reaction during the initial stage of hightemperature (3000 K) NVT-MD simulation. Combustion Chemistry. As mentioned before, the H2−O2 chemistry starts at a much lower temperature for the CHO2008 force field. Almost 80% of hydrogen molecules convert to water during the temperature ramp stage, which is a final product of combustion. The intermediate species of the H2−O2 subsystem, i.e., OH and HO2 radicals, are frequently produced and consumed during this stage. However, the number of CO molecules remained constant during this period, indicating that no CO has been converted to CO2. In the case of the CHO2016 force field, due to the high-energy barrier of the initiation reaction, no reaction is observed during the temperature ramp process. Additionally, during the equilibration stage at 3000 K, different dynamics is observed for different force fields. Figure 10a shows that for the CHO-2008 force field CO2 production starts only after a significant water pool is build up, indicating that a high density of background water plays an important role in CO2 production. However, for the CHO-2016 force field, both CO2 and water start producing at a similar time, though water has a larger rate of production (Figure 10b). This indicates that though water helps in CO2 productiona wellknown phenomenon in syngas combustion98a dense water background is not required for the CHO-2016 description. This reveals the probable presence of additional reaction pathways of CO2 conversion from CO in the CHO-2016 description. In the case of hydrocarbon oxidation, it is universally observed that water production takes place much earlier than CO2. Almost all of the hydrocarbon combustion initiates either by C−C bond cleavage or by abstraction of a terminal hydrogen by O2. Once the system is initiated, the reactions are driven by intermediate radicals like HO2 and OH, which then finally convert to water. In contrast to water production, CO2

production is dependent upon production of C1 species like CO, HCO, or CH2O, which only happens after several reaction steps, especially for large hydrocarbons. However, in the case of syngas, CO is already in the system; thus, as soon as H2−O2 chemistry kicks off, CO can react with the intermediates OH and HO2. The reaction CO + OH is typically barrierless, while the reaction CO + HO2 has a small barrier (∼23 kcal/mol). Thus, both reactions can easily occur at high temperature and produce CO2. Accordingly, we expect to see the production of CO2 in syngas simulation soon after the system initiates through H2−O2 chemistry. The CHO-2016 description could capture this phenomenon, while the CHO-2008 description fails and shows a strong dependency on the presence of water molecules in the system. Though the initiation takes place at a much later time for CHO-2016, almost all H2 molecules are converted to H2O molecules after 500 ps of simulation. However, the final amount of CO2 molecules in the system is slightly higher for CHO-2016. This indicates the presence of additional pathways of CO2 formation in CHO-2016, which has converted more CO to CO2. One of the major aims for this force field reparameterization is to describe C1 chemistry properly. In the CHO-2008 description, the main reaction channel of CO2 formation is the barrierless reaction CO + OH = CO2 + H. Apart from this reaction, CO + HO2 = CO2 + OH and CO + O2 = CO2 + O reactions are also observed but are much less frequent. Most importantly, no HCO radical is formed during any stage of simulation, which is a very important intermediate in syngas combustion.91 This is quite expected as neither the geometry of the HCO radical nor the energy barrier of HCO formation was correct in CHO-2008. Consequently, HCO radical and associated reaction pathways are completely absent in the CHO-2008 description. However, when simulation is performed using the CHO2016 description, it showed much richer CO chemistry. Apart from capturing all three pathways that were captured by CHO2008, it also captures a whole new branch of reaction pathways involving HCO radicals. Also, the HCO radical is observed from time to time, though the maximum number of it never went over one during the simulation. This is due to the fact that the only way to generate a HCO radical in this simulation is the reaction H + CO, which has a small barrier (∼4 kcal/mol); however, H is a very short-lived radical, and whenever it produces during the simulation, it prefers to react with O2 to produce HO2 which is a barrierless reaction. A fuel-rich condition with more CO and less O2 might produce more HCO radicals, which we plan to test in the future. 1061

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A

CH3 and HO2 then recombine and produce CH3OOH, which then further produces CH3O and OH. This will only occur if the two radicals remain close to each other. However, due to their high kinetic energy at 3000 K, these two radicals tend to drift away from each other, and we see this recombination reaction only a couple of times with CHO-2016. Instead, in this description, the O−O bond of HO2 prefers to break at high temperature and produces an O and OH radical. This O radical then combines with CH3 and produces CH3O. Thus, the CHO-2016 description gives another route of CH 3 O formation. Once CH3O is formed, the next steps are similar in both force fields where conversion of species is as follows: CH3O → CH2O → HCO → CO, which is also in good agreement with the existing kinetic models.100−102 However, in CHO-2008, few species like [HCO(OOH)] and H−O−CO were observed as intermediates to convert CO from HCO, though in CHO-2016, HCO directly produces H and CO. As HCO is an unstable radical species, at 3000 K, the probability of the unimolecular reaction HCO → H + CO occurring is higher than the bimolecular reaction HCO + HO2 → [HCO(OOH)], which requires an effective collision. Figure 12 shows a probable reaction network of methane oxidation

In CHO-2016, once a HCO radical is produced, either it quickly falls apart to H + CO or it reacts with OH or HO2 radicals to produce CO2. These reactions are well established in the literature;13,94,99 however, the CHO-2008 description fails to reproduce them. Figure 11 compares the reaction network of CO2 formation from CO observed during the simulation using both force fields. It clearly demonstrates that the CHO-2016 description can capture a completely new branch of reactions involving HCO radical. Table 1 lists all reactions involving HCO radicals and CO2 production in the literature. CHO-2008 can capture only three Table 1. Reactions Listed in the Literature vs Reactions Observed during the Simulation

Figure 12. Methane oxidation reaction pathways as observed from CHO-2016. Final products from the simulation are marked as bold.

of them, while the CHO-2016 description can capture most of them during the simulation. Reactions that are listed in the table but are not observed during the simulation with the CHO-2016 description mostly involves two radicals. As the radicals are very short lived and the number of atoms is very small in our system, it is highly unlikely that two radicals will exist at same time and collide with one another. Thus, these reactions might be observed with a bigger system where a large number of atoms would lead to radical pool build up. 3.3.2. Oxidation of Methane. It has been described previously that the CHO-2008 description has a favorable low-energy pathway of the initiation reaction R1 for methane oxidation which leads to faster than expected H abstraction by O2. Chenoweth et al. performed 40 independent NVT-MD simulations, a system where 1 CH4 molecule is placed with 100 O2 molecules at 2500 K, and the initiation time of the reaction R1 was below 1 ns for all cases they considered. They even observed this reaction at 2250 K, but the initiation time was much higher. To compare, we perform 10 independent NVTMD simulations of the same system with CHO-2016 at 2500 K, and surprisingly, not a one of them is initiated through reaction R1 or any other reaction during 1 ns of simulation. This confirms that the low-energy pathway of H abstraction by O2 does not exist in the CHO-2016 description. On the basis of the above discussion, we increase our simulation temperature to 3000 K and run 10 independent NVT-MD simulations of 1 CH4 molecule exposed in 100 O2 molecules for 900 ps. Seven out of those 10 cases were initiated, and all of them were initiated through reaction R1. In the CHO-2008 description,

using the CHO-2016 description. Additionally, in only 1 case out of the 10 we studied, CO converted to CO2 by reacting with OH, which is a barrierless reaction. Since we observed this reaction frequently in syngas simulation, we can conclude than this would be the most probable path of CO2 formation. Next, to compare the dynamics of both force fields, we perform ReaxFF NVT-MD simulation of 20 CH4 and 100 O2 molecules at 3000 K. Figure 13 shows the distribution of fuel, intermediates, and products of the first 0.8 ns of simulation using both descriptions. Figure 13a clearly indicates that for CHO-2008, all CH4 are consumed within the first 0.15 ns of simulation. However, Figure 13b shows that for the CHO-2016 description the CH4 depletion rate is much slower and there are still a couple of CH4 molecules left after 0.8 ns. Since reaction R1 is the only initiation reaction observed in both descriptions, we can conclude that O2 is less reactive in the CHO-2016 description, which matches with our previous observation. Moreover, faster depletion of CH4 in the CHO2008 description leads to a CH2O peak (Figure 13a), which is absent in the CHO-2016 description (Figure 13b). The CH2O peak is a feature of low-temperature combustion of larger alkanes (starting from propane), which shows cool flame characteristics,103 and at this high temperature (3000 K) it should easily decompose to products like CO and H2.104 Thus, the formaldehyde peak should not be evident at hightemperature methane combustion, and the CHO-2016 description confirms this observation where formaldehyde formation and depletion is observed time to time. Additionally, 1062

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A

Figure 13. Distribution of fuel (CH4), intermediates, and products (CH2O, CO, CO2, and H2O) observed in ReaxFF NVT-MD simulation of 20 CH4 and 100 O2 molecules at 3000 K using (a) CHO-2008 and (b) CHO-2016 description.

Figure 14. (a) Rate of consumption of JP-10 molecules as a function of time at different temperatures. Results are obtained by averaging 10 independent simulations at each temperature. (b) Arrhenius plot for JP-10 pyrolysis. Rate constants at each temperature are obtained by using a linear fit for each curve in a. Best fit line with the equation is also shown.

the initial results obtained using the CHO-2016 description and compare that with previous results. Using the CHO-2016 description we observe two dominant primary decomposition pathways of JP-10 molecules, especially at lower temperatures. In the first pathway, ethylene decomposes from JP-10 molecules and produces a C8 intermediate, while in the second pathway, one JP-10 molecule decomposes to two C5 hydrocarbons. These two initial pathways are in good agreement with CHO-2008 results. However, at higher temperatures (>2300 K), a couple of additional decomposition pathways are observed. These pathways include either C−H bond scission to produce H radical or C7/C3 intermediates. In general, though the reaction network might be different depending on the choice of force field, both descriptions show similar primary decomposition pathways. Figure 14a shows the consumption of JP-10 molecules as a function of simulation time at various temperatures ranging from 2000 to 2600 K. To construct the curve at each temperature, we take the average of the number of JP-10 molecules consumed at each data point over 10 independent simulations. Linear fit of these data gives us the rate parameter at each temperature. Since JP-10 decomposition requires only a single molecule, these data can be used to study the first-order kinetics of JP-10 pyrolysis. In order to do so we calculate the Arrhenius parameters from the Arrhenius plot shown in Figure 14b and find the values of the activation energy (Ea) and preexponential factor (A) to be 56.80 kcal/mol and 25 × 1014 1/s, respectively. The values of these Arrhenius parameters are in

faster consumption of CH4 molecules in the CHO-2008 description leads to faster product formation, and all of the product (H2O and CO2) molecules are formed within the first 0.4 ns of simulation. However, in the CHO-2016 description, only 40% of CH4 molecules converted to CO2 within our simulation time limit. This demonstrates that less reactive O2 molecules in CHO-2016 not only reduces the CH4 depletion rate but also slow down the entire methane oxidation chemistry and delays the product formation. Since methane is less reactive due to a thermodynamically stable noble gas-like electron configuration and has no functional group and often catalysts are required to perform methane oxidation105,106 to get any useful product, we can conclude that it is a relatively slower process which is well captured in the ReaxFF CHO-2016 description. 3.3.3. Pyrolysis of JP-10. Previously, the CHO-2008 parameter set has been applied to study the initiation mechanism and kinetic parameters of large hydrocarbons. For JP-10, which is a single-component jet fuel, initiation mechanisms and Arrhenius parameters were obtained by Chenoweth et al.,43 which were in good agreement with experimental107 results. Since the development of a new ReaxFF C/H/O description also focuses on reproducing the success stories of the CHO-2008 description, JP-10 pyrolysis would be an obvious test case to consider. Thus, we perform the exact same set of simulations mentioned by Chenoweth et al.43 Though a complete reaction mechanism with statistical data is subjected to future study, here we would like to report 1063

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A

Experimental studies109,110 on pyrolysis of n-butylbenzene indicate that decomposition starts with Cα−Cβ bond scission, which produces a benzyl and a propyl radical. In fact, at low temperature, this is the only decomposition mechanism they observed. At higher temperature, additional decomposition pathways may occur involving other C−C bonds in the side chain. Since the C−H bond in stronger than the C−C bond, initiation through C−H bond scission is unlikely but not impossible at higher temperature. Since Cα−Cβ bond scission produces benzyl radical, which in turn produces acetylenea soot precursor111it is important for ReaxFF to capture this decomposition pathway to correctly describe soot formation chemistry. Thus, we perform 10 independent simulations of n-butylbenzene pyrolysis at 2200 K and identified the mechanism of initiation for each case. Sixty percent of the cases were initiated through Cα−Cβ bond scission, while Cβ−Cγ and Cγ−Cδ bond scissions are observed in 20% of the cases each (Figure 15b). Thus, Cα−Cβ bond scission is the dominant pathway of n-butylbenzene decomposition in the ReaxFF CHO-2016 description, which is in good accordance with experimental observations.109,110 This further validates the quality of this description to simulate large aromatic hydrocarbons.

very good agreement with experimental and previous ReaxFF calculations (Table 2). Table 2. Comparison of Arrhenius Parameters between Experiment and ReaxFF equation of best fit curve experimental107 CHO-200843 CHO-2016

y= + y= + y= +

−31443x 31.734 −29394x 36.041 −28532x 35.44

activation energy (Ea) (kcal/mol)

pre-exponential factor (A) (s−1)

62.4

0.57 × 1014

58.4

45 × 1014

56.8

25 × 1014

Table 2 shows that activation energy calculated using both versions of thr ReaxFF description lies within an acceptable range when compared with experiment. However, there is a 2 orders of magnitude difference between ReaxFF and experimental results when a pre-exponential factor is considered. Both ReaxFF descriptions predict a much higher value than the experimental one. The reason behind this discrepancy is probably the high pressure and temperature that we use to simulate the system. In a first-order reaction, the preexponential factor quantifies the number of molecules which attain enough energy to cross the energy barrier and decompose. High pressure means a higher number of molecules compressed at a smaller region; thus, more molecules will have enough energy to react. Also, the preexponential factor is considered constant over a small temperature range only; however, our simulation temperature is significantly above the experimental temperature. Combination of these two issues might be responsible for this higher value of the pre-exponential factor in ReaxFF. Finally, the decomposition pathways and Arrhenius parameters for JP-10 pyrolysis using the ReaxFF CHO-2016 description are in good agreement with the CHO-2008 version. This demonstrates the transferability of the newer version from smaller to larger hydrocarbons. 3.3.4. Pyrolysis of n-Butylbenzene. To further test the transferability of the recently developed C/H/O description toward larger hydrocarbons, we simulate n-butylbenzene, which is an important component of jet fuel surrogates.108 Additionally, this alkyl benzene simulation would lead us to evaluate the quality of the force field for aromatic hydrocarbons. Figure 15a shows a n-butylbenzene molecule; the carbon atoms in the side chain are named α-, β-, γ-, and δ-carbons, respectively, where Cα is the closest to the phenyl radical.

4. CONCLUSION We reparametrized the ReaxFF C/H/O force field in order to ensure its transferability to study the oxidation and pyrolysis for a wider range of hydrocarbon fuel and fuel mixtures. Since the previous C/H/O force field could only predict the combustion chemistry of larger hydrocarbons (>C3), we extended its capability to smaller hydrocarbons so that a single force field can be used to study any hydrocarbon fuel irrespective of their molecular size or structure. We validated the quality of the recently developed parameters by performing high-temperature oxidation and pyrolysis simulation of four different hydrocarbon fuels. Syngas simulation results revealed that the CHO2016 description provides much richer chemistry of CO2 conversion from CO, which is a set of important reactions in combustion chemistry to generate products of combustion. Previously, ReaxFF was not successful to capture all relevant pathways of this conversion reaction, therefore providing a less accurate distribution of species. Syngas simulation using the CHO-2016 description led to more CO2 formation than the CHO-2008 description within the same simulation time, and our analysis indicated that this was due to the presence of an additional intermediate −HCO in the system, which introduced a couple of additional pathways of CO2 formation. Though the experimental literature supports the presence of HCO radical and associated pathways, the CHO-2008 parameters were unable to capture them. Additionally, simulation of methane oxidation revealed the presence of a more favorable pathway of H abstraction by O2 during combustion using the CHO-2008 description; however, the experimental barrier for this type of reaction is typically much higher, which led to faster than expected H abstraction by O2 from a fuel molecule in ReaxFF. Since O2 was too reactive in CHO-2008, H abstraction by O2 was the dominate initiation mechanism for any hydrocarbon oxidation in ReaxFF. However, since the C−C bond is much weaker than the C−H bond, typically the initiation of a larger hydrocarbon starts by C−C bond cleavage and at lower temperature pyrolysis dominates over oxidation. As this issue is fixed in the recent development and now methane requires a much higher temperature to initiate in typical MD time scale

Figure 15. (a) n-Butylbenzene molecule. (b) Products of nbutylbenzene pyrolysis observed in ReaxFF simulation using CHO2016 description; values in parentheses represent the percentage of their occurrence. 1064

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A

(2) Pitz, W. J.; Westbrook, C. K.; Proscia, W. M.; Dryer, F. L. A Comprehensive Chemical Kinetic Reaction Mechanism for the Oxidation of N-Butane. Symp. Combust., [Proc.] 1985, 20 (1), 831− 843. (3) Westbrook, C. K.; Pitz, W. J.; Herbinet, O.; Curran, H. J.; Silke, E. J. A Comprehensive Detailed Chemical Kinetic Reaction Mechanism for Combustion of N-Alkane Hydrocarbons from NOctane to N-Hexadecane. Combust. Flame 2009, 156 (1), 181−199. (4) Pitz, W. J.; Westbrook, C. K. Chemical Kinetics of the High Pressure Oxidation of N-Butane and Its Relation to Engine Knock. Combust. Flame 1986, 63 (1−2), 113−133. (5) Mueller, M. A.; Yetter, R. A.; Dryer, F. L. Flow Reactor Studies and Kinetic Modeling of the H2/O2/NOX and CO/H2O/O2/NOX Reactions. Int. J. Chem. Kinet. 1999, 31 (10), 705−724. (6) SLOANE, T. M. Ignition and Flame Propagation Modeling with an Improved Propane Oxidation Mechanism. Combust. Sci. Technol. 1992, 83 (1−3), 77−96. (7) Dagaut, P.; Cathonnet, M.; Boettner, J.-C. Kinetic Modeling of Propane Oxidation and Pyrolysis. Int. J. Chem. Kinet. 1992, 24 (9), 813−837. (8) Chakir, A.; Bellimam, M.; Boettner, J. C.; Cathonnet, M. Kinetic Study of N-Heptane Oxidation. Int. J. Chem. Kinet. 1992, 24 (4), 385− 410. (9) Wilk, R. D.; Cernansky, N. P.; Pitz, W. J.; Westbrook, C. K. Propene Oxidation at Low and Intermediate Temperatures: A Detailed Chemical Kinetic Study. Combust. Flame 1989, 77 (2), 145−170. (10) Davis, S. G.; Law, C. K.; Wang, H. Propene Pyrolysis and Oxidation Kinetics in a Flow Reactor and Laminar Flames. Combust. Flame 1999, 119 (4), 375−399. (11) Shaw, R. Semi-empirical Extrapolation and Estimation of Rate Constants for Abstraction of H from Methane by H, O, HO, and O2. J. Phys. Chem. Ref. Data 1978, 7 (3), 1179−1190. (12) Chaos, M.; Dryer, F. L. Syngas Combustion Kinetics and Applications. Combust. Sci. Technol. 2008, 180 (6), 1053−1096. (13) Starik, A. M.; Titova, N. S.; Sharipov, A. S.; Kozlov, V. E. Syngas Oxidation Mechanism. Combust., Explos. Shock Waves 2010, 46 (5), 491−506. (14) Curran, H. J.; Gaffuri, P.; Pitz, W. J.; Westbrook, C. K. A Comprehensive Modeling Study of Iso-Octane Oxidation. Combust. Flame 2002, 129 (3), 253−280. (15) Allara, D. L.; Edelson, D. A Computational Modeling Study of the Low-Temperature Pyrolysis of N-Alkanes; Mechanisms of Propane, N-Butane, and N-Pentane Pyrolyses. Int. J. Chem. Kinet. 1975, 7 (4), 479−507. (16) Jahangirian, S.; Dooley, S.; Haas, F. M.; Dryer, F. L. A Detailed Experimental and Kinetic Modeling Study of N-Decane Oxidation at Elevated Pressures. Combust. Flame 2012, 159 (1), 30−43. (17) Autoignition Chemistry of the Hexane Isomers: An Experimental and Kinetic Modeling Study; http://papers.sae.org/ 952406. (18) Axelsson, E. I.; Brezinsky, K.; Dryer, F. L.; Pitz, W. J.; Westbrook, C. K. Chemical Kinetic Modeling of the Oxidation of Large Alkane Fuels: N-Octane and Iso-Octane. Symp. Combust., [Proc.] 1988, 21 (1), 783−793. (19) Wang, L.-P.; Titov, A.; McGibbon, R.; Liu, F.; Pande, V. S.; Martínez, T. J. Discovering Chemistry with an Ab Initio Nanoreactor. Nat. Chem. 2014, 6 (12), 1044−1048. (20) Miller, S. L.; Urey, H. C. Organic Compound Synthes on the Primitive Earth. Science 1959, 130 (3370), 245−251. (21) Gao, C. W.; Allen, J. W.; Green, W. H.; West, R. H. Reaction Mechanism Generator: Automatic Construction of Chemical Kinetic Mechanisms. Comput. Phys. Commun. 2016, 203, 212−225. (22) Harper, M. R.; Van Geem, K. M.; Pyl, S. P.; Merchant, S. S.; Marin, G. B.; Green, W. H. Comprehensive Reaction Mechanism for N-Butanol Pyrolysis and Combustion. Combust. Flame 2011, 158 (10), 2075−2075. (23) Gao, C. W.; Vandeputte, A. G.; Yee, N. W.; Green, W. H.; Bonomi, R. E.; Magoon, G. R.; Wong, H.-W.; Oluwole, O. O.; Lewis,

through H abstraction, it is now possible with CHO-2016 to define a pyrolysis-dominated or oxidation-dominated temperature region for different hydrocarbons and their mixtures. We plan to address this topic in the future. Additionally, we demonstrated that the new parameters are as good as the previous ones for large hydrocarbons by simulating the pyrolysis of JP-10 molecules. Arrhenius parameters and initiation mechanism of decomposition are similar for both descriptions, indicating the transferability of the recently developed parameters. Furthermore, CHO-2016 parameters have successfully predicted the initiation mechanism of n-butylbenzene pyrolysis involving Cα−Cβ bond scission, which further indicates its capability to simulate large aromatic hydrocarbons. Since we started our training using carbon parameters from C-2013 and only did a minor training for Cthe −C−C angle and C−C−C−C torsion angle parameters, this new parameter set can accurately predict the equation of state of graphite and diamond, similar to C-2013 parameters. This further enhances the transferability of this force field as the same set of parameters can be used to study both the solid and the gas phase of carbon. Results presented in this article indicate a significant increase of ReaxFF capability to study the combustion dynamics of varieties of hydrocarbon fuels along with mechanical properties of condensed-phase carbon, all with a single force field. The previous CHO-2008 ReaxFF version was good for studying large hydrocarbon oxidation; however, since the recent development includes a better description of smaller hydrocarbon oxidation chemistry, a complete and accurate mechanism leading to formation of final products such as CO2 and H2O would now be possible. Thus, ReaxFF could be used to generate a detailed kinetic mechanism of both simple and complex fuels from atomistic-scale simulations. Additionally, coupled with acceleration mechanisms like bond boost63 or parallel replica dynamics,112 it will open the doors to access areas like soot formation and soot oxidation, which was previously inaccessible using existing computational chemistry methods.

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.6b12429. Additional comparison of QM and ReaxFF energies (PDF) Force field parameters (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*Phone: 814 863 6277. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS This research was funded by the Air Force Office of Scientific Research (AFOSR) under grant no. FA9550-14-1-0355. REFERENCES

(1) Westbrook, C. K.; Pitz, W. J. A Comprehensive Chemical Kinetic Reaction Mechanism for Oxidation and Pyrolysis of Propane and Propene. Combust. Sci. Technol. 1984, 37 (3−4), 117−152. 1065

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A D. K.; Vandewiele, N. M.; et al. JP-10 Combustion Studied with Shock Tube Experiments and Modeled with Automatic Reaction Mechanism Generation. Combust. Flame 2015, 162 (8), 3115−3129. (24) Allen, J. W.; Scheer, A. M.; Gao, C. W.; Merchant, S. S.; Vasu, S. S.; Welz, O.; Savee, J. D.; Osborn, D. L.; Lee, C.; Vranckx, S.; et al. A Coordinated Investigation of the Combustion Chemistry of Diisopropyl Ketone, a Prototype for Biofuels Produced by Endophytic Fungi. Combust. Flame 2014, 161 (3), 711−724. (25) Petway, S. V.; Ismail, H.; Green, W. H.; Estupiñań , E. G.; Jusinski, L. E.; Taatjes, C. A. Measurements and Automated Mechanism Generation Modeling of OH Production in Photolytically Initiated Oxidation of the Neopentyl Radical. J. Phys. Chem. A 2007, 111 (19), 3891−3900. (26) Slater, J. C.; Koster, G. F. Simplified LCAO Method for the Periodic Potential Problem. Phys. Rev. 1954, 94 (6), 1498−1524. (27) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. DREIDING: A Generic Force Field for Molecular Simulations. J. Phys. Chem. 1990, 94 (26), 8897−8909. (28) Allinger, N. L.; Yuh, Y. H.; Lii, J. H. Molecular Mechanics. The MM3 Force Field for Hydrocarbons. 1. J. Am. Chem. Soc. 1989, 111 (23), 8551−8566. (29) Lii, J. H.; Allinger, N. L. Molecular Mechanics. The MM3 Force Field for Hydrocarbons. 2. Vibrational Frequencies and Thermodynamics. J. Am. Chem. Soc. 1989, 111 (23), 8566−8575. (30) Lii, J. H.; Allinger, N. L. Molecular Mechanics. The MM3 Force Field for Hydrocarbons. 3. The van Der Waals’ Potentials and Crystal Data for Aliphatic and Aromatic Hydrocarbons. J. Am. Chem. Soc. 1989, 111 (23), 8576−8582. (31) Allinger, N. L.; Chen, K.; Lii, J.-H. An Improved Force Field (MM4) for Saturated Hydrocarbons. J. Comput. Chem. 1996, 17 (5− 6), 642−668. (32) Sun, H. COMPASS: An Ab Initio Force-Field Optimized for Condensed-Phase ApplicationsOverview with Details on Alkane and Benzene Compounds. J. Phys. Chem. B 1998, 102 (38), 7338−7364. (33) Brenner, D. W. Empirical Potential for Hydrocarbons for Use in Simulating the Chemical Vapor Deposition of Diamond Films. Phys. Rev. B: Condens. Matter Mater. Phys. 1990, 42 (15), 9458−9471. (34) Brenner, D. W.; Shenderova, O. A.; Harrison, J. A.; Stuart, S. J.; Ni, B.; Sinnott, S. B. A Second-Generation Reactive Empirical Bond Order (REBO) Potential Energy Expression for Hydrocarbons. J. Phys.: Condens. Matter 2002, 14 (4), 783. (35) Yu, J.; Sinnott, S. B.; Phillpot, S. R. Charge Optimized ManyBody Potential for the &\mathrm{Si}/{\mathrm{SiO}}_{2}$ System. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75 (8), 85311. (36) Shan, T.-R.; Devine, B. D.; Hawkins, J. M.; Asthagiri, A.; Phillpot, S. R.; Sinnott, S. B. Second-Generation Charge-Optimized Many-Body Potential for $\text{Si}/{\text{SiO}}_{2}$ and Amorphous Silica. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82 (23), 235302. (37) Nouranian, S.; Tschopp, M. A.; Gwaltney, S. R.; Baskes, M. I.; Horstemeyer, M. F. An Interatomic Potential for Saturated Hydrocarbons Based on the Modified Embedded-Atom Method. Phys. Chem. Chem. Phys. 2014, 16 (13), 6233. (38) 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 (41), 9396−9409. (39) Castro-Marcano, F.; van Duin, A. C. T. Comparison of Thermal and Catalytic Cracking of 1-Heptene from ReaxFF Reactive Molecular Dynamics Simulations. Combust. Flame 2013, 160 (4), 766−775. (40) Agrawalla, S.; van Duin, A. C. T. Development and Application of a ReaxFF Reactive Force Field for Hydrogen Combustion. J. Phys. Chem. A 2011, 115 (6), 960−972. (41) Castro-Marcano, F.; Kamat, A. M.; Russo, M. F., Jr.; van Duin, A. C. T.; Mathews, J. P. Combustion of an Illinois No. 6 Coal Char Simulated Using an Atomistic Char Representation and the ReaxFF Reactive Force Field. Combust. Flame 2012, 159 (3), 1272−1285. (42) Cheng, X.-M.; Wang, Q.-D.; Li, J.-Q.; Wang, J.-B.; Li, X.-Y. ReaxFF Molecular Dynamics Simulations of Oxidation of Toluene at High Temperatures. J. Phys. Chem. A 2012, 116 (40), 9811−9818.

(43) Chenoweth, K.; van Duin, A. C. T.; Dasgupta, S.; Goddard, W. A., III Initiation Mechanisms and Kinetics of Pyrolysis and Combustion of JP-10 Hydrocarbon Jet Fuel. J. Phys. Chem. A 2009, 113 (9), 1740−1746. (44) Chenoweth, 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 (5), 1040−1053. (45) Zheng, M.; Li, X.; Liu, J.; Guo, L. Initial Chemical Reaction Simulation of Coal Pyrolysis via ReaxFF Molecular Dynamics. Energy Fuels 2013, 27 (6), 2942−2951. (46) Qian, H.-J.; van Duin, A. C. T.; Morokuma, K.; Irle, S. Reactive Molecular Dynamics Simulation of Fullerene Combustion Synthesis: ReaxFF vs DFTB Potentials. J. Chem. Theory Comput. 2011, 7 (7), 2040−2048. (47) LaBrosse, M. R.; Johnson, J. K.; van Duin, A. C. T. Development of a Transferable Reactive Force Field for Cobalt. J. Phys. Chem. A 2010, 114 (18), 5855−5861. (48) Nomura, K.; Kalia, R. K.; Nakano, A.; Vashishta, P.; van Duin, A. C. T.; Goddard, W. A. Dynamic Transition in the Structure of an Energetic Crystal during Chemical Reactions at Shock Front Prior to Detonation. Phys. Rev. Lett. 2007, 99 (14), 148303. (49) Ostadhossein, A.; Cubuk, E. D.; Tritsaris, G. A.; Kaxiras, E.; Zhang, S.; van Duin, A. C. Stress Effects on the Initial Lithiation of Crystalline Silicon Nanowires: Reactive Molecular Dynamics Simulations Using ReaxFF. Phys. Chem. Chem. Phys. 2015, 17 (5), 3832− 3840. (50) Islam, M. M.; Ostadhossein, A.; Borodin, O.; Yeates, A. T.; Tipton, W. W.; Hennig, R. G.; Kumar, N.; van Duin, A. C. ReaxFF Molecular Dynamics Simulations on Lithiated Sulfur Cathode Materials. Phys. Chem. Chem. Phys. 2015, 17 (5), 3383−3393. (51) Strachan, A.; van Duin, A. C. T.; Chakraborty, D.; Dasgupta, S.; Goddard, W. A. Shock Waves in High-Energy Materials: The Initial Chemical Events in Nitramine RDX. Phys. Rev. Lett. 2003, 91 (9), 98301. (52) Verners, O.; van Duin, A. C. T. Comparative Molecular Dynamics Study of Fcc-Ni Nanoplate Stress Corrosion in Water. Surf. Sci. 2015, 633, 94−101. (53) Wood, M. A.; Cherukara, M. J.; Kober, E. M.; Strachan, A. Ultrafast Chemistry under Nonequilibrium Conditions and the Shock to Deflagration Transition at the Nanoscale. J. Phys. Chem. C 2015, 119 (38), 22008−22015. (54) Zou, C.; Van Duin, A. Investigation of Complex Iron Surface Catalytic Chemistry Using the ReaxFF Reactive Force Field Method. JOM 2012, 64 (12), 1426−1437. (55) Shin, Y. K.; Kwak, H.; Vasenkov, A. V.; Sengupta, D.; van Duin, A. C. T. Development of a ReaxFF Reactive Force Field for Fe/Cr/O/ S and Application to Oxidation of Butane over a Pyrite-Covered Cr2O3 Catalyst. ACS Catal. 2015, 5 (12), 7226−7236. (56) Hong, S.; van Duin, A. C. T. Molecular Dynamics Simulations of the Oxidation of Aluminum Nanoparticles Using the ReaxFF Reactive Force Field. J. Phys. Chem. C 2015, 119 (31), 17876−17886. (57) Verlackt, C. C. W.; Neyts, E. C.; Jacob, T.; Fantauzzi, D.; Golkaram, M.; Shin, Y.-K.; van Duin, A. C. T.; Bogaerts, A. AtomicScale Insight into the Interactions between Hydroxyl Radicals and DNA in Solution Using the ReaxFF Reactive Force Field. New J. Phys. 2015, 17 (10), 103005. (58) Yeon, J.; van Duin, A. C. T. ReaxFF Molecular Dynamics Simulations of Hydroxylation Kinetics for Amorphous and Nano-Silica Structure, and Its Relations with Atomic Strain Energy. J. Phys. Chem. C 2016, 120 (1), 305−317. (59) Yue, D.-C.; Ma, T.-B.; Hu, Y.-Z.; Yeon, J.; van Duin, A. C.; Wang, H.; Luo, J. Tribochemical Mechanism of Amorphous Silica Asperities in Aqueous Environment: A Reactive Molecular Dynamics Study. Langmuir 2015, 31 (4), 1429−1436. (60) Zou, C.; Shin, Y. K.; van Duin, A. C. T.; Fang, H.; Liu, Z.-K. Molecular Dynamics Simulations of the Effects of Vacancies on Nickel Self-Diffusion, Oxygen Diffusion and Oxidation Initiation in Nickel, Using the ReaxFF Reactive Force Field. Acta Mater. 2015, 83, 102− 112. 1066

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A (61) Beste, A. ReaxFF Study of the Oxidation of Lignin Model Compounds for the Most Common Linkages in Softwood in View of Carbon Fiber Production. J. Phys. Chem. A 2014, 118 (5), 803−814. (62) Bharti, A.; Banerjee, T. Reactive Force Field Simulation Studies on the Combustion Behavior of N-Octanol. Fuel Process. Technol. 2016, 152, 132−139. (63) Cheng, T.; Jaramillo-Botero, A.; Goddard, W. A.; Sun, H. Adaptive Accelerated ReaxFF Reactive Dynamics with Validation from Simulating Hydrogen Combustion. J. Am. Chem. Soc. 2014, 136 (26), 9434−9442. (64) Guo, F.; Cheng, X.; Zhang, H. ReaxFF Molecular Dynamics Study of Initial Mechanism of JP-10 Combustion. Combust. Sci. Technol. 2012, 184 (9), 1233−1243. (65) He, Z.; Li, X.-B.; Liu, L.-M.; Zhu, W. The Intrinsic Mechanism of Methane Oxidation under Explosion Condition: A Combined ReaxFF and DFT Study. Fuel 2014, 124, 85−90. (66) Liu, L.; Bai, C.; Sun, H.; Goddard, W. A. Mechanism and Kinetics for the Initial Steps of Pyrolysis and Combustion of 1,6Dicyclopropane-2,4-Hexyne from ReaxFF Reactive Dynamics. J. Phys. Chem. A 2011, 115 (19), 4941−4950. (67) Wang, Q.-D.; Wang, J.-B.; Li, J.-Q.; Tan, N.-X.; Li, X.-Y. Reactive Molecular Dynamics Simulation and Chemical Kinetic Modeling of Pyrolysis and Combustion of N-Dodecane. Combust. Flame 2011, 158 (2), 217−226. (68) Zhang, T.; Li, X.; Qiao, X.; Zheng, M.; Guo, L.; Song, W.; Lin, W. Initial Mechanisms for an Overall Behavior of Lignin Pyrolysis through Large-Scale ReaxFF Molecular Dynamics Simulations. Energy Fuels 2016, 30 (4), 3140−3150. (69) Weismiller, M. R.; Duin, A. C. T. v.; Lee, J.; Yetter, R. A. ReaxFF Reactive Force Field Development and Applications for Molecular Dynamics Simulations of Ammonia Borane Dehydrogenation and Combustion. J. Phys. Chem. A 2010, 114 (17), 5485−5492. (70) Mueller, J. E.; van Duin, A. C. T.; Goddard, W. A. Development and Validation of ReaxFF Reactive Force Field for Hydrocarbon Chemistry Catalyzed by Nickel. J. Phys. Chem. C 2010, 114 (11), 4939−4949. (71) Kylasa, S. B.; Aktulga, H. M.; Grama, A. Y. PuReMD-GPU: A Reactive Molecular Dynamics Simulation Package for GPUs. J. Comput. Phys. 2014, 272, 343−359. (72) Zheng, M.; Li, X.; Liu, J.; Wang, Z.; Gong, X.; Guo, L.; Song, W. Pyrolysis of Liulin Coal Simulated by GPU-Based ReaxFF MD with Cheminformatics Analysis. Energy Fuels 2014, 28 (1), 522−534. (73) Goverapet Srinivasan, S.; van Duin, A. C. T. MolecularDynamics-Based Study of the Collisions of Hyperthermal Atomic Oxygen with Graphene Using the ReaxFF Reactive Force Field. J. Phys. Chem. A 2011, 115 (46), 13269−13280. (74) Bagri, A.; Mattevi, C.; Acik, M.; Chabal, Y. J.; Chhowalla, M.; Shenoy, V. B. Structural Evolution during the Reduction of Chemically Derived Graphene Oxide. Nat. Chem. 2010, 2 (7), 581−587. (75) Huang, X.; Yang, H.; van Duin, A. C. T.; Hsia, K. J.; Zhang, S. Chemomechanics Control of Tearing Paths in Graphene. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 85 (19), 195453. (76) Ganesh, P.; Kent, P. R. C.; Mochalin, V. Formation, Characterization, and Dynamics of Onion-like Carbon Structures for Electrical Energy Storage from Nanodiamonds Using Reactive Force Fields. J. Appl. Phys. 2011, 110 (7), 73506. (77) Döntgen, M.; Przybylski-Freund, M.-D.; Kröger, L. C.; Kopp, W. A.; Ismail, A. E.; Leonhard, K. Automated Discovery of Reaction Pathways, Rate Constants, and Transition States Using Reactive Molecular Dynamics Simulations. J. Chem. Theory Comput. 2015, 11 (6), 2517−2524. (78) Srinivasan, S. G.; van Duin, A. C. T.; Ganesh, P. Development of a ReaxFF Potential for Carbon Condensed Phases and Its Application to the Thermal Fragmentation of a Large Fullerene. J. Phys. Chem. A 2015, 119 (4), 571−580. (79) Tersoff, J. Empirical Interatomic Potential for Carbon, with Applications to Amorphous Carbon. Phys. Rev. Lett. 1988, 61 (25), 2879−2882.

(80) Mortier, W. J.; Ghosh, S. K.; Shankar, S. ElectronegativityEqualization Method for the Calculation of Atomic Charges in Molecules. J. Am. Chem. Soc. 1986, 108 (15), 4315−4320. (81) Senftle, T. P.; Hong, S.; Islam, M. M.; Kylasa, S. B.; Zheng, Y.; Shin, Y. K.; Junkermeier, C.; Engel-Herbert, R.; Janik, M. J.; Aktulga, H. M.; et al. The ReaxFF Reactive Force-Field: Development, Applications and Future Directions. Npj Comput. Mater. 2016, 2, 15011. (82) Becke, A. D. A New Mixing of Hartree−Fock and Local Density-functional Theories. J. Chem. Phys. 1993, 98 (2), 1372−1377. (83) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37 (2), 785−789. (84) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. Selfconsistent Molecular Orbital Methods. XX. A Basis Set for Correlated Wave Functions. J. Chem. Phys. 1980, 72 (1), 650−654. (85) Jaguar; Schrö dinger, https://www.schrodinger.com/jaguar (accessed Sep 30, 2016). (86) Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54 (16), 11169. (87) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77 (18), 3865− 3868. (88) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27 (15), 1787−1799. (89) Hsu, C.-C.; Mebel, A. M.; Lin, M. C. Abinitio Molecular Orbital Study of the HCO+O2 Reaction: Direct versus Indirect Abstraction Channels. J. Chem. Phys. 1996, 105 (6), 2346−2352. (90) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81 (8), 3684−3690. (91) Boivin, P.; Jiménez, C.; Sánchez, A. L.; Williams, F. A. A FourStep Reduced Mechanism for Syngas Combustion. Combust. Flame 2011, 158 (6), 1059−1063. (92) van Duin, A. C.; Baas, J. M.; van de Graaf, B. Delft Molecular Mechanics: A New Approach to Hydrocarbon Force Fields. Inclusion of a Geometry-Dependent Charge Calculation. J. Chem. Soc., Faraday Trans. 1994, 90 (19), 2881−2895. (93) Baulch, D. L.; Bowman, C. T.; Cobos, C. J.; Cox, R. A.; Just, T.; Kerr, J. A.; Pilling, M. J.; Stocker, D.; Troe, J.; Tsang, W.; et al. Evaluated Kinetic Data for Combustion Modeling: Supplement II. J. Phys. Chem. Ref. Data 2005, 34 (3), 757−1397. (94) Li, J.; Zhao, Z.; Kazakov, A.; Chaos, M.; Dryer, F. L.; Scire, J. J. A Comprehensive Kinetic Mechanism for CO, CH2O, and CH3OH Combustion. Int. J. Chem. Kinet. 2007, 39 (3), 109−136. (95) Shaw, R. Semi-empirical Extrapolation and Estimation of Rate Constants for Abstraction of H from Methane by H, O, HO, and O2. J. Phys. Chem. Ref. Data 1978, 7 (3), 1179−1190. (96) Zhu, R.; Lin, C. The CH3 + HO2 Reaction: First-Principles Prediction of Its Rate Constant and Product Branching Probabilities. J. Phys. Chem. A 2001, 105 (25), 6243−6248. (97) Verstraelen, T.; Vandenbrande, S.; Ayers, P. W. Direct Computation of Parameters for Accurate Polarizable Force Fields. J. Chem. Phys. 2014, 141 (19), 194114. (98) Kim, T. J.; Yetter, R. A.; Dryer, F. L. New Results on Moist CO Oxidation: High Pressure, High Temperature Experiments and Comprehensive Kinetic Modeling. Symp. Combust., [Proc.] 1994, 25 (1), 759−766. (99) Mueller, M. A.; Yetter, R. A.; Dryer, F. L. Flow Reactor Studies and Kinetic Modeling of the H2/O2/NOX and CO/H2O/O2/NOX Reactions. Int. J. Chem. Kinet. 1999, 31 (10), 705−724. (100) de Joannon, M.; Sabia, P.; Tregrossi, A.; Cavaliere, A. Dynamic Behavior of Methane Oxidation in Premixed Flow Reactor. Combust. Sci. Technol. 2004, 176 (5−6), 769−783. (101) Hidaka, Y.; Sato, K.; Henmi, Y.; Tanaka, H.; Inami, K. ShockTube and Modeling Study of Methane Pyrolysis and Oxidation. Combust. Flame 1999, 118 (3), 340−358. 1067

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068

Article

The Journal of Physical Chemistry A (102) Petersen, E. L.; Davidson, D. F.; Hanson, R. K. Kinetics Modeling of Shock-Induced Ignition in Low-Dilution CH4/ O2Mixtures at High Pressures and Intermediate Temperatures. Combust. Flame 1999, 117 (1−2), 272−290. (103) Sheinson, R. S.; Williams, F. W. Chemiluminescence Spectra from Cool and Blue Flames: Electronically Excited Formaldehyde. Combust. Flame 1973, 21 (2), 221−230. (104) Skinner, G. B.; Lifshitz, A.; Scheller, K.; Burcat, A. Kinetics of Methane Oxidation. J. Chem. Phys. 1972, 56 (8), 3853−3861. (105) Periana, R. A.; Taube, D. J.; Gamble, S.; Taube, H.; Satoh, T.; Fujii, H. Platinum Catalysts for the High-Yield Oxidation of Methane to a Methanol Derivative. Science 1998, 280 (5363), 560−564. (106) Fujimoto, K.; Ribeiro, F. H.; Avalos-Borja, M.; Iglesia, E. Structure and Reactivity of PdOx/ZrO2Catalysts for Methane Oxidation at Low Temperatures. J. Catal. 1998, 179 (2), 431−442. (107) Nageswara Rao, P.; Kunzru, D. Thermal Cracking of JP-10: Kinetics and Product Distribution. J. Anal. Appl. Pyrolysis 2006, 76 (1− 2), 154−160. (108) Colket, M. B.; Seery, D. J. Reaction Mechanisms for Toluene Pyrolysis. Symp. Combust., [Proc.] 1994, 25 (1), 883−891. (109) Diévart, P.; Dagaut, P. The Oxidation of N-Butylbenzene: Experimental Study in a JSR at 10 Atm and Detailed Chemical Kinetic Modeling. Proc. Combust. Inst. 2011, 33 (1), 209−216. (110) Freund, H.; Olmstead, W. N. Detailed Chemical Kinetic Modeling of Butylbenzene Pyrolysis. Int. J. Chem. Kinet. 1989, 21 (7), 561−574. (111) Sivaramakrishnan, R.; Tranter, R. S.; Brezinsky, K. High Pressure Pyrolysis of Toluene. 2. Modeling Benzyl Decomposition and Formation of Soot Precursors. J. Phys. Chem. A 2006, 110 (30), 9400− 9404. (112) Joshi, K. L.; Raman, S.; van Duin, A. C. T. Connectivity-Based Parallel Replica Dynamics for Chemically Reactive Systems: From Femtoseconds to Microseconds. J. Phys. Chem. Lett. 2013, 4 (21), 3792−3797.

1068

DOI: 10.1021/acs.jpca.6b12429 J. Phys. Chem. A 2017, 121, 1051−1068