The impact of molecular configuration on the bond breaking rates of hydrocarbons: a computational study†

Jiang Wang *^a, Zhiling Li ^a and Wenli Zhang *^b
aCollege 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
^bSchool of Transportation Engineering, Guizhou Institute of Technology, Boshi Road, Dangwu Town, Gui’an New District, Guizhou 550025, China. E-mail: wenlizhang@git.edu.cn

First published on 16th August 2024

---

---

1 Introduction

Hydrocarbons are a crucial form of energy resource for humanity. The combustion (oxidation) of these organic compounds, including alkanes present in natural gas or emerging clean fuels like ethanol, yields substantial energy and heat.^1,2 Simultaneously, in the realm of petroleum refinement, the transformation of large hydrocarbon molecules from crude oil involves their breakdown into smaller components through pyrolysis processes.^3–5 Central to both energy extraction and industrial modification is the dissociation of covalent bonds within the molecular backbones of these compounds.^6,7 A comprehensive understanding of bond-breaking energies and dissociation rates becomes imperative for unraveling the intricate mechanisms underpinning these processes. Whether through combustion for energy release or thermal cracking for industrial utilization, a deeper insight into these molecular transformations is not only crucial for fundamental comprehension but also essential for the practical application of hydrocarbons.

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.²⁰ On the other hand, quantum mechanical approaches (QM) can determine the properties of electrons in molecules and is widely used to simulate chemical reactions²¹ 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.⁴¹ This observation aligns with similar results reported by Ding et al.³¹ 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.⁴³ 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.

2 Theory and methods

In this section, we will elucidate key elements of the theory and methodologies that underpin the eventual computation of bond dissociation rates, which depend on the molecular configuration and the specific positioning of bonds within the molecule.

2.1 The exponential distribution of molecular lifetime

For a molecule containing N_b bonds in its backbone, dissociation occurs when any one of these bonds undergoes breaking. We explicitly exclude scenarios where multiple bonds dissociate simultaneously, as our reactive molecular simulation employs a small time resolution—specifically, an integration time step of 0.1 fs. In such fine-grained temporal settings, the dissociation probability for a single bond at one step becomes infinitesimal. Consequently, the likelihood of multiple bonds breaking concurrently is negligible, given that it involves the multiplication of multiple infinitesimal probabilities.

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

	(1)

where N(t) represents the count of all remaining molecules at time t that have not undergone dissociation. In this context, when dealing with a single molecule, λ can be interpreted as follows: λdt signifies the probability that this singular molecule will dissociate within the next time interval of dt.

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)

and the lifetime t is defined as the time interval between a designated starting point and the moment of dissociation. It is noteworthy that the choice of the starting time is arbitrary: if t is shifted by any translation Δt, the exponential distribution remains invariant after normalization.

The average lifetime can be computed as follows:

	(3)

which implies that the average lifetime of a molecule, denoted as , is the reciprocal of its dissociation rate. Measurement of the average lifetime can be easily achieved through reactive molecular simulations. Consequently, the dissociation rate can be straightforwardly obtained as:

	(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.

2.2 Molecular dissociation rate subject to configurations

Several experimental studies reveal that molecular conformation influences the reaction rate.^13–16 Because the stretching of bonds, angles, torsions, and how atoms are interacting with the environment all depend on the configuration of the molecule, this suggests that the dissociation rate constant λ can be expressed as a function of molecular configuration, denoted as λ(r). Here, r ∈ ^3N represents the molecule's conformation in the 3N dimensional configurational space, where N denotes the number of atoms in the molecule.

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

	(5)

Here, ρ_e(r) denotes the equilibrium probability density distribution of r in the configurational space, with . The subscript ‘e’ signifies ‘equilibrium’. Q serves as the normalization denominator, representing the canonical partition function, while F(r) stands for the free energy of the molecule. Additionally, k_B denotes the Boltzmann constant.

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

	(6)

A concise proof can be presented as follows: considering a single molecule, it has the probability of ρ_e(r)dr to exist within the domain dr at location r in the configurational space. The probability that it dissociates within the future time interval dt is λ(r)dt. Therefore, the probability for a molecule, at any configuration, to dissociate within the future dt is given by the integration of ρ_e(r)drλ(r)dt over all configurations in ^3N, and by definition, this probability equals λdt, leading to the relationship:

	(7)

upon canceling dt and reordering terms on the left-hand side of the equation, we establish the validity of eqn (6).

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)

a concise proof is as follows: ρ_b(r)dr signifies the probability of a molecule dissociating and being within the domain dr at location r under the condition that the molecule dissociates in any configuration. Utilizing the definition of conditional probability, we have:

	(9)

where P(r,b) represents the probability of a molecule being at r + dr and dissociating within the future dt. This is given by ρ_e(r)drλ(r)dt. The term P(b) denotes the probability of the molecule dissociating in any configuration within the future dt, which is equivalent to λdt. By utilizing eqn (9), it is straightforward to demonstrate eqn (8) by canceling out dr and dt.

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

	(10)

and all terms on the right-hand side can be obtained through reactive molecular simulations. This results in the molecular dissociation rate constant, λ(r), as a function of molecular configuration.

2.3 Bond breaking rate subject to positions

The molecular dissociation rate, λ(r), encapsulates the collective breaking properties of all bonds. However, to delve into the breaking dynamics of individual bonds, we aim to calculate the breaking rate constant for each bond, denoted as λ(r,i), where i is the bond index, and i = 1, 2, …, N_b, with N_b being the total number of bonds within a molecule. This leads us to:

	(11)

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

	(12)

Here, P_b(i|r) represents the probability that the broken bond is i under the condition that the molecule dissociates at the configuration r. On the right-hand side, ρ_b(r,i) signifies, under the condition of molecular dissociation, the probability density for the molecule with configuration r to break its i-th bond. Utilizing eqn (8), we obtain:

	(13)

substituting eqn (13) into eqn (12), we get:

	(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)

and from eqn (12) we have:

ρb(r,i) = ρb(r)Pb(i|r)	(17)

substituting eqn (16) and (17) into eqn (15), we have:

	(18)

which implies that the overall bond dissociation rate for bond i, denoted as λ(i), is obtained by multiplying λ with the overall breaking probability for bond i, represented as P(i), and

	(19)

where N_k denotes the number of molecule dissociations involving the breaking of bond k as measured from ReaxFF simulations. Then we can have an alternative expression of λ(i) as:

	(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),⁵⁴ 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.

2.4 Reactive molecular dynamics simulation

In this study, we investigated the dissociation processes of three distinct molecules: n-tridecane (CH₃(CH₂)₁₁CH₃), n-pentane (CH₃(CH₂)₃CH₃), and 1,3-propanediol (CH₂(CH₂OH)₂). n-Tridecane possesses 13 backbone atoms and 12 backbone bonds, while n-pentane and 1,3-propanediol have 5 backbone atoms and 4 backbone bonds. The chemical structures of these molecules are depicted in Fig. 1. Initially, one of each target molecule was placed within a simulation box under periodic boundary conditions, accompanied by varying numbers of O₂ molecules to ensure comparable densities across simulations. The setup parameters for the simulations are detailed in Table 1.

---

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.²⁸ 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:

Esystem = Ebond + Eover + Eunder + Elp + Eval + Etor + ECoulomb + EvdWaals	(21)

where, the first six terms on the right-hand side of the formula represent specific energy contributions related to bond formation and breaking: bond energy, overcoordination energy penalty, undercoordination stability, long pair energy, valence angle energy, and torsion angle energy. These terms capture the intricacies of bond dynamics in the system. The last two terms, E_Coulomb and E_vdWaals, account for electrostatic and van der Waals potentials, respectively.^28,29

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,⁶⁶ 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,⁶⁷ 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.

3 Results and discussions

In this section, we will unveil the outcomes pertaining to three typical hydrocarbon molecules and delineate the methodology employed for calculating their bond dissociation rates.

3.1 n-Tridecane

In this part, we will use n-tridecane as our initial example, providing detailed results from key steps leading to the ultimate calculation of the bond dissociation rate.

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.

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

3.2 n-Pentane

In this part, we will illustrate the calculation of the bond dissociation rate for the n-pentane molecule. The discussion will be more concise compared to the n-tridecane part, focusing on notable features specific to n-pentane while assuming familiarity with the overall procedures outlined earlier.

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₂ molecules.

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⁻¹. 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¹¹ s⁻¹. (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⁻¹ 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₃) 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.

3.3 1,3-Propanediol

In this part, we will present the results for the last example: 1,3-propanediol. Its 3D configuration, together with bond indexing, is shown in Fig. 8a. Unlike the previous two examples whose backbones consist solely of C–C bonds, in 1,3-propanediol, the terminal bond is C–O, which will exhibit different properties. Fig. 8b displays the simulation system containing one 1,3-propanediol and 9 O₂ molecules.

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⁻¹. 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.

4 Conclusions

In this research, we introduced a ‘top-down’ method to calculate bond dissociation rates over different molecular configuration, and applied it to the ReaxFF molecular dynamics simulation of three different molecules: n-tridecane, n-pentane, and 1,3-propanediol. The following conclusions were drawn.

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.

Data availability

Codes for this article, including the script for generating simulation trajectories and the code for conducting bond breaking rate calculation, are available at Github repository REAX_BRR at https://github.com/cwangjiang/REAX_BBR.git.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This material is based upon work supported by Guizhou Provincial Science and Technology Projects (Grant No. Qian Ke He Ji Chu-ZK [2022] Yi Ban 183), and the High-Level Talent Startup Fund at Guizhou Institute of Technology (2023GCC054).

References

1. G. A. Olah and Á. Molnár, Hydrocarbon chemistry, John Wiley & Sons, 2003 Search PubMed .
2. I. Glassman, R. A. Yetter and N. G. Glumac, Combustion, Academic press, 2014 Search PubMed .
3. M. A. Fahim, T. A. Al-Sahhaf and A. Elkilani, Fundamentals of petroleum refining, Elsevier, 2009 Search PubMed .
4. F. D. Mango, The light hydrocarbons in petroleum: a critical review, Org. Geochem., 1997, 26, 417–440 CrossRef CAS .
5. P. E. Savage, Mechanisms and kinetics models for hydrocarbon pyrolysis, J. Anal. Appl. Pyrolysis, 2000, 54, 109–126 CrossRef CAS .
6. S. R. Turns, Introduction to combustion, McGraw-Hill Companies New York, NY, USA, 1996, vol. 287 Search PubMed .
7. G. Badger, Pyrolysis of hydrocarbons, Prog. Phys. Org. Chem., 1965, 3, 1–38 CrossRef CAS .
8. 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 .
9. 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 .
10. K. Han and T. Chu, Reaction rate constant computations: theories and applications, The Royal Society of Chemistry, 2013 Search PubMed .
11. H. DaCosta and F. Maohong, Rate constant calculation for thermal reactions, Wiley Online Library, 2011 Search PubMed .
12. J. L. Magee, Theory of the chemical reaction rate constant, Proc. Natl. Acad. Sci. U. S. A., 1952, 38, 764–770 CrossRef CAS PubMed .
13. 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 .
14. 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 .
15. 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 .
16. 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₃CHOO, Science, 2013, 340, 177–180 CrossRef CAS PubMed .
17. M. Dell’Angela, et al., Real-Time Observation of Surface Bond Breaking with an X-ray Laser, Science, 2013, 339, 1302–1305 CrossRef PubMed .
18. 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 .
19. 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 .
20. D. Frenkel and B. Smit, Understanding molecular simulation: from algorithms to applications, Elsevier, 2023 Search PubMed .
21. M. Orio, D. A. Pantazis and F. Neese, Density functional theory, Photosynth. Res., 2009, 102, 443–453 CrossRef CAS PubMed .
22. 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 .
23. 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 .
24. 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 .
25. 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 .
26. 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 .
27. 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 .
28. 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 .
29. 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 .
30. 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 .
31. Y. Liu and J. Ding, Combustion chemistry of n-heptane/ethanol blends: a ReaxFF study, Mol. Simul., 2021, 47, 37–45 CrossRef .
32. 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 .
33. 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 .
34. M. Fan and Y. Lu, Insights into carbon monoxide oxidation in supercritical H₂O/CO₂ mixtures using reactive molecular dynamics simulations, J. Supercrit. Fluids, 2022, 189, 105727 CrossRef CAS .
35. 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 .
36. 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 .
37. 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 .
38. 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 .
39. 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 .
40. 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 .
41. 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 .
42. 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 .
43. K. J. Laidler, The development of the Arrhenius equation, J. Chem. Educ., 1984, 61, 494 CrossRef CAS .
44. 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 .
45. D. G. Truhlar and B. C. Garrett, Variational transition state theory, Annu. Rev. Phys. Chem., 1984, 35, 159–189 CrossRef CAS .
46. J. L. Bao and D. G. Truhlar, Variational transition state theory: theoretical framework and recent developments, Chem. Soc. Rev., 2017, 46, 7548–7596 RSC .
47. D. G. Truhlar and B. C. Garrett, Variational transition-state theory, Acc. Chem. Res., 1980, 13, 440–448 CrossRef CAS .
48. 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 .
49. 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 .
50. J. L. Bao, X. Zhang and D. G. Truhlar, Barrierless association of CF₂ and dissociation of C₂F₄ 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 .
51. R. Gao, L. Zhu, Q. Zhang and W. Wang, Atmospheric oxidation mechanism and kinetic studies for OH and NO₃ radical-initiated reaction of methyl methacrylate, Int. J. Mol. Sci., 2014, 15, 5032–5044 CrossRef PubMed .
52. 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 .
53. 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 .
54. A. Mackiewicz and W. Ratajczak, Principal components analysis (PCA), Comput. Geosci., 1993, 19, 303–342 CrossRef .
55. 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 .
56. 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 .
57. 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 .
58. R. R. Coifman and S. Lafon, Diffusion maps, Appl. Comput. Harmon. Anal., 2006, 21, 5–30 CrossRef .
59. 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 .
60. 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 .
61. 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 .
62. S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys., 1995, 117, 1–19 CrossRef CAS .
63. 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 .
64. A. K. Rappe and W. A. Goddard III, Charge equilibration for molecular dynamics simulations, J. Phys. Chem., 1991, 95, 3358–3363 CrossRef CAS .
65. A. Nakano, Parallel multilevel preconditioned conjugate-gradient approach to variablecharge molecular dynamics, Comput. Phys. Commun., 1997, 104, 59–69 CrossRef CAS .
66. The Engineering ToolBox ( 2005). Adiabatic Flame Temperatures. https://www.engineeringtoolbox.com/adiabatic-flame-temperature-d_996.html.
67. The Engineering ToolBox ( 2003). Fuel Gases – Flame Temperatures. https://www.engineeringtoolbox.com/flame-temperatures-gases-d_422.html.
68. D. Hong and X. Guo, A reactive molecular dynamics study of CH₄ combustion in O₂/CO₂/H₂O environments, Fuel Process. Technol., 2017, 167, 416–424 CrossRef CAS .
69. M. Fixman, Radius of gyration of polymer chains, J. Chem. Phys., 1962, 36, 306–310 CrossRef CAS .
70. A. R. Khokhlov, A. Y. Grosberg and V. S. Pande, Statistical physics of macromolecules, Springer, 1994, vol. 1 Search PubMed .
71. P. J. Flory, Principles of polymer chemistry, Cornell University Press, 1953, pp. 428–429 Search PubMed .
72. W. L. Mattice and U. Suter, Conformational theory of large molecules: the rotational isomeric state model in macromolecular systems, 1994 Search PubMed .
73. D. N. Theodorou and U. W. Suter, Shape of unperturbed linear polymers: polypropylene, Macromolecules, 1985, 18, 1206–1214 CrossRef CAS .
74. 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 .
75. 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 .
76. E. V. Anslyn and D. A. Dougherty, Modern physical organic chemistry, University Science Books, 2006, p. 72 Search PubMed .
77. A. Jacobs, Understanding organic reaction mechanisms, Cambridge University Press, 1997, p. 68 Search PubMed .
78. J. Hutter, Car–Parrinello molecular dynamics, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2012, 2, 604–612 CAS .
79. 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