Automated Chemical Kinetic Modeling via Hybrid Reactive Molecular

Article pubs.acs.org/jcim

Cite This: J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Automated Chemical Kinetic Modeling via Hybrid Reactive Molecular Dynamics and Quantum Chemistry Simulations Malte Döntgen,†,‡ Felix Schmalz,† Wassja A. Kopp,† Leif C. Kröger,† and Kai Leonhard*,† †

Chair of Technical Thermodynamics, RWTH Aachen University, 52062 Aachen, Germany Molecular Science, Department of Chemistry, University of Helsinki, 00560 Helsinki, Finland

‡

Downloaded via KAOHSIUNG MEDICAL UNIV on July 12, 2018 at 12:10:16 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: An automated scheme for obtaining chemical kinetic models from scratch using reactive molecular dynamics and quantum chemistry simulations is presented. This methodology combines the phase space sampling of reactive molecular dynamics with the thermochemistry and kinetics prediction capabilities of quantum mechanics. This scheme provides the NASA polynomial and modified Arrhenius equation parameters for all species and reactions that are observed during the simulation and supplies them in the ChemKin format. The ab initio level of theory for predictions is easily exchangeable, and the presently used G3MP2 level of theory is found to reliably reproduce hydrogen and methane oxidation thermochemistry and kinetics data. Chemical kinetic models obtained with this approach are ready to use for, e.g., ignition delay time simulations, as shown for hydrogen combustion. The presented extension of the ChemTraYzer approach can be used as a basis for methodological advancement of chemical kinetic modeling schemes and as a black-box approach to generate chemical kinetic models.

■

INTRODUCTION The development of reliable kinetic models for complex chemical processes can take several decades of combined experimental and theoretical research.1,2 In order to reduce this time span, various black-box methodologies have been proposed for accelerating the tedious kick-off phase of chemical kinetic modeling.3−5 One of the most widely used chemical kinetic modeling approaches6 is the combination of empirical low-cost/high-uncertainty methods (e.g., rate rules,7 analogies,8 or group additivity9) with ab initio high-cost/lowuncertainty methods (e.g., high-level quantum mechanics (QM)10). Although this combination avoids some of the disadvantages of the separate methods, it still fails at a central aspect of chemical kinetic modeling: making sure that all relevant reactions are included. In the present work, this shortcoming of established chemical kinetic modeling approaches is remediated by using reactive molecular dynamics (MD) simulations in combination with QM calculations to provide a sound basis for detailed chemical kinetic modeling. The long-standing history of detailed chemical kinetic modeling enabled researchers to develop intuition and rely on their experience. As stated by Levenspiel, “Intuition and past practices were the prime guides for the practitioner”.11 As a result of these “past practices”, combustion scientists started using analogies and developed reaction classes for approaching fuel candidates more systematically.12 The more systematic the procedures became, the simpler it was to let a computer do the tedious parts of the job.13 Nowadays, software for reaction mechanism generation, such as RMG,3 EXGAS,4 or RING,5 is capable of estimating chemical kinetic models for a large variety of fuel candidates.14,15 The established rule- and classbased approaches, however, lack the flexibility to provide © XXXX American Chemical Society

elementary reactions beyond the predefined domain of chemical space. In order to elude the rigid nature of established approaches, various chemistry exploration approaches were proposed in recent studies.16−23 In 2013, Zádor and Najm16 developed the KinBot software package for automatically exploring potential energy surfaces (PESs). The underlying algorithms used predefined patterns in molecular connectivities to discover possible reactions. Ab initio methods were used to obtain reliable temperature- and pressure-dependent rate constants from Rice−Ramsperger−Kassel−Marcus (RRKM)/master equation (ME) simulations. In 2014, Rappoport et al.17 reported a heuristic approach for discovering the chemistry of complex molecular structures. The reaction pathways of these structures were predicted using rule-based transformations, and the kinetics was predicted as flux through high-energy stable stationary points as pseudo transition states (TSs) connecting low-energy species. In 2015, Bergeler et al.,18 Suleimanov and Green,19 and Zimmerman20 proposed empirical and heuristic rule-based chemistry exploration schemes. These schemes were used to explore the chemistry of catalytic ammonia production,18 combustion and atmospheric chemistry,19 and organic chemistry and hydrogen storage.20 Similar approaches are being developed for automatic prediction of catalysis mechanisms using, e.g., a global reaction route mapping strategy24 or graph-based reaction path sampling.25 Parallel to the above rule-based exploration approaches, recent studies made use of the information gained from Received: February 9, 2018 Published: June 13, 2018 A

DOI: 10.1021/acs.jcim.8b00078 J. Chem. Inf. Model. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Information and Modeling reactive molecular dynamics (MD) simulations.21−23 In 2014, Wang et al.21 used Born−Oppenheimer MD (BOMD) at the Hartree−Fock level to study the formation of complex molecular structures from small molecules in a so-called ab initio nanoreactor, accelerated using a virtual piston periodically squeezing the simulation domain.26 This approach was used to simulate the Urey−Miller experiment, in which the formation of complex biomolecules from simple amino acids is studied.27 In 2014, Zheng and Pfaendtner22 proposed a reaction exploration scheme based on Car−Parrinello MD (CPMD) at the BLYP level of theory, accelerated via metadynamics. This scheme was used to study the oxidation chemistry of single methanol molecules. In 2015, Dö ntgen et al.23 developed a scheme to simultaneously obtain pathways and rate constants from reactive MD simulations. This so-called Chemical Trajectory Analyzer (ChemTraYzer) exploration approach was successfully applied to derive a reaction mechanism for hightemperature methane combustion. In their work, Döntgen et al.23 utilized simulations with the reactive force field ReaxFF28 for chemistry exploration, identifying reaction events by evaluating changes in the molecular connectivity, similar to the work of Wang et al.21 and Zheng and Pfaendtner.22 Compared to ab initio MD,21,22 however, 104 times longer time spans could be simulated using the ReaxFF force field, yielding significantly improved phase space sampling.23 The ReaxFF force field has been used to study a wide range of chemical processes, but one focus of present research is the use of ReaxFF in combustion sciences.23,29−34 Although current parameter sets for ReaxFF exhibit some shortcomings,23 the ability to actually explore novel chemistry was proven recently.35 Here the previous ChemTraYzer version will be extended toward chemical kinetic modeling by addressing the following points:

present methodology because of its limited size as well as an interesting target for in silico experiments using reactive MD. On the other hand, we study the CH4 + O2 system, which is substantially more complex and suffers from considerable inconsistency among different models.39 The presently generated H2 and CH4 oxidation models will be tested against literature thermochemistry, kinetics, and ignition delay times. The presented kinetic modeling approach is intended to produce initial chemical kinetic models on a level comparable to the established approaches, such as RMG,3 but with higher flexibility in terms of exploring the chemical space. To the best of our knowledge, this is the first time that MD and QM have been combined in an automated scheme to provide detailed chemical kinetic models.

■

METHODOLOGY A chemical kinetic model is composed of three parts: elementary reactions, thermodynamic data, and reaction rate coef f icients. For each of these three “pillars of detailed chemical kinetic modeling”,40 the models and methods used in the present work will be discussed and put into context in the following. Elementary Reactions. The first part of a chemical kinetic model is the set of elementary reactions. Since the current research in this field has been reviewed in the Introduction, this section will focus on how the present work extends the previous version of the ChemTraYzer software package. The first step toward chemical kinetic modeling is the recalculation of TSs and species via QM, for which the respective molecular geometries must be extracted from the MD simulations. Here an on-the-fly hybrid MD and QM simulation concept is developed for the following reason. Storing molecular geometries, i.e., all atomic positions, for a long MD simulation exceeds typical hard-disk capacities. In order to drastically reduce hard-disk requirements and access significantly more information, the ReaxFF trajectory simulations are cut into batches. Each of these batches is processed directly after it has been simulated and is discarded after processing, leading to significantly reduced hard-disk requirements. This kind of “batched” simulation allows the rate of writing to the hard disk to be increased (each 0.5 fs in the present work). This not only improves the precision of reaction event detection but further allows for extraction of molecular geometries directly from the trajectories. Once extracted from the trajectory, the molecular structures are preoptimized at the ReaxFF level of theory. For stable stationary points, molecular structures are extracted and processed periodically in order to screen the conformational space of the respective species (each 5 ps in the present work). Different conformers of a species are distinguished by comparing the energies and principal moments of inertia of the energy-minimized molecular structures at the ReaxFF level of theory. For unstable stationary points, i.e. TSs, preoptimization is done by freezing the atoms “actively” involved in the bond cleavage and formation and energy-minimizing the remaining molecular structure. Again, different conformers of a TS are distinguished by comparing the energies and principal moments of inertia of the optimized TS structures at the ReaxFF level of theory. For preoptimized TSs, however, the positions of the “active” atoms may differ for the exact same conformer. As a consequence, rather than being defined by a unique set of properties at the ReaxFF level of theory, a TS conformer is defined by ranges of energies and moments of

(i) So far, TS structures have been extracted and recalculated via QM manually.23 The extraction and recalculation via QM of molecular geometries of species and TSs will be automated. (ii) So far, the manually obtained QM results have not been processed further, and solely the QM energies have been used. The postprocessing of QM results will be conducted in an automated manner to provide thermochemistry and kinetics data. (iii) Thermochemistry and kinetics data will be provided in the ChemKin format36 and thus be ready-to-use for zero-dimensional numerical simulations of, e.g., shock tube experiments. For validation, we choose two fundamental systems of combustion kinetics. On the one hand, we study the H2 + O2 system, which still suffers from lack of consistent experimental characterization.37 In 2008, Konnov stated that “in spite of apparent simplicity, contemporary reaction mechanisms for the H2−O2 system developed by different authors often differ by the number of reactions and their rate constants” and that further analysis and more accurate experiments would be necessary.37 Accordingly, Burke et al.38 found that “present uncertainties in the temperature and pressure dependence of rate constants for HO2 formation and consumption reactions ... substantially affect predictive capabilities”. Therefore, this system provides a good starting point for validating the B


Article

Journal of Chemical Information and Modeling inertia (cf. the Supporting Information for details on how these ranges are obtained). Only if the energy and moments of inertia of a molecular structure are included in the property ranges of a TS conformer are the two considered to be equal. The species’ and TSs’ molecular structures are stored in a database for later comparison with other molecular structures and to be used as input for QM calculations. These QM calculations produce results used for predicting thermodynamic and kinetics data as described in the following sections. Thermodynamic Data. Since the work of Benson and Buss,41 the systematic prediction of thermodynamic data has developed strongly42 and is on the verge of reaching experimental accuracy today,43 at the cost of computational power. While the highest levels of accuracy can be achieved for single compounds, the prediction of thermodynamic data for all species of, e.g., a detailed chemical kinetic model requires computationally less costly levels of theory, at the cost of higher uncertainties. The total uncertainty of thermodynamic data prediction originates from the uncertainties stemming from the electronic structure method used to calculate the PES of a species’ internal degrees of freedom (DoFs) and from the description of the motions of the respective atomic nuclei. While the lowest uncertainties would be achieved by simultaneously solving the electronic and nuclear parts of the Schrödinger equation, or at least by including the coupling of internal DoFs in the nuclear part of the Schrödinger equation, these approaches are hardly feasible44,45 and have minor effects on heavy atoms46 or are limited to small molecules,47,48 respectively. Currently, the internal DoFs are typically treated as decoupled,49,50 although this assumption is known to fail in some cases, making manual intervention inevitable.51 When the level of accuracy is reduced even further, we gain a more robust, automatable approach, which is described in the following. In the present work, thermodynamic data are predicted using the rigid rotor−harmonic oscillator (RRHO) model with approximate one-dimensional hindered rotor (1DHR) corrections. Although this is not the most accurate approach, it is the best trade-off among computational cost, accuracy, and automatability to our understanding. Truhlar and co-workers49,52−54 reviewed and developed numerous approximations for the anharmonic partition function of bond torsional modes. Here the anharmonicity of bond torsional modes of stable and unstable stationary structures is estimated as a smooth transition from the harmonic oscillator (HO) to the free rotor (FR) partition function.52−54 This transition is described as the hyperbolic tangent of the ratio of the free rotor and classical (clas) harmonic oscillator partition functions. Chuang and Truhlar52−54 reported that the maximum deviation from the true true 1DHR partition function amounts to a factor of 2.5 at high temperatures for C2 species. Here this transition model is applied to the sum of N conformers and the product of M bond torsional modes and yields an anharmonic correction Qcorr anh to the harmonic oscillator partition function:

corr Q anh =

Q anh

) (

E

)

i ∑i = 1 ∏j = 1 Q HO, exp − k Ti j B ÄÅ ÅÅ M GM ÅÅ ∏j = 1 Q FR, ÅÅ j Å ≈ tanhÅÅ E ÅÅ N M i ÅÅ ∑i = 1 ∏j = 1 Q clas, j exp − k Ti ÅÅÇ B N

(

M

(

) (

)

ÑÉÑ ÑÑ ÑÑ ÑÑ ÑÑ ÑÑ ÑÑ ÑÑÖ

(1)

where Qanh is the total anharmonic partition function, QiHO,j is the harmonic partition function for the ith conformer and jth mode, Qiclas,j is the respective classical partition function, QGM FR,j is the free rotor partition function for the jth mode of the global minimum (GM), kB is the Boltzmann constant, and T is the temperature. In this work, the molecular geometry of the GM is used to determine the moments of inertia of rotational modes. While this approximate 1DHR treatment cannot account for coupling of internal DoFs, its low- and high-temperature limits52−54 and functional form55 make it physically reasonable. Reaction Rate Coefficients. While the representation of rate constants has only slightly changed since the early work of Arrhenius,56 the prediction of rate constants has developed strongly.43,57 State-of-the-art kinetics predictions are based on Rice−Ramsperger−Kassel−Marcus (RRKM) theory and master equation (ME) simulations,58−60 which depend not only on the species’ and TSs’ molecular properties but also on the collisional and energy transfer dynamics.58 Modeling and quantification of these dynamic processes is an issue in current research,61,62 the methods of which are computationally too costly to be applied to all species of, e.g., a detailed chemical kinetic model. Therefore, we consider the canonical TST formulation as the most robust kinetic theory in terms of automatability today and use it in the present work for predicting high-pressure rate constants:57 k(T ) =

⧧ kBT ji ΔG (T ) zyz zz expjjj− j h kBT z{ k

(2)

where h is Planck’s constant and ΔG⧧(T) is the free energy barrier. For barrierless reactions, no TS can be found in potential energy, and performing QM calculations for such a reaction would be a waste of computational resources. In order to avoid this waste of computational power, rule-based detection of barrierless reactions is conducted whenever a reaction is detected in the MD simulations. This rule-based detection is based on the assumption that any bond formation or bond cleavage from/to unpaired electrons is barrierless. From this, predefined patterns of reactants’ and products’ spin multiplicities can be defined. If the spin multiplicities of a reaction’s reactants and products match one of these patterns, the reaction is identified as barrierless. These patterns include the formation of a singlet from two doublets, the formation of a doublet from a doublet and a triplet, and the formation of a singlet from two doublets and a triplet. Reactions are checked against these patterns in the forward and backward directions. Once identified, the rate constants of the barrierless reactions are calculated from the MD simulations as described by Kröger et al.63 Simulation Details. The chemical composition of the fuel/oxidizer mixture is defined by the desired number of molecules of each compound. For each species, the database is checked for an initial molecular structure. If no structure is C


Article

Journal of Chemical Information and Modeling

Figure 1. Program structure illustrated as a flowchart. Each colored box represents one object/class in the program structure.

conducted at constant temperature, controlled via a Nosé− Hoover thermostat68,69 with a damping constant of 10 fs. In order to reduce the amount of QM calculations, preliminary simulations are conducted to produce QM properties of H2, CH4, O2, and prominent species showing up during the inception stage of hydrogen and methane high-temperature oxidation. The molecular geometries at the ReaxFF level of theory are used as initial geometries for QM optimization of stable and unstable stationary points. Here the G3MP2 level of theory70 is employed to obtain molecular geometries, harmonic frequencies, hindered rotor information, and single-point energies using the Gaussian software package.71 TSs are validated using intrinsic reaction coordinate (IRC) scans at the MP2/631G(d) level of theory, which is also used in the G3MP2 compound method. Although it is not the most recent method, the G3MP2 level of theory is based on MP2 geometry optimization, which allows optimization of several HXOY TSs that cannot be found at the B3LYP level of theory. The level of theory applied for QM, however, can be changed easily without constraining the above-described procedure. It should be noted, however, that “The One” QM method with universal accuracy does not yet exist. Eventually one will encounter complex electronic structures with, e.g., multi-

available, its molecular geometry is estimated using the Open Babel software package.64 The initial configuration is generated using the Packmol software package65 so that the desired density in a cubic simulation box is matched. Initial velocities are distributed randomly so that the Maxwell−Boltzmann distribution at the desired temperature is resembled. In the present work, the inception stage of hydrogen and methane oxidation is studied using H2 + O2, CH4 + H2 + O2, and CH4 + O2 trajectories at a density of 1000 mol/m3 and a temperature of 3000 K, corresponding to an ideal gas pressure of roughly 246 atm. Reactive molecular dynamics simulations with the ReaxFF28 force field parameters for hydrogen combustion of Agrawalla et al.66 are performed using the LAMMPS software package.67 The chemical composition for each trajectory is stoichiometric, with 10 H2 molecules in the hydrogen oxidation trajectories, a single additional CH4 molecule in the CH4 + H2 oxidation simulations, and 10 CH4 molecules in the methane oxidation trajectories. Although the force field of Agrawalla et al.66 is trained for H2 + O2, simulations with methane molecules are conducted in order to broaden the database for comparison with the literature. It should be noted that for pure methane oxidation simulations, the recent CHO force field of Ashraf and van Duin33 would be better suited (cf. the Supporting Information). Simulations are D


Article

Journal of Chemical Information and Modeling reference effects, which are clearly beyond the capability of simple black-box QM methods. As a consequence, refining a fraction of the present QM results manually would be inevitable to improve the initial chemical kinetic model. Program Structure. From a technical point of view, the presented methodology wraps around the previous ChemTraYzer code and implements a framework for input, intermediate, and output data. In order to illustrate the functionalities and how they are interconnected, the process flow is compiled in Figure 1. The top part of the flowchart represents the simulation initialization. After initialization, the simulation starts and the aforementioned batched simulation/processing scheme is handled by the processes in the main body of the flowchart. Therein, the simulation and processing steps are included in the “Simulation” and “Processing” classes, respectively. At the end of the batched simulation/processing run, the QM data are received from the Gaussian output files. This information is stored in the database and used to predict thermochemistry and kinetics data, which are exported in the ChemKin format. The present extension of the ChemTraYzer software package is available online at www.sourceforge.net/projects/ chemtrayzer/.

Figure 2. Correlation between equilibrium constants obtained from the high-temperature hydrogen and methane oxidation trajectories and the GRI-Mech 3.0 model. The red dashed line represents identity of the equilibrium constants. The value of the correlation coefficient rXY is also shown.

■

RESULTS AND DISCUSSION Validation of ReaxFF Equilibrium Constants. The ReaxFF rate constants and the QM properties obtained from the above-described high-temperature hydrogen and methane oxidation simulations are used to compute equilibrium constants Keq for reactions observed during the simulations. Literature reference values are obtained from GRI-Mech 3.0,72 denoted as GRI3 below, and from Aramco Mech 1.3,39 denoted as ARAMCO below. These are state-of-the-art models for describing gas-phase chemical processes involving methane. The major difference between these models is the inclusion of larger hydrocarbons (which are not required here) and formic acid-based chemistry in the ARAMCO model. Previously, formic acid and its radicals were found to be important during high-temperature methane combustion simulations via ReaxFF.23 Thus, the ARAMCO model allows for validation of formic acid-based chemistry. Both reference models are largely based on group additivity thermochemistry data, which is typically accurate for the presently studied species in the temperature range for which the group additivity values have been proposed. ReaxFF equilibrium constants Keq,ReaxFF are computed as the ratios of the forward and backward rate constants at 3000 K. In Figure 2, the correlation between the ReaxFF equilibrium constants and the equilibrium constants computed via the GRI3 thermochemistry is shown. The statistical uncertainty of the ReaxFF equilibrium constants results from the uncertainties of the rate constants used to calculate the equilibrium constants. The worst case is the prediction of rate constants for reactions with only one observed reaction event, for which the statistical uncertainties amount to roughly 1 order of magnitude for a 95% confidence level.63 This rate constant uncertainty would cause a deviation of about 2 orders of magnitude for the respective equilibrium constant. The correlation between the ReaxFF and literature equilibrium constants is close to unity, indicating reliable reproduction of the literature thermochemistry in the ReaxFF simulations. When interpreting these results, however, one should keep in mind that the force field used for the present

reactive MD simulations was trained for the presently studied system. Since for most reactions the backward pathway could not be observed during simulation, the data set for comparing ReaxFF with the literature is too small to make a solid statement on the overall performance of the presently used force field. Validation of QM Equilibrium Constants. QM equilibrium constants Keq,QM can be obtained for any temperature and any reaction if all of the reactants and products are available in the present database. Figure 3 shows the correlations between the QM equilibrium constants and the equilibrium constants computed via (a, b) the GRI3 thermochemistry and (c, d) the ARAMCO thermochemistry for 1000 and 2000 K. In these correlations, any reaction is considered if its reactants and products can be found in the present database and in the GRI3 or ARAMCO model. The correlations between the present QM and literature equilibrium constants are very high, irrespective of the reference data (either the GRI3 or ARAMCO model). Especially at 1000 K, the correlation coefficients are almost unity, indicating very reliable reproduction of the reaction free energies reported in the literature. With increasing temperature, however, the correlation coefficients slightly decrease for both reference cases, and the decrease becomes more pronounced as the temperature increases further (cf. Figure S1). As will be discussed later (cf. Table 2), this trend is due to the combination of an approximate 1DHR model and the lack of a model to describe bond stretching and angle bending anharmonicity at high temperatures. Chuang and Truhlar52−54 indeed reported that the uncertainty of the approximate 1DHR model increases with increasing temperature. It has to be noted that the correlation between the GRI3 and ARAMCO models decreases with increasing temperature. In fact, the present ab initio model correlates better with each model than the two models correlate with each other above 1000 K (cf. Figure S1). Thermochemistry Uncertainties. Uncertainty quantification plays a major role in chemical kinetic modeling, as E


Article


Figure 3. Correlations between equilibrium constants obtained via quantum mechanics and those from (a, b) GRI-Mech 3.0 (GRI3) and (c, d) Aramco Mech 1.3 (ARAMCO). The red dashed lines represent identity of the equilibrium constants. The values of the correlation coefficient rXY are also shown.

Table 1. Standard Enthalpies of Formation and Their Deviations from Tabulated Literature Values (ATcT)74 for Three Oxygen-Related Species at the G3MP2 and CCSD(T)/aug-cc-pVXZ Levels of Theory in kJ/mola species O2 error HO2 error H2O2 error

exptl74 0.00 12.26 −135.46

G3MP2

X=T

X=Q

X=5

X=6

9.49 9.49 16.51 4.25 −127.91 7.55

23.06 23.06 38.97 26.71 −104.12 31.34

10.01 10.01 23.18 10.92 −122.88 12.57

6.46 6.46 18.92 6.66 −127.83 7.62

4.79 4.79 17.06 4.80 −129.90 5.56

a

B2PLYPD3 and aug-cc-PVTZ energy-minimized geometries and harmonic frequencies were used as input for CCSD(T) calculations and to obtain the thermal contributions, respectively.

in future studies. Proppe et al.73 moreover discussed the uncertainties stemming from anharmonic contributions to the partition function. Besides using an approximate 1DHR model and in addition to the points raised by Proppe et al.,73 the present QM-based thermochemistry prediction is subject to three main sources of error:

uncertainties are essential for merging data from different sources in the same chemical kinetic model. Proppe et al.73 applied a fundamental approach for quantifying uncertainties in the context of complex chemical networks. In there, the authors used Bayesian statistics to assess the uncertainty of the electronic structure method used.73 While this approach is beyond the scope of the present work, it should be considered F


Article


and 1500 K while having a negligibly small effect on the values at around 3000 K (cf. Table 2: Δcp,B3LYP/R). The combination of errors introduced by anharmonicity and the level of theory covers most parts of erroneous heat capacity prediction at the G3MP2 level of theory. In summary, the equilibrium constants obtained via the presented hydrogen trajectory and quantum chemistry approach are reliable. The increasing uncertainty with increasing temperature is related first to the approximate 1DHR model and second to modeling of bond stretching and angle bending modes as harmonic oscillators, which fails to describe the strong anharmonicity of these modes at high temperatures. Validation of Rate Constants. The rate constants obtained from the hybrid trajectory and quantum chemistry simulations are partly based on TST and partly based on analysis of the trajectories. The trajectory-based rate constants are used for barrierless reactions exclusively, since TST is not applicable to that type of reaction. Figure 4 shows the correlations between the presently predicted rate constants and the respective rate constants obtained from (a) the GRI3 model and (b) the ARAMCO model. In this figure, the rate constants at 1000 and 2000 K are presented in a single correlation plot for each reference, and a single correlation coefficient is given for both temperatures. Although rate constants for barrierless reactions are available only for the simulation temperature of 3000 K, these rate constants are assumed to be temperature-independent for simplicity. For both reference models, the presently predicted ab initio rate constants are in good correlation with the literature rate constants. Similar to the equilibrium constants, the rate constants correlate better at low temperature than at high temperature (cf. Figure S2). For the rate constants, however, this difference is much more pronounced: low-temperature rate constants have a correlation coefficient of roughly 90% with the literature values, whereas high-temperature rate constants partly have a correlation coefficient of roughly 55% (cf. Figure S2 for details). Again, this is likely due to the approximate hindered rotor modeling and the missing model for bond stretching and angle bending anharmonicity. In fact, the description of torsional anharmonicity is even more approximative for TSs, as their hindered rotor description is estimated from the reactants’ hindered rotors. Comparing the ReaxFF-based rate constants for barrierless reactions (circles) to those calculated via QM-based TST (diamonds) shows similar scattering for the two methods. This leads to the conclusion that barrierless reactions described via ReaxFF are about as reliable as reactions with barriers described via G3MP2 (cf. Figure 7). The reasonable performance of ReaxFF for barrierless reactions can be explained by the fact that barrierless reactions are mostly dominated by the relative translational motion of the reactants, which is described well by ReaxFF. In summary, the rate constants obtained via the presented hydrogen trajectory and quantum chemistry approach are still reliable in comparison to the GRI3 and ARAMCO kinetics. The low correlation at high temperatures should be addressed in future studies and is identified as a major drawback of the present kinetics predictions. With the thermochemistry and kinetics data combined, we obtain an initial chemical kinetic model that is ready to use for numerical simulations. Now, performing sensitivity analyses based on this initial chemical kinetic model would reveal which reactions require further

(1) Standard enthalpies of formation from atomization energies for species with oxygen−oxygen bonds suffer from slow basis set convergence of the respective single-point energy calculations. Table 1 shows the enthalpies of formation for O2, HO2, and H2O2 from the literature and at the G3MP2 and CCSD(T)/aug-cc-pVXZ levels of theory. In the latter, the standard enthalpy of formation for compound i is calculated according to ΔfHi° = (Hi,QM ° − Hatoms,QM ° ) + ΔfHatoms ° , in which “atoms” refers to the atomization of compound i. The elements’ standard enthalpies of formation are taken from the ATcT databse.74 The G3MP2 method is designed to extrapolate to the complete basis set limit, which appears to perform differently for the different oxygen-containing species. While for HO2 and H2O2 G3MP2 yields results similar to those of CCSD(T)/aug-cc-pVXZ with X = 6 and X = 5, respectively, the G3MP2 result for O2 is closest to that of CCSD(T)/aug-cc-pVQZ. The largest impact on thermochemistry prediction can be observed for O2, for which the uncertainty is on the same order as the change in enthalpy over temperature. In the present work, the enthalpies of formation of small oxygen−oxygen bond containing species, i.e., O2, HO2, and H2O2, are taken from tabulated literature values instead of being calculated from atomization energies. This avoids unphysical thermochemistry stemming from slow basis set convergence of the methods used for these species. (2) Anharmonic bond stretching and angle bending modes are approximated using the harmonic oscillator model. This harmonic approximation is well-known to cause errors in heat capacity prediction, which was studied in detail, e.g., for methane by Nguyen and Barker.75 The authors found that although at 500 K the anharmonicity is negligibly small, at 3000 K the anharmonic contribution Δcp,anh amounts to R, the universal gas constant (cf. Table 2: Δcp,anh/R). Thus, inclusion Table 2. Methane Heat Capacity Errors Due to Anharmonicity and Level of Theory T/K

Δcp,anh/Ra

Δcp,B3LYP/Rb

Δcp,G3MP2,res/Rc

500 1000 1500 2000 2500 3000

0.06 0.28 0.50 0.69 0.86 1.00

0.31 0.36 0.27 0.19 0.14 0.09

0.05 0.09 0.09 0.13 0.19 0.28

a

Error resulting from comparison of harmonic and anharmonic bond stretching and angle bending (anh). bError resulting from comparison of of MP2 to B3LYP harmonic frequencies (B3LYP). cRemaining error of the G3MP2 method (res), given by Δcp,G3MP2,res = cp,G3MP2 − cp,lit − cp,anh − cp,B3LYP.

of anharmonicity of bond stretching and angle bending would yield larger heat capacities compared with using the pure harmonic oscillator model. Two potential approaches for including anharmonic bond stretching would be (i) to calculate all of the bond dissociation energies (BDEs) or (ii) to carry out a second-order perturbation calculation.76 (3) The level of theory applied to calculate harmonic frequencies causes an error in partition function prediction. For comparison, methane harmonic frequencies are calculated at the B3LYP/TZVP level of theory, which was shown to be slighly more accurate than MP2.77 In contrast to anharmonicity, improving the harmonic frequency calculation mostly affects the heat capacity values at temperatures between 500 G


Article


Figure 4. Correlations between rate constants obtained from the high-temperature methane oxidation trajectories and from (a) GRI-Mech 3.0 (GRI3) and (b) Aramco Mech 1.3 (ARAMCO). The red dashed lines represent identity of the rate constants. QM/TST-based rate constants are marked as diamonds, while ReaxFF-based rate constants for barrierless reactions are marked as circles. The values of the correlation coefficient rXY are also shown.

Moreover, the former model builds on the latter and thus inherits any uncertainty embodied in it. Furthermore, the high temperatures of the methane oxidation simulations mostly lead to CH4 consumption via dissociation to Ċ H3 and Ḣ rather than hydrogen abstraction. At lower temperatures, however, hydrogen abstraction pathways dominate ignition and have to be included in the detailed chemical kinetic model. Therefore, we added hydrogen abstraction via Ö , Ȯ H, and O2 from the GRI3 model to the present ab initio model, and hereafter we refer to this as the modified ab initio model. For future reactive MD simulations, the initial chemical composition should include a variety of radicals to enable hydrogen abstraction prior to fuel dissociation. Ignition delay times calculated with the modified ab initio model are compared to those for methane combustion obtained using the GRI3 model72 and an RMG-generated model.3 To allow for fair comparison, the pressure-dependent reactions in the GRI3 model are fixed at the high-pressure limit and RMG is constrained to use high-pressure-limit rate constants. Figure 6 shows the ignition delay times calculated using the three models based on the aforementioned 95% fuel consumption criterion using the Cantera software package.79 As expected, the uncertainty of the predicted ignition delay times with respect to the literature is much higher for methane oxidation than for hydrogen oxidation. In contrast to hydrogen oxidation, the predictions appear to worsen with increasing pressure. While this is partly due to uncertainties attached to the modified ab initio model, the GRI3 model itself predicts too-fast ignition at high pressures.82 The methane oxidation model provided by the well-established RMG approach differs from the GRI3 model in the activation energy and absolute value of τig. Comparing the RMG model to the present ab initio model is more just, since both models are based on automated black-box reaction mechanism generation approaches. From a mechanistic point of view, the present ab initio model includes reactions involving the Criegee intermediate (H2Ċ OȮ ), ozone (O3), and formic acid-based chemistry.23,35

attention, and their rate constants could be recalculated at higher levels of theory. Ignition Delay Times. The ignition delay time (τig) calculations with the present ab initio model are compared to hydrogen combustion as modeled by Ó Connaire et al.78 Figure 5 shows the ignition delay times calculated with the chemical kinetic model of Ó Connaire et al.,78 denoted as Connaire below, and with the present ab initio model. To allow for fair comparison, the pressure-dependent reactions in the Connaire model are fixed at the high-pressure limit. Ignition delay times are obtained as the time of 95% fuel consumption from numerical simulations conducted with the Cantera software package.79 For most temperatures and pressures shown in Figure 5, the ignition delay times are overestimated by the present model compared with those predicted with the Connaire model. As can be seen from the deviation plot (cf. Figure 5b), the present predictions overestimate the ignition delay times by a factor of 6 at most. This deviation decreases with increasing pressure and approximately amounts to a factor of 2 over the whole temperature range for p = 1000 atm. Under high-temperature and low-pressure conditions, the ignition delay times are underpredicted by a factor of 8 by the present model compared with the Connaire model. The deviations between the present model and the Connaire model appear to arise mostly from deviations in the underlying rate constants, namely, inaccuracies in describing O + OH. As discussed in the next section, ignition delay times have the highest sensitivity to the reaction of these reactants irrespective of the chemical kinetic model used for the simulations. Since the H2 mechanism is typically used as a basis for larger mechanisms of hydrocarbon fuels, the uncertainties reported here would add up to further uncertainties in the fuel-specific mechanisms. However, the H2 base mechanism is well-studied80,81 and can be taken from the literature rather from the presently proposed hybrid trajectory and quantum chemistry approach. For methane oxidation, the underlying reaction mechanism is significantly more complex than that for hydrogen oxidation. H


Article


Figure 5. (a) Comparison of H2 ignition delay times calculated with the Connaire model78 (black) and the present ab initio model (red) for different constant pressures (line styles apply to both colors) and (b) the deviations between the predictions of the two models. The fuel (xfuel = 10%) was mixed with a 21% O2/79% Ar mixture to resemble Φ = 0.26. Pressures range from 1 to 1000 atm.

Figure 6. (a) Comparison of CH4 ignition delay times calculated with the GRI3 model72 (black), an RMG-generated model3 (blue), and the modified ab initio model (red) for different constant pressures (line styles apply to all colors) and (b) the deviations between the predictions of the GRI3 model and the modified ab initio model. The fuel (xfuel = 10%) was mixed with a 21% O2/79% Ar mixture to resemble Φ = 1.06. Pressures range from 1 to 100 atm.

These pathways are not included in the GRI3 model and mostly kick in at high pressures, at which recombination reactions are entropically favored. Compared with the work of Döntgen et al.23 on methane oxidation, the only notable difference between the methane oxidation mechanisms is the presence of ozone chemistry in the present work. Here O3 mainly forms via the addition of O2 and Ö radicals. The Ö radicals in turn are mostly produced via H2 + O2 chemistry at high temperatures. The formation of an ozone-related compound is even more pronounced under high-pressure conditions: O2 and Ȯ H radicals add to O3H. From a sensitivity point of view, however, these reactions are 2 orders of magnitude less significant compared with the reaction of Ö and Ȯ H to form HOȮ , to which ignition delay times are most sensitive (cf. the Supporting Information).

H2 + O2 Chemical Kinetic Model Rate Constants. In Figure 7, the rate constants of reactions involved in the above hydrogen combustion simulations are compared to the rate constants used in the Connaire model.78 To judge the deviations between the present ab initio rate constants and those used in the Connaire model, the ratios of these rate constants are plotted. The present ab initio rate constants deviate by up to roughly 2 orders of magnitude over the temperature regime from 1000 to 3000 K at most. Most of the rate constants, however, agree with the literature rate constants within 1 order of magnitude (cf. the shaded area in Figure 7). Interestingly, most of the ReaxFF-based rate constants used to estimate the kinetics of barrierless reactions are similarly accurate as the G3MP2-based rate constants for reactions with barriers (as mentioned above already). For each approach I


Article


chemistry simulations to obtain initial chemical kinetic models. This methodology is a more flexible alternative to the established rule- and class-based approaches for accelerating the tedious kick-off phase of detailed chemical kinetic modeling. An ab initio model for high-temperature hydrogen and methane oxidation was obtained from hybrid trajectory and quantum chemistry simulations in a fully automated manner. In contrast to the pure ReaxFF thermochemistry, the quantumchemical updates allow for the prediction of reliable equilibrium constants. Kinetics predictions from the hybrid simulations were found to be reasonably reliable, and minor improvements would be necessary in order to reproduce hydrogen ignition delay times from the literature. In the current state, the presently predicted hydrogen ignition delay times differ from literature predictions by a factor of 8 at most (while considering high-pressure-limit rate constants exclusively). For methane oxidation, the present ab initio model would provide a sound basis for further detailed chemical kinetic modeling, comparable to chemical kinetic models obtained via, e.g., RMG. The QM extension of the ChemTraYzer software package23 can be used to obtain purely ab initio models without any prior knowledge about the studied compounds. These initial chemical kinetic models can directly be used to numerically simulate gas-phase chemical processes. In the current state, however, simulations are constrained to the high-pressure limit because of the use of TST for kinetics predictions. For finite pressures, the balancing of reactions must be done microscopically, thus requiring prediction of microcanonical properties.55 These properties are readily available from the presented QM calculations, and the inclusion of pressure dependence will be the subject of future methodological advances. Furthermore, future work will focus on quantifying the uncertainties of thermochemistry and kinetics predictions based on the presented methodology. The presented approach is intended to be for the methodologist, who is interested in exploiting the potential of MD and QM, and for the practitioner, who can use it as a black-box methodology for chemical kinetic modeling. While the present methodology is by no means a replacement for detailed chemical kinetic modeling studies, it provides a sound basis on which detailed models can be built.

Figure 7. Temperature-dependent ratios of the present rate constants from the H2 + O2 ab initio model to those taken from the Connaire model.78 The red dashed line corresponds to identity of the rate constants, and the shaded area surrounding it represents any ratio within 1 order of magnitude deviation.

(ReaxFF-based for barrierless reactions and G3MP2-based for reactions with barriers), the rate constant for one reaction is highly inaccurate compared with the rate constant used in the Connaire model (HO2 + H → H2O2 and O2 + H2 → HO2 + H, respectively). While for the ReaxFF-based rate constants the uncertainty of the underlying force field could explain this deviation, the G3MP2-based rate constants deviate because of an erroneous spin. In the present work, spin multiplicities are derived from the number of electrons using the lowest spin multiplicity possible. However, the lowest spin multiplicity does not always correspond to the lowest-energy configuration, and there is a list of well-known exceptions, including molecular oxygen (triplet rather than singlet). Fast detection of the correct spin multiplicity will be addressed in a future study by using simple on-the-fly Hartree−Fock calculations of the lowest and next-highest spin multiplicities. The two above-discussed strongly deviating reactions, however, have no significant impact on hydrogen ignition delay time predictions because of their low sensitivities (

Automated Chemical Kinetic Modeling via Hybrid Reactive Molecular

Recommend Documents