Thermodynamic response functions and Stokes-Einstein breakdown in superheated water under gigapascal pressure

Liquid water is the most intriguing liquid in nature, both because of its importance to every known form of life, and its numerous anomalous properties, largely magnified under supercooled conditions. Among the anomalous properties of water is the seeming divergence of the thermodynamic response functions and dynamic properties below the homogenous nucleation temperature (~ 232 K). Furthermore, water exhibits an increasingly decoupling of the viscosity and diffusion, upon cooling, resulting in the breakdown of the Stokes-Einstein relationship (SER). At high temperatures and pressures, however, water behaves more like a “simple” liquid. Nonetheless, experiments at 400 K and GPa pressures (Bove et al. (2011) Phys. Rev. Lett., 111:185,901) showed that although the diffusion decreases monotonically with the pressure, opposite to pressurized supercooled water, a decoupling of the viscosity and diffusion, larger than that found in supercooled water at normal pressure, is observed. Here, we studied the validity of SER and different pressure-dependent thermodynamic response functions, known to exhibit an abnormal behavior upon cooling, including the density, isothermal compressibility, and the thermal expansion coefficient along the 400 K isotherm up to 3 GPa through molecular dynamics simulations. Seven different water models were investigated. A monotonic increase of the density (~ 50%) and decrease of the isothermal compressibility (~ 90%) and thermal expansion (~ 65%) is found. Our results also show that compressed hot water has various resemblances to cool water at normal pressure, with pressure inducing the formation of a new second coordination sphere and a monotonic decrease of the diffusion and viscosity coefficients. Whereas all water models provide a good account of the viscosity, the magnitude of the violation of the SER at high pressures (> ~ 1 GPa) is significantly smaller than that found through experiments. Thus, violation of the SER in simulations is comparable to that observed for liquid supercooled water, indicating possible limitations of the water models to account for the local structure and self-diffusion of superheated water above ~ 1 GPa.


Introduction
Liquid water exhibits numerous thermodynamic and dynamic anomalies, most of which intensify upon cooling below the freezing point [1][2][3]. Among the most prominent anomalies of liquid water at normal pressure are the increase of the magnitude of its thermodynamic response functions, density (ρ), the isothermal compressibility (K T ), and thermal expansion coefficient (α P ). These anomalies, related to an increase of the entropy and volume fluctuations at low temperatures, opposite to simple liquids, prevent yet, a comprehensive understanding of liquid water, including its temperature and pressure dependence [2,4,5].
Several explanations were put forward to provide a coherent picture of water's anomalous behaviour. These include the ''stability limit conjecture'' (SLC) [6], the singularityfree hypothesis [7] or the liquid-liquid critical point (LLCP) [8]. The LLCP picture, the most widely investigated, foresees the co-existence of high-density and low-density liquid (HDL/LDL) forms of water below a second critical point located already below the homogeneous nucleation temperature (T h ~ 232 K). The existence of HDL and LDL forms of liquid water would then relate to the abovementioned anomalies, including the temperature of maximum density (ρ) at 4 °C, the minimum of the isothermal compressibility (K T ) at 46 °C, and a negative thermal expansion coefficient (α P ) below 4 °C [2]. Furthermore, in addition to the seeming divergence of the above thermodynamic response functions, the LLCP would explain a similar apparent divergence of transport properties such as the viscosity, in deep supercooled water.
The LLCP has been observed in molecular simulations of several water models [8][9][10][11] but since its estimated locus ( T ~ 170-190 K; P ~ 0.17-0.19 GPa) [10] is already within the so-called "no man's land" (101-231 K), this poses serious experimental limitations. Several experiments attempted to overcome this kinetic limitation by investigating nanoconfined water [12] and water droplets [13,14], for which T h is below that of bulk water. However, no experimental study could yet refute or confirm the existence of a LLCP in liquid water.
Among the dynamic anomalies of liquid water are the decrease of the viscosity, η, upon compression up to ~ 1 kbar, below ~ 10 °C (see Fig. 4 of ref. [2]), and the increase of the self-diffusion, D, at low pressures, below room temperature [2,15]. Further, an increasingly decoupling of diffusion and viscosity, upon cooling below 350 K, at normal pressure, is observed, resulting in a breakdown of the Stokes-Einstein relationship (SER) [16][17][18][19][20][21][22][23][24]. The breakdown of the SER at low temperatures occurs because the viscosity increases faster than 1/D. Thus, such a violation is magnified when the self-diffusion of low tetrahedrality populations is considered, since the latter exhibit a less pronounced decrease of the self-diffusion coefficient with the temperature [24].
Recently, three of us showed that the breakdown of SER can be quantitatively explained using the translational jumpdiffusion (TJD) approach [22,[25][26][27][28][29][30]. The TJD splits the total diffusion coefficient on a jump-diffusion term, emanating from the translational jump movements, and a residual diffusion term, the remaining part of the diffusion contributed by small step displacements. Interestingly, the residual diffusion remains almost strongly coupled to the viscosity in supercooled water. Lowering the temperature, pressure, or both, facilitates a dramatic growth of the contribution of jump-diffusion in supercooled water leading to a gradual collapse of the SER.
Unlike supercooled water, superheated water behaves more like a simple liquid and the LLCP hypothesis should not have to be invoked to interpret its thermodynamic and structural properties. Nevertheless, the dynamics of superheated water at GPa pressure has been shown [31] to exhibit much larger deviations to the SER than those observed for supercooled water, at least, down to 240 K where both viscosity [16] and diffusion [32,33] data is available. Bove et al. [31], measured, through quasielastic incoherent neutron scattering, the translational and rotational diffusion coefficients of water along the 400 K isotherm up to 3 GPa, which is the melting point of ice VII. The translational diffusion D was found to strongly decrease with pressure although the slowdown rate diminishes above 1 GPa. However, this reduction of D is slower by a factor of nearly two than the increase of viscosity η of water at the same pressures. Therefore, the product Dη increases with the pressure, violating the SER. This breakdown of the SER is particularly strong above 1 GPa and its relationship with the putative LLCP is not known.
The main goal of this work is twofold: to investigate whether the abovementioned water models can quantitatively predict the SER breakdown at high-temperatures and high-pressures and assess putative structural and thermodynamic anomalies in superheated water, that might accompany the SER violation.

Molecular dynamics
MD simulations were performed for 1500 water molecules in a cubic box with periodic boundary conditions, for the distinct water models (SPC/E, SWM4-NDP, TIP3P/FB, TIP4P/FB, TIP5P, TIP4P-Ew, and TIP4P/2005). The initial configuration was generated with the Packmol software [40], and the MD simulations were carried out with the GROMACS program [41]. The systems were simulated at 400 K and twelve different pressures: 0.01, 0.1, 0.2, 0.3, 0.5, 0.8, 1.2, 1.5, 1.9, 2.4, 2.7, and 3.0 GPa. Following steepest descent energy minimization and a 100 ps MD in the NVT ensemble, the systems were further equilibrated for 10 ns in the NPT ensemble. The trajectories were then propagated for 40 ns in the NPT ensemble. Electrostatic interactions were computed via the particle-mesh Ewald (PME) method [42]. A cut-off of 1 nm was used for non-bonded van der Waals and the PME real space electrostatic interactions. The temperature was controlled via the Nosé-Hoover thermostat [43,44] with a coupling constant of 0.5 ps. The pressure was controlled isotropically with the Parrinello-Rahman barostat [45] with a coupling constant of 1 ps. The Verlet-Leapfrog algorithm was used to solve the equation of motion with a time-step of 2 fs.

Dynamic properties
The diffusion coefficient was estimated from the mean square displacement through the Einstein relation [40], where the < > indicate an ensemble average of the mean square displacement (MSD), r i (t) and r i (0) are the positions of the ith oxygen atom at time t and the origin, respectively, and D is the diffusion coefficient.
The shear viscosity, η, was calculated through integration of the stress (or pressure) tensor time correlation function ( P ) [46,47], where P is the symmetrized traceless portion of the stress tensor, , given by [42,48], is the Kronecker delta and there are six, out of the nine, distinct P elements. This expression improves the statistics over the original Green-Kubo formula which involves the average over only the three different off-diagonal elements, = , of the stress tensor [40,49], also used here for comparison purposes; no significant differences were found between the viscosity calculated with Eqs. (2) and (4). For the definition of the stress tensor elements for periodic systems treated with the Ewald or the particle-mesh Ewald methods see refs [50][51][52].

Thermodynamic properties
Various thermodynamic response functions, known to exhibit an abnormal behavior upon cooling, were calculated for the distinct models, to assess their pressure dependence in superheated water. These included the density ρ, the isothermal compressibility K T , and the thermal expansion coefficient P . The density ( = m∕V ) is the amount of mass (m) of a substance per unit volume (V) and is sensitive to temperature and pressure, especially for fluid systems. For simple liquids, a pressure increase, or temperature decrease, results in a volume reduction and, therefore, an increased density. The density of water, although increasing with pressure, decreases upon cooling below 4 °C; hence, a volume increase is accompanied by an entropy decrease, a paradoxical behaviour not found in simple liquids. The isothermal compressibility,K T = −( ln V∕ P) T , provides a measure of the change of volume as a response to a change in pressure, at constant temperature. The thermal expansion coefficient, P = ( ln V∕ T) P , provides a measure of the change of volume of a fluid or solid as a response to a change in temperature, at constant pressure. These properties were calculated using the respective fluctuation formulas as discussed by Allen and Tildesley [53]. The pressure dependence of the above properties for the distinct water models is displayed in Fig. 1 and compared with the available experimental data [54]. The density of water increases monotonically (as expected) with the pressure increase (Fig. 1a), consistent with the experimental trend. Moreover, the results are similar for the distinct models, except for SWM4-NDP and TIP5P, which predict, respectively, a lower and higher density at high pressures, compared to the experimental data. For the remaining five water models the density is in very good agreement with the experiments. For most models, ρ is almost coincident up to ~ 1 GPa. Thus, the density only diverges slightly, among the distinct water models, at high pressures (> 1 GPa), and an increase ~ 50% is found between the lowest and highest pressures studied, for most models.
The isothermal compressibility, K T (Fig. 1b), and thermal expansion coefficient, P (Fig. 1c), decay monotonically with the pressure for all water models. The largest difference among the water models is found for TIP5P, which estimates significantly larger K T and α P , especially at low pressures. Except for TIP5P and SWM4-NDP the simulated P matches well the experimental data at low pressures, whereas at high pressures deviations are also observed for TIP4P-Ew and TIP4P-fb. The TIP3P-fb, TIP4P/2005, and SPC/E models outperform the remaining 44 Page 4 of 9 ones for the entire pressure range. For most models a decrease of K T and α P around 90% and 65%, respectively, is found.

Diffusion and viscosity
The self-diffusion coefficient was estimated from the MSD (see Eq. 1). The MSD comprises three characteristic domains: a sort-time ballistic regime ( MSD ∝ t 2 ), an intermediate-time sub-diffusive regime ( MSD ∝ t ;0 < α < 1 ), and a long-time linear diffusive regime (MSD ∝ t ), from which the diffusion can be estimated (see Eq. 1).
The self-diffusion and viscosity coefficients are displayed in Fig. 2a and b, respectively. The self-diffusion coefficient at pressures below ~ 2.5 GPa is consistent with the experimental data [31], except for the TIP5P water model. At higher pressures (> 2.5 GPa), a faster decay of the self-diffusion coefficient, compared with the experiments, can be seen. Opposite to D, the experimental viscosity coefficient [31,55,56] does not exhibit such a marked downtrend variation at high pressures; "experimental" viscosity data are the same used by Bove et al. [31] and has been interpolated from experimental data at 373 K and 437 K [56]. Again, except for TIP5P the simulation viscosity is in good agreement with the experimental data for every water model (see Fig. 2b). The faster decay of D at high pressures, and the viscosity accord, compared with the experiments, indicates that model water should obey the SER more closely than real water.
The SER, D ∕T = c , where the constant c = k B ∕C r s ( k B is the Boltzmann constant, C is a constant, and r s is the hydrodynamic radius of the "molecule"), predicts a constant Dη product at constant temperature. Figure 2c  of the SER. However, the magnitude of this violation is significantly lower than that found from the product of the experimental D and η. The experimental Dη increases slowly up to 1 GPa, after which a sharp increase can be seen up to ~ 2.5 GPa. Above 2.5 GPa an even steeper increase is observed, denoting a major decoupling between the diffusion and viscosity. Thus, although the water models predict the decoupling phenomenon, this is moderate compared to that found for real water for pressures above ~ 1 GPa. This possibly reveals the limitations of the water models to accurately capture the diffusion-viscosity decoupling of pressurized hot water. The reason for this behavior is mostly related with the prediction of a faster decay of the diffusion coefficient at high pressures, compared to real water. The latter, in turn, could be associated with a poor description of the structure of water at high temperatures and pressures. We now turn attention to the analysis of the structure of water.

Structure
Bove et al. attributed the SER breakdown to the invariance of the number of hydrogen bonds (H-bonds) of water under pressure, associated with a rigid first coordination sphere, also based on MD simulations with the TIP4P/2005 water model [31]. Herein, we also calculated the average number of H-bonds per water molecule. A widely employed geometric H-bond definition based on two distances and an angle, was adopted: r ODOA < 3.5 Å, r OAHD < 2.45 Å, and the angle H D O D O A < 30°; where "D" and "A" denote the H-bond donor and acceptor, respectively [57][58][59][60][61]. Here, r ODOA and r OAHD are the positions of the first minimum of the O-O and O-H radial distributions functions (rdfs), which will shift slightly with pressure (see below). Furthermore, these values are different for the distinct water models. Nonetheless, because the H-bond breaking/forming dynamics is not controlled by translational [62], but rather by rotational [59,60] motions, the number of H-bonds does not change significantly with the abovementioned shift. Thus, we chose to use the same H-bond definition for all water models and across the pressure interval. Figure 3 shows that the average number of geometric H-bonds (n HB ) increases monotonically with the pressure. A relatively rapid increase is observed at low pressures whereas a moderate increase occurs at high pressures. This is consistent with the decrease of the self-diffusion coefficient. The fact that n HB does not plateau out explains the rapid decrease of the simulated self-diffusion coefficients at high pressures, compared to real water, resulting in a moderate violation of the SER. Notice that the self-diffusion coefficient of real water (see Fig. 2a) decreases more slowly at high pressures, nearly exhibiting a plateau above ~ 2.4 GPa. This "plateau", not observed in the simulations (see Fig. 2b), is responsible for the prominent violation of the SER by real water, depicted in Fig. 3b, and suggests that the local structure of real water, in particular the number of H-bonds, should be nearly insensitive to the pressure above ~ 2.4 GPa.
The n HB for all water models (except TIP5P) are nearly coincindent, implying only minor differences among these water models with respect to H-bond formation. The TIP5P water model predicts a much lower number of H-bonds compared to the other models. However, a larger growth upon compression is observed. Thus, the increase in the number of H-bonds in moving from 0.001 to 3 GPa is about 30% for TIP5P, while for the other water models is ~ 12%. This is consistent with the more steep decline of the self-diffusion coefficient of TIP5P water with the pressure.
These results show that while H-bonds are gradually broken upon compression in supercooled water, a pressure increase in superheated water increases geometric H-bonding. A nearly constant average number of geometric H-bonds is found in pressurized polarizable model water at room temperature [63].
Further insight into the structure of water can be gauged from water's partial rdfs. The rdf measures the probability of finding a particle at a distance r from some reference particle, revealing the local number density variation along the distance from the reference particle. The oxygen-oxygen, g OO (r), and oxygen-hydrogen, g OH (r), partial rdfs, are displayed in Figs. 4 and 5, respectively. Figure 4 shows that as pressure increases up to ~ 0.5 GPa, the second coordination sphere vanishes, with a new one starting to form at a higher distance upon further compression. This is accompanied by a broadening of the first coordination sphere. As the second coordination sphere becomes more and more prominent upon compression, it also starts shifting to shorter distances. A compression of the first coordination sphere concerning the position of the first minimum can also be seen above 0.5 GPa. The depth of the first minimum, in turn, starts to increase (i.e., deeper minimum) significantly above 0.5 GPa denoting a lower diffusion coefficient as water molecules are found less frequently crossing between the first and second coordination spheres. The g OH (r) (see Fig. 5) shows a first peak around 0.18 nm, corresponding to the nearest H atoms of neighbouring water molecules, which can form H-bonds with the reference O atom, and a second peak at around 0.32 nm. The first and second peaks undergo a shift to lower distances upon compression. However, while the height of the first peak slightly decreases that of the second increases. Moreover, the depth of the first minimum decreases (i.e., less deep) while that of the second exhibits the opposite behaviour.
The concomitant decrease of the height of the first peak and depth of the first minimum for all water models, explain the near constancy of the number of H neighbours across the low-and high-pressure ranges (see Fig. 6). The emergence of a clear third peak is also visible at high pressures, again revealing some resemblances with the structure of liquid water upon cooling. The first shell coordination number (CN) of water (O-O and O-H), calculated by integrating the rdfs are shown in Fig. 6. Clearly there is a sharp increase of CN as the pressure increases to 0.5 GPa, which is caused by the sudden shift to longer distances of the first minimum in the rdf. However, the O-H coordination number remains almost insensitive to the pressure, for the reasons previously noted.
We have also calculated the tetrahedral order parameter q for the TIP4P/2005 model at all pressures using the following equation [64,65].  , the distributions of q, at six representative pressures. At 0.01 GPa the tetrahedral order distribution peaks at around q ~ 0.48 and exhibits a shoulder around q ~ 0.7. With increasing pressure, the peak height increases, and the shoulder gradually vanishes such that the shoulder of the distribution is almost invisible at the highest pressure P = 3 GPa. Figure 7b plots the average tetrahedral order parameter < q > as a function of pressure. The tetrahedrality < q > remains almost constant up to 0.5 GPa followed by a steep decrease upon further compression. Therefore, although the number of H-bonds increases with increasing pressure, the tetrahedrality decreases, indicating the formation of more strained and possibly higher energy H-bonds. These structural changes should impact the diffusion-viscosity decoupling of highly pressurized water at 400 K, possibly inducing a faster decline of the diffusion coefficient, compared to the experiments.

Conclusions
Liquid water remains one of the most studied liquids both because of its importance to multiple chemical processes and its peculiar behavior, especially in supercooled metastable conditions. Opposite to supercooled water, superheated water behaves more like a simple liquid. Nonetheless, experimental neutron diffraction experiments revealed a large decoupling of the viscosity and diffusion upon compression. The possible connection between this decoupling and the putative second critical point of liquid water is unknown. Herein, to gain insight into the behaviour of superheated model water, molecular dynamics simulations at 400 K up to 3 GPa were performed with seven water models: SPC/E, TIP3P-FB, TIP4P-FB, TIP5P, SWM4-NDP, TIP4P-Ew, and TIP4P/2005.
Our results show that the TIP5P water model predicts the most disparate results among the models investigated, seemingly associated with a lower number of H-bonds. This model predicts a higher isothermal compressibility and thermal expansion coefficient as well as a lower viscosity and a higher self-diffusion coefficient, especially at low pressures. Upon compression the thermodynamic response functions and transport coefficients of TIP5P water become less dissimilar from the remaining water models. This is linked with a steeper increase of the number of the H-bonds with pressure, making TIP5P water's H-bond network more similar to that of the remaining models, at high pressures.
Although the number of H-bonds increases, the tetrahedrality of the H-bond network decreases with increasing pressure indicating a more strained network.
The distinct water models estimate a more moderate violation of the SER, than that found by Bove et al. [31]. Thus, a diffusion-viscosity decoupling similar to that found in real supercooled water, down to 242 K, is observed. The reason is that water models predict a more pronounced decrease of the self-diffusion coefficient at high pressures, than experiments. This suggests that water models commonly used to study liquid water, including metastable supercooled water, cannot accurately describe the structural and dynamic responses to compression, of superheated water. This problem should be revisited using other water models and/or ab initio MD to gain further insight into the structural source of the strong diffusion-viscosity decoupling observed in real superheated water.
Author contributions SD contributed to the methodology, analysis, and writing original draft. AM contributed to the analysis, and writing