Cluster Monte Carlo and numerical mean field analysis for the water liquid--liquid phase transition
Marco G. Mazza, Kevin Stokely, Elena Strekalova, H. Eugene Stanley, Giancarlo Franzese
aa r X i v : . [ c ond - m a t . s o f t ] J a n Cluster Monte Carlo and numerical mean field analysis for the water liquid–liquidphase transition
Marco G. Mazza, Kevin Stokely, Elena Strekalova, H. Eugene Stanley, and Giancarlo Franzese Center for Polymer Studies and Department of Physics,Boston University, Boston, Massachusetts 02215, USA Departament 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 fieldapproach, we study the low temperature phase diagram of water, adopting a cell model that repro-duces the known properties of water in its fluid phases. Both methods allows us to study the waterthermodynamic behavior at temperatures where other numerical approaches –both Monte Carloand molecular dynamics– are seriously hampered by the large increase of the correlation times. Thecluster algorithm also allows us to emphasize that the liquid–liquid phase transition corresponds tothe percolation transition of tetrahedrally ordered water molecules.
PACS numbers: 61.20.Ja, 61.20.Gy
INTRODUCTION
Water is possibly the most important liquid for life[1] and, at the same time, is a very peculiar liquid [2].In the stable liquid regime its thermodynamic responsefunctions behave qualitatively differently than a typi-cal liquid. The isothermal compressibility K T , for ex-ample, has a minimum as a function of temperature at T = 46 ◦ C, while for a typical liquid K T monotonicallydecreases upon cooling. Water’s anomalies become evenmore pronounced as the system is cooled below the melt-ing point and enters the metastable supercooled regime[3].Different hypothesis have been proposed to rational-ize the anomalies of water [4]. All these interpretations,but one, predict the existence of a liquid–liquid phasetransition in the supercooled state, consistent with theexperiments to date [4] and supported by different mod-els [2].To discriminate among the different interpretations,many experiments have been performed [5]. However,the freezing in the temperature-range of interest can beavoided only for water in confined geometries or on thesurface of macromolecules [4, 6]. Since experiments inthe supercooled region are difficult to perform, numeri-cal simulations have played an important role in recentyears to help interpret the data. However, also the sim-ulations at very low temperature T are hampered by theglassy dynamics of the empirical models of water [7, 8].For these reasons is important to implement more effi-cient numerical simulations for simple models, able tocapture the fundamental physics of water but also lesscomputationally expensive. Here we introduce the im-plementation of a Wolff’s cluster algorithm [9] for theMonte Carlo (MC) simulations of a cell model for wa-ter [10]. The model is able to reproduce all the differ-ent scenarios proposed to interpret the behavior of wa-ter [11] and has been analyzed (i) with mean field (MF) [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] haveshown that very large times are needed to equilibratethe system as T →
0, as a consequence of the onset ofthe glassy dynamics. The implementation of the Wolff’sclusters MC dynamics, presented here, allows us to (i)drastically reduce the equilibration times of the model atvery low T and (ii) give a geometrical characterization ofthe regions of correlated water molecules (clusters) at low T and show that the liquid–liquid phase transition canbe interpreted as a percolation transition of the tetrahe-drally ordered clusters. THE MODEL
The system consists of N particles distributed within avolume V in d dimensions. The volume is divided into N cells of volume v i with i ∈ [1 , N ]. For sake of simplicity,these cells are chosen of the same size, v i = V /N , but thegeneralization to the case in which the volume can changewithout changes in the topology of the nearest–neighbor(n.n.) is straightforward. By definition, v i ≥ v , where v is the molecule hard-core volume. Each cell has avariable n i = 0 for a gas-like or n i = 1 for a liquid-likecell. We partition the total volume in a way such thateach cell has at least four n.n. cells, e.g. as in a cubiclattice in 3 d or a square lattice in 2 d . Periodic boundaryconditions are used to limit finite–size effects.The system is described by the Hamiltonian [10] H = − ǫ X h i,j i n i n j − J X h ij i n i n j δ σ ij ,σ ji + − J σ X i n i X ( k,l ) i δ σ ik ,σ il , (1)where ǫ > J >
FIG. 1: A pictorial representation of five water molecules in3 d . Two hydrogen bonds (grey links) connect the hydrogens(in blue) of the central molecule with the lone electrons (smallgray lines) of two nearest neighbor (n.n.) molecules. A bondindex (arm) with q = 6 possible values is associated to eachhydrogen and lone electron, giving rise to q possible orien-tational states for each molecule. A hydrogen bond can beformed only if the two facing arms of the n.n. molecules arein the same state. Arms on the same molecule interact amongthemselves to mimic the O-O-O interaction that drives themolecules toward a tetrahedral local structure. four (Potts) variables σ ij = 1 , . . . , q representing bond in-dices of molecule i with respect to the four n.n. molecules j , δ a,b = 1 if a = b and δ a,b = 0 otherwise, and h i, j i de-notes that i and j are n.n. The model does not assume aprivileged state for bond formation. Any time two facingbond indices (arms) are in the same (Potts) state, a bondis formed. The third term represents an intramolecular(IM) interaction accounting for the O–O–O correlation[16], locally driving the molecules toward a tetrahedralconfiguration. When the bond indices of a molecule arein the same state, the energy is decreased by an amount J σ > k, l ) i indicates one of the six different pairs of the fourbond indices of molecule i (Fig.1).Experiments show that the formation of a hydrogenbond leads to a local volume expansion [2]. Thus in oursystem the total volume is V = N v + N HB v HB , (2)where N HB ≡ X n i n j δ σ ij ,σ ji (3)is the total number of hydrogen bonds, and v HB is the constant specific volume increase due to the hydrogenbond formation. MEAN–FIELD ANALYSIS
In the mean–field (MF) analysis the macrostate ofthe system in equilibrium at constant pressure P andtemperature T ( N P T ensemble) may be determined bya minimization of the Gibbs free energy per molecule, g ≡ ( h H i + P V − T S ) /N w , where N w = X i n i (4)is the total number of liquid-like cells, and S = S n + S σ is the sum of the entropy S n over the variables n i andthe entropy S σ over the variables σ ij .A MF approach consists of writing g explicitly usingthe approximations X
1. There-fore, the number density n σ of bond indices σ ij is in thetetrahedral state is n σ = 1 + ( q − m σ q . (11)Since an appropriate form for h is [10] h = 3 J σ n σ , (12)we obtain that J σ q ≤ h ≤ J σ .The MF expressions for the entropies S n of the N vari-ables n i , and S σ of the 4 N n variables σ ij , are then [12] S n = − k B N ( n log( n ) + (1 − n ) log(1 − n )) (13) S σ = − k B N n [ n σ log( n σ )+(1 − n σ ) log(1 − n σ ) + log( q − , (14)where k B is the Boltzmann constant.Equating p σ ≡ n σ + (1 − n σ ) q − , (15)with the approximate expression in Eq. (10), allows forsolution of n σ , and hence g , in terms of the order param-eter m σ and n .By minimizing numerically the MF expression of g with respect to n and m σ , we find the equilibrium values n ( eq ) and m ( eq ) σ and, with Eqs. (4) and (2), we calcu-late the density ρ at any ( T, P ) and the full equationof state. An example of minimization of g is presentedin Fig. 2 where, for the model’s parameters J/ǫ = 0 . J σ /ǫ = 0 . v HB /v = 0 . q = 6, a discontinuityin m ( eq ) σ is observed for P v /ǫ > .
8. As discussed inRef.s [10, 14] this discontinuity corresponds to a first or-der phase transition between two liquid phases with dif-ferent degree of tetrahedral order and, as a consequence,different density. The higher P at which the changein m ( eq ) σ is continuous, corresponds to the pressure of aliquid–liquid critical point (LLCP). The occurrence of theLLCP is consistent with one of the possible interpreta-tions of the anomalies of water, as discussed in Ref. [12].However, for different choices of parameters, the modelreproduces also the other proposed scenarios [11]. THE SIMULATION WITH THE WOLFF’SCLUSTERS MONTE CARLO ALGORITHM
To perform MC simulations in the
N P T ensemble, weconsider a modified version of the model in which weallow for continuous volume fluctuations. To this goal,(i) we assume that the system is homogeneous with allthe variables n i set to 1 and all the cells with volume v = V /N ; (ii) we consider that V ≡ V MC + N HB v HB ,where V MC > N v is a dynamical variable allowed tofluctuate in the simulations; (iii) we replace the first (vander Waals) term of the Hamiltonian in Eq. (1) with aLennard-Jones potential with attractive energy ǫ > J and truncated at the hard-core distance U W ( r ) ≡ ( ∞ if r r , ǫ h(cid:0) r r (cid:1) − (cid:0) r r (cid:1) i if r > r . (16) Tetrahedral Order Paramter m σ -1.95-1.9 M o l a r G i bb s F r ee E n e r gy g [ ε ] Pv / ε = 0.7 k B T/ ε = 0.06 k B T/ ε = 0.08 Tetrahedral Order Paramter m σ -1.75-1.7 M o l a r G i bb s F r ee E n e r gy g [ ε ] Pv / ε = 0.8 k B T/ ε = 0.05 k B T/ ε = 0.07 Tetrahedral Order Paramter m σ -1.55-1.5 M o l a r G i bb s F r ee E n e r gy g [ ε ] Pv / ε = 0.9 k B T/ ε = 0.04 k B T/ ε = 0.06 FIG. 2: Numerical minimization of the molar Gibbs free en-ergy g in the mean field approach. The model’s parametersare J/ǫ = 0 . J σ /ǫ = 0 . v HB /v = 0 . q = 6. In eachpanel we present g (dashed lines) calculated at constant P and different values of T . The thick line crossing the dashedlines connects the minima m ( eq ) σ of g at different T . Upperpanel: P v /ǫ = 0 .
7, for T going from k B T /ǫ = 0 .
06 (top) to k B T /ǫ = 0 .
08 (bottom). Middle panel:
P v /ǫ = 0 .
8, for T going from k B T /ǫ = 0 .
05 (top) to k B T /ǫ = 0 .
07 (bottom).Lower panel:
P v /ǫ = 0 .
9, for T going from k B T /ǫ = 0 . k B T /ǫ = 0 .
06 (bottom). In each panel dashed linesare separated by k B δT /ǫ = 0 . m ( eq ) σ increases when T decreases, being 0 (marking the absence oftetrahedral order) at the higher temperatures and ≃ . T , m ( eq ) σ changes in a continuous way for P v /ǫ = 0 . . P v /ǫ = 0 . P . where r ≡ ( v ) /d ; the distance between two n.n.molecules is ( V /N ) /d , and the distance r between twogeneric molecules is the Cartesian distance between thecenter of the cells in which they are included.The simplification (i) could be removed, by allowingthe cells to assume different volumes v i and keeping fixedthe number of possible n.n. cells. However, the resultsof the model under the simplification (i) compares wellwith experiments [12]. Furthermore, the simplification(i) allows to drastically reduce the computational costof the evaluation of the U W ( r ) term from N ( N −
1) to N − N = 2500 and N = 10000 molecules, each with four n.n. molecules,at constant P and T , in 2d, and with the same parame-ters used for the mean field analysis. To each moleculeswe associate a cell on a square lattice. The Wolff’s algo-rithm is based on the definition of a cluster of variableschosen in such a way to be thermodynamically correlated[18, 19]. To define the Wolff’s cluster, a bond index (arm)of a molecule is randomly selected; this is the initial el-ement of a stack. The cluster is grown by first checkingthe remaining arms of the same initial molecule: if theyare in the same Potts state, then they are added to thestack with probability p same ≡ min [1 , − exp( − βJ σ )][9], where β ≡ ( k B T ) − . This choice for the proba-bility p same depends on the interaction J σ between twoarms on the same molecule and guarantees that theconnected arms are thermodynamically correlated [19].Next, the arm of a new molecule, facing the initiallychosen arm, is considered. To guarantee that connectedfacing arms correspond to thermodynamically correlatedvariables, is necessary [18] to link them with the probabil-ity p facing ≡ min [1 , − exp( − βJ ′ )] where J ′ ≡ J − P v HB is the P –dependent effective coupling between two facingarms as results from the enthalpy H + P V of the sys-tem. It is important to note that J ′ can be positive ornegative depending on P . If J ′ > p facing ; if J ′ < p facing [20]. Only after everypossible direction of growth for the cluster has been con-sidered the values of the arms are changed in a stochasticway; again we need to consider two cases: (i) if J ′ > σ new = (cid:0) σ old + φ (cid:1) mod q (17)where φ is a random integer between 1 and q ; (ii) if J ′ <
0, the state of every single arm is changed (rotated) by the same random constant φ ∈ [1 , . . . q ] σ new i = (cid:0) σ old i + φ (cid:1) mod q. (18)In order to implement a constant P ensemble we letthe volume fluctuate. A small increment ∆ r/r = 0 . V ≡ V new − V old and van der Waals energy ∆ E W is computed and the move is accepted with probabilitymin (1 , exp [ − β (∆ E W + P ∆ V − T ∆ S )]), where ∆ S ≡− N k B ln( V new /V old ) is the entropic contribution. MONTE CARLO CORRELATION TIMES
The cluster MC algorithm described in the previoussection turns out to be very efficient at low T , allow-ing to study the thermodynamics of deeply supercooledwater with quite intriguing results [21]. To estimate theefficiency of the cluster MC dynamics with respect to thestandard Metropolis MC dynamics, we evaluate in bothdynamics, and compare, the autocorrelation function ofthe average magnetization per site M i ≡ P j σ ij , wherethe sum is over the four bonding arms of molecule i . C M ( t ) ≡ N X i h M i ( t + t ) M i ( t ) i − h M i i h M i i − h M i i . (19)For sake of simplicity, we define the MC dynamics au-tocorrelation time τ as the time, measured in MC steps,when C M ( τ ) = 1 /e . Here we define a MC step as 4 N up-dates of the bond indices followed by a volume update,i.e. as 4 N + 1 steps of the algorithm.In Fig. 3 we show a comparison of C M ( t ) for theMetropolis and Wolff algorithm implementations of thismodel for a system with N = 50 ×
50, at three tempera-tures along an isobar below the LLCP, and approachingthe line of the maximum, but finite, correlation length,also known as Widom line T W ( P ) [12]. In the top panel,at T ≫ T W ( P ) ( k B T /ǫ = 0 . P v /ǫ = 0 . τ W ≈ × , and for the Metropolis dynamics τ M ≈ .In the middle panel, at T > T W ( P ) ( k B T /ǫ = 0 . P v /ǫ = 0 .
6) the difference between the two correlationtimes is larger: τ W ≈ . × , τ M ≈ × . The bot-tom panel, at T ≃ T W ( P ) ( k B T /ǫ = 0 . P v /ǫ = 0 . τ W ≈ . × , while τ M is beyond the accessibletime window ( τ M > ).Since as T → τ M /τ W grows at lower T allowing the evalua-tion of thermodynamics averages even at T ≪ T C [21].In particular, the cluster MC algorithm turns out to bevery efficient when approaching the Widom line in thevicinity of the LLCP, with an efficiency of the order of10 . We plan to analyze in a systematic way how the Time t [MC steps] C o rr e l a ti on f un c ti on C M ( t ) MetropolisWolff k B T/ ε =0.11 Time t [MC steps] C o rr e l a ti on f un c ti on C M ( t ) MetropolisWolff k B T/ ε =0.09 Time t [MC steps] C o rr e l a ti on f un c ti on C M ( t ) MetropolisWolff k B T/ ε =0.06 FIG. 3: Comparison of the autocorrelation function C M ( t ) forthe Metropolis (circles) and Wolff (squares) implementation ofthe present model. We show the temperatures k B T /ǫ = 0 . k B T /ǫ = 0 .
09 (middle panel), k B T /ǫ = 0 . P v /ǫ = 0 . N = 50 × efficiency τ M /τ W grows on approaching the LLCP. Thisresult is well known for the standard liquid-gas criticalpoint [9] and, on the basis of our results, could be ex-tended also to the LLCP. However, this analysis is veryexpensive in terms of CPU time and goes beyond the goalof the present work. Nevertheless, the percolation analy-sis, presented in the next section, helps in understandingthe physical reason for this large efficiency.The efficiency is a consequence of the fact that the av-erage size of Wolff’s clusters changes with T and P in the same way as the average size of the regions of cor-related molecules [19], i.e. a Wolff’s cluster statisticallyrepresents a region of correlated molecules. Moreover,the mean cluster size diverges at the critical point withthe same exponent of the Potts magnetic susceptibility[19], and the clusters percolate at the critical point, aswe will discuss in the next section. PERCOLATING CLUSTERS OF CORRELATEDMOLECULES
The efficiency of the Wolff’s cluster algorithm is a con-sequence of the exact relation between the average sizeof the finite clusters and the average size of the regionsof thermodynamically correlated molecules. The proofof this relation at any T derives straightforward from theproof for the case of Potts variables [19]. This relationallows to identify the clusters built during the MC dy-namics with the correlated regions and emphasizes (i) theappearance of heterogeneities in the structural correla-tions [22], and (ii) the onset of percolation of the clustersof tetrahedrally ordered molecules at the LLCP [23], asshown in Fig. 4.A systematic percolation analysis [18] is beyond thegoal of this report, however configurations such as thosein Fig. 4 allow the following qualitative considerations.At T > T C the average cluster size is much smaller thanthe system size. Hence, the structural correlations amongthe molecules extends only to short distances. This sug-gests that the correlation time of a local dynamics, suchas Metropolis MC or molecular dynamics, would be shorton average at this temperature and pressure. Neverthe-less, the system appears strongly heterogeneous with thecoexistence of large and small clusters, suggesting thatthe distribution of correlation times evaluated amongmolecules at a given distance could be strongly heteroge-neous. The clusters appear mostly compact but with afractal surface, suggesting that borders between clusterscan rapidly change.At T ≃ T C there is one large cluster, in red on the rightof the middle panel of Fig. 4, with a linear size compa-rable to the system linear extension and spanning in thevertical direction. The appearance of spanning clustersshows the onset of the percolation geometrical transition.At this state point the correlation time of local, such asMetropolis MC dynamics or molecular dynamics wouldbe very slow as a consequence of the large extension ofthe structurally correlated region. On the other hand,the correlation time of the Wolff’s cluster dynamics isshort because it changes in one single MC step the stateof all the molecules in clusters, some of them with verylarge size. Once the spanning cluster is formed, it breaksthe symmetry of the system and a strong effective fieldacts on the molecules near its border to induce their reori-entation toward a tetrahedral configuration with respect FIG. 4: Three snapshots of the system with N = 100 × P C ( P v /ǫ =0 . ≃ P C v /ǫ ) and T > T C (top panel, k B T /ǫ = 0 . T ≃ T C (middle panel, k B T /ǫ = 0 . T < T C (bottompanel, k B T /ǫ = 0 . T ≃ T C . the molecules in the spanning cluster.As shown in Fig.3, the spanning cluster appears as afractal object, with holes of any size. The same large dis-tribution of sizes characterizes also the finite clusters inthe system. The absence of a characteristic size for theclusters (or the holes of the spanning cluster) is the con-sequence of the fluctuations at any length-scale, typicalof a critical point.At T < T C the majority of the molecules belongs toa single percolating cluster that represents the networkof tetrahedrally ordered molecules. All the other clus-ters are small, with a finite size that corresponds to theregions of correlated molecules. The presence of manysmall clusters gives a qualitative idea of the heterogene-ity of the dynamics at these temperatures. SUMMARY AND CONCLUSIONS
We describe the numerical solution of mean field equa-tions and the implementation of the Wolff’s cluster MCalgorithm for a cell model for liquid water. The meanfield approach allows us to estimate in an approximateway the phase diagram of the model at any state pointpredicting intriguing new results at very low T [21].To explore the state points of interest for these predic-tions the use of standard simulations, such as moleculardynamics or Metropolis MC, is not effective due to theonset of the glassy dynamics [8]. To overcome this prob-lem and access the deeply supercooled region of liquidwater, we adopt the Wolff’s cluster MC algorithm. Thismethod, indeed, allows to greatly accelerate the autocor-relation time of the system. Direct comparison of Wolff’sdynamics with Metropolis dynamics in the vicinity of theliquid-liquid critical point shows a reduction of the auto-correlation time of a factor at least 10 .Furthermore, the analysis of the clusters generatedduring the Wolff’s MC dynamics allows to emphasize howthe regions of tetrahedrally ordered molecules build up onapproaching the liquid–liquid critical point, giving rise tothe backbone of the tetrahedral hydrogen bond networkat the phase transition [23]. The coexistence of clustersof correlated molecules with sizes that change with thestate point gives a rationale for the heterogeneous dy-namics observed in supercooled water [22]. ACKNOWLEDGMENTS
We thank Andrew Inglis for introducing one of theauthors (MGM) to VPython, Francesco Mallamace fordiscussions, NSF grant CHE0616489 and Spanish MECgrant FIS2007-61433 for support. [1]
Aspects of Physical Biology: Biological Water, ProteinSolutions, Transport and Replication , G. Franzese andM. Rubi eds. (Springer, Berlin, 2008).[2] P. G. Debenedetti, J. Phys.: Condens. Matter 15 (2003)R1669.[3] P. G. Debenedetti and H. E. Stanley, “The Physics ofSupercooled and Glassy Water,” Physics Today 56 [issue6] (2003) 40.[4] G. Franzese, K. Stokely, X.-Q. Chu, P .Kumar,M. G. Mazza, S.-H. Chen, and H. E. Stanley, J. Phys.:Condens. Matter 20 (2008) 494210.[5] C. A. Angell, Science 319 (2008) 582.[6] H.E. Stanley, P. Kumar, G. Franzese, L.M. Xu, Z.Y. Yan,M.G. Mazza, S.-H. Chen, F.Mallamace, S. V. Buldyrev,“Liquid polyamorphism: Some unsolved puzzles of wa-ter in bulk, nanoconfined, and biological environments”,in
Complex Systems , M. Tokuyama, I. Oppenheim, H.Nishiyama H, eds. AIP Conference Proceedings, 982(2008) 251.[7] H. E. Stanley, S. V. Buldyrev, G. Franzese, N. Giovam-battista, F. W. Starr, Phil. Trans. Royal Soc. 363, 509(2005); P. Kumar, G. Franzese, S. V. Buldyrev, and H.E. Stanley, Phys. Rev. E 73 (2006) 041505.[8] P. Kumar, G. Franzese and H. E. Stanley, Phys. Rev.Lett. 100 (2008) 105701.[9] U. Wolff, Phys. Rev. Lett. 62 (1989) 361.[10] G. Franzese and H. E. Stanley, J. Phys.: Condens. Matter14 (2002) 2201; Physica A 314 (2002) 508.[11] K. Stokely, M. G. Mazza, H. E. Stanley, G. Franzese,“Effect of hydrogen bond cooperativity on the behaviorof water” arXiv:0805.3468v1 (2008).[12] G. Franzese and H. E. Stanley, J. Phys.: Condens. Matter19 (2007) 205126. [13] P. Kumar, G. Franzese and H. E. Stanley, J. Phys.: Con-dens. Matter 20 (2008) 244114.[14] G. Franzese, M. I. Marqu´es, and H. Eugene Stanley,Phys. Rev. E 67 (2003) 011103.[15] M. I. Marqu´es, Phys. Rev. E 76 (2007) 021503.[16] M.A. Ricci, F. Bruni, and A. Giuliani