A tetrahedral entropy for water
 ^{a}Center for Studies in Physics and Biology, The Rockefeller University, 1230 York Avenue, New York, NY 10021;
 ^{b}Department of Physics, Yeshiva University, 500 West 185th Street, New York, NY 10033; and
 ^{c}Center for Polymer Studies and Department of Physics, Boston University, Boston, MA 02215
See allHide authors and affiliations

Contributed by H. Eugene Stanley, October 5, 2009 (received for review May 15, 2009)
Abstract
We introduce the spacedependent correlation function C _{Q}(r) and timedependent autocorrelation function C _{Q}(t) of the local tetrahedral order parameter Q ≡ Q(r,t). By using computer simulations of 512 waterlike particles interacting through the transferable interaction potential with five points (TIP5 potential), we investigate C _{Q}(r) in a broad region of the phase diagram. We find that at low temperatures C _{Q}(t) exhibits a twostep timedependent decay similar to the selfintermediate scattering function and that the corresponding correlation time τ_{Q} displays a dynamic crossover from nonArrhenius behavior for T > T _{W} to Arrhenius behavior for T < T _{W}, where T _{W} denotes the Widom temperature where the correlation length has a maximum as T is decreased along a constantpressure path. We define a tetrahedral entropy S _{Q} associated with the local tetrahedral order of water molecules and find that it produces a major contribution to the specific heat maximum at the Widom line. Finally, we show that τ_{Q} can be extracted from S _{Q} by using an analog of the Adam–Gibbs relation.
 specific heat of water
 orientational entropy of tetrahedral liquids
 anomalies of liquid water
 Widom line
It has long been appreciated that the local structure around a molecule of liquid water arising from the vertices formed by four nearest neighbors is approximately tetrahedral at ambient pressure and that the degree of tetrahedrality increases when water is cooled (1–5). An important advance occurred in the past 10 years when computer simulations allowed the quantification of the degree of tetrahedrality (6–10) by assigning to each molecule a local tetrahedral order parameter Q (11–16). At high temperatures, the probability distribution P(Q,T) is bimodal, with one peak corresponding to a high degree of tetrahedrality and the other to a less tetrahedral environment (Fig. 1). Upon decreasing temperature, the peak associated with a high degree of tetrahedrality grows, suggesting that the local structure of water becomes much more tetrahedral at lower temperatures (13, 14, 17).
Introduction
Water has been hypothesized to belong to the class of polymorphic liquids, phase separating—at sufficiently low temperatures and high pressures—into two distinct liquid phases: a high density liquid (HDL) with smaller Q and a low density liquid (LDL) with larger Q (18). The coexistence line separating these two phases may terminate at a liquid–liquid (LL) critical point, above which (in the LL supercritical region) appears a line of correlation length maximum in the pressure–temperature plane. The locus of maximum correlation length in the onephase region is called the Widom line T _{W} ≡ T _{W}(P) (19), near which different response functions, such as isobaric specific heat C _{P} and isothermal compressibility K _{T}, display maxima. Recent neutronscattering experiments (20) and computer simulations (19) show that the dynamics of water gradually cross over from being nonArrhenius for T > T _{W} to Arrhenius for T < T _{W}.
Although dynamic heterogeneities have received considerable attention (21, 22), static heterogeneities have attracted less interest. In this paper, we ask how the tetrahedral order and its spatiotemporal correlations change upon crossing the Widom line. To this end, we introduce three quantities of physical interest: (i) the spacedependent correlation function C _{Q}(r), (ii) the timedependent autocorrelation function C _{Q}(t), and (iii) the tetrahedral entropy S _{Q}. These functions can also be usefully applied to study other locally tetrahedral liquids, such as silicon (23–25), silica (26), and phosphorus (27).
Model
We perform molecular dynamics (MD) simulations of N = 512 waterlike molecules interacting via the transferable interaction potential with five points (TIP5 potential), (28), which exhibits a LL critical point at T _{C} ≈ 217 K and P _{C} ≈ 340 MPa (29, 30). We carry out simulations in the NPT ensemble at atmospheric pressure (P = 1 atm) for temperatures ranging from 320 K down to 230 K.
Results
Space and TimeDependent Correlations of the Tetrahedral Order Parameter.
To quantify the local degree of tetrahedrality, we calculate the local tetrahedral order parameter (13) where ψ_{ikj} is the angle formed by the molecule k and its nearest neighbors i and j. The average value 〈Q〉 ≡(1/N) ∑_{k} Q _{k} increases with decreasing T and saturates at lower T, whereas d〈Q〉/dT has a maximum at the Widom line T _{W} ≈ 250 K (31–33).
To characterize the spacedependent correlations of the local order parameter, we find all the pairs of molecules i and j whose oxygens are separated by distances belonging to the interval [r −Δr/2,r + Δr/2]. The number of such pairs is where r _{ij} is the distance between the oxygens of molecules i and j. The sum is taken over all molecules in the system and For Δr ← 0, δ(r _{ij} − r, Δr) ← δ(r _{ij} − r)Δr where δ(r _{ij} − r) is the Dirac δfunction. N(r, Δr) can be approximated as where N is the total number of molecules in the system, g _{OO}(r) is the oxygen–oxygenpair correlation function, and ρ is the number density.
Next, we find the local order parameters Q _{k} and Q _{l} for all such pairs of molecules and compute their mean, their variance, and covariance, Finally, we introduce the spacedependent correlation function of the local tetrahedrality Q,
In Figs. 2 A and 2 B, we show Q(r) and σ_{Q} ^{2}(r) for different temperatures and atmospheric pressures. The behavior of Q(r) and its variance σ_{Q} ^{2}(r) as a function of the distance r has a clear physical meaning. For molecules separated by 0.32 nm, Q(r) has a deep minimum characterizing the distortion of the first tetrahedral coordination shell of the four nearest neighbors by the intrusion of a “fifth neighbor” (34). Conversely, the quantity σ_{Q} ^{2}(r = 0.32 nm) −σ_{Q} ^{2}(∞) for such molecules dramatically increases upon decreasing temperature (see Fig. 2 C), and has a maximum at T ≈ 246 K, which is approximately equal to the ambient pressure value of T _{W} (≈ 250 K for the TIP5P model) reported in refs. 19 and 31–33. A measure qualitatively similar to Q(r) has been studied in case of solvation structure of water around carbohydrate molecules (35). In Fig. 2 D, we show C _{Q}(r) for different temperatures for P = 1 atm. C _{Q}(r) has positive maxima at positions of the first and second peaks in the oxygen–oxygen pair correlation function. Water molecules separated by r ≈ 0.32 nm exhibit weak anticorrelation in local tetrahedral order at high temperatures, which changes to positive correlation upon decreasing the temperature below T _{W}.
To study the time development of the local tetrahedral order parameter, we introduce the timedependent autocorrelation function where Q _{i}(t) is the tetrahedral order parameter of each molecule in the system and 〈Q _{i}〉 is the ensemble average. In Fig. 3 A we show C _{Q}(t) for different temperatures. The decay of C _{Q}(t) is reminiscent of the decay of the selfintermediate scattering function. The longtime behavior of C _{Q}(t) is exponential at high T but at low T can be fit with a stretched exponential exp[−(t/τ)^{β}], where 0 < β < 1 (36).
For simplicity, we define the correlation time τ_{Q} as the time required for C _{Q}(t) to decay by a factor e. Fig. 3 B shows the values of τ_{Q} as function of 1/T on an Arrhenius plot. The behavior of τ_{Q} is nonArrehenius at high temperatures and can be fit by a power law B(T − T _{MCT})^{−γ} where B = 25.39, the mode coupling temperature (37) T _{MCT} ≈ 246.186, and γ = 1.17722 are the fitting parameters. At low T, τ_{Q} deviates from the powerlaw fit and becomes Arrhenius. This crossover in relaxation of local tetrahedral order occurs near T _{W}, similar to the crossover found in density relaxation (19, 20, 38).
Tetrahedral Entropy and Tetrahedral Specific Heat.
Because Q measures the local tetrahedral order, it must contribute to the entropy of the system. We next derive an expression for this “tetrahedral entropy” S(Q _{1},Q _{2},…,Q _{N}), which we define to be the logarithm of the number of states corresponding to the interval between (Q _{1},Q _{2},….,Q _{N}) and (Q _{1} + ΔQ _{1},Q _{2} +Δ Q _{2},…,Q _{N} + ΔQ _{N}) in the hypercubic space formed by Qs. According to Eq. 1,
We further define a measure of “tetrahedral specific heat” as C _{P} ^{Q} ≡ T(∂S _{Q}(T)/∂T)_{P}. Fig. 4 B shows the temperature dependence of C _{P} ^{Q} and total specific heat C _{P} ^{Total} for P = 1 atm. C _{P} ^{Q} shows a maximum at T ≈ 250 ± 10 K ≈ T _{W}, where the total specific heat C _{P} ^{Total} also has a maximum (19, 31, 32, 33), suggesting that the fluctuations in tetrahedral order reach a maximum at T _{W}(P). Moreover, comparing the tetrahedral specific heat in Fig. 4 B with the total specific heat, we find that the C _{P} ^{Q} is a major contribution to the C _{P} ^{Total} and, particularly, the difference ΔC _{P} of the two specific heats remains a constant within the error bar for all temperatures studied, hence suggesting that C _{P} ^{Q} is responsible for the Widomline transformations.
To relate the tetrahedral entropy S _{Q}(T) to the tetrahedral relaxation time τ_{Q}(T) associated with the local tetrahedral ordering, we propose the following generalization of the Adam–Gibbs relation between the translational relaxation time and configurational entropy (39, 40), where τ_{Q}(0) is the tetrahedral relaxation time at very large T, and A is a parameter playing the role of activation energy. Accordingly, we calculate τ_{Q}(T) (Fig. 3 B) from S _{Q}(T) by using Eq. 13 with three free parameters: A, τ(0), and S _{0}. We find that S _{0} ≈ 4.21k _{B}. The values of τ_{Q}(T) calculated by using Eq. 13 are within error bars of the values of τ_{Q}(T) calculated from the autocorrelation function (see Fig. 3 B).
Fig. 3B is an Arrhenius plot of τ_{Q}, calculated from Eq. 13 . The temperature dependence of τ_{Q} is different at low and high temperatures, changing from nonArrhenius (a Tdependent slope on the Arrhenius plot), which can be fit by a power law at high temperatures, to Arrhenius (a constant slope) at low temperatures (20, 19, 31–33).
Relation of Tetrahedral Entropy to Translational Order Parameter and Translational Entropy.
In this section, we investigate the relation of S _{Q} with (i) the translational order parameter t* (41) defined as and (ii) the translational entropy S ^{2} (42, 43) obtained from the twopoint translational correlations, which is defined as where ρ is density and g(r) is the oxygen–oxygen radial distribution function.
The excess entropy S ^{2} describes the contribution of twopoint correlations to the total entropy of the system. It has been found that S ^{2} agrees very well with the total entropy in different systems including some models of water (43, 44).
In Fig. 5 A, we show S _{Q} as a function of average tetrahedral order parameter Q̄ for different temperatures and atmospheric pressure. S _{Q} varies approximately as log(1 − Q̄) for the range of temperature studied. To justify this dependence of S _{Q} on Q̄, we expand Eq. 12, where F is a function of moments of the fluctuations (Q − Q̄). Hence, if we ignore F compared with the first term, then S _{Q} ∼ log(1 − Q̄). In Fig. 5 A, we plot the above function (dashed line), assuming a constant F = −0.40 for all temperatures studied (230 – 320 K) and find that it agrees well with the values of S _{Q} computed by using Eq. 12 at low T but deviates slightly at high T.
Furthermore, Fig. 5 B shows that S _{Q} varies linearly with t*. Combining the results of Figs. 5 A and B, we form Fig. 5 C, which demonstrates s simple relation between tetrahedral and translational oder parameters. We next find that, at low T, S _{Q} changes linearly with S ^{2}; however it is nonlinear for high T. Note that a simple nonlinear relation between Q and t* was found for the SPC/E model of water in the structurally anomalous region (13). Moreover, Fig. 5 D shows that the temperature where this change in behavior occurs roughly coincides with the temperature of the Widom line T _{W}.
Discussion and Summary
In summary, we have studied the space and time correlations of local tetrahedral order Q, presumably related to local tetrahedral heterogeneities. We find that the spacedependent correlation of the the local tetrahedral order is anticorrelated for the molecules separated by 3.2 Å at high temperatures. This negative correlation changes to positive correlation upon decreasing T below the Widom temperature T _{W}(P). Further, we define a measure of the tetrahedral entropy S _{Q} and find that the specific heat associated with the tetrahedral ordering is responsible for the Widomline transformations. Finally, we find that S _{Q} well describes tetrahedral relaxation using the Adam–Gibbs relation. Moreover, we find a simple approximate relation between tetrahedral and translational order parameters Q̄ and t*. On the completion of this work, we were made aware by a referee that Lazardis and Karplus have studied the orientational entropy calculated from the twopoint angular correlation function (45). Although they find that the orientational entropy of water has a major contribution to the total entropy, a contribution of orientational fluctuations to the total specific heat was not studied. The tetrahedral entropy S(Q) defined in this paper captures the twopoint orientational correlations but is much easier to compute in contrast to the computation of orientational entropy from the twopoint angular correlation function and hence offers a rather simpler way to investigate the entropy associated with local structures.
Acknowledgments
We thank J. R. Errington, P. H. Poole, S. Sastry, and F. W. Starr for fruitful discussions and the National Science Foundation Chemistry Division for support. P.K. acknowledges support from Keck Foundation National Academies Keck Futures Initiatives Award, and S.V.B. acknowledges partial support through the Dr. Bernard W. Gamson Computational Science Center at Yeshiva College.
Footnotes
 ^{1}To whom correspondence should be addressed. Email: pradeep.kumar{at}rockefeller.edu or hes{at}bu.edu

Author contributions: P.K., S.V.B., and H.E.S. designed research, performed research, contributed new reagents/analytic tools, analyzed data, and wrote the paper.

The authors declare no conflict of interest.

Freely available online through the PNAS open access option.
References
 ↵
 ↵
 Pople JA
 ↵
 ↵
 Eisenberg D,
 Kauzmann W
 ↵
 Debenedetti PG
 ↵
 ↵
 ↵
 ↵
 ↵
 Kusalik PG,
 Svishchev IM
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Xu L,
 et al.
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 Katayama Y,
 et al.
 ↵
 ↵
 ↵
 ↵
 Kumar P,
 et al.
 ↵
 Kumar P
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
Citation Manager Formats
Article Classifications
 Physical Sciences
 Physics