Refinement of the OPLSAA Force-Field for Liquid Alcohols.

We employ the popular all-atom optimized potential for liquid simulations, OPLSAA, force-field to model 17 different alcohols in the liquid state. Using the standard simulation protocol for few hundred nanosecond time periods, we find that 1-octanol, 1-nonanol, and 1-decanol undergo spontaneous transition to a crystalline state at temperatures which are 35-55 K higher than the experimental melting temperatures. Nevertheless, the crystal structures obtained from the simulations are very similar to those determined by X-ray powder diffraction data for several n-alcohols. Although some degree of deviations from the experimental freezing points are to be expected, for 1-nonanol and 1-decanol, the elevation of the freezing temperature warrants special attention because at room temperature, these alcohols are liquids; however, if simulated by the OPLSAA force-field, they will crystallize. This behavior is likely a consequence of exaggerated attractive interactions between the alkane chains of the alcohols. To circumvent this problem, we combined the OPLSAA model with the L-OPLS force-field. We adopted the L-OPLS parameters to model the hydrocarbon tail of the alcohols, whereas the hydroxyl head group remained as in the original OPLSAA force-field. The resulting alcohols stayed in the liquid state at temperatures above their experimental melting points, thus, resolving the enhanced freezing observed with the OPLSAA force-field. In fact, the mixed-model alcohols did not exhibit any spontaneous freezing even at temperatures much lower than the experimental values. However, a series of simulations in which these mixed-OPLSAA alcohols started from a coexistence configuration of the liquid and solid phases resulted in freezing transitions at temperatures 14-25 K lower than the experimental values, confirming the validity of the proposed model. For all of the other alcohols, the mixed model yields results very similar to the OPLSAA force-field and is in good agreement with the experimental data. Thus, for simulating alcohols in the liquid phase, the mixed-OPLSAA model is necessary for large (7 carbons and above) hydrocarbon chains.

Introduction

Alcohols are an important class of molecules with many scientific, medical, and commercial applications. Being associative liquids, they function as good solvents for chemical reactions in the lab as well as in industrial productions. Their importance has also been augmented by their utilization in the development of alternative energy sources.[1] Therefore, the ability to describe accurately their properties in computer simulations is highly desired.

The optimized potential for liquid simulation (OPLS),[2,3] and especially its all-atom version OPLSAA,[4,5] is arguably considered to be the best overall force-field to described small organic molecules in the liquid state. The main reason for its success is likely the choice to determine the partial atomic charges in a molecule empirically by reproducing the thermodynamic properties of the corresponding liquid phase, such as density, enthalpy of vaporization, and in some cases, free energy of hydration. It is interesting to mention that more elegant, and seemingly more rigorous, procedures to assign charges to a molecule, such as the restrained electrostatic potential[6] of the Amber force-field, often suffer from difficulties to reproduce thermodynamic properties of the neat liquid[7] or kinetic properties of ionic liquids.[8,9]

As with few other force-fields, it is remarkable that the OPLSAA has survived two decades of hardware and software computer improvements which resulted in few orders of magnitude longer simulation times, much larger system size, and a more accurate calculation of the long-range electrostatic forces. For example, an extensive study on the conformational populations of a homologous series of alkane chains from n-butane to n-dodecane reported a close agreement with the experimental data.[10] Nevertheless, for n-heptane and longer alkanes, the calculated heat of vaporization was substantially larger than the experimental value, suggesting overestimation of the intermolecular attractions in the liquid state. Another evidence supporting the assumption of exaggerated attractive forces between long alkane molecules comes from the transition temperature between liquid and gel phases of pentadecane. This temperature is predicted by the OPLSAA force-field to be much higher than the experimental temperature.[11]

These shortcomings and their consequence in impeding adequate description of lipids with the OPLSAA force-field have prompted Siu et al. to propose modifications to the partial charges and the Lennard-Jones parameters of the atoms in hydrocarbons.[11] In addition, the coefficients defining the potential function of several dihedral angles were also changed. The resulting set of parameters, referred to as the L-OPLS, have been shown to fix the above-mentioned problems encountered with the original force-field.

Recently, we investigated the structure and self-assembly of methanol molecules confined between two, or adsorbed on, graphene sheets.[12] In the process of extending this study for larger molecules, we tested the behavior of the bulk liquid phase (in three dimensions) at ambient conditions of a series of larger alcohols. However in this case, we observed that within few hundred nanoseconds, 1-octanol, 1-nonanol, and 1-decanol, when described by the OPLSAA force-field, spontaneously undergo a sharp first-order transition from liquid to solid. It is reasonable to assume that the excessive intermolecular attractive forces observed for long alkane chains is likely the reason for elevation in the freezing temperature. Therefore, we combined the original OPLSAA with the L-OPLS force-fields to model alcohols. The former was considered for the description of the hydroxyl head-group, whereas the latter was utilized to describe the hydrocarbon tail. We checked 17 alcohols of different types (primary and secondary, as well as, mono-, di-, and tri-alcohols) and compared their physical properties to those obtained using the original OPLSAA parameters. In addition, for 1-heptanol, 1-octanol, 1-nonanol, and 1-decanol, we studied the temperature-induced transition to the solid or glassy state for both alcohol models.

We note that there are few studies in the literature reporting refinements of the OPLSAA force-field for alkanes[13] and alcohols.[14] However, they focused on shorter aliphatic chains and apparently have not resolved the problems faced when simulating long-chain hydrocarbons.[11] In another study, Kulschewski and Pleiss reparametrized the OPLSAA force-field for alcohols, aiming to minimize the discrepancies between the calculated and experimentally determined diffusion constants.[15] To this end, they modified the partial charges of the oxygen and hydrogen atoms of the hydroxyl group and assigned different values for each alcohol molecule studied. However, besides abolishing transferability, the enthalpies of vaporization that their new models would exhibit have not been assessed. This and our inference that the elevation of the freezing temperature originates from inaccurate parameters of the alkane chain have led us to adopt the L-OPLS refinement.

Results and Discussion

In Table , we present the bulk densities of 17 alcohols considered in this study (see Figure ) for the original OPLSAA and the mixed-OPLSAA force-fields at ambient conditions of T = 298.15 K and P = 1 bar. Note that the description of methanol is identical in both models; nevertheless, we still considered it independently twice as an additional check for the convergence of the results. Examining the series of primary mono-alcohols (alkanols), as expected, the difference between the two models increases with the length of the hydrocarbon chain. Up to 1-hexanol, the difference relative to the experimental value is small, lower than 0.4%. This difference continues to increase for heptanol and octanol, but still lower than 1%. For all of the primary mono-alcohols, the OPLSAA densities are larger than their corresponding mixed-OPLSAA values, and from butanol to octanol, they are also closer to the experimental values. For nonanol and decanol, the difference between the two models jumps to 13%, where now the values derived from the mixed-OPLSAA model are by far the ones closer to the experimental values. As shown below, these large deviations of the OPLSAA model for 1-nonanol and 1-decanol are due to a phase transition to a crystalline state already at a temperature of 298.15 K. For all other alcohols, the difference between the densities of the two force-fields (relative to the experimental value) is again very small, up to 0.3% and, in most cases, the OPLSAA model is closer to the experimental data. For 2-propanol, ethylene glycol, and glycerol, both models exhibit large deviations in the density compared with the corresponding experimental values. Interestingly, the magnitude of these deviations are all around 3.5%.

Table 1

Bulk Density of the Different Alcohols Used in Our Study in kg/m3 for OPLSAA and Mixed-OPLSAA Models at T = 298.15 K and P = 1 bara

alcohol	OPLSAA	mixed-OPLSAA	expt.
methanol	775.9 ± 0.1	776.0 ± 0.1	791.4
ethanol	795.5 ± 0.1	794.3 ± 0.1	789.3
1-propanol	800.9 ± 0.1	799.0 ± 0.1	799.7
1-butanol	802.0 ± 0.1	799.3 ± 0.1	809.5
1-pentanol	807.8 ± 0.1	805.6 ± 0.1	814.4
1-hexanol	813.6 ± 0.1	810.4 ± 0.1	813.6
1-heptanol	819.5 ± 0.1	814.3 ± 0.1	821.9
1-octanol	825.4 ± 0.1	817.6 ± 0.1	826.2
1-nonanol	925.0 ± 0.1	820.6 ± 0.1	828.0
1-decanol	927.4 ± 0.3	823.2 ± 0.1	829.7
2-propanol	806.9 ± 0.1	808.1 ± 0.1	780.9
3-pentanol	814.6 ± 0.1	813.5 ± 0.1	820.3
4-heptanol	818.7 ± 0.1	815.0 ± 0.1	818.3
ethylene glycol	1077.3 ± 0.1	1075.3 ± 0.2	1113.5
1,3-propanediol	1047.1 ± 0.1	1045.3 ± 0.3	1053.8
1,2-propanediol	1034.9 ± 0.3	1035.4 ± 0.3	1036.1
glycerol	1215.9 ± 0.3	1216.4 ± 0.1	1261.3

The experimental values taken from the CRC Handbook of Chemistry and Physics[16] are also given. Error estimates lower than 0.1 kg/m3 were not considered.

Figure 1

Model alcohols investigated in this study. (a) Series of primary mono-alcohols starting from methanol until 1-decanol. (b) Secondary alcohols: 2-propanol, 3-pentanol, and 4-heptanol. (c) Primary diols: ethylene-glycol and 1,3-propanediol. (d) 1,2-Propanediol, and (e) glycerol.

The enthalpies of vaporization calculated by eq are displayed in Table (the values of Ugas and Uliquid are given in Table S4). For all alcohols studied, the value obtained by the OPLSAA force-field is larger than that of the mixed-OPLSAA model, except for 2-propanol for which both values are similar. Inspecting the series of the primary mono-alcohols up to 1-octanol, the difference in the heat of vaporization between two models is smaller than 2 kJ/mol. For 1-propanol and 1-butanol, the OPLSAA values are closer to experiments, nonetheless, even with the mixed-OPLSAA model, the relative deviations of the computed heats of vaporization from the experimental values are equal or less than 5%. For 1-pentanol and above, the mixed-OPLSAA model performs better when compared to the experimental data. We note that for alcohols larger than pentanol, the enthalpies of vaporization calculated from the simulations, using either of the force-fields, are larger than the experimental values. However, as was the case for the densities, the striking observations are the values of the OPLSAA for 1-nonanol and 1-decanol which deviates from the experimental data by 26 and 29 kJ/mol, respectively. This is easy to understand because, as mentioned above, these alcohols described by the OPLSAA force-field are solids at T = 298.15 K and therefore, their calculated “heat of vaporization” is actually the heat of sublimation, or the sum of the heat of vaporization and the heat of fusion. Experimentally,[20] the heats of fusion of 1-nonanol and 1-decanol were found to be 28.80 and 37.66 kJ/mol. Thus, if we subtract these values from the calculated “heats of vaporization” for 1-nonanol and 1-decanol, as shown in Table , we obtain 73.81 and 73.22 kJ/mol, respectively. For 1-nonanol, this vaporization enthalpy is quite similar to the experimental value; however, for 1-decanol, it differs by 8 kJ/mol. We note however that a different study[21] determined the heat of sublimation of 1-decanol at 298 K to be 112.5 kJ/mol which is very close to the value of 110.88 kJ/mol obtained by the OPLSAA force field. The OPLSAA and the mixed models display a similar performance (within a range of 3 kJ/mol) for the calculated enthalpies of vaporization for the secondary alcohols, diols, and glycerol. Although, the mixed-OPLSAA model is slightly closer to the experimental values, in most cases, the difference between the two models is much smaller than the deviation from the experimental data. Nevertheless, the agreement with the values obtained from experiments is quite good, except for 1,2-propanediol, glycerol, and 1,3-propanediol which display relative deviations of about 8%, 10%, and 17%, respectively.

Table 2

Enthalpy of Vaporization of Different Alcohols at 298.15 K and 1 bar in kJ/mol for OPLSAA and Mixed-OPLSAA modelsa

alcohol	OPLSAA	mixed-OPLSAA	expt.
methanol	36.65 ± 0.01	36.63 ± 0.02	37.43
ethanol	42.68 ± 0.05	42.51 ± 0.04	42.32
1-propanol	46.04 ± 0.06	44.76 ± 0.08	47.45
1-butanol	52.01 ± 0.09	50.69 ± 0.08	52.35
1-pentanol	57.52 ± 0.05	56.96 ± 0.04	57.02
1-hexanol	63.49 ± 0.12	62.75 ± 0.10	61.61
1-heptanol	69.85 ± 0.11	68.75 ± 0.13	66.81
1-octanol	76.31 ± 0.11	74.52 ± 0.20	70.98
1-nonanol	102.61 ± 0.14	80.63 ± 0.19	76.86
1-decanol	110.88 ± 0.24	86.31 ± 0.12	81.50
2-propanol	47.55 ± 0.05	47.70 ± 0.02	45.39
3-pentanol	54.08 ± 0.13	51.97 ± 0.08	54.0
4-heptanol	66.90 ± 0.26	64.16 ± 0.12	62.4[17]
ethylene glycol	67.08 ± 0.06	66.37 ± 0.09	63.9
1,3-propanediol	82.00 ± 0.09	80.73 ± 0.09	69.8
1,2-propanediol	73.19 ± 0.13	72.40 ± 0.12	67.5[18]
glycerol	83.90 ± 0.13	82.70 ± 0.16	91.7[19]

The experimental values are taken from the CRC Handbook of Chemistry and Physics[16] unless otherwise indicated.

Table shows the self-diffusion coefficient of different alcohols. For most of the alcohols, the value obtained from the OPLSAA force-field is smaller than that of the mixed-OPLSAA model. For the rest, the diffusion coefficient is the same (within the estimated errors) and only for 2-propanol, it is larger for the OPLSAA force-field. In comparison with the experimental results and considering the range of values reported in the literature, both models perform reasonably well for the primary mono-alcohols up to 1-pentanol, where the OPLSAA exhibits even closer values to the experimental data. From 1-hexanol and above, the mixed-OPLSAA continues to perform well; however, the OPLSAA displays retarded dynamics which is intensified with the length of the hydrocarbon tail of the alcohol. Obviously, this is expected for 1-nonanol and 1-decanol because of the transitions to crystals; however, the OPLSAA values for 1-hexanol, 1-heptanol, and 1-octanol suggest that the onset of the effect operates already at shorter alcohols. Although, we were not able to find experimental data for 1-nonanol and 1-decanol, the magnitudes of diffusion coefficients of the mixed-OPLSAA model characterize a liquid phase, as is known to be the case experimentally at T = 298.15 K, much more than the values found for the OPLSAA model. For the other alcohols, the experimental data indicate that the OPLSAA is performing slightly better than, or similar to, the mixed-OPLSAA model. Large difference between the two models is observed for 4-heptanol; however, also in this case, we could not find experimental data that prevent assessment. The alcohols displaying the largest, and substantial, disagreement with the experimental data are 3-pentanol, glycerol, and in particular, 1,2-propanediol.

Table 3

Self-Diffusion Coefficient of Different Alcohols in 10–9 m2/s at T = 298 K and P = 1 bara

alcohol	OPLSAA	mixed-OPLSAA	expt.
methanol	2.87 ± 0.07	2.80 ± 0.03	2.27;[22] 2.44[23]
ethanol	1.15 ± 0.02	1.17 ± 0.04	1.01;[22] 1.16;[23] 1.08[24]
1-propanol	0.75 ± 0.02	0.83 ± 0.01	0.646;[22] 0.590;[23] 0.627[24]
1-butanol	0.47 ± 0.01	0.56 ± 0.05	0.504;[22] 0.426;[25] 0.456[24]
1-pentanol	0.31 ± 0.01	0.36 ± 0.02	0.296;[26] 0.286[27]
1-hexanol	0.18 ± 0.01	0.24 ± 0.01	0.218[26]
1-heptanol	0.110 ± 0.004	0.162 ± 0.006	0.172[26]
1-octanol	0.066 ± 0.002	0.120 ± 0.004	0.138;[25] 0.142[26]
1-nonanol	0.0006 ± 0.0004	0.089 ± 0.001
1-decanol	0.0002 ± 0.0001	0.069 ± 0.001
2-propanol	0.558 ± 0.004	0.527 ± 0.002	0.649;[22] 0.582[23]
3-pentanol	0.331 ± 0.008	0.446 ± 0.006	0.232[28]
4-heptanol	0.080 ± 0.007	0.146 ± 0.003
ethylene glycol	0.125 ± 0.002	0.119 ± 0.009	0.0961;[29] 0.083[30]
1,3-propanediol	0.0126 ± 0.0003	0.0148 ± 0.0006
1,2-propanediol	0.0094 ± 0.0005	0.0086 ± 0.0003	0.045;[31] 0.041[31]
glycerol	0.0072 ± 0.0006	0.0091 ± 0.0001	0.0047;[31] 0.0025[32]

For few alcohols, no experimental values were found.

We also calculated the relative permittivity of the alcohols described by both models and summarized the results in Table S5. The mixed-OPLSAA exhibits very similar values (within the estimated errors) to the original OPLSAA model. This is not surprising given that the polar hydroxyl group(s) of the alcohols are unaltered between the two models. As reported previously in the literature,[33] the OPLSAA model deviates substantially relative to the experimental results. In all cases, the computed values are smaller than those determined experimentally, except for 4-heptanol and 1,3-propanediol. For the former, the modeled values reproduce well, whereas for the latter, they are larger than the experimental data.

The augmentation of the retarded diffusion with the length of the hydrocarbon chain displayed in Table for the OPLSAA force-field led us to speculate that shorter alcohols, for example, 1-octanol and 1-heptanol, might also exhibit melting or freezing temperatures which are much higher than in experiments. To examine this conjecture and to determine the freezing temperatures for 1-nonanol and 1-decanol, we performed additional series of simulations in which we started from the liquid state at high temperatures and gradually cooled down the systems. A freezing transition in three-dimensional liquid is known to occur via a discontinuous, first-order transition. This is a sharp transition in which the densities of the two phases differ substantially; therefore, in Figure , we plot the density as a function of temperature for the original OPLSAA and the mixed-OPLSAA models. For 1-decanol, 1-nonanol, and 1-octanol described by the OPLSAA force-field, there is a clear and significant jump in the density at, respectively, T ≃ 327.5, 322.5, and 294 K, suggesting a transition to a dense solid phase. At these freezing temperatures, there are also concomitant drops of the diffusion coefficients (displayed in Figure ) to very small values demonstrating transitions of the systems from fluids to solid states. The experimental melting temperatures of these three alcohols are[16] 280, 268, and 258.5 K, respectively, which means the melting temperatures predicted by the OPLSAA is about 35–55 K higher than the experimental values. For decanol and nonanol, this is critical because simulating these alcohols at ambient temperatures in which they are supposed to be liquids would transform them to solids. In contrast, for the mixed-OPLSAA model, the change of the density is linear without any signature of a phase transition. In addition, the corresponding diffusion constants also exhibit only gradual monotonic decreases with decreasing the temperature without any sharp drop. The values of these diffusion constants, as well as the shape of the mean-squared displacement (MSD) as a function of time (not shown), just above the experimental melting temperatures indicate that these mixed-OPLSAA alcohols are in the liquid state. As the temperature is decreased further, the dynamics continues to slow down; however, no first-order transition to a solid state is detected even when cooled down 60 K below the corresponding experimental melting temperatures. Therefore, it is likely that at some temperature, these alcohols transform into a supercooled liquid or a glassy state as can also be observed experimentally.[34]

Figure 2

Density of 1-decanol, 1-nonanol, 1-octanol, and 1-heptanol, as a function of temperature for the OPLSAA and the mixed-OPLSAA force-fields. In this series of simulations, adjacent temperatures were decreased or increased gradually from the simulation at T = 298 K (see Methods). The light blue vertical line denotes the experimental melting temperature of each alcohol.[16]

Figure 3

Diffusion constant of 1-decanol, 1-nonanol, 1-octanol, and 1-heptanol, as a function of temperature for the OPLSAA and the mixed OPLSAA force-fields. The light blue vertical line denotes the experimental melting temperature of each alcohol.[16]

Note that for the OPLSAA model, 1-decanol and 1-nonanol at T = 298.15 K started to display an increase in the density and a growth of the crystal nuclei already after the first 10 ns, whereas for 1-octanol at T = 280 K, it required approximately 200 ns. Obviously, one can argue that the time scale of the simulations is much smaller than the relaxation time required for a freezing transition of the mixed-OPLSAA alcohols, and that the relaxation times observed for the OPLSAA model is much faster than the actual experimental relaxation times. A first-order transition is known to exhibit hysteresis (ergodicity breaking). Thus, the transition temperature of freezing a liquid is lower than that of melting a crystal. Therefore, to assess the extent of hysteresis in these systems, we performed an additional series of simulations in which we started from the crystal structure at a low temperature and gradually (after 50 ns) heated up the system to the next sampled temperature. At each temperature, a minimum of 100 ns simulation time was conducted. Longer simulations, up to 400 ns, were performed at temperatures just below those in which a melting transition was observed. The starting structure for both models was taken from the lowest temperature configuration of the spontaneous crystallization of the OPLSAA model. For 1-heptanol where no such structure existed, the starting structure was adapted from that of 1-octanol by replacing the terminal methyl group with a hydrogen. The results are displayed in Figure indicating that the alcohols described by the mixed-OPLSAA melt at temperatures near (nonetheless, always higher than) those of the experimental melting points, with deviations in the range of 5–17 K. However, the alcohols described by the OPLSAA model melt at significantly higher temperatures with deviations in the range of 64–107 K relative to the experimental data.

Figure 4

Similar to Figure , however, in this case, the series of simulations started from a crystal structure at the lowest temperature and then gradually heated up (see text).

Figures and 4 indicate these systems display large hysteresis, and there is a need for a more accurate determination of the freezing/melting temperatures. To this end, we conducted additional series of simulations for 1-octanol, 1-nonanol, and 1-decanol where the starting configuration of the system at each temperature was a coexistence of the crystal and liquid phases.[35,36] This starting configuration, which was the same at each temperature, was taken from a corresponding simulation of the spontaneous crystallization observed above for the OPLSAA model (about halfway of the freezing process). The results are summarized in Figure , where the density of the phase formed is plotted as a function of the temperature. Furthermore, the evolution of the instantaneous density above and below the transition temperature in the case of the OPLSAA model is shown in Figure S1. These figures indicate that the freezing/melting temperature of 1-decanol, 1-nonanol, and 1-octanol using the OPLSAA model is around 335, 325, and 304 K, respectively, that is, values that are up to 10° higher than those determined based on spontaneous crystallization from the liquid state (Figure ). In contrast, the freezing/melting temperature of 1-decanol, 1-nonanol, and 1-octanol using the mixed-OPLSAA model is around 255, 255, and 245 K which are only 14–25 K lower than the experimental temperatures.

Figure 5

Similar to Figures and 4, however, in this case, the starting configuration of the system at each temperature was a coexistence (the same configuration for all temperatures) of the liquid and solids states.

Although, from Figures and 3, it is clear that 1-decanol, 1-nonanol, and 1-octanol of the OPLSAA model undergo spontaneous temperature-induced transition to a solid state, and the nature of this solid has not been revealed yet. In Figure , we present snapshots of the last configuration for these alcohols at a temperature just below the freezing transition. In these three cases, the alcohols transformed into an ordered crystal in which the long axis of the molecule is perpendicular to a layered structure. Within each layer, the alcohols form a hexagonal lattice. This is clearly displayed in Figure where the in-plane structure of only a slice of the simulation box that includes a single layer is shown. The same hexagonal structure can also be seen in Figure S2 which is a view on the entire simulation box along the long-axis of the molecules (thus, a rotation by 90° of the structure shown in Figure ). These features are in agreement with the crystal structures found by X-ray powder diffraction data for several n-alcohols including 1-butanol,[34] 1-pentanol,[37] 1-octanol,[38] and even for much longer n-alcohol chains with 11–37 carbon atoms.[39,40] In most cases, an alcohol can exhibit few crystal phases (depending on the temperature) that mainly differ in the orientational order, and the tilting degree, of the long-axis of the molecules, as well as in trans or gauche bond conformations involving the hydroxyl group atoms. Nevertheless, in all cases, a layered structure is formed by hydrogen bonds between the hydroxyl groups of different alcohols, and the long-axis of the molecules is either perpendicular or tilted with respect to the plane of the layers. In the latter, the long-axes orientations of the molecules can then form a zigzag pattern, as obtained in Figure , for 1-octanol and 1-nonanol, or a parallel arrangement. Experimentally 1-octanol[38] at 150 K and 1-pentanol[37] at 183 K exhibit a much stronger zigzag pattern than that displayed in Figure . Besides the difference in the temperatures, a plausible reason for the discrepancy is that the cubic shape of the simulation box suppresses larger tilting between the layers.

Figure 6

Snapshot of the last frame of the simulation of 1-octanol, 1-nonanol, and 1-decanol for the two investigated force-field models. In the OPLSAA case (left column), the chosen simulation is for the highest temperature for which crystallization occurs, whereas for the mixed-OPLSAA case (right column), it is for the lowest temperature studied.

Figure 7

Slice of the snapshot of 1-octanol shown in Figure that includes (for clarity) only a single layer of alcohols. The plane of the layer coincides with the plane of the image. The drawings of the two perfect hexagons are a guide for identifying the hexagonal in-plane structure within each layer. The green hexagon includes alcohols in which their hydroxyl groups point downward, whereas for the yellow hexagon, the hydroxyl groups point upward.

Nonetheless, the most substantial discrepancy between the simulation results and the experimental structure is that in the former within each layer, the hydroxyl group is pointing up and down as displayed, for example, by the yellow and green hexagons in Figure . It is clear that this does not correspond to the lowest (free) energy minimum. When all of the hydroxyl groups within each layer would point in one direction, then the adjacent layers can interact with one another via only the hydroxyl groups (head–head) and via only the alkane chains (tail–tail) in an alternate way. The relaxation times for this type of ordering is likely to be much longer than the time period for which we are able to perform the simulations, and therefore, the simulations did not result in an equilibrated crystal structure (although it is possible to identify regions in which one side of a layer has excess, whereas the opposite side is depleted, of hydroxyl groups). In Figure a, we present the radial distribution function between the carbon atoms of the alcohols indicating that the phases formed by 1-octanol, 1-nonanol, and 1-decanol exhibit long-range order and that the structures of these crystals are very similar to one another. The curves for the OPLSAA were calculated from the simulations of the spontaneous crystallization from the disordered liquid state (see Figure S3), whereas those for the mixed-OPLSAA were calculated from the simulations starting from a coexistence of liquid/solid phases. Nevertheless, the structures of the crystals obtained using the OPLSAA model are more ordered than those of the mixed-OPLSAA.

Figure 8

Radial distribution function between the carbon atoms (excluding intramolecular correlations) of 1-octanol, 1-nonanol, and 1-decanol for the two force-field models. The OPLSAA distributions were obtained from simulations in which the liquid states were spontaneously crystallized, whereas the mixed-OPLSAA distributions were taken from the simulations in which the starting conformations included a coexistence of the liquid and solid phases. In all cases, the temperature considered is the highest temperature supporting the solid phase.

Qualitatively, both models exhibit the same behavior with a temperature decrease for 1-heptanol. The curves in Figures and 3 lack any evidence for a transition to a solid state. The OPLSAA model exhibits slightly larger densities and smaller diffusion constants compared to the mixed-OPLSAA model. The difference between two models is largest around T = 220 K, and we do not know whether longer simulations would have allowed a transformation to solid. In Figure S4, we plot a snapshot of the last configuration at the lowest temperature for the two models confirming that both systems are disordered and no transition to a crystalline phase occurred. The short-range order is also evident by the fast decay of the correlations in the radial distribution function shown in Figure S5. On the basis of the results obtained for 1-heptanol, we assume that also shorter OPLSAA alcohols will not crystallize at lower temperatures and the behavior would be similar to the mixed-OPLSAA model.

To address concerns that the spontaneous crystallizations of the OPLSAA model reported above arise because of the type of barostat used[41] or, alternatively, because the random distribution of the molecules in starting configurations we prepared contained some degree of order, we performed two additional series of simulations for the OPLSAA model: the first with the Parrinello–Rahman barostat[42] (referenced to 1 bar with a time constant of 2 ps) and the second with the Berendsen barostat (as before). However, the starting configurations were taken from simulations of the corresponding gas phase simulated for 25 ns at 1000 K. In both cases, we considered the 1-nonanol and 1-decanol systems at T = 298.15 K and 1-octanol at T = 280 K. These two independent series of simulations reproduced the three spontaneous crystallization transitions from the liquid states reported above. In Table S6, we provide the densities obtained for the crystals. These densities are slightly lower than those obtained from the simulations reported in Figure (which conducted by slowly reducing the temperature), nevertheless, the crystal structures formed are very similar, as evidenced by the radial distribution function shown in Figure S6.

Conclusions

In this article, we find that the OPLSAA force-field for 1-octanol, 1-nonanol, and 1-decanol leads to spontaneous crystallization from the liquid phase at temperatures which are about 35–55 K higher than the experimental melting temperatures. Estimation of the melting temperatures starting from a coexistence of the solid and liquid phases resulted with even higher values, approximately 50 K above those determined experimentally. This overestimation of the melting temperature is especially important for 1-nonanol and 1-decanol because at ambient conditions (T ≈ 300 K), these alcohols are liquids; however, if simulated using the OPLSAA model, they will spontaneously crystallize. It is likely that this deficiency, which in the series of primary mono-alcohols intensifies with the length of the alkane chain, is closely related to the reports in the literature of the large deviations in the behavior of OPLSAA hydrocarbons from experiments. To rectify the description of these hydrocarbons, small changes have been introduced to several nonbonded parameters and dihedral angles in what is known as the L-OPLS force-field. We therefore incorporated this L-OPLS description into the original OPLSAA force-field to model various types of alcohols. These mixed-OPLSAA alcohols, including 1-octanol, 1-nonanol, and 1-decanol, remained in the liquid state at temperatures above their experimental melting points, thus, resolving the crystallization issue of the original OPLSAA model. Furthermore, simulations starting from a coexistence of the liquid and solid phases enabled the determination of their melting temperatures to be 14–25 K lower than the experimental values.

For simulation time periods considered here, a similar temperature-induced study for 1-heptanol described by the OPLSAA model did not exhibit any transition, and we speculate that shorter alcohols will not crystallize as well. For all other alcohols studied, which include a relatively short hydrocarbon chain, both models give similar results and reproduce reasonably well experimental data. Nevertheless, for short alkanols (1-pentanol or shorter) as well as for 3-pentanol and glycerol, the original OPLSAA version actually slightly outperforms the mixed model.

Except for a few cases, the OPLSAA model yields densities and enthalpies of vaporization larger, whereas diffusion constants smaller than the mixed-OPLSAA model. This is because the changes introduced in the L-OPLS parameters lead to weaker attractions between the alkane chains. The effect is demonstrated in Table S7, where we calculated the nonbonded energy between (as well as within) the alcohol molecules. In all cases, except for three, the energy derived from the mixed-OPLSAA model is larger (i.e., weaker attractive forces) than the OPLSAA model. This shifts the mixed-OPLSAA freezing temperatures of 1-octanol, 1-nonanol, and 1-decanol to lower values. In fact, we were not able to observe the spontaneous crystallization from a pure liquid phase in any of the mixed-OPLSAA alcohols examined, even when the temperature dropped to 60 K below the corresponding experimental freezing point. However, a series simulations in which these mixed-model alcohols were gradually heated up from the crystal structures at low temperatures resulted in melting transitions at temperatures 5–17 K higher than the experimentally determined melting temperatures.

We performed also simulations of only the tail segment of the alcohols, that is, the corresponding alkanes from n-heptane to n-decane. The results, as shown in Table , indicate that although the heat of vaporization is better described by the L-OPLS force field, the densities, except for decane, are closer to the experimental values when the alkanes are described by the OPLSAA force-field. Therefore, in the absence of the hydrophilic head group, for example, for alkanes, the exaggeration of the attractive interactions between the hydrocarbon molecules is diminished compared to that found in alkanols, and an obvious out-performance of the L-OPLS force field is evident only for decane. In fact, given the trends in the performance of both force-fields for the primary mono-alcohols and for n-alkanes, it might be that developing a transferable nonpolarizable force-field describing the gas, liquid, and solid states with high-fidelity is extremely difficult, if not impossible. One obvious obstacle is the fixed partial charges of the atoms in the molecule despite the different dielectric constants of its surrounding because of the different structures of these three states. In addition, the incorporation of polar groups into the molecule (such as when constructing an alcohol model from an alkane) can result in significant changes in the fixed partial charges or in induced dipoles challenging the transferability for nonpolarizable force-fields. Here, it is interesting to point out that a refinement of the OPLSAA force-field for carbohydrates found it necessary to apply scaling factors for the 1,5 and 1,6 electrostatic interactions.[43]

Table 4

Bulk Density (in kg/m3) and Enthalpy of Vaporization (in kJ/mol) for Several n-alkanes for the OPLSAA and Mixed-OPLSAA Models at T = 298.15 K and P = 1 bara

	density			ΔHvap
	OPLSAA	mixed-OPLSAA	expt.	OPLSAA	mixed-OPLSAA	expt.
heptane	679.3	673.9	679.5	40.53 ± 0.08	38.04 ± 0.04	36.57
octane	700.0	693.8	698.6	46.49 ± 0.11	43.58 ± 0.11	41.49
nonane	717.0	709.8	719.2	52.87 ± 0.14	48.59 ± 0.09	46.55
decane	731.6	722.7	726.6	59.46 ± 0.15	54.13 ± 0.12	51.42

These results were obtained from simulations with 200 ns equilibration time and 40 ns data-collection time using the same protocol as for the simulations of the alcohols. Experimental values were taken from the CRC Handbook of Chemistry and Physics.[16] Error estimates for the densities were around 0.05 kg/m3 or lower.

Interestingly, the structures of the alcohol crystals obtained in the simulations are very similar to those deducted by X-ray powder diffraction data of various n-alcohols. In particular, a layered structure was formed by hydrogen bonds between the hydroxyl groups, perpendicular or tilted orientations of the long-axis of the molecule with respect to the layer’s plane, and hexagonal arrangement of the molecules within the layers. The major discrepancy observed is that within each layer, the hydroxyl groups are pointing in both directions, whereas experimentally they point only toward one direction.

Methods

We simulated 17 different alcohols in the bulk liquid state. These model alcohols (see Figure ) include primary mono-alcohols (methanol to 1-decanol), secondary mono-alcohols (2-propanol, 3-pentanol, and 4-heptanol), diols (ethylene glycol, 1,3-propanediol, and 1,2-propanediol) and a triol (glycerol). Each system was composed of 864 molecules in a cubic-shaped box with a size ranging from 3.9 nm (for the smallest system of methanol) to 6.3 nm (for the largest system of 1-decanol). The starting configuration for each system was obtained by placing the alcohol molecules randomly in a large box and then compressing the system toward its bulk liquid density via the application of a barostat. Upon reaching a density in the vicinity of the bulk value, each system was then equilibrated for 60 ns, and subsequently data were collected for an additional 40 ns simulation time.

The molecular dynamics package GROMACS version 4.6.5[44] was utilized to perform all simulations, employing the leap-frog algorithm to integrate Newton’s equations of motion with a time step of 2 fs. Periodic boundary conditions were applied in all three dimensions. Electrostatic interactions were calculated by the PME[45] method. The cut-off distance defining the real space was 1.2 nm, and the grid spacing for the reciprocal space was 0.12 nm with quadratic interpolation. Lennard-Jones interactions were evaluated by a single cut-off distance of 1.2 nm with long-range dispersion corrections for the energy and pressure. The system was maintained at a temperature of 298.15 K by the velocity-rescaling thermostat[46] with τ = 0.1 ps and at a pressure of 1.0 bar by the Berendsen thermostat[47] with τ = 1.0 ps and a compressibility of 1 × 10–5 1/bar. All covalent bonds were described by harmonic potentials except those involving hydrogen atoms. In the latter case, the bonds were constrained using the LINCS algorithm.[48]

To calculate the enthalpy of vaporization per mole of liquidthe difference in the energy, U, between the gas and liquid states were evaluated (R is the gas constant, T is the temperature, P is the pressure, and V is the volume). Therefore, we also conducted simulations of one molecule of the corresponding alcohol in vacuum. In this case, simulations were performed at constant volume in which the cubic box had a length of 12.0 nm. Temperature thermostat was applied as above. The Lennard-Jones and electrostatic interactions were calculated by a cut-off distance of 2.8 nm, which in all cases, was larger than the size of the alcohol molecule in the box. This means that all pair interactions were calculated exactly and larger cut-off distances yielded identical results. Each system was equilibrated for 50 ns, and data were collected for additional 200 ns.

In a second series of simulations, we investigated the temperature-induced (keeping the pressure at 1 bar) phase transition from the liquid to the solid state for four primary mono-alcohols (1-decanol, 1-nonanol, 1-octanol, and 1-heptanol). In these cases, the last frames obtained from the simulations at T = 298.15 K were taken as the starting configurations for the simulations at the adjacent higher (310 K) and lower (290 K) temperatures. Then, the configurations obtained from these simulations after a minimum propagation of 100 ns were served as starting structures for the subsequent higher and lower temperatures, and so on. At each temperature, each system was equilibrated for at least 260 ns, and then data were collected for additional 40 ns. For systems in which the density did not exhibit convergence, we extended the equilibration time up to 510 ns. The errors in thermodynamic quantities obtained from the simulations were estimated using the block averaging method.[49]

The alcohol molecules were described by the all-atom OPLSAA force-field.[4] The relevant nonbonded interaction parameters are presented in Table S1. Recently, the OPLSAA parameters for hydrocarbons were refined aiming to improve the representation of long aliphatic chains.[11] This L-OPLS force-field modified several of the nonbonded parameters and dihedral angles. The modified parameters relevant to the description of the alcohol molecules considered in this study are shown in Tables S2 and S3. Therefore for comparison, all simulations were also conducted with the description of the alkane chain taken from the L-OPLS force-field, keeping the parameters of the hydroxyl group, as well as the partial charges of the first methylene group (which is covalently bonded to the hydroxyl group), the same as in the original OPLSAA force-field. Hereafter, we refer to this combined set of parameters as the “mixed-OPLSAA” force-field.

The self-diffusion coefficient, D, was obtained from a linear regression of the plot of the single-particle MSD as a function of time using Einstein’s relationwhere t is a time interval and Δr2(t) is the MSD of the particles during this time interval, calculated by the following expressionwhere N is the number of particles in the system, and the brackets on the left-hand side indicate an average over different time origins of the trajectory. The linear regression is fitted using the least-squares method ignoring the 10% segments at the beginning and the end of the MSD plot. To estimate the error of D, the fit interval is divided into two halves. The difference between the values of D obtained from each of these half fit intervals is taken as the error estimate of D.