Ready-to-Use Polymerization Simulations Combining Universal Machine Learning Interatomic Potential with Time-Dependent Bond Boosting for Polymer and Interface Design

Abstract

Although polymerization and curing reactions govern the performance of advanced materials, their simulation remains challenging owing to the need for accurate, transferable potentials and rarity of chemical events. Conventional reactive force fields such as ReaxFF require system-specific parametrization, while universal machine learning interatomic potentials (uMLIPs) exhibit limited sampling efficiency. This paper introduces a novel simulation framework integrating a uMLIP with a time-dependent bond-boost scheme. The bias potential increases monotonically with time, and the use of a unified parameter set across reaction classes enables consistent acceleration without system-specific tuning. For radical polymerization of vinyl monomers, the proposed framework reproduces characteristic trends, such as linear molecular-weight growth with conversion, initiator-concentration scaling, and relative monomer reactivity trends. For step-growth polycondensation of nylon-6,6, it captures the characteristic sharp increase in molecular weight at high conversion rates, consistent with experimental behavior. For epoxy curing at a copper substrate, it reveals interfacial ring-opening and cross-linking events, consistent with spectroscopic evidence of Cu—O–C bond formation. Overall, coupling uMLIPs with time-dependent bond boost enables practical and transferable simulations of polymerization and curing processes. The proposed framework reliably resolves mechanistic pathways and relative reactivity, offering molecular-level insights into polymer growth and interfacial adhesion.

Preferred Networks]Preferred Networks, Inc., Tokyo 100-0004, Japan Preferred Networks]Preferred Networks, Inc., Tokyo 100-0004, Japan Preferred Networks]Preferred Networks, Inc., Tokyo 100-0004, Japan Preferred Networks]Preferred Networks, Inc., Tokyo 100-0004, Japan

1 Introduction

Polymeric materials play an indispensable role across a wide range of industrial applications, including electronics, medical devices, energy storage, and structural components^{1, 2, 3, 4}. Their performance strongly depends on the molecular structures formed by microscopic chemical processes such as polymerization and cross-linking, highlighting the importance of understanding structural changes and reaction mechanisms at the molecular level^{5, 6}. Thus, to design high-performance materials and develop novel polymeric systems, the corresponding reaction pathways, intermediates, and transition states must be comprehensively explored.

Experimental techniques such as time-resolved spectroscopy and electron microscopy have been used to analyze polymer reactions^{7, 8}. However, direct observation of reaction intermediates and detailed mechanisms at the atomic level remains challenging^{9, 10, 11, 12}. Molecular simulations serve as powerful complementary tools that allow the identification of reaction intermediates and transition states that are difficult to capture experimentally^{13, 14, 15, 16}. By providing atomistic insights, these simulations offer valuable guidelines for the rational design of advanced materials.

Two methods are widely used for simulating cross-linking and polymerization reactions: pseudoreaction models and reactive potentials. Pseudoreaction models stochastically rearrange chemical bonds when specified distance or angular criteria are met^{17, 18, 19, 20, 21}. Although these models can be combined with classical molecular dynamics (MD) to handle large systems, they necessitate complex parameter tuning, including the definition of activation energies and frequency factors in advance. Quantum chemical approaches such as global reaction route mapping can accurately calculate activation energies. However, their computational cost remains prohibitive for large or structurally diverse condensed-phase systems^{5, 22, 23}.

A well-known reactive force field approach is ReaxFF, which describes the formation and breaking of chemical bonds during MD simulations²⁴. However, many chemical reactions in polymerization processes are rare events, meaning that straightforward nanosecond-scale MD simulations often fail to capture sufficient reaction progress. This problem becomes particularly pronounced in polymerization, which involves a large number of sequential reactions, necessitating unrealistically long simulation times²⁵. To address this, acceleration techniques such as bond-boost MD have been developed^{26, 27}. However, the creation and parametrization of specialized force fields for each target system require substantial expertise and effort. Moreover, the bond-boost method involves hyperparameters that, if set improperly, may compromise the reliability of the observed reactions. Another widely used approach for enhancing the sampling of rare events is metadynamics, which applies history-dependent bias potentials to efficiently explore free-energy landscapes ^{28, 29}. However, by design, metadynamics promotes both forward and backward reactions, making it less suitable for directly accelerating one-way polymerization processes.

Given these limitations, there exists a strong demand for simulation methods that are simple, accurate, and broadly applicable to diverse polymerization processes. Recently, universal machine learning interatomic potentials (uMLIPs) trained on density functional theory data have attracted increasing attention in this domain.^{30, 31, 32} In particular, uMLIPs such as PreFerred Potential (PFP), M3GNet, CHGNet, and MACE have demonstrated remarkable generalization across various atomic species and structures,^{33, 34, 35, 36, 37} enabling their application to catalysis, metal–organic frameworks, and solid–liquid interfaces.^{38, 39, 40, 41} The incorporation of these uMLIPs largely eliminates the need to construct system-specific force fields, thereby broadening the scope of atomistic simulations.

Considering these aspects, this paper introduces a time-dependent reaction acceleration strategy that combines the broad transferability of a uMLIP with efficient sampling. This approach enables polymerization simulations that are both accurate and straightforward to implement, even for chemically complex systems. To validate the method, we first benchmark the radical polymerization of vinyl monomers, focusing on kinetics and relative reactivity. We then demonstrate its broad applicability through step-growth polycondensation of nylon-6,6 and the curing of an epoxy resin at a copper oxide (CuO) interface, a system known to be particularly challenging to treat using conventional reactive force fields such as ReaxFF. Collectively, these applications highlight the potential of the proposed framework to provide accurate and convenient polymerization simulations across diverse chemical processes.

2 Computational Details

2.1 Time-dependent Bond-Boost Potential

To accelerate reactions, we use a modified version of the bond-boost method originally proposed by Vashisth et al. ²⁶. The original bond-boost potential is expressed as

$$V(r)=f_{1}\Bigl\{1-\exp\bigl[-f_{2}(r-r_{0})^{2}\bigr]\Bigr\}.$$ (1)

Here, $f_{1}$, $f_{2}$, $r_{0}$, and $r$ represent the depth of the boost potential, its range, the target distance, and the interatomic distance, respectively. Appropriate values for $f_{1}$, $f_{2}$, and $r_{0}$ must be specified prior to the simulation. Vashisth et al. ²⁶ and Dasgupta et al. ²⁷ examined multiple parameter sets and adopted the acceleration conditions that best reproduced the activation energies obtained using ReaxFF. Typically, determining $f_{1}$ requires estimating the activation energy in advance and then fine-tuning the parameters through trial and error. Furthermore, for bond cleavage, an appropriate distance $r_{0}$ must be defined as the threshold for bond breaking. As the number of target reactions increases, such preliminary parameter studies become increasingly cumbersome, complicating the overall simulation setup.

To address this, in this work, we introduce a time-dependent bond-boost (TDBB) method as a modification of the original approach. The key innovation is that the bias potential increases with time, as shown in Eqs. 2, 3, 4 and 5. Here, $V^{f}(r,t)$ and $V^{d}(r,t)$ represent the bias potentials applied to promote bond formation and bond dissociation, respectively:

$$V^{f}(r,t)=f_{1}(t)\Bigl\{1-\exp\bigl[-f_{2}(r-r_{0})^{2}\bigr]\Bigr\},$$ (2)

$$V^{d}(r,t)=f_{1}(t)\exp\bigl(-f_{2}r^{2}\bigr),$$ (3)

$$r_{0}=\lambda\sum_{a\in\{\text{pair}\}}r_{a}^{\mathrm{vdw}},$$ (4)

$$f_{1}(t)=\begin{cases}\gamma t,&f_{1}(t)<f_{1}^{\mathrm{max}},\\ f_{1}^{\mathrm{max}},&f_{1}(t)\geq f_{1}^{\mathrm{max}},\end{cases}$$ (5)

where $\gamma$ is the boost factor, $r_{a}^{\mathrm{vdw}}$ is the van der Waals radius of atom $a$ (summed over the two atoms forming the reactive pair), $\lambda$ is a scale factor defining the bonding threshold, and $f_{1}^{\mathrm{max}}$ is the maximum boost amplitude. Because $f_{1}(t)$ increases with time, reactions are gradually promoted as the simulation proceeds. Moreover, applying a single set of boost parameters to all possible reactions ensures consistent acceleration and the capturing of relative reactivities without system-specific parameterization. The overall acceleration level can be tuned by $\gamma$, with larger values yielding higher reaction rates.

Notably, we employ distinct functional forms for bond formation, $V^{f}(r_{ij})$, and bond breaking, $V^{d}(r_{ij})$, to ensure that once the product or dissociated state is reached, the bias potential vanishes, and normal MD resumes. In the bond-formation case, the potential drives two atoms to approximately 60% of the sum of their van der Waals radii. Once that distance is reached, the potential effectively becomes zero. In the bond-breaking case, once the bond extends beyond a critical length, the boost potential similarly disappears. This design eliminates the need to define separate parameters for each reaction, thereby simplifying the simulation setup.

According to Eqs. 2, 3, 4 and 5, the parameters $f_{1}$, $f_{2}$, and $r_{0}$ correspond closely to those in the original bond-boost method. Although these parameters can be fine-tuned if desired, $f_{1}^{\mathrm{max}}$ primarily serves to cap the boost strength and prevent excessively large acceleration. Similarly, $f_{2}$ controls how the boost varies with distance: Smaller values concentrate the effect at very short range (strong near contact), whereas larger values extend the influence over longer distances. By setting $r_{0,\mathrm{create}}$ to the sum of the van der Waals radii of the involved atoms, this parameter becomes independent of the reaction pathway. In practice, the same settings have been successfully applied to multiple reactions.

The basic functional form follows the original bond-boost design, ensuring that once a reaction occurs, the boost potential decays to zero and the system behaves as in a standard MD simulation without acceleration. An additional advantage of this straightforward Morse-like potential is that it consistently exerts force along the direction of the intended reaction coordinate. In contrast, adaptive bias methods such as metadynamics can, in some cases, drive the system in the opposite direction of the desired pathway. Although adding a high-wall restraint can prevent the bias from diverging, such methods typically fill energy minima before effectively accelerating the system toward the reaction coordinate, potentially delaying the onset of meaningful acceleration. Moreover, in metadynamics, the bias can continue to accumulate even after a reaction has occurred, necessitating additional steps (e.g., detecting the occurrence of the reaction) to halt bias deposition.

Lastly, the proposed method does not allow time rescaling, unlike conventional hyperdynamics ^{42, 43, 44}. According to transition-state theory, time rescaling in biased MD is valid only if (1) the added potential does not exceed the activation energy and (2) no potential is applied at the dividing surface. In contrast, in the proposed approach, the bias grows monotonically with time and can surpass the activation barrier regardless of the potential energy landscape. Thus, these conditions are generally not satisfied. Notably, while this method simplifies parameter setup and effectively drives chemical reactions, the time evolution of the trajectory no longer strictly corresponds to physical time once the bias becomes substantial.

2.2 Polymerization Procedure

The proposed simulation method involves four main steps (Figure 1): (1) definition of reactive atomic groups, (2) selection and acceleration of reaction-targeted groups, (3) reaction induction and structural relaxation, and (4) evaluation of reactions and updating of atomic groups.

Figure 1: Schematic of the proposed simulation workflow, including four main stages: (1) definition of reactive atomic groups, where potential reactive species are identified from chemical functionalities; (2) selection and acceleration of reaction-targeted groups, where candidate atomic combinations satisfying geometric and chemical criteria are chosen and subjected to time-dependent bias potentials; (3) reaction induction and structural relaxation, in which accelerated molecular dynamics simulations promote bond rearrangements followed by unbiased relaxation; and (4) reaction evaluation and group update, where the bonding states are analyzed and reactive groups are updated for subsequent iterations.

1. Definition of Reactive Atomic Groups: In the initial stage, reactive atomic groups are identified and defined as sets $G_{I}$, $G_{J}$, $G_{K}$, and $G_{L}$, each corresponding to a specific group (e.g., oxygen atoms in epoxy groups, nitrogen atoms in amine groups, or $\alpha$-carbon atoms in vinyl groups for radical polymerization). These sets are defined as

$$G_{X}=\{a\mid a\in X\},\quad X=I,J,K,L,$$ (6)

with $a$ representing an individual atom belonging to the group. Chemical reactions are modeled as rearrangements of atomic bonds, represented by bond formation and breakage.

2. Selection and Acceleration of Reaction-Targeted Groups: To efficiently identify possible reaction candidates, a two-stage selection process is employed. First, neighboring atomic groups are preselected based on geometric proximity to identify potential reacting partners. Subsequently, based on reaction rules, candidate groups are formed by selecting tuples $(i,j,k,l)$ such that $i\in G_{I}$, $j\in G_{J}$, $k\in G_{K}$, and $l\in G_{L}$, with interatomic distances constrained to the window $r_{ab}^{\min}\leq r_{ab}\leq r_{ab}^{\max}$. Defining the set of target pairs as $\mathcal{P}=\{(i,j),(i,k),(j,l)\}$, the candidate set is expressed as

$$G_{ijkl}=\{(i,j,k,l)\mid \,i\in G_{I},\,j\in G_{J},$$ (7)

$$\,k\in G_{K},\,l\in G_{L},$$

$$\,r_{ab}\in[r_{ab}^{\min},r_{ab}^{\max}]$$

$$\,\forall(a,b)\in\mathcal{P}\}.$$

Moreover, we allow multiple reaction types to coexist. For each reaction type $r$, a candidate set $G_{ijkl}^{(r)}$ is prepared with its own group assignments $(G_{I},G_{J},G_{K},G_{L})$, target-pair set $\mathcal{P}^{(r)}$, and distance windows $\bigl[r_{ab}^{\min},r_{ab}^{\max}\bigr]^{(r)}$. This setup enables systems in which different reactions proceed simultaneously within a single simulation box (e.g., Dehydration–condensation and radical polymerization occurring in parallel). Unless otherwise noted, candidates from all reaction types are pooled and subjected to the same nonoverlapping selection and biasing procedure described below.

As candidate groups may share atoms (e.g., $(0,1,2,3)$ and $(0,4,5,6)$), a nonoverlapping selection is performed by computing the score $d_{ijkl}=r_{ij}+r_{ik}+r_{jl}$ for each group, sorting the obtained scores in ascending order, and sequentially selecting groups while skipping those with overlapping atoms. This yields the refined set $G_{ijkl}^{\mathrm{acc}}$. For each group in $G_{ijkl}^{\mathrm{acc}}$, the set of reactive atomic pairs $P$ is then defined, consisting of one or more pairs selected from $(i,j)$, $(i,k)$, $(j,l)$, and $(k,l)$ that participate in bond formation or cleavage. A time-dependent bias potential is subsequently applied to each pair in $P$:

$$\Delta V(t) =\sum_{(i,j,k,l)\in G_{ijkl}^{\mathrm{acc}}}\sum_{p\in P}\Bigl[$$ (8)

$$\quad f_{p}\,V^{f}(r_{p},t)$$

$$\quad+\bigl(1-f_{p}\bigr)\,V^{d}(r_{p},t)\Bigr],$$

where $f_{p}=1$ denotes bond formation, $f_{p}=0$ denotes bond breakage, and $r_{p}$ is the interatomic distance of pair $p$.

3. Reaction Induction and Structural Relaxation: The system is subjected to an accelerated MD simulation under the applied bias $\Delta V(t)$ for either a predetermined period or until a reaction event is observed. A reaction event is defined to occur when the interatomic distances of a specific set of pairs ($ij$, $ik$, and $jl$) simultaneously satisfy the prescribed bonding conditions. For each pair, the bonding state is determined based on the van der Waals criterion: Two atoms are considered bonded if their separation falls below 60% of the sum of their van der Waals radii. If $f_{p}=1$, the event is triggered when these pairs newly satisfy the bonding condition (bond formation). If $f_{p}=0$, the event is triggered when the pairs cease to satisfy the condition (bond dissociation). Once such an event is detected, or the maximum simulation time is reached, the acceleration phase terminates. Subsequently, a standard MD simulation (without the bias) is performed to allow structural relaxation and equilibration of the resulting configurations.

4. Reaction Evaluation and Group Update: After structural relaxation, the system is analyzed to compare bonding states before and after the biased simulation. Atoms that have participated in new bond formations or bond cleavages are marked as “reacted” and removed from the reactive set. In contrast, atoms that experienced the bias but did not undergo any bond changes remain in the reactive set for subsequent iterations. If the resulting species remain reactive (e.g., in epoxy–amine curing or vinyl radical polymerization), these atoms are reassigned to Group $X$ and considered in subsequent iterations. Repeating these steps iteratively drives the overall polymerization process.

2.3 Comparison with existing acceleration methods

To highlight the novelty and advantages of the proposed framework, Table 1 summarizes the key features of widely used rare-event acceleration techniques. Conventional bond boost requires system-specific parameter tuning based on activation barriers, whereas hyperdynamics and metadynamics depend on carefully chosen collective variables and typically accelerate both forward and backward events indiscriminately. In particular, even infrequent metadynamics—the only variant that, in principle, allows approximate recovery of physical timescales—is intrinsically bidirectional: Gaussian hills are deposited regardless of whether the system proceeds forward or backward, which makes the method inefficient and physically inconsistent for inherently unidirectional processes such as radical polymerization. In contrast, bond-boost schemes apply a forward-only bias that selectively promotes bond formation, aligning with the irreversible nature of polymer growth. Our TDBB variant retains this property while eliminating the need for prior parametrization, allowing a single parameter set to be used across reaction types. Combined with uMLIP such as PFP, this approach offers a unique balance of simplicity, accuracy, and transferability, and is particularly well suited for chemically complex and heterogeneous polymerization systems.

Table 1: Comparison of acceleration methods, focusing on bias form and setup requirements.

Method	Setup requirements	Bias function	Reaction directionality	Time rescaling
Bond Boost 45	Reactive groups, bias height, target distance	Fixed bias	Unidirectional	No
Metadynamics	Collective variables (CV)	Adaptive hill bias	Bidirectional	No ‡
CVHD	CV and dividing surface	Adaptive hill bias	Bidirectional	Yes*
TDBB (this work)	Reactive groups	Monotonic, time-growing bias	Unidirectional	No

^* Hyperdynamics time rescaling is possible if conditions are met (zero bias at dividing surface, bias $<$ barrier, TST assumptions).
^‡  Except infrequent metadynamics.

3 Methods

We investigated three representative systems: radical polymerization of vinyl monomers, step-growth polycondensation of nylon-6,6, and epoxy curing at a CuO interface. All systems were prepared through a two-stage equilibration procedure: initial relaxation with a classical force field, followed by refinement using the PFP uMLIP under the NPT ensemble. Detailed compositions and equilibration protocols are provided in the Supporting Information.

All simulations were conducted using PFP (version v6.0.0) in the CRYSTAL_PLUS_D3 mode. PFP is a uMLIP trained on Perdew—Burke–Ernzerhof reference data and augmented with the D3 dispersion correction, which improves the reproducibility of condensed-phase properties such as density. Owing to its transferability across 96 elements, PFP has been applied to a wide range of materials, including inorganic and organic compounds, polymers, and catalytic systems.

Reactive acceleration MD was performed using the boost potential described in Section 2.1. Unless otherwise noted, simulations employed a time step of 0.25 fs, alternating biased and unbiased dynamics every 2000 steps (500 fs), and bias parameters of $f_{1}^{\mathrm{max}}=250$ kcal/mol for $V^{f}$, $f_{1}^{\mathrm{max}}=125$ kcal/mol for $V^{d}$, $f_{2}=10~\text{\AA }^{-2}$, and $\gamma=1.0$. Reaction events were identified based on van der Waals radii as detailed in Section 2.1. For each condition, three independent simulations with randomized initial configurations were performed, and the reported results represent averages over these runs. System-specific simulation parameters, including system size, composition, and equilibration protocol, are provided in the Supporting Information.

To evaluate the robustness of the acceleration scheme, we examined the sensitivity of the cumulative number of reactions to variations in $f_{1}^{\mathrm{max}}$, $f_{2}$, and $\gamma$, using styrene radical polymerization as a representative system (Fig. S4). The influence of $f_{1}^{\mathrm{max}}$ was negligible once it exceeded a moderate threshold: extremely low values suppressed reactions within the time window, whereas sufficiently large values simply shifted the reaction to occur later within the same window, after which the reactive-pair list was updated. System behavior was also insensitive to variations in $f_{2}$. Although $f_{2}$ formally defines the spatial range of the bias, propagation is governed by near-contact events, and $f_{2}\approx 10$ was sufficient to ensure consistent behavior. In contrast, increasing $\gamma$ accelerated all reaction channels uniformly, acting as a global rate scaling factor. These analyses demonstrate that the default parameter set is broadly representative and does not distort relative reactivity trends.

All MD simulations were carried out using OpenMM ver. 8.1.1 ⁴⁶. Forces and energies were evaluated using PFP, with GPU selection depending on the system size: an NVIDIA Tesla V100 (16 GB memory) was used for medium-sized systems, whereas larger systems were computed on an NVIDIA A100 (80 GB memory). A custom interface coupled OpenMM with PFP. The bias potential was implemented via OpenMM–Torch, enabling TorchScript-defined energy and force terms to be directly incorporated. Force derivatives were evaluated using PyTorch’s automatic differentiation, eliminating the need for manual implementation. The workflow of reactive acceleration MD was managed using an in-house Python script.

4 Results and Discussion

4.1 Radical polymerization of vinyl monomers

The simulations were designed with two primary objectives: (i) reproduce the kinetics of radical polymerization under solution conditions, and (ii) reproduce the relative reactivity order of representative vinyl monomers under bulk conditions. The monomer set included methyl acrylate, methyl methacrylate, styrene, vinyl acetate, diphenylethylene, and dimethyl itaconate. Polymerizations were performed both in toluene solution and in the bulk phase. Toluene was selected as a representative nonpolar solvent commonly used in experimental radical polymerizations, whereas bulk simulations enabled direct assessment of solvent-free reactivity. Together, these two environments and the diverse monomer set provided a stringent benchmark for evaluating the predictive capability of the proposed simulation framework.

Figure 2 illustrates three key aspects of the styrene radical polymerization under varying concentrations: (i) time–conversion profiles, (ii) relationship between monomer conversion and number-average molecular weight ($M_{n}$), and (iii) dependence of the polymerization rate on the initiator-to-monomer ratio (I/M). As shown in Figure 2 (b), systems with I/M ratios of 5/100 and 10/100 exhibit markedly higher conversions at earlier times compared to those with lower initiator loadings. This trend demonstrates that increasing the initiator and/or monomer concentrations accelerates polymerization, consistent with experimental observations and classical kinetic predictions. These results highlight the ability of the proposed uMLIP–accelerated MD framework to accurately capture the early-stage kinetics of radical polymerization.

The plots of conversion versus $M_{n}$ in Figure 2 (c) exhibit a clear linear relationship across all simulated conditions. This linear trend strongly suggests that, in the absence of termination and chain-transfer reactions (see Figure. S1), polymerization proceeds in a “living” manner, characterized by continuous chain growth without significant loss of active chain ends. Moreover, the slope of the conversion–molecular-weight curve for the I/M = 10/100 system is approximately half of that for the lower-initiator systems. This behavior is characteristic of living polymerization, where the number-average degree of polymerization is inversely proportional to the initiator concentration. This finding corroborates the fidelity of the proposed simulations in reproducing mechanistic features observed experimentally in controlled polymerization systems.

Figure 2 (d) illustrates the dependence of the polymerization rate on the I/M ratio. In the absence of termination reactions, the polymerization rate $R_{p}$ is expressed as

$$R_{p}=-\frac{d[M]}{dt}=k_{p}[M][P^{*}]$$ (9)

$$R_{p}\propto[M][I],$$ (10)

where $[M]$ and $[I]$ denote the monomer and initiator concentrations, respectively; $[P^{*}]$ is the concentration of propagating radicals; and $k_{p}$ is the propagation rate constant. In this idealized case, all radicals generated from the initiator remain active, resulting in a linear dependence of the polymerization rate on initiator concentration. This is unlike the square-root dependence observed in conventional radical polymerizations.

The simulation results closely follow this theoretical proportionality, confirming that the proposed approach reproduces the linear dependence of the polymerization rate on initiator concentration and its proportionality to monomer concentration under the no-termination assumption. This agreement between simulation and theory not only validates the physical accuracy of the uMLIP but also underscores the robustness of the acceleration scheme in capturing fundamental kinetic behavior across a range of concentrations. The linear dependence of $R_{p}$ on initiator concentration is a hallmark of living radical polymerization, where active chain ends persist without termination.

Figure 2: Structural and kinetic characteristics of styrene radical polymerization in toluene solvent. (a) Snapshot after polymerization and a representative polymer chain. (b) Time–conversion curves at different initiator-to-monomer (I/M) ratios. The vertical axis $\alpha$ represents the monomer conversion, defined as $\alpha=1-[M]/[M]_{0}$, where $[M]_{0}$ and $[M]$ are the initial and instantaneous monomer concentrations, respectively. (c) $M_{\mathrm{n}}$ as a function of conversion. (d) Polymerization rate plotted against $[I][M]$. The horizontal axis represents the product of the initiator and monomer concentrations, $[I][M]$, converted from number densities to molar units (1 $\text{\AA }^{-6}$ = 2.76 × 10^-4 mol² L^-2). Here, I, M, and S denote the initiator (AIBN), monomer (styrene), and solvent (toluene), respectively.

Collectively, these results demonstrate that combining machine-learned potentials with accelerated dynamics can quantitatively reproduce key mechanistic and kinetic trends in radical polymerization. By accurately capturing the dependence of reaction rate and molecular weight evolution on initiator and monomer concentrations, the proposed framework provides a solid foundation for predictive simulations of polymerization and offers valuable insights into the molecular-level dynamics governing polymer growth.

Figure 3 presents the apparent propagation rate constants ($k_{p}^{*}$) for various vinyl monomers. The apparent rate constants were obtained by fitting the monomer conversion profiles to a single-exponential function,

$$\alpha(t)=1-\exp(-k_{p}^{\ast}t),$$ (11)

where $\alpha(t)$ is the monomer conversion at time $t$, and $k_{p}^{\ast}$ represents the effective propagation rate under the simulated conditions. Representative fits are shown in Figure 3 (a), together with the extracted $k_{p}$ values for each monomer. These fitted values were then ranked and compared against experimental propagation rate constants (Figure 3 (b)). Rank correlations between the simulated and experimentally reported propagation rates yield a Spearman correlation coefficient ($\rho$) of 0.66, Kendall rank correlation ($\tau$) of 0.60, and mean absolute deviation of 1.0. While the agreement is not perfect, the simulations provide a reasonable basis for ranking relative monomer reactivities.

To probe the origin of discrepancies in ranking, explicit transition-state searches were performed using the PFP to estimate activation energies for representative propagation reactions (Table S4). Although experimental activation energies place styrene in the middle of the monomer series, the barriers computed using PFP rank styrene as the lowest among the studied monomers. This inconsistency in energetic ordering parallels the deviation observed in the simulated propagation rate constants, indicating that the adopted potential systematically underestimates the propagation barrier of styrene. In contrast, bulk densities are accurately reproduced (Table S4), suggesting that the discrepancy originates from the description of reaction barriers rather than from the condensed-phase structural properties. Improving the barrier fidelity of the PFP is therefore expected to directly enhance the reliability of kinetic rankings.

Nonetheless, rank inversions are observed for certain monomers, highlighting areas where the proposed model struggles to capture subtle differences in propagation kinetics. The range of propagation rate constants examined corresponds to activation energy variations of only a few kcal/mol. Accurately resolving such fine energetic differences will likely require higher-fidelity uMLIP trained on more comprehensive and accurate reference data, or advanced acceleration schemes capable of capturing the delicate energetics of reactive events. Notably, because the applied bias alters the effective timescale of the dynamics, the proposed approach is not intended to reproduce absolute propagation rate constants. Rather, its primary utility lies in capturing relative trends and reactivity rankings across monomers without requiring prior knowledge of activation energies. Incorporating higher-fidelity potentials or refined acceleration strategies is expected to further improve quantitative agreement and expand the applicability of this framework to complex polymerization systems.

Figure 3: (a) Time–conversion profiles (blue) and exponential fits (red) for representative vinyl monomers, with extracted apparent propagation rate constants $k_{p}$. (b) Ranking of the apparent $k_{p}$ values compared with experimental data. Rank correlations indicate Kendall’s $\tau=0.60$, Spearman’s $\rho=0.66$, and mean absolute deviation = 1.0.

4.2 Step-Growth Polycondensation of Nylon-6,6

To further evaluate the capability of the proposed framework in modeling step-growth polymerization, we examined the polymerization of nylon-6,6, a prototypical condensation polymer. The results (Figure 4) show the time evolution of monomer conversion and the relationship between conversion and DP_n. Here, DP_n denotes the number-average degree of polymerization, defined as the average number of monomeric units per polymer chain. DP_n is related to the number-average molar mass $M_{\mathrm{n}}$ by $\mathrm{DP}_{n}=M_{\mathrm{n}}/M_{0}$, where $M_{0}$ is the molar mass of the repeating unit. Consistent with classical theories of step-growth polymerization, DP_n increases monotonically with conversion, particularly in the later stages dominated by oligomer–oligomer coupling. This trend, a hallmark feature of step-growth kinetics, is reproduced with high fidelity across all runs, underscoring the robustness of the methodology. The quantitative agreement between simulation and theoretical expectations indicates that the uMLIP–accelerated framework accurately captures the essential thermodynamic and kinetic characteristics of condensation polymerizations.

Our results also align with experimental observations of nylon polycondensation. In melt polymerization without continuous water removal, conversion typically levels off at 60–80%.⁴⁷ Our simulations reproduce this incomplete conversion and limited chain length, in contrast to idealized models that assume perfect removal of the condensate. This correspondence highlights the ability of the proposed framework to capture the equilibrium-limited nature of step-growth polymerizations without requiring predefined activation barriers or frequency factors for each reaction pathway, demonstrating a clear advantage over heuristic methods.

Figure 4: Structural and kinetic characteristics of nylon-6,6 step-growth polymerization under melt conditions without water removal. (a) Representative snapshot after polymerization, showing polymer chains (tan) and generated water molecules (blue). (b) Time evolution of monomer conversion ($\alpha$), which levels off below unity owing to equilibrium limitation. (c) Number-average degree of polymerization (DP_n) as a function of conversion, compared with the theoretical Carothers relation. The simulations reproduce the characteristic nonlinear DP_n–conversion relationship and incomplete conversion observed experimentally in condensation polymerizations without continuous removal of condensate.

4.3 Epoxy Curing at the CuO Interface

Reliable adhesion between epoxy resins and metallic wiring layers is essential for improving the performance and reliability of semiconductor devices. Surface treatments such as oxygen plasma are commonly applied to enhance adhesion, and spectroscopic studies have revealed the presence of Cu–O–C linkages at treated interfaces.⁴⁸ These observations suggest that covalent bond formation across the epoxy–copper interface is a key factor underlying the improved interfacial strength. However, molecular-level understanding of these interfacial curing reactions remains challenging to achieve owing to their chemical complexity and large system sizes.

This section describes the assessment of the ability of the proposed accelerated simulation framework to capture interfacial curing chemistry. The considered case is particularly demanding: No publicly available ReaxFF parametrization simultaneously covers Cu surfaces and epoxy–amine chemistry, and full ab initio dynamics of interfacial systems of this scale are computationally infeasible. Consequently, conventional approaches are not directly applicable, and uMLIPs provide the only practical route to modeling these interfaces. Even so, straightforward uMLIP MD suffers from the rare-event bottleneck, failing to capture curing reactions within accessible time windows.

To analyze the curing reaction, Figure 5(a) shows the relative concentration change $c(t)/c(0)$ of key atomic species as a function of simulation time:

$$c(t)/c(0)=\frac{N_{\mathrm{species}}(t)}{N_{\mathrm{species}}(0)}.$$

Here, $N_{\mathrm{species}}(t)$ is the number of target atoms at time $t$. Both epoxy oxygen ($O_{\mathrm{e}}$) and primary amine nitrogen ($N_{\mathrm{p}}$) decrease over time, while secondary amine nitrogen ($N_{\mathrm{s}}$) increases. This indicates that epoxy ring-opening by primary amines is the dominant reaction in the early stage. The increase in $N_{\mathrm{s}}$ confirms the formation of secondary amines, consistent with the expected initial step of epoxy curing. Experimental studies of diglycidyl ether of bisphenol A/diethylenetriamine systems using near-infrared spectroscopy reveal that primary amines are nearly consumed before secondary amine conversion begins. Our simulations capture this early-stage behavior. The number of hydroxyl hydrogens on the CuO surface ($H_{\mathrm{CuOH}}$) decreases slightly, suggesting limited reaction at the interface.

To further characterize the spatial distribution of reaction events, we evaluated the depth-resolved reaction density $\rho_{\mathrm{rxn}}(z)$:

$$\rho_{\mathrm{rxn}}(z)=\frac{N_{\mathrm{rxn}}(z)}{A\,\Delta z\,N_{\mathrm{frames}}},$$

where $N_{\mathrm{rxn}}(z)$ is the number of reactions observed in bin $z$, $A$ is the cross-sectional area, and $\Delta z$ is the bin width, set to 8 Å. Figure 5(b) shows $\rho_{\mathrm{rxn}}(z)$ for epoxy–primary amine, epoxy–secondary amine, and epoxy–surface reactions, together with the atomic number density profiles $\rho(z)$ of epoxy, amine, and substrate. Interestingly, the reaction density near the interface is lower than in the bulk, even though the local atomic densities of both epoxy and amine are higher near the interface. This finding is consistent with prior studies reporting densification near solid substrates, where restricted segmental mobility leads to reduced reaction rates despite high local concentrations. In other words, our simulations confirm that interfacial reactions are suppressed in the near-surface region, in agreement with previous experimental and theoretical reports of reduced reactivity at polymer/solid interfaces. ⁴⁹

A representative snapshot is shown in Figure 5(c). Direct visualization confirms that epoxy ring-opening reactions occur at the CuO–OH surface, leading to covalent bond formation between the epoxy network and substrate. Even though the number of such interfacial bonds is small, their formation connects the polymer network to the substrate and is expected to significantly enhance adhesion, consistent with experimental observations.

Figure 5: Epoxy curing at the hydroxylated CuO interface. (a) Time evolution of relative concentrations $c(t)/c(0)$ of epoxy oxygen ($O_{\mathrm{e}}$), primary amine nitrogen ($N_{\mathrm{p}}$), secondary amine nitrogen ($N_{\mathrm{s}}$), and surface hydroxyl hydrogen ($H_{\mathrm{CuOH}}$). (b) Depth-resolved reaction density $\rho_{\mathrm{rxn}}(z)$ for epoxy–primary amine (top), epoxy–secondary amine (middle), and epoxy–surface reactions (bottom), together with atomic number density profiles $\rho(z)$ of epoxy (blue), amine (orange), and CuO substrate (red). (c) Representative snapshot of the epoxy/CuO interface, with magnified view showing covalent bond formation between epoxy and surface hydroxyl groups.

5 Conclusion

The proposed framework successfully captures relative reactivity trends and qualitative mechanistic features. However, several limitations remain to be addressed. First, applying a bias potential disrupts the correspondence between simulation and physical time, preventing direct reproduction of absolute rate constants. Second, although propagation events are systematically accelerated, the current implementation excludes other relevant pathways, such as chain transfer and termination, which are essential for describing molecular-weight distributions in radical polymerization. Third, the accuracy of predicted barriers and kinetics is ultimately bounded by the fidelity of the underlying uMLIP, which, in our case, does not fully resolve subtle energetic differences across monomers.

Thus, the following directions are recommended for future work: (i) development of higher-fidelity uMLIPs trained on more diverse and accurate reference data, including reaction barriers and transition-state information, to improve quantitative reliability; (ii) extension of the reaction catalog to incorporate side reactions and termination events, enabling a more comprehensive description of experimentally observed polymerization kinetics; and (iii) exploration of hybrid approaches that combine acceleration with approximate time rescaling or kinetic reconstruction, thereby bridging the gap between qualitative trend capture and quantitative rate prediction. In addition, coupling the present scheme with enhanced conformational sampling (e.g., replica exchange or basin-hopping strategies) may further improve the robustness of predicted mechanistic pathways.

Overall, the TDBB approach, combined with a transferable uMLIP, provides a practical and versatile platform for simulating polymerization and curing processes. Its demonstrated ability to reproduce key features across diverse systems without system-specific parametrization highlights its promise as a foundation for developing more quantitatively predictive and chemically comprehensive methodologies.

Additional figures, reaction pathway analyses, and simulation details are provided in the Supporting Information (PDF) available at the ACS Publications website.

We express our sincere gratitude to So Takamoto of Preferred Networks, Inc., and Ryosuke Kuwabara of MITSUI & CO., LTD., and Sannohe Kunio of IBLC Co., Ltd. as well as to all their colleagues, for their valuable comments, insightful discussions, and continuous support.