Cut-off Scale and Complex Formation in Density Functional Theory Computations of Epoxy-Amine Reactivity

Most of the properties of epoxy resins are tied to their degree of cross-linking, making understanding the reactivity of different epoxy systems a crucial aspect of their utilization. Here, epoxy-amine reactivity is studied with density functional theory (DFT) at various cut-off levels to explore the suitability of the method for estimating the reactivity of specific epoxy systems. Although it is common to use minimal structures in DFT to reduce computational cost, the results of this study highlight the important role of hydrogen bonding and other noncovalent interactions in the reactivity. This is a promising result for differentiating the most probable reactive paths for different resin systems. The significance of amine groups as a potential source of catalyzing H-bonds was also explored and, while not quite as effective as a catalyst as a hydroxyl group, a clear catalyzing effect was observed in the transition state energies. Unfortunately, the added complexity of a more representative reactive system also results in increased computational cost, highlighting the need for proper selection of structural cutoffs.


INTRODUCTION
For many epoxy applications, in-depth knowledge of the resin system chemistry is essential. The basic reactions of epoxy resins can be considered well understood, but the kinetics of curing, 1−3 ageing behavior, 4,5 polymer network formation, 3,6,7 and options for tailoring the structure and properties 8−10 remain subjects of active study in the field. The studies of the reactivity at the scale of individual molecules, bonds, and small differences between spatial isomers are only possible through ab initio or the computationally less arduous density functional theory (DFT) computations. 11−13 The performance and approximations offered by modern DFT implementations enable approaching the scale of full resin molecules, expanding the potential user base of the methods also toward materials science and chemical engineering.
Research on epoxy-amine reactivity commonly relates to different modifications of the base reaction, for example, the catalysis provided by hydrogen bonding. 14 The possible effects of hydrogen bonds (H-bonds) are twofold: H-bonding can lower the reactive barrier (catalysis) or stabilize the intermediate structure, causing a retardation effect. 13 Raman and Palmese 13 studied the latter by adding tetrahydrofuran into an epoxyamine mixture. A similar retardation could also result from Hbond acceptors in the structure of the reagents, such as oxygen bridges in the epoxy and/or amine backbone. The catalyzing Hbonds are mostly considered to form with the hydroxyl groups formed during the addition reaction between the primary amine and the epoxide ring. 11−15 Most studies ignore the possible role of another H-bond donorthe amine groups themselves likely due to the lower reactive barriers reported for hydroxyl group-catalyzed reactive paths. 11 However, the abundance of available amine end-groups in epoxy-amine systems, especially in the early stages of the curing reactions, makes the amine groups a promising source of H-bonds.
In addition to H-bonding, epoxy-amine reactions are affected by inductive effects from the various chemical groups connected to the reactive center. 14 One example is the α-effect. 16 The αeffect is strong in methylated amines 17 and can, therefore, contribute to the results of computations of simplified structures. At least, if correctly predicted by the computation. The inductive effects are much weaker for larger structures 18 but, nevertheless, need to be considered in the total effect substituents have on the energy of the reaction and when estimating how changing the scale of the model system affects the computational results.
In this study, epoxy-amine reactions are studied computationally using modern DFT methods. The energy levels of the structures along the reactive path are computed at varying structural cutoffs. The aim is to explore the accuracy versus computational cost relation for a model epoxy-amine system and to see the effects various neighboring chemical groups can have on the reactive path. This should allow the qualitative prediction of the most probable reactive path for a given cutoff and offer insights into the complex role of noncovalent bonding in stabilizing the structures along the reactive path. The different cut-off structures are manually built to represent commonly encountered epoxy-amine reactive systems with increasing complexity. For example, all epoxy structures are simplified cutoffs of diglycidylether of Bisphenol-A (DGEBA).

COMPUTATIONAL METHODS
The numerical computation was performed using the Schrodinger Materials Science Suite (Schrodinger Inc., New York, USA) software package (version 2020-4)most prominently the Jaguar (version 11.0) DFT program. 19 The longrange corrected hybrid nonlocal B3LYP-D3 theory level was used for all DFT computations. The functional has been shown to give a good combination of accuracy and computational efficiency with many types of chemical structures and noncovalent bonding cases. 20,21 The M06-2X and ωB97X-D functionals along with B3LYP-D3 with Becke−Johnson damping 22 were tested in single point energy (SPE) computations for comparison. This part of the discussion is presented in the Supporting Information. The basis set was varied depending on the computational task. The 6-31G** basis set was used for the transition state searches and other initial structural optimizations.
The role of the structural cutoffs was explored using reaction workflows recently introduced to the Schrodinger Materials Science Suite. These simulations involve freezing the primary reactive groups of the structures and swapping fragments from the rest of the structure based on pregenerated examples shown in Figure 1. Ethylamine was selected as A1, instead of the simpler methylamine, to reduce the otherwise expected significant contribution of the α-effect for the reference structure. Ethylamine was also used as the catalyzing amine for all relevant simulations. Similarly, isopropyl alcohol was used as the source of the catalyzing hydroxyl group. In the figure (Figure 1), carbons are colored green, nitrogens blue, oxygens red, and hydrogens white. The same coloring will be used throughout this study.
In order to find the most probable reactive paths, a conformational searchusing the MacroModel program with the OPLS3e force-field 23 for Monte Carlo Multiple Minimum searcheswas added to each transition state search and reaction workflow. For the transition state searches only the lowest energy conformers are available as outputs, whereas in the reaction workflows ten lowest energy conformers are requested from separate conformational search runs the same starting structure. The final structural optimizations in DFT were performed using the LACVP** basis set, whichas the workflow defaultresulted in significantly less convergence issues in the workflow than the 6-31G** basis set.
The final energies of each structureincluding the outputs of the transition state searcheswere calculated as SPE computations using the cc-pVTZ(-F) triple zeta basis set for improved accuracy. Other relevant parameters of the Jaguar code include grid density for the DFT computations (fine for SPE, medium for all other computations) and the accuracy of the pseudospectral computations (ultrafine grids with tight cutoffs for SPE, mixed grid with loose cutoffs for all other computations).
All results are reported as averages and standard deviations of grouped sets of these computational results. The conformers are grouped based on the used fragments ( Figure 1) and subgrouped based on conformational similarities, for example, similar H-bonding states. The computation times are reported from SPE computations. Because the actual duration of the computation depends on the computational performance of the hardware and, for example, the use of parallelization to speed up the computation, the values are reported relative to the fastest computation (A1,E1 uncatalyzed primary amine reaction) and are based on the total CPU time.

RESULTS AND DISCUSSION
The computations begin with a transition state search for the simple A1,E1 system. This was needed both to find the reference energy levels for the simplified system and for use as an input for the reaction workflow computations. Figure 2 presents the results of the transition state searches used to find the primary amine addition reactive paths, including the A1,E1 version of the amine-catalyzed reactive path. These act as starting points for the rest of the computations.
The results for uncatalyzed reactions are presented in Table 1 as transition state energies relative to the lowest reactant structure energy. This relative energy is used as a measure for the activation energy (E a ) or reactive barrier height.
The β-methylated hardeners (A2−A4) are known to result in slightly slower reactive systems compared to their linear counterparts. This is caused by the steric hindrance from the methyl group next to the amine. 14 The effect is overlooked by the simplified reactive system (A1,E1), but clearly evident as higher activation energies for each larger amine structure (A2− A4,E1). This already highlights the crucial role of the selection of the cutoffs when studying specific reactive systems.
For the more complex structures, certain specific noncovalent interactions emerge anddespite the vast number of possible conformations for the rest of the moleculethe structures involving these interactions share structural and energetic similarities. To highlight these specific interactions, they are presented in Schemes 1−3.
Contributions to the total energy can also be noted from conformations where the presence of H-bonding is uncertain but the proximity creates a more favorable energetic state. This is most common, for example, between an amine hydrogen and the nearby backbone oxygen (see, e.g., Figure 3a). Although similar to Scheme 1b, for larger structures, these often present as conformers with multiple oxygen groups coordinating around the amine hydrogen, as presented in 3e, Schemes 2a and 3. This type of complex is presented in Scheme 4. The hydroxyl group in the structure E2 manifested mostly conformations where the hydroxyl group interacts with the epoxide oxygen (See Figure 3c). Although accurate for the structure in question, these conformations are considered bad for the purposes of exploring the reactivity of a DGEBA-based epoxy system, in which the oxygen is part of the epoxy backbone (see e.g. Figure 3d,e). The A3 hydroxyl group also manifested some conformations that would be very unfavorable in a larger system (Figure 3b) at least with the assumption that this oxygen is similarly part of the amine backbone. Due to these considerations the choice was made to exclude these systems from the computations of the catalyzed systems.
Interactions are also observed between the aromatic ring and the polar functional groups. First examples of this type of behavior are observable from the A1,E3 secondary amine reaction, where the average height of barrier changes approximately 1−2 kcal/mol depending on the type of interaction. The highest reactive barriers (32.00 ± 0.48 kcal/ mol) are observed for the conformers with no interactions between the hydroxyl group from the primary amine reaction and the aromatic ring in the epoxy structure. For a PhH → OH type interaction (Scheme 3b), the average barrier height decreases to approximately 31.13 ± 0.21 kcal/mol. Aromatic H-bonds are more commonly reported between an H-bond donor and the aromatic ring. 24 It is worth noting that the aromatic hydrogen in ortho-position to an oxygen is expected to have a slightly more acidic nature compared to a benzene hydrogen, which makes such interactions plausible. The lowest energy conformers (average of 29.01 ± 0.77 kcal/mol) present a OH → Ph type interaction (Scheme 3a) resembling the aromatic H-bond discussed by Brinkley and Gupta. 24 The examples of such structures are presented in Figure 4. Similar interactions were observed in the uncatalyzed A4,E4 secondary amine system.
When functional groups (amine or hydroxyl) are added to catalyze the reaction, the H-bonding of these added functional groups clearly changes the energetics of the reaction, which is in line with existing studies on the topic. 11,15 The catalysis is based on the extra functional group interacting with the opening epoxide ring, as presented in Scheme 5.
Mobility around the reactive center is likely overestimated in such a small model system, that is, it is improbable that every reactive site would have sufficient free volume surrounding it to allow this type of catalysis. Nevertheless, the model still helps highlight the differences between different options for H-bond formation, namely the hydroxyl and amine groups, and helps estimate the role of such catalysis in total reactivity. The transition state energies for the catalyzed reaction paths are presented in Table 2.
The H-bond between the amine group and a backbone oxygen (Scheme 1) appears a favorable conformation for systems large enough to orient suitably. Figure 5 shows examples of the effect this interaction has on the energy in the A1,E4 system. The effect on barrier height is unexpectedly minor. A possible explanation is that the orientation forces some bond angles to deviate from their equilibrium values increasing the energy. Figure 6 shows the lowest and highest energy conformers for the amine-catalyzed A2,E1 secondary amine transition states. A complex H-bond network formed between the hydroxyl groupcreated in the primary amine reactionand the amine-catalyzed reactive core (Scheme 5b) significantly lowers the reactive barrier for this path. It is also worth noting that the required orientation of the hydroxyl group also brings it close to the remaining amine hydrogen. The combined effect is quite significant (29.21 ± 0.90 kcal/mol VS 35.19 ± 0.64 kcal/mol).
In general, for the larger systems, the conformations and their interpretation become challenging due to the number of available interactions including possible complex cases such as the one presented in Scheme 4. For example, in the catalyzed A4,E4 systems, multiple noncovalent interactions are noted for most conformers (see Figure 7). The lowest energy conformations show oxygens coordinating around the amine hydrogen as in Scheme 4. In the higher energy conformers, however, either the amine backbone oxygen or the hydroxyl group from the primary amine is oriented away from the amine hydrogen. The difference in energies between these conformations appears to be in the range of 2−4 kcal/mol but in all cases other interactions contribute. A similar coordination around the amine hydrogen can also be observed in the uncatalyzed system presented in Figure 3h, and a similar difference in energy is observed when one of the oxygens is oriented less favorably (see Figure 3g). Also, again additional interactions contribute, such as the OH → O H-bond presented in Figure 3h.
The aromatic interactions contribute significantly also in the catalyzed models. With the reactive center largely the same for a majority of the structures, the trend of increasing deviations for larger systems can be often attributed to how, for example, the hydroxyl group formed in the primary amine reaction is oriented in relation to the aromatic ring in the secondary amine reactions. The strength of the interactionfor cases resembling Scheme 3ais in these results similar to what is observed for H-bonding, giving validity to the term of aromatic H-bonding.
Predicting the actual importance of the catalyzed paths using the DFT computation is difficult as the systems are effectively considered infinitely diluted, whereas molecular motion in the real epoxy-amine reactive systems is far more restricted. The observed catalyzing effect partly explains the accelerated reactivity when amine hardener is added in excess of the stoichiometric ratio. It is also very likely that whatever importance the catalysis has, it is mitigated as the curing process progresses, and the required molecular motions are restricted until itafter vitrificationbecome practically nonexistent. Table 3 presents an overview of the computational times the easiest measure for the computational costfor the different systems. Some interesting observations can be made from these results. The most obvious one is as follows: increasing the complexity of the system increases the computational cost but also variation of the computational time. In the SPE computation, the biggest source of variation between similar structures appears to be the dispersion correction steps in the computation. The slightly longer computational times for hydroxyl-catalyzed compared to amine-catalyzed systems are likely a result of the slightly bigger system (isopropanol vs ethylamine).   The reactant conformations of the reaction workflows provide an interesting point of discussion. Especially for the catalyzed reactive paths the lowest energy conformers show significant interactions between the added catalyzing functional group (OH or NH) and either one or both of the amine and epoxy functional sites. These minimal energy conformers tend to orient unfavorably considering the reactive path of the epoxyamine addition reaction. Whether the needed reorientation manifests as an initial smaller reactive barrier, could be an interesting point of further study. However, the effect is likely small compared to the advantage offered by the lowering of the reactive barrier in the rate determining step studied here.
No conclusive retardation effect was observed based on the computation results. However, some hints at possible metastable intermediate structures can be drawn from the results of the conformational searches of the intermediate structures. For example, the uncatalyzed primary amine reaction workflows resulted in intermediate structures close to the final productindicating no local minima near the transition state. Whereas for the secondary amine reaction, the conformational search converged to local minima structurally and energetically near the transition state for all systems larger than the A1,E1 system. This could indicate the secondary amine reaction for our model system is in fact slower, even though the activation energies are largely similar. Here, as the core reactive path was kept largely constant in the set-up of the computations, possible key considerations, such as a backbone oxygen interacting and "deactivating" the catalyzing functional group, were not explored.

CONCLUSIONS
In this study, the epoxy-amine reaction was studied using the DFT computations at different cut-off scales to study the role various noncovalent interactions of the neighboring functional groups. A secondary goal was to test the potential of such computations in predicting the reactivity of a specific epoxyamine system, at least qualitatively, with promising results. The selected cutoffs significantly affect the final energy levels. Most   Scheme 5. Catalysis of the Epoxy-Amine Reaction by an Extra functional Group (OH or NH 2 = R 4 ) a a R 1 = epoxy backbone, R 2 = amine backbone, R 3 = either hydrogen or the epoxy from primary amine reaction. In the models presented here R 5 is connected to one of R 1-3 . consistent results were achieved with the largest and the smallest cut-off levels. The energy levels for the smallest and the largest structures are also very similar, which should be coincidental. Various important interactions were included in the computa-tions as a result of the increased complexity of the systems. Based on our computations, the order of importance for different types of interactions to epoxy-amine reactivity is as follows: The conformers are subgrouped based on notable common interactions in Schemes 1−5. "No catalysis" indicates conformations with no predicted interaction between the catalyzing functional group and the opening epoxide ring.  intermolecular H-bonding, intramolecular H-bonding, aromatic interactions, and inductive effects, from most important to least important. Computations with the smallest cutoffs (largest structures) unveil the possible effects of the chemical groups near the reactive site that can stabilize the transition state or lower the energy level and accelerate the reaction. However, these overlapping phenomena proved to be quite challenging to analyze. Expanding the size of a reactive system in DFT computations must be approached with care. A simulation more representative of your real system can offer insights into the behavior of the system. The drawbacks are losing generality and, of course, a significant increase in computational cost..   The reference time is approximately 2332 cpu seconds.