# Water-like anomalies as a function of tetrahedrality

See allHide authors and affiliations

Edited by Pablo G. Debenedetti, Princeton University, Princeton, NJ, and approved February 27, 2018 (received for review December 21, 2017)

## Significance

Water is the most common and yet least understood material on Earth. Despite its simplicity, water tends to form tetrahedral order locally by directional hydrogen bonding. This structuring is known to be responsible for a vast array of unusual properties, e.g., the density maximum at 4 °C, which play a fundamental role in countless natural and technological processes, with the Earth’s climate being one of the most important examples. By systematically tuning the degree of tetrahedrality, we succeed in continuously interpolating between water-like behavior and simple liquid-like behavior. Our approach reveals what physical factors make water so anomalous and special even compared with other tetrahedral liquids.

## Abstract

Tetrahedral interactions describe the behavior of the most abundant and technologically important materials on Earth, such as water, silicon, carbon, germanium, and countless others. Despite their differences, these materials share unique common physical behaviors, such as liquid anomalies, open crystalline structures, and extremely poor glass-forming ability at ambient pressure. To reveal the physical origin of these anomalies and their link to the shape of the phase diagram, we systematically study the properties of the Stillinger–Weber potential as a function of the strength of the tetrahedral interaction λ. We uncover a unique transition to a reentrant spinodal line at low values of λ, accompanied with a change in the dynamical behavior, from non-Arrhenius to Arrhenius. We then show that a two-state model can provide a comprehensive understanding on how the thermodynamic and dynamic anomalies of this important class of materials depend on the strength of the tetrahedral interaction. Our work establishes a deep link between the shape of the phase diagram and the thermodynamic and dynamic properties through local structural ordering in liquids and hints at why water is so special among all substances.

- tetrahedral liquids
- water’s anomalies
- water-like liquids
- two-state model
- modified Stillinger–Weber potential

Liquids do not possess long-range order but often have short-range order. For example, water, silicon, germanium, and carbon are known to form tetrahedral order locally because of the directional nature of hydrogen or covalent bonding. These liquids commonly exhibit anomalous thermodynamic and dynamic anomalies, which are absent in ordinary liquids, e.g., van der Waals liquids. Liquid anomalies include the density maximum as a function of temperature T, the steep increase in the isothermal compressibility and heat capacity upon cooling, the non-Arrhenius behavior of viscosity and diffusion constant at low pressures, and the minimum of viscosity and the maximum of diffusion constant as a function of pressure P (see, e.g., refs. 1⇓⇓⇓–5 for water anomalies). Furthermore, all these liquids commonly have V-shaped P–T solid–liquid phase diagrams, in which the melting point has a minimum at a positive pressure

As a coarse-grained classical model for tetrahedral materials, the Stillinger–Weber (SW) potential has emerged as an effective model potential capable of capturing all of the relevant physical properties that stem from the tetrahedrality of the interactions. The original parameterization of the SW potential was targeted to the bulk properties of silicon (7), but has also found widespread applicability in the modeling of other group XIV elements (8). Apart from atomic fluids, the SW potential has also found application in the coarse-grained description of complex molecular fluids. The most notable example is water, whose SW representation is intermediate between that of silicon and carbon and is known as monoatomic water (mW model) (8). While retaining a high degree of structural accuracy, the mW model has proved to be very efficient from a computational point of view and has played a big role in the study of water crystallization (9⇓⇓⇓–13), which otherwise requires advanced techniques (14, 15). As a good model of water, the mW water exhibits a vast array of thermodynamic and dynamic anomalies (16⇓⇓–19), and recently the behavior of the anomalies for different values of λ was considered in refs. 20 and 21. Ref. 20 focused on the location of the second critical point, studied by means of the isochore crossing technique, showing that changing λ can decrease the critical pressure to ambient conditions and down to the liquid–vapor spinodal, as in the critical point-free scenario. In ref. 21, the full hierarchy of anomalies was considered for three different values of λ, showing that they follow a silica-like hierarchy, which becomes a water-like hierarchy if the excess entropy and Rosenfeld scaling are considered. These works have shown the richness of the behavior of the SW model and opened the question of whether we can rationalize the anomalous behavior of tetrahedral liquids and whether we can connect their behavior to the underlying phase diagram.

In this article we consider the anomalous behavior of the liquid phase, focusing in particular on the liquid anomalies that occur at negative pressures. Our goal is to connect the behavior of liquid anomalies with the change of thermodynamic properties as a function of λ. We start by computing the full phase diagrams in the extended (T, P, λ) thermodynamic space, extending the results of ref. 22 to negative pressures, where clathrate structures are the stable crystals. The region at negative pressure is crucial for unveiling the origin of the anomalous behavior, and by measuring the density fluctuations, we track the stability limit of the liquid, i.e., the liquid-to-gas spinodal line. We find evidence for a transition from a positively sloped spinodal line to a reentrant spinodal as a function of the λ parameter, providing a unique example of such a transition in a water-like model. We show how this result is connected to the anomalous phase behavior of water and argue that a two-state modeling of the liquid phase (6, 16, 23⇓⇓⇓⇓⇓⇓–30) provides a simple theoretical framework which rationalizes the anomalous behavior of tetrahedral liquids as a function of the strength of the tetrahedral interaction. The model provides a deep link between the anomalies and the shape of the phase diagram, as a consequence of the fact that locally favored structures in liquids have the same local symmetry as the low-pressure diamond crystal. More precisely, we reveal that the value of λ corresponding to water maximizes two-state features and the resulting anomalies, providing structural flexibility to water: Water can change its physical and chemical properties by changing an extra structural degree of freedom, i.e., the fraction of the two states, in response to external perturbations.

## Results

### Two-State Model.

Liquid anomalies can be divided into two categories: thermodynamic and dynamic anomalies. Thermodynamic ones originate from the anomalous temperature dependence of a thermodynamic response function: Unlike the ordinary behavior of simple liquids, in water, thermodynamic fluctuations show an increase with lowering the temperature. An example is given by the isothermal compressibility

To rationalize both thermodynamic and dynamic anomalies, we use a two-state model. The history of two-state models of water dates back to Röntgen (31). The basic idea is that the anomalies of water can be understood if water is described as a mixture of two components in thermodynamic equilibrium, such that the concentration of the mixture is state dependent. Until recently, however, water was described as a mixture of distinct structural components, whose number is two (32⇓–34) to four (35).

Only recently, the importance of the degeneracy of states (or the large entropic loss upon the formation of locally favored tetrahedral structures) was properly recognized (23⇓–25). Furthermore, unlike previous approaches, where the order parameter is only density, it was proposed (36, 37) that we need at least two order parameters to understand the phenomena: one is the density ρ and the other is bond order s, which represents the local breakdown of rotational symmetry due to directional bonding. This bond order is also associated with the rotational symmetry that is broken upon crystallization, which is the key to a link between the two-state behavior and the phase diagram. The order parameter s is defined as the fraction of locally favored structures. The importance of the two-order parameter description was verified for model water by numerical simulations (38). Note that the density order parameter is conserved, but the bond order parameter is not since locally favored structures can be created and annihilated locally. This idea has been formalized by writing the free energy of water as that of a regular mixture of two components with very different degeneracy of states, under the additional equilibrium condition between the two components. This has produced a family of models that are often used to fit water’s equation of state with high precision (6, 16, 23⇓⇓⇓⇓–28, 30). Recent approaches go beyond the phenomenological use of a two-state equation of states and attempt to derive a two-state description starting from microscopic structural information (29, 39, 40). In our approach (27), we identify the two states according to the degree of translational order up to the second shell. By introducing a structural parameter that measures translational order (that we call ζ), we divide the population of water molecules into two collections of states: the S state, comprising highly ordered states, where there is a clear separation between the first and the second shell of nearest neighbors, and the ρ state, low-ordered states characterized by disordered arrangements of second-shell molecules, including configurations with shell interpenetration. We define s as the fraction of S state, which at any given T and P can be written, provided that there is little cooperativity in formation of locally favored structures, as (24, 25)

One compelling feature of our two-state model is that it can describe thermodynamic and dynamic anomalies (24⇓–26). In the case of dynamic anomalies, the predictions of two-state models are remarkably different from alternative explanations of dynamic anomalies. The major contender for the description of dynamic anomalies is based on glassy phenomenology, which is known as the fragile-to-strong transition (41⇓⇓⇓⇓–46).

In the case of our two-state model (24⇓–26), instead, the two different states have different activation energies, **1** can be approximated as (24⇓–26)**5**, and expanding to second order in β (high-T expansion), we get**5** predicts a full strong-to-strong transition, from activation energy

We note that very recently a different type of two-state model has been proposed, in which the ρ state behaves like a fragile liquid (47). We do not follow this route in the present article and instead give evidence that a pure ρ state behaves as a strong liquid in our model.

### Generalized SW Model.

The SW potential is composed of the sum of a pairwise term *Materials and Methods* for the definition of these terms):

By tuning λ one can continuously interpolate between the behavior of water-like materials and the behavior of simple fluids. To demonstrate this, in Fig. 1 we plot the melting line of the stable crystalline phase [the diamond cubic (dc) crystal] as a function of λ. The dc crystal is the only phase with a negatively sloped coexistence line, and the slope increases with increasing λ. Eventually at high λ the slope at

In *SI Appendix* we plot the full phase diagram of the model, extending the results of ref. 22 to negative pressures. Negative pressures are of great interest for at least two important reasons: (*i*) They stabilize clathrate lattices, which are crystalline structures with voids that can accommodate guest molecules and are studied for energy storage, carbon dioxide sequestration, separation, and natural gas storage (48⇓⇓–51); and (*ii*) contrasting theories of the thermodynamic anomalies (in particular for the case of water) can be tested in the negative pressure region, both numerically and experimentally (52⇓⇓–55). In *SI Appendix* we show that, at negative pressure, the body-centered cubic (BCC) phase is stable at lower λ and the Si34 phase is stable at higher λ. In the following sections we focus extensively on the line of liquid stability at negative pressures (the so-called spinodal). As a preliminary calculation, in *SI Appendix* we have mapped the location of the critical point (from which the spinodal emanates) for a large range of values of λ and reveal that increasing the tetrahedral parameter λ results in a lowering of both the critical temperature and pressure. As we will see later, this gives rise to a retracing spinodal (56) at low values of λ, when the spinodal line meets the line of density maxima.

### Thermodynamic Anomalies.

We have run extensive computer simulations to map the specific volume and compressibility anomalies in the **1**–**3**, which allows us to obtain the two-state model parameters *Left* column) and compressibility minima (*Right* column) anomalies for *Top* row to *Bottom* row). All anomalies shift to higher temperature with increasing λ, while also becoming more pronounced. The two-state model (solid lines) provides an excellent description of the anomalous behavior.

In Fig. 3 we plot the two-state model parameters obtained by fitting the thermodynamic anomalies of Fig. 2. Fig. 3*A* shows the increase of the fraction of the S state with decreasing T and for different values of λ. As λ is decreased from *B*–*D* we plot the variation with λ of the parameters *D* is less conclusive, but its decrease at high values of λ is in agreement with the change of the slope of the melting line at high λ displayed in Fig. 1. The increase in structural order in the ρ state with increasing λ, which is seen in the λ dependence of g, may be responsible for the decrease in

### Dynamic Anomalies.

The assumption about the two-state nature of water poses strong constraints on the nature of dynamic anomalies. As explained in *Two-State Model*, our two-state model predicts a strong-to-strong transition, contrary to the fragile-to-strong transition predicted by scenarios based on the glass-transition phenomenology. Here we emphasize that the strong-to-strong transition is the Arrhenius-to-Arrhenius transition and is independent from the glass transition. This is evident from the fact that the transition takes place far above the glass-transition point (*B* and *C*). Note that smaller λ means weaker directional bonds, resulting in the smaller energy difference between ρ and S states as well as the weaker constraint on particle configuration for the ρ state, which leads to the large degeneracy of the ρ state.

To study the dynamic behavior as a function of tetrahedrality, we run molecular dynamic simulations covering almost all of the accessible region of the

Fig. 4*A* shows the T dependence of the diffusion coefficient for selected values of λ. Our range goes from very high T [*Two-State Model*), which predicts the opposite transition, from super-Arrhenius to Arrhenius (i.e., the fragile-to-strong transition). The observed behavior is instead fully compatible with the two-state scenario for dynamic anomalies. Fig. 3*A*, *Inset* shows the amount of S state in the same range of **7** this transition can be fitted quadratically in *B* we plot the quadratic coefficient

From Eq. **7** we know that the quadratic term of the high-T expansion is *B* we superimpose the static term *B* clearly shows that the increase in

Fig. 4 *C* and *D* shows how the Arrhenius behavior of the ρ state changes with λ. We observe in particular that the activation energy

### Anomalies and Apparent Divergences.

We have seen that changing tetrahedrality is an effective tool to understand how both thermodynamic and dynamic anomalies emerge from ordinary liquid behavior. Here we show that altering λ can change the behavior of a water-like liquid at extreme conditions and affect its stability limit. We focus in particular on the liquid–gas spinodal line, or more precisely the line of liquid stability, below which the liquid becomes unstable to gas cavitation. We point out that simulation studies cannot access a true line of instability, as the cavitation of vapor is strongly system-size dependent. We nevertheless use the word “spinodal” to refer to this instability, as it is commonly used in the water literature (20, 21, 54). To determine this line we use two different procedures. First, we calculate the density dependence of the inverse of the isothermal compressibility at each temperature and obtain the spinodal points as the density where the inverse of isothermal compressibility sharply changes. The isothermal compressibility is computed in the

In Fig. 5 we summarize the loci of thermodynamic anomalies and liquid stability for *A*), 20.75, (Fig. 5*B*), 22.75 (Fig. 5*C*), and 23.15 (Fig. 5*D*). The most notable change occurs to the liquid spinodal (solid purple line) that emerges from the liquid–gas critical point (solid blue circle): While at high values of λ (Fig. 5 *C* and *D*) the spinodal displays usual monotonous behavior, and for low λ (Fig. 5 *A* and *B*) the spinodal intersects the line of density maxima (purple open circles) and retraces. Comparing these results with ref. 62, where the line of density maxima for silicon (

Next we focus on the apparent divergences in thermodynamic and dynamic properties of water at low T. In Fig. 5 we plot as open squares the estimated location of the apparent spinodal divergence **3**. We also plot as a solid red diamond in Fig. 5 the apparent dynamic divergence at *A*). Finally, we also plot the location of the Widom line (Fig. 5, solid green line) obtained from the two-state model (Eq. **4**), which is not a divergence, but the line along which *C* and *D*), where there is no retracing spinodal, and also the dynamic divergence does not occur at low values of λ, where the relaxation is more consistently fitted as Arrhenius (Fig. 4). These observations strongly hint at the fact that these divergences are only apparent. Their coincidence with the Widom (or Schottky) line indicates that the apparent divergences simply point to the loci of maximum change in the behavior of water, predicted as the Schottky anomaly of the two-state model: from the ρ state to the S state in the case of thermodynamic anomalies and from the high-T Arrhenius (strong) regime to the low-T Arrhenius (strong) regime. We thus believe that the two-state model can rationalize all of the observations of both thermodynamic and dynamic behavior across all values of λ, interpolating between simple liquid behavior found at low λ and the rich interplay of anomalies found at high λ.

## Discussion

In this article we have exploited the strategy of varying the tetrahedrality of the SW model to gain insights into the anomalous behavior of water and other tetrahedral materials in their liquid state as well as the phase behavior including all gas, liquid, and solid phases. The SW model has found widespread applicability in the study of thermodynamic anomalies of tetrahedral liquids, most notably silicon (62) and water (16). The first study to consider variations of λ as a means to change continuously the property of the materials was the seminal study of ref. 65, where the glass-forming ability was considered. Very recently the same idea was also applied to study the change in the anomalous properties of the liquid phase (20, 21).

In our work, we have computed the full phase diagram of the SW model as a function of the tetrahedral parameter λ. We have determined the phase diagram at negative pressures and also computed the λ dependence of the critical point. We then focused on liquid anomalies, with a special focus on the negative pressure region.

To rationalize the behavior of the anomalies we then applied a two-state model, fitting both the density and compressibility anomalies. The two-state model predicts an increase of the driving force toward the more ordered S state with increasing λ: Both the differences in energy

We then analyzed the behavior of dynamic anomalies, focusing on diffusion. We have shown that at small λ the dynamics are Arrhenius, while at large λ the dynamics cross to super-Arrhenius. The emergence of super-Arrhenius behavior from Arrhenius behavior, in coincidence with the increase in the fraction of the S state, is in line with the predictions of the two-state model, i.e., strong (Arrhenius)-to-strong (Arrhenius) transition, while it is at odds with interpretations based on the glass-transition singularity, i.e., fragile-to-strong transition. From a quadratic fit of the T dependence of the diffusion coefficient, we have also found that the activation energy difference

We also considered the location of the anomalies and apparent divergences in relation to the phase behavior. We found that by lowering λ the phase diagram changes to a retracing spinodal scenario, which occurs when the line of density maxima crosses the spinodal line. Increasing λ, the landscape changes from a retracing spinodal to a monotonous one, and the dynamic relaxation changes from Arrhenius to apparently super-Arrhenius. Despite these changes, all extrapolations based on singular behavior (spinodal divergences for thermodynamic anomalies and glass divergence for dynamic ones) always fall on top of the predicted two-state Schottky (or Widom) line. Starting from locally favored structures (the S state), the two-state model provides a unified description of water anomalies that is independent of singularities, while still being compatible with them. The similar scenario may be applied to water-like tetrahedral liquids such as silicon, germanium, and carbon. Here we also note that a unified scenario for other tetrahedral liquids, such as silica, was recently presented in ref. 66, and the connections to polyamorphism were discussed in ref. 67.

Finally, our study reveals that water is the material where tetrahedrality plays the bigger role: If tetrahedrality is weaker than water, the two-state feature becomes weaker, while if it is stronger than water, on the other hand, the volume difference between the two states becomes smaller, leading to a weaker density anomaly. On noting that the two-state feature is the origin of the flexibility of water properties or the large susceptibility of the properties to physical and chemical perturbations, our finding highlights the exceptional nature of water, which makes it so special compared with any other substances.

## Materials and Methods

### SW Potential.

Here, the pairwise term

### Numerical Methods.

To compute solid–liquid and liquid–gas coexistence lines, we run Monte Carlo simulations in the isothermal–isobaric

To obtain liquid–gas critical points, we run Monte Carlo simulations in the grand canonical ensemble. We compute the distribution functions of the mixing-order parameter M (

To compute liquid–gas spinodal points, we follow two strategies. In the first strategy we compute isothermal compressibilities, dividing the simulation box into smaller boxes to evaluate the size dependence of the compressibility; we then define the spinodal points as where the inverse of the compressibility vanishes. At lower temperatures, we instead run simulations in the

## Acknowledgments

This work was partially supported by Grants-in-Aid for Specially Promoted Research (25000002) from the Japan Society of the Promotion of Science. J.R. acknowledges support from the European Research Council Grant DLV-759187 and the Royal Society University Research Fellowship.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: tanaka{at}iis.u-tokyo.ac.jp or john.russo{at}bristol.ac.uk.

Author contributions: J.R. and H.T. designed research; J.R. and K.A. performed research; J.R., K.A., and H.T. analyzed data; and J.R., K.A., and H.T. wrote the paper.

The authors declare no conflict of interest.

↵*Tanaka H, Shi R, Russo J, The microscopic structural origin of water’s anomalies. Meeting of the American Physical Society, March 5–9, 2018, Los Angeles, CA, K57.00007 (abstr).

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1722339115/-/DCSupplemental.

- Copyright © 2018 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

## References

- ↵
- Debenedetti PG

- ↵
- ↵
- ↵
- ↵
- ↵
- Tanaka H

- ↵
- Stillinger FH,
- Weber TA

- ↵
- Molinero V,
- Moore EB

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Sosso GC,
- Li T,
- Donadio D,
- Tribello GA,
- Michaelides A

- ↵
- Pipolo S, et al.

- ↵
- ↵
- Sengupta S,
- Vasisht VV,
- Sastry S

- ↵
- Singh M,
- Dhabal D,
- Nguyen AH,
- Molinero V,
- Chakravarty C

- ↵
- Dhabal D, et al.

- ↵
- Angell CA,
- Kapko V

- ↵
- Dhabal D,
- Chakravarty C,
- Molinero V,
- Kashyap HK

- ↵
- Akahane K,
- Russo J,
- Tanaka H

- ↵
- ↵
- Tanaka H

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Holten V,
- Palmer JC,
- Poole PH,
- Debenedetti PG,
- Anisimov MA

- ↵
- Röntgen WC

- ↵
- Angell CA

- ↵
- ↵
- Ponyatovsky EG,
- Sinitsyn VV,
- Pozdnyakova TA

- ↵
- ↵
- ↵
- Tanaka H

- ↵
- ↵
- ↵
- ↵
- Xu L, et al.

- ↵
- Cerveny S,
- Mallamace F,
- Swenson J,
- Vogel M,
- Xu L

- ↵
- ↵
- Zhang Y, et al.

- ↵
- ↵
- Wang Z, et al.

- ↵
- Singh LP,
- Issenmann B,
- Caupin F

- ↵
- Florusse LJ, et al.

- ↵
- ↵
- Chatti I,
- Delahaye A,
- Fournaison L,
- Petitet JP

- ↵
- ↵
- Azouzi MEM,
- Ramboz C,
- Lenain JF,
- Caupin F

- ↵
- Pallares G, et al.

- ↵
- González MA,
- Valeriani C,
- Caupin F,
- Abascal JL

- ↵
- Holten V, et al.

- ↵
- ↵
- Romano F,
- Russo J,
- Tanaka H

- ↵
- Kennett MP,
- Chamon C,
- Cugliandolo LF

- ↵
- Rovere M,
- Hermann D,
- Binder K

- ↵
- ↵
- Prestipino S,
- Caccamo C,
- Costa D,
- Malescio G,
- Munaò G

- ↵
- ↵
- Stokely K,
- Mazza MG,
- Stanley HE,
- Franzese G

- ↵
- Rovigatti L,
- Bianco V,
- Tavares JM,
- Sciortino F

- ↵
- ↵
- Shi R,
- Tanaka H

- ↵
- Anisimov MA

- ↵
- ↵
- de Graaf J,
- Filion L,
- Marechal M,
- van Roij R,
- Dijkstra M

- ↵
- ↵
- Vega C,
- Sanz E,
- Abascal J,
- Noya E

- ↵
- Panagiotopoulos AZ

- ↵
- Tsypin M,
- Blöte H

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Applied Physical Sciences