Block Chemistry for Accurate Modeling of Epoxy Resins

Accurate molecular modeling of the physical and chemical behavior of highly cross-linked epoxy resins at the atomistic scale is important for the design of new property-optimized materials. However, a systematic approach to parametrizing and characterizing these systems in molecular dynamics is missing. We therefore present a unified scheme to derive atomic charges for amine-based epoxy resins, in agreement with the AMBER force field, based on defining reactive fragments—blocks—building the network. The approach is applicable to all stages of curing from pure liquid to gelation to fully cured glass. We utilize this approach to study DGEBA/DDS epoxy systems, incorporating dynamic topology changes into atomistic molecular dynamics simulations of the curing reaction with 127,000 atoms. We study size effects in our simulations and predict the gel point utilizing a rigorous percolation theory to recover accurately the experimental data. Furthermore, we observe excellent agreement between the estimated and the experimentally determined glass transition temperatures as a function of curing rate. Finally, we demonstrate the quality of our model by the prediction of the elastic modulus based on uniaxial tensile tests. The presented scheme paves the way for a broadly consistent approach for modeling and characterizing all amine-based epoxy resins.

script.Atomic charges were equivalenced as dictated by geometrical symmetry.Intramolecular charge constraints were used to set the total charge of capping groups so that any polymerization state of DGEBA, resulting from appropriate linking of the EPX, PHE, LIN and EPO fragments, was charge neutral.Charge constraints for OPO were set so that any possible cured or partially cured configuration of the DGEBA -DDS network was also charge neutral.For example, the total charge of DGEBA with n=1 is EPX + PHE + LIN + PHE + EPO = -C1 -2C2 -2C3 -2C2 -C1 = C2 -2C2 + 2C2 -2C2 + C2 = 0. S7 Figure S2.An example of ring spearing observed in the simulations after a smaller interatomic distance than 2.5 Å was used during the system preparation.

The Curing Reaction
The possible underlying chemical reaction mechanisms of the curing reaction have been analyzed by Ehlers et al. 1 Applying density functional theory on small model complexes of epoxy and aliphatic amines, they found that the reaction favors a step-wise pathway, which includes the formation of a zwitterionic intermediate IM, over a concerted reaction involving a cyclic transition state (Figure 2a).Both mechanisms can benefit from hydrogen bonding interactions either via unreacted amine hydrogens or even better via the hydroxyl groups of the product species.In line with the analysis by Ehlers et al., we predict similar reaction kinetics for model systems involving primary and secondary aromatic amines as shown in Figure 2b.All structures were fully optimized at the wBP97-xD 2 /def2-TZVP 3 level of theory and characterized by subsequent frequency calculations at the same level of theory.The transition states were further confirmed by IRC calculations. 4,5Final energies were obtained at the DLPNO-CCSD(T) 6 /def2-QZVPP 7 level of theory with Orca 8 .All optimized structures are summarized in Figure S3 and the respective energetic barriers are stated in Table S1.Furthermore, we estimated the influence of solvent effects on the barrier, which should better reflect the situation in the liquid epoxy mixture (Table S2).Our curing procedure intersperses 100-ps equilibration periods among actual cross-linking cycles.

S9
We also performed curing on a large system with twice-as-long equilibration intervals of 200 ps (Figure S11).It can be readily seen that the differences between the 100 and the 200 ps profiles are negligible, indicating that 100 ps is a sufficient, as well as affordable, equilibration time (the disparity between the final density values associated with the two curves is of 0.8%).This shows that after 100 ps, the deviation from the average equilibrium density drops below 1%.The average density was evaluated after a 10 ns NPT run applied to one 75 % cured large system.The graph was built by sampling the density over a set of observation windows of increasing duration, from 100 to 1000 ps, and calculating for each of them the maximum deviation from the mean.This suggests that 100 ps are indeed sufficient for equilibration during curing.

Computational efficiency of fix bond/react
Table S3.Benchmarking results for an example crosslinking simulation using fix bond/react and a standard NPT equilibration run.Both simulations were run at 503 K for 1 ns and otherwise identical setups.All simulations were run on the Meggie cluster at the FAU, with each node equipped with two Intel Xeon E5-2630v4 "Broadwell" processors (10 cores per chip).interval for the center of the hyperbola fitted to the mean 95%-curing plot.Finally, the uncertainty on T g was computed by error propagation; the greatest contribution to the error on T g comes from the estimated uncertainties on the b parameters.The 95% data set shows an excellent hyperbolic trend, as confirmed by the overlap of the hyperbola asymptotes with the analytic tangents to the lowest and highest temperature points of the hyperbola.Conversely, the 15% data set clearly deviates from any hyperbolic trend, as testified by the disagreement between asymptotes and tangents.In this case, neither the intersection of the analytic tangents nor the center of the hyperbola are reliable estimates of T g .

Young's modulus
Table S5.All simulation parameters and results regarding post-curing annealing, equilibration and uniaxial tensile test.

Impact of partial atomic charges on curing and physical properties.
To investigate the impact of Coulomb interactions on the cross-linking process and on the deriving physical properties of the cured system, we carried out an additional simulation leaving out all Coulomb interactions (i.e.without charges), but still following the same procedure adopted in the presence of our block-chemistry charges.Specifically, we equilibrated and cured a system of large size (127000 atoms) and then evaluated its glass transition temperature, elastic modulus and percolation point.

Percolation point
The percolation point of the material cured without partial atomic charges was located at 58%, in alignment with our findings in the presence of charges.This agreement is to be expected, since the Flory-Stockmayer theory of polymer growth can describe gelation solely in terms of the functionalities of the reactants, predicting gelation exactly at 58% in the present case of a tetrafunctional hardener (DDS) and a bi-functional resin (DGEBA).

Glass transition temperature
We investigated samples at 35 % and 75 % curing as examples for curing degrees below and above the gel point.We determined the glass transition temperature using the same analysis procedure as for the block-chemistry approach and found that neglecting charges leads to drastic consequences.Specifically, at 35 % curing, T g drops from 78 ± 11°C in the system with charges, to 28 ± 9 °C in the one without charges.Similarly, at 75 % curing, the decrease is from 145 ± 11 °C to 94 ± 10 °C.This is far from the excellent experimental agreement achieved by our blockchemistry approach at all investigated degrees of cross-linking (Figure 4 of the manuscript).
This difference between the systems with and without charges is easy to assess by observing the density as a function of temperature at 35 % and 75 % curing (Fig. 6b in the manuscript).First, the slopes of the low-temperature tangents differ from each other, suggesting a scaling behavior different from the one presented in the manuscript for charged systems.Second, without charges the glass transition region is shifted to lower temperatures.Additionally, the density of the material polymerizing without electrostatic interactions is, at temperatures higher than the transition, systematically lower than the one obtained by using charges.

Elastic modulus
Finally, we performed in-silico tensile tests on our system cured without charges.By analyzing the linear response region, we determined the Young's modulus to be 2.015 ± 0.002 GPa.This constitutes more than a 30% drop in the ability of the material to withstand tensile stress, compared to the result obtained with charges (Young's modulus of 3.176 ± 0.002 GPa).

Figure S1 .
Figure S1.The charge constraints imposed during RESP charge derivation using the RED.III.5

Figure S3 .
Figure S3.Optimized model reactants, transition states and products for the curing reaction with primary and secondary amines following the scheme in Figure S3.The concerted mechanism is shown in a and b, the step-wise reaction in c and d and the step-wise/alcohol-catalyzed reaction in e and f.

S11 3 .
Figure S4.Degree of curing over 9.5 ns for one of our systems, when no annealing is applied.Samples collected every 100 fs.

Figure S5 .
Figure S5.The analysis of the effect of cut-offs on the distribution of reactive groups in the system as a function of curing.Cut-offs are presented with dashed red and blue vertical lines.Radial distribution functions at 0% (red lines), 40% (blue lines), and 80% (green lines) between reactive C and a) primary, b) secondary, and c) tertiary amines are shown.The insets represent the number of reactive C atoms as a function of the separation from the corresponding amine type.d) Numberof amine conversion during primary to secondary reaction (black line) and secondary to tertiary reaction (gray line) as a function of curing steps.No secondary amine is discovered up until the primary cutoff (2.8 A) and no tertiary amine are seen up until the secondary cutoff of 3.1 A (Figure1c).Furthermore, the conversion of primary amines clearly preceeds the conversion of secondary ones.

Figure S6 . 4 .
Figure S6.Evolution of the curing extent as a function of time for one of the five epoxy model systems.The slope of the curing curve decreases considerably after the first 4 ns of simulated cross-linking, due to a corresponding decrease in the diffusivity of the reactants.It rises again, however, after an unreactive annealing period at 800 K (marked in red).As the curing extent approaches 90%, longer annealing cycles are needed to obtain an appreciable increase in reaction rate.

Figure S7 .
Figure S7.Initial equilibration of the liquid precursor before curing, for all five large systems.The density equilibrated within 250 ps.

Figure S9 .
Figure S9.NPT equilibration of the same system as in Figure S4, 75% cured, at 503 K. Maximum % distance from the overall 10-ns average density value as a function of the time width of the window of observation.This shows that after 100 ps, the deviation from the average equilibrium

(
a) T g was obtained by inverting the linear transformation on the temperature axis as.  =  95% -Uncertainties on a and b were obtained by visually estimating which set of scaling parameters S17 yielded an acceptable overlay onto the 95% mean plot, keeping into account the standard deviation in density at each temperature data point.The uncertainty on represents a 95% confidence  95%

Figure S10 .
Figure S10.Density as a function of temperature at 95% (left) and at 15% (right) curing extents.
average over all simulations up to 0.02 strain.

Table S4 .
The parameters a, b and c used to rescale each density data set onto the mean 95% plot, the temperature corresponding to the center of the hyperbola fitted to the mean 95% plot,  95%and the resulting values of glass transition temperature T g , as a function of curing percentage.