Reax-Lump: an automated method for generating lumped kinetics of polymer plastic pyrolysis from reactive molecular dynamics simulation

Plastic products have become an essential component in our modern society, with wide applications such as construction, packaging, healthcare, and industrial production^[1−3]. According to Plastics Europe^[4], global plastic production reached 430.9 million tons in 2024, of which polypropylene (PP) and polyethylene (PE) accounted for 19% and 26%, respectively. Due to their favorable properties, these two plastics are widely used across various fields^[5−7]. However, high production of plastic products stands in stark contrast to the low recycling rate^[8], as most are incinerated or landfilled, resulting in significant resource waste and severe ecological pollution^[9,10]. In recent years, pyrolysis has been proposed as a promising method for valorizing polymer plastics^[11−13]. Pyrolysis thermally decomposes long-chain polymers into smaller molecules under controlled temperatures and pressures, yielding oil, gas, and char. Studies indicate that at a moderate temperature of around 500 °C, the mass fraction of liquid oil can reach up to 80%^[14]. Moreover, pyrolysis offers operational flexibility, allowing product yields to be optimized by adjusting the process parameters to meet different industrial needs. But product distribution and conversion pathways of pyrolysis are influenced by multiple factors, including temperature, reaction time, catalyst, reactor, etc.^[15]. Revealing reaction mechanism and kinetics is essential for optimizing reaction conditions and designing efficient reactors.

Extensive experimental investigations have been carried out to explore the pyrolysis process of pure PP and PE, as well as their mixtures. Conventional techniques such as X-ray diffraction (XRD), transmission electron microscopy (TEM), and thermogravimetric analysis (TG)^[16−18] can only provide characterization of stable products after pyrolysis, failing to capture dynamic reaction evolution. Based on the results of thermogravimetric analysis (TGA), Schubert et al.^[19,20] proposed a lumped reaction mechanism involving four species (coke, heavy oil, light oil, and gas) and six reactions, which has been used in numerical simulations to predict the product distribution. However, as a simplified lumped mechanism derived purely from regression of experimental data, it only captures inter-species mass transfer, and is insufficiently detailed for more refined and in-depth investigations. Reactive force field (ReaxFF) is an empirical bond-order-based reactive force field introduced by van Duin et al.^[21] in 2001. By integrating ReaxFF with molecular dynamics (MD), the method enables continuous tracking of bond formation and breakage during the pyrolysis process without predefined reaction pathways. Compared with density functional theory (DFT) methods, ReaxFF MD achieves reasonable accuracy at a significantly lower computational cost. With the advantage of the unique strength in modeling condensed-phase systems, ReaxFF MD has been widely applied to investigate the pyrolysis processes of large and complex reactive systems for polymer plastics. Zhang et al.^[22] constructed large-scale atomic models containing over 80,000 atoms and employed ReaxFF MD simulations with experiments to investigate the co-pyrolysis of lignin and polyethylene. Their results indicated that adding polyethylene inhibited the formation of hydrocarbon tar and char, while promoting the production of oxygenated tar and light phenols.

Recently, some post-processing analysis tools for ReaxFF MD simulation have been developed^[23,24]. Among them, VARxMD^[25] has been created in our previous work, a reaction analysis and visualization tool that allows the extraction of a detailed list of reactions with reaction sites from ReaxFF MD simulations.

It is well recognized that DFT methods provide relatively high accuracy, but are limited to small systems^[26] (~100 atoms) and cannot fully capture the complexity of polymer pyrolysis, while experimental techniques such as TGA offer macroscopic kinetic data, but lack atomistic reaction mechanisms. Machine learning models fitted to experimental data can yield compact kinetic models, but often miss chemical interpretability. ReaxFF MD simulations bridge these gaps by enabling large-scale atomistic simulations with reasonable accuracy and providing detailed reaction trajectories. However, a significant challenge currently faced by ReaxFF MD in modeling complex systems is the excessive number of reactions and species in the simulation data, which complicates the extraction of useful information and skeleton kinetic models. For example, Li et al.^[27] simulated PP pyrolysis at 3,000 K for 1 ns and obtained 395,929 reactions, which are too extensive to be directly used in downstream computational fluid dynamics (CFD) numerical simulations. In general, a simple and chemically interpretable kinetic mechanism is needed to significantly reduce the computational cost of downstream CFD simulations. Moreover, it facilitates the interpretation of simulations from a fundamental reaction perspective, and supports the selection of appropriate kinetic models for reactor design or industrial optimization^[28]. The conventional algorithms, such as sensitivity and rate analyses, directed relation graph (DRG), and directed relation graph with error propagation (DRGEP)^[29−32], have been proposed to simplify complex reaction mechanisms, especially for hydrocarbon fuel systems. But the highly scattered data and diverse species in polymer plastics pyrolysis generated by ReaxFF MD pose great difficulties in convergence, statistical reliability, and computational cost, limiting the capability of conventional reduction methods.

Some studies have adopted novel approaches to investigate the reaction mechanisms and kinetic models of plastic pyrolysis. Shen et al.^[33] creatively adopted a chemical reaction neural network (CRNN) approach to directly fitting species concentration profiles from experiments. This method yielded locally optimal solutions expressed as reaction equations. But the number of resulting mechanisms depends on the number of hidden layers in the neural network, which lacks chemical interpretability. Leveraging the fact that the number of chemical reaction categories is finite, Yang et al.^[34,35] developed a machine learning-based reaction classification model (SRG-Reax), providing valuable new insights into the simplification of very complex reaction mechanisms. Despite the progress, the approach does not yet reduce species numbers on a network scale, and lacks a scheme for kinetic parameter extraction.

Reaction classification is an important strategy for extracting reaction information and constructing kinetic models. However, modeling the polymer plastic pyrolysis system is challenging due to its significant reaction diversity and pronounced low-frequency characteristics. Particularly, no chemically interpretable and statistically robust reduced kinetic model directly derived from ReaxFF MD simulation data is currently available for polymer pyrolysis. Thus, the objective of this work is to develop a method for constructing chemically interpretable and reliable lumped kinetic models for polymer pyrolysis directly from ReaxFF MD simulation data, which shows potential for application in downstream numerical simulations. The strategy of species clustering combined with further lump reactions is proposed to improve their statistical availability, which would serve as a sound basis for kinetic model development. Leveraging the detailed chemical pathways and reaction categories obtained from our previous work of VARxMD^[25] and SRG-Reax^[34,35] method, this work introduces a first attempt of Reax-Lump to construct a lumped kinetic model of plastics directly from the ReaxFF MD simulation. A key advantage of the mechanism generated with SRG-Reax is its inherent chemical interpretability. Importantly, a collaborative optimization framework integrating genetic algorithms (GA) and trust region optimization algorithms was proposed innovatively during the mechanism lumping process, which effectively addresses the dual challenges of large-scale reaction mechanism reduction and kinetic parameter estimation. It should be noted that compared with traditional chemical reaction mechanism reduction methods, Reax-Lump proposed in this work exhibits good applicability, satisfactory chemical interpretability, and high-efficiency mechanism reduction capability for plastic pyrolysis systems (Supplementary Table S1). The details of the relevant algorithms underlying this methodology are described in Methods. In Results and discussion, PP is used as a case study to perform the proposed strategy application. The final conclusions and discussions are presented in the last section.

The Reax-Lump approach for obtaining lumped mechanisms consists of four steps, which are shown in Fig. 1. First, a large-scale reaction space, including detailed reactions and species for the pyrolysis of polymer plastics, is obtained by using ReaxFF MD simulation and VARxMD analysis^[25]. To endow the final mechanism with chemical interpretability, the SRG-Reax method is subsequently employed, which provides a new insight into the mechanism from the perspective of reaction classification. In the third step, the detailed reactions are mapped onto lumped reactions by defining rules that assign detailed species to lumped species. This further reduces the reaction network from tens of thousands of individual reactions to several hundred, significantly enhancing its statistical reliability. Finally, the kinetic parameters are optimized using a collaborative training method that combines the genetic algorithm^[36] and the trust region^[37] optimization strategy.

Figure 1. 

Processing scheme of Reax-Lump for generating lumped mechanisms of polymer plastics based on ReaxFF MD simulations.

ReaxFF MD simulation and reaction recognition

This work adopted PP as the representative model for polymer plastics. The PP molecular model was constructed by Li et al.^[27] , which consists of 24 PP molecular chains (each with n = 150), with a total atomic count of 32,448 and a density of 0.74 g/cm³. It should be noted that energy minimization preprocessing was performed during the geometry optimization step to ensure a stable initial configuration for the subsequent ReaxFF MD simulations. According to the Arrhenius equation, shown in Eq. (1), the determination of the two kinetic parameters activation energy (E_a) and pre-exponential factor (lnA) requires fitting data from multiple temperatures that can be from isothermal ReaxFF MD simulations, or a single heat-up simulation with reasonable a heating rate^[33]. However, conducting numerous isothermal simulations to guarantee reliable statistics for each lumped reaction is computationally prohibitive. To avoid the high computational cost for extensive isothermal simulations, a series of heat-up ReaxFF MD simulations were conducted under the heating rates of 1, 2, 4, and 8 K/ps from 300 to 3,000 K. In addition, our previous work also indicated that 'heat-up' simulations with a sufficiently low heating rate (e.g., 2 K/ps for coal and PE)^[22,38] proved effective in capturing the initial pyrolysis mechanisms and the overall evolution stages. Simulations were performed using the Berendsen thermostat to control temperature, with a time-step of 0.25 fs under the NVT ensemble, and employing the CHO-2008 reactive force field^[39]. Post-simulation processing, as well as the identification and acquisition of the global reaction list, were performed using the software VARxMD^[25].

Improvement of machine learning classification method in plastic pyrolysis

It should be noted that the reliability of SRG-Reax has been previously verified in systems such as n-dodecane and jet fuels of RP-1 and RP-3^[40,41], where the reaction types are provided in Supplementary Table S2. However, two limitations arise when applying SRG-Reax to polymer plastic pyrolysis. First, the absence of reaction labels specific to plastic pyrolysis leads to poor generalization and unsatisfactory performance. Second, the difficulty in observing pyrolysis reactions of large molecules poses a challenge in the annotation of reaction labels. To overcome these issues, this work developed a reaction center labeling platform, which is shown in Supplementary Fig. S1. The platform enabled the extraction of local reaction information from the pyrolysis reactions of large molecular fragments. Obviously, the reaction category is determined by the reaction center, adjacent atoms, and their neighboring functional groups. In this work, the reaction center extraction algorithm is adopted from Yang et al.^[35], which extracts the reaction structural information within four topological distances around the reactive atom, where each functional group is treated as one topological distance. The structures of functional groups considered in the extraction of reaction centers are shown in Supplementary Table S3 and annotation rules for reaction categories are provided in Supplementary Table S2. The model's performance was evaluated using the F₁-score for multi-class classification, while leveraging a confusion matrix was constructed to guide the labeling strategy and error correction scheme.

Using the developed platform, 1,320 reaction data entries from the pyrolysis of PP at a heating rate of 2 K/ps were labeled. Among these entries, 924 were used for model training and 396 for validation. As shown in Fig. 2a, when trained solely based on the original dataset, the model achieved a mere 51% accuracy with an F₁-Score of 0.52. Performance was particularly low for reaction categories 1, 2, 5, and 23, which might be attributed to the differences caused by molecular size. The original dataset is derived from fuel systems of n-dodecane, RP-1, and RP-3^[40,41], where initial species contain no more than 20 carbon atoms, whereas the polymer plastic model constructed in this work contains about 300 carbon atoms. As shown in Table 1, there is a significant difference in the predictive performance between the initial models containing high and low carbon number reactions. After incorporating the PP pyrolysis reaction dataset with 924 entries, the prediction accuracy of the model and the F₁-Score of the model on the validation set has significantly improved, which is shown in Fig. 2b. Thus, SRG-Reax embedded with the new PP dataset was employed for the reaction classification of PP pyrolysis simulation results under different heating rates, to investigate the entire PP pyrolysis process from the perspective of reaction classification.

Figure 2. 

Confusion matrices of the reaction classification model: (a) the model training set consists of the original dataset (n-dodecane, RP-1, RP-3); (b) the model training set is augmented with an additional 924 data points from the PP pyrolysis dataset.

Reaction lumping and reaction vertical merging

It is noteworthy that a key advantage of ReaxFF MD simulation is the access to detailed atomistic species information, which gives Reax-Lump the ability of flexible adjustment of these lumping rules based on specific target products. Specifically, detailed species are classified into categories including polypropylene (PP), primary pyrolysis products of PP, heavy oil, light oil, wax, and gas (shown in Table 2), which is consistent with those commonly adopted in literature, and compatible with downstream simulation software such as Aspen^[33,42]. In addition, species with major typical gases (C₂H₆, C₂H₄, CH₄, H₂) were analyzed separately.

Following species lumping, the reaction network is further simplified through reaction mapping and merging. For example, the detailed reaction C₃₀₀H₆₀₂ → C₁₀₀H₃₀₂ + C₁₀₀H₃₀₀ is mapped to [C₂₀₀₊] → 2[C₄₁~C₂₀₀]. To further reduce the reaction scale and decrease the complexity of subsequent kinetic parameter optimization, we performed time-scale merging on the reaction list, that is also called vertical merging. If a product at this step acts as a reactant at the next step, these two sequential reactions are merged into a single overall reaction.

Lumping kinetic parameter optimization

According to the lumping rules described in Table 2, lumping treatment was performed on the PP pyrolysis system with heating rates of 1, 2, 4, and 8 K/ps. Among them, under the conditions of a heating rate of 2 K/ps from 300 to 3,000 K and an analysis interval of 2.5 ps, a total of 19,831 reactions were detected, and these reactions were lumped into 1,881 lumped reactions. Although the dimension of the reaction space has been reduced by an order of magnitude, it is still quite large, which is neither conforming to the reaction scale required for lumped kinetics, nor avoiding posing a significant challenge to the calculation of kinetic parameters. This study achieves reaction mechanism reduction and parameter optimization via optimization methods based on species concentration evolution. A genetic algorithm is adopted for mechanism reduction, as it exhibits strong search capability in combinatorial optimization problems^[43]. For parameter optimization, the trust region optimization algorithm is employed, which does not require calculating gradients of the original optimization function and demonstrates excellent convergence performance for complex optimization problems. As shown in Eq. (2), the optimization objective L is defined as the sum of the mean squared errors of each species between the predicted concentrations from the final mechanism and the species concentrations obtained from ReaxFF MD simulations, where activation energy of Eₐ and pre-exponential factor logarithm of lnA of each reaction serve as the optimization parameters.

where, C₀ is the initial species concentration vector, $ C_{\exp}\left(t_i;E_a,\mathrm{ln}A\right) $ is the predicted concentration at time t_i, $ {C}_{\exp }\left({t}_{i}\right) $ is the observed concentration in ReaxFF MD simulation at time t_i. $ f(C\left(t\right),{E}_{a},\ln A) $ represents the rate law describing the reaction rate as a function of concentration, E_a, and $ \ln A $. Here, m represents the number of reactions, and N represents the number of moments.

Notably, this optimization problem is characterized by three major difficulties. Firstly, the space of initial optimization parameters involves thousands of variables, which results in a high dimensionality, leading to high computational costs and unstable convergence. Secondly, the loss function evaluation is highly complex, which requires integrating over the entire reaction process to determine species concentrations. Conventional gradient-based optimization algorithms often suffer from gradient instability in gradient computation. Lastly, a trade-off exists between the scale and the accuracy of the lumped mechanism. While combinatorial optimization algorithms can handle a large-scale mechanism, further simplification is constrained by accuracy requirements. Moreover, the large-scale mechanism also leads to difficulties in convergence during numerical optimization.

By leveraging the combinatorial optimization capability of GA and the convergence stability of the trust region optimization algorithm, this study proposes a collaborative training framework, which is illustrated in Fig. 3. An initial estimation for kinetic parameters is required after obtaining the lumped reaction list. A reaction kinetic calculation formula proposed by Kröger et al.^[44], which is used successfully in elementary reactions, is introduced to estimate the reaction rate constant (k) in this work. Subsequently, the E_a of the reactions are initialized randomly between 100–300 kJ/mol, and the corresponding $ \ln A $. is calculated by the Arrhenius equation, shown in Eq. (2). Subsequently, shown as Step II in Fig. 3, the original mechanism is optimized using an iterative loop combining GA and the trust-region method. The mechanism reduced by GA serves as the input for the trust-region optimization, and the output from the trust-region optimization is fed back into the GA. This cycle ultimately achieves the dual goals of reducing the mechanism's size and enhancing its accuracy. In this work, GA is applied for mechanism simplification. Each individual vector in the population represents the selection or non-selection of a particular reaction from the full reaction list. The fitness function of Eq. (3) includes the error caused by concentration prediction and the penalty term for the mechanism size.

Figure 3. 

Schematic diagram of the collaborative training framework of Reax-Lump.

Here, $ {\theta }_{i} $ denotes the individual vector, and ξ represents the size penalty factor. A larger ξ poses a stronger penalty on the mechanism size, leading to a more reduced final mechanism. Conversely, a smaller ξ relaxes the size constraint, resulting in a larger mechanism. After a round of combinatorial optimization using genetic algorithms, the kinetic parameters are optimized using the trust-region algorithm. It iteratively solves subproblems to converge toward the optimal solution. The key advantage of this algorithm is that it does not require compution of the gradient of the original objective function, ensuring good stability for this question. Its loss function employs the least squares method, defined as the sum of mean squared errors between the predicted concentrations and the ReaxFF MD simulated concentrations, which is shown in Eq. (4).

It should be noted that both GA and trust region optimization have inherent limitations. When used alone, the ability of GA to reduce the reaction mechanism is constrained by the accuracy requirement, while trust region algorithms struggle to achieve deep convergence with a large mechanism scale. The novel collaborative strategy proposed alternates cyclically between the GA and trust region optimization, which enables simultaneous progress toward both small mechanism size and high predictive accuracy. During the practical optimization of the reaction mechanism model, a maximum of 20 training epochs was adopted. This threshold was selected based on empirical observations that both the kinetics scale and the loss for concentration prediction no longer decreases after 20 epochs, as shown in Supplementary Fig. S2. Additional key parameter values used in the optimization process of this work are available in Supplementary Table S4.

In general, the trade-off between mechanism scale and accuracy in lumped reaction mechanism, optimization was achieved by proposing a collaborative training framework integrating genetic algorithm and trust region optimization. It is a promising and simple method, which can refine thousands of reactions obtained from ReaxFF MD simulation into a reasonable lumped kinetic model.

The Reax-Lump method proposed in this work enables the direct extraction of lumped kinetic mechanisms from ReaxFF MD simulations of polymer plastics. By integrating the SRG-Reax tool, it incorporates detailed reaction classification into the mechanism, which was applied to PP pyrolysis at four heating rates of 1, 2, 4, and 8 K/ps. A detailed analysis of this mechanism is also conducted from the perspective of reaction classification.

Mechanism training and validation for Reax-Lump method

Using the lumping rules in Table 2, each detailed reaction obtained from ReaxFF MD was mapped to its corresponding lumped reaction. Identical lumped reactions are subsequently merged, substantially reducing the total number of unique reactions. Further reduction was achieved by vertical merging on the reaction list along the time-step dimension. As shown in Fig. 4a and Table 3, under the reaction condition of 4 K/ps, the original list of 11,539 reactions were condensed into 801 lumped reactions after lumping and vertical merging steps. Moreover, as shown in Table 4, the initial reaction exhibits a low reaction frequency of approximately 1.0. In contrast, the reaction frequency increased significantly after the application of the lumping procedure, which improved the reliability of the data.

Figure 4. 

(a) Number of reactions after different steps for PP heating at 4 K/ps. (b) The optimized loss function and mechanism scale curve for PP heating at 4 K/ps.

Figure 4b presents the variation trends of the loss function and the reaction mechanism size during the training process, which intuitively reflects the entire training progress. During the training process, both the scale of the reaction mechanism and the loss function for concentration prediction progressively reduced. In each training cycle, the mechanism was first optimized by the genetic algorithm and then further refined using the trust region algorithm. The final lumped results of PP pyrolysis obtained at a heating rate of 4 K/ps are presented in Table 5, while the lumped mechanisms trained at other heating rates are provided in Supplementary Tables S5−S7.

Figure 5a−j compares the number evolution of the 10 lumped species during PP pyrolysis as predicted by the developed lumped kinetic model against ReaxFF MD simulation results at heating rates of 1, 2, 4, and 8 K/ps. The embedded heat-map shown in Fig. 5k indicates the coefficient of determination (R²) between the developed model and the ReaxFF MD data for each simulation condition. Only the R² value of the [C₄₁~C₂₀₀] species at 2 K/ps falls below 0.7, which can be improved by increasing the sampling rate of the lumped species, enlarging the initial population size, or optimizing the lumping rules. All other R² values are above 0.85, indicating significant agreement between the lumped kinetic model and the original ReaxFF MD simulations across most species and heating rates. The lumped kinetic model successfully reproduces the overall PP pyrolysis process obtained from the simulations and captures the distribution of major pyrolyzes. By incorporating target products into the species lumping rules, this approach also enables the development of focused lumped kinetic models tailored to desired products.

Figure 5. 

(a)−(j) Comparison of the predicted number of fragments between ReaxFF MD and lumped mechanisms under 1, 2, 4, and 8 K/ps heating rates. (k) heatmap of R² for fitting between predicted and experimental concentrations of various species at different heating rates.

Furthermore, the lumped kinetic mechanism was employed to predict the yields of small-molecular products at various temperatures under different heating rates. As shown in Fig. 6, C₁, C₂, and C₃~C₅ fragments are methane, ethane plus ethylene, and alkanes/alkenes containing 3–5 carbon atoms, respectively. The molecular weight of the C₃~C₅ fraction is taken as its average value. The yield of C₃~C₅ species dominates the gaseous products, consistent with the findings of Das & Tiwari^[45]. The gas yield increases continuously with temperature, indicating that decomposition reactions are promoted at high temperature to enhance formation of small molecules. At a given temperature, the total gas yield under low heating rates is notably higher than that under high heating rates, with the C₃~C₅ fragments showing the most significant difference. Meanwhile, the yields of total gases and the C₃~C₅ gases under low heating rates begin to decline when T ≥ 2,600 K. This reversal is attributed to the secondary cracking of C₃~C₅ species into smaller C₁ and C₂ fragments, which simultaneously promote their yields. This phenomenon also couples with partial rearrangement reactions to produce higher-carbon species of C₆~C₁₀ and C₁₁~C₂₀.

Figure 6. 

The yields of small-molecule gases under different temperatures and heating rates.

Investigating PP cracking behavior through reaction classification

In addition to obtaining lumped reaction kinetics shown in Table 5 and Supplementary Tables S5−S7, more importantly, when combined with machine learning-driven reaction classification method of SRG-Reax, the specific reaction classes involved in each category of lumped reactions can be identified. For example, Fig. 7 presents the distribution of reaction classes in the PP pyrolysis system under various heating rates. The dominant reaction classes in PP pyrolysis include C–C bond homolysis (RxC1), β-scission (RxC2), intramolecular chain isomerization (RxC5), recombination of C-centered radicals (RxC12), intramolecular H-shift (RxC4), intermolecular H-abstraction by carbon radicals (RxC9), chain cyclization (RxC31), and β-ring opening of branched carbon radicals (RxC23). Obviously, the cleavage of PP and its intermediates into smaller fragments is governed primarily by C–C bond homolysis and β-scission, with the latter being overwhelmingly dominant. The proportion of C–C bond breaking reactions increases with the heating rate, which is attributed to the relatively long reaction times at low heating rates, leading to small-molecule fragments to accumulate and initiates numerous recombinations of C-centered radicals. Among H-transfer-related reactions, the intermolecular H-abstraction by carbon radicals (RxC9) reaction occurs most frequently, which facilitates the unsaturated environment for subsequent C-chain elongation reactions. Correspondingly, at low heating rates, the relatively long reaction residence time leads to the generation of many short-chain carbon-centered radicals, which promotes the occurrence of carbon-centered radical recombination (RxC12) and chain cyclization (RxC9) reactions. Consequently, these types of reactions all exhibit an increasing trends with decreasing heating rates.

Figure 7. 

The distribution of the total 29 reaction classes obtained from ReaxFF MD simulation by using SRG-Reax at different heating rates.

As shown in Fig. 5f, g, the species [C₆~C₁₀] and [C₁₁~C₂₀] play a dominant role and exhibit an increasing trend in the late PP pyrolysis stage at low heating rates of 1 and 2 K/ps. By integrating the SRG-Reax reaction classification, the analysis of the lumped reaction types derived from fitting under different heating rates was performed. Table 6 shows the reaction classes involved in the lumped reaction of [C₃~C₅] + C₂H₆ → [C₆~C₁₀] under different heating rates, with the numbers indicating the count of detailed reactions in each category. It can be observed that at a low heating rate, the total number of lumped reactions from [C₃~C₅] to [C₆~C₁₀] is significantly higher than that at high heating rates, particularly for the recombination of C-centered radicals (RxC12). This explains the reason that the yield of [C₆~C₁₀] does not decrease markedly in late stages at 1 and 2 K/ps. The number of cyclization reactions (RxC31) at heating rates of 1, 2, and 4 K/ps are 15, 2, and 4, respectively, which indicates that cyclization is more favorable under slow heating conditions.

Additionally, the increased radical pool at a low heating rate also promotes H-transfer related reactions of RxC4, RxC9, and RxC11. In terms of reaction proportion, as the rate decreases from 4 to 2 K/ps, the proportion of RxC12 reaction increases with the rise in C-radical concentration. When the rate further drops to 1 K/ps, the concentrations of both C- radicals and H-related active intermediates continue to increase, resulting in high proportions of cyclization reactions and H-transfer related reactions compared to 2 K/ps, which confirms that cyclization is more favorable under slow heating conditions, and the cyclization reaction is accompanied by a large number of H-transfer reactions.

As shown in Fig. 5g, the yield of [C₁₁~C₂₀] species shows a significant increasing trend during the late pyrolysis stage (T ≥ 2,400 K) at a heating rate of 1 K/ps. This phenomenon was analyzed from the perspective of the variation of reaction categories with the stages of the reaction process by using the SRG-Reax method. Most of [C₁₁~C₂₀] species are olefins and alkanes in light oils, and a small number of cycloalkanes and aromatic hydrocarbons appear at high temperatures, which is provided in Supplementary Fig. S3.

In addition, the dynamic profiles of a specific species during PP pyrolysis are analyzed by examining the evolution of the proportion of reaction classes with temperature. Figure 8a shows the number evolution of C₂H₄ with temperature under different heating rates, and Fig. 8b shows the ratio of β-scission (RxC 2) reactions to C-C bond homolysis reactions across temperatures and heating rates. It is observed that the initial formation temperature of C₂H₄ decreases with decreasing heating rate, which aligns well with the temperature range in which the ratio of β-scission to C–C bond homolysis reactions increases significantly in Fig. 8b. Moreover, Fig. 8a indicates that low heating rates promotes C₂H₄ formation, which is also corroborated by the phenomenon that a lower heating rate promotes a higher ratio of β-scission to C–C bond homolysis. It is inferred that a lower heating rate prolongs the reaction residence time per unit time at a given temperature, and results in more complete cracking.

Figure 8. 

Evolution trends of (a) C₂H₄ yields, and (b) RxC₂/RxC₁ ratios with temperature.

Figure 9 displays the proportion of reaction types among all reactions related to [C₁₁~C₂₀] at different temperature stages within 1 K/ps. As shown in Fig. 9a, the initial stage of 1,750–2,100 K is dominated by C-C bond homolysis (RxC1) and β-scission (RxC2). These two reaction classes cause PP and its primary cracking products ([C₂₀₀₊] fragments) to decompose into low-molecular-weight fragments of [C₄₁~C₂₀₀], and [C₂₁~C₄₀], that initially accumulate, and then decrease as they continue to undergo RxC1 and RxC2 reactions. Meanwhile, the formation of [C₁₁~C₂₀] fragments and a small amount of [C₁~C₅] components are observed, leading to the increasing trend of [C₁₁~C₂₀] species.

Figure 9. 

Comparison of the proportion of reaction categories related to [C₁₁~C₂₀] at various temperature stages under the heating rate of 1 K/ps. (a) 1,750–2,100 K; (b) 2,100–2,500 K; (c) 2,500–3,000 K.

Accordingly, during the middle stages of 2,100–2,500 K (Fig. 9b), cracking reactions of RxC1 and RxC2 remain dominant by reactions of. Particularly, a notable increase in hydrogen abstraction reactions (RxC9) raises the concentration of carbon-centered radicals, favoring a higher proportion of RxC2 and a lower proportion of RxC1. Correspondingly, the yields of [C₁₁~C₂₀] and [C₃~C₅] components decrease with temperature during this period (Fig. 5e, g). In the late stage (Fig. 9c), the proportion of RxC12 reactions increase significantly alongside a greater diversity of reaction types. Moreover, the increased proportions of hydrogen addition, abstraction, and transfer reactions (RxC9, RxC4, RxC5, RxC11) as well as cyclization reactions (RxC31) contribute to the conversion of low-carbon fragments to [C₁₁~C₂₀]. The increase in the proportion of cyclization reactions can be corroborated by the rise in the number of rings in the late stage (Supplementary Fig. S3). Therefore, in summary, the enhancement of recombination and cyclization reactions collectively contribute to the increase in [C₁₁~C₂₀] yield in the late PP pyrolysis stage. Therefore, the reaction lumping method, when integrated with the reaction classification strategy, not only can lump thousands of reactions and species into a manageable kinetic model suitable for downstream reactor simulations, but also reveals the specific reaction types encapsulated within each lumped reaction. This provides critical support for identifying the dominant reaction mechanisms responsible for the formation of different desired products.

This work proposes the Reax-Lump method, an automated approach for the extraction of lumped reaction mechanisms for polymer plastic pyrolysis from large-scale ReaxFF molecular dynamics simulation results. To address the critical challenge of balancing mechanism reduction and predictive accuracy, a collaborative training optimization scheme was developed based on genetic algorithm and trust-region method. This strategy iteratively refines the mechanisms iteratively refines the mechanism by alternating between scale reduction and kinetic parameters optimization, significantly reducing the dimensionality of the reaction mechanism while ensuring its prediction accuracy for polymer pyrolysis. In addition, by building a labeled database for PP pyrolysis reaction categories, this work extends the applicability of the SRG-Reax method to the PP systems, and successfully realized the combination of the Reax-Lump method and the SRG-Reax method.

The effectiveness of Reax-Lump was evaluated by comparing the predicted results of 10 lumped species against the original ReaxFF MD simulation results across four heating rates. The average coefficient of determination (R²) for different species of 0.956 demonstrates the reliability of the lumped mechanism extraction scheme proposed. Moreover, the classification accuracy of the SRG-Reax for PP pyrolysis increases from its original 51% to 97.7%, indicating its promising role in classifying complex reactions. By integrating species lumping, reaction classification, and kinetic collaborative optimization, the lumped kinetic models for PP pyrolysis at 1, 2, 4, and 8 K/ps were constructed. These models comprising 10 species and approximately 20 lumped reactions are chemically interpretable, with each lumped reaction assigned to its specific reaction classes, which explains the evolving trends of [C₃~C₅] and [C₁₁~C₂₀] yields.

The Reax-Lump method provides an automated, integrated workflow for constructing interpretable and compact pyrolysis mechanisms from ReaxFF MD data. Although demonstrated here for PP pyrolysis, the method is readily extendable to other single- or multi-component polymer plastic systems. Importantly, once calibrated with experimental data, the pyrolysis mechanisms generated by Reax-Lump are promising for use in downstream reactor modeling, which is promising for practical applications in predicting product distributions for polymer plastic pyrolysis.

It is notewothy that although validated for successful olefin polymer systems, the current method cannot be directly extended to systems containing heteroatoms beyond C, H, and O, as SRG-Reax lacks the capability to accurately classify the corresponding reaction classes. Furthermore, due to inherent differences between ReaxFF MD simulations and experiments, the lumped kinetic mechanism obtained here is not directly applicable to numerical simulations. The current lumping strategy also relies solely on carbon and hydrogen number, without incorporating species reactivity or industrial requirements. Future work will focus on redefining SRG-Reax reaction labels and constructing dedicated training datasets to extend classification capability to more complex polymer systems. Prior to numerical simulation integration, the lumped mechanisms will be calibrated against experimental data or literature. A universal template library will also be developed to enable flexible, targeted species lumping based on industrial needs or target species of interest.