Liquid-Liquid Phase Transition in Supercooled H2O and D2O: A Path-Integral Computer Simulation Study

We perform path-integral molecular dynamics (PIMD) and classical MD simulations of H2O and D2O using the q-TIP4P/F water model over a wide range of temperatures and pressures. The density ρ(T ), isothermal compressibility κT (T ), and self-diffusion coefficients D(T ) of H2O and D2O are in excellent agreement with available experimental data; the isobaric heat capacity CP (T ) obtained from PIMD and MD simulations agree qualitatively well with the experiments. Some of these thermodynamic properties exhibit anomalous maxima upon isobaric cooling, consistent with recent experiments and with the possibility that H2O and D2O exhibit a liquid-liquid critical point (LLCP) at low temperatures and positive pressures. The data from PIMD/MD for H2O and D2O can be fitted remarkably well using the Two-State-Equation-of-State (TSEOS). Using the TSEOS, we estimate that the LLCP for q-TIP4P/F H2O, from PIMD simulations, is located at Pc = 167±9 MPa, Tc = 159±6 K, and ρc = 1.02±0.01 g/cm . Isotope substitution effects are important; the LLCP location in q-TIP4P/F D2O is estimated to be Pc = 176± 4 MPa, Tc = 177± 2 K, and ρc = 1.13±0.01 g/cm . Interestingly, for the water model studied, differences in the LLCP location from PIMD and MD simulations suggest that nuclear quantum effects (i.e., atoms delocalization) play an important role in the thermodynamics of water around the LLCP (from the MD simulations of q-TIP4P/F water, Pc = 203±4 MPa, Tc = 175±2 K, and ρc = 1.03±0.01 g/cm ). Overall, our results strongly support the LLPT scenario to explain water anomalous behavior, independently of the fundamental differences between classical MD and PIMD techniques. The reported values of Tc for D2O and, particularly, H2O suggest that improved water models are needed for the study of supercooled water. ∗ ali.eltareb@brooklyn.cuny.edu, gustavo.lopez1@lehman.cuny.edu, ngiovambattista@brooklyn.cuny.edu


I. INTRODUCTION
temperature is overestimated (T c = 237 K, P c = 167 MPa, ρ c = 0.99 g/cm 3 ) [41]; while in the TIP4P/2005 and TIP4P/Ice water models it is underestimated (T c = 172 K, P c = 186 MPa, and ρ c = 1.03 g/cm 3 for TIP4P/2005; T c = 188 K, P c = 175 MPa, ρ c = 1.01 g/cm 3 for TIP4P/Ice) [18]. In all these cases, the LLCP pressure is overestimated by approximately 100 MPa [6]. The computer simulation studies that find a LLCP in water are based on classical models where the atoms delocalization due to nuclear quantum effects (NQE) are neglected. This can be troublesome because water is a light molecule and delocalization effects of its H atoms occur even at standard temperatures and pressures [45,46]. For example, the temperature of maximum density and the glass transition temperature (T g ) of H 2 O and D 2 O differ by 7 − 10 K, a clear sign of non-negligible NQE. Path-integral computer simulations of water-like models that have a LLCP clearly show that NQE can indeed shift the location of the LLCP as well as alter the shape and slope of the C P (T ) and κ T (T ) maxima lines [47,48]. Interestingly, while experiments in glassy water have estimated differences in the location of LLCP in H 2 O and D 2 O (∆T c ≈ 10 K, ∆P c ≈ 50 MPa) [49,50], the issue of isotope effects on the location of the LLCP has not been explored in computational and theoretical studies.
In this work we perform extensive path-integral molecular dynamics (PIMD) and classical MD simulations of light and heavy water using the q-TIP4P/F model and explore the corresponding phase diagram and thermodynamic/dynamical anomalous properties. One goal of this work is to determine the NQE (due to atoms delocalization) on the location of the LLCP, LLPT, and supercritical anomalous lines (such as maxima lines in ρ, C P , and κ T ) in q-TIP4P/F water (H 2 O). The second goal of this work is to study isotope substitution effects in water, i.e., whether PIMD simulations of H 2 O and D 2 O can reproduce the subtle differences in the phase diagram and anomalous properties of H 2 O and D 2 O observed in experiments. In a previous study, we performed PIMD simulations of q-TIP4P/F water at P = 0.1 MPa and showed that this model reproduces some signatures of the LLPT scenario, specifically, a maximum in κ T (T ) was found in H 2 O and D 2 O at T ≈ 230 − 235 K (P = 0.1) MPa [46]. Here, we extend our previous study to a wide range of temperatures and pressures in order to explore the possible existence of a LLCP in H 2 O and D 2 O. By combining the PIMD/MD results and the two-state equation of state (TSEOS) [14,22,23], we are able to identify a LLCP in both H 2 O and D 2 O. The TSEOS is based on the assumption that liquid water is a mixture of two interconvertible (liquid) states. The TSEOS has been shown to fit remarkably well the computer simulations results obtained from classical MD simulations of ST2 and TIP4P/2005 water as well as a water model based on DFT and machine learning techniques [14,[19][20][21]23]; it has also been applied to the case of real water [20,51].
Our paper is organized as follows. In Sec. II, we present the computer simulation details.
In Sec. III, we discuss the results from our PIMD and classical MD simulations of H 2 O using the q-TIP4P/F water model. The phase diagram of D 2 O is briefly discussed and compared with the phase diagram of H 2 O. A summary and discussions are included in Sec. IV.

II. SIMULATION METHOD
Our results are based on PIMD and classical MD simulations of a system composed of N = 512 water molecules in a cubic box with periodic boundary conditions. H 2 O and D 2 O molecules are represented using the non-rigid q-TIP4P/F model [52]. This model is based on the TIP4P/2005 model for water [53], commonly used in classical MD simulations.
The q-TIP4P/F water model was optimized for path integral computer simulations and has been shown to be able to reproduce remarkably well the properties of liquid water at P = 0.1 MPa [46,52]. Here, we perform PIMD and MD simulations at constant N , P , and T over a wide range of temperatures and pressures, 180 ≤ T ≤ 375 K and −250 ≤ P ≤ 500 MPa; see Fig. S1 of the Supplementary Information (SI). The temperature of the system is maintained constant using a stochastic (local) path integral Langevin equation (PILE) thermostat [54] while the pressure of the system is controlled by using a Monte Carlo Barostat (additional computational details can be found in Ref. [46]). In the PIMD simulations, the time step dt is set to 0.25 fs and the number of beads per ring-polymer/atom was set to n b = 32; in Ref. [46], it is shown that this value of n b is large enough to obtain well-converged dynamical, thermodynamic, and structural properties of q-TIP4P/F water at P = 0.1 MPa and T ≥ 210 K. Short-range (Lennard-Jones pair potential) interactions are calculated using a cutoff r c = 1.0 nm and long range electrostatic interactions are computed using the Particle Mesh Ewald (PME) method with the same cutoff r c . In the classical MD simulations, we employ a time step dt = 0.50 fs and set n b = 1. All PIMD and classical MD simulations are performed using the OpenMM software package (version 7.4.0) [55].
In all PIMD and classical MD simulations, the system is equilibrated for a time interval t eq followed by a production run of time length t prod . The values of t eq and t prod depend on the state point simulated. Typical simulation times for PIMD range from 2.5 ns to 100 ns.
Simulation times for classical MD simulations range from 2.5 ns to 2 − 3 µs. To confirm that the system reaches equilibrium, we monitor the mean-square displacement (MSD) of the water molecules in the system as a function of time and confirm that the PIMD and classical MD simulations satisfy the requirement that t eq , t prod > τ , where τ is the time it takes for the MSD of water molecules to reach 1 nm 2 .

III. RESULTS
The results are presented as follows. In Sec. III A, we show that the phase diagrams of H 2 O from MD and PIMD simulations are consistent with the existence of a LLPT and LLCP at low temperatures. Since a LLCP generates anomalous loci of maxima in C P and κ T , the behavior of C P (T ) and κ T (T ) are discussed in Sec. III B. The self-diffusion coefficient of  (similar values of ∆ρ hold for the case of MD simulations). It follows that both MD and PIMD simulations of q-TIP4P/F water predict the correct location (T and P) for the density maximum of water. In these and other cases, the computer simulation results can be fitted remarkably well using the TSEOS [14,[20][21][22][23]56]. In Figs. 1(c) and 1(d), we compare the ρ(T ) isobars obtained from the MD and PIMD simulations of q-TIP4P/F water with the corresponding fit using the TSEOS (a brief explanation of the methodology used to obtain the TSEOS can be found in Refs. [14,21]). The TSEOS is fitted using only the PIMD and MD data for 180 ≤ T ≤ 325 K and −50 ≤ P ≤ 350 MPa. As shown in Figs. 1(c) and 1(d), the TSEOS isobars are in excellent agreement with the simulation results over a majority of the state points simulated. Interestingly, the TSEOS also predicts a minimum in the ρ(T ) isobars of q-TIP4P/F water. While at the studied temperatures we do not observe density minima in the classical MD and PIMD simulations, density minima were reported at different pressures in TIP4P/2005 and ST2 water [15,42]. The optimized parameters for the TSEOS are given in Table S1 of the SI.
The TSEOS also provides a good estimation of the LLCP location. For example, in the case of TIP4P/2005 water, the TSEOS predicts that T c = 182 K, P c = 170 MPa, and ρ c = 1.02 g/cm 3 [14], which is in good agreement with recent MD simulations that were able to access the LLCP, T c = 172 K, P c = 186 MPa, and ρ c = 1.03 g/cm 3 [18]; similar results were found in a MD simulation study combined with the potential energy landscape theoretical approach [43]. Using the TSEOS, one can estimate the LLCP location (ρ c , T c , P c ). The values of (ρ c , T c , P c ) for the case of q-TIP4P/F H 2 O, based on classical MD and PIMD simulations, are given in Table I  Interestingly, recent studies based on water-like monoatomic model liquids that exhibit a LLCP, show that including NQE has the effects of lowering T c and increasing P c , while leaving ρ c unaffected [47,48]. We obtain the isothermal compressibility of q-TIP4P/F water by calculating the density fluctuations of the system [61], where ... indicates average over time and k B is the Boltzmann's constant. Fig. 3 and 100 MPa. This is an anomalous property that was originally predicted by the LLPT hypothesis scenario and later confirmed by experiments [5,17]. The experimental data from Ref. [5] is included in Fig. 3 This is fully consistent with the PIMD simulations results shown in Fig. 3(a). Specifically, as the pressure increases from P = 0.1 MPa to P = 100 MPa, the κ T -maximum shifts to lower temperatures and increases in height. The same behavior of κ T is found in classical MD simulations of water models that exhibit a LLCP [14,15,21,62].
We also calculate κ T (T ) using the TSEOS. The TSEOS provides an expression for the Gibbs free energy of the system from which the isothermal compressibility can easily be obtained, Next, we discuss the isobaric heat capacity, In our previous work (at P = 0.1 MPa) [46], the enthalpy was calculated directly from MD and PIMD simulations at selected temperatures and then, the values of H(T ) were fitted using a fourth-order polynomial. The resulting analytic expression for H(T ) was then used in Eq. 3 to calculate C P (T ). The use of a fourth-order polynomial in the fitting procedure is rather arbitrary. It captures the qualitative increase of C P (T ) upon cooling at low pressures but it may play a relevant role in identifying a C P -maximum, which is known to occur in experiments [10]. Accordingly, in this work, we take advantage of the TSEOS and use it to calculate H(T ) at selected pressures; after all, the TSEOS reproduces very well the behavior (and maxima) of ρ(T ) (see Fig. 1) and κ T (T ) (see Fig. 3). Specifically, for a given pressure, we use the polynomial expression of G(T ) given by the TSEOS and obtain an analytical expression for H(T ) using the Gibbs-Helmholtz equation, The

C. Diffusion Coefficient
We also calculate the self-diffusion coefficient of q-TIP4P/F water D(T ) as function of temperature at selected pressures. To obtain D(T ), we employ the same methodology used in Ref. [46]. Briefly, we first calculate the mean-square displacement (MSD) of the oxygen atoms/ring-polymer centroids. D(T ) is then evaluated from the slope of the MSD at long times, in the so-called diffusive regime. Fig. 6(a)   and hence, this temperature is underestimated in both MD and PIMD simulations.
Overall, the results from classical MD simulations in Fig. 7(b) are consistent with the TSEOS. Accordingly, one may conclude that the reported location of the LLCP based on the TSEOS is robust in the case of classical MD simulations. For example, the κ T -maxima from MD simulations and from the TSEOS (blue triangles and blue solid lines) overlap; similarly, the corresponding C P -maxima lines (red triangles and red solid lines) also overlap.
In addition, we find that the T M CT (P ) line shows a sudden change in slope at the intersection with the Widom line. This is consistent with the view that the Widom line separates the LDL-like liquid at low pressures from the HDL-like liquid at high pressures [72].

IV. SUMMARY AND DISCUSSION
We performed classical MD and PIMD simulations of H 2 O using the q-TIP4P/F model over a wide range of temperatures and pressures. At supercritical temperatures, most properties studied are practically insensitive to whether one employs classical MD and PIMD simulations (−100 ≤ P ≤ 500 MPa). Specifically, the ρ(T ) and κ T (T ) obtained from classical MD or PIMD simulations overlap (within error bars) down to T ≈ 225 K and T ≈ 200 K, respectively ( Fig. 1 and Fig. S2 of the SI). In the case of C P (T ), the quantitative values vary for classical MD and PIMD results, but this is expected to occur due to NQE [66]. Nonetheless, the qualitative behavior of C P (T ) is independent on whether NQE are included or excluded. Similarly, the behavior of D(T ) is not affected whether one employs classical MD or PIMD simulations down to T ≈ 260 K (Fig. 6a). Relative to the classical MD simulations, including NQE (PIMD simulations) increases the values of D(T ) at T < 260 K (inset of Fig. 6a). In both cases, the D(T ) from computer simulations are in good agreement with experiments where data is available (Fig. 6b and 6c). Consistent with the existence of a LLCP, our study shows the presence of loci of maxima in C P and κ T in the P-T phase diagram of q-TIP4P/F water. These anomalous maxima lines, together with the loci of maxima in D and ρ are included in Fig. 7. We stress that the location of the LLCPs reported in this work are estimations provided by the TSEOS.
While our results conclusively show that the LLCP in H 2 O shifts to lower T when NQE are included, obtaining the exact values of (ρ c , T c , P c ) may require additional data at T < 180 K (particularly for the case of PIMD simulations of H 2 O where T c is low).
Overall, our results for H 2 O (e.g., Fig. 7) are consistent with previous classical computer simulations of water using the (rigid) ST2 [42] and TIP4P/2005 [15] water models. It follows that the present study validates the LLPT hypothesis for water to the case where NQE are included. We note, however, that the LLCP in q-TIP4P/F water, as well as in ST2, TIP4P/2005, and TIP4P/Ice water, is located at pressures and temperatures that are off compared to the experimental predictions [6,46]. This explains why these water models are unable to reproduce the sharp increase in C P (T ) and κ T (T ) observed in experiments at P = 0.1 MPa [5,10]. Since C P (T ) and κ T (T ) quantify the fluctuations in entropy and volume, respectively, it follows that, at least for the q-TIP4P/F model studied, the inclusion of NQE (atoms delocalization) is not sufficient to reproduce the anomalously large fluctuations (in entropy and volume) of real water in the supercooled regime. Accordingly, additional sources of fluctuations may be missed in rigid (e.g., TIP4P/2005, ST2) as well as flexible water models, such as q-TIP4P/F model.
We also performed extensive PIMD simulations of heavy water using the q-TIP4P/F model. The results are summarized in the phase diagram of Fig. 8 (see also the SI). The shift considerably the location of the LLCP [47,48]. In particular, the differences in the relative values of (ρ c , P c , T c ) between q-TIP4P/F H 2 O and D 2 O are somewhat consistent with the predictions from experiments in glassy water (∆P c = 50 MPa, ∆T c = 10 K) [49,50] and with the relative location of the κ T -maximum at 1 bar [10].

V. DATA AVAILABILITY
All study data are included in the article and/or SI. G. analyzed data; and A. E., G. E. L., and N. G. wrote the paper.

VIII. ADDITIONAL INFORMATION
The authors declare no conflict of interest.