Determination of interaction parameters in a bottom-up approach employed in reactive dissipative particle dynamics simulations for thermosetting polymers

Abstract

The limitations in previous dissipative particle dynamics (DPD) studies confined simulations to a narrow resin range. This study refines DPD parameter calculation methodology, extending its application to diverse polymer materials. Using a bottom-up approach with molecular dynamics (MD) simulations, we evaluated solubility parameters and bead number density governing nonbonded interactions via the Flory–Huggins parameter and covalent-bonded interactions. Two solubility parameter methods, Hildebrand and Krevelen–Hoftyzer, were compared for DPD simulations. The Hildebrand method, utilizing MD simulations, demonstrates higher consistency and broader applicability in determining solubility parameters for all DPD particles. The DPD/MD curing reaction process was examined in three epoxy systems: DGEBA/4,4′-DDS, DGEBA/MPDA and DGEBA/DETA. Calculations for the curing profile, gelation point, radial distribution function and branch ratio were performed. Compared to MD data for DGEBA/4,4′-DDS, the maximum deviation in secondary reactions between epoxy and amine groups according to DPD simulations with Krevelen–Hoftyzer was 14.8%, while with the Hildebrand method, it was 1.7%. The accuracy of the DPD curing reaction in reproducing the structural properties verifies its expanded application to general polymeric material simulations. The proposed curing DPD simulations, with a short run time and minimal computational resources, contributes to high-throughput screening for optimal resins and investigates mesoscopic inhomogeneous structures in large resin systems.

---

1. Introduction

The characteristics of carbon fiber reinforced plastics (CFRPs), especially their static strength and fatigue durability, strongly depend on those of the matrix resin.¹ Generally, thermosetting resin is utilized for the matrix resin in CFRPs, and the thermosetting resin consists of a basic resin and a curing agent. In experiments, the base resin and liquid curing agent are usually used; after a course of hours of cross-linking reactions, a cross-linked structure is formed. Besides the type of basic resin and curing agent, the curing reaction process will also largely determine the properties of the matrix resin. To select the appropriate resin and suitable reaction conditions to improve properties such as the toughness of CFRPs, a lot of experiments have been conducted in the past.^2–7 Maurice et al.² emphasized the influence of the conversion rate on material properties, particularly the tensile modulus. Zhang et al.³ investigated the impact of the heating rate on the curing behavior and phase separation in thermoplastic-toughened epoxy systems. Morancho et al.⁵ discussed the impact of excess base epoxy resin under non-stoichiometric conditions on the cross-linking process and the properties of the resulting materials.

Since there are numerous types of epoxy resin, the experimental selection of matrix resin is expensive and time-consuming. For this reason, molecular-scale curing simulations are expected to help reproduce and elucidate the thermoset resin reaction process. Molecular dynamics (MD) simulations reproduce molecular motions based on classical Newtonian mechanics. This has been proven to be an efficient way to perform curing reaction processes in various epoxy systems.^8–12 However, it also requires considerable computational resources and has high time costs due to the large number of atoms and complex force field calculations in the MD system. Molecular-scale curing simulations are required for high-throughput screening of resins to assist in experimental material exploration. Therefore, saving computational time is a crucial consideration. Furthermore, most of the mesoscale structure-forming phenomena that take place in polymeric materials influence the material properties, such as reaction-induced phase separation.¹³ However, such mesoscopic structures cannot be captured in MD simulations. Hence, to meet the demands of high-throughput screening and accurately replicate mesoscopic inhomogeneous structures, there is a pressing need for a curing simulation capable of faster calculations for larger systems.

The dissipative particle dynamics (DPD) simulation method has gained increasing popularity in polymer simulations in recent years. DPD treats multiple atoms as coarse-grained particles, referred to as beads. These beads interact with each other through various forces, typically including dissipative, random, and conservative forces. In DPD simulations, the internal degrees of freedom of the beads are neglected, and only the motion of the beads is considered. Extensive studies have demonstrated that the DPD simulation method is an efficient coarse-grained simulation technique capable of capturing longer time and length scales.^14–17 To enhance the accuracy of DPD simulations in reproducing the curing reaction process, the calculation methods and selection of interaction parameters play a crucial role in the overall DPD model. This paper focuses on precisely addressing these aspects, aiming to improve the fidelity of DPD simulations for capturing the intricate curing dynamics of epoxy materials.

DPD simulation is a simulation method aimed at enhancing computational efficiency and system scalability by consolidating atomic groups into coarse-grained particles. In 1992, Hoogerbrugge and Koelman introduced DPD as a computational approach to studying fluid flow and hydrodynamics,^18,19 with the theoretical foundation later established by Español and Warren.²⁰ Groot and Warren²¹ developed the conventional parameterization scheme for DPD interaction, linking the Flory–Huggins theory to interaction parameters. Building upon the aforementioned DPD force calculation method, Lísal et al. proposed a reaction ensemble DPD method to explore polydisperse homopolymer systems and supramolecular di-block copolymer melt.^22,23 While their focus was on replicating the reaction equilibrium of polymer systems rather than the polymerization process itself, their work illustrated the applicability of DPD for simulating polymerization. Liu et al. developed a curing reaction model using the same DPD method and applied it to surface-initiated polymerization.²⁴ However, in realistic systems, there are more than two kinds of beads, and the atomic groups represented in beads will influence their properties; specifically, the densities of different types of beads and the “bead–bead” interactions between different kinds of beads are variable.

There are several important approaches to parameterize DPD interactions for systems with variable bead–bead repulsions. Travis et al.²⁵ constructed bead–bead DPD interactions based on bead solubility parameters and defined pair–wise interactions using the Hildebrand regular solution theory. However, in their work, bead–bead interactions are computed independently without considering how to recover correct bead densities. Alternatively, Backer et al.²⁶ and Spaeth et al.²⁷ introduced beads with variable masses at different coarse graining levels. In their study, the system consisted of three regions: one comprising “normal” DPD particles, another comprising larger DPD particles, and a small region at the intersection where both types of particles coexist. Interactions between the beads, known as “bead–bead” interactions, vary depending on the composition of these two types of beads. However, the densities of the “normal” beads remain constant throughout the system. Kacar et al.¹⁵ introduced a density-dependent parameterization allowing the recovery of experimental bead densities from simulations. In their method, they kept most of the Groot–Warren framework intact, including the Flory–Huggins theory, but extended it by allowing bead-type dependent densities and like–like parameters. They combined this DPD simulation method with MD simulation in a realistic epoxy system, using MD simulation for parameterization.²⁸ In Kacar's framework, the calculation of solubility parameters utilized the Krevelen–Hoftyzer method,²⁹ relying on directly deriving molar contributions from functional groups. However, for certain functional groups, data in Krevelen's dataset is unavailable, as is the case with the sulfur-containing group in 4,4′-DDS. Therefore, there is a need to expand the methodology for calculating DPD parameters proposed by Kacar et al. and apply it to a wider range of polymer materials for research.

For the curing reaction model, Komarov et al.³⁰ proposed a coarse-grained curing molecular dynamics simulation to replicate the crosslinking process of thermoset resins. Li and Strachan³¹ applied this approach in MD simulation and used a novel method to update the charge after forming the crosslinked structure. Numerous studies have simulated curing processes using a distance-based criterion and constant reaction probability.^32,33 However, these methods are not very accurate in reproducing the reaction process, as the reaction probability should vary depending on the type of reaction substance and reaction conditions. Therefore, a new reaction model was proposed by Okabe et al., which was applied to MD simulation.⁸ There are two criteria for the reaction process: the distance between the reaction sites and the Arrhenius-type reaction probability determined from the activation energy and local temperature. Then, Kawagoe et al.³⁴ performed the curing reaction process by combining Okabe's curing reaction model with DPD simulation and introduced a reverse mapping procedure to reconstruct an all-atom MD system from a DPD system to verify the validity of the proposed curing DPD simulation and to evaluate its structural and thermomechanical properties. However, they employed the DPD scheme only for the single epoxy system given by Kacar et al.¹⁵ without discussing the calculation method of the interaction parameters. Consequently, to extend the application of the DPD method to a broader range of epoxy systems, it becomes essential to undertake a bottom-up exploration of DPD parameter calculation methods.

In order to effectively capture mesoscale structures and reduce computational costs, the DPD method has gained significant traction in the field of crosslinked polymers. Building upon previous research, the primary objective of this study is to develop a bottom-up approach for calculating DPD parameters, enabling their broader application in extensive epoxy systems. The DPD parameters are derived through MD simulations, primarily providing solubility parameters and bead densities, which are crucial for nonbonded interactions. Additionally, parameters related to covalent-bonded interactions, such as equilibrium distance and equilibrium angle, are also determined. To validate the accuracy of the developed DPD parametrization, structural properties are calculated in the on-reaction and post crosslinking process and compared against results obtained from MD simulations. By ensuring the fidelity of the DPD parametrization, this research aims to facilitate the widespread utilization of the DPD method in investigating complex epoxy systems and accurately reproducing their structural properties.

2. Simulation methods

2.1. DPD simulation

2.1.1. DPD theory. Different from the force field of Hoogerbrugge and Koelman,¹⁸ this study considers the bonded force between beads, which includes bond force and angle force as shown in eqn (1). Beads interact with a pair-wise conservative linear repulsive force. Additionally, a dashpot, named dissipative force, reduces the relative velocity between the particles, and a stochastic force, named random force, gives repulsions to the particles. Both of these forces are applied and are correlated with the fluctuation–dissipation theorem.Here, F^C_ij is the conservative force between beads i and j, which is represented by the following soft repulsion interaction:where a_ij is the repulsion parameter which is equal to e_ij/r_DPD. e_ij is the interaction parameter in units of energy. r_DPD is the cutoff distance. r̂_ij = r_ij/r_ij is the unit vector from beads j to i where r_ij = r_i − r_j and r_ij = |r_ij|. H(r) is the Heaviside step function. The dissipative force, F^D_ij, and the random force, F^R_ij, represent the viscous drag and thermal noise, respectively. The bond force, F^B_ij, and the angle force, F^A_ijk, represent the bonded force between beads. The calculation formula of interaction parameters e_ij¹⁵ is as follows:where ê_ij denotes the neutral repulsion parameter for every bead pair given byIn eqn (4), ρ^pure_i is the density of every bead, k is the Boltzmann constant, and T is the temperature. α is equal to 0.101 in this paper. p is the pressure, which is set to p = 40r_DPD⁻³kT.¹⁵χ_ij is the Flory–Huggins parameter, quantifying the excess repulsion between beads of different species. The calculation formula for Flory–Huggins interaction parameters χ_ij is as follows:where V_bead is the averaged molar volume of the system. δ_i and δ_j are the solubility parameters of beads i and j, respectively. The value of V_bead can be obtained from the following equation:where N_i is the number of beads in the simulation box for each type of bead, and ρ is the averaged molar density of the system. The cutoff distance r_DPD can be calculated using the following equation, taking into account a non-dimensional density of 3:²¹

r_DPD³ = 3ρ⁻¹. (8)

In this equation, δ can be evaluated directly from Hildebrand's theory,³⁵ referred to as the Hildebrand solubility parameter hereafter. In the present study, we also employed the Krevelen–Hoftyzer's²⁹ method, which offers an easier calculation based on the group contribution method to determine the solubility parameters. The two calculation methods for the solubility parameter mentioned above will be illustrated in the next section.

The remaining two forces F^D_ij and F^R_ij are given by

F^D_ij = −γw^D(r_ij)(r̂_ij·v_ij)r̂_ij, (9)

F^R_ij = σw^R(r_ij)θ_ijr̂_ij, (10)

where w^D(r_ij) and w^R(r_ij) are weight functions that vanish for r > r_c. γ is the friction coefficient, σ is the noise amplitude, v_ij = v_i − v_j is the relative velocity, and θ_ij is a randomly fluctuating variable with Gaussian statistics.

Español and Warren²⁰ demonstrated that the system samples the canonical ensemble and adheres to the fluctuation–dissipation theorem if the following relationships hold:

w^D(r_ij) = [w^R(r_ij)]², (11)

σ² = 2γkT, (12)

where T is the temperature and k is Boltzmann's constant. Typical expressions for w^D(r_ij) and w^R(r_ij) are chosen²¹ as

2.1.2. Parametrization of conservative forces. As described in the previous section, the calculation of conservative forces relies on the selection of solubility parameters. In this study, two methods were employed to calculate the solubility parameters: Krevelen–Hoftyzer²⁹ and Hildebrand³⁵ solubility parameters.

Travis et al.²⁵ constructed bead–bead DPD interactions using bead solubility parameters and defined pair-wise interactions based on Hildebrand's regular solution theory. In Hildebrand's theories, the solubility parameter can be obtained using the following equations:

where V is the molar volume and E_coh denotes the cohesive energy given by the equation as follows:

E_coh = ΔH_vap − kT, (15)

Here, ΔH_vap denotes the heat of vaporization given by the equation as follows:

ΔH_vap = −E_nonb + kT. (16)

Here, E_nonb is the intermolecular potential energy. From eqn (14)–(16), the calculation formula of solubility parameters is as follows:

The above describes the calculation method for the Hildebrand solubility parameter. In numerous DPD studies, the Krevelen–Hoftyzer solubility parameter is also commonly employed to calculate the conservative force.¹⁵ In the estimation of the solubility parameter, the cohesive energy is decomposed into three components: the contribution of dispersion forces, E_d, the contribution of polar forces, E_p, and the contribution of hydrogen bonding, E_h. The cohesive energy is calculated as follows:

E_coh = E_d + E_p + E_h. (18)

Accordingly, the calculation formula for the solubility parameter is as follows:

δ² = δ_d² + δ_p² + δ_h², (19)

where δ_d is the contribution of dispersion forces to the solubility parameter, δ_p is the contribution of polar forces to the solubility parameter, and δ_h is the contribution of hydrogen bonding to the solubility parameter.

The available experimental data, however, prove that it is impossible to derive a simple system for an accurate prediction of solubility parameter components from the chemical structure. In particular, the interaction of different structural groups in producing overall polar and hydrogen-bonding properties is so complicated that the solubility parameter components do not obey simple rules. Therefore, the approach of Hoftyzer and Van Krevelen²⁹ was used in this study as follows:

where F_di represents the contribution from the dispersion force. F_pi represents the contribution from the polarity force, and E_hi stands for the contribution from the hydrogen bond interaction energy. These three values can be obtained from a table provided by Hoftyzer and Van Krevelen.²⁹ Some of the table values used in this study are shown in Table 1.

Table 1 Group contributions to the solubility parameter component

Structural group	F di (MJ m^−3)^0.5 mol^−1	F pi (MJ m^−3)^0.5 mol^−1	E hi (J mol^−1)
–CH3	420	0	0
–CH2–	270	0	0
[double bond splayed left]CH–	80	0	0
–NH2	280	0	8400
–NH–	160	210	3100
[double bond splayed left]N–	20	800	5000
–OH	210	500	20 000
–O–	100	400	3000
–C6H5	1430	110	0

---

2.1.2. Parametrization of other forces. The calculation formula of random force becomes as follows:Here, ζ_ij is a Gaussian random number with zero mean and unit variance that is chosen independently for each pair of interacting particles, and Δt is the time step. r̂_ij = r_ij/r_ij is the unit vector from beads j to i where r_ij = r_i − r_j and r_ij = |r_ij|.

The friction coefficient γ of dissipative force is related to noise amplitude of random force as shown in eqn (23). Therefore, we only need to select an appropriate noise amplitude σ and time step Δt to determine the magnitude of dissipative force and random force.

The conversion from dimensionless to dimensional quantities in temperature in the DPD simulation utilizes the following formula:

where φ_DPD is the energy base value equal to kT. k is the Boltzmann constant. T is the real temperature value in the present simulation condition. T* is the dimensionless value of temperature which becomes 1.

The time step Δt in the simulations is selected as 0.02t_DPD. For our purpose the exact mapping of the time scale is not important since our main interest is the final crosslinked structure. Therefore, the unit DPD time, t_DPD, can be calculated using the following dimensional calculation.

where φ_DPD is the energy scale which equals kT. r_DPD is the length scales of DPD simulation as given below. Here, we determined r_DPD = 7.7 × 10⁻¹⁰ m first. m_estimated is the estimated DPD bead mass where m_estimated = 50 g mol⁻¹. Then, the value of t_DPD can be calculated as 3.465 ps which is similar with those employed by Kawagoe et al.³⁴

We chose the DPD length scale such that the overall dimensionless DPD density is approximately ρr_DPD³ = 3. This is the most used value in DPD simulations.²¹ To compute r_DPD, the ‘real’ density ρ of the system needs to be known. This value can be determined from the experimental value of the number density. Here, we estimate the average molecular volume by assuming that this volume is the weighted average of the pure species volumes. The volume is calculated as follows:

where N_i is the number of i beads, V_i is the actual volume of bead i, and ρ_i,pure is the number density of bead i. The volume of each bead will be obtained by MD simulation.

Referring to Groot,²¹ the dimensionless value of noise amplitude σ* was set to 3. The noise amplitude used in this study can be obtained by the following dimensional calculation.

After determining r_DPD, we can calculate the value of noise amplitude σ in the three systems which will greatly influence the dissipative force and random force.

For intramolecular interactions, F^B_ij and F^A_ijk are the two-body interaction of a bond and the three-body interaction of an angle, respectively.

V^B_ij = k^B_ij(r_ij − r₀)², (28)

V^A_ijk = k^A_ijk(θ_ijk − θ₀)², (29)

where V^B_ij and V^A_ijk are the bond and angle potentials, respectively. k^B_ij and k^A_ijk are the force constants. r₀ and θ₀ are the equilibrium distance and angle, respectively.

Kacar et al.¹⁵ explored the complete construction method of DPD interaction parameters, including the discussion of the above bond and angle force constants. They found that all bead sequences form highly stiff bonds and stiff angles with different spring constants for different bond and angle species. However, these values were impractically high to be used in the DPD simulations. Therefore, relatively larger force constants were employed for all bond and angle species in this study, as follows, referring to Kacar et al.¹⁵

k^B_ij = 500(kT/r_DPD²), (30)

k^A_ijk = 50kT. (31)

2.2 DPD modeling and MD simulation details for DPD parameters

In this study, DGEBA was used as the base resin, and 4,4′-DDS, MPDA, and DETA were employed as curing agents. Three epoxy systems were investigated: DGEBA/4,4′-DDS, DGEBA/MPDA, and DGEBA/DETA. In DPD simulations, the base resin and curing agents were coarsened into several beads corresponding to atomic groups in MD simulations. DGEBA was coarse-grained with five beads, denoted as “A, B, C, B, A” as shown in Fig. 1(a). Fig. 1(a) illustrates that the A bead contains the reactive epoxy group, the B bead contains a benzene ring with an oxygen atom bonded to it, while the C bead is composed solely of the hydrocarbon segment. Curing agents were coarse-grained with three beads: “D, E, D” beads for 4,4′-DDS, “F, G, F” beads for MPDA, and “H, I, H” beads for DETA as shown in Fig. 1(b). Fig. 1(b) demonstrates that the amino group, capable of participating in the curing reaction, is contained in the D, F and H beads, while the I bead also contains an amino group for potential reaction. In the three curing reaction systems, the A bead in DGEBA can undergo two crosslinking processes with the D beads of 4,4′-DDS, the F beads of MPDA, and the H beads of DETA. After the crosslinking process, new A′, D′, D′′, F′, F′′, H′, and H′′ beads are generated, as shown in Fig. 1(c). Additionally, the I bead in the middle of DETA can perform one crosslinking process with the A bead, resulting in the generation of a new I′ bead.

Fig. 1 Chemical structures of (a) base resin, (b) curing agents, and (c) products after reaction. 4,4′-DDS, MPDA and DETA are divided into DPD beads in a similar way. The two terminals are divided into one kind of beads for amino groups, represented by blue and purple circles, respectively, and the middle is divided into other kinds of beads, represented by orange and brown circles.

In this study, DPD parameters are obtained using MD simulation. This section introduces the simulation processes for the solubility parameter and number density of different types of beads, contributing to conservative force. Furthermore, the simulation process for the equilibrium distance and angle, which contribute to bonded force, is also described here.

First, the Winmostar³⁶ software was utilized to visually construct both the full atomistic models corresponding to DPD beads and bead sequences. Subsequently, the Gaussian 16 C.01³⁷ software was employed for the optimization of the structure (using the OPTIMIZE method) to obtain the molecular structure with the minimum energy. After breaking up an original molecule into fragments, the charges of the outermost atoms in the beads are retained as the charges of the original molecules. To ensure a net charge of 0, any excess charge is evenly distributed among the inner atoms, primarily the carbon atoms. To maintain the accuracy of simulating chemical reactions, the charges of oxygen and nitrogen atoms in the beads involved in the reaction groups are kept consistent with the charges of the original molecules. The assigned charges to each atom in the beads are shown in Fig. 2.

Fig. 2 Assignment of charges in (a) DGEBA particles, (b) 4,4′-DDS particles, (c) MPDA particles and (d) DETA particles. The charges of oxygen and nitrogen atoms involved in the reaction groups remain unchanged compared to the original molecule.

Due to the investigation of a large system and a shorter calculation time, the DREIDING force field³⁸ was adopted. In the following procedures, we utilized the Nosé–Hoover chains thermostat³⁹ to control the temperature and employed the Martyna–Tobias–Klein method⁴⁰ in the NPT ensemble. The cut-off distance of 12 Å was adopted. Periodic boundary conditions (PBCs) were applied in all three dimensions. The r-RESPA multiple-time-step method⁴¹ was employed, with 1 fs for intermolecular forces and 0.2 fs for intramolecular forces. The smooth particle mesh Ewald (SPME) method⁴² was utilized to calculate the long-range Coulomb force in the in-house code.

To calculate solubility parameters and number density, firstly, the initial simulated box size was set at 50.0 × 50.0 × 50.0 ų, with approximately 5000 atoms depending on the species of atomic groups. The temperature was set at 300 K, and the pressure at 1 atm, creating a liquid phase simulation system. Notably, C, H′′, I′, and F particles exist in the gaseous state under normal pressure, and to transition them to the liquid state, a minimum pressure was applied. Subsequently, 2000 step energy minimization using the steepest descent method was carried out. Following the energy minimization, a 1 ns NVT (constant number of particles, volume, and temperature) ensemble simulation and a 2 ns NPT (constant number of particles, pressure, and temperature) ensemble simulation were performed. The volume and nonbonded intermolecular energy could be directly obtained in the MD simulation.

To calculate the equilibrium distance and angle, we initially set the simulated initial box size at 230.0 × 230.0 × 230.0 ų, with 30 bead sequences. The temperature was set to 300 K, and the pressure was maintained at 1 atm, establishing the simulation system as a gas phase. Subsequently, we initiated the simulation process with 2000 step energy minimization using the steepest descent method. Following this, a 2 ns NVT ensemble simulation was conducted. The equilibrium distance and angle were then determined by calculating the center of mass of the beads. By obtaining the probability distribution of bonds and angles, we traced the peak values, representing the most probable values, to ascertain the equilibrium distance and angle. These equilibrium values are utilized for initial state setup and in calculations of angle and bond potentials as shown in the formulas of eqn (28) and (29).

2.3 Curing DPD/MD simulation details

In the DPD/MD simulation system, the molecular numbers of the basic resin and curing agents were determined by the stoichiometric ratio. For DGEBA/4,4′-DDS, DGEBA/MPDA, and DGEBA/DETA, they were 2 : 1, 2 : 1, and 5 : 2, respectively. The box size of the DPD simulation was the same as that of the uncured MD system, which was averaged in the NPT simulation for 1 atm and 453 K. The number of molecules, average volumes of beads, r_DPD and noise level σ in the three epoxy systems are shown in Table 2. Due to the different volume acquisition methods used in the two methods, the average volume and r_DPD will also be different. The noise level is a crucial parameter in calculating the random force and dissipative force within the DPD forces. The calculation method can be found in the ESI.† Periodic boundary conditions were applied in all directions. The time step of the DPD simulation was 69.3 fs. The system temperature was controlled at 453 K.

Table 2 The average volumes of beads, r_DPD and noise level σ in the three epoxy system

Epoxy system	Number of molecules	Average volumes (Hildebrand/Krevelen) (Å^3)	r DPD (Hildebrand/Krevelen) (Å)	Noise level σ (Hildebrand/Krevelen) (10^−17/J s^0.5 m^−1)
DGEBA/4,4′-DDS	600/300	167.2/126.9	7.951/7.243	4.547/4.760
DGEBA/MPDA	600/240	152.2/123.6	7.702/7.178	4.618/4.781
DGEBA/DETA	600/300	164.9/120.5	7.902/7.121	4.557/4.801

---

To obtain the initial non-crosslinked structure, a relaxation step is necessary before the curing reaction. In the DPD simulation system, a DPD simulation of 0.693 ns which is equal to 10 000 DPD simulation steps at 300 K was conducted prior to crosslinking. In the MD system, a series of simulations, including 2 ns of NPT (300 K, 1 atm) and 0.5 ns of NVT (300 K), were performed before the crosslinking process. Both DPD and MD simulations in this study underwent sufficient mixing and achieved equilibrium states before crosslinking. During the crosslinking process, the distances between all epoxy groups and amino groups in the system need to be calculated. Based on these distances, a determination is made regarding whether a chemical reaction occurs between all pairs of functional groups within a threshold value, denoted as R_c. The value of R_c is set to 5.64 Å, four times the equilibrium C–N bond length,³¹ in both the MD and DPD curing reaction processes.³⁴ The settings for equilibrium angles and equilibrium distances used in this study are based on calculations from an all-atom MD model. Therefore, the corresponding value of R_c remains consistent. In the reaction determination, the temperature around the functional group should be calculated first. Fig. 3 is a detailed flow chart of the curing reaction process. In this procedure, the reaction probability was calculated based on the following Arrhenius equation.⁴³

where E_a is the activation energy, R is the gas constant, and T represents the local temperature, which is calculated from neighboring atoms between an epoxy group and an amino group. This range of neighboring atoms is nearly equivalent to the coarse-graining approach used in DPD simulations. The fluctuation of temperature comes from the small number of atoms or beads to be engaged in the temperature evaluation. But at least for reaction cycles the acceptance of the reaction takes place when the temperature occasionally gets higher relative to the activation energy. This will give the kinetic process which realize the reaction process in our model, even though there is an arbitrariness for the definition of local temperature. Here, the setting of the acceleration constant A is to bridge the timescale difference between the experiment and MD simulation. Therefore, based on previous research validation,⁴⁴ the acceleration constant A was set as 10¹⁹ for all reaction pairs to ensure that the reactions in the three systems can be reproduced within a reasonable computational time frame. The number of reaction sites are equivalent between MD and DPD simulation, whether atoms are coarse-grained or not, so we used the same pre-exponential factor for both systems.

Fig. 3 Flow chart of the curing reaction process.

To determine the reaction probability in the reaction calculation, the Monte Carlo method⁴⁵ was used here. In the reaction judgement, the reaction probability is compared with a random number in the range of 0 to 1. If the reaction probability is greater than the number, a new bond is formed.

For the MD system, the following is the equation for calculating the kinetic energy after reaction.

K_after = K_before + H_f, (33)

where K_after represents the kinetic energy after the reaction, K_before indicates the kinetic energy before the reaction and H_f stands for the heat of formation.

In the DPD system, an adjustment factor β was added such that the localized temperature increment in a coarse-grained system due to H_f is almost identical to that of an MD system. The following is the equation for calculating the kinetic energy after the reaction.

K_after = K_before + βH_f, (34)

Here, β is the adjustment factor such that the localized temperature increment in a coarse-grained system due to H_f is almost identical to that of an MD system. β was set as 0.069 referring to Kawagoe et al.³⁴ The H_f values for all three systems generally range from 20 to 25. The value of H_f fluctuates very little, and it has been verified that the coefficient can control the increase in temperature of the three systems.³⁴ Therefore, in this research we consider using the same coefficient.

The activation energy and the heat of formation determined using formulas from ref. 32–34 have been referenced from data in the previous publication.⁴⁴ These values were the activation energy and heat of formation precisely calculated using the global reaction route mapping (GRRM) method⁴⁶ on the basis of ab initio quantum chemical (QC) calculations.

Adding H_f induces chain reactions as T increases, i.e., k increases in eqn (32). During the relaxation phase of the simulation, multiple-step relaxation cycles were conducted to relax and equilibrate the new topology after cross-linking, and preparing for the next cross-linking step. For the MD curing reaction relaxation, the DREIDING force field was selected. Structural relaxation was carried out using the steepest descent method for 1000 steps. Subsequently, a 15 ps NPT (453 K, 1 atm) ensemble simulation and a 1 ps NVT (453 K) ensemble simulation were performed. In the relaxation process following DPD crosslinking, a 2.5 ps DPD simulation was executed. The crosslinking structure was formed using an external curing code written in the Python language, while the relaxation calculations were conducted using the in-house DPD/MD package. Following all cross-linking cycles, we conduct a relaxation process lasting 10 ns NPT (300 K, 1 atm) ensemble to obtain a stable structure.

3. Results and discussion

3.1 DPD interaction parameters

The number density of beads was obtained through MD simulation. One molecule is divided into several groups corresponding to the DPD particles, and this division results in the formation of broken covalent bonds. The broken covalent bond was attached to hydrogen atoms to form a closed-shell molecular structure. Subsequently, we conducted MD simulation to obtain the density of beads, including the newly formed beads.

Tables 3 and 4 present the number density of various bead types. The number density of newly generated beads is smaller than that of the corresponding original beads, indicating that the volume of the newly generated beads is larger. For example, the volumes of beads D′ and D′′ are larger than that of bead D. This expansion in volume is attributed to the crosslinking process, where carbon atoms are added, leading to a reduction in hydrogen atoms in beads D, F, and H. The relatively larger volume of carbon atoms contributed to the overall volume increase. The calculation of the Hildebrand solubility parameter involves nonbonded energy and system volume. Conversely, for the Krevelen–Hoftyzer solubility parameter, the volume was sourced from the SciFinder-n database,⁴⁷ which is smaller compared to the volume obtained from MD simulations for the Hildebrand solubility parameter. This may be attributed to the Dreiding force field, which often leads to a larger volume due to the selection of LJ parameters. The LJ parameters in DREIDING were mainly derived from Buckingham potential (exponential-6, X6) parameters but without direct derivation from physical properties like density. The three components of the molar attraction function were determined by referencing the table provided by Hoftyzer and Van Krevelen.²⁹ This study explores the differences between Krevelen–Hoftyzer and Hildebrand solubility parameters.

Table 3 The number density of DPD beads obtained from MD simulation for the Hildebrand solubility parameter. The table displays the number density of DGEBA and curing agent beads, along with the number density of the new beads produced after the DPD curing reaction

Bead type (DGEBA)	Number density (Å^−3)	Bead type (DDS)	Number density (Å^−3)	Bead type (MPDA)	Number density (Å^−3)	Bead type (DETA)	Number density (Å^−3)
A	0.0061	D	0.0057	F	0.0175	H	0.0078
B	0.0058	E	0.0126	G	0.0065	I	0.0065
C	0.0056	D′	0.0048	F′	0.0078	H′	0.0065
A′	0.0067	D′′	0.0044	F′′	0.0065	H′′	0.0063
						I′	0.0063

---

Table 4 The number density of DPD beads for the Krevelen–Hoftyzer solubility parameter⁴⁷

Bead type (DGEBA)	Number density (Å^−3)	Bead type (DDS)	Number density (Å^−3)	Bead type (MPDA)	Number density (Å^−3)	Bead type (DETA)	Number density (Å^−3)
A	0.0088	D	0.0066	F	0.0201	H	0.0128
B	0.0069	E	0.0139	G	0.0073	I	0.0085
C	0.0083	D′	0.0055	F′	0.0128	H′	0.0085
A′	0.0092	D′′	0.0048	F′′	0.0085	H′′	0.0071
						I′	0.0071

---

The evaluated solubility parameters are presented in Table 5. It is observed that the Hildebrand solubility parameters in this study are consistently smaller than the Krevelen–Hoftyzer solubility parameters. This deviation is partly attributed to the difference in volume of beads between the two methods. The volume calculated using the Hildebrand method is generally larger than that calculated using the Krevelen–Hoftyzer method, as indicated in Tables 3 and 4, resulting in a smaller calculated solubility parameter. Based on the DPD results compared with MD simulations, the Hildebrand method demonstrates a higher level of consistency. Furthermore, due to the limited data available for the Krevelen method, it cannot comprehensively calculate the solubility parameters for all DPD particles. In contrast, the Hildebrand method theoretically has the potential to calculate the solubility parameters for all DPD particles, making it more promising for widespread use. It is worth emphasizing this point in Tables 3–5: H and F′, H′ and F′′, and I, H′′, and I′ are the same type of particles with identical compositions, only differing in naming convention. Therefore, the volumes and solubility parameters of these particles are the same.

Table 5 The solubility parameter of DPD beads comprises two types: the Hildebrand solubility parameter and the Krevelen–Hoftyzer solubility parameter. This table documents the solubility parameters for DGEBA and the curing agent, encompassing parameters for newly generated beads in the DPD curing process. The Krevelen–Hoftyzer parameter data inside the parenthesis is substituted with the Hildebrand solubility parameter due to an inability to query the decomposed energy

Bead type (DGEBA, DETA)	Hildebrand solubility parameter ((J cm^3)^0.5)	Krevelen–Hoftyzer solubility parameter ((J cm^3)^0.5)	Bead type (4,4′-DDS, MPDA)	Hildebrand solubility parameter ((J cm^3)^0.5)	Krevelen–Hoftyzer solubility parameter^15 ((J cm^3)^0.5)
A	15.34	16.44	H	18.61	19.99
B	18.58	24.54	I	14.37	15.87
C	9.37	12.35	H′	14.37	15.87
A′	13.19	23.54	H′′	18.38	19.38
D	21.01	18.51	I′	18.38	19.38
E	25.35	(25.35)	F	10.97	(10.97)
D′	18.96	18.22	G	22.34	(22.34)
D′′	18.78	19.95	F′	18.61	19.99
			F′′	14.37	15.87

---

The equilibrium distance is a parameter used to calculate the bond force between two beads, while the equilibrium angle is used to calculate the potential for angle bending between three atoms. In this study, the equilibrium distances and angles between different bead sequences were obtained from MD simulations, as presented in Tables 6 and 7, respectively. The equilibrium distance and angle curves are provided in the ESI.† A total of 300 sample data points were obtained from the MD simulation, and the equilibrium distance was determined as the peak value of the maximum occurrence probability. These equilibrium values are utilized for the initial state setup and in calculations of angle and bond potentials as shown in the formulas of eqn (28) and (29).

Table 6 This table displays the equilibrium distances between beads, the equilibrium distances between the curing agent beads and the equilibrium distances of the new sequences generated after crosslinking. The term “new sequences” here pertains to the generated new sequences, such as A–D and A–F. The equilibrium distance remains consistent between A′–D and A–D before and after the reaction

Bead sequence (DGEBA)	Equilibrium distance (Å)	Bead sequence (curing agent)	Equilibrium distance (Å)	Bead sequence (between DGEBA and curing agent)	Equilibrium distance (Å)
A–B	4.713	D–E	3.725	A–D	4.841
B–C	3.883	F–G	2.809	A–F	3.885
		H–I	3.302	A–H	4.247

---

Table 7 This table displays the equilibrium angle between beads. The term “new sequences” here pertains to the generated new sequences, such as A′–B–C and D′–E–D. The equilibrium distance remains consistent between A′–B–C and A–B–C before and after the reaction

Bead sequence (DGEBA)	Equilibrium angle (Å)	Bead sequence (curing agent)	Equilibrium angle (Å)
A–B–C	112.4	D–E–D	153.2
B–C–B	143.5	F–G–F	178.2
		H–I–H	154.8

---

3.2 Curing process

In the curing reaction process, the ratios of unreacted beads, primary reaction beads, secondary reaction beads and middle reaction beads in DETA at different curing rates were monitored. Fig. 4 shows the ratio of unreacted beads, primary reaction beads, secondary reaction beads and middle reaction beads as a function of curing rate for the three systems.

Fig. 4 Reacted/unreacted bead fractions as a function of curing rate in the (a) DGEBA/4,4′-DDS, (b) DGEBA/MPDA, and (c) DGEBA/DETA. Grey lines represent the ratio of unreacted beads such as D in the DDS system, blue lines represent the ratio of beads via the 1st reaction such as D1, and green lines represent the ratio of beads via the 2nd reaction such as D2. Circles represent MD simulations, blocks represent DPD simulations (Hildebrand), and triangles represent DPD simulations (Krevelen–Hoftyzer).

For the total simulation duration, DPD simulations require less time, primarily because they involve fewer simulation steps compared to MD simulations and fewer particles in the system. The MD simulation ran for approximately 50 hours on 200 cores, whereas the DPD simulation ran for only about 20 minutes on 20 cores, which is roughly 1/150 of the MD simulation.

The calculation formula for the deviation value is |N_DPD − N_MD|/N_MD where N_DPD represents the number of reactions obtained by a DPD curing reaction, and N_MD represents the number of reactions obtained by an MD curing reaction. For the DGEBA/4,4′-DDS system, the DPD curing reaction profiles with Hildebrand parameters are more consistent with the MD curing reaction curve. In the DPD curing process with the Hildebrand solubility parameter, there is excellent agreement with MD simulations, in terms of both primary and secondary reactions. However, in the DPD curing process with the Krevelen–Hoftyzer solubility parameter, the ratio of primary reaction beads consistently falls below that of the MD curing reaction and the DPD curing process with the Hildebrand solubility parameter. Conversely, the ratio of secondary reaction beads is consistently higher. At the maximum deviation point, occurring near a 40% curing rate for the DPD curing reaction (Hildebrand) and MD curing reaction, the deviation of the primary reaction is 2.2%, and the deviation of the secondary reaction is 1.7%. In contrast, for the DPD curing reaction with Krevelen–Hoftyzer parameters, the maximum deviation point occurs near a 50% curing rate, with a deviation of the primary reaction of 25.8% and a deviation of the secondary reaction of 14.8%.

For the DGEBA/MPDA system, both types of DPD curing reactions yield similar results, although they deviate significantly from the MD curing reaction. This discrepancy can be attributed to the notably high size ratio of DPD particles, which influences the accuracy of the DPD simulation. For instance, considering the reactive particle F in MPDA, it comprises only four atoms, leading to a substantial size difference between bead F and the B particle. In the DGEBA/DETA system, the profiles of DPD curing reactions with the Hildebrand solubility parameters closely resemble those of the MD curing reaction compared to the Krevelen–Hoftyzer approach. Considering all three systems, the results obtained from the Hildebrand method are closer to those of MD simulations. Hildebrand solubility parameters are parameterized based on energy and volume from MD simulations using the Dreiding force field, which is the same force field applied in the MD reaction simulations of this study. The correspondence between MD and DPD (Krevelen) methods is not good, because, in this case, solubility is parameterized based on literature values rather than MD simulations with the Dreiding force field.

To elucidate the impact of variations in the curing reaction process on physical properties, the gelation point was examined. Criteria for identifying the gelation point included the largest molecular weight (LMW) and the second largest molecular weight (SMW). In this study, a more analytical mass-based criterion, denoted as the reduced molecular weight (RMW), was employed. The RMW is defined as the molecular weight average of all the molecules in the system. It increases until the largest molecule starts to dominate, after which it declines.

The gelation point, corresponding to the maximum point of the RMW profile concerning the curing rate, is utilized to assess the gelation process. From Fig. 5, gelation points for different systems can be determined. For the DPD curing reaction with the Krevelen–Hoftyzer solubility parameter, the gelation points are as follows: DGEBA/4,4′-DDS: 60%, DGEBA/MPDA: 65%, DGEBA/DETA: 55%. In the case of the DPD curing reaction with the Hildebrand solubility parameter, the gelation points are as follows: DGEBA/4,4′-DDS: 55%, DGEBA/MPDA: 65%, DGEBA/DETA: 50%. For the MD curing reaction, the gelation points are observed as follows: DGEBA/4,4′-DDS: 55%, DGEBA/MPDA: 60%, DGEBA/DETA: 50%. In the DGEBA/4,4′-DDS system, the gelation point obtained from the DPD curing reaction with the Hildebrand solubility parameter aligns well with that obtained from the MD curing reaction. Furthermore, it also shows good agreement with the theoretical estimation by M. Livraghi et al., suggesting a gelation point of about 57 ± 6% based on percolation theory.⁴⁸

Fig. 5 RMW, SMW LMW curve in (a) DGEBA/4,4′-DDS, (b) DGEBA/MPDA, and (c) DGEBA/DETA. Blue lines represent LMW curves, green lines represent SMW curves, and red lines represent RMW curves.

The ratio of branches at the molecular scale represents the internal molecular connectivity of the crosslinking network structure within a material. This ratio is determined by the number of monomers with each branching structure, divided by the total number of monomers. It has implications for thermophysical properties such as glass transition temperature and Young's modulus. In a crosslinked system, each basic resin can undergo a maximum of two reactions, and each curing agent can undergo a maximum of four reactions. Consequently, resin/curing agent molecules will form varying numbers of branches, including 0, 1, 2, 3, or 4 branches. The ratios of branches can be compared between DPD simulation and MD simulation. In the DPD system, the reaction times of A, D, F, and H beads in a basic resin/curing agent molecule are calculated, while in the MD system, the reaction times of the amino group and epoxy group in a basic resin/curing agent molecule are counted. The number of branches changes continuously during the curing reaction. Before the gel point of 55%, branch 1 and branch 0 dominate the overall crosslinking network. From Fig. 6(a), it can be observed that branch 1 is located at the end of the crosslinked polymer chain. Therefore, during the period from the onset of cross-linking to the gel point, the increase in the proportion of branch 1 divides the system into smaller blocks of polymer. After exceeding approximately 40% conversion, the proportion of branch 1 starts to decrease, indicating that the dispersed polymer blocks begin to coalesce into larger polymer structures. Until reaching the gel point of 55%, the largest weight polymers begin to form and continue to expand with increasing conversion.

Fig. 6 Branch ratios depending on the curing conversion for the DGEBA/4,4′-DDS system. (a) Calculation diagram of the branch quantity. The colors correspond to the number of branches: gray for a zero-branching structure, yellow for a one-branching structure, orange for a two-branching structure, green for a three-branching structure, and purple for a four-branching structure. (b) Gray, red, green, blue and purple curves represent the ratio of 0, 1, 2, 3 and 4 branches, respectively. For both figures, circles represent MD simulations, the blocks represent DPD simulations (Hildebrand), and triangles represent DPD simulations (Krevelen–Hoftyzer).

In the DGEBA/4,4′-DDS system, the ratio of branches in the DPD curing process, using the Hildebrand solubility parameter, closely aligns with that observed in the MD curing reaction process. This similarity arises from the comparable ratios of primary and secondary reactions observed in the MD system, as compared to the DPD curing process using the Hildebrand solubility parameter. As a result, the branch structures remain relatively similar during the curing reaction process. However, in the DPD curing process with the Krevelen–Hoftyzer solubility parameter, the growth rate of primary reactions increases significantly after reaching a 20% curing rate. Consequently, during this stage, secondary reactions occur slowly following the primary reactions. For example, the generation of 1 branch from the primary reaction quickly occurs, but the change to 2 branches is slow. As a result, the curve representing the presence of 1 branch reaches a peak in the DPD curing process, whereas the curve for 1 branch in the MD curing process represents a plateau.

3.2 Polymer structure after crosslinking

After conducting DPD and MD curing reactions, two types of coarse-grained structures were obtained using DPD methods with the Krevelen–Hoftyzer solubility parameter and the Hildebrand solubility parameter. Throughout the curing process, the formation of new bonds between beads led to changes in inter-bead distances, impacting the structural properties of the system. In contrast, the crosslinking structure in MD simulations is represented using full atomistic models. To compare structures between the DPD curing process and the MD curing process, it is essential to calculate the center of mass for atomic groups corresponding to each DPD bead in the MD system. The radial distribution function (RDF) results are presented in Fig. 7 and 8 for the DGEBA/4,4′-DDS and DGEBA/DETA systems, respectively. The distance between the epoxy group and SO₂ varies depending on whether the cross-linking point is a secondary amine or a tertiary amine, leading to two peaks in A′–E. However, in DPD simulations, the softness of DPD interactions results in RDFs showing clearer exclusion volumes compared to MD simulations. Moreover, DPD simulations lack the ability to capture forces like hydrogen bonding in full atomic detail, contributing to differences in the RDF curves. The RDF curves obtained from the Krevelen–Hoftyzer solubility parameter and the Hildebrand solubility parameter are similar, but some deviations are observed in the curves for C–D′′, C–E, C–H′′, and C–I′ pairs. Overall, the RDF curves in both types of DPD systems show good agreement with those obtained from the MD system.

Fig. 7 RDF between each group pair in the DGEBA/4,4′-DDS system. The red line represents the MD simulation, the blue dashed line represents the DPD (Hildebrand) simulation, and the green dashed line represents the DPD (Krevelen–Hoftyzer) simulation.

Fig. 8 RDF between each group pair in the DGEBA/DETA system. The red line represents the MD simulation, the blue dashed line represents the DPD (Hildebrand) simulation, and the green dashed line represents the DPD (Krevelen–Hoftyzer) simulation.

4. Conclusions

In this paper, we explored DPD interaction parameters using a bottom-up approach based on MD simulations and analyzed the accuracy of DPD simulations in terms of reproducing the curing process and structural properties of epoxy crosslinking polymers. The density and solubility parameters of the coarse-grained beads were evaluated by MD simulations, which were then used to calculate interaction parameters. Moreover, the equilibrium bond length and bond angle between covalently bonded DPD beads were also determined by MD simulations. Herein, two types of solubility parameters, i.e., the Hildebrand solubility parameter and the Krevelen–Hoftyzer solubility parameter, were calculated in this study.

Afterwards, the DPD/MD curing reaction process was conducted in three epoxy systems: DGEBA/4,4′-DDS, DGEBA/MPDA and DGEBA/DETA. Subsequently, we compared the number of reactions for different types. We observed that the DPD simulation with the Hildebrand solubility parameter is more consistent with MD simulations compared to the DPD simulation with the Krevelen–Hoftyzer solubility parameter. The Curing properties in the reaction and the structure properties were compared between these methods. As a result, the DPD simulation can reproduce the structure of the full atomistic model to a certain degree. From a comprehensive perspective of the above structural parameters, the Hildebrand parameter is more suitable for reproducing the curing process in MD simulations. In comparison to the Krevelen–Hoftyzer method, the Hildebrand method, based on comprehensive simulations to calculate solubility parameters, theoretically enables the calculation of solubility parameters for any type of DPD particles. This conclusion significantly expands the applicability of DPD simulations for polymeric materials.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors gratefully acknowledge the financial support from the Council for Science, Technology, and Innovation (CSTI), Cross-ministerial Strategic Innovation Promotion Program (SIP) (Funding agency: JST). Numerical simulations were performed on the Supercomputer system “AFI-NITY” at the Advanced Fluid Information Research Center, Institute of Fluid Science, Tohoku University. This work was supported by JST, the establishment of university fellowships towards the Creation of Science Technology Innovation, Grant Number JPMJFS2102.

References

1. D. D. R. Cartié and P. E. Irving, Effect of resin and fibre properties on impact and compression after impact performance of CFRP, Composites, Part A, 2002, 33(4), 483–493,  DOI:10.1016/S1359-835X(01)00141-5.
2. M. J. Marks and R. V. Snelgrove, Effect of conversion on the structure-property relationships of amine-cured epoxy thermosets, ACS Appl. Mater. Interfaces, 2009, 1(4), 921–926,  DOI:10.1021/am900030u.
3. J. Zhang, Q. Guo and B. L. Fox, Study on thermoplastic-modified multifunctional epoxies: influence of heating rate on cure behaviour and phase separation, Compos. Sci. Technol., 2009, 69(7–8), 1172–1179,  DOI:10.1016/j.compscitech.2009.02.016.
4. Z. Yang, H. Peng, W. Wang and T. Liu, Crystallization behavior of poly(ε-caprolactone)/layered double hydroxide nanocomposites, J. Appl. Polym. Sci., 2010, 116(5), 2658–2667,  DOI:10.1002/app.31787.
5. J. M. Morancho, X. Ramis, X. Fernández-Francos, J. M. Salla, A. O. Konuray and À. Serra, Curing of off-stoichiometric amine–epoxy thermosets, J. Therm. Anal. Calorim., 2018, 133(1), 519–527,  DOI:10.1007/s10973-018-7158-2.
6. J. Hu, J. Shan, D. Wen, X. Liu, J. Zhao and Z. Tong, Flame retardant, mechanical properties and curing kinetics of DOPO-based epoxy resins, Polym. Degrad. Stab., 2014, 109, 218–225,  DOI:10.1016/j.polymdegradstab.2014.07.026.
7. N. Odagiri, K. Shirasu, Y. Kawagoe, G. Kikugawa, Y. Oya and N. Kishimoto, et al., Amine/epoxy stoichiometric ratio dependence of crosslinked structure and ductility in amine-cured epoxy thermosetting resins, J. Appl. Polym. Sci., 2021, 138(23), 50542,  DOI:10.1002/app.50542.
8. T. Okabe, T. Takehara, K. Inose, N. Hirano, M. Nishikawa and T. Uehara, Curing reaction of epoxy resin composed of mixed base resin and curing agent: Experiments and molecular simulation, Polymer, 2013, 54(17), 4660–4668,  DOI:10.1016/j.polymer.2013.06.026.
9. H. B. Fan and M. M. F. Yuen, Material properties of the cross-linked epoxy resin compound predicted by molecular dynamics simulation, Polymer, 2007, 48(7), 2174–2178,  DOI:10.1016/j.polymer.2007.02.007.
10. T. Okabe, Y. Oya, K. Tanabe, G. Kikugawa and K. Yoshioka, Molecular dynamics simulation of crosslinked epoxy resins: Curing and mechanical properties, Eur. Polym. J., 2016, 80, 78–88,  DOI:10.1016/j.eurpolymj.2016.04.019.
11. V. Varshney, S. S. Patnaik, A. K. Roy and B. L. Farmer, A molecular dynamics study of epoxy-based networks: cross-linking procedure and prediction of molecular and material properties, Macromolecules, 2008, 41(18), 6837–6842,  DOI:10.1021/ma801153e.
12. Y. Kawagoe, K. Kawai, Y. Kumagai, K. Shirasu, G. Kikugawa and T. Okabe, Multiscale modeling of process-induced residual deformation on carbon-fiber-reinforced plastic laminate from quantum calculation to laminate scale finite-element analysis, Mech Mater., 2022, 170, 104332,  DOI:10.1016/j.mechmat.2022.104332.
13. H. Kishi, Y. Kunimitsu, Y. Nakashima, T. Abe, J. Imade and S. Oshita, et al., Control of nanostructures generated in epoxy matrices blended with PMMA-b-PnBA-b-PMMA triblock copolymers, eXPRESS Polym. Lett., 2015, 9(1), 23–35,  DOI:10.3144/expresspolymlett.2015.4.
14. H. Liu, H. Li and Z. Y. Lu, Incorporating chemical reactions in dissipative particle dynamics simulations, Proc. Comput. Sci., 2011, 4, 1021–1028,  DOI:10.1016/j.procs.2011.04.108.
15. G. Kacar, E. A. J. F. Peters and G. De With, Mesoscopic simulations for the molecular and network structure of a thermoset polymer, Soft Matter, 2013, 9(24), 5785–5793,  10.1039/C3SM50304F.
16. S. Thomas, M. Alberts, M. M. Henry, C. E. Estridge and E. Jankowski, Routine million-particle simulations of epoxy curing with dissipative particle dynamics, J. Theor. Comput. Chem., 2018, 17(3), 1–21.
17. H. Liu, Y. L. Zhu, Z. Y. Lu and F. Müller-Plathe, A kinetic chain growth algorithm in coarse-grained simulations, J. Comput. Chem., 2016, 37(30), 2634–2646,  DOI:10.1142/S0219633618400047.
18. P. J. Hoogerbrugge and J. M. V. A. Koelman, Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics, Europhys. Lett., 1992, 19(3), 155–160,  DOI:10.1209/0295-5075/19/3/001.
19. J. M. V. A. Koelman and P. J. Hoogerbrugge, Dynamic simulations of hard-sphere suspensions under steady shear, Europhys. Lett., 1993, 21(3), 363–368,  DOI:10.1209/0295-5075/21/3/018.
20. P. Espanol and P. Warren, Statistical mechanics of dissipative particle dynamics, Europhys. Lett., 1995, 30(4), 191–196,  DOI:10.1209/0295-5075/30/4/001.
21. R. D. Groot and P. B. Warren, Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation, J. Chem. Phys., 1997, 107(11), 4423–4435,  DOI:10.1063/1.474784.
22. M. Lísal, J. K. Brennan and W. R. Smith, Mesoscale simulation of polymer reaction equilibrium: Combining dissipative particle dynamics with reaction ensemble Monte Carlo. I. Polydispersed polymer systems, J. Chem. Phys., 2006, 125(16), 164905,  DOI:10.1063/1.2359441.
23. M. Lísal, J. K. Brennan and W. R. Smith, Mesoscale simulation of polymer reaction equilibrium: Combining dissipative particle dynamics with reaction ensemble Monte Carlo. II. Supramolecular diblock copolymers, J. Chem. Phys., 2009, 130(10), 104902,  DOI:10.1063/1.3079139.
24. H. Liu, M. Li, Z. Y. Lu, Z. G. Zhang, C. C. Sun and T. Cui, Multiscale simulation study on the curing reaction and the network structure in a typical epoxy system, Macromolecules, 2011, 44(21), 8650–8660,  DOI:10.1021/ma201390k.
25. K. P. Travis, M. Bankhead, K. Good and S. L. Owens, New parametrization method for dissipative particle dynamics, J. Chem. Phys., 2007, 127(1), 014109,  DOI:10.1063/1.2746325.
26. J. A. Backer, C. P. Lowe, H. C. J. Hoefsloot and P. D. Iedema, Combined length scales in dissipative particle dynamics, J. Chem. Phys., 2005, 123(11), 114905,  DOI:10.1063/1.2013208.
27. J. R. Spaeth, T. Dale, I. G. Kevrekidis and A. Z. Panagiotopoulos, Coarse-graining of chain models in dissipative particle dynamics simulations, Ind. Eng. Chem. Res., 2011, 50(1), 69–77,  DOI:10.1021/ie100337r.
28. G. Kacar, E. A. J. F. Peters and G. De With, Multi-scale simulations for predicting material properties of a cross-linked polymer, Comput. Mater. Sci., 2015, 102, 68–77,  DOI:10.1016/j.commatsci.2015.02.021.
29. D. W. Van Krevelen, Properties of Polymers: Their Correlation with Chemical Structure; their Numerical Estimation and Prediction from Additive Group Contributions. Properties of Polymers: Their Correlation with Chemical Structure; their Numerical Estimation and Prediction from Additive Group Contributions, Elsevier Science, 4th edn, 2009, pp. 1–1004.
30. P. V. Komarov, Y. T. Chiu, S. M. Chen, P. G. Khalatur and P. Reineker, Highly cross-linked epoxy resins: An atomistic molecular dynamics simulation combined with a mapping/reverse mapping procedure, Macromolecules., 2007, 40(22), 8104–8113,  DOI:10.1021/ma070702.
31. C. Li and A. Strachan, Molecular simulations of crosslinking process of thermosetting polymers, Polymer, 2010, 51(25), 6058–6070,  DOI:10.1016/j.polymer.2010.10.033.
32. C. Li, G. A. Medvedev, E. W. Lee, J. Kim, J. M. Caruthers and A. Strachan, Molecular dynamics simulations and experimental studies of the thermomechanical response of an epoxy thermoset polymer, Polymer, 2012, 53(19), 4222–4230,  DOI:10.1016/j.polymer.2012.07.026.
33. J. R. Gissinger, B. D. Jensen and K. E. Wise, Modeling chemical reactions in classical molecular dynamics simulations, Polymer, 2017, 128, 211–217,  DOI:10.1016/j.polymer.2017.09.038.
34. Y. Kawagoe, G. Kikugawa, K. Shirasu and T. Okabe, Thermoset resin curing simulation using quantum-chemical reaction path calculation and dissipative particle dynamics, Soft Matter, 2021, 17(28), 6707–6717,  10.1039/d1sm00600b.
35. J. H. Hildebrand, Solubility, J. Am. Chem. Soc., 1916, 38(8), 1452–1473,  DOI:10.1021/ja02265a002.
36. T. M. Winmostar and T. M. Winmostar Available from: https://winmostar.com/en/.
37. O. Buneman, Computer Simulation Using Particles (R. W. Hockney and J. W. Eastwood), SIAM Rev., 1983, 25(3), 425–426,  DOI:10.1137/1025102.
38. S. Mayo, B. Olafson and W. G. P. chemistry, DREIDING: a generic force field for molecular simulations, ACS Publ., 1990, 94, 91101,  DOI:10.1021/j100389a010.
39. S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys., 1984, 81(1), 511–519,  DOI:10.1063/1.447334.
40. G. J. Martyna, M. E. Tuckerman, D. J. Tobias and M. L. Klein, Explicit reversible integrators for extended systems dynamics, Mol. Phys., 1996, 87(5), 1117–1157,  DOI:10.1080/00268979600100761.
41. M. Tuckerman, B. J. Berne and G. J. Martyna, Reversible multiple time scale molecular dynamics, J. Chem. Phys., 1992, 97(3), 1990–2001,  DOI:10.1063/1.463137.
42. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, A smooth particle mesh Ewald method, J. Chem. Phys., 1995, 103(19), 8577–8593,  DOI:10.1063/1.470117.
43. K. J. Laidler, The development of the Arrhenius equation, J. Chem. Educ., 1984, 61(6), 494–498,  DOI:10.1021/ed061p494.
44. Y. Zhao, G. Kikugawa, Y. Kawagoe, K. Shirasu, N. Kishimoto and Y. Xi, et al., Uncovering the Mechanism of Size Effect on the Thermomechanical Properties of Highly Cross-Linked Epoxy Resins, J. Phys. Chem. B, 2022, 126(13), 2593–2607,  DOI:10.1021/acs.jpcb.1c10827.
45. N. Metropolis and S. Ulam, The Monte Carlo Method, J. Am. Stat. Assoc., 1949, 44(247), 335–341,  DOI:10.1080/01621459.1949.10483310.
46. S. Maeda, K. Ohno and K. Morokuma, Systematic exploration of the mechanism of chemical reactions: the global reaction route mapping (GRRM) strategy using the ADDF and AFIR methods, Phys. Chem. Chem. Phys., 2013, 15(11), 3683–3701,  10.1039/c3cp44063j.
47. S. W. Gabrielson, SciFinder, J. Med. Libr. Assoc., 2018, 106(4), 588,  DOI:10.5195/jmla.2018.515.
48. M. Livraghi, K. Höllring, C. R. Wick, D. M. Smith and A. S. Smith, An Exact Algorithm to Detect the Percolation Transition in Molecular Dynamics Simulations of Cross-Linking Polymer Networks, J. Chem. Theory Comput., 2021, 17(10), 6449–6457,  DOI:10.1021/acs.jctc.1c00423.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d3sm01743e

---