{'fraction_non_alphanumeric': 0.035168028068636586, 'fraction_numerical': 0.013211994956416864, 'mean_word_length': 4.6127692307692305, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Water-like anomalies as a function of tetrahedrality

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.

L iquids do not possess long-range order but often have shortrange 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 Px . 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)(10)(11)(12)(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)(17)(18)(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 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.
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)(24)(25)(26)(27)(28)(29)(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 κT , which is proportional to volume fluctuations and displays a minimum at around T = 319 K, below which it shows a rapid increase. Similar anomalies are shown by the density ρ (which has a maximum at T = 277.15 K), and the specific heat Cp (which has a minimum around T = 308 K).
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)(33)(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)(24)(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)(24)(25)(26)(27)(28)30). Recent approaches go beyond the phenomenological use of a two-state equation of states and attempt to derive a twostate 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) where ∆E = Eρ − ES is the energy difference between the S and ρ states, ∆v = vS − vρ is their specific volume difference, and g is a measure of the degeneracy of the S state compared with the degeneracy of the ρ state (∆σ = k B ln g, where ∆σ is the entropy difference between the two states). The fraction of the S state controls the degree of anomalous behavior of the mixture. Following the notation of ref. 36, the specific volume is then given by v (T , P ) = a(P )T + b(P ) + s∆v , [2] and the isothermal compressibility by where the first two terms in each equation [a(P ), b(P ) and k (P ), n(P )] 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 s = 1/2 or the line of the Schottky anomaly (36) and can be written as Note 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 (J = 0 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)(25)(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)(42)(43)(44)(45)(46).
In the case of our two-state model (24)(25)(26), instead, the two different states have different activation energies, E a ρ and , and the diffusion process can be written as wheres is the fraction of the dynamic S state. It was assumed (36) thats = s, i.e., that the dynamic and static fractions of the S 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 g 1, corresponding to a much lower degeneracy of the S state compared with the ρ state, the expression Eq. 1 can be approximated as (24-26) Substituting this expression in Eq. 5, and expanding to second order in β (high-T expansion), we get [7] Eq. 5 predicts a full strong-to-strong transition, from activation energy E a ρ to E a S , instead of a fragile-to-strong transition. For small values of s, the crossover between the two strong behaviors can be fitted quadratically, with a coefficient that is proportional to ∆E a , i.e., the difference in the activation energy between the two states, S 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 U2 and a three-body interaction term U3 (see Materials and Methods 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 λ = 23.15, while for silicon the original parameter is λ = 21.0 (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 P = 0 becomes positive, when the diamond phase becomes more dense than the liquid. Fig.  1 shows that the change of slope occurs around λ ∼ 25. Thus, the V-shape feature of the phase diagram with ∂Tm/∂P |P=0 < 0 exists only in a limited range of λ; i.e., 16 < λ < 25. As shown later, this range roughly correspond to the region where we see water-like anomalies.
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)(49)(50)(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)(53)(54)(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 (T , P ) plane, for the values of λ = 19.55, 20.75, 22.75, and 23.15. For each value of λ we perform a multiparameter fit, where all simulation results are fitted against Eqs. 1-3, which allows us to obtain the two-state model parameters ∆E , ∆v , and g. In Fig. 2 we plot both the density maxima (Left column) and compressibility minima (Right column) anomalies for λ = 19.55, 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.
In Fig. 3 we plot the two-state model parameters obtained by fitting the thermodynamic anomalies of Fig. 2. Fig. 3A shows the increase of the fraction of the S state with decreasing T and for different values of λ. As λ is decreased from λ = 23.15 (the value of mW water), the fraction s decreases, and T W (the  Widom temperature, where s = 1/2) moves to lower temperatures. Also the variation of s with T becomes steeper at lower values of λ, meaning that the anomalies become more localized at lower T . To understand these results, in Fig. 3 B-D we plot the variation with λ of the parameters ∆E , g, and ∆v , respectively. ∆E , the energy difference between the S and ρ states, has the strongest dependence with λ, increasing by almost a factor of 5 going from λ = 19.55 to λ = 23.15. Similarly to ∆E , also g, the ratio between the degeneracies of the S state and the ρ state, increases rapidly with λ. This rapid increase in g 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 ∆E and g 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 ∆v in Fig. 3D 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 ∆v . 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 fragileto-strong transition predicted by scenarios based on the glasstransition phenomenology. Here we emphasize that the strongto-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 (∼2Tg ). 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 S -state rich liquid can be followed only up to s 0.5, as crystallization intervenes at high values of λ, while at low values of λ the increase of s is very weak in the observable T window, due to the small energy and large entropy difference (Fig. 3 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 T − λ parameter space and keeping P = 0. The simulations are limited at high T by the location of the liquid-gas spinodal (beyond which there is cavitation) and at low T either by crystallization (for λ 18 and λ 19) or by dynamical slowing down (for 18 λ 19). 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. 4A shows the T dependence of the diffusion coefficient for selected values of λ. Our range goes from very high T [Tmax = 0.2, which in mW units (57) corresponds to ∼Tmax = 622 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 lowtemperature behavior, we note that deviations from the Arrhenius behavior appear only at high values of λ, while for λ 21 the relaxation appears to be Arrhenius down to the lowest temperatures. For λ 21 the behavior changes from Arrhenius to super-Arrhenius with lowering T . 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. 3A, Inset shows the amount of S state in the same range of 1/T where deviations from Arrhenius behavior appear. For λ 21 the fraction s is negligible, and thus we expect the dynamics to display the strong (Arrhenius) behavior of the ρ state. For λ 21, instead, the fraction of s 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 S state. According to Eq. 7 this transition can be fitted quadratically in β = 1/k B T , and in Fig. 4B we plot the quadratic coefficient A2 (black circles) as a function of λ. The value of A2 confirms that the quadratic term is negligible for λ 21 and increases considerably at higher λ. The connection of fragile behavior at λ 21 with the increase of S 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.
From Eq. 7 we know that the quadratic term of the high-T expansion is A2 = ∆E a g(∆E − P ∆v ), where we can distinguish a dynamical term ∆E a , which is the difference in the activation energy between the S and ρ states, and a static term g∆E , where, without loss of generality, we used the fact that we are working at P = 0. In Fig. 4B we superimpose the static term g∆E (red squares), showing that it has a much weaker λ dependence than the quadratic term A2. This implies that also the dynamic term ∆E a is a strongly increasing function of λ (note that g is a constant). So the effect of tetrahedrality is to increase not only the energy (∆E ) and entropy difference (g) between the S and ρ states, but also the difference in their activation energies, ∆E a . The comparison of A2 and g∆E in Fig. 4B clearly shows that the increase in ∆E a (i.e., the latter) is the main cause of the non-Arrhenius behavior. The λ dependences of ∆E and ∆Ea explain why static and dynamic anomalies emerge from ordinary fluid behavior at high λ, respectively. Fig. 4 C and D shows how the Arrhenius behavior of the ρ state changes with λ. We observe in particular that the activation energy E a ρ has a minimum around 18 λ 19, which explains why the diffusion constant has a maximum in this region. This minimum in E a ρ may be a consequence of the competition between density and bond orderings (6): Small λ (λ ≤ 18) leads to a higher density and weaker directional bonds, whereas large λ (λ ≥ 18) 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 NVT ensemble via block analysis (59), where the distribution of the densityorder parameter is computed in blocks of different sizes. In the second procedure, which we use at lower T , we run extensive NVT simulations at size N = 512 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 ∂P ∂ρ |T = 0 and d 2 P d ρ 2 > 0. 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 λ = 19.55 (Fig. 5A), 20.75, (Fig. 5B ), 22.75 (Fig. 5C ), and 23.15 (Fig. 5D). 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 (λ = 21) was shown to just miss the spinodal line, we can estimate that the reentrant behavior of the spinodal starts at ∼ λ < 21. 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.
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 Tsp(P ), obtained from fitting the increase of the isothermal compressibility with the relation κT (T , P ) = k (P )T 2 + n(P ) + K (P )(T − Tsp(P )) −γ , [9] where the first two terms are the background behavior of the compressibility and whose coefficients are the same as the ones used in Eq. 3. We also plot as a solid red diamond in Fig. 5 the apparent dynamic divergence at P = 0, obtained by a Vogel-Fulcher-Tammann (VFT) fit of the diffusion data (Fig. 4A). 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 s = 1/2. 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 P = 0 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 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 ∆E and in degeneracy g 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 anoma-lies emerge continuously by increasing λ, moving progressively toward higher T and P .
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 twostate 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 ∆E a ρ 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 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 U 2 models a steep repulsion at short distances and a short-range attraction, while the three-body interaction term U 3 is a directional repulsive interaction which promotes tetrahedral angles between triplets of particles, The parameters for the models in this work are A = 7.049556277, B = 0.6022245584, p = 4, q = 0, cos θ 0 = −1/3, γ = 1.2, and a = 1.8. 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 NPT ensembles. The size and shape of the simulation box can fluctuate to allow crys-talline phases to change their structures (68,69). A volume-change attempt occurs every N translation attempts. The number of particles in the box is N = 1,024. 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 M (M = ρ + mu; ρ is density, u is internal energy per particle, and m 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 NVT ensemble and construct the whole equation of state (60,61), detecting the spinodal points as the points where ∂P ∂ρ | T = 0. To equilibrate simulations in the unstable region, we reduced the number of particles to N = 512. Both techniques gave similar results in the region of overlap. To obtain lines of density maxima and compressibility minima, we run NPT Monte Carlo simulations and compute averages and fluctuations of densities.