A Molecular Dynamics Study of Monomer Melt Properties of Cyanate Ester Monomer Melt Properties

The objective of this work was to computationally predict the melting temperature and melt properties of thermosetting monomers used in aerospace applications. In this study, we applied an existing voids method by Solca. to examine four cyanate ester monomers with a wide range of melting temperatures. Voids were introduced into some simulations by removal of molecules from lattice positions to lower the free-energy barrier to melting to directly simulate the transition from a stable crystal to amorphous solid and capture the melting temperature. We validated model predictions by comparing melting temperature against previously reported literature values. Additionally, the torsion and orientational order parameters were used to examine the monomers’ freedom of motion to investigate structure–property relationships. Ultimately, the voids method provided reasonable estimates of melting temperature while the torsion and order parameter analysis provided insight into sources of the differing melt properties between the thermosetting monomers. As a whole, the results shed light on how freedom of molecular motions in the monomer melt state may affect melting temperature and can be utilized to inspire the development of thermosetting monomers with optimal monomer melt properties for demanding applications.


Introduction
Cyanate esters that polymerize into polycyanurates have excellent processability, low water uptake, low dielectric loss, and excellent thermal stability [1]. These thermosetting polymers occupy a niche area between high-glass-transition-temperature bismaleimides and tetrafunctional epoxy resins [2]. In this respect, they have been utilized for applications such as filament wound Carbon-Fiber-Reinforced Polymers (CFRPs) in heat shields used for atmospheric reentry and airframes [3]. The combination of high performance with good processability is a unique asset of polycyanurates compared to alternative high-temperature thermosetting resins. Research and development of these polymer matrices in the aerospace and defense industries can profit from the production of novel advanced thermosets and a better understanding of the underlying dynamics of these systems at an atomistic level. Accordingly, examining industrially relevant cyanate ester thermosetting monomers with Molecular Dynamics (MD) simulation is of particular interest in this study.
The melting temperature (T m ) of cyanate ester monomers is a key parameter in the fabrication of CFRPs with techniques such as filament winding and resin transfer molding [4]. Optimally, the cyanate ester exhibits a T m slightly above room temperature (315-340 K) yet is readily solidified for storage [5]. High T m (above 400 K) makes it difficult to maintain a monomer melt for adequate processing times without advancing the polymerization. Generally, measurement of monomers by differential scanning calorimetry is a common

Molecular Dynamics Simulations
MD simulations were performed to simulate the solid-to-liquid phase transition of cyanate ester monomers, calculate T m , and thereby appraise the quality of the molecular model relative to available experimental values ( Figure 2). The Schrödinger Materials Science Suite [21] and Desmond [22] were used to for model construction, analysis, and simulation. Simulations were performed in the canonical ensemble with a constant number of atoms, pressure, and temperature (NPT) with a Nosé-Hoover thermostat [23,24] and the Martyna−Tobias−Klein barostat [25] unless otherwise noted. The OPLS3e forcefield was used for all molecules [26][27][28].

Molecular Dynamics Simulations
MD simulations were performed to simulate the solid-to-liquid phase transition of cyanate ester monomers, calculate Tm, and thereby appraise the quality of the molecular model relative to available experimental values ( Figure 2). The Schrödinger Materials Science Suite [21] and Desmond [22] were used to for model construction, analysis, and simulation. Simulations were performed in the canonical ensemble with a constant number of atoms, pressure, and temperature (NPT) with a Nosé-Hoover thermostat [23,24] and the Martyna−Tobias−Klein barostat [25] unless otherwise noted. The OPLS3e forcefield was used for all molecules [26][27][28]. We examined four different cyanate ester monomers, with previously reported experimental melting temperatures: (1) 2′-(4-cyanatophenyl)propane (BADCy, Tm = 355.15 K) [29], (2) 4,4′-(ethane-1,1-diyl)bis(cyanatobenzene) (LECy, Tm = 302.15 K) [30], (3) tris(4cyantophenyl)methylsilane (SiCy-3, Tm = 390.62 K) [5], and (4) bis(4-cyanatophenyl)dimethylsilane, (SiMCy, Tm = 333.05 K) [31]. The crystal structures of BADCy [32], LECy [30], SiCy-3 [5], and SiMCy [5] have also been previously published. The structures for BADCy We examined four different cyanate ester monomers, with previously reported experimental melting temperatures: (1) 2 -(4-cyanatophenyl)propane (BADCy, T m = 355.15 K) [29], (2) 4,4 -(ethane-1,1-diyl)bis(cyanatobenzene) (LECy, T m = 302.15 K) [30], (3) tris(4-cyantophenyl) methylsilane (SiCy-3, T m = 390.62 K) [5], and (4) bis(4-cyanatophenyl) dimethylsilane, (SiMCy, T m = 333.05 K) [31]. The crystal structures of BADCy [32], LECy [30], SiCy-3 [5], and SiMCy [5] have also been previously published. The structures for BADCy (Refcode: TACHAG), SiCy-3 (Refcode: WAGWEJ), and SiMCy (Refcode: YERHAF) are available from the Cambridge Structural Database as .cif files [33]. For LECy, the crystallographic data were only available as unit cell vectors, so the .cif for BADCy was modified due to structural similarities. In this instance, a methyl group was deleted from each BADCy molecule in the unit cell, and subsequently a simulation cell was constructed from the edited structure. The supplied unit cells were replicated in all three directions to create supercells with 246 molecules (4 × 4 × 4) and 400 molecules (5 × 5 × 4). For 256-molecule ensembles (0% voids), LECy had the lowest number of atoms (8192), whereas SiCy-3 had the largest number of atoms (11,264). For larger systems with 400 molecules (0% voids), the number of atoms ranged from 12,800 to 17,600. The system sizes were chosen since they are large enough to ensure the boxes were greater than the minimum boundary condition (>20 nm), yet small enough to minimize the computational cost. For example, simulations containing 100,000 atoms take approximately 24 h to run on 8 cores. In view of this high computational cost of performing simulations with 100,000 atoms, the use of smaller system sizes was particularly attractive. The procedure to simulate the solid-liquid phase transition imitates the process of heating a crystalline solid from low temperature while tracking the change in volume of the system in tandem. To calculate T m from MD simulations, researchers have typically used the volume (or density) as a function of temperature [12,13,[34][35][36][37][38]. Within this thermodynamic theory of melting, at the melting temperature (constant temperature and pressure), the solid-liquid equilibrium is defined as: where s represents solid properties and l indicates liquid properties. For example, upon the first-order melting phase transition there is a discontinuous increase in volume. We used a void method, put forward by Solca whereby voids are "added" into the system by deleting a certain percentage of monomers from the source structure [20]. Molecules were selected at random and deleted to prepare supercells with defect concentration ranges from 1 to 10% to allow us to assess the effect of voids on the quality of T m calculations with regard to experimental values. The 10% voids upper limit was selected to prevent mechanical destabilization and collapse of simulation cells, which has been previously attributed to an excessive number of defects, shear instability, or excessive vibrational motion [39]. An energy minimization step was used to generate optimized geometrical configurations and to provide molecule systems with realistic densities. MD simulations were run stepwise from 150 K specified increments until 600 K was reached. The simulation cells were subjected to an equilibration period of 20 ns at each temperature step, during which properties were elicited [40].
It should be mentioned that the heating rate influences thermodynamic outcomes [41]; therefore, initially two different heating rates (25 K/20 ns and 50 K/20 ns) were tested to observe the reproducibility of calculated T m .

Orientational Order Parameter
The orientational order parameter was utilized to describe the solid-to-liquid phase transition. For a given configuration, this parameter is defined by: where cosθ i is derived from the scalar product of a vector for the descriptor (principal moment of inertia) and the vector for the director, and P 2 (cosθ i ) is the second-order Legendre polynomial, (3cos2θ i − 1)/2. The sum is calculated from all descriptor vectors defined for all molecules, and the order parameter is taken as the average value of the Legendre polynomial over all descriptors. P 2 is used so that the order parameter does not distinguish between descriptor vectors in opposite directions [21]. The order parameter equals 1 when molecules have perfect alignment along the director, which signifies a crystalline lattice in the solid state. This value decreases as the system becomes disordered and approaches 0 when molecules have isotropic orientations with respect to the director, which indicates liquid-like behavior [37,42]. The procedure of the calculation consisted of running an annealing simulation below the T m (T m − 50 K) and above the T m (T m + 50 K). The simulation cells were subjected to an NVT Brownian dynamics stage for 1200 ps, an NVT MD stage for 1.2 ns at 10 K, and an NPT MD stage for 1.2 ns at the desired temperature to reach an equilibrium value using the MD Multistage Workflow module [43].

Torsion
The torsion angle describing the orientation of the aromatic rings, relative to the bridging methyl group, was calculated using the Torsion Profile Analysis module. Trajectories were obtained in an analogous manner to the order parameter calculations. For all monomers, trajectory files for systems below and above the T m (T m − 50 K or T m + 50 K) were used as input for analysis over the simulation time. The SMARTS [44] pattern for the dihedral angle between the methyl carbon on bridging groups and aromatic carbons on the benzene ring was used to define the torsional angle.

Melting of Monomers with 0% Voids
The calculation of monomer melt properties relies on simulations that represent the existing physical systems and that are supported by empirical data. It is fundamental to our approach that simulating the systems should take less time compared to synthesizing and measuring properties experimentally. In an effort to reach a good compromise between simulation time and accuracy, we decided to apply a well-known void method [20] to our thermosetting monomer systems. Three parameters that may affect the calculated T m of simulated cyanate ester monomers were tested: system size, heating rate, and percentage voids. Simulations were performed as described above. The phase transition temperature of the simulated system was estimated by observing a discontinuity, to represent melting, or inflection, to represent a glass transition, from the temperature-volume curve. This is analogous to a T m calculation performed on a differential scanning calorimeter performed on themosetting monomers. The effect of system size on the calculated T m was examined ( Figure S1). We observed that increasing the number of molecules did not cause a significant change in the calculated T m except as noted in Supplementary Material, consistent with previous findings [20]. Based on this result, subsequent simulations were conducted with 256 molecules for efficiency. In addition, the optimal heating rate was tested on the 256-monomer systems. Simulations were run with stepwise temperature ramps of 25 K/20 ns and 50 K/20 ns ( Figure S2). It was determined that there was no significant difference in the phase transition temperature, so the 50 K/20 ns (and 256 mol systems) was a good compromise between the simulation size and heating rate. Figure 3 shows the evolution of volume-temperature curves (50 K/20 ns, 256 mol) for systems without voids. Due to the nature of the temperature increments, the T m was expected to be resolved within a 50 K temperature range, rather than a true discontinuity typical for a first-order phase transition. The intent was not to obtain an exact number, but rather a range sufficient to resolve a T m by 50 K, which would be useful for the purpose of ranking. Figure 3b shows the volume-temperature curve for LECy. The volume increases linearly up to 450 K, after which an abrupt increase in volume observed from 450 to 500 K. The T m is assigned to the onset of this discontinuity, which is evidence of the melting of a crystalline solid [35]. The volume-temperature graphs for BADCy and SiCy-3, shown in Figure 3a,c, do not exhibit a discontinuous T m . The discontinuous T m is likely not captured in our temperature range, so 600 K is taken as the calculated T m . For all systems with 0% voids, the predicted phase T m s are significantly greater than the experimental values. This finding is in good agreement with earlier observations of superheating for perfect crystalline lattices [13,45,46]. [13]. In the absence of voids, another important consideration is that melting is a fundamentally a kinetic phenomenon, and varying the heating rate is known to affect the observed Tm [47]. Heating rates achieved in MD simulations are intrinsically much faster than heating rates accessible for experimental systems. In this context, significant superheating was observed to occur since molecular motions are hindered below the Tm at short timescales accessible in MD simulations.

Dependence of Monomer Melting Temperature on Percentage of Voids
The following part of this study involved removal of 1 to 10% of molecules from simulation cells. This was to assess the effect of the void percent on the quality of the Tm calculation with regard to experimental values and to observe the monomer melt behavior on an atomistic scale. The 10% voids upper limit was established to prevent mechanical instability and collapse of the simulated crystalline solid into an amorphous solid [20]. The calculated Tms for all of the simulated monomer systems were overestimated, which was expected and has previously been attributed to the high-free-energy barrier to melting [13]. The barrier to melting for simulations without voids is largely on account that homogenous nucleation is the sole mechanism to initiate melting, which is determined by the probability of spontaneously forming liquid-like droplets in the solid phase [13]. In the absence of voids, another important consideration is that melting is a fundamentally a kinetic phenomenon, and varying the heating rate is known to affect the observed T m [47]. Heating rates achieved in MD simulations are intrinsically much faster than heating rates accessible for experimental systems. In this context, significant superheating was observed to occur since molecular motions are hindered below the T m at short timescales accessible in MD simulations.

Dependence of Monomer Melting Temperature on Percentage of Voids
The following part of this study involved removal of 1 to 10% of molecules from simulation cells. This was to assess the effect of the void percent on the quality of the T m calculation with regard to experimental values and to observe the monomer melt behavior on an atomistic scale. The 10% voids upper limit was established to prevent mechanical instability and collapse of the simulated crystalline solid into an amorphous solid [20]. The evolution of volume as a function of stepwise heating for monomers is shown in Figure S3. The variation in the calculated T m as a function of percentage of total number of molecules removed from simulation cells is shown in Figure 4a. The most striking and immediate observation for the plots is the decrease in the predicted T m with an increase in percent voids, complementing previous studies [11,20]. The lower T m reflects a decrease in the superheating effect in the direct heating to melting on account of a decrease in the freeenergy barrier to the formation of a solid-liquid interface to nucleate melting. In general, monomers LECy and SiMCy have the lowest T m s, BADCy an intermediate value, whereas SiCy-3 has the highest T m . sible inaccuracy related with superheating due to the rapid heating of simulated monomers. On the other hand, the current study does support the observations of Solca that the Tm decreases with an increase in the percentage of voids, then remains constant in a "plateau region" that corresponds to the thermodynamic Tm [20]. Alvares used the voids method to calculate Tm for CaO, and the plateau value was reported to span from 9 to 27% voids [40]. For Alavi, the plateau spanned from 5 to 10% voids, despite the similar approach [19]. These findings emphasize the difficulty in predicting, a priori, (and calculating computationally) the Tm of monomers. A challenge that becomes formidable when targeting additional constraints on reactivity and performance to satisfy composite applications.  While the voids were confirmed to be necessary to nucleate melting of the simulated thermosetting monomers, there were some issues with the application of the void method to these polyatomic monomers. In the traditional void method, the estimated T m has been shown to reach a constant value that is independent of the percentage. This was not observed in all thermosetting monomers simulated; therefore, an important consideration is the discrepancy from the traditional voids method, with regard to the absence of a flat region. It is evident that LECy and SiMCy are in best agreement with the conventional void method as they reach an appreciable plateau from 6 to 10% voids (Figure 4a). In the case of BPACy and SiCy-3, the systems do reach a narrow plateau where T m is independent of percent voids; however, it spans from 9 to 10% and is more narrow than expected. The authors recognize that an important consideration regarding the calculations is the lack of a theoretical foundation in the assumption that the thermodynamic Tm for a crystal coincides with a plateau region for a crystal lattice with a particular void size. Furthermore, a meaningful variable that may affect the T m values in the plateau region is the possible inaccuracy related with superheating due to the rapid heating of simulated monomers. On the other hand, the current study does support the observations of Solca that the T m decreases with an increase in the percentage of voids, then remains constant in a "plateau region" that corresponds to the thermodynamic T m [20]. Alvares used the voids method to calculate T m for CaO, and the plateau value was reported to span from 9 to 27% voids [40]. For Alavi, the plateau spanned from 5 to 10% voids, despite the similar approach [19]. These findings emphasize the difficulty in predicting, a priori, (and calculating computationally) the T m of monomers. A challenge that becomes formidable when targeting additional constraints on reactivity and performance to satisfy composite applications. Figure S4 compared the accuracy of calculated T m s against the experimental values, as a function of the percentage of voids. Most of the calculated values are well above the experimental T m partly due to the lack of heterogeneous nucleation and also due to the measurement technique (heating rate). Of note is the good agreement of calculated T m s with experimental values from 9 to 10% voids for most monomer systems. Since our simulations were run in 50 K steps, we expected the technique would be capable of ranking modeled systems, rather than precise numeric agreement. With this consideration, the 9% void systems were established for the comparison analysis between simulated monomer systems and are summarized in Table 1 and Figure 4b. Additionally, there is good agreement of the results on BPACy, LECy, and SiCy-3, and to a lesser degree SiMCy, with the experiment.
A point to note from Figure S4 is that the T m of SiMCy was most similar to the experimental value from 3 to 5% voids, but the accuracy diminishes further with increasing percentage of voids. This behavior has been observed elsewhere and is related to the mechanical instability and collapse of the simulated crystalline solid into an amorphous solid [20]. While this was unexpected when considering results from other monomers, when the volume versus temperature curves for SiMCy were plotted (Figure 5a), important differences between simulated behavior became evident. It can be seen that simulations of SiMCy with 9% and 10% voids exhibited an inflection, which suggests a glass transition to a rubbery regime. The appearance of the T g suggests the incorporation of high void percentages may cause the collapse of the solid lattice into an amorphous solid at the start of the MD simulation because the system is unable to maintain solid configuration [48]. It should be mentioned for the simulated systems with 256 molecules that the voids are located relatively close to each other. In the instance of high percent voids, these interfaces may overlap, transforming the solid-liquid phase into supercooling amorphous liquid. Considering the effect of size, it can be seen that for the system with 400 molecules, all the simulations exhibited a discontinuous solid-to-liquid melting phase transition (Figure 5b). In a comparison of the four monomers with 9% voids, it can be seen that the correct rank ordering is retained.
experimental Tm partly due to the lack of heterogeneous nucleation and also due to the measurement technique (heating rate). Of note is the good agreement of calculated Tms with experimental values from 9 to 10% voids for most monomer systems. Since our simulations were run in 50 K steps, we expected the technique would be capable of ranking modeled systems, rather than precise numeric agreement. With this consideration, the 9% void systems were established for the comparison analysis between simulated monomer systems and are summarized in Table 1 and Figure 4b. Additionally, there is good agreement of the results on BPACy, LECy, and SiCy-3, and to a lesser degree SiMCy, with the experiment.
A point to note from Figure S4 is that the Tm of SiMCy was most similar to the experimental value from 3 to 5% voids, but the accuracy diminishes further with increasing percentage of voids. This behavior has been observed elsewhere and is related to the mechanical instability and collapse of the simulated crystalline solid into an amorphous solid [20]. While this was unexpected when considering results from other monomers, when the volume versus temperature curves for SiMCy were plotted (Figure 5a), important differences between simulated behavior became evident. It can be seen that simulations of SiMCy with 9% and 10% voids exhibited an inflection, which suggests a glass transition to a rubbery regime. The appearance of the Tg suggests the incorporation of high void percentages may cause the collapse of the solid lattice into an amorphous solid at the start of the MD simulation because the system is unable to maintain solid configuration [48]. It should be mentioned for the simulated systems with 256 molecules that the voids are located relatively close to each other. In the instance of high percent voids, these interfaces may overlap, transforming the solid-liquid phase into supercooling amorphous liquid. Considering the effect of size, it can be seen that for the system with 400 molecules, all the simulations exhibited a discontinuous solid-to-liquid melting phase transition (Figure 5b). In a comparison of the four monomers with 9% voids, it can be seen that the correct rank ordering is retained.

Structural Characterization
Following the examination and validation of system construction with T m calculations, the structural parameters were also explored. Order parameter studies were performed to determine the possible correlation of local short-ranged ordering to T m and observe monomer melt behavior on the molecular level. The order parameter in MS suite was utilized and applied to monomer systems with 9% voids. The calculations involved the use of simulation cells slightly below (T m − 50 K) and above (T m + 50 K) the T m . The low-and high-temperature order parameters for each class of monomer were repeatable across multiple simulations and are summarized in Table 2, and snapshots of this order-to-disorder transition are shown in Figure S6. The simulated values were expected to be at a maximum for monomers with high T m s, and at a minimum for monomers with comparatively lower T m s. The orientational order parameters for the monomers below the T m vary from~0.6 to~0.3 and are ranked from highest to lowest in the following order: SiCy-3 > BADCy > LECy > SiMCy. SiCy-3 has the highest order parameter and the highest experimental T m of the monomers characterized in this study. Analysis shows a comparable ordering between BADCy and LECy systems, with BADCy being slightly higher. It is well-known that the BPA-derived cyanate ester monomer with two methyl groups on the bridgehead tend to have increased rigidity. Compared to BADCy, LECy has one less methyl group on the bridgehead group, resulting in an increase in flexibility and decrease in close contact ordering ( Figure 6). The correlation between order parameter and T m seen here lends support to the hypothesis that an increase in ordering between close contacts is the driving factor in the melting temperatures. However, it was surprising the same trend does not hold with SiMCy, as it has the lowest order parameter, yet the second lowest Tm. This result may be rationalized by the monomer chemical structure and that the increased flexibility from the silicon substituted bridgehead increases isotropy in the ordering, resulting in the lowest order parameter value. Work by Guenthner has shown that the rotational barriers experienced by the silicon containing monomer is low compared to BADCy (quaternary carbon) due to "extra" conformational freedom of SiMCy due to the quaternary silicon group [5]. For SiMCy, the two phenyl rings are not "locked in" and can occupy a variety of twist angles, which may rationalize the low-order parameter. It should also be mentioned that the low-order parameter value signifies the presence of amorphous domains in SiMCy simulation cells, which is in agreement with observation of a Tg for SiMCy systems with 9 to 10% voids. When considering the previous section where the calculated Tm value (9% voids) was compared to the experimental Tm value, the presence of amorphous domains, rather than the crystal structure, may result in the lack of congruous trends formerly observed with molecular simulations using voids method (underpredicted Tm).
For simulations run above Tm, all monomers exhibited observable drops in order parameter that approached zero over time, which was expected. As seen in Figure 7, BADCy However, it was surprising the same trend does not hold with SiMCy, as it has the lowest order parameter, yet the second lowest T m . This result may be rationalized by the monomer chemical structure and that the increased flexibility from the silicon substituted bridgehead increases isotropy in the ordering, resulting in the lowest order parameter value. Work by Guenthner has shown that the rotational barriers experienced by the silicon containing monomer is low compared to BADCy (quaternary carbon) due to "extra" conformational freedom of SiMCy due to the quaternary silicon group [5]. For SiMCy, the two phenyl rings are not "locked in" and can occupy a variety of twist angles, which may rationalize the low-order parameter. It should also be mentioned that the low-order parameter value signifies the presence of amorphous domains in SiMCy simulation cells, which is in agreement with observation of a T g for SiMCy systems with 9 to 10% voids. When considering the previous section where the calculated T m value (9% voids) was compared to the experimental T m value, the presence of amorphous domains, rather than the crystal structure, may result in the lack of congruous trends formerly observed with molecular simulations using voids method (underpredicted T m ).
For simulations run above T m , all monomers exhibited observable drops in order parameter that approached zero over time, which was expected. As seen in Figure 7, BADCy was observed to reach an equilibrium value of 0.09, indicating the crystalline lattice structure completely dematerializes [50]. Similar trends were observed for the other monomers ( Figure S5). These observations indicate an order-to-disorder transition as a result of the melting process, specifically an increase in mobility and decrease in associations [51,52]. To further investigate structure-property relationships and examine possible correlations between chemical structure (freedom of motion) and monomer Tm, Torsion Profile Analysis was performed on monomer systems. The module allows for the investigation of a property that is not readily accessible experimentally: phenyl ring torsion of monomer melts. It was hypothesized that the difference in Tm between monomers was influenced by the chemical structure and assumed flexibility of the bridgehead group [5]. Furthermore, phenylene flips would indicate the solid-to-liquid phase transition, with high-melting monomers having a narrow distribution of dihedral angles, indicating high rigidity. Conversely, it was expected that low-melting monomers would have broader distributions of dihedral angles, indicating increased monomer flexibility.
In order to examine this, equilibrated systems slightly below (Tm − 50 K) and above (Tm + 50 K) the Tm were examined. For most systems, 9% voids was used. For SiMCy, 7% voids was used, since the system maintained a crystalline structure prior to melting, as evidenced by a first-order solid-to-liquid phase transition. Figure 8 shows examples of changes in the dihedral distribution (phenyl group orientation compared to methyl bridging group) upon the solid-liquid phase transitions. In the case of BADCy, below the Tm, the mean values fluctuate from around −174 to −154°; −124° and −94 to −84°; −54° and −30 to 25°; 55° and 85 to 95°; and 125° and 155 to 175°. The peak positions of the BADCy dihedral angle shift as the temperature is increased from 300 to 450 K. The pronounced peaks in the distributions merge and range from −175 to −125°, −95 to 55°, and 95 to 175°. Not only do the peak widths in the distribution increase, but the peak heights also decrease upon melting, signifying a decrease in the barrier to motion as the relative populations of conformers broaden [53]. The torsion of phenyl rings with respect to the two methyl groups at the bridge is able to assume a broad range of conformations, and the ring acquires more freedom to rotate. A similar trend is seen with LECy (Figure 8b). Below the Tm, the torsion has maxima at −153 to −143°, −113 to −64°, −34 to −24°, 15 to 35°, 65 to 75°, 94°, 114°, and 144 to 154°. Above the Tm, the probability greatly changes, and the maxima occur from −142 to −123°, −63 to −24°, 26 to 56°, and 105 to 145°. Compared to BADCy, it can be seen that LECy has more broad peaks below the Tm. This broadening indicates a Figure 7. The orientational order parameter for BADCy 50 K below the calculated melting temperature (black line) and 50 K above the calculated melting temperature (blue line). Included is a snapshot of monomers in the ordered state below the calculated melting temperature and disordered state above the calculated melting temperature.
To further investigate structure-property relationships and examine possible correlations between chemical structure (freedom of motion) and monomer T m , Torsion Profile Analysis was performed on monomer systems. The module allows for the investigation of a property that is not readily accessible experimentally: phenyl ring torsion of monomer melts. It was hypothesized that the difference in T m between monomers was influenced by the chemical structure and assumed flexibility of the bridgehead group [5]. Furthermore, phenylene flips would indicate the solid-to-liquid phase transition, with high-melting monomers having a narrow distribution of dihedral angles, indicating high rigidity. Conversely, it was expected that low-melting monomers would have broader distributions of dihedral angles, indicating increased monomer flexibility.
In order to examine this, equilibrated systems slightly below (T m − 50 K) and above (T m + 50 K) the T m were examined. For most systems, 9% voids was used. For SiMCy, 7% voids was used, since the system maintained a crystalline structure prior to melting, as evidenced by a first-order solid-to-liquid phase transition. Figure 8 shows examples of changes in the dihedral distribution (phenyl group orientation compared to methyl bridging group) upon the solid-liquid phase transitions. In the case of BADCy, below the T m , the mean values fluctuate from around −174 to −154 • ; −124 • and −94 to −84 • ; −54 • and −30 to 25 • ; 55 • and 85 to 95 • ; and 125 • and 155 to 175 • . The peak positions of the BADCy dihedral angle shift as the temperature is increased from 300 to 450 K. The pronounced peaks in the distributions merge and range from −175 to −125 • , −95 to 55 • , and 95 to 175 • . Not only do the peak widths in the distribution increase, but the peak heights also decrease upon melting, signifying a decrease in the barrier to motion as the relative populations of conformers broaden [53]. The torsion of phenyl rings with respect to the two methyl groups at the bridge is able to assume a broad range of conformations, and the ring acquires more freedom to rotate. A similar trend is seen with LECy (Figure 8b).
Below the T m , the torsion has maxima at −153 to −143 • , −113 to −64 • , −34 to −24 • , 15 to 35 • , 65 to 75 • , 94 • , 114 • , and 144 to 154 • . Above the T m , the probability greatly changes, and the maxima occur from −142 to −123 • , −63 to −24 • , 26 to 56 • , and 105 to 145 • . Compared to BADCy, it can be seen that LECy has more broad peaks below the T m . This broadening indicates a lower barrier to motion for LECy compared to BADCy, which would be expected on the basis of the decrease in rotational steric hinderance upon replacing one methyl bridging group with a hydrogen atom. Ultimately, the lower barrier to motion results in a concomitant decrease in T m, and the expectedly lower order parameter seen for LECy in the previous section. However, when considering the distribution of dihedral angles for all of the monomers in the data set, important differences between monomer freedom of motion become apparent. Figure 8 shows the broadest peaks in the melt (above the T m ) are in the order of: BADCy > SiMCy > SiCy-3 > LECy. However, the T m values occur in the following order: SiCy-3 > BADCy > SiMCy > LECy. While the computed results appear contradictory to previously reported experimental value, this suggests that the barrier to motion (internal entropy) is not the sole contributor to the T m , and other factors should be considered, such as interactions with neighboring molecules (intermolecular interactions). These results align with those of Ghiassi, who used X-ray crystallography and thermal and computational simulations to study the structure-property relationships that determine the melting characteristics for a series of cyanate ester monomers. The authors similarly concluded the expected correlation between restricted motion and entropy of melting (melting characteristics) was not observed for the examined monomers [7]. basis of the decrease in rotational steric hinderance upon replacing one methyl bridging group with a hydrogen atom. Ultimately, the lower barrier to motion results in a concomitant decrease in Tm, and the expectedly lower order parameter seen for LECy in the previous section. However, when considering the distribution of dihedral angles for all of the monomers in the data set, important differences between monomer freedom of motion become apparent. Figure 8 shows the broadest peaks in the melt (above the Tm) are in the order of: BADCy > SiMCy > SiCy-3 > LECy. However, the Tm values occur in the following order: SiCy-3 > BADCy > SiMCy > LECy. While the computed results appear contradictory to previously reported experimental value, this suggests that the barrier to motion (internal entropy) is not the sole contributor to the Tm, and other factors should be considered, such as interactions with neighboring molecules (intermolecular interactions). These results align with those of Ghiassi, who used X-ray crystallography and thermal and computational simulations to study the structure-property relationships that determine the melting characteristics for a series of cyanate ester monomers. The authors similarly concluded the expected correlation between restricted motion and entropy of melting (melting characteristics) was not observed for the examined monomers [7]. The dihedrals between aromatic rings and a bridging methyl group are of particular interest because they define the phenyl ring positions in the molecule. In order to explore this further, the dihedral angle was plotted as a function of time during equilibration for BADCy ( Figure 9). The sequence of images shows the time-dependent distributions of dihedral angles that demonstrate the phase transition of BADCy monomers upon heating above the Tm. In the crystal state the torsion angles were relatively stagnant and distributed without deviation along a gradient, as expected. At 450 K, a fluctuation in the dihedral angle of phenyl group orientation compared to methyl bridging group suggests ring flips occur, likely due to the increase in thermal energy when the melting temperature is approached. The same trend can be seen for all systems ( Figure S7), which show fluctuating dihedral angles upon melting, in accordance with the increase in thermal energy and freedom of motion. The dihedrals between aromatic rings and a bridging methyl group are of particular interest because they define the phenyl ring positions in the molecule. In order to explore this further, the dihedral angle was plotted as a function of time during equilibration for BADCy ( Figure 9). The sequence of images shows the time-dependent distributions of dihedral angles that demonstrate the phase transition of BADCy monomers upon heating above the T m . In the crystal state the torsion angles were relatively stagnant and distributed without deviation along a gradient, as expected. At 450 K, a fluctuation in the dihedral angle of phenyl group orientation compared to methyl bridging group suggests ring flips occur, likely due to the increase in thermal energy when the melting temperature is approached. The same trend can be seen for all systems ( Figure S7), which show fluctuating dihedral angles upon melting, in accordance with the increase in thermal energy and freedom of motion. Polymers 2022, 14, x FOR PEER REVIEW 13 of 16 Figure 9. Display of angles for each individual torsion as a bar plot for the simulation of a time range with colors denoting the size of the angle for BADCy.

Conclusions
This research was aimed at studying the effect of chemical structure on the melting behavior of cyanate ester thermosetting monomers and investigating the performance of the voids method introduced by Solca in capturing the melting temperature [20]. MD simulations were carried out on the crystalline monomers with void percentages up to 10%. A computational framework and a compared analysis of simulated Tm is presented and validated through comparison to previously reported empirical data. Additionally, the torsion and orientational order parameter were used to examine the monomers' freedom of motion in order to determine structure-property relationships. It was found that for crystal lattices with up to 400 molecules, 9% voids provided the best estimate of Tm for most systems. The orientational order parameter and torsion simulations gave insight on the local short-ranged ordering of monomers and monomer flexibility. The monomers with lower melting temperatures tended to have lower orientational order parameters than monomers with higher values, as expected [54], based on molecular symmetry and packing. Furthermore, these results, in tandem with torsion simulations, indicate that the internal barrier to motion is not the sole contributor to Tm, and intermolecular interactions must also be considered. As a whole, the results shed light on how the chemical structure of cyanate ester monomers may affect Tm, which was both computationally and experimentally determined, and can be utilized to inspire the development of thermosetting monomers with optimal monomer melt properties for demanding applications.  Figure S3: The evolution of volume as a function of temperature starting at T = 150 K for (a) BADCy, (b) LECy, (c) SiCy-3, and (d) SiMCy with 1 to 10% voids, Figure S4: Percent error of the calculated melting temperature, compared to the experimental values, for simulations with 0 to 10% voids (in % of the total monomers removed from simulation cell) for: (a) BADCy, (b) LECy, (c) SiCy-3, and (d) SiMCy, Figure S5: The orientational order parameter for monomers: (a) LECy, (b) SiCy-3, and (c) SiMCy 50 K below the calculated melting temperature (black line) and 50 K above the calculated melting temperature (blue line), Figure S6: Snapshots of the monomers at different stages of the melting process with 9% voids for the following monomers:

Conclusions
This research was aimed at studying the effect of chemical structure on the melting behavior of cyanate ester thermosetting monomers and investigating the performance of the voids method introduced by Solca in capturing the melting temperature [20]. MD simulations were carried out on the crystalline monomers with void percentages up to 10%. A computational framework and a compared analysis of simulated T m is presented and validated through comparison to previously reported empirical data. Additionally, the torsion and orientational order parameter were used to examine the monomers' freedom of motion in order to determine structure-property relationships. It was found that for crystal lattices with up to 400 molecules, 9% voids provided the best estimate of T m for most systems. The orientational order parameter and torsion simulations gave insight on the local short-ranged ordering of monomers and monomer flexibility. The monomers with lower melting temperatures tended to have lower orientational order parameters than monomers with higher values, as expected [54], based on molecular symmetry and packing. Furthermore, these results, in tandem with torsion simulations, indicate that the internal barrier to motion is not the sole contributor to T m , and intermolecular interactions must also be considered. As a whole, the results shed light on how the chemical structure of cyanate ester monomers may affect T m , which was both computationally and experimentally determined, and can be utilized to inspire the development of thermosetting monomers with optimal monomer melt properties for demanding applications.  Figure S3: The evolution of volume as a function of temperature starting at T = 150 K for (a) BADCy, (b) LECy, (c) SiCy-3, and (d) SiMCy with 1 to 10% voids, Figure S4: Percent error of the calculated melting temperature, compared to the experimental values, for simulations with 0 to 10% voids (in % of the total monomers removed from simulation cell) for: (a) BADCy, (b) LECy, (c) SiCy-3, and (d) SiMCy, Figure S5: The orientational order parameter for monomers: (a) LECy, (b) SiCy-3, and (c) SiMCy 50 K below the calculated melting temperature (black line) and 50 K above the calculated melting temperature (blue line), Figure S6: Snapshots of the monomers at different stages of the melting process with 9% voids for the following monomers: