Evidence of a liquid–liquid phase transition in H\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2O and D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2O from path-integral molecular dynamics simulations

We perform path-integral molecular dynamics (PIMD), ring-polymer MD (RPMD), and classical MD simulations of H\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2O and D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2O using the q-TIP4P/F water model over a wide range of temperatures and pressures. The density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (T)$$\end{document}ρ(T), isothermal compressibility \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa _T(T)$$\end{document}κT(T), and self-diffusion coefficients D(T) of H\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2O and D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2O are in excellent agreement with available experimental data; the isobaric heat capacity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_P(T)$$\end{document}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 H\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2O and D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2O exhibit a liquid-liquid critical point (LLCP) at low temperatures and positive pressures. The data from PIMD/MD for H\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2O and D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2O 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 H\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2O, from PIMD simulations, is located at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_c = 167 \pm 9$$\end{document}Pc=167±9 MPa, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_c = 159 \pm 6$$\end{document}Tc=159±6 K, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _c = 1.02 \pm 0.01$$\end{document}ρc=1.02±0.01 g/cm\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^3$$\end{document}3. Isotope substitution effects are important; the LLCP location in q-TIP4P/F D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2O is estimated to be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_c = 176 \pm 4$$\end{document}Pc=176±4 MPa, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_c = 177 \pm 2$$\end{document}Tc=177±2 K, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _c = 1.13 \pm 0.01$$\end{document}ρc=1.13±0.01 g/cm\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^3$$\end{document}3. 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, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_c = 203 \pm 4$$\end{document}Pc=203±4 MPa, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_c = 175 \pm 2$$\end{document}Tc=175±2 K, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _c = 1.03 \pm 0.01$$\end{document}ρc=1.03±0.01 g/cm\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^3$$\end{document}3). 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_c$$\end{document}Tc for D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2O and, particularly, H\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2O suggest that improved water models are needed for the study of supercooled water.


Simulation method
Our results are based on PIMD/RPMD 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 54 . This model is based on the TIP4P/2005 model for water 55 , 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,54 . 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 Supplementary 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 56 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. In order to ensure that the conclusions in Ref. 46 applied to the pressures we considered in this work, we have also performed additional PIMD simulations using n b = 72 beads per ring-polymer (see SI). Consistent with Ref. 46 , we found that most of the thermodynamic and dynamical properties converged with n b = 32 , with the enthalpy being the only expected exception. 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) 57 . The OpenMM software package is also used to perform the RPMD simulations which are used for the calculation of the diffusion coefficients of H 2 O and D 2 O. Details on the calculation of the diffusion coefficients can be found in Ref. 46 . We note that in the OpenMM software package, the RPMD application sets the mass of the ring-polymer beads to the physical mass of the corresponding atom. When used Scientific Reports | (2022) 12:6004 | https://doi.org/10.1038/s41598-022-09525-x www.nature.com/scientificreports/ to calculate static equilibrium properties (energies, density, and RDF), the RPMD simulations reduce to PIMD simulations. In all PIMD/RPMD 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/RPMD 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/RPMD 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 .

Results
The results are presented as follows. In "Liquid-liquid phase transition", 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 "Thermodynamic response functions: κ T and C p ". The self-diffusion coefficient of H 2 O and D 2 O are the topic of "Diffusion coefficient" where we identify the anomalous locus of diffusivity maxima. A complete phase diagram for H 2 O is presented in "Phase diagram" where we also discuss similar results obtained for D 2 O. Figure 1a shows the density of H 2 O from both classical MD (open circles) and PIMD simulations (solid circles) along the isobars P = −100, 0.1, 100, . . . , 500 MPa and at temperatures in the range T = 180-375 K. The densities from MD and PIMD simulations overlap practically throughout the entire temperature and pressure range considered with some deviations being noticeable only at P = 100-200 MPa and T < 240 K. As we will show below, these T-P conditions are in the proximity of the LLCP. We note that the densities of q-TIP4P/F H 2 O are in remarkable good agreement with the corresponding experimental values. To show this, we include in Fig. 1b the densities from experiments and PIMD simulations of q-TIP4P/F H 2 O. Deviations between experiments and PIMD simulations are small, �ρ < 0.02-0.03 g/cm 3 , and are present only at P > 200 MPa and T < 250 K (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]58 . In Fig. 1c,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 Fig. 1c,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.

Liquid-liquid phase transition.
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 314 , 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 318 ; similar results were found in an 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 1 and are indicated in Fig. 1c,d by a red star. It follows that including NQE can shift the location of the LLCP. Specifically, relative to the classical case (MD simulations), adding NQE (PIMD simulations) lowers T c and P c by 16 ± 6 K, 36 ± 10 MPa, respectively; ρ c is not affected by the inclusion of NQE. 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 also compare the volumes predicted by the TSEOS with the corresponding values obtained from our MD and PIMD simulations. Fig. 2a,b show P(V) along isotherms based on the classical MD and PIMD simulations, respectively. In both cases, MD and PIMD simulations, the values of P(V) obtained from the TSEOS are in excellent agreement with our simulations. This strongly indicates that the TSEOS is reliable in predicting the properties of q-TIP4P/F water from both MD and PIMD simulations. We note that the P(V) isotherms shown in Fig. 2a,b seem to develop an inflection point as the temperature decreases, consistent with the existence of a LLCP at T < 200 K. Similarly, as shown in Fig. 2c, the potential energy PE(V) along isotherms is an increasing function of V at high temperatures but it develops a concave region (i.e., ∂ 2 PE/∂V 2 N,T < 0 ) at low Thermodynamic response functions: κ T and C p . We obtain the isothermal compressibility of q-TIP4P/F water by calculating the density fluctuations of the system 63 , where . . . indicates average over time and k B is the Boltzmann's constant. Fig. 3a,b show the κ T (T) for H 2 O obtained from PIMD simulations (solid circles) together with available experimental data (open symbols) at low and high pressures, respectively. At P ≥ 200 MPa (Fig. 3b), the experimental and PIMD simulation values of κ T (T) practically overlap; a similar agreement is found at P = 100 MPa (Fig. 3a). However, at P = 0.1 MPa, where more experimental data is available, the experimental κ T (T) increases more rapidly upon cooling than found in PIMD simulations. Hence, relative to real water, the density fluctuation in q-TIP4P/F water are underestimated at P = 0.1 MPa and in the supercooled regime. We note that the values of κ T (T) obtained from classical MD 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 Fig. 3a; the experimental κ T -maximum occurs at T = 228 K ( P = 0.1 MPa; open right triangles) and it is very sharp. While the κ T -maximum in PIMD simulations occurs at a similar temperature ( T = 230 K), this maximum is much smaller and wider relative to the experiments. Within the LLCP hypothesis scenario, the κ T -maximum is expected to increase as one approaches the LLCP and it should diverge at the LLCP. This is fully consistent with the PIMD simulations results shown in Fig. 3a. 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,64 .

is included in
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, A comparison of the values of κ T (T) obtained from the TSEOS and our MD/PIMD simulations are presented in Fig. 3c,d. The predictions from the TSEOS agree rather well with the MD simulation results [inset of Fig. 3d]. In the case of PIMD simulations [inset of Fig. 3c], the TSEOS provides compressibility values that are in good www.nature.com/scientificreports/ agreement with the simulation results at high temperatures. However, at lower temperatures, the TSEOS predicts slightly larger maxima in κ T that are shifted to lower temperatures relative to the simulations. This suggests that, the location of the LLCP in q-TIP4P/F water from PIMD may be located at slightly lower T c and/or higher P c relative to the corresponding estimated values resulting from the TSEOS. 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 TSEOS predictions are in excellent agreement with our simulations throughout the entire temperature and pressure range considered in this work. Figure 5a,b show, respectively, the C P (T) of q-TIP4P/F water from PIMD and classical MD simulations at selected pressures, above and below the estimated LLCP pressure; open symbols are experimental values. The C P (T) from classical MD and PIMD simulations are qualitatively similar. Specifically, at the temperature studied, C P (T) exhibits a maximum at approximately P ≤ 200 MPa. This C P -maximum increases and shifts to lower temperatures as the pressure increases towards the LLCP pressure. At P > 200 MPa > P c , C P (T) is a monotonic decreasing function of T. Note that classical MD simulations predict much larger values of C P (T) than found in PIMD simulations (which is known to occur when NQE are omitted 68 ).
Differences between the experimental data and MD/PIMD simulations are noticeable. For example, as shown in Fig. 5a, at P ≥ 100 MPa, PIMD simulations predict that C P (T) decreases upon heating while experiments show the opposite behavior. In particular, at P = 0.1 MPa, the C P (T) of q-TIP4P/F water is in semiquantitative   www.nature.com/scientificreports/ agreement with experiments down to T ≈ 240 K. The maximum in C P (T) , at P = 0.1 MPa, occurs at 228 K and 216 K in experiments and q-TIP4P/F water, respectively. However, the maxima of C P (T) in q-TIP4P/F water is much smaller and wider than found in the experiments of Pathak et al. 10 . This is consistent with the estimated location of the LLCP in experiments and in our simulations. The LLCP in real water is estimated to be located at P C ≈ 50-100 MPa and T C ≈ 220 K, while in q-TIP4P/F water we find T c =159 K and P c = 167 MPa. Accordingly, at P = 0.1 MPa, experiments are much closer to the corresponding LLCP than in the case of q-TIP4P/F water 6 . Our results for C P (T) also imply that PIMD and MD simulations of q-TIP4P/F water cannot reproduce the entropy fluctuations observed in real water at low temperatures and pressures.
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, using the RPMD simulation technique, 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. Figure Fig. 6b,c along different isotherms. At high temperatures, approximately T > 300 K (Fig. 6c), the values of D(P) from classical MD and PIMD simulations practically overlap with the available experimental data at all pressures studied. Instead, at T < 300 K (Fig. 6b), our simulation results predict values of D(P) that deviate by up to a factor of 4 from the corresponding experimental data, depending on pressure. Interestingly, overall, classical MD simulations are in slightly better agreement with experiments compared to the RPMD simulation results. Not surprising, as shown in the inset of Fig. 6a, the inclusion of NQE increases water diffusivities, particularly upon cooling within the supercooled regime.
One of the main points of Fig. 6b is the presence of an anomalous maximum in the diffusion coefficient of q-TIP4P/F water. Such a D-maximum is consistent with experiments and implies that there is a range of temperatures at which D increases (anomalously) with increasing P. Phase diagram. Figure 7a,b show, respectively, the phase diagram of q-TIP4P/F water obtained from PIMD and classical MD simulations. These phase diagrams include the LLCP, coexistence line, C P -maxima, κ T -maxima, and Widom line calculated from the TSEOS. Also included are the lines of ρ-maxima, D-maxima, C P -maxima, κ T -maxima, and κ T -minima obtained from our MD/PIMD simulations. We also show the maxima/ minima of these properties reported in experiments, where available. The liquid-vapor boundary lines shown in Fig. 7a,b correspond to the conditions at which spontaneous cavitation occurs during our computer simulation.
The phase diagrams of q-TIP4P/F water resulting from the classical MD and PIMD simulations are qualitatively similar. As found previously in ST2 42 and TIP4P/2005 water 15 , the C P and κ T -maxima lines obtained from the TSEOS originate at the LLCP and deviate from each other at higher temperatures. The κ T -maxima line (blue up-triangles) connects smoothly with the κ T -minima line (blue down-triangles). In addition, as shown in Ref. 35 , the point in Figs. 7 and 8 where the ρ-maxima line has infinite slope is located on the κ T -extrema line, at which (∂κ T /∂T) P = 0 . The ρ-maxima line has a nose shape. In particular, our simulations suggest that, at low pressures, the ρ-maxima line is re-entrant and deviates from the liquid-vapor boundary line [see, in particular, www.nature.com/scientificreports/ Fig. 7b]. This implies that the 're-entrant spinodal' hypothesis scenario proposed to explain water anomalous behavior is not supported by MD/PIMD simulations of q-TIP4P/F water. Hence, our results are consistent with Refs. 7,72 , where computer simulations performed using the ST2 and TIP4P water model found that the spinodal line is also not re-entrant. We note that Fig. 7a,b also include available experimental data. When compared with experiments, both classical MD and PIMD reproduce correctly the location of the κ T -maxima and minima lines; the PIMD simulation results reproduce slightly better the location of the ρ-maxima line.
Regarding the D-maxima line, both MD and RPMD simulations overestimate the corresponding pressures relative to the experiments, with the RPMD simulations performing slightly better. We also include in Fig. 7a,b the MCT temperatures extracted from Fig. 6a. The experimental MCT temperature at P = 0.1 MPa is T MCT = 221 K 73 and hence, this temperature is underestimated in both MD and RPMD simulations.
Overall, the results from classical MD simulations in Fig. 7b 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 MCT (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. The sharp crossover in T MCT (P) is reminiscent of the glass transition temperature of water as a function of pressure which shows a sudden change as the system evolves from LDL (low pressure) to HDL (high pressure) 31,74,75 .
In the case of PIMD simulations (Fig. 7a), the κ T -maxima line obtained from the TSEOS (blue solid line) is located at slightly lower pressure relative to the PIMD simulation results (blue up-triangles). Similarly, the Widom line predicted by the TSEOS is located at a pressure slightly lower than the pressure at which the slope of T MCT (P) suddenly changes. Hence, in the case of PIMD simulations for water, the reported location of the LLCP may shift slightly if additional data points at T < 180 K are considered in the TSEOS calculation. Additional PIMD simulations are also needed to improve the determination of the κ T -maxima line at low temperatures.
The similarities in the phase diagrams of Fig. 7a,b imply that, at least for the q-TIP4P/water model, the LLPT hypothesis scenario remains a solid, viable explanation of water anomalous behavior even if NQE (i.e., atoms delocalization) are included. This is important since (i) the LLCP scenario has been tested only in classical MD simulations and, mostly, rigid water models, and (ii) the location of the LLCP is extremely sensitive to small variations in the water model considered (e.g., small changes in the partial charges of the H and O atoms can easily shift the location of the LLCP to (P, T) conditions that are physically inaccessible to the liquid state; see, e.g., Refs. 76,77 ). Overall, including NQE shifts the location of the LLCP, LLPT line, and maxima/minima lines towards lower temperatures (see also Ref. 47,48 ).

Summary and discussion
We performed classical MD, PIMD, and RPMD 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 68 . 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 RPMD simulations down to T ≈ 260 K (Fig. 6a). Relative to the classical MD simulations, including NQE (RPMD 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,c). Deviations between MD and PIMD simulations become noticeable at approximately P = 100-200 MPa and T < 225 K. Our results strongly indicate that at these thermodynamic conditions, q-TIP4P/F water exhibits a LLCP. Using a two-state equation of state, we estimate that the LLCP is located at (ρ c = 1.03 g/cm 3 , P c = 203 MPa, T c = 175 K) when NQE are not included (classical MD); including NQE (PIMD simulations) shifts the location of the LLCP to (ρ c = 1.02 g/cm 3 , P c = 167 MPa, T c = 159 K); see Fig. 1c,d. 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. Our estimation of the LLCP for H 2 O is based on PIMD simulations using n b = 32 beads per ring-polymer. While PIMD simulations with n b > 32 are computationally expensive, additional PIMD simulations employing n b > 32 beads per ringpolymer are desirable at low temperatures in order to obtain a more precise estimation of the LLCP location in H 2 O (q-TIP4P/F water). 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 provides a thermodynamic explanation of 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 . Specifically, these water models predict that P c > 150 MPa, while P c ≈ 50 − 100 MPa estimated from experiments 1,50 . Accordingly, computer simulations show a mild increase in κ T and C P , relative to experiments, upon isobaric cooling at P = 0.1 MPa. We note that the C P (T) and κ T (T) quantify the fluctuations in entropy and volume, respectively. Hence, from a microscopic point of view, the weak increase of κ T and C P upon isobaric cooling at P = 0.1 MPa is due to the inability of current water models (ST2, TIP4P/2005, TIP4P/Ice, etc) to reproduce the anomalously large fluctuations (in entropy and volume) of real water in the supercooled regime. At least for the q-TIP4P/F model studied, the inclusion of NQE (quantum fluctuations due to atoms delocalization) is not sufficient to reproduce the anomalously large fluctuations (in entropy and volume) of real water at low temperatures. 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 PIMD simulations confirm that isotope substitution has minor effects on the properties of water. While the phase diagram of D 2 O is qualitatively identical to the phase diagram of H 2 O, the location of the corresponding LLCP differ. Specifically, calculations based on the TSEOS, applied to PIMD data at T ≥ 190 K, predict that in the case of D 2 O, ( ρ c = 1.13 g/ cm 3 , P c = 176 MPa, T c = 177 K) which represents a non-negligible shift relative to H 2 O ( �ρ c = 0.11 g/cm 3 , P c ≈ 9 MPa, T c ≈ 18 K). This is important since computer simulations of water-like models show that the introduction of NQE can indeed 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 . The present study shows that many questions previously addressed in computational studies of supercooled water at low temperatures, using classical water models, are accessible via PIMD simulations. For example, it would be interesting to explore the relationship between dynamics and structure of water at low temperatures and how the Stokes-Einstein and Stokes-Einstein-Debye relationships are affected by isotope substitution [82][83][84][85] .

Data availibility
All study data are included in the article and/or SI.