Multiscale prediction of polymer relaxation dynamics via computational and data-driven methods

Abstract

We present a multiscale modeling approach that integrates molecular dynamics simulations, machine learning, and the elastically collective nonlinear Langevin equation (ECNLE) theory to investigate the glass transition dynamics of polymer systems. The glass transition temperatures (T_g) of four representative polymers are estimated using simulations and a machine learning model trained on experimental datasets. These predicted T_g values are used as inputs to the ECNLE theory to compute the temperature dependence of structural relaxation times and diffusion coefficients, and the dynamic fragility. The T_g values predicted from simulations show good quantitative agreement with experimental data. While machine learning tends to slightly overestimate T_g, the resulting dynamic fragility values remain close to experimental fragilities. Overall, ECNLE calculations using these inputs agree well with broadband dielectric spectroscopy results. Our integrated approach provides a practical and scalable tool for predicting the dynamic behavior of polymers, particularly in systems where experimental data are limited.

---

1 Introduction

Organic polymers such as poly(phenylene sulfide) (PPS), poly(isoprene) (PI), poly(propylene glycol) (PPG), and poly(butadiene) (PB) have been widely studied due to their valuable thermal, mechanical, and chemical properties. PI and PB are elastomeric polymers with double bonds in their backbone, which provide high chain flexibility. This makes them suitable for applications requiring elasticity such as tires, seals, and conveyor belts,^1,2 as well as biomedical devices.^3,4 PPG is a flexible polymer capable of forming hydrogen bonds and shows good chemical stability across a range of environments. It is commonly used in coatings, adhesives, and biocompatible products.^5–8 In contrast, PPS has a rigid backbone containing sulfur atoms. Since this structural feature improves its thermal stability and chemical resistance, PPS can be used under demanding engineering conditions.⁹ To improve the performance of such polymers, accurate determination and control of their glass transition temperature and related properties is essential. Despite extensive research, the fundamental nature of the glass transition remains a central challenge in the study of amorphous materials.

Several experimental techniques are commonly used to determine the glass transition temperature (T_g), including differential scanning calorimetry (DSC),¹⁰ dynamic mechanical analysis (DMA),^11,12 and broadband dielectric spectroscopy (BDS).¹³ Each method has distinct advantages and disadvantages. DSC has been widely adopted because of its simplicity and ability to detect heat flow associated with thermal transitions. However, when the glass transition overlaps with melting or crystallization events, it is hard to determine T_g. DMA has been widely used to measure the temperature dependence of viscoelastic behaviors of polymers with significant mechanical changes typically observed near the T_g. However, it requires precise preparation of sample dimensions and careful mechanical calibration to ensure accurate results. BDS data provide detailed insight into molecular relaxation dynamics over a broad frequency range. It allows us to calculate the temperature and pressure dependence of structural relaxation times, T_g, and the dynamic fragility. Despite its capabilities, BDS requires specialized instrumentation, can be costly, and is not universally applicable to all types of materials. These limitations have led to the development of computational approaches that include molecular dynamics (MD) simulations and theoretical models to complement or substitute the experimental determination of T_g and structural relaxation time.

MD simulations provide atomistic insight into the structural and thermal behavior of polymers by tracking the time evolution of particles under realistic interatomic interactions.¹⁴ The glass transition temperature can be estimated from such simulations by analyzing the temperature dependence of volume or density. However, MD simulations are limited to short timescales, typically up to 10⁵ picoseconds. Meanwhile, experimental observations often span milliseconds to hundreds of seconds. This means that MD simulations cannot fully capture the slow relaxation dynamics characteristic of systems near and below the T_g.

To address the timescale limitations of MD simulations, the elastically collective nonlinear Langevin equation theory was developed to describe structural relaxation in glass-forming systems. The ECNLE theory describes how local cage-scale hopping cooperates with long-range elastic displacements to govern structural relaxation. It predicts the temperature dependence of relaxation times and diffusion coefficients across a wide dynamic range from picoseconds to hundreds of seconds, which cover both simulation and experimental timescales. This theory has been successfully applied to various amorphous materials including polymers,^15,16 molecular liquids, organic compounds,^17,18 metallic glasses,¹⁹ and confined systems. ECNLE predictions have shown strong agreement with experimental data and/or simulation.^15–19 A key input to ECNLE calculations is T_g, which provides the thermal mapping from density to temperature. However, for newly synthesized or computationally designed materials, experimental T_g values are often unavailable and this restricts the applicability of the theory.

To overcome this challenge, in this study, we present an integrated approach that combines MD simulations, machine learning, and ECNLE theory to investigate the glassy dynamics of representative polymer systems. First, MD simulations and trained ML models are used to estimate the T_g values of each polymer system based on structural or compositional information. These predicted T_g values and experimental counterparts are then used as inputs for the ECNLE theory through the thermal mapping. We subsequently compute the temperature dependence of the structural relaxation time and diffusion coefficient, and the dynamic fragility. Finally, we validate our predictions by comparing them against available experimental data. This integrated MD–ML–ECNLE approach provides a predictive and scalable platform for understanding and designing amorphous materials, particularly when experimental data are limited or unavailable.

2 Theoretical background

In this section, we present the theoretical and computational methods used in this study. We first describe how the ECNLE theory is applied to predict the temperature dependence of structural relaxation dynamics. We then introduce MD simulations and ML models for estimating the T_g of polymers. By integrating these techniques, we establish a multiscale framework capable of predicting both the T_g and relaxation behavior for polymer systems without relying on adjustable or experimental parameters. Each approach is described in the following subsections.

2.1 The ECNLE theory

In the ECNLE theory, activated events in amorphous systems are analyzed by describing the material as a dense hard-sphere fluid system.^15–30 The model treats particles as hard spheres with a number density, ρ, and a characteristic diameter, d. The dynamics of a particle are influenced by short-range interactions with neighboring particles, thermal noise, and frictional damping. These effects are collectively described by a nonlinear Langevin equation. Solving this equation provides the dynamic free energy, F_dyn(r), which quantifies the effective confinement or caging constraints imposed on the tagged particle by its surroundings^15–30where k_B is the Boltzmann constant, T is the ambient temperature, r is the displacement, q is the wavevector, S(q) is the static structure factor which can be obtained using the Percus–Yevick theory,³² and Φ = ρπd³/6 is the volume fraction. The first term on the right-hand side of eqn (1) represents the ideal fluid-like contribution, whereas the second term captures the caging constraints arising from interactions with neighboring particles.

The dynamic free energy provides key insights into the local particle dynamics in dense fluids. At low volume fractions (Φ < 0.432), the dynamic free energy F_dyn(r) decreases monotonically with increasing r and has no local barrier. In this regime, particles are not dynamically localized and diffuse freely due to weak interparticle interactions. This corresponds to fluid-like behavior.¹⁵ As the system becomes denser, a local barrier emerges in the dynamic free energy due to stronger interactions between the tagged particle and its nearest neighbors. These interactions give rise to transient cages that temporarily confine the particle's motion. The radius of the particle cage, r_cage, is determined as the position of the first minimum in the radial distribution function, g(r). This function is related to the static structure factor S(q) via an inverse Fourier transform. As density increases, the local barrier in the dynamic free energy becomes higher. The height of this barrier is given by F_B = F_dyn(r_B) − F_dyn(r_L), where r_L is the localization length and r_B is the position of the barrier maximum in F_dyn(r) as illustrated in Fig. 1. The jump displacement Δr = r_B − r_L represents the distance that the particle must move to escape its cage.

Fig. 1 Schematic illustration of the dynamic free energy in the ECNLE theory. Key length and energy quantities are defined.

As a particle attempts to escape its cage, the surrounding particles must undergo cooperative rearrangements. This process leads to a slight expansion of the cage surface and generates a radially symmetric displacement field, u(r), that propagates outward from the particle cage. Based on Lifshitz's continuous mechanics analysis,³³ the displacement field is

where Δr_eff is the amplitude of the cage expansion, which is

The propagation of the displacement field induces vibrations in the surrounding particles. Since u(r) is small, these vibrations can be approximated as harmonic and the elastic energy associated with a single particle is given by with K₀ = |∂²F_dyn(r)/∂r²|_r=rL being the spring constant. To quantify the effect of collective motion on relaxation, the total elastic energy is obtained by summing the contributions from all particles outside the cage^15–30

Based on Kramers’ theory, we can calculate the structural relaxation time via^15–19

where τ_s is the short-time relaxation and is the magnitude of the curvature of the dynamic free energy at the peak of the barrier. The analytical expression for τ_s has been discussed in detail in prior studies.^24–27

Note that the calculations above provide τ_α as a function of volume fraction Φ. To quantitatively compare with experimental data, it is necessary to map Φ to temperature. To achieve this, Mirigian and Schweizer developed a thermodynamic mapping scheme for rigid molecular liquids.^20–22 This method relies on matching the dimensionless amplitude of long-wavelength density fluctuations, defined as S^expt₀(T) = ρk_BTκ_T, where κ_T is the isothermal compressibility. This quantity reflects nanoscale thermal density fluctuations and can be obtained from the experimental equation of state. In the hard-sphere fluid model with the Percus–Yevick approximation, the corresponding zero-wavenumber structure factor is given by . By equating S^HS₀ with S^expt₀(T),^15–30 we can determine the volume fraction as a function of temperature (Φ(T)). This thermal mapping allows us to calculate the temperature dependence of τ_α within the ECNLE theory.

Although ECNLE calculations associated with Schweizer's thermal mapping have been successful in predicting the dynamics of many glass-forming materials^21,24 without any adjustable parameter, its assumption of a universal link between local and long-range dynamics has a limitation. In reality, variations in microstructure and chemical interactions lead to significant differences in relaxation behavior, particularly in the dynamic fragility of polymers with different chemistries. To address this, Xie and Schweizer extended the ECNLE theory by introducing a non-universal coupling between local hopping and collective elastic barriers.³¹ They argued that the microscopic hopping distance needed for a particle to escape its cage is specific to the chemical structure and nanoscale rearrangements of the material. This chemical specificity is encoded into the theory through a single material-dependent parameter, a_c, which scales the collective elastic barrier. The new elastic barrier becomes F_e,new = a_cF_e, and the structural relaxation time in eqn (5) is modified as^15–30

The extended theory significantly improves the predictions for the glass transition temperature, dynamic fragility, and temperature dependence of τ_α simultaneously. However, calculations for the extended theory are highly dependent on equation-of-state data and the adjustable parameter a_c, which limits its practical application.

To address the limitations of previous thermal mapping schemes, we introduced an alternative approach based on the principle of thermal expansion

where T_g is the glass transition temperature defined at the structural relaxation time of 100 s, Φ_g is the corresponding volume fraction where τ_α = 100 s, Φ₀ = 0.5 is a characteristic volume fraction, and β ≈ 12 × 10⁻⁴ is the effective thermal expansion coefficient assumed constant for amorphous materials. Eqn (7) shows that Φ can be directly mapped to temperature using only T_g as the input. Importantly, Φ_g depends on the material-specific parameter a_c through eqn (6). For a_c = 1, the value of Φ_g is 0.6157, which is consistent with prior studies. This mapping significantly simplifies input requirements and is well suited for high-throughput application to diverse systems. The T_g value can be obtained experimentally via DSC or BDS,^10,13 or alternatively estimated using computational simulations³⁴ or machine learning methods.^19,28

The relationship between the diffusion coefficient (D) and the structural relaxation time is described by^19,35,36

In several previous studies,^19,35,36 diffusion coefficients calculated using the ECNLE theory have shown good agreement with experimental measurements for various glass-forming systems.^19,35,36

The ECNLE theory has several distinct advantages over other theoretical models such as mode-coupling theory (MCT),^37,38 random first order transition (RFOT) theory,³⁹ and the Adam–Gibbs (AG) model⁴⁰ to describe the thermal dependence of structural relaxation in supercooled liquids and polymers. The MCT captures short-time caging and two-step relaxation behavior and breaks down below the crossover temperature because it neglects activated hopping processes. The RFOT assumes the existence of cooperatively rearranging regions and a free-energy landscape without a precise microscopic formulation for the energy barrier. Similarly, although the AG model links relaxation time to configurational entropy, it remains largely phenomenological and lacks a clear microscopic basis. In contrast, ECNLE explicitly accounts for both local hopping out of transient cages and long-range collective elastic distortions. This dual contribution allows it to accurately model relaxation dynamics over a wide temperature range, including near and below the glass transition. In general, the ECNLE theory provides a more physical and predictive framework to understand relaxation time, fragility, and glass formation in a broad range of amorphous materials.^15–30

2.2 Molecular dynamics simulations

To estimate the glass transition temperature of PPS, PI, PPG, and PB, molecular dynamics simulations were conducted using the large-scale atomic/molecular massively parallel simulator (LAMMPS).⁴¹ Although the MD cooling rates are significantly faster than those in experiments, our aim is not to exactly reproduce experimental T_g values. Instead, we use the MD-derived T_g values as system-specific inputs to ECNLE theory to model temperature-dependent dynamic properties without relying on experimental data. This allows us to build a self-contained predictive framework. The calculation procedure involves three key steps: (1) constructing and minimizing initial configurations, (2) equilibrating and cooling the system under constant pressure and temperature, and (3) determining T_g based on the simulated temperature dependence of the system volume.

For each polymer, simulation cells containing either 10 or 20 chains with a polymerization degree of N = 10 were generated. These structures were packed into a periodic cubic box of dimensions 60 × 60 × 60 ų to provide an initial mass density close to 1.0 g cm⁻³. The Dreiding force field was employed due to its applicability to both organic molecules and polymeric systems.⁴² This force field accounts for bonded interactions (bond, angle, dihedral) and nonbonded interactions (van der Waals and hydrogen bonding), enabling realistic modeling of polymer conformations and intermolecular forces.

The chosen polymerization degree of N = 10 corresponds to molecular weights of approximately 631, 683, 656, and 570 g mol⁻¹ for PPS, PI, PPG, and PB, respectively. While these values lie in the low-to-intermediate molecular weight regime, they are comparable to the experimental molecular weights of PI (1040 g mol⁻¹), PPG (192 g mol⁻¹), and PB (777 g mol⁻¹).⁴³ This allows for meaningful comparisons of their T_g values and relaxation dynamics.⁴³ Prior studies⁴³ have shown that the glass transition temperature of PPG remains relatively insensitive across a broad molecular weight range (134–18 000 g mol⁻¹), and the T_g values of PB with molecular weights near 570 and 777 g mol⁻¹ are nearly identical. Thus, for these polymers, although the molecular weight has significant effects on the glass transition temperature and fragility,^43–45 our simulation conditions are representative of experimentally accessible behavior. However, for PPS, the simulated molecular weight is significantly lower than the experimental value (∼44 000 g mol⁻¹),⁴³ and quantitative comparisons of T_g and fragility in this case are cautiously studied.

Following energy minimization, each polymer system was equilibrated at a high temperature selected to ensure molecular mobility and thermal stability. The equilibration temperatures were set to 340 K for PPS, 300 K for PI and PPG, and 330 K for PB. Simulations were performed in the NPT ensemble using the Nosé–Hoover thermostat and barostat for 250 ps to maintain constant thermodynamic conditions. Subsequently, the systems were cooled in a stepwise manner to target temperatures of 150 K for PPS, 100 K for PI and PPG, and 50 K for PB. The cooling process was executed over multiple temperature intervals of 20–25 K with each interval involving 200 ps of equilibration. At each step, the specific volume was averaged over the final 100 ps of simulation. To determine T_g, the temperature–volume data were analyzed by fitting two linear regimes corresponding to high- and low-temperature regimes. The intersection of these two linear fits was taken as the glass transition temperature for each polymer.

2.3 Machine learning models

In addition to MD simulation approaches, machine learning can be exploited to predict the T_g value of polymers. Among various ML techniques, Gaussian process regression (GPR) has recently been shown to provide highly accurate T_g predictions for polymeric systems.²⁸ Based on this, we selected a GPR model using a custom Tanimoto kernel to construct the predictive model in the present study. Since the experimental T_g values of PPS, PI, PPG, and PB are all below 320 K, we selected a subset of 1586 polymers with experimental T_g values under 320 K from a larger dataset comprising 7174 polymers with T_g values ranging from 134 K to 768 K.⁴⁶

The chemical structures of the polymers were first represented using the simplified molecular input line entry system (SMILES) format,⁴⁷ which encodes how atoms are connected within each molecule. These representations were then processed using RDKit⁴⁸ to generate quantitative molecular fingerprints that capture key features of the molecular topology, such as the arrangement of atoms and bonds. Each fingerprint consists of 2048 binary values (0 and 1), where each bit indicates the presence or absence of a specific structural pattern. These molecular descriptors were used as the input features for training and testing the machine learning model. All dataset and code used in this study are openly available in the GitHub repository⁴⁹ to ensure full reproducibility and transparency. The dataset was randomly divided into training (80%) and testing (20%) sets. To evaluate the predictive accuracy, we used two common performance metrics: the coefficient of determination (R²), which indicates how well the model captures variance in the data, and the root mean square error (RMSE), which measures the average deviation between predicted and actual T_g values.

The histograms in Fig. 2 show the frequency distributions of actual and predicted T_g from our machine learning model. For the training set (Fig. 2a), the distribution of the predicted T_g values shows excellent overlap and is nearly indistinguishable from the actual T_g distribution. Crucially, for the test set (Fig. 2b), our predictions have a strong coherence between the actual and predicted distributions. Specifically, the predicted T_g distribution in the test set successfully captures the overall shape, range, and modal features of the experimental data. This strong agreement in both training and test sets suggests that our model is reliable for predicting the T_g values.

Fig. 2 Distribution of actual and predicted glass transition temperatures of polymers for (a) training data and (b) test data.

3 Results and discussion

Fig. 3 shows the temperature dependence of volume (in ų) obtained from MD simulations for four polymer systems. Each system was equilibrated at a high temperature and subsequently cooled in discrete steps to a lower target temperature. This cooling process resulted in a gradual volume decrease. This change indicates the transition from a supercooled liquid with a higher thermal expansivity to a glassy state with significantly reduced molecular mobility. Linear fits were applied to the averaged volume data in both high- and low-temperature regions, and the intersection point of these lines was taken as the predicted T_g. The calculated T_g values are 227 K for PPS, 221 K for PI, and 192 K for PPG.

Fig. 3 Temperature dependence of the simulated volume from MD simulations for (a) PPS, (b) PI, (c) PPG, and (d) PB. Blue squares indicate raw data at each temperature step, while black circles indicate the averaged values used for analysis. Linear fits to the low-temperature and high-temperature regimes are shown as brown and black lines, respectively. The intersection of these two lines corresponds to the predicted T_g of the polymer.

Our machine learning model, which is based on GPR, predicts T_g with R² of 0.78 and RMSE of 19.5 K. Fig. S2 in the SI shows the predicted versus experimental T_g values for the training and test sets of the low-T_g subset to confirm comparable errors for both. This performance indicates reasonable predictive accuracy within the low-T_g range relevant to our selected polymers. Compared to a previous GPR study,²⁸ which reported R² of over 0.89 and an RMSE below 37 K for a full dataset spanning T_g values from 134 to 768 K, our model has a lower overall R² but a smaller RMSE. This difference arises because we focus our training on a subset of 1586 polymers with T_g values below 320 K to better represent the temperature range of interest. As detailed in Tables SI–SIII of the SI, cross-validation shows that this choice yields slightly lower global R² due to the narrower T_g range, but improves accuracy (lower RMSE) for low-T_g polymers and avoids the systematic upward shift in predicted T_g that occurs when higher-T_g polymers are included in the training set.

Several recent studies have applied more complex machine-learning architectures to predict polymer T_g. A graph neural network (GNN) model was specifically developed for polyimides, focusing on high-T_g materials within a relatively narrow chemical space.⁵⁰ However, its maximum R² of 0.86 implies a non-negligible amount of unexplained variance and does not provide higher accuracy than our GPR model. Park et al. combined a graph convolutional network (GCN) with a neural-network (NN) regressor to determine T_g and other properties over a dataset with a similar T_g range and distribution to our full database.⁵¹ Although their deep-learning architecture is more complex, the resulting prediction accuracy is only comparable to that of our simpler GPR approach. Cheminformatics models coupled to coarse-grained MD simulations⁵² can identify physical descriptors, but these simulations are time-consuming and do not always reproduce experimental behavior accurately. As a result, this approach yields only 0.47 ≤ R² ≤ 0.74 for T_g prediction across diverse polymers and is difficult to deploy in a high-throughput setting. More recently, quantum-chemical-augmented NN/GNN frameworks⁵³ that rely on quantum-chemistry-derived local-cluster descriptors report R² ≈ 0.74 for T_g and require additional quantum-chemical calculations for each polymer. In contrast, our GPR model operates on simple binary fingerprints generated from SMILES strings, avoids any MD- or QC-based feature generation, and still achieves higher accuracy in the low-T_g regime of interest. This balance of accuracy, simplicity, and computational efficiency makes the GPR framework the most appropriate choice for the present work.

To further evaluate the accuracy of both ML and MD approaches, we compare their predicted T_g values with experimental data in ref. 43. As summarized in Table 1, the MD results show good quantitative agreement with experimental values and deviations are typically within 10–15 K. The GPR model predicts the T_g values for PPS and PI with moderate accuracy. The ML predictions overestimate the experimental counterparts by 20 K and 36 K, respectively. For PPG and PB, the deviations are larger, on the order of 80 K. This likely indicates both the limited number of structurally similar low-T_g polymers in the training data and the intrinsic scatter of the experimental T_g database. Hence, in this work we use the GPR model mainly to capture overall trends in the low-T_g regime and to provide an additional input for testing ECNLE predictions, rather than as a highly accurate predictor for each individual flexible polymer. While MD simulation provides more accurate absolute T_g predictions, the GPR model remains a useful and efficient tool for trend prediction and rapid screening, especially when experimental data are unavailable. A systematic exploration of more expressive architectures and additional descriptors as shown in ref. 50–53 to further reduce these errors is an important direction for future work, but is beyond the scope of this study, which focuses on testing ECNLE predictions against different T_g inputs.

Table 1 The T_g values of selected polymers obtained from MD simulations and ML calculations in this work and in a prior experimental study⁴³

Material	Tg,expt (K)	Tg,MD (K)	Tg,ML (K)
PPS	229	227	249
PI	185	221	240
PPG	191	192	277
PB	166	181	249

---

Fig. 4 shows the structural relaxation time as a function of 1000/T for PPS, PI, PPG, and PB. The relaxation times were first calculated as a function of volume fraction using eqn (1)–(5), and then converted to temperature using the thermal mapping described in eqn (7), where T_g is a key input. The T_g values used in this mapping were obtained from MD simulations, ML predictions, and experimental data as listed in Table 1. ECNLE predictions using experimental T_g values (T_g,expt) quantitatively agree with BDS data in ref. 43 for all polymers. When using MD-predicted T_g values (T_g,MD), the predicted relaxation times remain consistent with experimental data for most systems, although a noticeable deviation is observed for PI. This finding is predictable since the difference between T_g,MD and T_g,expt is largest for PI. In contrast, using ML-predicted T_g values (T_g,ML) leads to greater discrepancies between ECNLE results and experimental data, particularly at lower temperatures.

Fig. 4 Temperature dependence of the structural relaxation time for (a) PPS, (b) PI, (c) PPG, and (d) PB calculated using eqn (1)–(5) and (7). The T_g used in the thermal mapping (eqn (7)) is obtained from MD simulations, ML predictions, and experimental data in ref. 43. Experimentally, the molecular weights obtained were 44 000 g mol⁻¹ for PPS, 1040 g mol⁻¹ for PI, 192 g mol⁻¹ for PPG, and 777 g mol⁻¹ for PB. In contrast, the corresponding values derived from MD simulation were 631, 683, 656, and 570 g mol⁻¹, respectively.

To facilitate direct comparison with experimental results, the relaxation time data from Fig. 4 are replotted as a function of T_g/T in Fig. 5. This normalized representation shows that the results from the MD simulations and the GPR model agrees more closely with the experimental data than when plotted in Fig. 4. Moreover, plotting log τ_α versus T_g/T allows us to calculate the dynamic fragility, which is defined as .⁵⁴ This quantity quantifies the steepness of the increase in the structural relaxation time near T_g. Based on the value of the dynamic fragility index m, glass-forming materials are typically categorized into three groups. Materials with m < 30 are considered “strong” glass formers and their τ_α(T) exhibits relatively Arrhenius-like behavior. Those with m > 100 are classified as “fragile” and the temperature dependence of the structural relaxation time shows highly non-Arrhenius dynamics near T_g. Systems with m values between 30 and 100 are regarded intermediate glass formers.

Fig. 5 The logarithm of structural relaxation times as a function of T_g/T for (a) PPS, (b) PI, (c) PPG, and (d) PB from the same data as in Fig. 4.

Table 2 summarizes the dynamic fragility calculated using ECNLE theory using different T_g inputs and experimental fragility values. In most cases, the calculated m values are lower than the experimental values, but the deviations remain within an acceptable range. When using experimental T_g as the input, the calculated m values for PPS and PI are within 10–15% of the experimental fragilities. Although the T_g,MD values are close to the experimental T_g, the fragilities calculated using T_g,ML actually show better agreement with experimental values across all polymers. This suggests that machine learning, even with slightly overestimated T_g values, can effectively predict the dynamic fragility. The observed differences may result from experimental uncertainties, the sensitivity of glassy dynamics to small variations in T_g, and approximations inherent in the theoretical model. We also train our ML model based on the full dataset of polymer and find that the calculated fragilities deviate more significantly from experimental values (see the SI). This further validates the importance of appropriate data selection. Overall, these findings reveal that combining ECNLE theory with T_g inputs from either MD or ML provides a reliable and practical approach for predicting the thermal behavior of polymer glasses.

Table 2 The dynamic fragility values of selected polymers computed using ECNLE theory with T_g inputs from experiments, MD simulations, and ML predictions as listed in Table 1, compared with experimental fragility

Material	mexpt	m (Tg,expt)	m (Tg,MD)	m (Tg,ML)
PPS	116	84	84	89
PI	80	75.4	82	89
PPG	79	67	67	103
PB	84	64	67	92

---

To examine how specific molecular characteristics such as intramolecular flexibility, side-chain dynamics, and intermolecular interactions affect the τ_α(T) calculations, we adopted the method proposed by Xie and Schweizer.³¹ This approach uses a material-specific parameter a_c to adjust the contribution of the collective elastic barrier to the total energy barrier (F_B + a_cF_e). All above calculations in this work use the default value a_c = 1. Fig. 4 presents the structural relaxation times of PPS, PI, PPG, and PB as a function of 1000/T computed using the ECNLE theory with varying values of the parameter a_c based on eqn (1)–(4), (6), and (7). Adjusting a_c improves quantitative agreement between numerical results and experimental data. This finding indicates that accurate prediction of T_g alone is insufficient to quantitatively capture the temperature dependence of structural relaxation time, as the dynamic behavior is also sensitive to how the collective elastic barrier is treated. However, even with optimal tuning, the model does not fully reproduce the relaxation behavior over the entire temperature range. It is also important to note that the predictions of this extended approach strongly depend on the selection of a_c, which currently lacks a first-principles basis and must be empirically adjusted. The dependence of this method on fitting limits its use to systems for which experimental data already exist (Fig. 6).

Fig. 6 Temperature dependence of the structural relaxation time for (a) PPG, (b) PI, (c) PPS, and (d) PB. Data points correspond to experimental data in ref. 43. Solid curves are ECNLE theoretical predictions computed using eqn (1)–(4), (6) and (7) with experimental T_g values and different a_c parameters.

The above results indicate that using a_c = 1 in the ECNLE theory provides a reasonable prediction of τ_α(T). Based on this, we apply eqn (8) to compute the diffusion coefficient as a function of T_g/T for PPS, PI, PPG, and PB. The particle diameter d, which is used to calculate the diffusion coefficient, was estimated for each polymer using the atomic structure of its monomer derived from SMILES codes and processed with the RDKit software. The computed diffusion coefficients are presented in Fig. 7. Using the T_g values obtained from MD, ML, and experimental data, we find that the diffusion coefficient decreases rapidly as temperature decreases. Notably, the ECNLE predictions using T_g,MD and T_g,expt are relatively close, despite the quantitative differences between these T_g values. Although the ECNLE theory is a microscopic statistical mechanical theory and does not explicitly account for spatially heterogeneous dynamics or the decoupling behavior between diffusion and relaxation time observed near the glass transition temperature,^55–57 our calculated diffusion coefficients show good agreement with experimental data for PPG and PB provided in ref. 58. This agreement reinforces the physical relevance of our predictions and suggests that, for the systems studied here, the ECNLE theory can provide a quantitatively accurate description of average translational dynamics.

Fig. 7 The predicted temperature dependence of the diffusion coefficient for (a) PPS, (b) PI, (c) PPG, and (d) PB calculated using eqn (1)–(5), (7) and (8). The T_g used in the thermal mapping (eqn (7)) is obtained from experimental, MD simulations, ML predictions, and experimental data.⁵⁸ The molecular weights for PPG were 770 g mol⁻¹ experimentally and 656 g mol⁻¹ via MD simulation, while the values for PB were 466 g mol⁻¹ and 570 g mol⁻¹, respectively.

4 Conclusion

In conclusion, we have developed a multiscale approach combining MD simulations, machine learning, and the ECNLE theory to study glass transition and relaxation dynamics in polymers. The glass transition temperatures of four representative polymers were estimated using both MD simulations and the GPR model trained on experimental data. These values were then used to calculate temperature-dependent structural relaxation times, diffusion coefficients, and dynamic fragility using the ECNLE theory. The MD predicted values T_g agree well with the experimental data within 10–15 K. Meanwhile, the T_g predictions based on machine learning are relatively higher than the experimental T_g values. ECNLE calculations using the T_g values from the MD simulations and experiments show good agreement with the BDS data when plotted as log τ_α versus 1000/T. In the Angell representation (log τ_α versus T_g/T), the ECNLE results based on T_g,MD, T_g,ML and T_g,expt are all quantitatively in agreement with the experimental relaxation data. This representation can reliably estimate the dynamic fragility. In particular, the fragilities calculated using T_g,ML show better agreement with the experimental values compared to those based on T_g,MD. This suggests that machine learning provides not only trend preservation but also quantitative accuracy in fragility prediction. Although the extended ECNLE model allows tuning through a material-specific parameter a_c to scale effects of collective dynamics on the glass transition, we found that the standard formulation (a_c = 1) already provides reasonable accuracy without any adjustable parameters. Overall, our study clearly validates the predictivity of our integrated approach, especially when experimental data are limited.

Conflicts of interest

The authors have no conflicts to disclose.

Data availability

Data and code for this article are available at GitHub at https://github.com/NgoQue/GPR_Model.

Acknowledgements

This research was funded by the Vietnam National Foundation for Science and Technology Development (NAFOSTED) under Grant No. 103.01-2023.62.

References

1. P. Zuo, A. Tcharkhtchi, M. Shirinbayan, J. Fitoussi and F. Bakir, Macromol. Mater. Eng., 2019, 304(5), 1800686.
2. P. S. Gopala Krishnan, K. Ayyaswamy and S. K. Nayak, J. Macromol. Sci., Part A: Pure Appl. Chem., 2013, 50(1), 128–138.
3. A. M. Díez-Pascual and A. L. Diez-Vicente, ACS Appl. Mater. Interfaces, 2014, 6(13), 10132–10145.
4. Y.-J. Li, K. H. Matthews, T.-M. Chen, Y.-F. Wang, M. Kodama and T. Nakaya, Chem. Mater., 1996, 8(7), 1441–1450.
5. J. A. Okolie, iScience, 2022, 25(9).
6. M. Sara, T. Rouissi, S. K. Brar and J. F. Blais, Plat. Chem. Biorefinery, 2016, 77–100.
7. R. K. Wang and J. B. Elder, Lasers Surg. Med., 2002, 30(3), 201–208.
8. Z. Li, Z. Zhang, K. L. Liu, X. Ni and J. Li, Biomacromolecules, 2012, 13(12), 3977–3989.
9. X. C. Tong, Advanced materials for thermal management of electronic packaging, Springer Science & Business Media, 2011, vol. 30.
10. M. S. Rahman, I. M. Al-Marhubi and A. Al-Mahrouqi, Chem. Phys. Lett., 2007, 440(4–6), 372–377.
11. W. Sun, A. P. Vassilopoulos and T. Keller, Int. J. Adhes. Adhes., 2014, 52, 31–39.
12. A. S. Cantor, J. Appl. Polym. Sci., 2000, 77(4), 826–832.
13. F. Kremer, J. Non-Cryst. Solids, 2002, 305, 1–9.
14. J.-L. Barrat, J. Baschnagel and A. Lyulin, Soft Matter, 2010, 6, 3430–3446.
15. A. Ghanekarade, A. D. Phan, K. S. Schweizer and D. S. Simmons, Proc. Natl. Acad. Sci. U. S. A., 2021, 118(31), e2104398118.
16. A. D. Phan and K. S. Schweizer, ACS Macro Lett., 2020, 9, 448–453.
17. A. D. Phan and K. S. Schweizer, Macromolecules, 2019, 52, 5192–5206.
18. A. D. Phan and K. S. Schweizer, J. Chem. Phys., 2024, 160, 074902.
19. A. D. Phan, D. T. Nga, N. T. Que, H. Peng, T. Norhourmour and L. M. Tu, Comput. Mater. Sci., 2025, 251, 113759.
20. S. Mirigian and K. S. Schweizer, J. Phys. Chem. Lett., 2013, 4(21), 3648–3653.
21. S. Mirigian and K. S. Schweizer, Macromolecules, 2015, 48(6), 1901–1913.
22. S. Mirigian and K. S. Schweizer, J. Chem. Phys., 2014, 141(16), 161103.
23. A. D. Phan and K. S. Schweizer, J. Chem. Phys., 2019, 150(4), 044508.
24. S. Mirigian and K. S. Schweizer, J. Chem. Phys., 2014, 140, 194506.
25. S. Mirigian and K. S. Schweizer, J. Chem. Phys., 2014, 140, 194507.
26. A. D. Phan and K. S. Schweizer, Macromolecules, 2018, 51, 6063–6075.
27. A. D. Phan and K. S. Schweizer, Macromolecules, 2019, 52, 5192–5206.
28. A. D. Phan, N. T. Que, N. T. T. Duyen, P. Thanh Viet, Q. K. Quang and B. Mei, J. Appl. Phys., 2025, 138, 44703.
29. S. Mirigian and K. S. Schweizer, J. Chem. Phys., 2015, 143, 244705.
30. A. D. Phan and K. S. Schweizer, J. Phys. Chem. B, 2018, 122(35), 8451–8461.
31. S.-J. Xie and K. S. Schweizer, Macromolecules, 2016, 49, 9655–9664.
32. J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, Academic Press, London, 2006.
33. L. D. Landau and E. M. Lifshitz, Theory of Elasticity, 3rd edn, Permagon Press, London, 1975.
34. J. Han, R. H. Gee and R. H. Boyd, Macromolecules, 1994, 27(26), 7781–7784.
35. A. D. Phan, K. Koperwas, M. Paluch and K. Wakabayashi, Phys. Chem. Chem. Phys., 2020, 22, 24365.
36. A. D. Phan, N. K. Ngan, D. T. Nga, N. B. Le and C. V. Ha, Phys. Status Solidi RRL, 2022, 16, 2100496.
37. W. Götze, Complex Dynamics of Glass-Forming Liquids:A Mode-Coupling Theory, Oxford University Press, 2009.
38. L. M. C. Janssen, Front. Phys., 2018, 6, 97.
39. T. R. Kirkpatrick, D. Thirumalai and P. G. Wolynes, Phys. Rev. A, 1989, 40, 1045.
40. G. Adam and J. H. Gibbs, J. Chem. Phys., 1965, 43, 139–146.
41. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
42. S. L. Mayo, B. D. Olafson and W. A. Goddard, J. Phys. Chem., 1990, 94, 8897–8909.
43. B. Schmidtke, M. Hofmann, A. Lichtinger and E. A. Rössler, Macromolecules, 2015, 48(9), 3005–3013.
44. X. Xu, J. F. Douglas and W.-S. Xu, Soft Matter, 2025, 21, 2664.
45. Q.-L. Yuan, X. Xu, J. F. Douglas and W.-S. Xu, Macromolecules, 2025, 58, 9528.
46. M. J. Uddin and J. Fan, Polymers, 2024, 16(8), 1049.
47. D. Weininger, J. Chem. Inf. Comput. Sci., 1988, 28, 31–36.
48. RDKit: Open-source cheminformatics software. Available at: https://www.rdkit.org/.
49. Data and code for this article are available at GitHub: https://github.com/NgoQue/GPR_Model.
50. H. Qiu, X. Qiu, X. Dai and Z.-Y. Sun, J. Mater. Chem. C, 2023, 11, 2930–2940.
51. J. Park, Y. Shim, F. Lee, A. Rammohan, S. Goyal, M. Shim, C. Jeong and D. S. Kim, ACS Polym. Au, 2022, 2, 213–222.
52. A. Karuth, A. Alesadi, W. Xia and B. Rasulev, Polymer, 2021, 218, 123495.
53. Y.-W. Zhang, Z. Ye, D. Hu, S. Qi, Z. Sun, J. Yang, Y. Ma, W. Zhang, J. Zhang and Z. Li, ACS Appl. Polym. Mater., 2025, 7, 12176–12186.
54. J. Wu, G. Huang, L. Qu and J. Zheng, J. Non-Cryst. Solids, 2009, 355(34–36), 1755–1759.
55. X. Xu, J. F. Douglas and W.-S. Xu, Macromolecules, 2023, 56, 4929.
56. Q.-L. Yuan, X. Xu, J. F. Douglas and W.-S. Xu, J. Phys. Chem. B, 2024, 128, 9889.
57. Z. Yang, X. Xu, J. F. Douglas and W.-S. Xu, J. Chem. Phys., 2024, 160, 041101.
58. N. Fatkullin, S. Stapf, M. Hofmann, R. Meier and E. A. Rössler, J. Non-Cryst. Solids, 2015, 407, 309–317.

---

Footnote

† Present address: Center for Materials Innovation and Technology, VinUniversity, Hanoi 100000, Vietnam.

---