Systematic coarse-graining of epoxy resins with machine learning-informed energy renormalization

The study addresses how to develop chemistry-specific, computationally efficient coarse-grained (CG) molecular models of epoxy resins that remain accurate across varying degrees of crosslinking (DC). All-atom (AA) molecular dynamics can predict DC-dependent properties but are too expensive for high-throughput design. Existing CG models for epoxies often focus on highly crosslinked networks and lack transferability across temperatures and curing states. The research question is how to systematically calibrate a CG force field that reproduces AA-predicted density, dynamics, and mechanical properties across DC for different epoxy chemistries. The authors propose combining energy renormalization (ER)—adjusting non-bonded interactions to account for entropic differences between AA and CG—with machine learning (Gaussian process surrogates) for multi-response, DC-dependent functional calibration. The work targets DGEBA-based epoxies cured with PACM or Jeffamine D400, reflecting practical interest in partially cured states and mixed hardeners that yield nanoscale heterogeneity and enhanced performance.

Prior AA simulations have effectively predicted DC effects on glass transition, thermal expansion, and elastic response, and have been used to study fracture in epoxy composites. However, AA simulations remain computationally expensive. Previous CG efforts for epoxies captured structural or thermo-mechanical properties mostly for highly crosslinked states and lacked transferability over temperature or curing degree. Energy renormalization has been used to achieve temperature transferability in simpler polymers and glass-formers by rescaling non-bonded interactions to match picosecond-scale dynamics (Debye-Waller factor). ML techniques have been explored in CG modeling, including particle swarm optimization for temperature-dependent calibration at select curing states, but a general framework for multi-property, DC-dependent calibration across different cure chemistries was missing. The present work builds on ER concepts and introduces Gaussian process-based functional calibration to address high-dimensional parameter spaces and multi-response targets with uncertainty quantification.

• Systems and chemistries: Epoxy networks based on DGEBA with two curing agents: PACM (4,4-diaminodicyclohexylmethane) or Jeffamine D400 (polyoxypropylene diamine) at stoichiometric ratios. AA networks prepared at DC from 0% up to 90% (PACM) and 95% (D400). CG mappings include 7 bead types, with bead assignments for DGEBA (beads 1–3), PACM (beads 4–5), and D400 (beads 5–7).
• All-atom (AA) simulations: LAMMPS with DREIDING force field. Crosslinking via Polymatic by creating epoxide–amine bonds within a 6 Å cutoff, 16 bonds per cycle, NPT relaxation at 600 K, 1 atm; two replicas for statistics. Annealing protocol includes heating to 600 K and 1000 atm, equilibration, quenching to 300 K, P=0 atm, and final equilibration. Properties measured: density, Debye–Waller factor ⟨u²⟩ from MSD at t ≥ 3 ps, Young’s modulus from tensile tests at 0.5 s⁻¹ using 0–2% strain slope, and yield stress via 0.2%-offset-like method (intersection of stress curve with modulus line shifted to 3% strain).
• Coarse-grained (CG) model: Bonds/angles from Boltzmann inversion (BI) of AA distributions; parameters listed (e.g., multiple bond and angle stiffness/equilibrium values). Non-bonded interactions via lj/gromacs with arithmetic mixing rules ε_ij = (ε_i ε_j)^1/2 and σ_ij = (σ_i σ_j)^1/2. Initial ε_i and σ_i from BI of AA radial distribution functions. Crosslinking in CG between DGEBA bead 3 and amine bead 5 (PACM/D400) within 15 Å, 10 bonds per cycle, with NPT relaxation at 1000 K, P=0 atm. CG annealing: equilibrations at 800 K/100 atm, then 500 K/0 atm, then 300 K/0 atm (each 2 ns). Dynamics and tensile tests performed similarly to AA (strain rate 0.5 s⁻¹).
• Parameter space and design of experiments: Treat non-bonded calibration as multi-objective optimization of 14 parameters [ε_i, σ_i] (i=1..7) across DC to match AA targets (density, ⟨u²⟩, Young’s modulus, yield stress). Parameter ranges for ε_i expanded based on matching AA dynamics at DC=0% or modulus at high DC, and ±20% around BI estimates for σ_i. Surrogate modeling required training data:
  • AA surrogates: 19 samples (PACM) and 20 samples (D400) across DC only.
  • CG surrogates: 700 simulations (PACM) in an 11D hypercube (DC plus 5 bead pairs), and 500 simulations (D400) in a 13D hypercube (DC plus 7 bead pairs), sampled via Sobol sequences.
• Gaussian process (GP) surrogates: Trained for both AA and CG systems to predict the four responses with predictive means and uncertainties. GP surrogates enable variance-based global sensitivity analysis (Sobol indices) and efficient calibration.
• Sensitivity analysis: Using GP surrogates to compute main and total sensitivity indices. Findings: σ_i dominate density; ε_i dominate dynamics and mechanical properties; DC strongly affects ⟨u²⟩ and yield stress; notable influence of σ_6 in DGEBA+D400 due to bead 6 prevalence and volume contribution.
• Functional calibration strategy: Initial attempt with sigmoids for DC dependence (by analogy to temperature dependence) proved too restrictive. Adopted radial basis functions (RBFs) with shared shape parameter ω and three centers at DC = 0%, 50%, 100% for each of 14 functions ε_i(DC), σ_i(DC). Calibration minimizes discrepancy between CG and AA GP predictions over DC for all responses. To avoid overfitting and ensure identifiability, uncertainty quantification was performed by evaluating how perturbations around optimal RBFs affect the probability of CG and AA producing the same target responses (likelihood-like objective), enabling approximate Bayesian computation to obtain means and uncertainty bands for each function.
• Model simplification via uncertainty-guided selection: Where uncertainty envelopes were narrow, function class was fixed (e.g., linear); where broad, simpler forms (constant or linear) were chosen without accuracy loss. Final simplified set reduced free parameters from 43 (full RBFs) to 21 by using constants or linear functions for most parameters and retaining an RBF only for ε_3.
• Final simplified parameterization (examples from Table 2):
  • ε_1(DC) = 1.07 (const), ε_2(DC) = 0.47 + 0.52·DC, ε_3(DC) = RBF with ω = −0.44 and centers (1.89, 2.21, 3.16), ε_4(DC) = 1.74, ε_5(DC) = 0.90 − 0.26·DC, ε_6(DC) = 0.40, ε_7(DC) = 0.12.
  • σ_1(DC) = 5.76 − 0.14·DC, σ_2(DC) = 5.71, σ_3(DC) = 3.65 − 0.45·DC, σ_4(DC) = 5.36, σ_5(DC) = 4.04, σ_6(DC) = 5.15, σ_7(DC) = 4.06.
• Validation: Compare GP-predicted CG responses and actual CG simulations against AA GP predictions across DC for both chemistries. Evaluate RMS percentage errors and confidence intervals. Also validate full MSD and stress–strain curves beyond the targeted scalar metrics.

• Multi-property, DC-transferable CG force field: The calibrated CG model accurately reproduces AA-predicted density, Debye–Waller factor ⟨u²⟩, Young’s modulus, and yield stress across DC from uncrosslinked to highly crosslinked networks for both DGEBA+PACM and DGEBA+D400.
• Accuracy: Average root mean squared percentage error (RMSPE) ≈ 10% across all targets and DC. High-DC Young’s moduli align with experiments (~2.5–3 GPa).
• Parameter simplification: Uncertainty-guided simplification reduced the calibration from 43 free parameters (RBFs for all 14 functions) to 21 parameters (constants/linear for 13 functions, RBF for ε_3) with negligible loss in accuracy.
• Sensitivity insights: σ parameters chiefly control density; ε parameters dominate dynamics and mechanics. DC is as influential as ε_i for ⟨u²⟩ and yield stress, less so for modulus. In DGEBA+D400, σ_6 strongly impacts all properties due to bead 6 prevalence and volume share.
• Computational efficiency: CG simulations are ~10^3 times faster than AA for the same system size, enabling larger-scale studies.
• Beyond target properties: The CG model reproduces overall tensile stress–strain curves and, for fully crosslinked networks, MSD behavior; a discrepancy appears at longer times for uncrosslinked systems where AA shows faster dynamics despite matched picosecond ⟨u²⟩.
• Generalizable framework: The ML-informed ER approach supports multi-objective calibration with uncertainty quantification and is extendable to other chemistries and property sets.

The work demonstrates that incorporating energy renormalization with Gaussian process-based functional calibration enables a CG force field whose non-bonded parameters vary with DC to reconcile entropic-resolution differences between AA and CG models. By simultaneously targeting density, dynamics (⟨u²⟩), stiffness, and yield strength, the CG model captures the complex, nonlinear dependence of properties on crosslinking and accommodates chemistry-specific effects (e.g., rigidity of PACM vs flexibility of D400). Sensitivity analysis elucidates how specific beads and parameters influence macroscopic responses, informing model parsimony. The uncertainty-aware simplification balances accuracy and identifiability, avoiding overfitting while maintaining multi-response fidelity. This methodology addresses the core challenge of transferability across curing states and provides a scalable path for high-throughput, physics-informed CG modeling of thermosets. The observed MSD discrepancy at long times for uncrosslinked systems highlights limits of matching only picosecond dynamics in chemically heterogeneous networks, guiding future inclusion of additional dynamical targets if needed.

The study establishes a systematic, ML-informed energy renormalization framework for DC-dependent coarse-graining of epoxy resins. For DGEBA+PACM and DGEBA+D400, the calibrated CG models match AA-derived density, Debye–Waller factor, Young’s modulus, and yield stress at room temperature across curing degrees with ~10% RMSPE, while reducing calibration complexity from 43 to 21 parameters. The approach yields interpretable sensitivity insights and 10^3-fold computational speedups, enabling simulations beyond the nanoscale and facilitating studies of heterogeneity, fracture, and composite behavior. Future directions include extending to multi-temperature calibration, incorporating bond-breaking for fracture and impact, exploring mixed hardener ratios more quantitatively, and expanding target sets to address longer-time dynamics in partially cured systems.

• Long-time dynamics: Despite matching picosecond ⟨u²⟩, the CG model exhibits slower long-time MSD than AA in uncrosslinked systems, suggesting that matching only short-time dynamics may not ensure full dynamical equivalence in chemically heterogeneous networks.
• AA crosslinking limits: AA networks did not reach full DC (max 90% for PACM, 95% for D400), constraining calibration at very high DC.
• Parameter non-uniqueness: Multiple functional parameterizations can achieve similar accuracy, especially for parameters with broad uncertainty bands; solutions depend on search space and training data.
• Trade-offs in multi-objective calibration: Perfectly matching one response (e.g., ⟨u²⟩ at high DC) can degrade others (e.g., modulus), requiring compromises.
• Scope: Calibration performed at fixed temperature (300 K) for two chemistries; transferability across temperature or to other chemistries would require additional data and calibration.
• Mapping choices: Although the protocol is likely robust, alternative CG mappings were not rigorously tested here.