Characterization of the Glass Transition of Water Predicted by

Molecular dynamics simulations allow detailed study of the experimentally inaccessible liquid state of supercooled water below its homogeneous nucleat...
0 downloads 0 Views 4MB Size
Download PDF
Article pubs.acs.org/JPCB

Characterization of the Glass Transition of Water Predicted by Molecular Dynamics Simulations Using Nonpolarizable Intermolecular Potentials Cara A. Kreck and Ricardo L. Mancera* School of Biomedical Sciences, CHIRI Biosciences, Curtin University, GPO Box U1987, Perth, Western Australia 6845, Australia S Supporting Information *

ABSTRACT: Molecular dynamics simulations allow detailed study of the experimentally inaccessible liquid state of supercooled water below its homogeneous nucleation temperature and the characterization of the glass transition. Simple, nonpolarizable intermolecular potentials are commonly used in classical molecular dynamics simulations of water and aqueous systems due to their lower computational cost and their ability to reproduce a wide range of properties. Because the quality of these predictions varies between the potentials, the predicted glass transition of water is likely to be influenced by the choice of potential. We have thus conducted an extensive comparative investigation of various three-, four-, five-, and six-point water potentials in both the NPT and NVT ensembles. The Tg predicted from NPT simulations is strongly correlated with the temperature of minimum density, whereas the maximum in the heat capacity plot corresponds to the minimum in the thermal expansion coefficient. In the NVT ensemble, these points are instead related to the maximum in the internal pressure and the minimum of its derivative, respectively. A detailed analysis of the hydrogen-bonding properties at the glass transition reveals that the extent of hydrogen-bonds lost upon the melting of the glassy state is related to the height of the heat capacity peak and varies between water potentials.

■

INTRODUCTION Water is the most abundant component of living organisms. Effective vitrification of water is essential for the successful long-term storage and protection of biological tissues through cryopreservation. This can be achieved by applying ultrafast cooling rates to induce the formation of glassy water at extremely low temperatures with the help of cryosolvents, thus avoiding the damaging formation and growth of ice crystals, the subsequent dehydration, and osmotic effects on cell components and membranes and indeed arresting all physical and biochemical processes.1−4 However, studying vitrified neat water experimentally is challenging, and its formation is poorly understood. As temperature decreases below the melting point (Tm) at 273 K, water enters the supercooled regime, where it undergoes anomalous increases in its isobaric heat capacity and isothermal compressibility and anomalous decreases in its density and expansivity coefficient.5 This unusual behavior has led to the theoretical proposals of a liquid−liquid critical point6 and a lower limit of metastability7 below the homogeneous nucleation temperature, TH (231 K), in supercooled water. TH is referred to as the supercooling “limit” of pure water, beyond which nucleation and thus crystallization is inevitable.5 The addition of ionic or hydrophilic solutes depresses TH about twice as much as it depresses the Tm.8 Water can reach a glassy state with the help of solute or surface effects if it is cooled very rapidly to below the glass-transition temperature (Tg). © 2014 American Chemical Society

In the glassy state the viscosity of water reaches very high values (close to 1013 poise), stopping molecular diffusion altogether and thus preventing nucleation. Because of the nonequilibrium nature of the glassy state, its physical properties depend on the method used to achieve vitrification.5 Consequently, the determination of the value of Tg has proved to be controversial. Differential scanning calorimetry (DSC) measurements of hyperquenched water estimated it to be 136 K by detecting the small increase in the heat capacity of water associated with vitrification.9,10 However, this has been reinterpreted as possibly being a prepeak preceding the true glass transition, which has been theorized to occur at the higher temperature of 165 K.9−11 Molecular dynamics (MD) annealing simulations using very high cooling rates have predicted Tg to be 188 K using the SPC/E water potential12 and 196 K using the TIP4P water potential13,14 and also confirmed the existence of the prepeak at a lower temperature, which depends on the cooling rate and aging length.12 The unique properties of water make it difficult to emulate. Nuclear quantum effects and polarization have a significant impact on many of the properties of water but are impedingly time-consuming to calculate.15 Instead, simple nonpolarizable intermolecular potentials are commonly used in classical MD Received: November 28, 2013 Revised: January 15, 2014 Published: January 27, 2014 1867

dx.doi.org/10.1021/jp411716y | J. Phys. Chem. B 2014, 118, 1867−1880

The Journal of Physical Chemistry B

Article

Figure 1. Relationship between quadrupole moment (QT), dipole moment (μ), and melting temperature (Tm) of ice Ih for all water potentials with the experimental value for the gas phase26,27 and values for liquid water calculated using density functional theory.28

Table 1. Summary of Predicted Properties by Different Water Potentialsa experiment SPCIII SPC/EIIIb TIP3PIII TIP4PIV TIP4P-EwIVb TIP4P/2005IVb TIP4P/IceIV TIP5PV TIP5P-EV six-siteVI

Tm (K)

Tmaxρ (K)

D298K (× 10−5 cm2/s)

σ300K (mN/m)

CP298K (J/mol K)

273.15 19129 215 146 232 24629 252 27229 274 27129 28929

277.15 22829 241 182 253 27329 278 29529 285 28129 28717

2.30 3.8530 2.54 5.49 3.23 2.433 2.06

71.73 54.731 63.6 52.3 59 65.231 69.3 80.131 52.6 ∼57.635 61.837

87 11332 99.6 89.9 96 92.433 101

2.78 2.834 2.236

139 13134 11717

Tm = melting point, Tmaxρ = temperature of maximum density, D = self-diffusion coefficient, σ = surface tension, and CP = heat capacity. bIndicates potentials with self-polarization correction. Roman numerals indicate number of points used in the potential. Data taken from Vega15 unless otherwise specified. a

The main weakness of the SPC/E potential was its inability to correctly predict equilibrium phase changes. The predicted phase diagram has been found to be dependent on the quadrupole moment (QT) and its ratio to the dipole moment (μ).18 With the exception of the TIP5P and TIP5P-E potentials, decreasing QT reduces the Tm of thermodynamically stable ice Ih (see Figure 1), while increasing the μ/QT ratio promotes the formation of ice II instead. The five-site potentials have the smallest QT due to the small charges used on the hydrogen atoms and lone pairs but were parametrized to reproduce the temperature of maximum density so nonetheless have a reasonable Tm.19 The QT affects molecular-orientationbased properties under both ambient conditions and in the supercooled regime but becomes more significant as thermal motion decreases.18 The SPC, SPC/E, TIP3P, and TIP5P potentials do not qualitatively reproduce the phase diagram or predict ice Ih as a stable phase, whereas the TIP4P family of potentials does.20,21 Additional investigations into water clusters have found that three-site potentials such as the SPC/E potential lead to the formation of unphysical planar trimers in both the gas22 and liquid23 phases. A comparative study of hydrophobic hydration found a correlation between the correctness of liquid density predictions and thermodynamic properties over a range of temperatures,24 with the TIP4P family of potentials performing better than the three- and five-site potentials. The

simulations due to their low computational cost and their ability to reproduce a wide range of thermodynamic, structural, and dynamical properties. 16 Parameterization of these water potentials is done using different approaches and target properties; therefore, different potentials are better at reproducing different properties of water under specific conditions. The most commonly used target properties for parametrization are the enthalpy of vaporization and liquid water density under ambient conditions. Most water potentials use between three and six sites to represent the structure and interactions of water. A Lennard-Jones interaction site is always placed on the oxygen atom, and positive charges are always placed on the hydrogen atoms. The structural difference between water potentials is the location of the negative charge(s) and the geometry of the water molecule. The six-site potential uniquely includes a LennardJones site on the hydrogen atoms.17 A recent review of the properties of the TIP3P, SPC/E, TIP4P, TIP4P/2005, and TIP5P water potentials concluded that the best overall predictions were obtained with potentials parametrized using a vaporization enthalpy that had been corrected for self-polarization (SPC/E and TIP4P/2005), rather than the standard experimental value, as well as with the four-site potentials (TIP4P family), where the negative charge is placed on the H−O− H bisector.15 The TIP4P/2005 potential was therefore found to be the best potential overall, and the TIP3P potential was the worst. 1868

dx.doi.org/10.1021/jp411716y | J. Phys. Chem. B 2014, 118, 1867−1880

The Journal of Physical Chemistry B

Article

Figure 2. Experimental49 and predicted radial distribution functions for all water potentials at 298 K for (a) O−O, (b) O−H, and (c) H−H. Insets show magnification of first peak.

given temperature.25 A summary of selected properties for these and additional water potentials of interest is provided in Table 1. A number of MD simulation studies have been conducted with the TIP4P potential to characterize the changes to the glass

TIP4P/2005 potential was again found to be the superior model. Accuracy in the prediction of the density was also found to be correlated with the ability of water potentials to form fully coordinated water clusters of at least two solvation shells at a 1869

dx.doi.org/10.1021/jp411716y | J. Phys. Chem. B 2014, 118, 1867−1880

The Journal of Physical Chemistry B

Article

Figure 3. H-bond distribution of all water potentials obtained from equilibrium simulations at 298 K. The average number of H-bonds is overlaid in black.

transition of water induced by small molecules such as DMSO14 and amides13 to try to provide a rationale for the molecular mechanism of solvent cryoprotection. However, because the quality of the predictions of structure, density, and heat capacity varies between the water potentials, the glass transition of water is likely to be influenced by the choice of potential. There is currently no evidence as to which water potential would better represent the vitrification of water under nonequilibrium simulation conditions. Therefore, we have undertaken an investigation of all water potentials listed in Table 1.

interactions, and a cutoff of 1.0 nm was used for both electrostatic and van der Waal interactions. No long-range corrections were applied to the van der Waal interactions. The SPC,41 SPC/E,42 TIP3P,32 TIP4P,32 TIP4P-Ew,33 TIP4P/2005,43 TIP4P/Ice,44 TIP5P,19 TIP5P-E,33 and six-site17 (also known as NvdE, NE6, and TIP6P) water potentials were used. Total energy, enthalpy, density, pressure, and temperature data were saved every 1 ps during heating, and a smoothing function45 was applied to the data after averaging over all 50 simulations for each water potential prior to drawing plots and conducting numerical differentiation. The heat capacity was calculated by numerical differentiation of the total energy or enthalpy versus temperature. The glass transition is indicated by a peak in the heat capacity upon heating.46−48 Further analyses of structural (radial distribution functions (RDFs), hydrogenbonding, tetrahedral order parameters) and dynamic (diffusion coefficients, viscosities) properties of the systems were undertaken using appropriate analysis modules in GROMACS 4.5.6.

■

COMPUTATIONAL METHODS Initial molecular configurations using cubic simulation boxes of 500 water molecules were generated using Packmol38 and subsequently energy-minimized and equilibrated at 360 K for 1.0 ns. MD annealing simulations in the isobaric−isothermal (NPT) ensemble (constant number, pressure, and temperature) were carried out at 1 atm of pressure. Simulations that were performed in the canonical (NVT) ensemble (constant number, volume, and temperature) used the average volume calculated from equilibrated systems using the NPT ensemble at 298 K. The equilibrium simulations at 298K were run for 5 ns with coordinates and energies saved every 25 fs and the last 4.5 ns used for analysis. The Berendsen thermostat and barostat39 were used for controlling the temperature and pressure, respectively. Temperature coupling was performed with a time constant of 0.1 ps, while pressure coupling was performed with a time constant of 1.0 ps and a compressibility of 4.5 × 10−5 bar−1. Simulated annealing was used to change the reference temperature of the thermostat for each temperature-coupling step (i.e., every 100 fs). A cooling rate of 3 × 1010 K/s was used to anneal the systems from the initial temperature of 360 down to 0 K. This was followed by heating back to 360 K at the same rate. All simulations were carried out using the MD simulation program GROMACS40 with 50 repetitions to ensure adequate sampling. Some simulations were performed using version 3.3.3 and others with 4.5.4. A comparison of identical systems run using both versions found no significant differences in results. A simulation time step of 1.0 fs was used throughout. The particlemesh Ewald method was used to calculate electrostatic

■

RESULTS AND DISCUSSION Structure and Properties at Room Temperature. To better understand the general differences between the water potentials, we performed analyses of various properties under equilibrium conditions at 298 K. Dynamic properties such as viscosities and diffusion coefficients require a large number of snapshots to be calculated accurately and were therefore analyzed under equilibrium conditions. Although such analyses are widely reported in the literature, they have not all been reported from simulations carried out under identical conditions for such a wide range of water potentials, allowing direct comparison. RDFs were calculated at 298 K for all water potentials and are presented in Figure 2, with all potentials overlaid on the same graphs. Traditionally comparative RDFs are presented offset from one another, but this makes it difficult to perceive subtle differences between them. All water potentials exaggerate the height of the first peak of the O−O and O−H RDFs with respect to experiment.49 The TIP4P/Ice potential in particular overestimates the heights of these peaks and the depth of most troughs and displays the highest degree of longest range order. The six-site potential displays the next highest degree of order. 1870

dx.doi.org/10.1021/jp411716y | J. Phys. Chem. B 2014, 118, 1867−1880

The Journal of Physical Chemistry B

Article

Figure 4. (a) Relationship between the angle (Sg) and distance (Sk) tetrahedral order parameters for all water potentials at 298 K. The rescaled q axis is also provided to facilitate comparison with literature values. (b) Relationship between the corrected self-diffusion coefficient (Dinfinity) and the viscosity (η) at 298 K for all water potentials. The experimental values for water were taken from refs 56 and 57.

Information. The main contribution to the average number of H-bonds is, as expected, the proportion of four-coordinated species, which increases primarily at the expense of three- and two-coordinated species. The quantity of five-coordinated species is reasonably consistent across all water potentials at 8−10%, with only the TIP5P and TIP5P-E potentials deviating significantly with values of 5 to 6%. This results in an unusually low average number of H-bonds for the five-site potentials. The tetrahedral order parameter of Chau and Hardwick54 was used to analyze the degree of tetrahedrality of water for each potential. This order parameter has two components: an angle order parameter (Sg), which measures the deviations of the cosines of all interbond angles from perfect tetrahedral values of −1/3, and a distance order parameter (Sk), which measures the variance between the radial distances, with both taking a value of 0 for a perfect tetrahedral arrangement. Sg has been rescaled to produce a value of 1 for a perfect tetrahedral arrangement in the alternative parameter q.55 The values of all three order parameters are presented graphically in Figure 4a and numerically in the Supporting Information. At 298 K, the TIP4P/Ice potential is the most ordered according to both Sk and Sg. Considering both Sg and Sk parameters, the level of order then follows the sequence six-site > TIP4P/2005 > TIP4P-Ew > TIP4P > SPC/E. TIP5P-E and TIP5P then follow in decreasing distance order, while the SPC and TIP3P potentials follow in decreasing angle order. The self-diffusion coefficient (D) at 298 K was calculated from the slope of the plot of the average mass-weighted mean square displacement against time using Einstein’s relation: D = (r − r0)2/6t. The measured DL is strongly affected by the simulation box length L. The corrected D∞ (for an infinitely sized box) was then calculated according to eq 1:58,59

This is to be expected as these two water potentials were parametrized to study the properties of ice and liquid water near the melting point. Most water potentials incorrectly predict the first peak in the O−H RDF to be higher than the second, with the exception of the SPC, TIP3P, TIP5P, and TIP5P-E potentials, and all potentials predict the first peak to be much narrower than experiment. The TIP3P potential has the lowest degree of order, particularly after a radius of 0.4 nm. Although it provides the closest match to the experimental height of the first O−O and first two O−H peaks, it significantly underestimates the heights of all H−H peaks and the depth of all troughs in the O−O, O−H, and H−H RDFs. SPC also predicts a low degree of long-range order in the O−O RDF. The TIP5P and TIP5P-E potentials produce very similar RDFs to each other, with the TIP5P potential displaying a slightly higher degree of order. They both exaggerate the first H−H peak greater than every other potential except the six-site and predict an early and very narrow second H−H peak. The SPC/E potential and the TIP4P family of potentials produce similarly shaped RDFs that follow the experimental plots reasonably well. The main deviation from experiment and between the potentials is the height of the first O−O and O−H peaks. However, the TIP4P-Ice potential exaggerates the size of all peaks and troughs in the O−O RDF markedly. Therefore the SPC/E potential and the TIP4P family of potentials (excluding TIP4P/Ice) appear to reproduce the overall RDFs of water better than the others, even though they incorrectly predict the height of the first O−O and O−H peaks. Hydrogen-bonding was analyzed using a hydrogen−acceptor cutoff of 2.5 Ǻ and a hydrogen−donor−acceptor angle cutoff of 37° to achieve a donor−hydrogen−acceptor angle range of 130−180°.50−53 The average number and the distribution of hydrogen-bonds (H-bonds) per water molecule for all water potentials are shown in Figure 3. The water potentials have been ordered on the basis of increasing proportion of four-coordinated water molecules. This order is similar to the order of the heights of the first peak in the O−O RDFs and is similar to the order of the magnitude of the quadrupole moments (Figure 1). It is interesting to note that the average number of H-bonds for all water potentials is correlated with QT (R2 = 0.8834), and the corresponding graph has been provided in the Supporting

D∞ = DL +

2.837kBT 6πηL

(1)

The shear viscosity (η) was calculated using the Einstein relation over 15 ps at the beginning of the plateau in the autocorrelation function of the pressure tensor, as the accuracy of this autocorrelation function decreases with time.60 D∞ and η are plotted in Figure 4b with numerical values provided in the Supporting Information. The viscosity values can only be regarded as estimates due to the effect of the pressure coupling. 1871

dx.doi.org/10.1021/jp411716y | J. Phys. Chem. B 2014, 118, 1867−1880

The Journal of Physical Chemistry B

Article

Figure 5. Heat capacity upon warming in the (a) NPT and (b) NVT ensembles for all water potentials.

Table 2. Predicted Properties Related to the Glass Transition for All Water Potentialsa

a

SPC

SPC/E

TIP3P

TIP4P

TIP4P-2005

TIP4P-Ew

TIP4P-Ice

TIP5P

TIP5P-E

six-site

TgNPT (K) TgNVT (K) ΔTg TmaxCP (K)

174 174 0.03% 204

196 192 −2.3% 226

155 162 4.6% 185

195 185 −5.6% 226

212 196 −8.3% 242

208 196 −6.4% 238

230 204 −13% 260

237 235 −0.65% 264

236 228 −3.4% 263

241 216 −12% 267

TmaxCV (K)

207

224

192

220

230

231

243

263

257

251

ΔTmaxC ΔmaxC Tminρ (K) Tmaxρ (K) ρTg (g/cm3)

1.2% −0.32% 183 227 0.980

−1.2% −4.6% 204 251 0.986

3.6% 7.5% 167 205 1.017

−2.6% −9.9% 196 258 0.961

−4.9% −20% 213 280 0.943

−3.0% −13% 209 275 0.953

−6.5% −28% 231 303 0.931

−0.31% −0.47% 244 285 0.957

−2.3% −9.5% 241 286 0.962

−6.0% −40% 240 299 0.920

ρmin (g/cm3) ρmax (g/cm3) ρ298K (g/cm3) QT (DÅ) μ/QT (Å−1)

0.979 0.999 0.969 1.969 1.155

0.985 1.003 0.991 2.035 1.155

1.016 1.027 0.976 1.721 1.363

0.961 0.997 0.985 2.147 1.014

0.943 0.989 0.987 2.297 1.004

0.953 0.989 0.986 2.164 1.073

0.931 0.980 0.981 2.434 0.996

0.957 0.976 0.974 1.560 1.460

0.963 0.993 0.992 1.560 1.460

0.920 0.990 0.991 2.331 0.811

Equilibrium density used in the NVT simulations is ρ298K. All other parameters were measured during heating.

from 54 to 55 J/mol K for the three- and four-site potentials and 56 to 57 J/mol K for the five- and six-site potentials. All potentials predict a CP greater than the expected equilibrium value of 75 J/mol K at 300 K, whereas all potentials except TIP3P overestimate CV at 300 K. When this is taken into account, the shape and position of the right-hand side of the CP peak for the TIP4P potential displays the most similarity to experimental results for supercooled water.62 The quantitatively incorrect heat capacity values are due to the classical nature of the water potentials.63 In the NPT ensemble, the number of interaction sites in the water potential appears to be related to position and height of the glass-transition peak. The glass transitions in the NVT ensemble occur in a similar order, but the heat capacity peaks have different heights and slopes, with the most striking changes occurring in the six-site and TIP4P family of potentials. The height of the peak for the TIP3P potential increased in the NVT ensemble, whereas all others decreased, although for the SPC and TIP5P potentials this decrease is insignificant. The slope of the peak for the TIP3P potential is also the only one to significantly increase in the NVT ensemble, with the others either

However, there is still reasonable agreement with literature values for both D and η. A summary of the wide degree of variation within the literature for the SPC, SPC/E, TIP4P, and TIP4P/2005 potentials is available.61 The six-site potential provides the closest match to the experimental water values of both D∞ and η, followed closely by TIP4P/2005. The TIP4P-Ew and SPC/E potentials also perform reasonably well, with the uncorrected DL corresponding with the experimental value; however, they slightly overestimate D∞ and underestimate η. The TIP5P, TIP5P-E, TIP4P, SPC, and TIP3P potentials increasingly overestimate D∞ and underestimate η, while the TIP4P-Ice potential underestimates D∞ and overestimates η. Glass Transitions. The heat capacity plots upon heating in both the NPT and NVT ensembles are presented in Figure 5. All potentials anomalously predict initial CV and CP values of 50 J/mol K at 0 K and very similar heat capacities at the beginning of the glass transition, with CP ranging from 56 to 57 J/mol K for the three- and four-site water potentials and 60 to 61 J/mol K for the five- and six-site potentials. The corresponding CV values range 1872

dx.doi.org/10.1021/jp411716y | J. Phys. Chem. B 2014, 118, 1867−1880

The Journal of Physical Chemistry B

Article

Figure 6. Plots of density against temperature on cooling (dashed lines) and heating (solid lines) for (a) three- and five-site potentials and (b) four- and six-site potentials. Experimental values originally from Hare and Sorensen72 and fitted as reported by Ashbaugh et al.24 are given as a black dotted line.

between room temperature and halfway between the Tminρ and Tmaxρ. A density minimum at ∼210 K has been measured by small-angle neutron scattering in confined supercooled D2O,66 whereas the experimental maximum of water occurs at 277 K. The TIP4P-Ew and TIP4P/2005 potentials reproduce both the experimental Tminρ and Tmaxρ well upon heating, although during cooling their Tminρ values occurred at ∼190 K. The cooling Tminρ and Tmaxρ values for the TIP4P/2005 potential are consistent with other simulation values measured at equilibrium;67 however, the density range is reduced. Tmaxρ values for all potentials are consistent with equilibrium values reported in the literature,24,25 with the exception of the six-site potential, for which a Tmaxρ value calculated using PME for long-range electrostatics could not be found. However, the ρmax values are slightly lower than most found in the literature,24,29 and we suspect that this is due to a lack of long-range van der Waals corrections in our simulations. Between the Tminρ and Tmaxρ, structural rearrangements occur, which have been interpreted as either a collapse in overall tetrahedral structure with a homogeneous model of water or the dominant structure shifting from predominantly more tetrahedral, low density liquid (LDL) local structures to predominantly less tetrahedral, high density liquid (HDL) local structures in an inhomogeneous mixture model of water.68,69 Only the six-site and the TIP4P family of potentials predict densities at 0 K to be lower than the maximum density, as predicted by Raman and FTIR spectroscopy of confined water.70 It is due to the use of classical statistical mechanics that the densities continue to increase at very low temperatures.71 To investigate the correlation between Tminρ and Tg further, the derivative of the density with respect to temperature (ρ′) was calculated. Figure 7a shows the relative positions of the CP, ρ, and ρ′ peaks for the TIP4P potential on an arbitrary scale. Similar plots for the other water potentials can be found in the Supporting Information. The maximum of ρ′ coincides with the maximum in CP. The thermal expansion coefficient (αP) can be calculated from ρ′ according to eq 2:

decreasing or not changing significantly as with the heights. The slope of the linear section preceding the glass transition decreased slightly in the NVT ensemble for all potentials. Approximate values of Tg were calculated by finding the intersection of two lines fitted to the linear sections before and after the increase in heat capacity. This is presented in Table 2. The order with which the glass transitions occur is very similar to the order of the melting temperatures provided in Table 1. Although the values vary between the NPT and NVT ensembles, the order is mostly unaffected, with the exception of the six-site potential. The change in the predicted Tg from the NVT to the NPT ensemble appears to be correlated to the change in the magnitude (R2 = 0.9057) and position (R2 = 0.9795) of the maximum in CP/CV, with only the TIP3P potential decreasing in all of these parameters. With the exception of the five-site potentials, all of these absolute and comparative parameters are also correlated with QT. The CV peak maximum occurs 29−39 K above the Tg, and the CP maximum occurs 27−32 K above Tg. The heat capacity and other response functions such as the isothermal compressibility are known to diverge around the same temperature at a given pressure known as the Widom line. At atmospheric pressure, the Widom line lies just below TH, and the experimental temperature of maximum CP (TmaxCP) is estimated to occur at ∼225 K.64 The Tg predicted in the NPT ensemble is also strongly correlated with the temperature of minimum density (Tminρ) upon heating (R2 = 0.9830), with the density at both temperatures being nearly equivalent for all potentials. Again, except for the five point potentials, there is a correlation between all of these properties and QT. There is also a correlation between QT and the difference between Tg and the temperature of minimum density (Tminρ). For the six-site potential and the TIP4P family of potentials, Tminρ and Tg differ by only ∼1 K. Equilibrium simulations of TIP4P/2005 have also noted that the system became glassy near the Tminρ.65 The densities upon cooling and heating for all water potentials are shown in Figure 6. Simulation with all potentials results in hysteresis between cooling and heating and small differences between local temperatures of maximum and minimum densities, with the general shape of the plots consistent between potentials with the same number of interaction points. Equilibrium simulations were also performed at selected temperatures for each water potential (results not shown) and were found to follow the nonequilibrium cooling curves closely

αP =

1 dV 1 dρ ≈− V dT ρ dT

(2)

Plots of αP for all water potentials are provided in Figure 7b. Although the positions of the peaks in αP consistently follow those of the CP peaks, there is no relationship between the heights of the peaks. A minimum in αP near the Widom line at 1873

dx.doi.org/10.1021/jp411716y | J. Phys. Chem. B 2014, 118, 1867−1880

The Journal of Physical Chemistry B

Article

Figure 7. Plots of (a) the derivative of the density with respect to temperature upon heating for the TIP4P potential and (b) the thermal expansion coefficient for all water potentials.

Figure 8. Plots of the derivative of the internal pressure with respect to temperature upon heating for (a) the TIP4P potential and (b) for all water potentials.

∼225 K has been measured experimentally in confined water.70,73 This minimum is thought to occur when the dominant local structure of water changes between LDL and HDL. This would then also be responsible for the behavior of other response functions such as the heat capacity, as a large amount of energy would be required to distort and break H-bonds as the structure changes. As density remains constant in the NVT ensemble simulations, the variation of the internal pressure (P) with respect to temperature was measured instead. Figure 8a shows that the P of the TIP4P potential in the NVT ensemble is almost the mirror image of the ρ predicted in the NPT ensemble; therefore, its derivative (P′) is inverted. The Tg and TmaxCv now approximately correspond to the maximum of P (TmaxP) and minimum of P′ (TminP′), respectively. Figure 8b presents the plots of P′ for all water potentials. Interestingly, this time, the TIP5P, TIP5P-E, and six-site potentials are separated from the rest. A comparison of the density and pressure data is provided in Table 3. Annealing simulations using the same cooling rate in the NVT ensemble have previously predicted the Tg of the SPC/E potential to be 188 K,12 in reasonable agreement with our results. Equilibrium simulations have predicted the maximum in the CP peak of TIP5P to occur at ∼245 K at 0 atm74 and the maximum in

the isothermal compressibility of TIP4P/2005 to occur at ∼232 K.65 These equilibrium response functions diverge about 10−20 K lower than those calculated during annealing. Effect of the Rate of Annealing on the Glass Transition. The TIP4P and TIP5P-E potentials were arbitrarily chosen to investigate the effect of the rate of annealing on the predicted glass transition. Simulations were repeated with speeds 10 times faster and 10 times slower than the original rate (3 × 1010 K/s) for both the cooling and heating steps. The heat capacity plots for the TIP5P-E and TIP4P potentials are presented in Figure 10a,b, respectively, and those for their densities are presented in Figure 10c,d, respectively. A higher annealing rate shifts the heat capacity peak to higher temperature, decreases the slope, and slightly increases the baseline heat capacity preceding the peak. As the annealing rate decreases by a factor of 10, Tg decreases by 6−10 K for the TIP4P and TIP5P-E potentials. A similar shift in Tg is also seen when comparing simulations in the NPT ensemble to the NVT ensemble for these potentials. This suggests that using the NVT ensemble, rather than the NPT ensemble, may be a computationally efficient method of improving the calculation of the glass-transition temperature equivalent to the much more computationally expensive method of decreasing the annealing rate by a factor of 10. 1874

dx.doi.org/10.1021/jp411716y | J. Phys. Chem. B 2014, 118, 1867−1880

The Journal of Physical Chemistry B

Article

Table 3. Comparative Analysis between Points on the Density/Pressure and Heat Capacity Plots for All Water Potentials SPC

SPC/E

TIP3P

TIP4P

TIP4P/2005

TIP4P-Ew

TIP4P/Ice

TIP5P

TIP5P-E

six-site

Tminρ TgNPT difference TminαP

183 174 9 201

204 196 8 225

167 155 13 186

196 195 1 225

213 212 1 241

209 208 1 238

231 230 1 259

244 237 7 265

241 236 5 263

240 241 −1 267

TmaxCP

205

225

187

225

241

238

260

264

263

267

difference TmaxP TgNVT difference TminP′ TmaxCV

−4 178 174 4 200 207

0 196 192 4 221 224

−1 167 162 5 185 192

0 186 185 1 212 220

0 197 196 1 225 230

0 194 196 −2 225 231

0 202 204 −2 236 243

0 239 235 4 259 263

0 227 228 −1 254 257

0 211 216 −5 245 251

avg = −0.4 R2 = 0.9814

difference

−7

−3

−7

−8

−6

−5

−7

−5

−4

−5

avg = −5.7

R2 = 0.9831 avg = 4.5 R2 = 0.9987

avg = 1.0 R2 = 0.9965

Figure 9. Plots of the heat capacity and density for different annealing rates for the TIP5P-E potential (a,c) and the TIP4P potential (b,d).

Hydrogen-Bonding Analysis in the Glass Transition. The average number and distribution of H-bonds were also calculated at specific temperatures during the heating step in the NPT ensemble and are presented in Figures 11 and 12, respectively. These temperatures were chosen to be 10 K below Tg (Tg−10), Tg, Tminρ, Tmaxρ′, and Tmaxρ. All snapshots within ±0.5 K of the desired temperature were used in the analysis. Figure 10a reveals that the number of H-bonds decreases with increasing temperature and that this is dependent on the number of interaction points in the water potential. Importantly, the height of the CP peaks appears to correlate with the net loss of H-bonds between Tg/Tminρ and Tmaxρ. This is to be expected as

Plots of the density upon both heating and cooling are provided in Figure 10c,d, revealing significant hysteresis. The hysteresis in the density and the final density in the NPT ensemble both increase with increasing speed. This reflects a greater inability for the water molecules to relax into tetrahedral conformations at higher annealing rates. While the Tminρ for the cooling and heating steps differ and are affected by the annealing speed, the Tmaxρ remains reasonably steady until very high speeds are reached. Equilibrium simulations (results not shown) confirm that the position and magnitude of the density maximum are unaffected by the standard annealing speed for all potentials studied. 1875

dx.doi.org/10.1021/jp411716y | J. Phys. Chem. B 2014, 118, 1867−1880

The Journal of Physical Chemistry B

Article

Figure 10. Hydrogen-bonding analysis at specific temperatures during the heating stage for all water potentials: (a) average number of H-bonds and (b) average number of O−O interactions between first neighbors (