Cluster Monte Carlo and numerical mean field analysis for the water liquid–liquid phase transition Marco G. Mazza,1 Kevin Stokely,1 Elena Strekalova,1 H. Eugene Stanley,1 and Giancarlo Franzese2 1Center for Polymer Studies and Department of Physics, Boston University, Boston, Massachusetts 02215, USA 2Departament de Fisica Fonamental, Universitat de Barcelona, Diagonal 647, 08028 Barcelona, Spain By the Wolff’s cluster Monte Carlo simulations and numerical minimization within a mean field approach,westudythelow temperaturephasediagram ofwater, adoptingacell modelthat repro- ducestheknownpropertiesofwaterinitsfluidphases. Both methodsallows ustostudythewater 9 thermodynamic behavior at temperatures where other numerical approaches –both Monte Carlo 0 andmoleculardynamics–areseriouslyhamperedbythelargeincreaseofthecorrelationtimes. The 0 clusteralgorithm also allows ustoemphasize thattheliquid–liquidphasetransition correspondsto 2 thepercolation transition of tetrahedrally ordered water molecules. n a PACSnumbers: 61.20.Ja,61.20.Gy J 0 2 INTRODUCTION [10, 12, 13], (ii) with Metropolis MC simulations [8, 14] and (iii) with Wang-Landau MC density of state algo- ] rithm [15]. Recent Metropolis MC simulations [8] have t Water is possibly the most important liquid for life f shown that very large times are needed to equilibrate o [1] and, at the same time, is a very peculiar liquid [2]. the system as T → 0, as a consequence of the onset of s In the stable liquid regime its thermodynamic response . the glassy dynamics. The implementation of the Wolff’s t functions behave qualitatively differently than a typi- a clusters MC dynamics, presented here, allows us to (i) m cal liquid. The isothermal compressibility K , for ex- T drasticallyreducetheequilibrationtimesofthemodelat ample, has a minimum as a function of temperature at - verylowT and(ii)giveageometricalcharacterizationof d T = 46 ◦C, while for a typical liquid K monotonically n T theregionsofcorrelatedwatermolecules(clusters)atlow decreases upon cooling. Water’s anomalies become even o T and show that the liquid–liquid phase transition can morepronouncedasthesystemiscooledbelowthemelt- c be interpreted as a percolation transition of the tetrahe- [ ing point and enters the metastable supercooled regime drally ordered clusters. [3]. 2 v Different hypothesis have been proposed to rational- 8 ize the anomalies of water [4]. All these interpretations, THE MODEL 8 but one, predict the existence of a liquid–liquid phase 6 transition in the supercooled state, consistent with the ThesystemconsistsofN particlesdistributedwithina 4 . experiments to date [4] and supported by different mod- volumeV inddimensions. ThevolumeisdividedintoN 0 els [2]. cells of volume v with i∈[1,N]. For sake of simplicity, 1 i 8 To discriminate among the different interpretations, thesecellsarechosenofthesamesize,vi =V/N,butthe 0 many experiments have been performed [5]. However, generalizationtothecaseinwhichthevolumecanchange v: the freezing in the temperature-range of interest can be without changes in the topology of the nearest–neighbor i avoided only for water in confined geometries or on the (n.n.) is straightforward. By definition, vi ≥ v0, where X surface of macromolecules [4, 6]. Since experiments in v0 is the molecule hard-core volume. Each cell has a ar the supercooled region are difficult to perform, numeri- variable ni = 0 for a gas-like or ni = 1 for a liquid-like cal simulations have played an important role in recent cell. We partition the total volume in a way such that years to help interpret the data. However, also the sim- each cell has at least four n.n. cells, e.g. as in a cubic ulations atverylow temperatureT arehamperedby the lattice in 3d or a squarelattice in 2d. Periodicboundary glassy dynamics of the empirical models of water [7, 8]. conditions are used to limit finite–size effects. For these reasons is important to implement more effi- The system is described by the Hamiltonian [10] cient numerical simulations for simple models, able to capture the fundamental physics of water but also less H =−ǫ ninj −J ninjδσij,σji+ computationally expensive. Here we introduce the im- hXi,ji Xhiji plementation of a Wolff’s cluster algorithm [9] for the −J n δ , (1) Monte Carlo (MC) simulations of a cell model for wa- σ i σik,σil ter [10]. The model is able to reproduce all the differ- Xi (Xk,l)i ent scenarios proposed to interpret the behavior of wa- where ǫ>0 is the strength of the van der Waals attrac- ter [11] and has been analyzed (i) with mean field (MF) tion, J >0 accounts for the hydrogenbond energy, with 2 constant specific volume increase due to the hydrogen bond formation. MEAN–FIELD ANALYSIS In the mean–field (MF) analysis the macrostate of the system in equilibrium at constant pressure P and temperature T (NPT ensemble) may be determined by a minimization of the Gibbs free energy per molecule, g ≡(hH i+PV −TS)/N , where w N = n (4) w i i X is the total number of liquid-like cells, and S =S +S n σ is the sum of the entropy S over the variables n and n i the entropy S over the variables σ . σ ij A MF approach consists of writing g explicitly using the approximations FIG. 1: A pictorial representation of five water molecules in n n −→2Nn2 (5) 3d. Two hydrogen bonds (grey links) connect the hydrogens i j (inblue)ofthecentralmoleculewiththeloneelectrons(small <Xij> gray lines) of two nearest neighbor (n.n.) molecules. A bond n n δ −→2Nn2p (6) index (arm) with q = 6 possible values is associated to each i j σij,σji σ <ij> hydrogen and lone electron, giving rise to q4 possible orien- X tational states for each molecule. A hydrogen bond can be ni δσik,σil −→6Nnpσ (7) formed only if the two facing arms of the n.n. molecules are Xi (Xk,l)i inthesamestate. Armsonthesamemoleculeinteractamong themselves to mimic the O-O-O interaction that drives the where n = Nw/N is the average of ni, and pσ is the molecules toward a tetrahedral local structure. probabilitythat two adjacentbond indices σ are in the ij appropriate state to form a hydrogen bond. Therefore, in this approximation we can write four(Potts)variablesσ =1,...,qrepresentingbondin- ij dicesofmoleculeiwithrespecttothefourn.n. molecules V =Nv +2Nn2p v , (8) 0 σ HB j, δ = 1 if a =b and δ =0 otherwise, and hi,ji de- a,b a,b hHi=−2[ǫn+(Jn+3J )p ]nN. (9) notesthatiandj aren.n. Themodeldoesnotassumea σ σ privilegedstateforbondformation. Anytimetwofacing The probabilityp , properlydefined as the thermody- σ bondindices(arms)areinthesame(Potts)state,abond namic average over the whole system, is approximated is formed. The third term represents an intramolecular astheaverageovertwoneighboringmolecules,underthe (IM) interaction accounting for the O–O–O correlation effect of the mean-field h of the surrounding molecules [16], locally driving the molecules toward a tetrahedral configuration. When the bond indices of a molecule are pσ = δσij,σji h. (10) in the same state, the energy is decreased by an amount The ground state of the(cid:10)system(cid:11)consists of all N vari- J >0 and we associate this local ordered configuration σ ables n = 1, and all σ in the same state. At low to a local tetrahedral arrangement [17]. The notation i ij temperatures, the symmetry will remain broken, with (k,l) indicates one of the six different pairs of the four i the majority of the σ in the preferred state. We as- bond indices of molecule i (Fig.1). ij sociate this preferred state to the tetrahedral order of Experiments show that the formation of a hydrogen the molecules and define m as the density of the bond bond leads to a local volume expansion [2]. Thus in our σ indicesinthetetrahedralstate,with0≤m ≤1. There- system the total volume is σ fore, the number density n of bond indices σ is in the σ ij V =Nv +N v , (2) tetrahedral state is 0 HB HB 1+(q−1)m where σ n = . (11) σ q N ≡ n n δ (3) HB i j σij,σji Since an appropriate form for h is [10] <i,j> X is the total number of hydrogen bonds, and v is the h=3J n , (12) HB σ σ 3 weTohbetaMinFtehxaptre3sJqsσio≤nshfo≤rt3hJeσe.ntropiesSn oftheN vari- [ε]y g-1.9 kBT/ε = 0.06 Pv0/ε = 0.7 ables n , and S of the 4Nn variables σ , are then [12] g i σ ij r e n S =−k N(nlog(n)+(1−n)log(1−n)) (13) E n B e e r F s Sσ =−kB4Nn[nσlog(nσ)+ b b (1−nσ)log(1−nσ)+log(q−1)], (14) Gi ar -1.95 k T/ε = 0.08 where kB is the Boltzmann constant. ol B M Equating 0 0.2 0.4 0.6 0.8 1 Tetrahedral Order Paramter m σ pσ ≡n2σ+ (1q−−n1σ)2, (15) [ε]g k T/ε = 0.05 Pv /ε = 0.8 y B 0 g with the approximate expression in Eq. (10), allows for r-1.7 e solutionofn ,andhenceg,intermsoftheorderparam- n σ E eter mσ and n. e e By minimizing numerically the MF expression of g r F withrespectton andmσ, wefindtheequilibriumvalues bs n(eq) and m(eq) and, with Eqs. (4) and (2), we calcu- b σ i late the density ρ at any (T,P) and the full equation G-1.75 of state. An example of minimization of g is presented ar k T/ε = 0.07 ol B in Fig. 2 where, for the model’s parameters J/ǫ = 0.5, M 0 0.2 0.4 0.6 0.8 1 Jσ/ǫ = 0.05, vHB/v0 = 0.5, q = 6, a discontinuity Tetrahedral Order Paramter m σ in mσ(eq) is observed for Pv0/ǫ > 0.8. As discussed in ε] Ref.s [10, 14] this discontinuity correspondsto a first or- [g k T/ε = 0.04 Pv /ε = 0.9 der phase transition between two liquid phases with dif- y B 0 g ferent degree of tetrahedralorderand, as a consequence, r e n-1.5 different density. The higher P at which the change E in m(eq) is continuous, corresponds to the pressure of a e σ e liquid–liquidcriticalpoint(LLCP).Theoccurrenceofthe Fr LLCP is consistent with one of the possible interpreta- s b tions of the anomaliesof water, as discussed in Ref. [12]. ib-1.55 G Hreopwroedveurc,esfoarlsdoifftehreenotthcehropicreospoosfepdasrcaemneatreiorss,[1th1e].model lar kBT/ε = 0.06 o M 0 0.2 0.4 0.6 0.8 1 Tetrahedral Order Paramter m THE SIMULATION WITH THE WOLFF’S σ CLUSTERS MONTE CARLO ALGORITHM FIG. 2: Numerical minimization of the molar Gibbs free en- ergy g in the mean field approach. The model’s parameters To performMC simulationsin the NPT ensemble, we are J/ǫ=0.5, Jσ/ǫ=0.05, vHB/v0 =0.5 and q=6. In each consider a modified version of the model in which we panel we present g (dashed lines) calculated at constant P allow for continuous volume fluctuations. To this goal, and different values of T. The thick line crossing the dashed (i) we assume that the system is homogeneous with all lines connects the minima m(σeq) of g at different T. Upper the variables ni set to 1 and all the cells with volume panel: Pv0/ǫ=0.7, for T going from kBT/ǫ=0.06 (top) to v = V/N; (ii) we consider that V ≡ V +N v , kBT/ǫ = 0.08 (bottom). Middle panel: Pv0/ǫ = 0.8, for T MC HB HB where V > Nv is a dynamical variable allowed to going from kBT/ǫ = 0.05 (top) to kBT/ǫ = 0.07 (bottom). MC 0 Lower panel: Pv0/ǫ = 0.9, for T going from kBT/ǫ = 0.04 fluctuateinthesimulations;(iii)wereplacethefirst(van (top) to kBT/ǫ = 0.06 (bottom). In each panel dashed lines der Waals) term of the Hamiltonian in Eq. (1) with a are separated by kBδT/ǫ = 0.001. In all the panels m(σeq) Lennard-Jones potential with attractive energy ǫ > J increases when T decreases, being 0 (marking theabsence of and truncated at the hard-core distance tetrahedralorder)atthehighertemperaturesand≃0.9(high tetrahedralorder)atthelowesttemperature. BychangingT, U (r)≡ ∞ if r6r0, (16) m(σeq) changes in a continuous way for Pv0/ǫ = 0.7 and 0.8, W (ǫ rr0 12− rr0 6 if r>r0. but discontinuousfor Pv0/ǫ=0.9 and higher P. h(cid:0) (cid:1) (cid:0) (cid:1) i 4 where r ≡ (v )1/d; the distance between two n.n. the same random constant φ∈[1,...q] 0 0 molecules is (V/N)1/d, and the distance r between two generic molecules is the Cartesian distance between the σinew = σiold+φ mod q. (18) center of the cells in which they are included. In order to impleme(cid:0)nt a con(cid:1)stant P ensemble we let The simplification (i) could be removed, by allowing the volume fluctuate. A small increment ∆r/r = 0.01 thecellstoassumedifferentvolumesv andkeepingfixed 0 i is chosen with uniform random probability and added the number of possible n.n. cells. However, the results to the current radius of a cell. The change in volume of the model under the simplification (i) compares well ∆V ≡ Vnew − Vold and van der Waals energy ∆E with experiments [12]. Furthermore, the simplification W is computed and the move is accepted with probability (i) allows to drastically reduce the computational cost min(1,exp[−β(∆E +P∆V −T∆S)]), where ∆S ≡ of the evaluation of the U (r) term from N(N −1) to W W −Nk ln(Vnew/Vold) is the entropic contribution. N−1operations. Thechanges(i)–(iii)modifythemodel B used for the mean field analysis and allow off-lattice MC simulations for a cell model in which the topology of the MONTE CARLO CORRELATION TIMES molecules (i.e. the number of n.n.) is preserved. The comparison of the mean field results with the MC simu- The cluster MC algorithm described in the previous lationsshowthatthesechangesdonotmodifythephysics section turns out to be very efficient at low T, allow- of the system. ing to study the thermodynamics of deeply supercooled We perform MC simulations with N = 2500 and water with quite intriguing results [21]. To estimate the N = 10000 molecules, each with four n.n. molecules, efficiencyoftheclusterMCdynamicswithrespecttothe at constant P and T, in 2d, and with the same parame- standard Metropolis MC dynamics, we evaluate in both ters used for the mean field analysis. To each molecules dynamics, and compare, the autocorrelation function of we associate a cell on a square lattice. The Wolff’s algo- theaveragemagnetizationpersiteM ≡ 1 σ ,where rithm is based on the definition of a cluster of variables i 4 j ij the sum is over the four bonding arms of molecule i. choseninsuchawaytobethermodynamicallycorrelated P [18,19]. TodefinetheWolff’scluster,abondindex(arm) of a molecule is randomly selected; this is the initial el- 1 hM (t +t)M (t )i−hM i2 i 0 i 0 i C (t)≡ . (19) ement of a stack. The cluster is grown by first checking M N hM2i−hM i2 the remaining arms of the same initial molecule: if they Xi i i are in the same Potts state, then they are added to the For sake of simplicity, we define the MC dynamics au- stack with probability psame ≡ min[1,1−exp(−βJσ)] tocorrelationtime τ as the time, measured in MC steps, [9], where β ≡ (kBT)−1. This choice for the proba- whenCM(τ)=1/e. HerewedefineaMCstepas4N up- bility psame depends on the interaction Jσ between two dates of the bond indices followed by a volume update, arms on the same molecule and guarantees that the i.e. as 4N +1 steps of the algorithm. connected arms are thermodynamically correlated [19]. In Fig. 3 we show a comparison of C (t) for the M Next, the arm of a new molecule, facing the initially Metropolis and Wolff algorithm implementations of this chosen arm, is considered. To guarantee that connected model for a system with N =50×50, at three tempera- facingarmscorrespondtothermodynamicallycorrelated tures along an isobar below the LLCP, and approaching variables,isnecessary[18]tolinkthemwiththeprobabil- the line of the maximum, but finite, correlation length, itypfacing ≡min[1,1−exp(−βJ′)]whereJ′ ≡J−PvHB alsoknownas Widomline TW(P) [12]. In the toppanel, istheP–dependenteffectivecouplingbetweentwofacing at T ≫ T (P) (k T/ǫ = 0.11, Pv /ǫ = 0.6), we find W B 0 arms as results from the enthalpy H +PV of the sys- a correlation time for the Wolff’s cluster MC dynamics tem. It is important to note that J′ can be positive or τ ≈3×103,andfortheMetropolisdynamicsτ ≈106. W M negative depending on P. If J′ > 0 and the two facing In the middle panel, at T > T (P) (k T/ǫ = 0.09, W B arms are in the same state, then the new arm is added Pv /ǫ = 0.6) the difference between the two correlation 0 to the stack with probability pfacing; if J′ < 0 and the times is larger: τW ≈2.5×103, τM ≈3×106. The bot- twofacingarmsareindifferentstates,thenthe new arm tom panel, at T ≃ T (P) (k T/ǫ = 0.06, Pv /ǫ = 0.6) W B 0 is added with probability pfacing [20]. Only after every shows τW ≈3.7×102, while τM is beyond the accessible possible directionofgrowthfor the clusterhas beencon- time window (τ >107). M sideredthevaluesofthearmsarechangedinastochastic SinceasT →0thesystementersaglassystate[8],the way; again we need to consider two cases: (i) if J′ > 0, efficiency τ /τ grows at lower T allowing the evalua- M W all arms are set to the same new value tion of thermodynamics averages even at T ≪ T [21]. C In particular, the cluster MC algorithm turns out to be σnew = σold+φ mod q (17) very efficient when approaching the Widom line in the where φis a randomin(cid:0)teger bet(cid:1)ween1 andq; (ii) ifJ′ < vicinity of the LLCP, with an efficiency of the order of 0, the state of every single arm is changed (rotated) by 104. We plan to analyze in a systematic way how the 5 1 the same way as the average size of the regions of cor- ) (t related molecules [19], i.e. a Wolff’s cluster statistically M C 0.8 represents a region of correlated molecules. Moreover, on k T/ε=0.11 the mean cluster size diverges at the critical point with cti 0.6 B the same exponent of the Potts magnetic susceptibility n u [19], and the clusters percolate at the critical point, as f Metropolis on 0.4 Wolff we will discuss in the next section. ati el rr 0.2 o PERCOLATING CLUSTERS OF CORRELATED C MOLECULES 0 100 101 102 103 104 105 106 107 Time t [MC steps] TheefficiencyoftheWolff’sclusteralgorithmisacon- 1 ) sequence of the exact relation between the average size (t M of the finite clusters and the average size of the regions C 0.8 on k T/ε=0.09 of thermodynamically correlated molecules. The proof cti 0.6 B ofthisrelationatanyT derivesstraightforwardfromthe n proof for the case of Potts variables [19]. This relation u n f Metropolis allows to identify the clusters built during the MC dy- atio 0.4 Wolff namicswiththecorrelatedregionsandemphasizes(i)the el appearance of heterogeneities in the structural correla- orr 0.2 tions[22],and(ii)theonsetofpercolationoftheclusters C of tetrahedrally ordered molecules at the LLCP [23], as 0100 101 102 103 104 105 106 107 shown in Fig. 4. Time t [MC steps] A systematic percolation analysis [18] is beyond the 1 goal of this report, however configurations such as those ) (t in Fig. 4 allow the following qualitative considerations. M C 0.8 At T >TC the averagecluster size is much smaller than on k T/ε=0.06 thesystemsize. Hence,thestructuralcorrelationsamong cti 0.6 B the molecules extends only to short distances. This sug- n u gests that the correlation time of a local dynamics, such f Metropolis on 0.4 Wolff asMetropolisMCormoleculardynamics,wouldbeshort ati on average at this temperature and pressure. Neverthe- rrel 0.2 less,thesystemappearsstronglyheterogeneouswiththe Co coexistence of large and small clusters, suggesting that the distribution of correlation times evaluated among 0 100 101 102 103 104 105 106 107 moleculesatagivendistancecouldbestronglyheteroge- Time t [MC steps] neous. The clusters appear mostly compact but with a fractal surface, suggesting that borders between clusters FIG.3: ComparisonoftheautocorrelationfunctionCM(t)for can rapidly change. theMetropolis(circles)andWolff(squares)implementationof the present model. We show the temperatures kBT/ǫ=0.11 AtT ≃TC thereisonelargecluster,inredontheright (top panel), kBT/ǫ = 0.09 (middle panel), kBT/ǫ = 0.06 of the middle panel of Fig. 4, with a linear size compa- (bottom panel), along the isobar Pv0/ǫ = 0.6 close to the rable to the system linear extension and spanning in the LLCP for N =50×50. vertical direction. The appearance of spanning clusters showstheonsetofthepercolationgeometricaltransition. At this state point the correlation time of local, such as efficiency τM/τW grows on approaching the LLCP. This Metropolis MC dynamics or molecular dynamics would result is well known for the standard liquid-gas critical be very slow as a consequence of the large extension of point [9] and, on the basis of our results, could be ex- the structurally correlated region. On the other hand, tended also to the LLCP. However, this analysis is very the correlation time of the Wolff’s cluster dynamics is expensiveintermsofCPUtimeandgoesbeyondthegoal short because it changes in one single MC step the state ofthe presentwork. Nevertheless,the percolationanaly- of all the molecules in clusters, some of them with very sis,presentedinthe nextsection,helpsinunderstanding largesize. Oncethe spanningclusterisformed,itbreaks the physical reason for this large efficiency. the symmetry of the system and a strong effective field The efficiency is a consequenceofthe factthatthe av- actsonthemoleculesnearitsbordertoinducetheirreori- erage size of Wolff’s clusters changes with T and P in entation towarda tetrahedralconfigurationwith respect 6 the molecules in the spanning cluster. As shown in Fig.3, the spanning cluster appears as a fractalobject,withholesofanysize. Thesamelargedis- tribution of sizes characterizes also the finite clusters in the system. The absence of a characteristic size for the clusters (or the holes of the spanning cluster) is the con- sequence of the fluctuations at any length-scale, typical of a critical point. At T < T the majority of the molecules belongs to C a single percolating cluster that represents the network of tetrahedrally ordered molecules. All the other clus- ters are small, with a finite size that corresponds to the regions of correlated molecules. The presence of many small clusters gives a qualitative idea of the heterogene- ity of the dynamics at these temperatures. SUMMARY AND CONCLUSIONS Wedescribethenumericalsolutionofmeanfieldequa- tions and the implementation of the Wolff’s cluster MC algorithm for a cell model for liquid water. The mean field approach allows us to estimate in an approximate way the phase diagram of the model at any state point predicting intriguing new results at very low T [21]. Toexplorethe statepointsofinterestforthese predic- tions the use of standard simulations, such as molecular dynamics or Metropolis MC, is not effective due to the onset of the glassy dynamics [8]. To overcomethis prob- lem and access the deeply supercooled region of liquid water, we adopt the Wolff’s cluster MC algorithm. This method,indeed,allowstogreatlyacceleratetheautocor- relationtimeofthesystem. DirectcomparisonofWolff’s dynamicswithMetropolisdynamicsinthevicinityofthe liquid-liquidcriticalpoint showsa reductionofthe auto- correlationtime of a factor at least 104. Furthermore, the analysis of the clusters generated duringtheWolff’sMCdynamicsallowstoemphasizehow theregionsoftetrahedrallyorderedmoleculesbuildupon approachingtheliquid–liquidcriticalpoint,givingriseto the backbone of the tetrahedral hydrogen bond network at the phase transition [23]. The coexistence of clusters of correlated molecules with sizes that change with the state point gives a rationale for the heterogeneous dy- namics observed in supercooled water [22]. FIG. 4: Three snapshots of the system with N =100×100, ACKNOWLEDGMENTS showing the Wolff’s clusters of correlated water molecules. For each molecule we show the states of the four arms and associate different colors to different arm’s states. The state We thank Andrew Inglis for introducing one of the pointsareat pressureclose tothecritical valuePC (Pv0/ǫ= authors (MGM) to VPython, Francesco Mallamace for 0.72 ≃ PCv0/ǫ) and T > TC (top panel, kBT/ǫ = 0.0530), discussions, NSF grant CHE0616489 and Spanish MEC T ≃ TC (middle panel, kBT/ǫ = 0.0528), T < TC (bottom grant FIS2007-61433for support. panel,kBT/ǫ=0.0520), showing theonsetofthepercolation at T ≃TC. 7 [13] P.Kumar,G.FranzeseandH.E.Stanley,J.Phys.: Con- dens. Matter 20 (2008) 244114. [14] G. Franzese, M. I. Marqu´es, and H. Eugene Stanley, [1] Aspects of Physical Biology: Biological Water, Protein Phys. Rev.E 67 (2003) 011103. Solutions, Transport and Replication, G. Franzese and [15] M. I.Marqu´es, Phys. Rev.E 76 (2007) 021503. M. Rubieds. (Springer, Berlin, 2008). [16] M.A. Ricci, F. Bruni, and A. Giuliani Similarities be- [2] P.G.Debenedetti, J.Phys.: Condens. Matter15 (2003) tween confinedandsupercooled water,toappearonFara- R1669. day Discussion (2008). M. Chaplin “Water’s Hydrogen [3] P. G. Debenedetti and H. E. Stanley, “The Physics of Bond Strength”, cond-mat/0706.1355 (2007). SupercooledandGlassyWater,”PhysicsToday56[issue [17] The model does not differentiate “donor” molecule and 6] (2003) 40. “acceptor” molecule in the hydrogen bond definition. [4] G. Franzese, K. Stokely, X.-Q. Chu, P .Kumar, This simplification increases the number of possible M. G. Mazza, S.-H. Chen, and H. E. Stanley, J. Phys.: bonded configurations, hence increases the entropy as- Condens. Matter 20 (2008) 494210. sociatedtothelocaltetrahedralconfigurations.Asimple [5] C. A. Angell, Science 319 (2008) 582. modification of the model could explicitly take into ac- [6] H.E.Stanley,P.Kumar,G.Franzese,L.M.Xu,Z.Y.Yan, countthisfeature,howeverthecomparisonoftheresults M.G. Mazza, S.-H. Chen, F.Mallamace, S. V. Buldyrev, from the present version of the model with experiments “Liquid polyamorphism: Some unsolved puzzles of wa- and simulations from more detailed models shows good terin bulk,nanoconfined, and biological environments”, qualitative agreement. in Complex Systems, M. Tokuyama, I. Oppenheim, H. [18] V. Cataudella, G. Franzese, M. Nicodemi, A. Scala, and Nishiyama H, eds. AIP Conference Proceedings, 982 A. Coniglio, Phys. Rev. E 54 (1996) 175; G.Franzese, J. (2008) 251. Phys. A 29 (1996) 7367. [7] H. E. Stanley, S. V. Buldyrev, G. Franzese, N. Giovam- [19] A. Coniglio and F. Peruggi, J. Phys.A 15 (1982) 1873. battista, F. W. Starr, Phil. Trans. Royal Soc. 363, 509 [20] The results of [18, 19] guarantee that the cluster algo- (2005); P. Kumar, G. Franzese, S. V. Buldyrev, and H. rithmdescribedheresatisfiesthedetailedbalanceandis E. Stanley,Phys. Rev.E 73 (2006) 041505. ergodic. Therefore, it is a valid Monte Carlo dynamics. [8] P. Kumar, G. Franzese and H. E. Stanley, Phys. Rev. [21] M.G.Mazza,K.Stokely,H.E.Stanley,andG.Franzese, Lett.100 (2008) 105701. “Anomalous specific heat of supercooled water”, cond- [9] U.Wolff, Phys.Rev.Lett. 62 (1989) 361. mat/0807.4267 (2008). [10] G.FranzeseandH.E.Stanley,J.Phys.: Condens.Matter [22] M.G. Mazza et al., Phys. Rev. Lett. 96 (2006) 057803; 14 (2002) 2201; Physica A 314 (2002) 508. N. Giovambattista et al., J. Phys. Chem. B 108 (2004) [11] K. Stokely, M. G. Mazza, H. E. Stanley, G. Franzese, 6655; M.G.Mazzaetal.,Phys.Rev.E76(2007)031203. “Effect of hydrogen bond cooperativity on the behavior [23] A.Oleinikova,I.Brovchenko,J.Phys.: Condens.Matter of water” arXiv:0805.3468v1 (2008). 18 (2006) S2247. [12] G.FranzeseandH.E.Stanley,J.Phys.: Condens.Matter 19 (2007) 205126.