当前位置:首页 >> 化学 >>

Metastable liquid–liquid transition in a molecular model of water


LETTER

doi:10.1038/nature13405

Metastable liquid–liquid transition in a molecular model of water
Jeremy C. Palmer1, Fausto Martelli2, Yang Liu1{, Roberto Car2, Athanassios Z. Panagiotopoulos1 & Pablo G. Debenedetti1

Liquid water’s isothermal compressibility1 and isobaric heat capacity2, and the magnitude of its thermal expansion coefficient3, increase sharply on cooling below the equilibrium freezing point. Many experimental4–8, theoretical9–11 and computational12,13 studies have sought to understand the molecular origin and implications of this anomalous behaviour. Of the different theoretical scenarios9,14,15 put forward, one posits the existence of a first-order phase transition that involves two forms of liquid water and terminates at a critical point located at deeply supercooled conditions9,12. Some experimental evidence is consistent with this hypothesis4,16, but no definitive proof of a liquid–liquid transition in water has been obtained to date: rapid ice crystallization has so far prevented decisive measurements on deeply supercooled water, although this challenge has been overcome recently16. Computer simulations are therefore crucial for exploring water’s structure and behaviour in this regime, and have shown13,17–21 that some water models exhibit liquid–liquid transitions and others do not. However, recent work22,23 has argued that the liquid–liquid transition has been mistakenly interpreted, and is in fact a liquid–crystal transition in all atomistic models of water. Here we show, by studying the liquid–liquid transition in the ST2 model of water24 with the use of six advanced sampling methods to compute the free-energy surface, that two metastable liquid phases and a stable crystal phase exist at the same deeply supercooled thermodynamic condition, and that the transition between the two liquids satisfies the thermodynamic criteria of a first-order transition25. We follow the rearrangement of water’s coordination shell and topological ring structure along a thermodynamically reversible path from the low-density liquid to cubic ice26. We also show that the system fluctuates freely between the two liquid phases rather than crystallizing. These findings provide unambiguous evidence for a liquid–liquid transition in the ST2 model of water, and point to the separation of time scales between crystallization and relaxation as being crucial for enabling it. Although several recent investigations using free-energy methods designed specifically to study phase transitions25 have shown that the ST2 model of water undergoes a liquid–liquid transition17–19, other investigations22,23 involving seemingly identical simulations using the same model found only a single liquid and a crystalline phase and concluded that what in reality is a crystallization transition had been mistaken for a liquid–liquid transition. Because there are stringent thermodynamic conditions that a first-order transition must satisfy, it is possible, albeit computationally expensive, to definitively verify or falsify the existence of a liquid–liquid transition. To this end we use six state-of-the-art freeenergy methods (four of which are documented in Methods) and scaling analysis, and we construct a thermodynamically reversible path between the liquid and crystalline phases of the ST2 model. Figure 1a, b shows perspective and orthographic projections of the freeenergy surface for ST2 water at 228.6 K and 2.4 kbar, as a function of density and an order parameter27, Q6, that can distinguish crystalline states from configurations lacking long-range order. It can be seen that two disordered (low-Q6) phases of different density are in equilibrium (same free energy) with each other, both of them being metastable with respect
1

to the crystal phase, the latter having much lower free energy. To our knowledge, this is the first time that two metastable liquid phases in equilibrium with each other and a third, stable crystalline phase have been identified in a pure substance at the same temperature and pressure in a computer simulation. The system-size dependence of the free-energy barrier separating coexisting phases is a stringent test of the presence of a true first-order transition in computer simulations22,23,25. We have calculated the free-energy surface in the low-Q6 region corresponding to the two liquids for system sizes N 5 192, 300, 400 and 600, with Fig. 2a showing that the corresponding barrier heights satisfy the N2/3 scaling characteristic of firstorder transitions. This scaling is a consequence of the surface free energy increasing as the interface between the liquids grows with system size22,23,25. For the range of system sizes examined, the interface manifests itself through the formation of water clusters with local environments characteristic of each distinct liquid phase. Figure 2b shows an example of this behaviour in a water configuration taken from a simulation performed near the barrier region for N 5 600. Because non-zero average Q6 values in an amorphous phase arise solely from fluctuations in finite systems, this quantity must also exhibit a system-size dependence22,27. Figure 2c shows that Q6 vanishes as N21/2, in agreement with the theoretical expectation22,27, thereby confirming the amorphous character of the low-density liquid (LDL) phase. Figure 3 shows the free-energy surface computed from standard Monte Carlo (MC) simulations at fixed temperature and pressure, during which each system was sampled for 100 relaxation times without the imposition of any constraint. Over the course of these long simulations, the systems fluctuate between the two coexisting phases enough times so as to allow the calculation of the free-energy surface, which is in excellent agreement with the results shown in Fig. 1b for the low-Q6 region and also with those obtained from the four other sampling techniques (see Methods and Extended Data Fig. 1). During this time, Q6 remains invariably in the amorphous region and the systems show no evidence of crystallization. The LDL phase exhibits slow dynamics, and proper scrutiny of a metastable state requires sampling to occur over times that comfortably exceed the system’s structural relaxation time, while being significantly shorter than the nucleation time of the stable ice phase. The latter’s density is very similar to that of LDL (Fig. 1). Ice nucleation, should it occur, takes place within LDL26, rather than from the high-density liquid (HDL). The inset to Fig. 3 shows the relaxation of fluctuations in density and structural order (Q6) in LDL. It can be seen that, after an intermediate transient period during which these processes are separated by as much as two orders of magnitude, fluctuations in density and Q6 decay on very similar time scales. As documented in Methods, this is a general result, but the transient behaviour is sensitive to the particular algorithm used to sample configurations (Extended Data Fig. 2). The results of Fig. 3 confirm that the LDL phase is a properly equilibrated liquid, and that under the conditions investigated here, the characteristic time for nucleation of the stable ice phase is much longer than the structural relaxation time for the LDL phase in the ST2 model of water. The key role of kinetics in stabilizing the liquid–liquid transition is further emphasized by the fact that

Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544, USA. 2Department of Chemistry, Princeton University, Princeton, New Jersey 08544, USA. {Present address: Air Products and Chemicals Inc., Allentown, Pennsylvania 18195, USA. 1 9 J U N E 2 0 1 4 | VO L 5 1 0 | N AT U R E | 3 8 5

?2014 Macmillan Publishers Limited. All rights reserved

RESEARCH LETTER
a a
10 8 βΔF( ρ) 6 4 2 192 300 600 400

b

80

0 20 30 40 50 60 70 80

60 βΔF

c 0.08
0.06 ?Q6? 0.04 0.02

N2/3

40 0.5 0.4 0.3 0.2 0.1
Q6

192 300 400 600

20 0 1.250 1.175 1.100 ρ( g c 1.025 0.950 m –3 0.875 )

0

0

0.02 0.04 0.06 0.08 0.10 N–1/2

b

1.250

1.175

1.100 ρ (g cm–3)

1.025

Figure 2 | Finite-size scaling and the liquid–liquid interface. a, Free-energy barrier separating the HDL and LDL basins computed at coexistence for systems containing N 5 192, 300, 400 and 600 ST2 water molecules. The barrier height increases with system size, obeying the N2/3 scaling law expected for a first-order phase transition. Error bars were computed using the bootstrap analysis described in Methods. b, Large clusters are formed near the barrier region by water molecules with local coordination environments characteristic of HDL and LDL (blue and red molecules, respectively). The local structure index order parameter, I, described in the text and Methods was used to characterize each molecule’s local environment, with blue molecules (HDL) ? 2. The green ? 2 and red molecules (LDL) having I . 0.12 A having I # 0.12 A simulation box containing 600 ST2 water molecules has been replicated across its periodic boundaries to illustrate that the clusters span the length of the unit cell. c, The mean value of Q6 averaged over the LDL basin decreases with system size, scaling as N21/2, and confirming the disordered nature of the liquid. The symbol size is larger than the estimated uncertainty for ?Q6?.

0.950

0.875

0.1

0.2

0.3 Q6

0.4

0.5

Figure 1 | Thermodynamic equilibrium between metastable liquid polymorphs. a, Reversible free-energy surface (F 5 free energy, b 5 1/kBT) described by density and the crystalline order parameter, Q6, for 192 ST2 water molecules at a point of liquid–liquid coexistence (228.6 K and 2.4 kbar). b, An orthographic projection of the free-energy surface shown in a. The HDL and LDL basins (r < 1.15 g cm23 and r < 0.90 g cm23, respectively) located at Q6 < 0.05 are separated by a ,4kBT free-energy barrier and are metastable with respect to cubic ice (Q6 < 0.52, r < 0.90 g cm23) by ,75kBT at this temperature and pressure. The average uncertainty in the free-energy surface is less than 1kBT. Contours are 1kBT apart.

the system fluctuates freely between the two liquid basins in unconstrained simulations (Fig. 3) without crystallization, even though the barrier separating the two liquids is comparable to that separating LDL and ice (Fig. 1). Figures 1 and 3 show that the metastable liquids are not distinguished by Q6 because of their amorphous nature, suggesting that other order parameters must be used to characterize their structure. The local structure index28 (LSI) is an order parameter that quantifies the extent to which a molecule possesses a tetrahedral environment with well-separated first and second coordination shells. Figure 4 shows the free-energy surface of
3 8 6 | N AT U R E | VO L 5 1 0 | 1 9 J U N E 2 0 1 4

ST2 water at 228.6 K and 2.4 kbar for N 5 192 projected onto the space parameterized by the first moment of the molecular LSI distribution,  I , and Q6. Water molecules within the HDL phase have a disordered coordination structure, resulting in  I < 0, because of the presence of interstitial molecules residing between the first and second neighbour shell that disrupt local tetrahedral order. The coordination structure of LDL ? 2) is more ordered, with two distinct neighbour shells that ( I < 0.15 A give rise to its ice-like density. The LSI parameter also distinguishes the ? 2) with its well-defined coordination structure that ice phase ( I < 0.25 A exhibits long-range, crystalline order. The inset to Fig. 4 shows that the changes in the coordination structure along the HDL–LDL and LDL– crystal paths are accompanied by large topological rearrangements des . The average cribed by the first moment of the ring size distribution, R ring size decreases monotonically along the HDL–LDL path, suggesting a continuous rearrangement process. In contrast, abrupt, non-monotonic behaviour is observed along the transition from LDL to the crystal in the vicinity of the saddle point in the  I –Q6 free-energy surface, which is consistent with structural rearrangements that have been observed in ice nucleation trajectories taken from long molecular-dynamics simulations of the TIP4P water model29. We note, however, that the system size examined here, although suitable for accurate free-energy calculations, may be insufficient to provide information about the mechanisms governing the ice nucleation and growth process. Such behaviour should therefore be investigated in future studies using larger systems. Our free-energy calculations demonstrate that the ST2 model of water exhibits a liquid–liquid phase transition under deeply supercooled conditions. An emerging question is to understand which aspects of intermolecular interactions cause some water models to undergo a liquid–liquid transition with well-separated relaxation and crystallization times, whereas

?2014 Macmillan Publishers Limited. All rights reserved

LETTER RESEARCH
1.250

1.175

other models do not show this behaviour. The present results suggest that constraints associated with the breaking and forming of hydrogen bonds, present in ST2 (ref. 24) but not in coarse-grained models13, have an important role. Further research using state-of-the-art free-energy methods, such as those employed here, can provide insights into this question and may thereby also improve our understanding of the phase behaviour of real water under deeply supercooled conditions.
1.0 0.8 C(t) 0.6 0.4 0.2 0 103 104 105 106 107 t (MC moves) 108

1.100 ρ (g cm–3)

METHODS SUMMARY
The reversible free-energy surface described by density, r, and the bond-orientational order parameter, Q6, was computed for the Ewald-compatible variant of the ST2 water model24 described by Liu et al.17 using MC simulations in the isothermal– isobaric ensemble, augmented with collective, N-particle rotational and translational MC moves and umbrella sampling30. A harmonic umbrella bias potential was used to restrict each MC simulation to a different window in r–Q6 parameter space. Each simulation used to generate Fig. 1 was equilibrated for ,104tQ6 , followed by a production phase of equal or greater duration, where tQ6 is the integrated autocorrelation time30 for Q6 in the sampling window. Two-dimensional r –Q6 histograms were generated from uncorrelated samples collected in each umbrella window. The histograms were subsequently combined30 to produce an unbiased estimate of the free energy. Special care was taken to ensure reversibility in the low-density region (r , 0.98 g cm23), enhancing sampling of degrees of freedom associated with structural order by performing Hamiltonian exchange MC moves30, in which umbrella restraint parameters were swapped between simulations in adjacent windows along Q6. Bi-directional sampling was also performed in this region, seeding two separate generations of simulations with initial configurations extracted from a freezing (LDL R Ice Ic) or melting (Ice Ic R LDL) trajectory. Reversibility was explicitly checked by comparing histograms from each generation of simulations to monitor for hysteresis (path dependence). Saved simulation trajectories were analysed to examine the structural and topological properties of each phase identified in the free-energy surface. The final data sets were subjected to critical scrutiny and were found to be insensitive to the sampling methodology and duration, yielding estimates for the ice Ic–HDL free-energy difference and the HDL–LDL surface tension in harmony with independent simulations and thermodynamic expectations (see Extended Data).
Online Content Any additional Methods, Extended Data display items and Source Data are available in the online version of the paper; references unique to these sections appear only in the online paper. Received 28 November 2013; accepted 24 April 2014.

1.025

0.950

0.875

0.1 Q6

0.2

0.3

Figure 3 | Free-energy surface from unconstrained simulations. The r–Q6 free-energy surface at 228.6 K and 2.4 kbar computed from 16 unconstrained MC simulations initialized in the low-Q6 region. Contours are separated by 1kBT. Because of the separation of timescales between structural relaxation in the liquid phase and ice nucleation, each simulation was run for more than ,100 relaxation times without exhibiting any sign of crystallization. The inset shows autocorrelation functions for density (blue line) and Q6 (red line) computed from the unconstrained MC simulations performed in the LDL region. Fluctuations in density and structural order (Q6) decay in tandem on timescales that are relevant to relaxations within the liquid phase, as demonstrated by both order parameters having mean autocorrelation times of ,106 MC moves.
0.275

0.225

1.

2.

0.175 I (?2)

3. 4.

6.0 0.125 R 6.3 6.6 6.9 0.075 7.2 7.5 HDL 0.1 0.2 0.3 Q6 0.025 0.1 0.2 Q6 0.3 0.4 0.5 0.4 0.5 Ice LDL

5. 6. 7. 8.

9.

10. 11. 12. 13. 14. 15.

Figure 4 | Structural and topological order in the metastable coexisting liquids and in cubic ice. The free-energy surface at 228.6 and 2.4 kbar described by the first moment of the molecular local structure index distribution,  I , and I the crystalline order parameter, Q6. Contours are 1kBT apart. Parameter  successfully distinguishes the three phases based on structural order, characterizing the extent to which molecules in each phase possess a tetrahedral environment with well-separated first and second coordination shells. The inset shows that the three phases have distinctive topological features characterized . by the first moment of the ring size distribution, R

16.

Speedy, R. J. & Angell, C. A. Isothermal compressibility of supercooled water and evidence for a thermodynamic singularity at 245 uC. J. Chem. Phys. 65, 851–858 (1976). Angell, C. A., Oguni, M. & Sichina, W. J. Heat capacity of water at extremes of supercooling and superheating. J. Phys. Chem. 86, 998–1002 (1982). Hare, D. E. & Sorensen, C. M. Densities of supercooled H2O and D2O in 25 mm glass capillaries. J. Chem. Phys. 84, 5085–5089 (1986). Mishima, O. & Stanley, H. E. Decompression-induced melting of ice IV and the liquid–liquid transition in water. Nature 392, 164–168 (1998). Soper, A. K. & Ricci, M. A. Structures of high-density and low-density water. Phys. Rev. Lett. 84, 2881–2884 (2000). Winkel, K., Elsaesser, M. S., Mayer, E. & Loerting, T. Water polyamorphism: reversibility and (dis)continuity. J. Chem. Phys. 128, 044510 (2008). Huang, C. et al. The inhomogeneous structure of water at ambient conditions. Proc. Natl Acad. Sci. USA 106, 15214–15218 (2009). Clark, G. N. I., Hura, G. L., Teixeira, J., Soper, A. K. & Head-Gordon, T. Small-angle scattering and the structure of ambient liquid water. Proc. Natl Acad. Sci. USA 107, 14003–14007 (2010). Poole, P. H., Sciortino, F., Grande, T., Stanley, H. E. & Angell, C. A. Effect of hydrogenbonds on the thermodynamic behavior of liquid water. Phys. Rev. Lett. 73, 1632–1635 (1994). Tanaka, H. Simple view of waterlike anomalies of atomic liquids with directional bonding. Phys. Rev. B 66, 064202 (2002). Holten, V. & Anisimov, M. A. Entropy-driven liquid–liquid separation in supercooled water. Sci. Rep. 2, 713 (2012). Poole, P. H., Sciortino, F., Essmann, U. & Stanley, H. E. Phase behavior of metastable water. Nature 360, 324–328 (1992). Moore, E. B. & Molinero, V. Structural transformation in supercooled water controls the crystallization rate of ice. Nature 479, 506–508 (2011). Speedy, R. J. Stability-limit conjecture. An interpretation of the properties of water. J. Phys. Chem. 86, 982–991 (1982). Sastry, S., Debenedetti, P. G., Sciortino, F. & Stanley, H. E. Singularity-free interpretation of the thermodynamics of supercooled water. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 53, 6144–6154 (1996). Amann-Winkel, K. et al. Water’s second glass transition. Proc. Natl Acad. Sci. USA 110, 17720–17725 (2013).
1 9 J U N E 2 0 1 4 | VO L 5 1 0 | N AT U R E | 3 8 7

?2014 Macmillan Publishers Limited. All rights reserved

RESEARCH LETTER
17. Liu, Y., Palmer, J. C., Panagiotopoulos, A. Z. & Debenedetti, P. G. Liquid–liquid transition in ST2 water. J. Chem. Phys. 137, 214505 (2012). 18. Poole, P. H., Bowles, R. K., Saika-Voivod, I. & Sciortino, F. Free energy surface of ST2 water near the liquid–liquid phase transition. J. Chem. Phys. 138, 034505 (2013). 19. Palmer, J. C., Car, R. & Debenedetti, P. G. The liquid–liquid transition in supercooled ST2 water: a comparison between umbrella sampling and well-tempered metadynamics. Faraday Discuss. 167, 77–94 (2013). 20. Li, Y. P., Li, J. C. & Wang, F. Liquid–liquid transition in supercooled water suggested by microsecond simulations. Proc. Natl Acad. Sci. USA 110, 12209–12212 (2013). 21. Overduin, S. D. & Patey, G. N. An analysis of fluctuations in supercooled TIP4P/ 2005 water. J. Chem. Phys. 138, 184502 (2013). 22. Limmer, D. T. & Chandler, D. The putative liquid–liquid transition is a liquid-solid transition in atomistic models of water. J. Chem. Phys. 135, 134503 (2011). 23. Limmer, D. T. & Chandler, D. The putative liquid–liquid transition is a liquid–solid transition in atomistic models of water. II. J. Chem. Phys. 138, 214504 (2013). 24. Stillinger, F. H. & Rahman, A. Improved simulation of liquid water by molecular dynamics. J. Chem. Phys. 60, 1545–1557 (1974). 25. Lee, J. Y. & Kosterlitz, J. M. New numerical method to study phase transitions. Phys. Rev. Lett. 65, 137–140 (1990). 26. Moore, E. B. & Molinero, V. Ice crystallization in water’s ‘no-man’s land’. J. Chem. Phys. 132, 244504 (2010). 27. Steinhardt, P. J., Nelson, D. R. & Ronchetti, M. Bond-orientational order in liquids and glasses. Phys. Rev. B 28, 784–805 (1983). 28. Shiratani, E. & Sasai, M. Molecular scale precursor of the liquid–liquid phase transition of water. J. Chem. Phys. 108, 3264–3276 (1998). 29. Matsumoto, M., Saito, S. & Ohmine, I. Molecular dynamics simulation of the ice nucleation and growth process leading to water freezing. Nature 416, 409–413 (2002). 30. Frenkel, D. & Smit, B. Understanding Molecular Simulation: From Algorithms to Applications 2nd edn (Academic, 2002). Acknowledgements Computations were performed at the Terascale Infrastructure for Groundbreaking Research in Engineering and Science (TIGRESS) facility at Princeton University. P.G.D. acknowledges support from the National Science Foundation (CHE 1213343), A.Z.P. acknowledges support from the US Department of Energy (DE-SC0002128), and R.C. acknowledges support from the US Department of Energy (DE-SC0008626). Author Contributions J.C.P., R.C., A.Z.P. and P.G.D. planned the study. J.C.P., Y.L. and F.M. performed the simulations and numerical data analysis. J.C.P. and P.G.D. wrote the main paper and methods information. All authors discussed the results and commented on the manuscript at each stage. Author Information Reprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interests. Readers are welcome to comment on the online version of the paper. Correspondence and requests for materials should be addressed to P.G.D. (pdebene@princeton.edu).

3 8 8 | N AT U R E | VO L 5 1 0 | 1 9 J U N E 2 0 1 4

?2014 Macmillan Publishers Limited. All rights reserved

LETTER RESEARCH
METHODS
General sampling protocol. MC simulations in the isothermal-isobaric ensemble employing collective, N-particle smart MC moves31 were used to investigate the lowtemperature phase behaviour of the ST2 water model24, modified for compatibility with the Ewald treatment for long-range electrostatic interactions17,32. The r–Q6 range relevant to each phase under consideration was explored systematically using windowed umbrella sampling30,33. The parameter space was divided into overlapping windows. Independent MC simulations were performed in each window, restricting sampling to the target region with a harmonic restraint: ? ? kr ? ? N ? ?2 kQ ? ? ? ?2 r r {r? z 6 Q6 r N {Q? ?1? W rN ~ 6 2 2 n where r is the vector describing the microscopic coordinates of the N-particle system, kr and kQ6 are spring constants, and parameters r* and Q? 6 specify the window’s centre. Values ranging from 5,000kBT to 10,000kBT (cm6 g22) and from 2,000kBT to 6,000kBT for kr and kQ6 , respectively, proved sufficient to ensure that the simulations sampled in the vicinity of their target window. Technical details regarding the basic MC algorithm, implementation of the ST2 water model, and definition and calculation of Q6 are described in ref. 19. Free-energy analysis. Time series data were collected for r and Q6 in each umbrella window during the post-equilibration, production phase of the MC simulations. The data were subsequently re-sampled with an interval ? equal ? to the maximum statistical inefficiency in each window, g :1z2|max tr ,tQ6 , where tr and tQ6 are the integrated autocorrelation times associated with r and Q6, respectively. The relaxation times for each observable were typically found to be comparable in magnitude (that is, tr <tQ6 ), including within sampling windows in the vicinity of the LDL basin. Two-dimensional r–Q6 histograms were generated from the uncorrelated time series data and subsequently combined using the weighted histogram analysis method of Kumar et al.34 to produce an unbiased estimate of the free energy, F ?r,Q6 ?~{kB T ln?2?r,Q6 ??, where 2 is the microstate probability distribution. Points of liquid–liquid coexistence, where the LDL and HDL basins have equal depths, were located by reweighting in pressure17: F ?r,Q6 ; pzDp,T ?~F ?r,Q6 ; p,T ?zDpN =r ?2? similar duration. Comparison of the results for N 5 192 with the more extensive calculations used to generate Fig. 1 provided verification that the sampling duration and explored range of Q6 were sufficient to accurately reproduce the low-Q6 portion of the free-energy surface. For each system size considered, the height of the barrier separating the liquids was computed at coexistence from the free energy profile along r: ?  ?3? F ?r?~{kB T ln exp?{bF ?r,Q6 ??dQ6 where b 5 (kBT)21. Unconstrained sampling. The r–Q6 free-energy surface in Fig. 3 was computed by performing long, unbiased MC simulations of 192 ST2 water molecules at the estimated point of liquid–liquid coexistence (228.6 K and 2.4 kbar). Equilibrated HDL and LDL configurations extracted from umbrella sampling calculations were used to initialize eight independent simulations in the vicinity of each liquid basin. The unbiased simulations were run for two orders of magnitude longer than the integrated Q6 autocorrelation time in the LDL region. Time series data collected over the duration of each MC simulation were analysed, as described above, to compute free energy. Analysis of local structure index and topological rings. The LSI28 is an order parameter sensitive to heterogeneity in water’s coordination shell capable of distinguishing between molecular configurations characteristic of HDL, LDL and ice. The free-energy surface parameterized by the first moment of the molecular LSI distribution,  I , and Q6 was computed from time series data generated by analysing saved trajectories from the long umbrella sampling simulations used to construct Fig. 1. The uncorrelated data were subsequently re-weighted38,39 to remove the bias, generating the final estimate of F ? I ,Q6 ? shown in Fig. 4. A 0.37-nm cutoff based on the O–O separation distance between neighbouring water molecules was used in the calculation of  I . Additional details regarding the definition and calculation of  I may be found in ref. 28. Ring statistics based on King’s criteria40 were computed from saved trajectories at selected points along the HDL–LDL and LDL–ice paths, applying a 0.35-nm oxygen-based cutoff to determine topological connectivity between adjacent water molecules. Consistency among sampling methods. To verify that our results withstand critical scrutiny, we have studied the dependence of the free-energy surface on sampling duration and methodology, computing the r–Q6 free-energy surface at 228.6 K and 2.4 kbar using several state-of-the-art computational techniques. Extended Data Table 1 lists the methods we have used, along with tr and tQ6 computed in the LDL basin, and the sampling duration in each umbrella window. Four and sixteen identical simulation replicas were used in the well-tempered metadynamics41 and unconstrained MC calculations, respectively, with each replica being run for the reported duration. The free-energy surfaces computed using the different sampling methods are shown in Extended Data Fig. 1. In each case we find two coexisting liquids separated by a ,4kBT free-energy barrier, demonstrating that such results are independent of sampling technique. Extended Data Fig. 1 also demonstrates that the results are devoid of non-equilibrium artefacts. Limmer and Chandler23 have suggested that the LDL basin is such an artefact associated with the sluggish dynamics of ice coarsening, and consequently it was posited that the LDL basin should progressively age as the sampling duration increases, until it eventually vanishes at ,103tQ6 (ref. 23). In contrast, we do not observe significant changes even when the sampling duration is increased by two orders of magnitude from 102 to 104tQ6 . As shown in Extended Data Fig. 1, the techniques that yield such satisfyingly consistent free-energy surfaces include the hybrid MC sampling method42 employed by Limmer and Chandler23. Figure 1 shows that consistent results are obtained even when reversible sampling is performed between the LDL and ice Ic basins. Finally, our results are qualitatively consistent with free-energy calculations employing different variants of the ST2 water model18, microsecond-long MD trajectories exhibiting abrupt and infrequent transitions between HDL and LDL17,43, and previous finite-size scaling studies43. Limmer and Chandler23 have proposed a theory of artificial polyamorphism, which posits that a purported separation of timescales between density and structural relaxations (that is, tQ6 ?tr ) gives rise to an illusory LDL basin associated with the coarsening of ice. To scrutinize this prediction, we examined the density and Q6 autocorrelation functions computed in the LDL region, using the various sampling techniques employed in our study. Extended Data Fig. 2 shows representative autocorrelation functions for three of the sampling techniques. Whereas the density and Q6 autocorrelation functions exhibit transient behaviour at short times where they are separated by more than one order of magnitude, Extended Data Fig. 2 clearly shows that such short-time behaviour is sensitive to the sampling technique and therefore does not provide a physically meaningful description of the coupling between density and structural relaxations in the system. For each sampling technique, we find that density and Q6 fluctuations decay in tandem at long times. It is this techniqueindependent, long-time behaviour that is relevant to sampling the physical properties of the system. Hence, our results demonstrate that Limmer and Chandler’s theory23

where Dp is the pressure shift. Uncertainties in F(r,Q6;p 1 Dp,T) were estimated from the variance computed from 500 resampled r–Q6 free-energy surfaces generated using the Bayesian bootstrap technique described by Hub et al.35. This technique has been shown to provide robust error estimates even in extreme cases where the sampling duration is limited to timescales on the order of the characteristic relaxation time of the biased observable35. Computing the three-phase diagram. Umbrella sampling MC simulations of 192 ST2 water molecules at 228.6 K and 2.2 kbar were used to compute the reversible free-energy surface in Fig. 1. The high-density region (r $ 0.98 g cm23) was sampled by performing independent simulations in 27 density windows in the range 0.98 g cm23 # r* # 1.24 g cm23 using a spacing of 0.01 g cm23 and Q? 6 5 0.05. Simulations in the low-density region (r , 0.98 g cm23) were carried out in four density windows, namely r* 5 0.91, 0.93, 0.95 and 0.97 g cm23. Sampling along Q6 was enhanced at each of the four target densities using Hamiltonian exchange MC moves36, in which attempts were made to swap parameters kQ6 and Q? 6 between replicas in neighbouring Q6 windows. Two independent sets of replicas were used for each value of r* in the low-density region. The first set comprised 16 replicas evenly distributed over the range 0.02 # Q? 6 # 0.17, and 32 replicas were used in the second group to span the interval 0.16 # Q? 6 # 0.625. Exchange attempts were made between even or odd numbered replica pairs with equal probability once every 200 MC moves on average. Simulations were equilibrated for ,104tQ6 in each sampling window. Density, Q6 and the configurational energy were carefully monitored for drift to verify that each simulation had completely equilibrated by the end of this period. Bi-directional sampling between the LDL and ice phase was also performed to serve as an additional check for equilibration in the low-density region. The first generation of simulations was seeded using initial configurations extracted from a trajectory of the LDL freezing into ice Ic, and the second generation was initialized with configurations taken from a melting trajectory. The freezing and melting trajectories were produced by applying a strong umbrella bias to accelerate the phase transition process37. After equilibration, data collection was performed in each window for ,104tQ6 . Histograms generated using data collected from the two generations of simulations in the low-density region were compared to explicitly check for reversibility. The absence of hysteresis confirmed that the simulations were properly equilibrated and sampling the reversible r–Q6 free-energy surface. Finite-size scaling. Umbrella sampling calculations for N 5 192, 300, 400 and 600 ST2 water molecules were performed in the low-Q6 region, using 35 evenly distributed density windows in the range 0.90 g cm23 # r* # 1.24 g cm23. Simulations in each window were equilibrated for at least 102tQ6 , followed by a production phase of

?2014 Macmillan Publishers Limited. All rights reserved

RESEARCH LETTER
can only be justified if the long-time behaviour is completely neglected by defining the relaxation time, for instance, as C(t) 5 e21. Although this definition can in general be used to estimate t, a more careful analysis is required when comparing correlation functions that are decoupled at shorter times but invariably decay together at long times. The physically relevant long-time behaviour may be captured by using a different metric such as the integrated autocorrelation time. Using this definition, we find that tr <tQ6 for each sampling method listed in Extended Data Table 1. Moreover, by re-sampling our data using an interval equal ? ? to the maximum statistical inefficiency in each window, g :1z2|max tr ,tQ6 , we have excluded the possibility that transient, short-time correlations are embedded in the freeenergy surfaces shown in Extended Data Fig. 1 and Figs 1 and 3. We also do not observe significant changes in the free-energy surface shown in Fig. 1 even when the data are re-sampled using an interval of 102g. Consequently, the presence of a LDL basin cannot be attributed to finite-time artefacts associated with transient behaviour occurring on timescales that are orders of magnitude shorter than the sampling interval. Thermodynamic consistency. The free-energy surface in Fig. 1 shows that the coexisting liquids are metastable with respect to ice Ic at 228.6 K and 2.4 kbar, with the ice phase being lower in free energy by ,75kBT (in extensive units for N 5 192) or, equivalently, DGIc–L 5 2742 J mol21. In contrast, Limmer and Chandler’s23 free energy calculations predict ice Ic–liquid coexistence at a nearby state condition (230 K and 2.6 kbar) for the same variant of the ST2 water model (see the middle column in Fig. 13 of ref. 23). To resolve this significant discrepancy, we have used thermodynamic integration (TI) along with an empirical equation of state (EEOS) parameterized to reproduce the experimental properties of water and ice44,45, to estimate DGIc–L under comparable state conditions for water (thus allowing us to subject both our results and those of Limmer and Chandler23 to thermodynamic scrutiny); and to estimate the melting temperature for the ST2 model at 2.6 kbar (thus allowing us to test the very different predictions for the equilibrium melting temperature of ice Ic at 2.6 kbar in the ST2 model). Thermodynamic integration was performed using the identity DGIc{L ?P,T ?:DGIh{L ?P,T ?zDGIc{Ih ?P,T ? ?4? the equation-of-state calculations. In contrast, Extended Data Table 3 shows that Limmer and Chandler’s simulations23, purportedly for the same ST2 variant and at a nearby state condition (2.2 kbar and 230 K), predict that DGIc–L is an order of magnitude smaller than the values calculated by TI using LE and the EEOS. In fact, we find similar disagreement between the TI calculations and the DGIc–L values estimated from the free-energy surfaces reported by Limmer and Chandler22,23, even for the other ST2 variants considered in their studies. Limmer and Chandler23 observed ice Ic–liquid coexistence (that is, DGST2 Ic{L <0) at 230 K and 2.6 kbar for the same variant of the ST2 water model examined in our study (see the middle column of Fig. 13 in ref. 23). Reweighting the free-energy surface shown in Fig. 1 in pressure and using the HDL as a reference, we find DGST2 Ic{L < 2705 J mol21 at 228.6 K at 2.6 kbar. This value for DGST2 Ic{L was used along with LE and the EEOS to predict the melting temperature of ice Ic for the ST2 water model m,ST2 ), providing an estimate of temperature at which our simulations should (TIc be performed to find ice Ic–liquid coexistence at 2.6 kbar. Starting from the initial temperatures (T1) listed in Extended Data Table 2, the LE and EEOS expressions for DGIc–L were integrated at 2.6 kbar to find the temperature, T2, satisfying   ? T m,ST2  ? T2  Ic LDGIc{L LDGST2 Ic{L dT ~ dT ~705 J mol{1 ?6? LT LT T1 228:6K P P We note that T1 and T2 are either defined with respect to ST2’s melting temperature m,ST2 sc for ice Ih at 1 bar (that is, T0 ,Ih ), or the supercooling, DT , as described above. Thus, m,ST2 m,ST2 m,ST2 m < T for calculations performed using T ~ T <T2 z T 2 Ic 0,Ih 0,Ih , whereas TIc  
m,ST2 m,ST2 m m T0 ,Ih {T0,Ih for the latter scenario, where T0,Ih {T0,Ih is the difference between

where subscripts Ic, Ih and L denote ice Ic, ice Ih and the liquid phase, respectively. Two levels of TI were considered for evaluating DGIh–L(P,T): (i) A simple linear extrapolation (LE) using experimental data for the specific m m volume (DV0 ,Ih{L ) and entropy (DS0,Ih{L ) change upon melting for ice Ih at 1 bar:     m m m DGIh{L ?P,T ?~D Vm ?5? 0,Ih{L P{P0,Ih {DS0,Ih{L T {T0,Ih ,
m m where P0 ,Ih 5 1 bar and T0,Ih is the melting temperature at 1 bar. (ii) The empirical equation of state (EEOS) developed in refs 44, 45, which is applicable over the ranges 0–22 kbar and 175–360 K and accurately describes the phase behaviour of liquid water and several ice polymorphs, including ice Ih. The difference in free energy between ices Ic and Ih, DGIc–Ih, was calculated from experimental vapour pressure data for these ice phases46 and the enthalpy difference, DHIc–Ih, measured by calorimetry47–50. The ice phases were assumed to be incompressible, which is justified by the fact that their specific volumes are relatively insensitive to pressure44,45. Because the ST2 water model is over-structured in comparison m,ST2 with real water, it has a melting temperature T0 ,Ih < 300 K for ice Ih at 1 bar (ref. 51), 44 m for water . Two different approaches were which is significantly higher than T0 ,Ih used to account for this behaviour: m,ST2 m was assumed for ice Ih at 1 bar. (i) A melting temperature of T0 ,Ih ~T0,Ih (ii) Thermodynamic integration calculations were performed at the same superm cooling, DT sc :T0 ,Ih {T , with respect to the melting temperature of ice Ih at 1 bar. Our simulations at 228.6 K, for example, were conducted at a supercooling of m,ST2 DT sc 5 71.4 K with respect to T0 ,Ih . In the second approach, the TI was therefore sc m m performed from T0 ,Ih 5 273.15 K to T 5 T0,Ih 2 DT 5 201.75 K to achieve the same supercooling for real water. Extended Data Table 2 shows the values of DGIc–L predicted by LE and the EEOS for water, along with our DGIc–L calculation for the ST2 model at 228.6 K and 2.4 kbar obtained from Fig. 1. Although LE predicts the largest DGIc–L due to the assumption of incompressibility, it provides a reasonable order-of-magnitude estimate for this quantity. The more accurate EEOS, which accounts for the changes in the thermodynamic response functions of the liquid as a function of T and P, predicts that DGIc–L is smaller by a factor of 2 than the estimate obtained using LE. Because the ice phase produced by freezing LDL contains natural imperfections, the predicted DGIc–L for ST2 underestimates the difference in free energy that would be computed using an ideal ice Ic crystal prepared by artificial means. Defects in the crystal may also arise because the number of molecules in our simulations is not commensurate with a cubic surpercell of ice Ic. Despite such defects, however, we find that our DGIc–L value for the ST2 model is in reasonable agreement with the thermodynamic analysis, regardless of the approach used to compute or to assign the reference temperature in

the melting temperature of ice Ih at 1 bar for the ST2 model and real water. m,ST2 at 2.6 kbar obtained using the Extended Data Table 4 lists the estimates of TIc same procedures and reference temperatures as those reported in Extended Data m,ST2 < 260 K, whereas calculations with the more accurTable 2. The LE predicts TIc m,ST2 in the range 272–276 K at 2.6 kbar. To confirm these ate EEOS estimate TIc m,ST2 directly from simulation, using two different techpredictions, we computed TIc m,ST2 niques. In the first approach, TIc was determined using two-phase, ice Ic–liquid (N,PZ,T) MC simulations30, imposing a pressure of 2.6 kbar in the direction perpendicular to the ice Ic–liquid interface. Extended Data Fig. 3 shows the time evolution of the crystalline order parameter, Q6, for simulations performed at different temm,ST2 value predicted by the EEOS. Below 270 K, the simulations peratures near the TIc exhibited a gradual drift towards higher values of Q6, indicating that the system was freezing. Similarly, Q6 decreased for simulations performed above 275 K because of the melting of ice. Our estimate of the melting temperature is therefore the average of m,ST2 < 273 6 3 K at 2.6 kbar, which is in excellent agreement these temperatures, TIc with the range 272–276 K predicted using the EEOS. We also computed the r–Q6 free-energy surface at 275 K and 2.2 kbar for N 5 216 ST2 water molecules using the umbrella sampling procedure described above. Extended Data Fig. 4 shows the resulting r–Q6 free-energy surface after reweighting in pressure using equation (2) to locate the point of ice Ic–liquid coexistence, 275 K and ,2.7 kbar. As Extended Data Table 4 shows, this result is in excellent agreement with our thermodynamic calculations using the EEOS and interfacial simulations. Such values are 30–46 K m,ST2 predicted by Limmer and Chandler at the same pressure23, higher than the TIc demonstrating that those free-energy calculations are inconsistent with reasonable thermodynamic expectations based on accurate equations of state for real water and the established physical properties of the ST2 water model. We have shown that the free-energy surface shown in Fig. 1 is consistent with expectations based on thermodynamic arguments. This is demonstrated by the fact that our estimate of DGIc–L for the ST2 model at 228.6 K and 2.4 kbar is in good agreement with calculations performed using the accurate EEOS for water. In addition, we have also demonstrated thermodynamic consistency by using the EEOS along m,ST2 < 272–276 K. This with our DGST2 Ic{L value at 228.6 and 2.6 kbar to predict TIc prediction was verified by performing simulations of the ice Ic–liquid interface and m,ST2 at umbrella sampling calculations. Such results demonstrate conclusively that TIc ,2.6 kbar is ,40–45 K higher than reported by Limmer and Chandler23. It therefore seems that their free-energy surface (middle column of Fig. 13 in ref. 23) is distorted to such an extent that the output of their simulations corresponds to an effectively higher temperature. To observe ice Ic–liquid coexistence at 2.6 kbar, as reported by Limmer and Chandler23, this effective temperature would have to be well above the estimated liquid–liquid critical temperature (Tc < 237 K for our model32) for any reasonable variant of the ST2 water model, explaining the absence of a LDL basin in their free-energy surfaces23. Because the two liquids are only separated by a ,4kBT barrier at 228.6 K and 2.4 kbar, the free-energy surface must be accurately computed to observe the LDL basin. At odds with this requirement, we find a ,70kBT discrepancy between our respective estimates for DGIc–L near 228.6 K and 2.6 kbar, which cannot simply be dismissed as non-equilibrium artefacts, as suggested by Limmer and Chandler22,23. Although the precise numerical origin of this discrepancy

?2014 Macmillan Publishers Limited. All rights reserved

LETTER RESEARCH
is still under investigation, we showed above (see the section on Consistency among sampling methods, and Extended Data Fig. 1) that liquid–liquid coexistence is observed when we perform free-energy calculations using the hybrid MC technique42 employed by Limmer and Chandler23. In our hybrid MC implementation we use the molecular dynamics integrator of Miller et al.52, whereas Limmer and Chandler23 employed the constraint algorithm SETTLE53 to simulate rigid ST2 water molecules. Although we have not yet implemented this integrator, Reinhardt et al.37 recently observed ‘catastrophic’ divergence from the well-established equation of state for the TIP4P/2005 water model when hybrid MC simulations were performed with SETTLE. A more comprehensive discussion of the different perspectives regarding the liquid–liquid phase transition in ST2 water, computational approaches and related studies has recently been published19,54. As a final check, we followed the procedure described by Hunter and Reinhardt55 to estimate the liquid–liquid surface tension, cL–L, from our finite-size scaling data. We find that cL–L < 2 mJ m22, which is comparable to vapour–liquid surface tensions for various water models56 at similar reduced temperatures near the vapour–liquid critical point (that is, cV–L < 5.6–1.5 mJ m22 for T/Tc < 0.95–0.98), and an order of magnitude smaller than the cIh–L < 23 mJ m22 reported by Handel et al.57 for the ice Ih–liquid surface tension in TIP4P. Thus, the small value of cL–L is thermodynamically consistent with our observation that two liquids are forming an interface, not a liquid and a coarsening crystal, and with the fact that our simulations are performed relatively close to the estimated liquid–liquid critical point at a reduced temperature of T/Tc < 0.96.
31. Rossky, P. J., Doll, J. D. & Friedman, H. L. Brownian dynamics as smart Monte Carlo simulation. J. Chem. Phys. 69, 4628–4633 (1978). 32. Liu, Y., Panagiotopoulos, A. Z. & Debenedetti, P. G. Low-temperature fluid-phase behavior of ST2 water. J. Chem. Phys. 131, 104508 (2009). 33. Torrie, G. M. & Valleau, J. P. Nonphysical sampling distributions in Monte Carlo free-energy estimation: umbrella sampling. J. Comput. Phys. 23, 187–199 (1977). 34. Kumar, S., Bouzida, D., Swendsen, R. H., Kollman, P. A. & Rosenberg, J. M. The weighted histogram analysis method for free-energy calculations on biomolecules. 1. The method. J. Comput. Chem. 13, 1011–1021 (1992). 35. Hub, J. S., de Groot, B. L. & van der Spoel, D. g_wham: a free weighted histogram analysis implementation including robust error and autocorrelation estimates. J. Chem. Theory Comput. 6, 3713–3720 (2010). 36. Sugita, Y. & Okamoto, Y. Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 314, 141–151 (1999). 37. Reinhardt, A., Doye, J. P. K., Noya, E. G. & Vega, C. Local order parameters for use in driving homogeneous ice nucleation with all-atom models of water. J. Chem. Phys. 137, 194504 (2012). 38. Bonomi, M., Barducci, A. & Parrinello, M. Reconstructing the equilibrium boltzmann distribution from well-tempered metadynamics. J. Comput. Chem. 30, 1615–1621 (2009). 39. Gee, J. & Shell, M. S. Two-dimensional replica exchange approach for peptidepeptide interactions. J. Chem. Phys. 134, 064112 (2011). 40. King, S. V. Ring configurations in a random network model of vitreous silica. Nature 213, 1112–1113 (1967). 41. Barducci, A., Bussi, G. & Parrinello, M. Well-tempered metadynamics: a smoothly converging and tunable free-energy method. Phys. Rev. Lett. 100, 020603 (2008). 42. Duane, S., Kennedy, A. D., Pendleton, B. J. & Roweth, D. Hybrid Monte Carlo. Phys. Lett. B 195, 216–222 (1987). 43. Kesselring, T. A. et al. Finite-size scaling investigation of the liquid–liquid critical point in ST2 water and its stability with respect to crystallization. J. Chem. Phys. 138, 244506 (2013). 44. Choukroun, M. & Grasset, O. Thermodynamic model for water and high-pressure ices up to 2.2 GPa and down to the metastable domain. J. Chem. Phys. 127, 124506 (2007). 45. Choukroun, M. & Grasset, O. Thermodynamic data and modeling of the water and ammonia-water phase diagrams up to 2.2 GPa for planetary geophysics. J. Chem. Phys. 133, 144502 (2010). 46. Shilling, J. E. et al. Measurements of the vapor pressure of cubic ice and their implications for atmospheric ice clouds. Geophys. Res. Lett. 33, L17801 (2006). 47. Mayer, E. & Hallbrucker, A. Cubic ice from liquid water. Nature 325, 601–602 (1987). 48. Yamamuro, O., Oguni, M., Matsuo, T. & Suga, H. Heat capacity and glass transition of pure and doped cubic ices. J. Phys. Chem. Solids 48, 935–942 (1987). 49. Handa, Y. P., Klug, D. D. & Whalley, E. Difference in energy between cubic and hexagonal ice. J. Chem. Phys. 84, 7009–7010 (1986). 50. Mcmillan, J. A. & Los, S. C. Vitreous ice: irreversible transformations during warm-up. Nature 206, 806–807 (1965). 51. Weber, T. A. & Stillinger, F. H. Pressure melting of ice. J. Chem. Phys. 80, 438–443 (1984). 52. Miller, T. F. et al. Symplectic quaternion scheme for biophysical molecular dynamics. J. Chem. Phys. 116, 8649–8659 (2002). 53. Miyamoto, S. & Kollman, P. A. SETTLE: an analytical version of the SHAKE and RATTLE algorithm for rigid water models. J. Comput. Chem. 13, 952–962 (1992). 54. Palmer, J. C. General discussion. Faraday Discuss. 167, 118–127 (2013). 55. Hunter, J. E. & Reinhardt, W. P. Finite-size scaling behavior of the free energy barrier between coexisting phases: determination of the critical temperature and interfacial tension of the Lennard–Jones fluid. J. Chem. Phys. 103, 8627–8637 (1995). 56. Vega, C. & de Miguel, E. Surface tension of the most popular models of water by using the test-area simulation method. J. Chem. Phys. 126, 154707 (2007). 57. Handel, R., Davidchack, R. L., Anwar, J. & Brukhno, A. Direct calculation of solid– liquid interfacial free energy for molecular systems: TIP4P ice–water interface. Phys. Rev. Lett. 100, 036104 (2008).

?2014 Macmillan Publishers Limited. All rights reserved

RESEARCH LETTER

Extended Data Figure 1 | Reversible free-energy surfaces at 228.6 K and 2.4 kbar computed using different sampling techniques. Surfaces on the top row were computed using (from left to right) umbrella sampling MC, welltempered metadynamics and unconstrained MC; the bottom row shows results

from hybrid MC, parallel tempering MC and Hamiltonian exchange MC simulations. The free-energy barrier separating the liquid basins is ,4kBT for all of the surfaces shown. Contours are 1kBT apart and uncertainties are estimated to be less than 0.5kBT.

?2014 Macmillan Publishers Limited. All rights reserved

LETTER RESEARCH

Extended Data Figure 2 | Autocorrelation functions for different sampling techniques. Autocorrelation functions for density (blue) and Q6 (red) computed in the LDL region using unconstrained MC (left), hybrid MC (centre) and Hamiltonian exchange MC (right). The correlation functions were

calculated by averaging results from at least 12 independent simulations. Density and Q6 fluctuations decay on very similar timescales, despite exhibiting technique-dependent transient behaviour where these processes may be separated by more than one order of magnitude.

?2014 Macmillan Publishers Limited. All rights reserved

RESEARCH LETTER

Extended Data Figure 3 | Time evolution of the crystalline order parameter in two-phase MC simulations of the ice Ic–liquid interface at 2.6 kbar. The MC simulations were initiated from configurations containing 512 and 670 ST2 water molecules in the ice Ic and liquid phases, respectively. The x and y dimensions of the simulation cells were fixed in accord with the lattice constant for ice Ic, which was determined at each temperature by performing a separate

calculation for the bulk ice phase, while the z dimension was allowed to fluctuate so as to impose a constant pressure of 2.6 kbar perpendicular to the ice–liquid interface. Drift of Q6 towards higher or lower values indicates that the system is freezing or melting. The melting temperature of 273 6 3 K at 2.6 kbar was estimated by averaging the lowest and highest temperatures, respectively, at which melting and freezing were observed.

?2014 Macmillan Publishers Limited. All rights reserved

LETTER RESEARCH

Extended Data Figure 4 | Reversible free-energy surface at 275 K and 2.7 kbar demonstrating ice Ic–liquid coexistence. The liquid and ice Ic basins have equal depths with respect to the saddle point, indicating that the reported state condition is a point of coexistence. Such results confirm the estimates of

the melting temperature for ice Ic at 2.6 kbar obtained from TI calculations using the EEOS and the two-phase MC simulations of the ice–liquid interface. Contours are 1kBT apart.

?2014 Macmillan Publishers Limited. All rights reserved

RESEARCH LETTER
Extended Data Table 1 | Sampling methods

* Collective, smart MC31 moves used. { Relaxation times estimated from unbiased simulations using the same types of MC moves. { Rigid body integrator of Miller et al.52; ,10 molecular dynamics integration steps per MC move. 1 Eight replicas spaced between 228.6 and 272 K. | | Bi-directional sampling performed between the LDL and crystal to ensure reversibility. State-of-the-art sampling methods used to perform free-energy analysis, along with integrated autocorrelation times for density and the crystalline order parameter Q6 (tr and tQ6 , respectively) computed in the LDL basin at 228.6 K and 2.4 kbar, and the sampling duration in each umbrella sampling window given in terms of tQ6 .

?2014 Macmillan Publishers Limited. All rights reserved

LETTER RESEARCH
Extended Data Table 2 | Comparison of ice Ic–liquid free-energy differences obtained from thermodynamic integration and from results presented in the text for the ST2 model

Ice Ic–liquid free-energy differences (DGIc–L) predicted by LE and the EEOS for water are in good agreement with the DGIc–L value calculated for the ST2 model at 228.6 K and 2.4 kbar from the data presented in Fig. 1. The TI calculations using LE and the EEOS were performed using two different reference temperatures (described in Methods) to account for ST2’s over-structured nature in comparison with real water.

?2014 Macmillan Publishers Limited. All rights reserved

RESEARCH LETTER
Extended Data Table 3 | Comparison of ice Ic–liquid free-energy differences obtained from thermodynamic integration and from results presented by Limmer and Chandler23 for the ST2 water model

* Estimated from Fig. 5(b) of ref. 23. Ice Ic–liquid free-energy differences (DGIc–L) predicted by LE and the EEOS for water are in poor agreement with the DGIc–L value obtained by Limmer and Chandler23 for the ST2 model at 230 K and 2.2 kbar. Such disagreement demonstrates that Limmer and Chandler’s results do not withstand thermodynamic scrutiny and fail to provide a reasonable description of ST2’s phase behaviour. The TI calculations using LE and the EEOS were performed using two different reference temperatures (described in Methods) to account for ST2’s over-structured nature in comparison with real water.

?2014 Macmillan Publishers Limited. All rights reserved

LETTER RESEARCH
Extended Data Table 4 | Estimates of the melting temperature for ice Ic at 2.6 kbar for the ST2 water model

* Coexistence pressure is 2.7 kbar. { Estimated from Fig. 13 (middle column) of ref. 23. Comparison of melting temperature estimates for ice Ic at 2.6 kbar for the ST2 water model calculated using the TI schemes and empirical equations of state for water described in Methods. The estimates of m,ST2 obtained from TI using the accurate EEOS of Choukroun and Grasset44,45 are in excellent TIc agreement with values computed directly from two-phase MC simulations of the ice Ic–liquid interface m,ST2 and umbrella sampling MC simulations. In contrast, the TIc at 2.6 kbar estimated from Limmer and Chandler’s23 umbrella sampling simulations with the ST2 water model is lower by more than 40 K, demonstrating severe thermodynamic inconsistencies with their free-energy calculations.

?2014 Macmillan Publishers Limited. All rights reserved


相关文章:
更多相关标签: