Triacylglyceride melting point determination using coarse-grained molecular dynamics

This study is carried out using the COGITO force field to determine whether the thermodynamic melting point of pure triacylglyceride crystals can be predicted using molecular dynamics simulations. The triacylglycerides used in this study are both saturated and unsaturated, as well as symmetrical and asymmetrical, to test the robust-ness of both the force field and the direct heating methodology described in this paper. Given the nonequilibrium nature of a melting system, a larger number of simulations are required to ensure that the results are sufficiently converged, that is, with little fluctuation and a small confidence interval. The study also highlights the impor-tance of the presence of defects, in this case as voids, to lower the melting nucleation energy barrier of the crystals and avoid superheating of the systems being tested. The size of these defects is much larger than what would be found in a physical crystal, however, the simple and robust procedure that was developed allows the accurate prediction of melting points of the different triacylglycerides.


| INTRODUCTION
The solid-to-liquid phase transition temperature is an important macroscopic property of any material, and the simulation of this phase transition has been the research focus of a number of groups, 1-8 who have attempted to simulate the melting of various metals or molecules. Determining the melting point computationally using molecular dynamics (MD), however, is not trivial. The force field (FF) being used, and how this has been parameterized, that is, whether melting point has been part of the parameterization efforts, or not, has a significant impact on being able to predict melting points. Having a FF which is able to predict liquid and solid macroscopic properties does not imply the FF's ability to predict phase transition, 3 and hence any such simulations need to be validated thoroughly.
A common problem encountered when trying to determine the melting point of a crystal using MD is overheating (or superheating) of the system. 2,4,7,9 This is due to the system being a perfect crystal, and thus having no surface, defect or interface, or alternatively, superheating may be caused by a very high heating rate compared to empirical measurements or, the influence of the FF potential. 7 The issue of having a defect, or void, within a crystal when carrying out MD simulations has been studied extensively 1,2,4,8,9 and shown to be critical to obtaining a good determination of the melting point using MD. Having such void defects in a crystal allow for the fact that physical crystals are never perfect, and that boundaries and lattice defects reduce the possibility of crystal superheating. 2 Superheating can be mitigated by a slower heating rate, however, completing a simulation at heating rates approximating those obtainable in a lab, over a similar temperature range, is not feasible as the computation time would be too long.
Different approaches have been proposed to determine the simulated melting point. The first is a direct approach method, that is, one where the system is heated until phase transition is observed, mimicking an experimental procedure. The second approach is to use free energy methods, with the latter being deemed to be more accurate than the former, however, it is also more complicated to carry out, 10 and can be very difficult if the solid-liquid coexistence conditions are not known in advance. 2  with the COGITO FF. 11 All equilibrations were done using an NPT ensemble, using a time-step of 25 fs and a pressure of 1.01325 bar using a v-rescale thermostat and a Berendsen barostat. Anisotropic pressure coupling, with a compressibility of 1 Â 10 À5 bar À1 in the x, y and z directions was used. Temperature coupling was set at 1 ps, while pressure coupling was set at 10 ps. The cut-off scheme was set to Verlet, with the Coulomb and vdW cut-off distances set to 1.1 nm. The vdW-modifier was set to Potential-shift and the electrostatics (coulomb type) was set to Particle-Mesh Ewald (PME) (PME order = 4). All linear heating rates were set to 0.5 K/ns (5 Â 10 8 K/s).   11 In the case of sn-COC, given that no crystal structure was available, a crystal was built using the CG representation of sn-POP as the starting point. The sn-COC unit cell was built by reducing the number of beads on the capric fatty acid chains, moving the molecules close to each other, using the sn-POP distances as a guideline, and reducing the unit cell dimensions accordingly.

| Crystals with void generation
The simulation box at any void size and type was varied between each individual simulation. These were built by starting with a perfect crystal of 800 molecules ( Figure 1, and exemplified in Figure 2A), from which a molecule was chosen at random. Further molecules were then chosen at random until the required number of molecules was reached. In the case of "pit" voids the chosen molecules could be in any position within the crystal, with no relational proximity of the missing molecules except for those occurring through the random choice of molecules ( Figure 2B).
In the case of "crack" voids a "seed" molecule was chosen at random. The subsequent choice of molecules to remove were chosen by randomly choosing a molecule which was up to one molecule away in all three dimensions from an "anchor" molecule already chosen to be removed from the crystal. The "anchor" molecule was again chosen at random each time the next molecule to be removed was to be determined, until the required number of molecules to remove were chosen. This resulted in a range of crack voids which could be anywhere in the crystal and with any shape, however, all the chosen molecules were in close relational proximity ( Figure 2C). Void generation was done in an automated fashion by using a custom Python script (see Data S1).

| Melting point determination
Melting point onset for any given simulation was determined as follows: 1. the Potential energy over the trajectory was extracted using the gmx energy program in GROMACS, 2. the best-fit line (gradient and intercept) was calculated using linear regression on the data from 5% to 30% of the trajectory ("the sample set") versus temperature ( C), assuming that the temperature increase is perfectly linear with time, 3. the predicted Potential energy values over the whole trajectory/ temperature range were calculated using the linear regression parameters found in the previous step, 4. the sum of the square of errors for just the sample set, that is,

| TAG/crystal melting temperature determination
The melting temperature for any given TAG and void size was deter-  Comparing the plots for β 2 sn-POSt, the melting onset for pit voids was observed to be higher than that for crack voids by approximately 10 C at any given void size ( Figure 3). In this case, it is clear that the free energy barrier is higher in the cases where pit voids where used. Given these results, with the very close agreement of the simulated melting point of β 2 sn-POSt when using crack voids to the empirical melting point, as well as the fact that a crack void mimics defects found in crystals more closely than missing individual molecules, the remainder of the simulations, and results reported, were all carried out using crack voids.

| Accuracy of the determined melting point by bootstrap statistics
In a melting simulation on any given crystal with a void, given the nonequilibrium nature of the system, the determined melting point onset varies between repeats, even when using the same starting configuration. Given this, a number of simulations are needed to determine an average melting point. As more simulations are carried out, the changing calculated melting point ( Figure 4A), and its standard error ( Figure 4B), at any given void size, can be determined by using bootstrap statistics. 16,17 This statistical methodology creates numerous sub-samples by repeated resampling from a data population (in this case, using 10,000 sub-samples), thus allowing for the calculation of a population mean, standard error and confidence interval. As can be seen in Figure 4A, the calculated melting point had sta- Given the fast heating rate of the MD simulations, most crystals did not exhibit a sharp discontinuity when going from solid to melt at any given temperature and instead showed a gradual melting, starting at the void nucleus and extending towards the rest of the crystal, observed as a gradual change in the potential energy ( Figure 5A).
Given this, the melting point of that crystal is assumed to be the same as the melting onset temperature (red line in Figure 5A), that is, the temperature at which the potential energy of the system started to deviate from the expected (predicted) value (upper orange dashed line in Figure 5A). For each TAG and void size the melting points showed a distribution ( Figure 5B) Figure 7). Not only was the general trend from low to high melting point reproduced, but the determined melting points were also close to the empirically determined values, as well as the predicted values from different models.
One thing of note from these results is that the plateau was reached at different void sizes for the different TAGs (Table 1).
While the exact value will be dependent on the specific void sizes being used to determine where the plateau is reached, a general trend can be observed, namely, the void size required to reach the plateau increases with the melting point of the crystalline TAG. This is due to the increased thermodynamic stability of the crystal, and, hence, the increased need for a lowering of the melting energy barrier. This also means that the ideal void size for any given TAG cannot be known a priori but needs to be determined

| CONCLUSIONS
This study has shown that the thermodynamic melting point of pure TAG crystals can be estimated accurately. This is achieved by introducing void defects into a perfect crystal, heating the resulting structure and determining the melting point onset. Due to the non-equilibrium nature of a melting simulation, an increased number of simulations per void size are required to estimate an accurate melting point, although we find that after 150 simulations the results are sufficiently converged. All simulations were carried out using the COGITO FF, 11 which has proven to be able to reproduce the thermodynamic melting point of various TAGs, despite not being specifically parameterized against this macroscopic property.

ACKNOWLEDGMENTS
The authors thank Mondel ez International for funding this work.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are openly available in University of Strathclyde KnowledgeBase at https://doi.org/10. 15129/3cbe0feb-386e-4337-b386-4226a34247fe.