Water-like anomalies as a function of tetrahedrality.

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 [Formula: see text] We uncover a unique transition to a reentrant spinodal line at low values of [Formula: see text], 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.

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 , 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 (see, e.g., refs. 1–5 for water anomalies). Furthermore, all these liquids commonly have V-shaped – solid–liquid phase diagrams, in which the melting point has a minimum at a positive pressure . It was argued (6) that there is a deep link between the shape of the phase diagram and these anomalous thermodynamic and kinetic behaviors, as a consequence of local tetrahedral ordering. However, it has remained elusive how the degree of tetrahedrality controls the shape of the phase diagram and the anomalies. To address this problem, we need a model where we can control tetrahedrality in a systematic manner.

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 (, , ) 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 , which is proportional to volume fluctuations and displays a minimum at around K, below which it shows a rapid increase. Similar anomalies are shown by the density (which has a maximum at ), and the specific heat (which has a minimum around ).

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 , 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 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 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 as the fraction of state, which at any given and can be written, provided that there is little cooperativity in formation of locally favored structures, as (24, 25)where is the energy difference between the and states, is their specific volume difference, and is a measure of the degeneracy of the state compared with the degeneracy of the state (, where is the entropy difference between the two states). The fraction of the state controls the degree of anomalous behavior of the mixture. Following the notation of ref. 36, the specific volume is then given byand the isothermal compressibility bywhere the first two terms in each equation [, and , ] represent the background behavior and are obtained by fitting the specific volume and the compressibility far from the anomalous region. In the framework of the two-state model, the Widom line is nothing but the equimolar line or the line of the Schottky anomaly (36) and can be written asNote that two-state model predictions can accommodate a liquid–liquid critical point through a positive free enthalpy of the mixing term, but the Schottky anomalies arise whether this term is present or not. In this work we set the enthalpy of mixing to zero ( in the notation of refs. 24, 25, 28, and 36), as it produces the best results, also in line with what was observed for the mW model in ref. 16.

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, and (with ), and the diffusion process can be written aswhere is the fraction of the dynamic state. It was assumed (36) that , i.e., that the dynamic and static fractions of the states coincide. We adopt this assumption here. However, we note that for an accurate description of dynamic anomalies a hierarchical two-state model has to be considered.* In the case of , corresponding to a much lower degeneracy of the state compared with the state, the expression Eq. can be approximated as (24–26)Substituting this expression in Eq. , and expanding to second order in (high- expansion), we getEq. predicts a full strong-to-strong transition, from activation energy to , instead of a fragile-to-strong transition. For small values of , the crossover between the two strong behaviors can be fitted quadratically, with a coefficient that is proportional to , i.e., the difference in the activation energy between the two states, and .

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 and a three-body interaction term (see for the definition of these terms):Therefore, is the only parameter which differentiates the models. is a dimensionless parameter controlling the relative strength between pairwise and three-body terms. For the mW model of water (8) the value is , while for silicon the original parameter is (7). Furthermore, germanium and carbon are described by = 20.0 and 26.2, respectively.

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 becomes positive, when the diamond phase becomes more dense than the liquid. Fig. 1 shows that the change of slope occurs around . Thus, the V-shape feature of the phase diagram with exists only in a limited range of ; i.e., . As shown later, this range roughly correspond to the region where we see water-like anomalies.

Fig. 1.

Melting lines of the dc crystal at different values of . With increasing the temperature of the melting line increases, and the slope at goes from negative to positive; i.e., the crystal becomes more dense than the liquid.

In 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 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 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 plane, for the values of , 20.75, 22.75, and 23.15. For each value of we perform a multiparameter fit, where all simulation results are fitted against Eqs. –, which allows us to obtain the two-state model parameters , , and . In Fig. 2 we plot both the density maxima (Left column) and compressibility minima (Right column) anomalies for , 20.75, 22.75, and 23.15 (from 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.

Fig. 2.

Thermodynamic anomalies as a function of . From Top row to Bottom row, . (Left column) Specific volume () as a function of and for different . (Right column) Same as Left column, but for the isothermal compressibility . Symbols are results from simulations, while lines are fits according to the two-state model, Eqs. and .

In Fig. 3 we plot the two-state model parameters obtained by fitting the thermodynamic anomalies of Fig. 2. Fig. 3 shows the increase of the fraction of the state with decreasing and for different values of . As is decreased from (the value of mW water), the fraction decreases, and (the Widom temperature, where ) moves to lower temperatures. Also the variation of with becomes steeper at lower values of , meaning that the anomalies become more localized at lower . To understand these results, in Fig. 3 we plot the variation with of the parameters , , and , respectively. , the energy difference between the and states, has the strongest dependence with , increasing by almost a factor of going from to . Similarly to , also , the ratio between the degeneracies of the state and the state, increases rapidly with . This rapid increase in is more likely due to a decrease in the degeneracy of the state: As is increased, the liquid becomes progressively more ordered. Taken together, the increase of both and at high causes the emergence of anomalous behavior at higher temperatures and can be understood as an increase in the tetrahedral ordering of the fluid with (which controls the strength of the three-body interaction). They are also responsible for the ease of crystallization of the systems at high and the high glass-forming ability at low . At lower the thermodynamic driving force to form locally favored structures decreases, as the energy gain strongly decreases and the entropy loss also increases. The behavior of in Fig. 3 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 , may be responsible for the decrease in .

Fig. 3.

Two-state model analysis. () Fraction of the state as a function of for and different values of . Inset displays as a function of inverse temperature. , , and depict the dependence of two-state model parameters , , and , respectively.

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 (). Thus, the term “strong” liquid simply means a simple liquid obeying an Arrhenius law in this context. On a practical level, the transition from a strong -state rich liquid to a strong -state rich liquid can be followed only up to , as crystallization intervenes at high values of , while at low values of the increase of is very weak in the observable window, due to the small energy and large entropy difference (Fig. 3 ). Note that smaller means weaker directional bonds, resulting in the smaller energy difference between and 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 parameter space and keeping . The simulations are limited at high by the location of the liquid–gas spinodal (beyond which there is cavitation) and at low either by crystallization (for and ) or by dynamical slowing down (for ). Simulations are equilibrated in two steps, with isobaric–isothermal Monte Carlo first and isothermal molecular dynamics second. After equilibration, simulations are conducted in the microcanonical ensemble.

Fig. 4 shows the dependence of the diffusion coefficient for selected values of . Our range goes from very high [, which in mW units (57) corresponds to ∼ K] down to the homogeneous nucleation temperature. We find that if we include large temperatures, the diffusion coefficient displays sub-Arrhenius behavior, but if we limit the fits to low temperatures, we recover Arrhenius behavior. Interestingly this is the same behavior observed in lattice models of 2D doped antiferromagnets without quenched disorder (58). Focusing on the low-temperature behavior, we note that deviations from the Arrhenius behavior appear only at high values of , while for the relaxation appears to be Arrhenius down to the lowest temperatures. For the behavior changes from Arrhenius to super-Arrhenius with lowering . We already note that this is in contradiction to the glass-transition scenario (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, Inset shows the amount of state in the same range of where deviations from Arrhenius behavior appear. For the fraction is negligible, and thus we expect the dynamics to display the strong (Arrhenius) behavior of the state. For , instead, the fraction of increases considerably, and we thus expect the system to display a transition from the strong (Arrhenius) behavior of the state to the strong (Arrhenius) behavior of the state. According to Eq. this transition can be fitted quadratically in , and in Fig. 4 we plot the quadratic coefficient (black circles) as a function of . The value of confirms that the quadratic term is negligible for and increases considerably at higher . The connection of fragile behavior at with the increase of state and the observation of super-Arrhenius behavior at high emerging continuously from a pure Arrhenius relaxation at low strongly support the two-state interpretation of the dynamic anomaly.

Fig. 4.

Dynamics and anomalous behavior. () Arrhenius plot for the inverse of the diffusion coefficient vs. the inverse temperature . Symbols are simulation results, while the dashed lines are the Arrhenius fit to the high- behavior. () Coefficient of the quadratic term in Eq. (black circles) and the product of (red squares). () The activation energy of the state, , as a function of . () Inverse of the diffusion prefactor, , as a function of .

From Eq. we know that the quadratic term of the high- expansion is , where we can distinguish a dynamical term , which is the difference in the activation energy between the and states, and a static term , where, without loss of generality, we used the fact that we are working at . In Fig. 4 we superimpose the static term (red squares), showing that it has a much weaker dependence than the quadratic term . This implies that also the dynamic term is a strongly increasing function of (note that is a constant). So the effect of tetrahedrality is to increase not only the energy () and entropy difference () between the and states, but also the difference in their activation energies, . The comparison of and in Fig. 4 clearly shows that the increase in (i.e., the latter) is the main cause of the non-Arrhenius behavior. The dependences of and explain why static and dynamic anomalies emerge from ordinary fluid behavior at high , respectively.

Fig. 4 shows how the Arrhenius behavior of the state changes with . We observe in particular that the activation energy has a minimum around , which explains why the diffusion constant has a maximum in this region. This minimum in may be a consequence of the competition between density and bond orderings (6): Small () leads to a higher density and weaker directional bonds, whereas large () leads to a lower density and stronger bonds. Note that both higher density and stronger bonds result in the higher 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 ensemble via block analysis (59), where the distribution of the density-order parameter is computed in blocks of different sizes. In the second procedure, which we use at lower , we run extensive simulations at size at different densities and equilibrate the equation of state in the unstable region (60, 61), where the spinodal point can be obtained from the condition and . The results of the two methods match in the region of intermediate temperatures. To get lines of density maxima, we compute the isobaric temperature dependence of densities and obtain temperatures of density maxima by polynomial fitting. The same procedure is applied to compute the line of compressibility minima.

In Fig. 5 we summarize the loci of thermodynamic anomalies and liquid stability for (Fig. 5), , (Fig. 5), (Fig. 5), and (Fig. 5). 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 ) the spinodal displays usual monotonous behavior, and for low (Fig. 5 ) 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 () was shown to just miss the spinodal line, we can estimate that the reentrant behavior of the spinodal starts at . This shows that a transition to a reentrant spinodal (56), a result which was predicted in terms of a mean-field cell model (63), can be achieved with a microscopic model of water. Very recently, in patchy particles colloidal models, the authors of ref. 64 also observed a retracing spinodal, which in their case extended to positive pressures.

Fig. 5.

Stability limits and anomalies. A–D differ in their value of : (), (), (), and (). In A–D, the solid blue circle is the liquid–gas critical point; the solid blue line is the liquid–gas coexistence line; the solid purple line is the liquid spinodal (or better stability limit of the liquid phase); the solid orange line is the line of compressibility minima; the open magenta circles are the line of density maxima; the open squares are the location of the spinodal line extrapolated from the apparent divergence of the compressibility, which may or may not be contiguous with the liquid spinodal; the solid red diamond is the ideal glass-transition temperature extrapolated by VFT fit of the diffusivity; the solid blue triangle is the homogeneous nucleation point at ; the solid green line is the Schottky (or Widom) line predicted by the two-state model; and the dashed black lines represent the liquid–solid coexistence line. A and B, Insets show an enlargement of the minimum in the vicinity of the spinodal line.

Next we focus on the apparent divergences in thermodynamic and dynamic properties of water at low . In Fig. 5 we plot as open squares the estimated location of the apparent spinodal divergence , obtained from fitting the increase of the isothermal compressibility with the relationwhere the first two terms are the background behavior of the compressibility and whose coefficients are the same as the ones used in Eq. . We also plot as a solid red diamond in Fig. 5 the apparent dynamic divergence at , obtained by a Vogel–Fulcher–Tammann (VFT) fit of the diffusion data (Fig. 4). Finally, we also plot the location of the Widom line (Fig. 5, solid green line) obtained from the two-state model (Eq. ), which is not a divergence, but the line along which . The spinodal line, the glass-transition point, and the Widom line cannot be accessed in equilibrium, as they lie below the homogeneous nucleation line, whose point is plotted as a solid blue triangle in Fig. 5. It is important to observe that the spinodal divergence and the dynamical divergence fall on top of each other within the errors and on top of the Widom line. We have shown that the spinodal divergence should not occur at high (Fig. 5 ), 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 state in the case of thermodynamic anomalies and from the high- Arrhenius (strong) regime to the low- 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 state with increasing : Both the differences in energy and in degeneracy increase considerably with , as a consequence of a tendency for the liquid to become more ordered locally as the strength of the tetrahedral interaction becomes stronger. This explains how the anomalies emerge continuously by increasing , moving progressively toward higher and .

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 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 dependence of the diffusion coefficient, we have also found that the activation energy difference is a strongly increasing function of .

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 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 models a steep repulsion at short distances and a short-range attraction,while the three-body interaction term is a directional repulsive interaction which promotes tetrahedral angles between triplets of particles,The parameters for the models in this work are , , , , , , and . The parameter sets the energy scale and the length scale. They correspond to the depth of the two-body interaction potential and the particle diameter, respectively, and are determined by materials for which the model is used. We use internal units where and are the units of energy and length, respectively.

Numerical Methods.

To compute solid–liquid and liquid–gas coexistence lines, we run Monte Carlo simulations in the isothermal–isobaric ensembles. The size and shape of the simulation box can fluctuate to allow crystalline phases to change their structures (68, 69). A volume-change attempt occurs every translation attempts. The number of particles in the box is . We perform Gibbs–Duhem integration (70) and Hamiltonian Gibbs–Duhem integration (71) to obtain coexistence lines along the pressure axis and along the axis, respectively. Triple lines are computed in the same way as in ref. 22.

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 (; is density, is internal energy per particle, and is mixing parameter) and use histogram reweighting methods (72) to fit them into the Ising universal curve (73). Liquid–gas coexistence lines are instead computed by locating a coexistence point close to the critical point with successive umbrella sampling simulations and then running Gibbs–Duhem integration to trace the coexistence line at lower temperatures.

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 ensemble and construct the whole equation of state (60, 61), detecting the spinodal points as the points where . To equilibrate simulations in the unstable region, we reduced the number of particles to . Both techniques gave similar results in the region of overlap. To obtain lines of density maxima and compressibility minima, we run Monte Carlo simulations and compute averages and fluctuations of densities.