We use molecular dynamics simulations to study relations between thermodymamic, structural, and dynamical properties of TIP4P/2005 water models with systematically reduced partial charges and, thus, weaker hydrogen bonds. Observing a crossing of isochores in the P–T diagram, we show that these water-like models have a readily accessible liquid–liquid critical point (LLCP) associated with a transition between high-density liquid (HDL) and low-density liquid (LDL) forms and determine the dependence of the critical temperature Tc, pressure Pc, and density ρc on the charge-scaling factor from fits to a two-structure equation of states. The results indicate that the water-like models exhibit liquid polyamorphism in a wide range of interaction parameters. Considering elongated systems, we observe a decomposition into extended and stable HDL-like and LDL-like regions at appropriate pressures and low temperatures and analyze the respective structural and dynamical properties. We show that the diverse local order results in very different correlation times of local dynamics, while the fragility is hardly changed. The results yield insights into the origin of a dynamical crossover, which is observed when lowering the temperature along isobars and was previously interpreted in terms of a fragile-to-strong transition. Our findings imply that the effect does not involve two liquid phases with an exceptionally large difference of the fragility but rather a high temperature dependence near the LLCP results from a rapid conversion from HDL-like environments with faster dynamics to LDL-like ones with slower dynamics.

The anomalies of water are remarkable examples of competing many-particle effects. Most researchers agree that they result from the competition between hydrogen bonding and other intermolecular interactions and originate in the supercooled regime of the liquid,1 where thermodynamic response functions appear to diverge upon cooling below the homogeneous nucleation temperature,2 which was recently specified as TH ≈ 231 K in studies on micrometer-sized water droplets.3 In particular, it was conjectured that the anomalies of water are caused by a second critical point,4 which terminates a liquid–liquid phase transition (LLPT) line at reduced temperatures T and elevated pressures P. The proposed liquid phases of water differ in density motivating the denotations as high-density liquid (HDL) and low-density liquid (LDL).5–9 In addition, it was argued that they show diverse structures, e.g., the tetrahedral order is lower in the HDL than that in the LDL phase, and distinguishable dynamics.10 

The rapid crystallization of water, however, interferes with a straightforward characterization of the structural and dynamical properties of the proposed HDL and LDL phases. The term no man’s land was coined to denote these difficulties pertaining to the studies on supercooled water in the temperature range of ∼150 K–231 K. Therefore, despite recent progress in the reduction of the no man’s land,11–14 the existence of the HDL and LDL states and, thus, of the LLPT line is still a subject of controversial discussions.15,16 Notwithstanding, experimental evidence for the liquid polyamorphism of water is provided, e.g., by the existence of amorphous ice phases, which are separated by a first-order-like transition and show glass transitions at different temperatures Tg.5–9,17

A liquid polyamorphism is not unique to water. A recent review article provides a comprehensive overview.18 Experimental evidence for this phenomenon was also reported for, e.g., carbon,19 phosphorous,20,21 and sulfur.22 In simulation work, liquid polyamorphism was observed for various systems ranging from patchy particles with tetrahedral functionality23 to the spherically symmetric two length scale Jagla model.24 Moreover, even simple theoretical approaches can produce this feature.25,26 These studies imply that a LLPT line and the associated liquid–liquid critical point (LLCP) can occur when spacious low-energy and low-entropy local structures such as tetrahedral order exist. In such a situation, cooling promotes the LDL phase consisting of these structural motifs, while pressurizing favors the HDL phase with more compact particle arrangements.

Molecular dynamics (MD) simulations are well suited for investigations of liquid polyamorphism.27–31 In studies on water, it is advantageous that high cooling rates allow one to avoid crystallization. A possible drawback is that the LLCP, if existent, often occurs on the edge of the temperature–pressure region accessible to current computational approaches. Nevertheless, MD simulations were successfully used to scrutinize the existence of a LLCP for various water models.32 While evidence for a LLPT was reported for the TIP4P/2005,33–37 TIP4P/Ice,37 TIP5P,36 WAIL,38 E3B3,39 and ST24,40,41 models, this phenomenon does not occur for the mW model, where the entropy of mixing is too high for a phase separation.42 For the SPC/E model, a LLCP was proposed to exist in a parameter range inaccessible to current MD simulations.43 Even if observed, appropriate equilibration of the systems may be an issue owing to the slow structural relaxation near the LLCP, and opposing views can exist relating to the location of the LLCP33–36 and the effects of ice-like clusters,44,45 as was reported for TIP4P/2005 water.

For the theoretical description of polyamorphism, the two-structure equation of states (TSEOS) formalism was successfully applied.36,41,42,46–52 The TSEOS approach assumes that the Gibbs energy of HDL–LDL mixtures is given by the contributions of the individual phases supplemented by a nonideal mixing term, which can be enthalpic35,41,53,54 or entropic42,48,55 in nature. Although it is based on a mean-field approximation, it often enables a good description of the data and a reliable determination of the location of the LLCP,52 i.e., of the critical temperature Tc, pressure Pc, and density ρc, including cases where a modest extrapolation to inaccessible state points is required.

This literature survey suggests that liquid polyamorphism is more abundant than anticipated. However, it is still elusive for which types of particle interactions the phenomenon occurs and at which state points the LLPT is located. A major advantage of MD simulations is that the interactions can be varied continuously and that temperature, pressure, and density changes are straightforward. These capabilities were exploited to ascertain the LLPT of various liquids. For the Stillinger–Weber potential of silicon,56 it was found that the LLPT is erased when the three-body repulsion parameter is reduced.57,58 For the ST2 model of water, it was observed that the LLPT becomes more prominent when the bond angle is varied.59 Such variations are of particular relevance for water, for which a unique model does not yet exist so that the dependence of the properties on the parameters is an important issue.

Another efficient way to alter the interactions in chemically realistic models is to vary the partial charges of the atoms. While the model of Woodcock, Angell, and Cheeseman (WAC) of silica did not show a LLCP in the accessible region of the phase diagram,28 a clear crossing of isochores in the P–T phase diagram and extensive free energy calculations indicated a LLCP when the partial charges of this model are decreased in magnitude.29,60 Moreover, charge scaling proved to be very useful to investigate the role of the particle interactions in the glass transition of water,61,62 silica,63 and an ionic liquid.64 In particular, this approach allows one to tune the strength of the hydrogen bonds and the degree of the tetrahedral order in water models.61,62

The conjectured liquid polyamorphism is assumed to result also in peculiar dynamical behaviors of water. For example, the dynamics of water speeds up when the pressure is increased.1 To rationalize this anomalous effect, one may argue that the HDL replaces the LDL upon pressurizing and that the former shows faster dynamics than the latter. Consistent with this conjecture, it was reported that the HDL has a lower glass transition temperature Tg than the LDL,17 but the interpretation of the experimental observations is still controversial.65 In an extensive simulation study on the TIP4P/2005 water model, different glass transition temperatures were clearly identified, whereby the pressure dependence Tg(P) was found to be normal for the HDL, but anomalous for the LDL.66 Furthermore, a LLPT may be accompanied by a fragile-to-strong (FS) dynamical transition.67 Explicitly, assuming that the HDL and LDL are fragile and strong liquids, i.e., deviate from and conform to Arrhenius behaviors, respectively, a FS transition is expected in temperature-dependent studies, when, depending on the pressure, the LLPT line or the Widom line are crossed, where the latter marks the position of the maxima of thermodynamic response functions in the one-phase region above the LLCP.68–71 However, alternative explanations for crossovers in the temperature-dependent dynamics of water were also put forward.72 For example, the effect was attributed to an onset of hopping motions or secondary relaxations.70,73 MD simulations proved to be very useful to investigate the dynamical behavior of the ST2 water model near the LLPT.74 

Here, we perform MD simulations to ascertain the thermodynamic, structural, and dynamical properties of TIP4P/2005 water models with reduced partial charges and, thus, weaker hydrogen bonds. Observing an isochore crossing in the P–T phase diagram, it is shown that the studied water-like models have LLCPs. Performing a TSEOS analysis, the dependence of the critical parameters on the magnitude of the partial charges is determined. The analysis reveals that charge scaling shifts the LLCP to computationally better accessible state points. This effect is a major advantage of our approach and allows us to follow the evolution of the structure and dynamics of the charge scaled models across the LLCP. To rationalize the findings, we consider systems, which decompose into extended HDL-like and LDL-like spatial regions and determine the structure and dynamics in these domains. The findings show that a FS-like dynamical crossover is mainly related to a rapid conversion from HDL-like environments with faster dynamics to LDL-like ones with slower dynamics.

We performed MD simulations for several charge scaled TIP4P/2005 water models. Specifically, we started from the original TIP4P/2005 model75 and scaled the partial charges of all atoms by a factor q ranging between 0.86 and 0.91, while all other interaction parameters remained unchanged. Thus, the water-like molecules of the charge scaled TIP4P/2005 models are less polar and form weaker hydrogen bonds.

The simulations were conducted using the GROMACS simulation package.76,77 Periodic boundary conditions were employed and the leapfrog integrator with a time step of 2 fs was used. The particle mesh Ewald method78 was utilized to calculate the Coulomb and Lennard-Jones interactions where the cutoff was set to 1.0 nm and the Fourier spacing to 0.16 nm−1. All systems were carefully equilibrated prior to data acquisition. Specifically, the equilibration runs and the production runs were at least several ten times longer than the rotational correlation time at the given temperature.

The majority of the studied systems were cubic and comprised N = 2000 molecules, leading to box lengths of ∼4 nm. Most simulations of these systems were performed in the isothermal–isochoric (NVT) ensemble, applying the Nosé–Hoover thermostat.79,80 For simulations in the isothermal–isobaric (NPT) ensemble, we additionally employed the Parrinello–Rahman barostat.81 Moreover, we utilized elongated rectangular simulation boxes, which facilitate a decomposition into two phases because they allow for smaller interface-to-volume ratios than their cubic counterparts.29,82 These coexistence simulations were performed in the NVT ensemble for q = 0.86 and q = 0.88. We chose a ratio of the box dimensions of Lx : Ly : Lz ≡ 1 : 1 : 3 for systems comprising N = 2000 molecules and carried out additional simulations for N = 12 000 molecules and a ratio of 1:1:6. With these choices, the smaller box dimensions were Lx = Ly ≥ 2.7 nm and Lx = Ly ≥ 4.0 nm for the smaller and larger elongated systems, respectively.

To evaluate the isochores and determine the LLCP, we follow previous studies35,36,41,51 and perform a TSEOS analysis. This approach assumes that the liquid consists of an interconverting mixture of a high-density structure A with a fraction 1 − x and a low-density structure B with a fraction x. The molar Gibbs energy is written as35,36,41,51,52

G=GA+xGBA+RTxlnx+(1x)ln(1x)+ωx(1x).
(1)

Here, GBA = GBGA is the difference of the Gibbs energies GA and GB of the pure states. The third summand accounts for the entropy of mixing the structures A and B, where the first two terms of this contribution describe ideal mixing and the last one nonideal contributions. Because the determination of the LLCP is of interest, it is convenient to express temperature and pressure as dimensionless variables, ΔT^=TTcTc and ΔP^=PPcρcRTc. GA is modeled as a third-order expansion around the LLCP,

GA=RTcm,nĉmnΔT^mΔP^n,
(2)

with ĉmn as adjustable coefficients. GBA determines the A ⇌ B equilibrium; explicitly, the equilibrium constant amounts to K(T, P) = exp(−GBA/RT). It is approximated as

GBART=λ(ΔT^+aΔP^+bΔT^ΔP^).
(3)

The meaning of the parameters λ, a, and b was discussed in detail in previous studies.35,41,51 The LLPT, LLCP, and the Widom line, i.e., the continuation of the LLPT line beyond the LLCP, are determined by the condition ln K(T, P) = 0.

The nonideal contribution to the entropy of mixing, ωx(1 − x), is assumed to be symmetric in x. The parameter ω is a measure for the relevance of this effect. While a temperature-independent value of ω indicates an entropic origin, ω1T is found when the nonideality is enthalpic in nature. In previous works on the TIP4P/2005 water model,35,51 it was shown that the data are well described by the form

ω=2+ω0ΔP^T^,
(4)

with T^=T/Tc and ω0 as the adjustable parameter. Therefore, we also apply this form to all modified TIP4P/2005 models.

The equilibrium fraction xe is defined by the condition

1RTGxP,T=GBART+lnx1x+ω(12x)=0.
(5)

This equation has to be treated numerically. When describing isochoric data, the molar volume can be obtained from the partial derivative of the dimensionless Gibbs energy Ĝ=GRTc with respect to P^=PρcRTc. In a recent work,51 the TSEOS analysis was augmented with contributions from the liquid–vapor spinodal to improve the performance in a wide range of state points. We found that this extension was not necessary in our study, which was restricted to the vicinity of the LLCP, in particular, to T < 1.25Tc. More details about our TSEOS analysis can be found in the supplementary material.

1. Structure

We use two local structure identifiers to distinguish between HDL-typical and LDL-typical environments. In the following, we will refer to these particle arrangements as high-density (HD) and low-density (LD) local structures, respectively. For their calculation, we consider a molecule i and obtain the distances di,j to its neighboring molecules j from the positions of the respective oxygen atoms. Moreover, we determine the number n(i) of neighbors within the cutoff distance di,j < 0.37 nm. Ordering the values of di,j according to increasing distance, the local structure index (LSI) of molecule i is defined as83 

I(i)=1n(i)j=1n(i)Δ(i,j)Δ¯(i)2.
(6)

Here, Δ(i, j) = di,j+1di,j and Δ¯(i) is the average of these differences over all n(i) neighbors. For tetrahedral LD local structures, the difference between the di,j values is small for the first four neighbors but large between the fourth and fifth neighbors, resulting in large LSI values. For HD local structures, by contrast, the distribution of the Δ(i, j) is narrower, and thus, the LSI values are smaller. Moreover, we delineate a simple structure identifier N4(i) based on the number n(i) of neighboring molecules. In detail, we define N4(i) = 0 if n(i) > 4 and N4(i) = 1 if n(i) ≤ 4. Thus, the average value of N4 is small in HDL-like regions with more than four neighbors in the next neighbor shell, while it is large in tetrahedral LDL-like neighborhoods. Likewise, the mere number of neighbors was used to define a structure identifier in a previous computational approach to water properties near the LLCP.74 In the supplementary material, we also analyze the distance to the fifth neighbor, d5di,5, and exploit that low and high values are indicative of HD and LD local structures, respectively.

2. Dynamics

To characterize the translational motion of the water-like molecules, we analyze the time-dependent positions ri(t) of their oxygen atoms. Specifically, we calculate the mean-square displacement (MSD)

r2(t)=|ri(t+t0)ri(t0)|2,
(7)

where the pointed brackets denote the averages over all molecules i and various time origins t0, and determine the self-diffusion coefficient D by fitting the MSD in the diffusive regime at long times to r2(t) = 6Dt.

For the analysis of the reorientation dynamics, we consider the O–H bond vectors of the water-like molecules and compute the rotational correlation function (RCF),

F1(t)=ei(t+t0)ei(t0).
(8)

Here, the unit vectors ei(t) describe the time dependence of the O–H bond orientations. To quantify the time scale of the molecular reorientation, we fit the RCF F1(t) beyond the vibrational regime to a scaled Kohlrausch–Williams–Watts (KWW) function,

Aexptτkwwβkww.
(9)

Here, the factor A is used to consider that vibrational motion causes a minor loss of correlation prior to the onset of structural relaxation. Moreover, τkww and βkww are the time constant and the stretching parameter of the KWW function, respectively. From the fit results, we calculate the mean correlation time τ=1A0F1(t)dt according to τ = (τkww/βkww)Γ(1/βkww), where Γ(x) denotes the gamma function.

3. Spatially resolved analyses

In the simulations for elongated boxes near the LLCP, we observe a time-dependent decomposition into HDL-like and LDL-like fractions along the z dimension. To characterize the space–time characteristics of the structural and dynamical properties in these regions, we follow the respective behaviors in layers, which are parallel to the xy plane and have a width of Δz = 1 nm, along the trajectories. To analyze the density fluctuations, we define the layer-averaged density ρL(z, t) based on the number of molecules in the layer centered about the position z at time t. Likewise, we average over the molecules in the respective layers to obtain z- and t-dependent values of layer-averaged structural identifiers, IL and N4,L. For N = 2000, these averages include about 250 molecules. To analyze the overall temporal aspect of the fluctuations, we calculate correlation functions of the layer-averaged quantities,84 

Cξ(t)=ξL(t+t0)ξL(t0)ξL2σξ2,
(10)

where ξL denotes the respective random variable and σξ is its standard deviation. Moreover, we perform KWW fits to obtain mean time constants τξ of these fluctuations. In these analyses, vibrational motions and numerical uncertainties lead to rapid changes in the numbers of particles within the layers, which can, however, be removed by averaging the coordinates over the time scale of the vibrational motion. To investigate the local dynamics–structure relations, we calculate F1(t) separately for molecules in structurally different environments. Specifically, we distinguish between molecules characterized by various layer-averaged densities or structures at the time origin t0, where the latter are obtained from the layers centered about the instant molecular positions, as described above.

Commonly, liquids show a monotonic decrease in the pressure P when the temperature T is reduced along an isochore. However, when two liquid phases with different densities compete, the pressure can increase upon cooling. A crossing of isochores in the P–T diagram indicates a LLPT.28,29,60,85 Therefore, we calculate the pressure along isochores of the charge scaled water models from the virial stress. The P–T diagrams for the lowest (q = 0.86) and highest (q = 0.91) studied partial charges are shown in Fig. 1. For both q values, we observe that the pressure decreases when the temperature is reduced along isochores for sufficiently high densities, while it increases along those with sufficiently low densities in the studied temperature range. More importantly, the isochores for appropriate intermediate densities cross, indicating the existence of a LLPT. At low temperatures, the isochores end where the equilibration times exceed 2 µs. We see that the metastable region and the expected spread of the crossing isochores below the LLPT can be explored to some degree. The highest temperature, for which an isochore crossing is observed, and the corresponding pressure and density mark the LLCP. Qualitatively, the same behavior is observed for all studied scaling factors between q = 0.86 and q = 0.91; see the supplementary material.

FIG. 1.

Isochores of charge scaled TIP4P/2005 water models in the P–T diagram: (a) q = 0.86 and (b) q = 0.91. The respective densities are indicated. The colored solid lines are the results of the TSEOS analysis and involve the data points, which are included in the fitting procedure. The black crossed circles mark the positions of the LLCP, the black solid lines are the LLPT lines, and the black dashed lines are the Widom lines.

FIG. 1.

Isochores of charge scaled TIP4P/2005 water models in the P–T diagram: (a) q = 0.86 and (b) q = 0.91. The respective densities are indicated. The colored solid lines are the results of the TSEOS analysis and involve the data points, which are included in the fitting procedure. The black crossed circles mark the positions of the LLCP, the black solid lines are the LLPT lines, and the black dashed lines are the Widom lines.

Close modal

Taking the isochore crossing as evidence for the existence of two competing liquid states, we use the TSEOS approach to further analyze the simulation data. The results of this analysis are included as solid lines in Fig. 1. We see that the TSEOS approach well describes the simulation results for q = 0.86 and q = 0.91 in the vicinity of the LLCP. Similar agreement is observed for the other studied scaling factors, in harmony with previous findings for various model liquids with the tetrahedral order.35,36,41,51 The parameters resulting from the TSEOS fits for all studied scaling factors q are compiled in the supplementary material.

In Fig. 2, we depict the location of the LLCP for the studied water-like models, as obtained from the TSEOS analysis. We see that the critical temperature Tc, pressure Pc, and density ρc continuously decrease when the scaling factor q is reduced, consistent with previous findings for the WAC model of silica.60 Qualitatively, our findings can be rationalized by the fact that the electrostatic interactions and, thus, the hydrogen bonds are weaker for small values of q. Because of this q dependence of the interactions, the relevant thermal energy and, thus, Tc decrease. Moreover, weaker hydrogen bonds lead to a weaker tetrahedral order so that less pressure is required to push a fifth molecule into the next neighbor shell. Specifically, the LLCP is located at elevated pressures when the magnitude of the partial charges is lowered by less than 10% but at negative pressures when the molecular polarity is further reduced. Finally, the electrostatic interactions pull the molecules closer together than expected based on the Lennard-Jones parameters so that the density is reduced at small scaling factors q. To validate the outcome of the TSEOS analysis, we follow the evolution of the density ρ upon isobaric cooling at atmospheric pressure in the supplementary material. As expected, we find that ρ decreases very strongly near Tc for the higher scaling factors, for which Pc ≈ 1 bar, whereas such drops are not observed for the lower q values.

FIG. 2.

Dependence of the LLCP obtained from the TSEOS analysis on the charge-scaling factor q: (a) critical temperature Tc, (b) critical pressure Pc, and (c) critical density ρc. The solid lines are linear fits. (d) Rotational correlation times τ and jump correlation times a2/6D at the LLCP, where the latter were calculated from the self-diffusion coefficients of the oxygen atoms using a jump length of a = 0.325 nm, similar to the oxygen–oxygen next neighbor distance. The results are obtained from the dynamical behaviors in the vicinity of the LLCP by interpolation; see the text for details.

FIG. 2.

Dependence of the LLCP obtained from the TSEOS analysis on the charge-scaling factor q: (a) critical temperature Tc, (b) critical pressure Pc, and (c) critical density ρc. The solid lines are linear fits. (d) Rotational correlation times τ and jump correlation times a2/6D at the LLCP, where the latter were calculated from the self-diffusion coefficients of the oxygen atoms using a jump length of a = 0.325 nm, similar to the oxygen–oxygen next neighbor distance. The results are obtained from the dynamical behaviors in the vicinity of the LLCP by interpolation; see the text for details.

Close modal

To estimate the location of the LLCP for the regular TIP4P/2005 model, we extrapolate the basically linear dependence of the critical parameters on the scaling factor to q = 1. We obtain Tc = 169 K, Pc = 136 MPa, and ρc = 1017 kg/m3. In Table I, we observe that these values of Tc and Pc are smaller but still close to those reported for the TIP4P/2005 model in the literature, while ρc is in the range of the reported critical densities.33–35,87 We expect that the former differences result because the q range of the extrapolation is nearly twice as large as that of the analysis and the q dependence is not necessarily linear in a broader parameter range.

TABLE I.

Critical parameters of the LLCP reported in different studies on the TIP4P/2005 water model.

StudyTc (K)Pc (MPa)ρc (kg/m3)
Abascal and Vega33  193 135 1012 
Sumi and Sekino34  182 158–162 1020 
Yagasaki et al.86  185 … 1020 
Singh et al.35  182 170 1017 
Handle and Sciortino87  175 175 997 
Debenedetti et al.37  177 175 … 
Present extrapolation 169 136 1017 
StudyTc (K)Pc (MPa)ρc (kg/m3)
Abascal and Vega33  193 135 1012 
Sumi and Sekino34  182 158–162 1020 
Yagasaki et al.86  185 … 1020 
Singh et al.35  182 170 1017 
Handle and Sciortino87  175 175 997 
Debenedetti et al.37  177 175 … 
Present extrapolation 169 136 1017 

Having determined the LLCP of the charge scaled TIP4P/2005 water models, we move onto their structure and dynamics. Figures 3 and 4 show results for q = 0.86 along the isochore with a density close to the critical one, ρ = 925 kg/m3ρc.

FIG. 3.

Temperature-dependent local structure of the charged scaled TIP4P/2005 water model with q = 0.86 along the isochore with a density of ρ = 925 kg/m3ρc. The critical temperature amounts to Tc = 124 K. (a) Distribution of the LSI p(I) at the indicated temperatures. The dashed line separates I values typical of HD and LD local structures. (b) Temperature-dependent fractions of the LD local structures xLD. The fractions determined from the criteria N4 = 1 and I ≥ 0.0006 nm2 are shown together with that from the TSEOS analysis. The results (open symbols) obtained from the isochore with ρ = 925 kg/m3 are compared with data (solid symbols) along the isobar with P = −65 MPa, which is close to the critical pressure Pc = −55 MPa.

FIG. 3.

Temperature-dependent local structure of the charged scaled TIP4P/2005 water model with q = 0.86 along the isochore with a density of ρ = 925 kg/m3ρc. The critical temperature amounts to Tc = 124 K. (a) Distribution of the LSI p(I) at the indicated temperatures. The dashed line separates I values typical of HD and LD local structures. (b) Temperature-dependent fractions of the LD local structures xLD. The fractions determined from the criteria N4 = 1 and I ≥ 0.0006 nm2 are shown together with that from the TSEOS analysis. The results (open symbols) obtained from the isochore with ρ = 925 kg/m3 are compared with data (solid symbols) along the isobar with P = −65 MPa, which is close to the critical pressure Pc = −55 MPa.

Close modal
FIG. 4.

Temperature-dependent dynamics of the charged scaled TIP4P/2005 water model with q = 0.86 along the isochore with a density of ρ = 925 kg/m3ρc: (a) MSD r2(t) and (b) RCF F1(t). In both panels, results for temperatures from 158 K (red) to 116 K (blue) are shown. The critical temperature amounts to Tc = 124 K.

FIG. 4.

Temperature-dependent dynamics of the charged scaled TIP4P/2005 water model with q = 0.86 along the isochore with a density of ρ = 925 kg/m3ρc: (a) MSD r2(t) and (b) RCF F1(t). In both panels, results for temperatures from 158 K (red) to 116 K (blue) are shown. The critical temperature amounts to Tc = 124 K.

Close modal

In Fig. 3(a), we observe a weakly bimodal LSI distribution, in particular, at lower temperatures, where a peak at larger values I ≈ 0.0011 nm2 emerges and grows with respect to that at smaller values I ≈ 0.0002 nm2 upon cooling. Although the effect is less prominent than in an inherent structure analysis for the original TIP4P/2005 model,88 these findings indicate that two local structural motifs become preferred in the vicinity of the LLCP and that the HDL-like environments are transformed into LDL-like arrangements upon cooling. However, at none of the studied temperatures, we observe a decomposition into extended HDL and LDL phases for the studied cubic system.

To further analyze the interconversion between the HD and LD local structures, we determine their temperature-dependent fractions. For this purpose, we assume that a water-like molecule is in a LD local structure if N4 = 1 or I ≥ 0.0006 nm2, where the latter value is roughly midway between the p(I) peaks. In Fig. 3(b), we see that the LD fractions xLD obtained from the different structural identifiers and from the TSEOS analysis, respectively, gradually grow when the temperature is decreased along the studied isochore with ρρc. By contrast, a sharp rise in xLD occurs near Tc when cooling under isobaric conditions with PPc. In the supplementary material, we present the pressure dependence of xLD at characteristic temperatures. We find that the LD fraction gradually decreases when the pressure is increased through Pc well above Tc, while a sharp drop occurs near the LLCP. Remarkably, well below Tc, the coexisting HD:LD compositions result in re-entrant xLD(P) behaviors for all studied structural identifiers, also including d5. This effect was not reported in an analogous study of the original TIP4P/2005 water model,35 where it is computationally more expensive to explore state points below the LLCP, once more demonstrating the added value of our charge-scaling approach.

In Fig. 4, the MSD and RCF indicate a prominent slowdown of the translational and rotational dynamics when the temperature is decreased. Typical of viscous liquids, the temporary caging of the molecules by their neighbors leads to a plateau of both quantities at times intermediate between the regimes of vibrational motion and viscous flow. To characterize the translational and rotational dynamics, we determine the self-diffusion coefficients D and the mean correlation times τ from r2(t) and F1(t), respectively. Analogous analyses are performed for various isochores and isobars. For none of the studied cubic systems, we find evidence for a clear bimodality of water dynamics due to coexisting HDL-like and LDL-like regions with different molecular mobilities.

Figures 5 and 6 present the diffusion coefficients and correlation times of the charge scaled TIP4P/2005 water model with q = 0.86 obtained along various isochores and isobars. We see that the translational and rotational dynamics strongly slow down when the density is decreased in the vicinity of ρc, in harmony with the dynamical anomalies reported for water and other tetrahedrally ordered liquids such as silica.1,60 When cooling along isochores, we do not find evidence for changes in the temperature dependence in the vicinity of the LLCP. In particular, there is no FS crossover for the isochore with the density ρ = 925 kg/m3ρc in the studied temperature range. Rather, the diffusion coefficients and correlation times can be described by the Vogel–Fulcher–Tammann (VFT) equation, e.g.,

τ=τ0expBTT0,
(11)

where the VFT parameters τ0, B, and T0 can be described as second order polynomial functions of the density. On the other hand, we observe an intriguing temperature dependence, when we cool under isobaric conditions at pressures in the vicinity of Pc = 55 MPa. Here, isobaric conditions mean that we perform NPT simulations and apply a barostat to set a desired pressure P and to adjust the box volume V at each studied temperature T along the cooling path. Specifically, the temperature dependence of the diffusion coefficients D(T) at P = const. resembles that at V = const. at sufficiently high or low temperatures, while it is much stronger under isobaric conditions in the vicinity of the LLCP, resulting in a sigmoidal shape; see Fig. 5. Similar effects are found for the correlation times τ(T), as is demonstrated in the supplementary material. When comparing the results for various P values, we observe that the onset of the enhanced temperature dependence is consistent with the location of the Widom line. In harmony with our findings, the temperature dependence of the diffusion coefficients was found to be higher along isobars than along isochores near the LLCP,74 and the enhanced temperature dependence was reported to occur upon crossing the Widom line for the ST2 and TIP4P/2005 water models.4,89 Thus, the gradual and rapid interconversions between the HD and LD local structures observed above upon isochoric and isobaric cooling across Tc, respectively, are accompanied by corresponding changes in the translational and rotational dynamics. This is a first hint at a structural origin of the FS dynamical crossover for the studied charge scaled water models.

FIG. 5.

Temperature-dependent reciprocal self-diffusion coefficients for the charge scaled TIP4P/2005 water model with q = 0.86 in the vicinity of the LLCP: (a) 1/D along isochores. The crosses indicate the results along the isochore with the density ρ = 925 kg/m3ρc. The solid lines are VFT fits. (b) 1/D along isobars. For comparison, the results for the isochore with the lowest density ρ = 870 kg/m3 are included. The dashed lines are 1/D values at the indicated pressures P calculated from the observed behaviors along the isochores using the information from the TSEOS analysis; see the text for details. The dotted lines are extrapolations of these behaviors to lower temperatures outside the studied ranges. In both panels, the 1/D value at the LLCP is marked by a circled cross and the values along the Widom line are shown as the dashed-dotted line.

FIG. 5.

Temperature-dependent reciprocal self-diffusion coefficients for the charge scaled TIP4P/2005 water model with q = 0.86 in the vicinity of the LLCP: (a) 1/D along isochores. The crosses indicate the results along the isochore with the density ρ = 925 kg/m3ρc. The solid lines are VFT fits. (b) 1/D along isobars. For comparison, the results for the isochore with the lowest density ρ = 870 kg/m3 are included. The dashed lines are 1/D values at the indicated pressures P calculated from the observed behaviors along the isochores using the information from the TSEOS analysis; see the text for details. The dotted lines are extrapolations of these behaviors to lower temperatures outside the studied ranges. In both panels, the 1/D value at the LLCP is marked by a circled cross and the values along the Widom line are shown as the dashed-dotted line.

Close modal
FIG. 6.

Temperature-dependent parameters of the rotational dynamics in the charge scaled TIP4P/2005 water model with q = 0.86 along isochores as obtained from KWW fits of F1 correlation functions: (a) mean correlation times τ and (b) stretching parameters βkww. The crosses indicate the results for the isochore with the density ρ = 925 kg/m3ρc. The legends refer to both panels. In (a), VFT fits are shown as solid lines, the τ value at the LLCP is marked by a circled cross, and the values along the Widom line are shown as the dashed-dotted line. These results for cubic systems (solid symbols) are compared with those for an elongated system with q = 0.86 (open symbols). For the latter “phase-separated” system, spatially resolved τL values are determined in layers with average densities ρL, which match the densities ρ of the studied cubic systems. Identical symbols are used for data with matching ρL and ρ values. The black stars characterize the time scale τρ of the fluctuations of the layer-averaged density ρL as extracted from the correlation functions Cρ(t) for the elongated system.

FIG. 6.

Temperature-dependent parameters of the rotational dynamics in the charge scaled TIP4P/2005 water model with q = 0.86 along isochores as obtained from KWW fits of F1 correlation functions: (a) mean correlation times τ and (b) stretching parameters βkww. The crosses indicate the results for the isochore with the density ρ = 925 kg/m3ρc. The legends refer to both panels. In (a), VFT fits are shown as solid lines, the τ value at the LLCP is marked by a circled cross, and the values along the Widom line are shown as the dashed-dotted line. These results for cubic systems (solid symbols) are compared with those for an elongated system with q = 0.86 (open symbols). For the latter “phase-separated” system, spatially resolved τL values are determined in layers with average densities ρL, which match the densities ρ of the studied cubic systems. Identical symbols are used for data with matching ρL and ρ values. The black stars characterize the time scale τρ of the fluctuations of the layer-averaged density ρL as extracted from the correlation functions Cρ(t) for the elongated system.

Close modal

Next, we discuss the shape of the F1(t) decays for the q = 0.86 system. For this purpose, we depict the KWW stretching parameters βkww obtained along various isochores in Fig. 6. To compare the situation at isokinetic states, these data are not plotted as a function of the inverse temperature but of the correlation time. Along all isochores, the stretching parameter continuously decreases, and hence, the nonexponentiality of the decays increases when the motion slows down. In particular, we do not observe a minimum of βkww(τ) near the LLCP, which would result if prominent structural inhomogeneity owing to coexisting HD and LD local structures at these state points caused enhanced dynamical heterogeneity, as was reported in recent work.90 At a given correlation time τ, βkww is smaller for lower values of ρ than for higher ones. This observation may be unexpected because it implies that the dynamics is more heterogeneous when the tetrahedral order is, on average, higher at lower densities. However, the βkww(ρ) dependence may be non-monotonous and show a weak minimum at ρρc owing to a moderate dynamical inhomogeneity associated with the LLCP.

The knowledge about the density dependence of the VFT parameters τ0(ρ), B(ρ), and T0(ρ) from the above analysis can be exploited to obtain the correlation time τc and, analogously, the diffusion coefficients Dc at the LLCP of the charge scaled water models from interpolations. In Fig. 2, we see that the translational and rotational dynamics at the critical point are significantly faster for smaller scaling factors q. This means that the LLCP is accessible more easily to simulations when the water-like molecules are less polar, which is an important advantage of our charge-scaling approach and will also be exploited in the following spatially resolved analysis. In addition, we find that the q dependence of the diffusion coefficients at the LLCP can be mapped onto that of the correlation times using the relation a2 = 6Dcτc and setting the jump length a to 0.325 nm, slightly larger than the next neighbor distance. Thus, short-range and long-range dynamics are related in the same way for all studied water-like models at the LLCP, while this is not necessarily true well inside the metastable region.

The findings from the isochoric and isobaric cooling schemes can be compared in more detail using the information from the above TSEOS and VFT analyses. Specifically, the TSEOS results yield the temperature-dependent densities ρ required to keep the pressure constant at a desired value upon cooling. Then, the knowledge about the density dependence of the VFT parameters can be employed to calculate the dynamics along an isobar from that along the isochores. In Fig. 5, we see that the, thus, calculated diffusion coefficients well agree with the above discussed NPT simulation results. Altogether, the results imply that the sigmoidal temperature dependence along isobars, which is particularly apparent when extrapolating the TSEOS–VFT based predictions to lower temperatures, cannot be interpreted in terms of a FS transition in the strict sense, i.e., as a result of a transition from an exceptionally fragile to a truly strong liquid. Rather, the diffusion coefficients and correlation times show a high temperature dependence in the vicinity of the LLCP because there is a rapid transition from a HDL-like liquid with fast dynamics to a LDL-like liquid with slow dynamics. Accordingly, when the pressure P of the isobar is reduced and the distance to the LLCP increases, the inflection points of 1/D(1/T) and τ(1/T) shift to higher temperatures and the sigmoidal shape becomes less prominent; see Fig. 5.

Similarly, previous studies used two-state models of water dynamics and argued that the FS transition results from a rapid change between HDL and LDL phases with different water mobilities. Specifically, Shi et al.91 proposed that there is no true FS transition between liquids with largely different fragilities but rather HDL and LDL waters obey Arrhenius laws albeit with different activation energies and a rapid crossover from one Arrhenius law to the other occurs. Likewise, Caupin and co-workers92,93 described experimental and computational results for the temperature- and pressure-dependent viscosity of supercooled water by a weighted superposition of two liquid behaviors and conjectured that the very high temperature dependence near the Widom line cannot be ascribed to either of the liquid phases but to their rapid interconversion. However, unlike Tanaka and co-workers, they reported that the high-temperature water phase is a fragile liquid and, hence, does not show Arrhenius behavior.

To scrutinize the liquid–liquid coexistence expected in the metastable region of the phase diagram and investigate the properties of both phases, we move from cubic to elongated systems and perform simulations at and below the LLCP. In doing so, we exploit that, for a given system size, the energetically unfavorable interface between the phases is minimized in elongated simulation boxes when it aligns parallel to the smallest box surface, in our case, the xy plane, facilitating phase separation. These simulations are restricted to the scaling factors q = 0.86 and, in particular, q = 0.88, where the LLCP is clearly detectable and the dynamics in the vicinity of this state point is relatively fast; see Fig. 2. It is ensured in the supplementary material that the pressure P at a given state point is the same within the error margins in the simulations for the cubic and elongated boxes.

In Fig. 7, we show the spatial variation of the density and structure for the charge-scaled model with q = 0.88 and ρ = 940 kg/m3 at three exemplary times during a trajectory at T = 125 K < Tc. Specifically, we depict the normalized density ρL/ρ and the structural identifiers IL and N4,L along the elongated z direction, as obtained from averaging over layers, which are Δz = 1 nm thick and centered about z. Evidently, the system separates into two phases with significantly different properties. In particular, the layer-averaged density ρL is ∼10% lower in the center than at the edges of the simulation box, while the layer-averaged IL and N4,L parameters are higher in the former than the latter regions. The interface between both regions weakly changes during the course of the simulation. These findings indicate that the system decomposes into extended and stable LDL-like and HDL-like regions in the central and outer box regions. In the supplementary material, we present the results from analogous analyses at T = 132 K > Tc, which reveal that corresponding fluctuations of the structural properties also exist somewhat above the LLCP but their length and time scales are smaller so that a stable decomposition into extended LDL-like and HDL-like regions does not occur. We note in passing that, because of the periodic boundary conditions, the system is comprised of only one HDL-like region and one LDL-like region at 125 K so that the extension of these regions is artificially determined by the finite size of the simulation box.

FIG. 7.

Spatially resolved structural properties of the charge scaled TIP4P/2005 water model with q = 0.88 at ρ = 940 kg/m3 and T = 125 K for exemplary times during the trajectory: (a) ρL/ρ, (b) IL, and (c) N4,L. The results are obtained for an elongated system (Lx : Ly : Lz ≡ 1 : 1 : 3) comprising N = 2000 molecules. The properties are averaged over Δz = 1 nm thick layers, which are parallel to the xy plane and centered about the respective z value. A low-pass filter was applied to reduce the noise of the data.

FIG. 7.

Spatially resolved structural properties of the charge scaled TIP4P/2005 water model with q = 0.88 at ρ = 940 kg/m3 and T = 125 K for exemplary times during the trajectory: (a) ρL/ρ, (b) IL, and (c) N4,L. The results are obtained for an elongated system (Lx : Ly : Lz ≡ 1 : 1 : 3) comprising N = 2000 molecules. The properties are averaged over Δz = 1 nm thick layers, which are parallel to the xy plane and centered about the respective z value. A low-pass filter was applied to reduce the noise of the data.

Close modal

Figure 8 compares the probability distribution of the layer-averaged local density ρL and structural identifiers IL and N4,L for the charge-scaled model with q = 0.88 at various temperatures. We see that all distributions have a Gaussian shape at T > Tc ≈ 131 K, while they are bimodal well below the LLCP. We emphasize that the considered distributions do not characterize the structural diversity of the immediate local environments, but rather the variation of the average properties in Δz = 1 nm thick layers. Thus, for the observation of bimodal behaviors, it is not sufficient that HD and LD local structures exist but it is additionally required that the respective structural identity is largely preserved throughout the respective layer. Therefore, the observed crossovers from unimodal to bimodal behaviors indicate in our case that the length scale of structural fluctuations becomes larger than the thickness of the layers, and hence, the extension of HDL-like and LDL-like regions reaches the nanometer regime, in harmony with the above findings from the visual inspection of Fig. 7.

FIG. 8.

Probability distributions of layer-averaged properties for the charge scaled TIP4P/2005 water model with q = 0.88 and ρ = 940 kg/m3 at the indicated temperatures: (a) ρL/ρ, (b) IL, and (c) N4,L. The results are obtained from an elongated system (Lx : Ly : Lz ≡ 1 : 1 : 3) comprising N = 2000 molecules. The critical temperature amounts to Tc ≈ 131 K. In the case of N4,L, a moving average was applied to reduce the noise caused by the fact that N4 can take only two discrete values.

FIG. 8.

Probability distributions of layer-averaged properties for the charge scaled TIP4P/2005 water model with q = 0.88 and ρ = 940 kg/m3 at the indicated temperatures: (a) ρL/ρ, (b) IL, and (c) N4,L. The results are obtained from an elongated system (Lx : Ly : Lz ≡ 1 : 1 : 3) comprising N = 2000 molecules. The critical temperature amounts to Tc ≈ 131 K. In the case of N4,L, a moving average was applied to reduce the noise caused by the fact that N4 can take only two discrete values.

Close modal

Next, we ascertain the relation between structure and dynamics for the charge scaled water models in more detail. To achieve this objective, one may tag the molecules with a structural identifier, e.g., the LSI, and characterize the dynamical behaviors of the, thus, defined subsets separately. If the used structural identifier changes rapidly in space and time, this approach will, however, not yield the true dynamics in HDL-like and LDL-like environments. Specifically, the occurring structural fluctuations affect the observed dynamical behaviors of the subsets, e.g., their correlation functions, to a degree, which depends on the relative length and time scales on which the structural identifier and the molecular dynamics are probed. Therefore, we exploit that extended and stable HDL-like and LDL-like regions exist in the elongated systems at lower temperatures. Moreover, we use layer-averaged quantities for tagging the molecules to ensure that the, thus, specified structural environments differ over sufficiently large regions, and we employ a local probe of the molecular dynamics, explicitly, the RCF F1(t), to minimize the changes in the structural environments while assessing the respective dynamics. Therefore, we expect that our findings are characteristic for the rotational motions in regions with clearly distinguishable structures.

Figure 9 presents the results for the elongated q = 0.88 system at T = 125 K < Tc. Comparing F1(t) for subsets of molecules characterized by different layer-averaged N4,L values, it is evident that the loss of correlation shifts to longer times when N4,L and, hence, the tetrahedral order increases. To quantify the dependence, the corresponding correlation times of the subsets τL are depicted as a function of N4,L. We find that the rotational motion is nearly two orders of magnitude slower in the LDL-like regions with N4,L > 0.7 than that in the HDL-like regions with N4,L < 0.3. Consistently, a similar slowdown is observed when the layer-averaged density ρL decreases from HDL-like to LDL-like values. As a consequence of this substantial dependence of the dynamics on the structure, the ensemble averaged RCF shows a two-step signature, which is, however, smeared out to some extent because the interfacial regions lead to deviations from a perfect bimodality of the N4 distribution; see Fig. 8. This bimodal nature is quantified in the supplementary material. Comparing these results for F1(t) with the decay of CN4(t), we find that the molecular orientation changes faster than the layer-averaged structural identifier. Hence, the correlation times at low N4,L (high ρL) and high N4,L (low ρL), indeed, reflect the dynamical behaviors in extended and stable LDL-like and HDL-like regions, respectively, consistent with the observation that the dependence of τL on N4,L or ρL levels off in these regions.

FIG. 9.

Relation between structure and dynamics for the charge scaled TIP4P/2005 water model with q = 0.88 at ρ = 940 kg/m3 and T = 125 K: (a) RCF F1(t) for molecules in regions characterized by various layer-averaged structural identifiers N4,L at the time origin t0. The legend specifies the midpoint of the N4,L bins with a width of 0.05. The colored solid lines are KWW fits. The black solid line is the RCF F1(t) of the whole ensemble. The black dashed line is the correlation function CN4(t) of the layer-averaged structural identifier N4,L; see Eq. (10). (b) Dependence of the corresponding correlation times τL on the layer-averaged properties N4,L (blue circles) and ρL/ρ (red squares).

FIG. 9.

Relation between structure and dynamics for the charge scaled TIP4P/2005 water model with q = 0.88 at ρ = 940 kg/m3 and T = 125 K: (a) RCF F1(t) for molecules in regions characterized by various layer-averaged structural identifiers N4,L at the time origin t0. The legend specifies the midpoint of the N4,L bins with a width of 0.05. The colored solid lines are KWW fits. The black solid line is the RCF F1(t) of the whole ensemble. The black dashed line is the correlation function CN4(t) of the layer-averaged structural identifier N4,L; see Eq. (10). (b) Dependence of the corresponding correlation times τL on the layer-averaged properties N4,L (blue circles) and ρL/ρ (red squares).

Close modal

Finally, we investigate the temperature dependence of the rotational dynamics in HDL-like and LDL-like regions. Figure 10 shows the rotational correlation times τL in layers with various average structural identifiers N4,L in an Arrhenius diagram. Independent of the value of N4,L, we observe weak deviations from Arrhenius behavior across the critical temperature Tc, whereby the temperature dependence seems to be higher for larger N4,L values. However, the time scale τN4 of the N4,L fluctuations becomes comparable to the correlation times τL for N4,L > 0.7 when increasing the temperature toward Tc. This finding implies that the τL values for large N4,L do no longer reflect the true LDL-related dynamics near the LLCP, because the correlation functions of these subsets decay prematurely when LDL-like environments are transformed into HDL-like environments with faster dynamics. In other words, the results for the LDL-related dynamics in the elongated system with N = 2000 are reliable only at sufficiently low temperatures, e.g., T = 125 K used in the above analysis.

FIG. 10.

Temperature-dependent correlation times τL of rotational dynamics in regions with the indicated layer-averaged structural identifiers N4,L. Solid and open symbols are results for elongated q = 0.88 systems with N = 2000 and N = 12 000 molecules and aspect ratios of 1:1:3 and 1:1:6, respectively. For comparison, the time scales τN4 of the N4,L fluctuations, as obtained from KWW fits of CN4(t), are shown as lines: (dashed) N = 2000 and (dotted) N = 12000. The arrow marks the critical temperature Tc ≈ 131 K. Reliable results for N4,L = 0.1 and N4,L = 0.8 can be obtained only well below the LLCP because these values are rarely found at higher temperatures, resulting in poor statistics.

FIG. 10.

Temperature-dependent correlation times τL of rotational dynamics in regions with the indicated layer-averaged structural identifiers N4,L. Solid and open symbols are results for elongated q = 0.88 systems with N = 2000 and N = 12 000 molecules and aspect ratios of 1:1:3 and 1:1:6, respectively. For comparison, the time scales τN4 of the N4,L fluctuations, as obtained from KWW fits of CN4(t), are shown as lines: (dashed) N = 2000 and (dotted) N = 12000. The arrow marks the critical temperature Tc ≈ 131 K. Reliable results for N4,L = 0.1 and N4,L = 0.8 can be obtained only well below the LLCP because these values are rarely found at higher temperatures, resulting in poor statistics.

Close modal

Therefore, we increase the system size to N = 12 000 and the aspect ratio to 1:1:6 to further reduce the size and relevance of a potential HDL–LDL interface and, hence, to stabilize the phase separation. In Fig. 10, we see that the fluctuations of the layer-averaged N4,L values are, indeed, slower for the larger than the smaller system. Moreover, the correlation times τL for high N4,L values mildly increase, while those for small N4,L remain essentially unchanged. On the other hand, long simulation times prevent us from exploring the temperature range well below the LLCP for N = 12 000. Thus, a reliable determination of the exact temperature dependence of the rotational dynamics in LDL-like regions is possible neither for the small nor for the large system. Still, the combined results for these systems indicate that distinct changes of LDL-like dynamics at Tc do not occur.

The correlation times τL extracted for q = 0.86 in regions with different layer-averaged densities ρL from the elongated N = 2000 system can be compared with those obtained along the isochores with the corresponding densities ρ from the cubic N = 2000 system. In Fig. 6, we observe that the temperature-dependent τL from the spatially resolved analysis of the “phase-separated” elongated system agrees with the correlation times τ extracted from the isochores for high densities ρ = ρL. This agreement confirms that these time constants characterize the dynamics in defined HDL-like regions. By contrast, the correlation times τL from the former analyses are smaller than the τ values from the latter for low densities ρ = ρL < ρc at temperatures above and near Tc. Again, this discrepancy can be attributed to the relatively fast HDL–LDL fluctuations in the elongated system at these temperatures. Accordingly, the τL values for low ρL tend to show a higher temperature dependence and approach the extrapolations of the isochore data with the corresponding densities ρ when the temperature is decreased well below the LLCP.

We performed MD simulations to ascertain the thermodynamic, structural, and dynamical properties of various water-like models in the supercooled temperature range. Specifically, we started from the well established TIP4P/2005 model and systematically reduced the magnitude of the partial charges of the atoms and, thus, the strength of the hydrogen bonds and the degree of the tetrahedral order.61 This charge-scaling approach allowed us to investigate the role of local order for liquid polyamorphism and to determine the consequences of this interplay for the dynamical behavior. This knowledge is of use for all liquids with spacious low-energy and low-entropy structural motifs and, in particular, for water, for which a unique model is still lacking.

For such liquids with a spacious local order, the existence of a LLPT and a related LLCP was proposed to cause anomalous properties but is still the subject of a controversial scientific debate. Here, we exploited that charge scaling allows us to shift these phenomena to state points, which are conveniently accessible in state-of-the-art computer simulations. Similar results were reported for the WAC model of silica.60 For all studied charge-scaling factors, q = 0.86–0.91, we observed isochore crossing in the P–T diagram indicative of a LLPT. Despite its mean-field nature, the TSEOS formalism well describes the isochore data and enables the determination of the LLCP parameters for the water-like models. These results indicate that a LLCP exists for a significant range of interaction parameters. The finding that the phenomenon is not restricted to very specific particle interactions suggests that it can also occur for various real liquids, including water.

The TSEOS analysis showed that the critical temperature Tc, pressure Pc, and density ρc decrease when the magnitude of the partial charges is reduced, leading even to negative values of Pc at q < 0.9. While the critical parameters showed an approximately linear q dependence in the studied range q = 0.86–0.91, we expect a mildly nonlinear behavior in a broader range, because extrapolation to q = 1.00 underestimates the Tc and Pc values known from studies of the regular TIP4P/2005 model.33–35,37,86,87 When the charge-scaling factor and the hydrogen-bond strength are decreased, the isochores flatten out in the P–T diagram, suggesting that the liquid polyamorphism vanishes not too far below the covered q range. Studies on silica and silicon reported that the LLCP disappears below the liquid–vapor spinodal.31,60 Because the critical pressure takes strongly negative values for q = 0.86, one may speculate that the liquid–vapor spinodal is also close for water-like models with substantially reduced partial charges.

As for structure and dynamics, we found prominent changes when cooling in the vicinity of the LLCP. Specifically, distributions of local structure identifiers become bimodal, indicating a decomposition into HDL-like and LDL-like environments. Furthermore, the temperature dependence of the diffusion coefficients D and correlation times τ differ significantly under isochoric and isobaric conditions, consistent with previous results for the ST2 water model.74 Along the isochores, we observed moderate deviations from Arrhenius behavior, and the temperature dependence does not change even when crossing the LLCP for ρρc. Along the isobars, by contrast, the temperature dependence of D and τ is much higher near the LLCP than at TTc and TTc. These distinct differences between isochoric and isobaric conditions are confirmed when we use the temperature-dependent correlation times along the isochores to successfully predict those along the isobars. In addition, we observed that the dynamics of the water-like models slows down when the density or the pressure is decreased in the vicinity of the LLCP, consistent with the dynamical anomaly of water.

To ascertain the origin of the remarkable dynamical behaviors, we analyzed the respective dynamics in regions rich in HD or LD local structures, which span over a few nanometers and are stable for a few microseconds. For this purpose, we exploited that a decomposition into such regions occurs at sufficiently low temperatures in elongated systems, where the interface between these structurally different regions can be minimized along the smallest (xy) cross section of the simulation box. We characterized the variation of the structural and dynamical properties along the elongated z dimension by averaging over Δz = 1 nm thick layers centered about the respective z value. In this way, we ensured that, except for the highest studied temperatures, the structural differences persist on the length and time scales on which the dynamical properties are probed. Contrarily, in previous approaches to the structure–dynamics relations of water,90,94 different subensembles were discriminated based on local structural identifiers, which fluctuate more rapidly in space and time so that partially averaged dynamical behaviors may be observed.

This spatially resolved analysis was restricted to the systems with the small scaling factors of q = 0.86 and q = 0.88, which exhibit relatively fast dynamics at the LLCP. It revealed that the rotational motion of the water-like molecules is about two orders of magnitude faster in HDL-like regions than that in LDL-like regions. Consistently, experimental work reported that the rotational correlation times differ by about this amount when heating high-density and low-density amorphous ice phases through their glass transition temperatures Tg.17 Considering that LDL-like configurations with slow dynamics are transformed into HDL-like ones with fast dynamics when the pressure is increased, these findings yield an explanation for the anomalous pressure dependence of the diffusion coefficients and correlation times. For neither the LDL-like nor the HDL-like regions, we found a change in the temperature dependence of the dynamics near the LLCP.

These findings yield insights into the nature of dynamical crossovers observed for water and other tetrahedral liquids near the LLCP upon cooling along isobars. This observation was often interpreted in terms of a FS transition related to a LLPT. Explicitly, one may argue that the HDL and LDL are very fragile and truly strong liquids, respectively, so that a FS transition occurs when the LLPT or the Widom line are crossed upon cooling. At variance with this conjecture, the present findings for water-like systems indicated that both HDL-like and LDL-like regions show only weak deviations from Arrhenius behavior. Thus, the dynamic crossover is not due to a transition between liquids with highly different fragilities. In other words, the HDL-like and LDL-like dynamic properties are observed only at a sufficient distance to the LLCP, and the effects in the vicinity of the LLCP cannot be attributed to either of the two liquid phases. Rather, the high temperature dependence near the LLCP is caused by the interconversion between fast HDL-like structures and slow LDL-like structures, i.e., it results from the crossover between the largely different D(T) or τ(T) curves of the respective regions. As a consequence, in the one-phase regime of the phase diagram, the temperature dependence along isobars becomes higher when the interconversion occurs more rapidly upon approaching Pc. While the dynamics in the NVT and NPT ensembles agree at the same state point, a clear dynamical crossover is not observed along the isochores where rapid HDL–LDL interconversions do not occur. This scenario has certain similarities with two-state models of water dynamics by Tanaka and co-workers91 and Caupin and co-workers,92,93 who proposed that the FS transition involves a crossover between water phases with dynamics on very different time scales.

See the supplementary material for details on the TSEOS analysis, isochores for further charge-scaling factors, temperature-dependent densities along isobars at atmospheric pressure, pressure-dependent fractions of the low-density local structures, additional data for the temperature-dependent dynamics along isochores and isobars, a comparison of pressures for cubic and elongated simulation boxes, and spatially resolved analyses above Tc.

The authors thank the Deutsche Forschungsgemeinschaft (DFG) for funding through Forschungsgruppe FOR 1583 (Project no. Vo 905/9-2).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
P. G.
Debenedetti
,
J. Phys.: Condens. Matter
15
,
R1669
(
2003
).
2.
R. J.
Speedy
and
C. A.
Angell
,
J. Chem. Phys.
65
,
851
(
1976
).
3.
C.
Goy
,
M. A. C.
Potenza
,
S.
Dedera
,
M.
Tomut
,
E.
Guillerm
,
A.
Kalinin
,
K.-O.
Voss
,
A.
Schottelius
,
N.
Petridis
,
A.
Prosvetov
,
G.
Tejeda
,
J. M.
Fernández
,
C.
Trautmann
,
F.
Caupin
,
U.
Glasmacher
, and
R. E.
Grisenti
,
Phys. Rev. Lett.
120
,
015501
(
2018
).
4.
P. H.
Poole
,
F.
Sciortino
,
U.
Essmann
, and
H. E.
Stanley
,
Nature
360
,
324
(
1992
).
5.
O.
Mishima
and
H. E.
Stanley
,
Nature
396
,
329
(
1998
).
6.
H. E.
Stanley
,
S. V.
Buldyrev
,
M.
Canpolat
,
O.
Mishima
,
M. R.
Sadr-Lahijany
,
A.
Scala
, and
F. W.
Starr
,
Phys. Chem. Chem. Phys.
2
,
1551
(
2000
).
7.
O.
Mishima
and
Y.
Suzuki
,
Nature
419
,
599
(
2002
).
8.
K.
Amann-Winkel
,
R.
Böhmer
,
F.
Fujara
,
C.
Gainaru
,
B.
Geil
, and
T.
Loerting
,
Rev. Mod. Phys.
88
,
011002
(
2016
).
9.
P.
Gallo
,
K.
Amann-Winkel
,
C. A.
Angell
,
M. A.
Anisimov
,
F.
Caupin
,
C.
Chakravarty
,
E.
Lascaris
,
T.
Loerting
,
A. Z.
Panagiotopoulos
,
J.
Russo
,
J. A.
Sellberg
,
H. E.
Stanley
,
H.
Tanaka
,
C.
Vega
,
L.
Xu
, and
L. G. M.
Pettersson
,
Chem. Rev.
116
,
7463
(
2016
).
10.
S.
Lemke
,
P. H.
Handle
,
L. J.
Plaga
,
J. N.
Stern
,
M.
Seidl
,
V.
Fuentes-Landete
,
K.
Amann-Winkel
,
K. W.
Köster
,
C.
Gainaru
,
T.
Loerting
, and
R.
Böhmer
,
J. Chem. Phys.
147
,
034506
(
2017
).
11.
F.
Caupin
,
J. Non-Cryst. Solids
407
,
441
(
2015
).
12.
K. H.
Kim
,
A.
Späh
,
H.
Pathak
,
F.
Perakis
,
D.
Mariedahl
,
K.
Amann-Winkel
,
J. A.
Sellberg
,
J. H.
Lee
,
S.
Kim
,
J.
Park
,
K. H.
Nam
,
T.
Katayama
, and
A.
Nilsson
,
Science
358
,
1589
(
2017
).
13.
K. H.
Kim
,
K.
Amann-Winkel
,
N.
Giovambattista
,
A.
Späh
,
F.
Perakis
,
H.
Pathak
,
M. L.
Parada
,
C.
Yang
,
D.
Mariedahl
,
T.
Eklund
 et al,
Science
370
,
978
(
2020
).
14.
L.
Kringle
,
W. A.
Thornley
,
B. D.
Kay
, and
G. A.
Kimmel
,
Science
369
,
1490
(
2020
).
15.
16.
J.
Swenson
,
Phys. Chem. Chem. Phys.
20
,
30095
(
2018
).
17.
K.
Amann-Winkel
,
C.
Gainaru
,
P. H.
Handle
,
M.
Seidl
,
H.
Nelson
,
R.
Böhmer
, and
T.
Loerting
,
Proc. Natl. Acad. Sci. U. S. A.
110
,
17720
(
2013
).
18.
H.
Tanaka
,
J. Chem. Phys.
153
,
130901
(
2020
).
19.
M.
Togaya
,
Phys. Rev. Lett.
79
,
2474
(
1997
).
20.
Y.
Katayama
,
T.
Mizutani
,
W.
Utsumi
,
O.
Shimomura
,
M.
Yamakata
, and
K.-i.
Funakoshi
,
Nature
403
,
170
(
2000
).
21.
Y.
Katayama
,
Y.
Inamura
,
T.
Mizutani
,
M.
Yamakata
,
W.
Utsumi
, and
O.
Shimomura
,
Science
306
,
848
(
2004
).
22.
L.
Henry
,
M.
Mezouar
,
G.
Garbarino
,
D.
Sifré
,
G.
Weck
, and
F.
Datchi
,
Nature
584
,
382
(
2020
).
23.
F.
Smallenburg
,
L.
Filion
, and
F.
Sciortino
,
Nat. Phys.
10
,
653
(
2014
).
24.
L.
Xu
,
S. V.
Buldyrev
,
C. A.
Angell
, and
H. E.
Stanley
,
Phys. Rev. E
74
,
031108
(
2006
).
25.
A.
Ben-Naim
,
J. Phys. Chem.
128
,
024505
(
2008
).
26.
L.
Heckmann
and
B.
Drossel
,
J. Chem. Phys.
138
,
234503
(
2013
).
27.
V. V.
Vasisht
,
S.
Saw
, and
S.
Sastry
,
Nat. Phys.
7
,
549
(
2011
).
28.
E.
Lascaris
,
M.
Hemmati
,
S. V.
Buldyrev
,
H. E.
Stanley
, and
C. A.
Angell
,
J. Chem. Phys.
140
,
224502
(
2014
).
29.
R.
Chen
,
E.
Lascaris
, and
J. C.
Palmer
,
J. Chem. Phys.
146
,
234503
(
2017
).
30.
D.
Dhabal
,
C.
Chakravarty
,
V.
Molinero
, and
H. K.
Kashyap
,
J. Chem. Phys.
145
,
214502
(
2016
).
31.
G.
Zhao
,
Y.
Yu
,
J.
Yan
,
M.
Ding
,
X.
Zhao
, and
H.
Wang
,
Phys. Rev. B
93
,
140203
(
2016
).
32.
J. C.
Palmer
,
P. H.
Poole
,
F.
Sciortino
, and
P. G.
Debenedetti
,
Chem. Rev.
118
,
9129
(
2018
).
33.
J. L. F.
Abascal
and
C.
Vega
,
J. Chem. Phys.
133
,
234502
(
2010
).
34.
T.
Sumi
and
H.
Sekino
,
RSC Adv.
3
,
12743
(
2013
).
35.
R. S.
Singh
,
J. W.
Biddle
,
P. G.
Debenedetti
, and
M. A.
Anisimov
,
J. Chem. Phys.
144
,
144504
(
2016
).
36.
J.
Russo
and
H.
Tanaka
,
Nat. Commun.
5
,
3556
(
2014
).
37.
P. G.
Debenedetti
,
F.
Sciortino
, and
G. H.
Zerze
,
Science
369
,
289
(
2020
).
38.
Y.
Li
,
J.
Li
, and
F.
Wang
,
Proc. Natl. Acad. Sci. U. S. A.
110
,
12209
(
2013
).
39.
Y.
Ni
and
J. L.
Skinner
,
J. Chem. Phys.
144
,
214501
(
2016
).
40.
Y.
Liu
,
J. C.
Palmer
,
A. Z.
Panagiotopoulos
, and
P. G.
Debenedetti
,
J. Chem. Phys.
137
,
214505
(
2012
).
41.
V.
Holten
,
J. C.
Palmer
,
P. H.
Poole
,
P. G.
Debenedetti
, and
M. A.
Anisimov
,
J. Chem. Phys.
140
,
104502
(
2014
).
42.
V.
Holten
,
D. T.
Limmer
,
V.
Molinero
, and
M. A.
Anisimov
,
J. Chem. Phys.
138
,
174501
(
2013
).
43.
F.
Sciortino
,
E. L.
Nave
, and
P.
Tartaglia
,
Phys. Rev. Lett.
91
,
155701
(
2003
).
44.
D. T.
Limmer
and
D.
Chandler
,
J. Chem. Phys.
138
,
214504
(
2013
).
45.
S. D.
Overduin
and
G. N.
Patey
,
J. Chem. Phys.
143
,
094504
(
2015
).
46.
J. C.
Palmer
,
F.
Martelli
,
Y.
Liu
,
R.
Car
,
A. Z.
Panagiotopoulos
, and
P. G.
Debenedetti
,
Nature
510
,
385
(
2014
).
47.
L. G. M.
Pettersson
and
A.
Nilsson
,
J. Non-Cryst. Solids
407
,
399
(
2015
).
48.
V.
Holten
and
M.
Anisimov
,
Sci. Rep.
2
,
713
(
2012
).
49.
C. E.
Bertrand
and
M. A.
Anisimov
,
J. Phys. Chem. B
115
,
14099
(
2011
).
50.
51.
J. W.
Biddle
,
R. S.
Singh
,
E. M.
Sparano
,
F.
Ricci
,
M. A.
González
,
C.
Valeriani
,
J. L. F.
Abascal
,
P. G.
Debenedetti
,
M. A.
Anisimov
, and
F.
Caupin
,
J. Chem. Phys.
146
,
034502
(
2017
).
52.
M. A.
Anisimov
,
M.
Duška
,
F.
Caupin
,
L. E.
Amrhein
,
A.
Rosenbaum
, and
R. J.
Sadus
,
Phys. Rev. X
8
,
011004
(
2018
).
53.
P.
Mausbach
,
H.-O.
May
, and
G.
Ruppeiner
,
J. Chem. Phys.
151
,
064503
(
2019
).
54.
F.
Caupin
and
M. A.
Anisimov
,
J. Chem. Phys.
151
,
034503
(
2019
).
55.
V.
Holten
,
J. V.
Sengers
, and
M. A.
Anisimov
,
J. Phys. Chem. Ref. Data
43
,
043101
(
2014
).
56.
F. H.
Stillinger
and
T. A.
Weber
,
Phys. Rev. B
31
,
5262
(
1985
).
57.
V.
Molinero
,
S.
Sastry
, and
C. A.
Angell
,
Phys. Rev. Lett.
97
,
075701
(
2006
).
58.
C. A.
Angell
and
V.
Kapko
,
J. Stat. Mech.: Theory Exp.
2016
,
094004
.
59.
F.
Smallenburg
and
F.
Sciortino
,
Phys. Rev. Lett.
115
,
015701
(
2015
).
60.
E.
Lascaris
,
M.
Hemmati
,
S. V.
Buldyrev
,
H. E.
Stanley
, and
C. A.
Angell
,
J. Chem. Phys.
142
,
104506
(
2015
).
61.
R.
Horstmann
and
M.
Vogel
,
J. Chem. Phys.
147
,
034505
(
2017
).
62.
J.
Geske
,
M.
Harrach
,
L.
Heckmann
,
R.
Horstmann
,
F.
Klameth
,
N.
Müller
,
E.
Pafong
,
T.
Wohlfromm
,
B.
Drossel
, and
M.
Vogel
,
Z. Phys. Chem.
232
,
1187
(
2018
).
63.
E.
Pafong Sanjon
,
B.
Drossel
, and
M.
Vogel
,
J. Chem. Phys.
148
,
104506
(
2018
).
64.
T.
Pal
and
M.
Vogel
,
J. Chem. Phys.
150
,
124501
(
2019
).
65.
J. J.
Shephard
and
C. G.
Salzmann
,
J. Phys. Chem. Lett.
7
,
2281
(
2016
).
66.
J.
Engstler
and
N.
Giovambattista
,
J. Chem. Phys.
147
,
074505
(
2017
).
67.
L.
Liu
,
S.-H.
Chen
,
A.
Faraone
,
C.-W.
Yen
, and
C.-Y.
Mou
,
Phys. Rev. Lett.
95
,
117802
(
2005
).
68.
K.
Ito
,
C. T.
Moynihan
, and
C. A.
Angell
,
Nature
398
,
492
(
1999
).
69.
F. W.
Starr
,
C. A.
Angell
, and
H. E.
Stanley
,
Physica A
323
,
51
(
2003
).
70.
M. D.
Marzio
,
G.
Camisasca
,
M.
Rovere
, and
P.
Gallo
,
J. Chem. Phys.
146
,
084502
(
2017
).
71.
I.
Saika-Voivod
,
P. H.
Poole
, and
F.
Sciortino
,
Nature
412
,
514
(
2001
).
72.
S.
Cerveny
,
F.
Mallamace
,
J.
Swenson
,
M.
Vogel
, and
L.
Xu
,
Chem. Rev.
116
,
7608
(
2016
).
73.
J.
Swenson
and
S.
Cerveny
,
J. Phys.: Condens. Matter
27
,
033102
(
2015
).
74.
P. H.
Poole
,
S. R.
Becker
,
F.
Sciortino
, and
F. W.
Starr
,
J. Phys. Chem. B
115
,
14176
(
2011
).
75.
J. L. F.
Abascal
and
C.
Vega
,
J. Chem. Phys.
123
,
234505
(
2005
).
76.
D.
Van Der Spoel
,
E.
Lindahl
,
B.
Hess
,
G.
Groenhof
,
A. E.
Mark
, and
H. J. C.
Berendsen
,
J. Comput. Chem.
26
,
1701
(
2005
).
77.
B.
Hess
,
C.
Kutzner
,
D.
van der Spoel
, and
E.
Lindahl
,
J. Chem. Theory Comput.
4
,
435
(
2008
).
78.
T.
Darden
,
D.
York
, and
L.
Pedersen
,
J. Chem. Phys.
98
,
10089
(
1993
).
79.
S.
Nosé
,
J. Chem. Phys.
81
,
511
(
1984
).
80.
W. G.
Hoover
,
Phys. Rev. A
31
,
1695
(
1985
).
81.
M.
Parrinello
and
A.
Rahman
,
J. Appl. Phys.
52
,
7182
(
1981
).
82.
J.
Guo
and
J. C.
Palmer
,
Phys. Chem. Chem. Phys.
20
,
25195
(
2018
).
83.
E.
Shiratani
and
M.
Sasai
,
J. Chem. Phys.
104
,
7671
(
1996
).
84.

Because ergodicity was not fully reached for the layer-averaged quantities at the lowest studied temperatures, it proved useful to calculate ξL2 as ⟨ξL(t + t0)⟩⟨ξL(t0)⟩ and to normalize the correlation functions by C(t = 0). In the ergodic limit, this is equivalent to Eq. (10).

85.
A.
Skibinsky
,
S. V.
Buldyrev
,
G.
Franzese
,
G.
Malescio
, and
H. E.
Stanley
,
Phys. Rev. E
69
,
061206
(
2004
).
86.
T.
Yagasaki
,
M.
Matsumoto
, and
H.
Tanaka
,
Phys. Rev. E
89
,
020301
(
2014
).
87.
P. H.
Handle
and
F.
Sciortino
,
J. Chem. Phys.
148
,
134505
(
2018
).
88.
K. T.
Wikfeldt
,
A.
Nilsson
, and
L. G. M.
Pettersson
,
Phys. Chem. Chem. Phys.
13
,
19918
(
2011
).
89.
M.
De Marzio
,
G.
Camisasca
,
M.
Rovere
, and
P.
Gallo
,
J. Chem. Phys.
144
,
074503
(
2016
).
90.
R.
Shi
and
H.
Tanaka
,
J. Chem. Phys.
148
,
124503
(
2018
).
91.
R.
Shi
,
J.
Russo
, and
H.
Tanaka
,
Proc. Natl. Acad. Sci. U. S. A.
115
,
9444
(
2018
).
92.
L. P.
Singh
,
B.
Issenmann
, and
F.
Caupin
,
Proc. Natl. Acad. Sci. U. S. A.
114
,
4312
(
2017
).
93.
P.
Montero de Hijes
,
E.
Sanz
,
L.
Joly
,
C.
Valeriani
, and
F.
Caupin
,
J. Chem. Phys.
149
,
094503
(
2018
).
94.
G.
Camisasca
,
N.
Galamba
,
K. T.
Wikfeldt
, and
L. G. M.
Pettersson
,
J. Chem. Phys.
150
,
224507
(
2019
).

Supplementary Material