Table Of ContentCluster 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.
