Jiang
Wang
*^{a},
Zhiling
Li
^{a} and
Wenli
Zhang
*^{b}
^{a}College of Science, Guizhou Institute of Technology, Boshi Road, Dangwu Town, Gui’an New District, Guizhou 550025, China. E-mail: cwangjiang@git.edu.cn; Tel: +86-0851-86871682
^{b}School of Transportation Engineering, Guizhou Institute of Technology, Boshi Road, Dangwu Town, Gui’an New District, Guizhou 550025, China. E-mail: wenlizhang@git.edu.cn

Received
4th June 2024
, Accepted 16th August 2024

First published on 16th August 2024

The dissociation of hydrocarbon bonds plays a pivotal role in their utilization, whether through fuel combustion or the thermo-cracking of large hydrocarbons in petroleum refinement. Previous studies have primarily focused on the effects of temperature, pressure, and chemical environment on hydrocarbon reactions. However, the influence of molecular configuration on bond breaking rates has not been thoroughly explored. In this study, we propose an approach to compute bond dissociation rates, and apply it to the reactive molecular dynamics simulation (ReaxFF) trajectories of three molecules: n-tridecane, n-pentane, and 1,3-propanediol. Our results reveal that the bond dissociation rate depends not only on the bond position in the chain, but also on the molecular configuration. Stretched configurations exhibit higher dissociation rates, particularly favoring the breaking of central bonds. Conversely, when the molecule is coiled, resulting in a reduced size, terminal bonds exhibit higher dissociation rates. This research contributes to a deeper understanding of molecular dissociation properties in the oxidation of hydrocarbons, and provides a way to explore the bond breaking properties of other molecules.

Taking n-alkane as a representative example, numerous theoretical studies have investigated their bond dissociation energy (BDE) using first-principle approaches or density functional theory (DFT).^{8,9} A consistent finding across these investigations is the non-uniform distribution of dissociation energy among C–C bonds at different positions. Specifically, the terminal C–C bond (α) exhibits the highest BDE, whereas the second C–C bond (β) demonstrates the lowest BDE. Remarkably, all other C–C bonds exhibit relatively similar BDE values. This variation in bond strengths sheds light on the distinctive reactivity of different carbon–carbon bonds within the alkane structure.

The characterization of bond fission properties is effectively captured through the dissociation rate constant.^{10–12} This constant, essential for understanding bond dynamics, is not solely contingent on the bond dissociation energy (BDE) or enthalpy. Instead, it is intricately influenced by the complex interactions between molecules. Moreover, the specific configuration of a molecule, as well as the location of the bond within the molecular structure, introduces additional variability into the bond-breaking rate constant. As shown by a number of experimental studies, the shape and conformation of molecules will affect their reaction rate.^{13–16} Therefore, the dissociation rate constant is dependent on both the molecular configuration and the site of the bond within the molecule.

While direct observation of chemical bond breaking is feasible for specific reactions,^{17–19} quantitatively exploring bond-breaking properties at the atomistic level in high-temperature processes such as combustion or pyrolysis remains challenging through experimental means. In this context, computer simulations offer an alternative avenue for studying bond dissociation behaviors. Traditional molecular dynamics (MD) simulations, however, fall short in capturing the formation or breaking of chemical bonds.^{20} On the other hand, quantum mechanical approaches (QM) can determine the properties of electrons in molecules and is widely used to simulate chemical reactions^{21} of small systems with high accuracy. In addition to QM approaches, various reactive force field MDs that could simulate chemical reactions have been proposed within the past decade.^{22–26}

To determine the bond breaking rate of a molecule in a specific configuration, a straightforward computational approach involves initially fixing the configuration and then simulating the reaction using either quantum mechanics (QM) or reactive molecular dynamics (MD) methods, and then measuring the breaking rate directly. However, the number of possible configurations of a molecule is large, especially when the molecule is long and has high degrees of freedom. This complexity makes it inaccessible to evaluate the bond dissociation rate for all possible configurations comprehensively. To circumvent this issue, we propose an alternative method. Instead of analyzing each configuration one after another using ‘bottom-up’ theories and methods, we conduct a large number of independent QM or reactive MD simulations. Each simulation samples the dissociation event of one bond in a possible configuration. By analyzing the ensemble of these simulation trajectories in a ‘top-down’ scheme, we can extract the bond breaking rate for specific configurations.

To implement the proposed approach, it is necessary to obtain a large number of samples of various molecular configurations at which the molecule dissociates. In this research, we will use MD based on the ReaxFF force field (FF)^{27–30} to generate this data, balancing simulation speed and accuracy.

Introduced by Van Duin et al., the ReaxFF force field (FF) is one of the most popular FF that is deliberately designed and parametrized by QM calculations, it is capable of accurately modeling bond formation and breaking in MD simulations, and maintains the computational efficiency close to classical MD,^{27–30} and it is widely adopted in research on chemical reactions to extract the most prominent information of larger systems. ReaxFF has proven particularly effective in studying phenomena such as fuel combustion,^{30–35} pyrolysis of n-alkanes,^{36,37} and other hydrocarbons.^{38–40}

Han et al. conducted a comprehensive investigation into the combustion of n-heptane using ReaxFF MD. Their findings revealed distinct properties among C–C bonds along the chain, highlighting a tendency for the central C–C bonds to dissociate prior to the terminal ones.^{41} This observation aligns with similar results reported by Ding et al.^{31} Moreover, studies focusing on the combustion or pyrolysis of hydrocarbons in diesel fuel consistently indicate that the dissociation properties of chemical bonds are contingent on their specific locations within the molecular structure.^{32,33,42}

Several methods are involved in the calculation of bond-breaking rate constants. The first method utilizes the Arrhenius equation, which provides a crucial link between the chemical reaction rate constant and the activation energy.^{43} In ReaxFF simulations, researchers often measure the overall reaction rate directly at various temperatures. Subsequently, activation energies can be deduced using the Arrhenius equation.^{28,31,35,44} However, the reverse process is not straightforward. Without prior knowledge of the activation energies for bond dissociation, calculating the bond-breaking rate constant using the Arrhenius equation becomes impractical, particularly for bonds in molecules with different configurations. Alternatively, variational transition state theory (VTST)^{45–47} and the POLYRATE program developed by Truhlar et al.^{48,49} are widely used to calculate and predict reaction rates for clearly defined reaction in specific configurations.^{50–53} However, it has not been directly applied to the calculation of single bond-breaking rates for molecules over different configurations. To be the best of our knowledge, to date, there is currently no established method for quantitatively calculating bond dissociation rate constants as a function of molecular configurations.

In this research, we present an approach to calculate bond dissociation rate constants in the oxidation and pyrolysis reactions of hydrocarbons from a ‘top-down’ scheme, with a specific emphasis on linear chain polymers. By applying this approach to ReaxFF simulation data, our study reveals that bond-breaking rates depend not only on the location of the bond within the molecule but also on the overall molecular configuration.

The findings of this research contribute to a deeper understanding of the intricate bond-breaking properties exhibited by hydrocarbons during the oxidation process. The proposed method is not limited to the studied hydrocarbons and ReaxFF simulations but is also extendable to different molecules in other types of reactive simulations or quantum mechanical (QM) simulations.

A molecule has a dissociation rate constant λ, which can be defined as:

(1) |

The dissociation rate constant, denoted as λ for a molecule, leads to the probability density of its lifetime t to obey an exponential distribution:

ρ(t,λ) = λe^{−λt}, t ∈ [0, +∞) | (2) |

The average lifetime can be computed as follows:

(3) |

(4) |

However, it is crucial to note that the dissociation rate calculated through this method represents the overall rate, encompassing contributions from all bonds under all possible molecular configurations.

The molecular configuration r follows a Boltzmann distribution when the system is in equilibrium:

(5) |

As λ(r) represents the distribution of the rate constant as a function of molecular configuration, a crucial relationship can be deduced:

(6) |

(7) |

Upon molecular dissociation, the configuration probability density distribution is denoted as ρ_{b}(r), where subscript ‘b’ signifies ‘breaking’, with . It is crucial to distinguish ρ_{b}(r) from ρ_{e}(r): ρ_{e}(r) is the equilibrium distribution, which can be estimated using data from classical molecular dynamics (MD) or reactive MD before dissociation occurs. On the other hand, ρ_{b}(r) is the configurational distribution under the condition of bond dissociation. Practical estimation of ρ_{b}(r) can be achieved through reactive MD simulations. By conducting numerous independent simulations, each resulting in a configuration r upon the moment of molecular dissociation, ρ_{b}(r) can be then estimated by collecting and analyzing all obtained configurations.

A crucial relationship connecting ρ_{e}(r) and ρ_{b}(r) can be derived as follows:

ρ_{e}(r)λ(r) = λρ_{b}(r) | (8) |

(9) |

From eqn (8), we can derive the expression for λ(r) as follows:

(10) |

(11) |

We could use the reactive MD data to calculate a conditional probability:

(12) |

(13) |

(14) |

By utilizing the highlighted part of eqn (14), we can compute the bond breaking rate constant as a function of both molecular configuration r and bond site position i.

After obtaining λ(r,i), we can also derive the overall bond breaking rate for a specific bond:

(15) |

From eqn (13) we have:

ρ_{e}(r)λ(r,i) = λρ_{b}(r,i) | (16) |

ρ_{b}(r,i) = ρ_{b}(r)P_{b}(i|r) | (17) |

(18) |

(19) |

(20) |

The calculation of λ(i) can be achieved using either eqn (15) and (18) or eqn (20), they all yield identical results. Notably, eqn (20) is relatively more straightforward.

The bond dissociation rate, λ(r,i), is a function of the molecular configuration r, existing in a 3N dimensional space. For molecules with a substantial number of atoms (where N is large), λ(r,i) becomes a high-dimensional function. This intricacy can pose challenges in terms of comprehension and practical application.

In practical scenarios, dealing directly with the 3N dimensional configurational vector r can be unwieldy. Instead, we often rely on lower-dimensional representations of molecular conformation using a set of collective variables denoted as ξ, where ξ ∈ ^{M} and M ≪ 3N. These collective variables, ξ, are chosen to capture the most pertinent and crucial information about the molecular configuration. Consequently, the bond dissociation rate constant is then expressed as λ(ξ,i).

These collective variables can take various forms. Some are straightforward and easily understood, such as the head-to-tail distance of a chain molecule, or they may be closely related to experimental measurements, such as the radius of gyration (R_{g}) for proteins. Apart from those with apparent meanings, modern approaches often employ linear or nonlinear dimension reduction techniques, such as principal component analysis (PCA),^{54} time-lagged independent component analysis (TICA),^{55–57} and diffusion maps (dMaps).^{58–61} These techniques are widely used to extract intrinsic low-dimensional manifolds and collective variables ξ from high-dimensional configurational data.

Molecule |
n (O_{2}) |
T (K) | L (nm) | Density (g cm^{−3}) |
No. of sim. |
---|---|---|---|---|---|

n-Tridecane | 17 | 3500 | 2.5 | 0.0773 | 5000 |

n-Pentane | 9 | 3500 | 2.0 | 0.0747 | 2000 |

1,3-Propanediol | 9 | 3500 | 2.0 | 0.0756 | 2000 |

Reactive molecular dynamics (MD) simulations were conducted utilizing the large-scale atomic/molecular massively parallel simulator (LAMMPS) package.^{62,63} The simulations employed the ReaxFF force field with the CHO-2016 parameter sets developed by van Duin et al.^{28} In ReaxFF, the potential energy of the system (E_{system}) is determined by the bond order, which is a function of the distance between a pair of atoms. The expression for E_{system} is given by the formula:

E_{system} = E_{bond} + E_{over} + E_{under} + E_{lp} + E_{val} + E_{tor} + E_{Coulomb} + E_{vdWaals} | (21) |

Initially, an equilibrium simulation for each molecule at a low temperature (1000 K) was conducted for 10 ps. During this period, no molecular dissociation occurred. Frames extracted from this equilibrium trajectory served as the starting configurations for subsequent oxidation simulations. Specifically, for n-tridecane, 5000 independent oxidation simulations were performed, each with a distinct starting configuration. Similarly, 2000 simulations were conducted for both n-pentane and 1,3-propanediol.

In principle, conducting more independent simulations provides additional sampling data, enhancing result accuracy. However, this also comes with an increased computational cost. In our approach, we opted for a balance between accuracy and computational efficiency. Specifically, we performed 5000 independent simulations for the larger molecule, n-tridecane, and 2000 simulations each for the smaller molecules, n-pentane and 1,3-propanediol. These numbers of simulations ensure accurate results while maintaining acceptable computational costs.

All simulations were carried out in the NVT ensemble using 4 × 5.2 GHz Intel i9-12900KF cores. The system was coupled with a Nosé–Hoover thermostat featuring a damping parameter of 0.01 ps, maintaining a temperature of 3500 K during oxidation simulations. A simulation time step of 0.1 fs was employed, and charge equilibrium calculations were performed at every step.^{64,65}

Given that the adiabatic flame temperature for n-pentane in pure oxygen can reach as high as 3848 K,^{66} we set the reaction temperature for all three molecules to 3500 K to ensure comparability, despite lacking specific combustion flame temperatures data for n-tridecane or 1,3-propanediol. Furthermore, it is worth noting that realistic hydrocarbon combustion temperatures typically fall below 3000 K, notably higher in pure oxygen than in air,^{67} and occur within a few seconds. In contrast, our ReaxFF simulations, like many others in the field, employ higher temperatures. This choice is made to facilitate rapid reactions occurring within a few picoseconds, enabling their capture by molecular dynamics simulations. Despite the disparity in time and temperature between ReaxFF simulations and experiments, the initial reaction mechanisms and kinetics observed in high-temperature ReaxFF simulations show good agreement with experimental data. This approach has been successfully employed by numerous researchers.^{27,42,68}

The criteria for determining C–C or C–O backbone bond breakage in this research are as follows: (a) the bond length must exceed a threshold of 1.9 Å (the average bond length is around 1.5 Å for both bonds), and (b) the bond length must remain greater than 1.9 Å for an entire consecutive 50 fs. The small threshold of 1.9 Å in criterion (a) ensures that we accurately pinpoint the moment the bond breaks. If the threshold is too large, the molecular configuration may significantly deform from the initial moment of bond breakage. The 50 fs duration in criterion (b) prevents ‘fake breaking’, where a bond temporarily exceeds 1.9 Å but quickly rebinds. Criterion (b) ensures that the bond remains broken and that atoms permanently separate.

Due to the exponential distribution governing the lifetime of molecules, the dissociation time in each simulation varies. Imposing a short time limit in the simulation may not guarantee bond breaking, while setting a long simulation time limit might be inefficient as many molecules may dissociate early. To balance these considerations, we implemented the simulation strategy by splitting the entire trajectory into iterative fragments and terminating the iteration once breaking events are identified. The details of this workflow are presented in Algorithm S1 in the ESI.†

Using ReaxFF-MD, we can generate sufficient data within a reasonable timeframe. For higher accuracy, we can utilize quantum mechanical (QM) approaches, such as ab initio or density functional theory (DFT) simulations. Although these methods incur higher computational costs, they produce more accurate trajectories.

In Fig. 2a, the 3D structure of n-tridecane is presented along with the corresponding bond indices. The molecule comprises 12 bonds, with symmetrically assigned bond indices—small indices near the terminals and larger indices at the center. Our focus is exclusively on the backbone, specifically limiting our consideration to C–C bonds within the backbone of n-tridecane.

Fig. 2 (a) The 3D visualization of the n-tridecane along with the bond index. (b) One n-tridecane alongside with 17 oxygen molecules in the simulation box. The dissociation of the n-tridecane occurs approximately at 5.00 ps in Movie S1 (ESI†). |

As detailed in Table 1, we conducted 5000 independent simulations for the dissociation of individual n-tridecane molecules within a box containing 17 O_{2} molecules, as illustrated in Fig. 2b. Each simulation provided data on the molecule's lifetime and the ultimate configuration at the point of dissociation.

The exponential distribution and the calculation of λ, λ(i).
In Fig. 3a, the probability density distribution of lifetimes is presented. Notably, the distribution deviates from an exponential form, featuring a distinct peak around t = 7 ps, indicating a non-equilibrium phase in each oxidation simulation that needs to be excluded, the reason of this non-equilibrium phase is the temperature difference between 3500 K, under which we conduct ReaxFF simulation, and 1000 K, during which, we generate initial configurations. To determine the optimal offset time (Δt) for removal, we process the data by eliminating data with t < Δt. The remaining data, with t ≥ Δt, is shifted leftwards to initiate from the origin by subtracting their lifetimes by Δt: t → t − Δt, which is called the modified data. Subsequently, we compute the average lifetime and the corresponding dissociation rate (λ = 1/) using the modified data. The remaining data is then normalized to obtain the probability density distribution. To assess the fit to an exponential distribution, we calculate the root mean square error (RMSE) between this distribution and the theoretical exponential distribution expressed in eqn (2). A small RMSE indicates a better following of an exponential distribution.

The RMSE for different offsets, Δt, is depicted in the inset of Fig. 3a. Notably, both small and large Δt values yield large errors. The minimum error is observed at Δt = 7.000 ps, as indicated by the red vertical dashed line in Fig. 3a. The probability density distribution for the modified data is presented in Fig. 3b, resulting in an average lifetime of 5.116 ps and λ = 1/ = 0.195 ps^{−1}. This distribution aligns well with the exponential distribution, as depicted by the red curve, and the vertical dashed line indicates the average lifetime.

Utilizing the modified data (all subsequent analyses are conducted using this modified data), we enumerate the number of dissociated molecules with the breaking bond i, where i = 1, 2, …, 6. We then employ eqn (19) to calculate the overall breaking probability for each bond and eqn (20) to determine the corresponding dissociation rate constant. Since the rate constant λ(i) is proportional to P(i), they are depicted together in Fig. 4a. The observed probabilities and bond dissociation rates vary among bond positions, with the terminal bond (bond 1) exhibiting the lowest dissociation rate. This finding aligns with previous DFT studies, indicating that the terminal bond has the highest bond dissociation energy (BDE),^{8,9} and is consistent with the results obtained by Ding et al., showing that the terminal bond dissociates last compared to central bonds,^{31,41} and our results show that the highest dissociation rate appeared between the terminal and the center for bond 3 and 4.

To streamline the analysis of dissociation properties among bonds, and based on the assumption that adjacent C–C bonds possess similar properties and breaking probabilities, which can be seen in Fig. 4a, we introduce a unified bond classification. Bond 1 and 2 are lumped as bond I, representing the ‘terminal’ bond and colored in pink in Fig. 4a and b. Bonds 3 and 4 are amalgamated as bond II, representing the ‘middle’ bond, and is colored in blue, and bonds 5 and 6 form bond III, colored in green, representing the ‘central’ bonds. The probability and dissociation rate for the unified bond are simply the summation of the probability and rate from individual bonds. In Fig. 4b, it is evident that the terminal bond I has the lowest dissociation rate, while bond II has the highest dissociation rate.

R
_{g} as a one-dimensional collective variable.
The radius of gyration, R_{g}, is widely used to describe molecules’ configurational properties, which is defined as the square root of the mass average of s_{i} over all atoms:

where s_{i} indicates the distance of the i-th atom from the center of mass of a molecule, N is the number of atoms and m_{i} is the mass of the i-th atom. For a physical system containing N identical particles, R_{g} can also be calculated as , where g_{i} is the i-th principal moments (eigenvalues) of the gyration tensor S, and

where r^{(k)}_{i} is the i-th cartesian coordinate of position vector r^{(k)} of the k-th particle, and the center of mass of the system is at the origin of the coordinate system, such that .^{69–74}R_{g} is also closely related to experimental measurement using techniques such as static light scattering and small-angle neutron scattering (SNAS).^{69,74,75} A small R_{g} indicates that the molecule is compact and coiled with a diminutive size, while a large R_{g} suggests that the molecule is in a stretched configuration, characterized by a larger size. In the ESI,† we show that R_{g} is well correlated with the top eigenvector of diffusion maps dimensionality reduction analysis,^{58–61} indicating that R_{g} can indeed effectively capture the overall molecular configuration, and span the configurational dependent bond breaking rate, so r could be replaced by R_{g} in all equations: r → R_{g}.

(22) |

(23) |

The calculation of λ(R_{g}) and λ(R_{g},i).
In Fig. 5, the procedures for calculating λ(R_{g},i) are illustrated. Fig. 5a depicts the number counting density of breaking each bond at various configurations, denoted as ρ_{i}(R_{g}), where i represents the bond index (I, II, or III). It is evident that the density for both small and large R_{g} is lower than that for medium R_{g}. Additionally, different bonds exhibit distinct density patterns, with bond I having the smallest density and bond II displaying the largest density for most values of R_{g}.

this probability, denoted as P_{i}(R_{g}), is equivalent to P_{b}(i|r) = P_{b}(i|R_{g}) in eqn (14). The depiction of P_{i}(R_{g}) is shown in Fig. 5b. It is evident that bond I consistently has the lowest probability for all values of R_{g}. When R_{g} is small, bond II exhibits the highest breaking probability, while at large R_{g}, bond III demonstrates the largest probability.

For any value of R_{g}, we could calculate the probability of breaking bond i as follows:

(24) |

The probability densities ρ_{e}(R_{g}) and ρ_{b}(R_{g}) measured from the simulation data are depicted in Fig. 5c, with representative molecular configurations corresponding to small, medium, and large R_{g}. Both ρ_{e}(R_{g}) and ρ_{b}(R_{g}) exhibit bell-shaped distributions, with the peak probability density occurring in the middle of the plot. This suggests that configurations with intermediate R_{g} have the lowest free energy. However, it is important to note that ρ_{e}(R_{g}) and ρ_{b}(R_{g}) are not identical. When R_{g} is small, ρ_{b}(R_{g}) is smaller than ρ_{e}(R_{g}), while for large R_{g}, ρ_{b}(R_{g}) exceeds ρ_{e}(R_{g}). It is noteworthy that ρ_{b}(R_{g}) may lack data points when R_{g} is small due to the rarity of breaking events and the absence of sampling for breaking. In such cases, we treat ρ_{b}(R_{g}) as 0 when R_{g} is very small.

The final step involves using eqn (10) to obtain λ(r) and eqn (14) to obtain λ(r,i), with r being replaced by the collective variable R_{g}. The overall bond dissociation rate as a function of R_{g} is depicted as the black curve in Fig. 5d, while the dissociation rates for each bond as a function of R_{g} are represented by the red, blue, and green curves. It is evident that λ(R_{g}) increases with the augmentation of R_{g}, signifying that the molecule is more likely to dissociate when it is stretched and possesses a larger size (R_{g} = 4.5 Å). Conversely, it tends to maintain stability when coiled up into a smaller-sized configuration (R_{g} = 2.7 Å).

It is important to highlight that although ρ_{b}(R_{g}) reaches its maximum at R_{g} = 3.5 Å (Fig. 5c), this doesn’t necessarily imply that the molecule has the maximum dissociation rate at this particular size. The reason ρ_{b}(R_{g}) is large in the middle is that the equilibrium distribution ρ_{e}(R_{g}) attains its maximum in this region. Even a small dissociation rate results in a significant number of observed molecule dissociations when R_{g} is intermediate. The dissociation rate constant (eqn (10)) is determined by the ratio of ρ_{b}(R_{g}) to ρ_{e}(R_{g}). Therefore, only a relatively smaller ρ_{e}(R_{g}) and a larger ρ_{b}(R_{g}) can imply a substantial dissociation rate.

In the final analysis, a comparison of λ(R_{g},i) among different bonds is presented in Fig. 5d. When R_{g} is small, all λ(R_{g},i) values are small and closely clustered. The zoomed-in plot in the inset illustrates that the dissociation rate λ(R_{g},I) for bond I consistently remains the minimum for all configurations throughout all R_{g}. For R_{g} values smaller than 4.00 Å, λ(R_{g},II) is the largest. However, when R_{g} exceeds 4.00 Å, bond III becomes the most vulnerable, exhibiting the largest dissociation rate.

These results underscore the variability in dissociation rates among different bonds, highlighting their sensitivity to changes in molecular configuration. Stretched molecular configurations favor the breaking of all bonds, with a pronounced preference for the central bond (III). In contrast, coiled configurations favor the breaking of bond II but leads to lower overall dissociation rates.

By obtaining λ, ρ_{e}(R_{g}), ρ_{b}(R_{g}), λ(R_{g},i), and P_{b}(i|R_{g}) using the above procedures, the overall dissociation rate for a bond, λ(i), can be obtained by using eqn (15) or eqn (18). The result is consistent with what is derived from eqn (20) and is presented in Fig. 4. Despite the significantly higher dissociation rate of bond III compared to bonds I and II when R_{g} is large, Fig. 4b indicates that bond II has the largest overall dissociation rate. This arises because ρ_{e}(R_{g}) in the integration (eqn (15)) is maximized in the range of 3.5 Å < R_{g} < 4.0 Å. In this region, even though λ(R_{g},II) is only slightly larger than λ(R_{g},III), the weighted integration amplifies this effect. Consequently, the overall λ(II) surpasses λ(III), despite λ(R_{g},III) being significantly larger when R_{g} is large.

We also conducted pyrolysis simulations of n-tridecane in the gas phase, where no oxygen molecules are present. Similar calculation procedures were carried out, and the results are shown in the ESI.† In the non-oxygen environment, the breaking rate constant for all C–C bonds is higher than in the oxygen environment, and the rate constants are more uniformly distributed over different configurations. A comparison with previous research using DFT and POLYRATE programs is provided,^{53} showing that ReaxFF simulation can qualitatively estimate the bond-breaking rate constant and capture bond-breaking behavior.

In Fig. 6a, the 3D structure of n-pentane and its bond indices are depicted. n-Pentane has four bonds symmetrically distributed along the molecule chain. Bond index 1 corresponds to the terminal bond, while index 2 represents the central bonds. The simulation system, as illustrated in Fig. 6b, consists of one n-pentane molecule placed in a box along with 9 O_{2} molecules.

Fig. 6 (a) The 3D visualization of n-pentane along with the bond index. (b) One n-tridecane alongside 9 oxygen molecules in the simulation box. The dissociation of n-pentane occurs approximately at 7.00 ps in Movie S2 (ESI†). |

Through 2000 independent simulation trajectories of n-pentane, 2000 configurations were generated at the point of molecular dissociation. Similar to the scenario with n-tridecane, the probability density distribution of n-pentane's lifetime is not precisely exponential, as illustrated in Fig. S8a in the ESI,† indicating a non-equilibrium phase at small values of t. Following the same approach used for n-tridecane, we modified the data by excluding the non-equilibrium phase. The optimal offset cut for n-pentane is determined to be 13.600 ps, as also indicated in Fig. S8a (ESI†).

In Fig. 7a, the probability density distribution of lifetime from the modified data is presented, showcasing an exponential distribution. The average lifetime, denoted by and marked by the red vertical dashed line, is determined to be 12.257 ps. The corresponding rate constant, calculated as λ = 1/, is found to be 0.082 ps^{−1}. The red curve represents the exponential distribution ρ(t,λ) = λe^{−λt}, aligning well with the simulated data. Moving to Fig. 7b, which displays the overall probability and dissociation rate for bonds 1 and 2. Unlike the n-tridecane case, n-pentane features only two bonds, obviating the need for defining a unified bond. Notably, bond 1 at the terminal exhibits a significantly lower dissociation rate compared to bond 2 at the center.

Similar to the n-tridecane scenario, we employ R_{g} as a straightforward collective variable to characterize the configuration of n-pentane. A small R_{g} corresponds to a coiled structure, while a large R_{g} indicates that the molecule is stretched, as illustrated in Fig. 7c. Both ρ_{e}(R_{g}) and ρ_{b}(R_{g}) are depicted in Fig. 7c, revealing a notable disparity between these two distributions. Specifically, when R_{g} is small, ρ_{b}(R_{g}) is lower than ρ_{e}(R_{g}), while it surpasses ρ_{e}(R_{g}) significantly when R_{g} is large.

Subsequently, we can derive the counting number density as a function of R_{g} for each bond, denoted as ρ_{i}(R_{g}), and the bond probability as a function of R_{g} for each bond, denoted as P_{i}(R_{g}). These values are illustrated in Fig. S8b and c in the ESI.† By utilizing eqn (10) and (14), we can obtain the ultimate bond dissociation rate for each bond across various configurations. As depicted in Fig. 7d, the black curve represents the dissociation rate for the entire molecule, considering all bonds as a function of R_{g}. The blue and red curves correspond to the dissociation rates for bonds 1 and 2, respectively, as functions of R_{g}. The graph illustrates that when R_{g} is large, the bond dissociation rate experiences a rapid increase, primarily driven by the significantly higher dissociation rate of bond 2 compared to bond 1. However, for extremely small R_{g} (e.g., R_{g} = 1.5 Å), as shown in the zoomed-in inset, bond 1 exhibits a higher dissociation rate. Due to the fact that, for the majority of R_{g} values, bond 2 possesses a higher dissociation rate than bond 1, and the overall rate for bond 2 is larger, as indicated in Fig. 7b.

Pyrolysis simulations of single n-pentane in the gas phase were also conducted, and the results are shown in the ESI.† Similar to the situation with n-tridecane, the bond-breaking rates are higher in the gas phase than when oxygen molecules are present. Together with the results from n-tridecane, both our findings and those from the referenced studies indicate the following: (1) C–C bond dissociation rates are similar among different positions in one alkane chain, with differences being less than one order of magnitude. (2) The bond-breaking rates are consistent in terms of their orders of magnitude, around 10^{11} s^{−1}. (3) Breaking of bond 1 is more difficult than other C–C bonds, with bond 1 having the smallest breaking rate constants. (4) As the alkane chain size increases, bond-breaking rates slightly decrease. These consistencies indicate that ReaxFF simulations can qualitatively estimate bond-breaking rate constants and capture bond-breaking kinetics.

We also explored the C–H bond breaking rate constant in pyrolysis and oxidation simulations. The results show that C–H bonds have significantly lower breaking rate constants compared to C–C bonds. This is because the bond dissociation energy (BDE) for C–H bonds is about 10 kcal mol^{−1} higher than that for C–C bonds, making C–H bonds more robust. Additionally, all C–H bonds are not highly sensitive to the configuration of the n-pentane molecule, and hydrogen atoms at the terminal alkyl (CH_{3}) group are slightly harder to abstract than those bonded to other carbon atoms. These findings are illustrated in the figures in the ESI.†

However, since the breaking of C–H bonds does not change the size of the alkane molecules, it is less significant for alkane oxidation and cracking compared to C–C bond breaking. Therefore, this paper primarily focuses on exploring C–C bond breaking.

Fig. 8 (a) The 3D visualization of 1,3-propanediol along with the bond index. (b) One 1,3-propanediol alongside 9 oxygen molecules in the simulation box. The dissociation of the 1,3-propanediol occurs approximately at 9.00 ps in Movie S3 (ESI†). |

In Fig. 9a, the probability density distribution of the lifetime in the modified data aligns well with the exponential distribution (the distribution of the full data without modification is shown in Fig. S9 in the ESI†). The calculated average lifetime is 12.339 ps, corresponding to a dissociation rate of λ = 1/ = 0.081 ps^{−1}. These values closely resemble those obtained for n-pentane. It is noteworthy that the molecular dissociation rate for 1,3-propanediol and n-pentane is smaller than that for n-tridecane. This discrepancy can be attributed to the smaller number of bonds in the backbone of smaller molecules, resulting in a reduced likelihood of bond breakage.

In Fig. 9b, the overall dissociation rates for bond 1 and bond 2 are illustrated, calculated using eqn (20). Unlike the previous two examples, the terminal C–O bond exhibits a higher dissociation rate, even though the C–O bond has a slightly higher bond dissociation energy (BDE) than the C–C bond.^{76,77} This suggests that the location of a bond and the configurational properties of the molecule play crucial roles in influencing the dissociation rate. The bond dissociation rate is positively correlated with bond dissociation energy (BDE) but is not directly determined by it. Various factors influence bond dissociation rates, including the chemical environment and the molecular configuration. For example, the edge C–O bond, despite having a higher BDE, can experience altered dissociation properties due to interactions with surrounding environmental molecules, such as oxygen. Factors such as bond stretching, torsion, angles, and molecular interactions with the environment depend on specific configurations. These factors collectively determine the overall observed bond dissociation rate.

In Fig. 9c, the probability density distributions ρ_{e}(R_{g}) and ρ_{b}(R_{g}) as functions of R_{g} are displayed. Similar to the results for n-pentane, these two distributions exhibit significant discrepancies. This difference arises from the heterogeneous distribution of dissociation rates across different configurations. By obtaining ρ_{i}(R_{g}) and P_{i}(R_{g}) (depicted in Fig. S9b and c in the ESI†) and utilizing eqn (10) and (14), we can derive the dissociation rate as a function of R_{g} for bond 1 and 2. The results are presented as the blue and red curves in Fig. 9d. The black curve represents the dissociation rate for the entire molecule, denoted as λ(R_{g}), which is the summation of the rate constants for bond 1 and 2: λ(R_{g}) = λ(R_{g},1) + λ(R_{g},2).

In Fig. 9d, it is observed that the dissociation constant for large R_{g} values is higher, and the dissociation of bond 2 is more favored. The zoomed-in inset reveals that small R_{g} values (R_{g} < 1.9 Å) favor the dissociation of bond 1. λ(R_{g},1) is consistently greater than λ(R_{g},2) in this region, leading to an overall λ(1) > λ(2) after integrating over all R_{g} through eqn (15).

Pyrolysis simulations of single 1,3-propanediol in the gas phase were also conducted, with the results shown in the ESI.† The measured bond dissociation rate constants are higher overall, and bond 1 has an even higher rate constant than bond 2.

The lifetime of molecules follows an exponential distribution, which can be extracted from ReaxFF simulation data. The molecular dissociation rate is then calculated as the reciprocal of the average lifetime. Larger molecules, such as n-tridecane, exhibit higher dissociation rates compared to smaller molecules. This is attributed to the larger number of bonds in n-tridecane, increasing the likelihood of bond breakage.

The bond dissociation rate depends not only on the bond position but also on the molecular configuration. Stretched configurations (characterized by a large radius of gyration R_{g}) are associated with higher dissociation rates.

Large R_{g} values are particularly favorable for the breaking of central bonds, which exhibit the highest dissociation rates. However, when R_{g} is small, central bonds do not necessarily have the highest dissociation rate. The terminal bond in n-pentane and 1,3-propanediol, as well as the type II bond in n-tridecane, which lies between the terminal and central bonds, surpass the central bond to possess the largest dissociation rate.

After integrating over all configurations, the overall dissociation rate for the terminal C–O bond in 1,3-propanediol is larger than that for the central C–C bond. Conversely, in n-pentane, the terminal C–C bond has a smaller dissociation rate than the central C–C bond. In n-tridecane, the type II bond, positioned between the terminal and central bonds, exhibits the highest dissociation rate.

This research presents a quantitative ‘top-down’ approach for calculating bond dissociation rates, bypassing the need for detailed reaction pathways, and instead requiring sufficient reactive simulation trajectories to sample all possible reactions effectively. It can be applied to ReaxFF simulation trajectories as demonstrated in this paper. It can also be applied to quantum mechanical (QM) simulation trajectories, potentially offering higher accuracy given sufficient simulation data. This approach provides insights into molecular properties during the oxidation of hydrocarbons in fuel combustion.

Future work could involve applying this approach to more accurate simulation data obtained from MD with ab initio quantum mechanical methods, such as Car–Parrinello molecular dynamics^{78,79} and QM molecular mechanical (QM/MM) methods, or to larger and more complex molecules, and the distribution of bond breaking rates over other aspects of the molecular configuration, such as the helicity, could be studied. Additionally, further exploration could be conducted using variational transition state theory (VTST) to calculate bond-breaking rates in various configurations.

- G. A. Olah and Á. Molnár, Hydrocarbon chemistry, John Wiley & Sons, 2003 Search PubMed .
- I. Glassman, R. A. Yetter and N. G. Glumac, Combustion, Academic press, 2014 Search PubMed .
- M. A. Fahim, T. A. Al-Sahhaf and A. Elkilani, Fundamentals of petroleum refining, Elsevier, 2009 Search PubMed .
- F. D. Mango, The light hydrocarbons in petroleum: a critical review, Org. Geochem., 1997, 26, 417–440 CrossRef CAS .
- P. E. Savage, Mechanisms and kinetics models for hydrocarbon pyrolysis, J. Anal. Appl. Pyrolysis, 2000, 54, 109–126 CrossRef CAS .
- S. R. Turns, Introduction to combustion, McGraw-Hill Companies New York, NY, USA, 1996, vol. 287 Search PubMed .
- G. Badger, Pyrolysis of hydrocarbons, Prog. Phys. Org. Chem., 1965, 3, 1–38 CrossRef CAS .
- J. Wu, H. Ning, L. Ma and W. Ren, Accurate prediction of bond dissociation energies of large n-alkanes using ONIOM-CCSD (T)/CBS methods, Chem. Phys. Lett., 2018, 699, 139–145 CrossRef CAS .
- K. C. Hunter and A. L. East, Properties of C–C bonds in n-alkanes: relevance to cracking mechanisms, J. Phys. Chem. A, 2002, 106, 1346–1356 CrossRef CAS .
- K. Han and T. Chu, Reaction rate constant computations: theories and applications, The Royal Society of Chemistry, 2013 Search PubMed .
- H. DaCosta and F. Maohong, Rate constant calculation for thermal reactions, Wiley Online Library, 2011 Search PubMed .
- J. L. Magee, Theory of the chemical reaction rate constant, Proc. Natl. Acad. Sci. U. S. A., 1952, 38, 764–770 CrossRef CAS PubMed .
- A. Shushin and A. Barzykin, Effect of local molecular shape and anisotropic reactivity on the rate of diffusion-controlled reactions, Biophys. J., 2001, 81, 3137–3145 CrossRef CAS PubMed .
- Y.-P. Chang, K. Długołecki, J. Küpper, D. Rösch, D. Wild and S. Willitsch, Specific chemical reactivities of spatially separated 3-aminophenol conformers with cold Ca
^{+}ions, Science, 2013, 342, 98–101 CrossRef CAS PubMed . - L. Khriachtchev, E. Maç, M. Pettersson and M. Rasanen, Conformational memory in photodissociation of formic acid, J. Am. Chem. Soc., 2002, 124, 10994–10995 CrossRef CAS PubMed .
- C. A. Taatjes, O. Welz, A. J. Eskola, J. D. Savee, A. M. Scheer, D. E. Shallcross, B. Rotavera, E. P. F. Lee, J. M. Dyke, D. K. W. Mok, D. L. Osborn and C. J. Percival, Direct measurements of conformer-dependent reactivity of the Criegee intermediate CH
_{3}CHOO, Science, 2013, 340, 177–180 CrossRef CAS PubMed . - M. Dell’Angela, et al., Real-Time Observation of Surface Bond Breaking with an X-ray Laser, Science, 2013, 339, 1302–1305 CrossRef PubMed .
- J. M. Clough, C. Creton, S. L. Craig and R. P. Sijbesma, Covalent bond scission in the Mullins effect of a filled elastomer: real-time visualization with mechanoluminescence, Adv. Funct. Mater., 2016, 26, 9063–9074 CrossRef CAS .
- K. H. Kim, et al., Direct observation of bond formation in solution with femtosecond X-ray scattering, Nature, 2015, 518, 385–389 CrossRef CAS PubMed .
- D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to applications, Elsevier, 2023 Search PubMed .
- M. Orio, D. A. Pantazis and F. Neese, Density functional theory, Photosynth. Res., 2009, 102, 443–453 CrossRef CAS PubMed .
- D. W. Brenner, Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films, Phys. Rev. B: Condens. Matter Mater. Phys., 1990, 42, 9458 CrossRef CAS PubMed .
- J. Tersoff, New empirical approach for the structure and energy of covalent systems, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 6991 CrossRef PubMed .
- H. S. Johnston and C. Parr, Activation energies from bond energies. I. Hydrogen transfer reactions, J. Am. Chem. Soc., 1963, 85, 2544–2551 CrossRef CAS .
- D. M. Root, C. R. Landis and T. Cleveland, Valence bond concepts applied to the molecular mechanics description of molecular shapes. 1. Application to nonhypervalent molecules of the P-block, J. Am. Chem. Soc., 1993, 115, 4201–4209 CrossRef CAS .
- C. R. Landis, T. Cleveland and T. K. Firman, Valence bond concepts applied to the molecular mechanics description of molecular shapes. 3. Applications to transition metal alkyls and hydrides, J. Am. Chem. Soc., 1998, 120, 2641–2649 CrossRef CAS .
- K. Chenoweth, A. C. Van Duin and W. A. Goddard, ReaxFF reactive force field for molecular dynamics simulations of hydrocarbon oxidation, J. Phys. Chem. A, 2008, 112, 1040–1053 CrossRef CAS PubMed .
- C. Ashraf and A. C. Van Duin, Extension of the ReaxFF combustion force field toward syngas combustion and initial oxidation kinetics, J. Phys. Chem. A, 2017, 121, 1051–1068 CrossRef CAS PubMed .
- M. Kowalik, C. Ashraf, B. Damirchi, D. Akbarian, S. Rajabpour and A. C. Van Duin, Atomistic scale analysis of the carbonization process for C/H/O/N-based polymers with the ReaxFF reactive force field, J. Phys. Chem. B, 2019, 123, 5357–5367 CrossRef CAS PubMed .
- E. M. Kritikos, A. Lele, A. C. van Duin and A. Giusti, A reactive molecular dynamics study of the effects of an electric field on n-dodecane combustion, Combust. Flame, 2022, 244, 112238 CrossRef CAS .
- Y. Liu and J. Ding, Combustion chemistry of n-heptane/ethanol blends: a ReaxFF study, Mol. Simul., 2021, 47, 37–45 CrossRef .
- Z. Chen, W. Sun and L. Zhao, Combustion mechanisms and kinetics of fuel additives: A reaxff molecular simulation, Energy Fuels, 2018, 32, 11852–11863 CrossRef CAS .
- Z. Chen, W. Sun and L. Zhao, Initial mechanism and kinetics of diesel incomplete combustion: ReaxFF molecular dynamics based on a multicomponent fuel model, J. Phys. Chem. C, 2019, 123, 8512–8521 CrossRef CAS .
- M. Fan and Y. Lu, Insights into carbon monoxide oxidation in supercritical H
_{2}O/CO_{2}mixtures using reactive molecular dynamics simulations, J. Supercrit. Fluids, 2022, 189, 105727 CrossRef CAS . - R. F. Goncalves, B. K. Iha, J. A. Rocco and A. E. Kuznetsov, Reactive molecular dynamics of pyrolysis and combustion of alternative jet fuels: A ReaxFF study, Fuel, 2022, 310, 122157 CrossRef CAS .
- Z. Chen, W. Sun and L. Zhao, High-temperature and high-pressure pyrolysis of hexadecane: molecular dynamic simulation based on reactive force field (ReaxFF), J. Phys. Chem. A, 2017, 121, 2069–2078 CrossRef CAS PubMed .
- Y. Wang and G. Liu, Inhomogeneity Effects on Reactions in Supercritical Fluids: A Computational Study on the Pyrolysis of n-Decane, JACS Au, 2022, 2, 2081–2088 CrossRef CAS PubMed .
- C. Ashraf, S. Shabnam, A. Jain, Y. Xuan and A. C. van Duin, Pyrolysis of binary fuel mixtures at supercritical conditions: A ReaxFF molecular dynamics study, Fuel, 2019, 235, 194–207 CrossRef CAS .
- X. Zhang, Y. Pan, Y. Ni, X. Shi and J. Jiang, Atomistic insights into the pyrolysis of methyl ethyl ketone peroxide via ReaxFF molecular dynamics simulation, Process Saf. Environ. Prot., 2022, 161, 316–324 CrossRef CAS .
- Y. Wang, S. Gong, H. Wang, L. Li and G. Liu, High-temperature pyrolysis of isoprenoid hydrocarbon p-menthane using ReaxFF molecular dynamics simulation, J. Anal. Appl. Pyrolysis, 2021, 155, 105045 CrossRef CAS .
- J. Ding, L. Zhang, Y. Zhang and K.-L. Han, A reactive molecular dynamics study of n-heptane pyrolysis at high temperature, J. Phys. Chem. A, 2013, 117, 3266–3278 CrossRef CAS PubMed .
- Q.-D. Wang, J.-B. Wang, J.-Q. Li, N.-X. Tan and X.-Y. Li, Reactive molecular dynamics simulation and chemical kinetic modeling of pyrolysis and combustion of n-dodecane, Combust. Flame, 2011, 158, 217–226 CrossRef CAS .
- K. J. Laidler, The development of the Arrhenius equation, J. Chem. Educ., 1984, 61, 494 CrossRef CAS .
- X. Z. Jiang and K. H. Luo, Reactive and electron force field molecular dynamics simulations of electric field assisted ethanol oxidation reactions, Proc. Combust. Inst., 2021, 38, 6605–6613 CrossRef CAS .
- D. G. Truhlar and B. C. Garrett, Variational transition state theory, Annu. Rev. Phys. Chem., 1984, 35, 159–189 CrossRef CAS .
- J. L. Bao and D. G. Truhlar, Variational transition state theory: theoretical framework and recent developments, Chem. Soc. Rev., 2017, 46, 7548–7596 RSC .
- D. G. Truhlar and B. C. Garrett, Variational transition-state theory, Acc. Chem. Res., 1980, 13, 440–448 CrossRef CAS .
- A. D. Isaacson, D. G. Truhlar, S. N. Rai, R. Steckler, G. C. Hancock, B. C. Garrett and M. J. Redmon, POLYRATE: A general computer program for variational transition state theory and semiclassical tunneling calculations of chemical reaction rates, Comput. Phys. Commun., 1987, 47, 91–102 CrossRef CAS .
- R. Meana-Pañeda, et al., Polyrate 2023: A computer program for the calculation of chemical reaction rates for polyatomics. New version announcement, Comput. Phys. Commun., 2024, 294, 108933 CrossRef .
- J. L. Bao, X. Zhang and D. G. Truhlar, Barrierless association of CF
_{2}and dissociation of C_{2}F_{4}by variational transition-state theory and system-specific quantum Rice–Ramsperger–Kassel theory, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 13606–13611 CrossRef CAS PubMed . - R. Gao, L. Zhu, Q. Zhang and W. Wang, Atmospheric oxidation mechanism and kinetic studies for OH and NO
_{3}radical-initiated reaction of methyl methacrylate, Int. J. Mol. Sci., 2014, 15, 5032–5044 CrossRef PubMed . - Z. Zhang, S. Wang, Y. Shang, J. Liu and S. Luo, Theoretical Study on Ethylamine Dissociation Reactions Using VRC-VTST and SS-QRRK Methods, J. Phys. Chem. A, 2024, 128, 2191–2199 CrossRef CAS PubMed .
- W. Xu and H. Fan, A DFT study on recombination of alkyl radicals to C2-C17 normal alkanes & branched C8 alkanes and corresponding CC bond pyrolysis reaction, Mol. Phys., 2020, 118, e1773002 CrossRef .
- A. Mackiewicz and W. Ratajczak, Principal components analysis (PCA), Comput. Geosci., 1993, 19, 303–342 CrossRef .
- C. R. Schwantes and V. S. Pande, Improvements in Markov state model construction reveal many non-native interactions in the folding of NTL9, J. Chem. Theory Comput., 2013, 9, 2000–2009 CrossRef CAS PubMed .
- L. Molgedey and H. G. Schuster, Separation of a mixture of independent signals using time delayed correlations, Phys. Rev. Lett., 1994, 72, 3634 CrossRef PubMed .
- G. Pérez-Hernández, F. Paul, T. Giorgino, G. De Fabritiis and F. Noé, Identification of slow molecular order parameters for Markov model construction, J. Chem. Phys., 2013, 139, 015102 CrossRef PubMed .
- R. R. Coifman and S. Lafon, Diffusion maps, Appl. Comput. Harmon. Anal., 2006, 21, 5–30 CrossRef .
- R. R. Coifman, S. Lafon, A. B. Lee, M. Maggioni, B. Nadler, F. Warner and S. W. Zucker, Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 7426–7431 CrossRef CAS PubMed .
- R. R. Coifman, I. G. Kevrekidis, S. Lafon, M. Maggioni and B. Nadler, Diffusion maps, reduction coordinates, and low dimensional representation of stochastic systems, Multiscale Model. Simul., 2008, 7, 842–864 CrossRef .
- M. A. Rohrdanz, W. Zheng, M. Maggioni and C. Clementi, Determination of reaction coordinates via locally scaled diffusion map, J. Chem. Phys., 2011, 134, 124116 CrossRef PubMed .
- S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys., 1995, 117, 1–19 CrossRef CAS .
- A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, LAMMPS-a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comput. Phys. Commun., 2022, 271, 108171 CrossRef CAS .
- A. K. Rappe and W. A. Goddard III, Charge equilibration for molecular dynamics simulations, J. Phys. Chem., 1991, 95, 3358–3363 CrossRef CAS .
- A. Nakano, Parallel multilevel preconditioned conjugate-gradient approach to variablecharge molecular dynamics, Comput. Phys. Commun., 1997, 104, 59–69 CrossRef CAS .
- The Engineering ToolBox ( 2005). Adiabatic Flame Temperatures. https://www.engineeringtoolbox.com/adiabatic-flame-temperature-d_996.html.
- The Engineering ToolBox ( 2003). Fuel Gases – Flame Temperatures. https://www.engineeringtoolbox.com/flame-temperatures-gases-d_422.html.
- D. Hong and X. Guo, A reactive molecular dynamics study of CH
_{4}combustion in O_{2}/CO_{2}/H_{2}O environments, Fuel Process. Technol., 2017, 167, 416–424 CrossRef CAS . - M. Fixman, Radius of gyration of polymer chains, J. Chem. Phys., 1962, 36, 306–310 CrossRef CAS .
- A. R. Khokhlov, A. Y. Grosberg and V. S. Pande, Statistical physics of macromolecules, Springer, 1994, vol. 1 Search PubMed .
- P. J. Flory, Principles of polymer chemistry, Cornell University Press, 1953, pp. 428–429 Search PubMed .
- W. L. Mattice and U. Suter, Conformational theory of large molecules: the rotational isomeric state model in macromolecular systems, 1994 Search PubMed .
- D. N. Theodorou and U. W. Suter, Shape of unperturbed linear polymers: polypropylene, Macromolecules, 1985, 18, 1206–1214 CrossRef CAS .
- M. Y. Lobanov, N. Bogatyreva and O. Galzitskaya, Radius of gyration as an indicator of protein structure compactness, Mol. Biol., 2008, 42, 623–628 CrossRef CAS .
- M. Nierlich, F. Boué, A. Lapp and R. Oberthur, Radius of gyration of a polyion in salt free polyelectrolyte solutions measured by SANS, J. Phys., 1985, 46, 649–655 CrossRef CAS .
- E. V. Anslyn and D. A. Dougherty, Modern physical organic chemistry, University Science Books, 2006, p. 72 Search PubMed .
- A. Jacobs, Understanding organic reaction mechanisms, Cambridge University Press, 1997, p. 68 Search PubMed .
- J. Hutter, Car–Parrinello molecular dynamics, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 604–612 CAS .
- T. D. Kühne, Second generation Car–Parrinello molecular dynamics, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 391–406 Search PubMed .

## Footnote |

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cp02271h |

This journal is © the Owner Societies 2024 |