Analytical Modeling of Thermodynamic and Transport Anomalies of Water

The structures and properties of biomolecules like proteins, nucleic acids, and membranes depend on water. Water is also very important in industry. Overall, water is an unusual substance with more than 70 anomalous properties. The understanding of water is advancing significantly due to the theoretical and computational modeling. There are different kinds of models, models with fine-scale properties and increasing structural detail with increasing computational expense, and simple models, which focus on global properties of water like thermodynamics, phase diagram and are less computationally expensive. Simplified models give a better understanding of water in ways that complement more complex models. Here, we review analytical modelling of properties of water on different levels, the two- and three-dimensional Mercedes– Benz (MB) models of water and experimental water.


Introduction
The Earth is a watery place by water being the most important fluid in nature for life and for humans in the industry. [1][2][3][4] About 71 percent of the Earth's surface is water-covered, and the oceans hold about 96.5 percent of all Earth's water. Water also exists in the air as water vapor, in rivers and lakes, in icecaps and glaciers, in the ground as soil moisture etc. Water controls the planets geochemical cycles; is a dominant driver of biomolecules, drug interactions, and biological actions; and central to green chemistry and many industrial processes. 5,6 Water is essential for our bodies. Every major system in our body depends on water to function since approximately two thirds of your body is water. Due to all these facts the structure and thermodynamics of water and aqueous solutions is of great importance for all sciences, especially chemistry and biology. Water exhibits many anomalous properties that affect life at a larger scale. Many animals benefit from the large latent heat of water to cool them down through sweating. The large heat capacity of water prevents local temperature fluctuations, facilitating thermal regulation of organisms. The density anomaly and lower ice (hexagonal ice I h ) density have a huge effect on surviving of organisms in frozen seas and lakes. Water is an almost universal solvent. 7 Nearly all known chemical substances will dissolve in water at least to a small extent. In comparison to other liquids, it has the most puzzling behavior. 7,8 It is said that water is an anomalous liquid. Anomalous liquids are liquids that exhibit unexpected behavior upon variations of the thermodynamic conditions in comparison to normal (argon-like) liquids. Water is the classic example of such anomalous liquids. Water's density maximum at 4 °C, the lower density of the solid phase compared to the liquid phase, high and nearly constant heat capacity in the liquid phase, negative expansion coefficient below the temperature of the density maximum, as well as high surface tension and viscosity are the most known examples of anomalous properties. If we continue, we have an anomalous increase in the compressibility and specific heat by cooling, unusually high boiling, freezing and critical points. The reason for water's complexity is due to its strong orientation-dependent hydrogen bonding and strong intermolecular associations.
An understanding of the hydrogen bonding is therefore crucial to understand the behavior and properties of water and aqueous solutions. Yet, despite extensive theory and simulations, the fact how water's properties come from its molecular structure remains poorly understood. Many models of varying complexity have been developed and analyzed to model water's extraordinary properties, Urbic: Analytical Modeling of Thermodynamic and Transport ...
for review see literature. [7][8][9][10][11][12][13][14][15][16][17][18][19] The rigid models are considered the simplest water models and rely on non-bonded interactions. The electrostatic interaction is modeled using Coulomb's law, and the dispersion and repulsion forces using the Lennard-Jones potential. Examples of such models are SPC (simple point-charge) 20 , TIP3P (transferable intermolecular potential with 3 points) 20 and TIP4P (transferable intermolecular potential with 4 points) 21 etc. In polarizable models we consider many-body energies which can be effectively accounted for by a single term representing classical many-body polarization. Several polarizable water models, with different degrees of sophistication, have been developed and used in molecular dynamics and Monte Carlo simulations of aqueous systems. 22 A key goal of the liquid-state statistical thermodynamics is to develop a quantitative theory for water and aqueous solutions. Theory and simulations have not yet been able to explain how water's molecular structure leads to its density, compressibility, expansion coefficient and heat capacity as functions of temperature and pressure, including its well-known anomalies. The properties of water can be determined with quantum-mechanical calculations. 22,23 These methods offer the highest degree of exactness, but a high computational cost of these approaches limits their use to small water systems, even though these insights allow the development and fine-tuning of simplified water models. [24][25][26] There have been two main approaches to modeling liquids. One approach is to perform computer simulations of atomically detailed models. Another way captures many properties of water and aqueous solutions by simpler models.
One of the simplest models for water is the socalled Mercedes-Benz (MB) model, 27 which is a 2-dimensional model that was originally proposed by Ben-Naim in 1971. 28,29 Each MB water particle is modeled as a disk that interacts with other particles through: (1) a Lennard-Jones (LJ) interaction and (2) an orientation-dependent hydrogen bonding interaction through three radial arms arranged as in the MB logo. Interest in simplified models is due to insights that are not obtainable from all-atom computer simulations. Simpler models are more flexible in providing insights and illuminating concepts, and they do not require big computer resources. The analytical models can also provide functional relationships for engineering applications and lead to improved models of greater computational efficiency. For the MB model, the NPT Monte Carlo simulations have shown that it predicts qualitatively the density anomaly, the minimum in the isothermal compressibility as a function of temperature, the large heat capacity, as well as the experimental trends for the thermodynamic properties of solvation of nonpolar solutes 27 and cold denaturation of proteins. 30 The MB model was also extensively studied with analytical methods like integral equation and thermodynamic perturbation theory [31][32][33][34][35][36] and statistical mechanic modeling [37][38][39] . Recently also phase diagram of liquid part and percolation curve of the model was calculated and reported. 40 The MB model has also been used to study systems with water molecules confined in partially quenched disordered matrix [41][42][43] and within small geometric spaces. 44,45 Nonequilibrium Monte Carlo and molecular dynamics simulations were used to study the effect of translational and rotational degrees of freedom on the structural and thermodynamic properties of this MB model. [46][47][48] By holding one of the temperatures constant and varying the other one, the effect of faster motion in the corresponding degrees of freedom on the properties of the simple water model was investigated. The situation where the rotational temperature exceeded the translational one is mimicking the effects of microwaves on the water model. A decrease of rotational temperature leads to the higher structural order while an increase causes the structure to be more Lennard-Jones fluid like. The 2D MB model was also extended to 3D by Dias et al. 49 and Bizjak et al. 50,51 Even though computer simulations play an important role in understanding the properties of liquids, they can be quite time consuming, even for simple two-dimensional 2D models. Due to this it is equally important to develop simplified, more analytical approaches. One such model is a statistical mechanical model, developed by Urbic and Dill 33 . The model is directly descendant from a treatment of Truskett and Dill, who developed a nearly analytical version of the 2D MB model. 52,53 In the model, each water particle interacts with its neighboring particle through a van der Waals interaction and an orientation-dependent interaction that models hydrogen bonds. Recently this theory was extended to 3D MB model 38 and later parametrized to describe properties of experimental water. 54 In this paper, we made review of analytical modeling for MB model of water, its properties in bulk which are starting point to develop the theory for solvation of polar and nonpolar solutes, important for example in self-association of surface-active compounds such as ionic liquids, 55 protein folding, etc. The outline of the paper is as follows. We present the 2D and 3D MB model in Sec. 2, and the details of the statistical mechanical methods are done in Sec. 3. In Sec. 4 we show and discuss the results and summarize everything in Sec. 5.

1. 2D MB Model
In 2D, the water particles are modelled as a two-dimensional disk with three bonding arms separated by an angle of 120°, which is fixed as in Mercedes-Benz logo (See Figure 1). 27 These arms mimic formation of hydrogen bonds. The interaction potential between particles i and j is a sum of a Lennard-Jones (LJ) and a hydrogen-bonding (HB) term (1) Where r ij is the distance between centers of particles i and j. X u i , X u j are the vectors representing the coordinates and the orientation of the particles i and j. The Lennard-Jones part has a standard form (2) s LJ and e LJ are the depth and the contact value of the LJ potential. The hydrogen bonding part is the sum of interactions between all pairs of the arms of different molecules (3) and is described by Gaussian function in distance and both angles (4) Here, e HB = -1 is a HB energy parameter and r HB = i is a characteristic length of HB, u u ij is the unit vector along r u ij and i u k is the unit vector of the k th arm of the particle I, and q i is the unit vector of the i th arm of the particle j. q i , q j are the orientations of the particle with respect to x axes. G(x) is the unnormalized Gaussian function (5) The strongest hydrogen bond occurs when an arm of one particle is co-linear with the arm of another particle and the two arms point in opposing directions. The LJ well-depth e LJ is 0.1 times the HB interaction energy |e HB | and the Lennard-Jones contact parameter s LJ is 0.7 r HB . The width of the Gaussian function for distances and angles (s = 0.085 r HB ) is small enough that a direct hydrogen bond is more favorable than a bifurcated one.

3D MB Model
In 3D, each water molecule is represented as a Lennard-Jones sphere (LJ) with four arms oriented tetrahedrally. 50 The angle between neighboring arms in a molecule is 109.47° (see Figure 2). Like in 2D, in 3D the inter-  Urbic: Analytical Modeling of Thermodynamic and Transport ... action potential between two water molecules is a sum of the Lennard-Jones potential and the hydrogen bond term (6) The Lennard-Jones part of the potential is the same as in 2D. The hydrogen bonding part of the interaction potential is (7) Where Ω u i . Ω u j are the orientational vectors of both particles and U kl describes the interaction between two HB arms of different molecules (8) Like in 2D, the strongest hydrogen bond occurs when an arm of one particle is colinear with the arm of another particle pointed towards each other. The model does not make a distinction between hydrogen bond donors and acceptors. Apart from the dimensionality, we want to keep the 3D MB model as similar as possible to the original 2D MB model. Hence, the parameters of our 3D model are the same as used in the 2D MB model calculations, except for the depth of the Lennard-Jones potential well e LJ . This change was made to maintain the same ratio between strength of the Lennard-Jones interaction and hydrogen bond interaction due to the different geometries; e LJ =1/35 e HB . These model parameters were not chosen or optimized to compare with experiments and can undoubtedly be improved for those purposes.

2D MB Model
In the theory, the system consists of N water molecules. 37 To keep track of the state of interaction of each possible hydrogen bonding arm of each water molecule we are using an underlying ice lattice as a bookkeeping tool. For the 2D water model, the underlying lattice is hexagonal (See Figure 3). We focus on a single water molecule in the hexagon and the relationship of that water to its clockwise neighbor. Figure 4 shows the three possible relationships: the test water can either form a hydrogen bond, a van der Waals contact, or no interaction at all. We compute the isothermal-isobaric statistical weights, Δ HB of the hydrogen-bonded molecules, Δ LJ of the van der Waals contacts, and Δ 0 of the unbonded population as functions of temperature, pressure, and interaction energies.
The hydrogen-bonded state. If the test water molecule points one of its three hydrogen bonding arms at an angle θ to within π/3 of the center of its clockwise neighbor water, it forms a hydrogen bond. The energy of interaction of the test water is (9) k is the angular spring constant that describes the weakening of the hydrogen bond as it becomes increasingly off-angle, and e HB and e LJ are the potential energy parameters. We regard this type of hydrogen bond as weak insofar as it is not cooperative with neighboring hydrogen bonds. We consider a more cooperative type of hydrogen bonding below. To compute the isothermal-isobaric partition func- tion, Δ HB , of this HB state, we integrate this Boltzmann factor over all the allowable angles and over all the allowable separations x and y of the test molecule relative to its clockwise neighbor, (10) Where b = 1/k B T is inverse temperature, c(T) is the 2D version of the kinetic energy contribution to the partition function and v HB is volume per molecule in this state. òò dxdy represents the volume over which the second molecule has translational freedom to form a hydrogen bond with the first water and is equal to the effective volume . The volume v HB of the hydrogen-bonded state is determined in the following way. First, we estimate an upper bound on the volume, from a simple geometric calculation. For the perfect hexagon crystal, representing low-pressure ice, the volume of the solid if the center of the hexagon is unoccupied is (11) Second, since liquid water is denser than ice, we estimate a lower bound on the volume using high-pressure ice, where another MB water occupies the center of each hexagonal cage 52,53 (12) Since the density of liquid water must lie between these limits, we estimate its volume as (13) Where x v = 1.01 is chosen empirically by fitting the density dependence vs. temperature. Using these definitions and performing the integration in Equation (10) gives (14) The van der Waals (vdW) state. Here, the test water molecule forms only a van der Waals contact with its clockwise neighboring water. The water molecule has an energy (15) The isothermal-isobaric partition function, Δ LJ of this state is given by integrating over angles and positions of the test particle relative to its clockwise neighbor as in case of the HB state (16) The integral òò dxfy represents the translation volume over which the second molecule forms a van der Waals contact with the first water and is equal to the effective volume v eff LJ = 0.104. The volume occupied by molecule in this state, v LJ , is volume of packed LJ disks (17) Integrating the partition function gives (18)  The non-interacting state. In this third possible state, the test water has no interaction with its clockwise neighbor (19) The same way as for the other two states, the isothermal-isobaric partition function is obtained by integrating over translational degrees of freedom (20) Urbic: Analytical Modeling of Thermodynamic and Transport ...
where v 0 is the volume available to the test molecule in this state and is calculated as the van der Waals gas approximation (21) Where v d = v s for MB water. Integrating of the partition function gives (22) These three expressions, Equations (14), (18) and (22), give the isobaric-isothermal ensemble Boltzmann weights of the three possible states of each water molecule. We assume a mean-field attractive energy, -Na/v, 52,53 among cages, where a is the van der Waals dispersion parameter (0.02, here) and v is the average molar volume, which we will define below. The partition function for a single full hexagon of 6 waters would be given by (23) Here, we treat hexagons a little differently instead. We define cooperative HB state or solid state that involves a higher degree of HB cooperativity than the hydrogen bonding that is just formed pairwise among nearest neighbor waters in the liquid state. So, the partition function for each hexagon will be given by (24) where δ = e -βe c is the Boltzmann factor for the cooperativity energy, e c , that applies only when 6 water molecules all connected into a full hexagonal cage. The terms on the right-side of this expression simply replace the statistical weight for each weakly hydrogen-bonded full hexagonal cage with the statistical weight for a cooperative strongly hydrogen-bonded hexagonal cage. Δ s is the Boltzmann factor for a cooperative HB or solid state. It differs from Δ HB only in that the former uses the hexagonal cage volume per molecule, v s , while the HB state uses the liquid water hydrogen bonding volume per molecule, v HB . Now we combine the Boltzmann factors for the individual water molecules to get the partition function Q for the whole system of N particles, (25) The factor N/6 accounts for the 3 possible interaction sites per water molecule and corrects for double counting the hydrogen bonds. We compute the populations of the states i = 1 (HB), 2 (LJ), 3(o), 4(solid) using (26) The chemical potential is given by (27) The molar volume is (28) and all the other thermodynamic properties below are obtained as described previously. 52,53 For all the model calculations, we used the parameters from potential function e HB , e LJ , r HB and s LJ . The parameter k = 10 was determined from the shape of the MB potential while e c = 0.03 was determined empirically.

2. 3D MB Model
Here, we will point out only the differences between the theory in 3D in comparison to 2D. 38 We consider a system of 3D MB model water molecules, modeled as three-dimensional spheres, and suppose that the structure of the liquid state of 3D model water is a perturbation from an underlying hexagonal (ice) lattice; (See Figure 5). Each molecule can be in one of the three possible orientational states like in 2D. These states are graphically presented in Fig. 6. The hydrogen-bonded state. If the test water molecule points one of its four hydrogen bonding arms at an angle θ to within π/3 of the center of its clockwise neighbor water, it forms a hydrogen bond. This is equivalent to about one fourth of the full solid angle. The energy of interaction of the test water is (29) k is the angular spring constant that describes the weakening of the hydrogen bond as it becomes increasingly off angle. To compute the isothermal-isobaric partition func-tion, Δ HB , of this HB state, we integrate this Boltzmann factor over all the allowable angles and over all the allowable separations of the test molecule relative to its clockwise neighbor, (30) Where c(T) is here the 3D version of the kinetic energy contribution to the partition function. òòò dxfydz represents the volume over which the second molecule has translational freedom to form a hydrogen bond with the first water and is equal to effective volume v eff HB = 0.242. The double integral òò dafy sums the orientations over which the test molecule has orientational freedom and is equal to 4π 2 . The volume v HB of the hydrogen-bonded state we determine similarly as for the 2D model. For the perfect hexagon crystal representing low-pressure ice, the volume of the solid is (31) We estimate volume v HB as perturbation of this state as (32) where x v = 1.12 is chosen empirically because density of the liquid state at room temperature is about 12% more dense than ice. Using these definitions and performing the integration in Equation (30) gives (33)

The van der Waals (vdW) state.
Here, the test water molecule forms only a van der Waals contact with its clockwise neighboring water. The water molecule has an energy (34) The isothermal-isobaric partition function, Δ LJ is (35) The triple integral òòò dxfydz represents the translation volume over which the second molecule forms a van der Waals contact with the first water and is equal to effective volume v eff LJ = 0.086. The integrals over angles are equal to 8π 2 . The volume v LJ of this state is approximated as a volume of cubic close-packed crystal where the closest molecules are at distance σ LJ √ 6 -2 we also tried other symmetries, but the results did not change much. Integrating of the partition function gives (37) The non-interacting state. In this third possible state, the test water has no interaction with its clockwise neighbor and the isothermal-isobaric partition function is equal to (38) Here, we also assume a mean-field attractive energy, -Na/v, 52,53 among cages, where a is the van der Waals dispersion parameter (0.045, here). The partition function for a single full hexagon of 6 waters and other properties are calculated in the same way as in 2D. For all the model calculations, we used the parameters from potential function e HB , e LJ , r HB and s LJ . The parameter k = 80 was determined from the shape of the 3D MB potential while e c = 0.18 was determined empirically.

The Real Water -CageWater
Here we made slight modification in comparison with 3D MB. 54 Two water molecules can interact through a hydrogen bond (which depends on their relative orientations), interact through a contact (which is orientation independent and occurs when they are close in space and no HB is present), or be noninteracting (when they are far apart, as in van der Waals gas). Hydrogen bonds are further parsed into two types: an HB can occur between 2 adjacent waters that have no higher order structure or can occur within a 12-water hexagonal unit cell (cage) forming 15 HBs. Parameters needed for calculations were obtained by getting good agreement with tempera-Urbic: Analytical Modeling of Thermodynamic and Transport ... ture dependence of density at normal pressure and of boiling point position and are presented in Table 1.

Results
We present our results below in dimensionless or reduced units for both MB models, normalized to the strength of the optimal hydrogen bond e HB and hydrogen bond separation r HB for 2D MB model and for 3D MB model).

1. 2D MB Model
Analytical theory has additional approximations compared to computer simulations, which is why we first  checked the quality of the predictions of the analytical theory. We calculated the temperature dependence of the density (ρ * ), heat capacity (c p * ), thermal expansion coefficient (α * ), and isothermal compressibility (κ * ) for different pressures. For a 2D MB model it was previously shown that the Mercedes-Benz water qualitatively correctly reproduces the anomalies of water 27,31,32 for these quantities. In Fig. 7 a comparison of predictions of the present theory (lines) for the density, the thermal expansion coefficient, the isothermal compressibility, and the heat capacity vs temperature to NPT Monte Carlo simulations (symbols) of the 2D MB model with the same parameters is shown. The calculations of the theory were performed at reduced pressure of 0.08, 0.12, 0.19 and 0.32. The theory is in good general agreement with the simulations, including the density maximum (minima in molar volume). The thermal expansion coefficient is negative at low temperatures, which is consistent with computer simulations and with experiments for water. The Monte Carlo simulations of MB water do not show an experimentally observed minimum in the isothermal compressibility versus temperature. On the other hand, the present theory predicts a minimum for higher pressures. At low temperatures, our present model shows a drop in heat capacity as the temperature is reduced. Being satisfied with the prediction of the model, we calculated non crystalline part of the phase diagram, shown in Fig. 8. The 2D MB model exhibits two critical points: the liquid-gas critical point (C 1 ) at temperature T * = 0.118 and pressure p * = 0.00035 which is slightly lower than obtained from computer simulations, 40 and the liquid-liquid critical point (C 2 ) at temperature 0.0212 and pressure 0.42. There exists also a region of pressures between both critical points where we have only one fluid phase, at higher pressures we have two liquid phases, and at lower pressures the liquid and the gas phases.  We have also developed the theory for dynamical properties. The diffusion processes occur in fluid or gas whenever a property is transported in a manner resembling a random walk. If we assume that the water molecules are doing random walk, we can say, for our 2D molecules in each state, that the diffusion is proportional to the step size, l, and a step frequency, n. The step frequency is proportional to the Boltzmann factor for the change of energy from bonded to free state. This means that the energy of interaction is negative of bonding energy of molecule. We assumed that the step size is equal to the characteristic length of interaction in each state (l HB = l s = r HB for HB and s state, l LJ = s LJ for LJ state and for 0 state for 0 state the average distance between molecules in this state l 0 = √ v 0 ). For our model, we calculated the diffusion constant as a sum of all states of individual contributions (39) Where D i = λ I n i for HB, s, LJ and 0 state of water. The step frequency is equal to Boltzmann factor of negative average bonding energy (u i ) (40) To model viscosity, h, we start from Stokes-Einstein relation between viscosity and diffusion coefficient, D, We can express viscosity from this equation as (42) We see that we can calculate viscosity from diffusion coefficient, temperature and diameter, d, of particle. In our case, we use averaged particle diameter which we calculated as a sum over all possible states of water (43) For water molecules in states HB, LJ and 0 we used diameter of molecule equal to r HB while for state s waters form hexagons and we use the diameter of water in hexagon state as equal to 2r HB . Next, we calculated the speed of sound c s as (44) and thermal conductivity l using modified Eyring's formula as (45) and thermal diffusivity l d as Figure 10: Pressure dependence of the diffusion (D * ), viscosity (h * ), thermal conductivity (l * ) and thermal diffusivity (l d * ) for 2D MB water for different temperatures calculated by the theory. (46) In Fig. 9 and 10, we have plotted temperature and pressure dependence of the dynamical properties (diffusivity, viscosity, thermal conductivity, and thermal diffusivity). All the quantities have similar anomalous nonmonotonic behavior as for experimental water.

2. 3D MB Model
As for 2D MB we first checked the quality of the predictions of the analytical theory also for 3D MB. We also calculated the temperature dependence of the molar volume, heat capacity, isothermal compressibility, and thermal expansion coefficient for different pressures. For a 3D MB model, it was previously shown that the Mercedes-Benz water qualitatively correctly reproduces the anomalies of water 49,50,51 for these quantities. In Fig. 11, a comparison of predictions of the present theory (lines) for the molar volume (V * /N), heat capacity (c p * ), thermal expansion coefficient (α * ), and isothermal compressibility (κ * ), vs temperature to NPT Monte Carlo simulations (symbols) of the 3D MB model with the same parameters is shown. The calculations of the theory were performed at reduced pressure of 0.12 and 0.19. The theory is in good general agreement with the simulations, including the minima in molar volume (density maximum). The thermal expansion coefficient is negative at low temper-atures, which is consistent with computer simulations and with experiments for water. The Monte Carlo simulations of MB water do not show an experimentally observed minimum in the isothermal compressibility versus temperature. On the other hand, the present theory predicts a minimum for higher pressures. At low temperatures, our present model shows a drop in heat capacity as the temperature is reduced. Being satisfied with the prediction of the model, we continued our research by calculating the density of 3D MB water as a function of temperature along isobars (up to p* = 0.25) and determine critical points of the model. Results are shown in Fig. 12. In this pressure range, upon increase of temperature density increases, reaches a maximum, and then decreases. We determined non crystalline part of the phase diagram, shown in Fig. 13. The 3D MB model exhibits two critical points: the liquid-gas critical point (C 1 ) at temperature T * = 0.117 and pressure p * = 0.0115, and the liquid-liquid critical point (C 2 ) at temperature 0.0779 and pressure 0.167. There exists also a region of pressures between both critical points where we have only one fluid phase, at higher pressures we have two liquid phases, and at lower pressures the liquid and the gas phases.

3. The Real Water -CageWater
Here, we compare the measured properties over water's liquid range to those predicted by parametrization for Figure 11: Temperature dependence of the molar volume (V * /N), heat capacity (c p *), thermal expansion coefficient (α*), and isothermal compressibility (κ*), for 3D MB water for pressures p* = 0.19 (orange) and p* = 0.12 (blue). Results from the theory are plotted with lines and from the computer simulations by symbols. 50,51 Urbic: Analytical Modeling of Thermodynamic and Transport ... experimental data, called CageWater. 54 In Fig. 14a we have plotted experimental data 8 and data by best practices water simulation models TIP4P/2005, 21 TIP3P, 20 SPC, 20 and mW 56 of the four main thermal and volumetric properties of water: the density, the thermal expansion coefficient, the isothermal compressibility, and the heat capacity. The comparison to experiments shows that the present model gives equal or better agreement than the simulation models over the normal and supercooled liquid temperature range and does not have the fluctuation errors that simulations have, but the theoretical model has more parameters The model allows us to parse the experimental observables into hydrogen-bonding, caging, van der Waals, and non-interacting molecular components. Water is known to have a high heat capacity (ability to absorb thermal energy upon heating) among liquids of similar molecular size. The main conclusions from Figure 14b. are the following. In the normal liquid range, the high heat capacity comes from the breaking of two types of bonds: pairwise H bonds and Lennard−Jones-like contacts. Heating hot water near the boiling point leads to lower density, as it would for any LJ fluid, because heating hot water changes the contact interactions more than the H bonds. Figure  14c shows the same bulk properties as in Figure 14a except now computed as a function of pressure, not temperature. As increasing pressure squeezes water to become more compact (density increases and compressibility decreases), it crumples the hexagonal water cages breaking them into component pieces that just have pairwise water−water hydrogen bonding with little change to LJ and noninteracting water populations. Pressure decreases the heat capacity (bond-breaking capability) because although it melts out some cages it is also "freezing in" some pairwise H bonds. The thermal expansion coefficient increases with pressure because pressure melts out the rigid cages into fragmented H-bond pairs, which can be more readily squeezed together by pressure. Figure 15 shows that CageWater accurately reproduces the anomalous hallmark thermal and volumetric signatures of the LLPT, namely, the divergent increasing heat capacity and compressibility with lowered temperature. Moreover, this model gives the microscopic components of those observables. We find that the large diverging heat capacity is due to the water cages, which have dominant populations in cold and supercooled water. The heat capacity is the sum of two contributions for each state: the population of that state multiplied by the individual heat capacity. We also find that the negative thermal expansion of supercooled water is dominated by the cage term. Heating supercooled water shrinks the average volume by melting the cages, which are voluminous, and converts them to smaller H-bonded fragments, like breaking a glass jar into shards that pack more compactly. This same physics is reflected in the peak of the compressibility at the supercooling peak temperature. Our model indicates that the two liquids that are in equilibrium around −50 °C are cage structures and broken H-bonded pieces.

Conclusion and Future Perspectives
We developed an analytical theory of water and applied it to 2D MB, 3D MB water models and parametrized it for experimental data. We used it for explaining how the pVT properties of liquid water arise from water's hydrogen bonding and contacts. The theory predicts volumetric and energetic properties rather well, for experimental data it is more accurate than explicit simulation models yet is much faster to compute. Its simplicity and predictive power come from representing water using only four factors in the partition function, hydrogen bonds, Lennard-Jones contacts, noninteracting terms and cooperative cages, rather than as a more extensive density expansion. The analytical theory advances our understanding of water's structure−property relations in showing that water's long-known 2-density behavior is encoded in relatively infrequent cages which melt out strongly with temperature and pressure. This un-   derstanding of water structure−property relations may aid in engineering filtration, osmosis, and desalination materials, in better solvation models for drugs and biomolecule actions, and for interpreting planetary geochemistry and hydrological cycles. The challenge of the current model is in calculating dielectric permittivity 57,58 and solvation of polar molecules and ions. To be able to do this the model will have to be upgraded to version to include also charges on the water test particles, but this will add new parameters. Another challenge in theoretical modelling lies in developing the same kind of theory for other compounds like ionic liquids 55 , alcohols etc. We are expecting that all this can be done.

Acknowledgments
The financial support of the Slovenian Research Agency through Grant P1-0201 as well as to projects J7-1816, J1-1708, N1-0186 and N2-0067 is acknowledged as well as National Institutes for Health RM1 award RM1GM135136.