Connection between water’s dynamical and structural properties: insights from ab initio simulations

Abstract

Among all fluids, water has always been of special concern for scientists from a broad variety of research fields due to its rich behavior. In particular, some questions remain unanswered nowadays concerning the temperature dependence of bulk and interfacial transport properties of supercooled and liquid water, e.g. regarding the fundamentals of the violation of the Stokes-Einstein relation in the supercooled regime or the subtle relation between structure and dynamical properties. Here we investigated the temperature dependence of the bulk transport properties from ab initio molecular dynamics based on density functional theory, down to the supercooled regime. We determined from a selection of functionals, that SCAN better describes the experimental viscosity and self-diffusion coefficient, although we found disagreements at the lowest temperatures. For a limited set of temperatures, we also explored the role of nuclear quantum effects on water dynamics using ab initio molecular dynamics that has been accelerated via a recently introduced machine learning approach. We then investigated the molecular mechanisms underlying the different functionals performance and assessed the validity of the Stokes-Einstein relation. We also explored the connection between structural properties and the transport coefficients, verifying the validity of the excess entropy scaling relations for all the functionals. These results pave the way to predict the transport coefficients from the radial distribution function, helping to develop better functionals. On this line, they indicate the importance of describing the long-range features of the radial distribution function.

Water is an ubiquitous liquid, essential for life on earth, and therefore constitutes one of the most important chemical substances. Despite the apparent simplicity of its chemical formula, water is a complex liquid, that after much effort, still evades our complete understanding at the molecular level [1]. Due to its critical relevance with regard to energy harvesting and water purification, several efforts have been carried in order to obtain molecular insights about water behavior under different thermodynamic conditions. Water molecular interactions arise from a balance between van der Waals and hydrogen bonding forces [2, 3], thus a complete description exclusively from classical force field (FF) simulations may hinder some critical mechanisms. Ab initio molecular dynamics (AIMD), where interatomic forces are computed from the electronic structure, may play a key role in understanding some important physical processes for bulk and confined water [4, 5, 6, 7] as, for instance, the controversial liquid-liquid critical point [8].

A very efficient approach to determine the electronic structure is density functional theory (DFT), based on a formulation of the many-body problem in terms of a functional of the electron density. So far however, many widely used approximations for the exchange-correlation functional do not provide a sufficiently accurate description of many of water properties [9, 10]. The main challenge of semi-local and hybrid density functional approximations in predicting the structure and energetics of water lies in their description of dispersion and exchange-overlap interactions [10, 11]. Additionally, nuclear quantum effects (NQEs) play an important role in determining water structure [12, 13, 14, 15]: while on the one hand, NQEs tend to strengthen the hydrogen bond, on the other hand, competing effects due to the spread of the protons in the normal direction tend to weaken it. Although NQEs can be modelled via ab initio path integral molecular dynamics (PIMD), accounting for this subtle competition adds an additional level of complexity [12]. Despite recent advances in describing the water structure and thermodynamic properties [16, 17], predicting transport properties from first principles represents a further challenge [18, 19, 20]. Limited work has been dedicated to the study of the temperature evolution of the diffusivity and of the shear viscosity with ab initio methods [8, 21, 3]. A clearer picture of the molecular mechanisms controlling the water viscosity and diffusion is needed, especially in the supercooled regime [22, 23], where water viscoelasticity is poorly understood [24, 25, 26] and the validity of the Stokes-Einstein (SE) relation at low temperatures remains an open question [27, 28, 29, 30, 31, 32].

Water dynamics is also crucial in the field of nanofluidics [33], where in particular the performances can be boosted by liquid-solid slip, arising from a competition between bulk liquid viscosity and interfacial friction [34, 35, 36, 37]. Further, reaching clearer insights into the molecular properties controlling water dynamics would enable to determine a relationship between the structural correlations and molecular transport. Establishing such a thermodynamic link between structure and dynamics would also be instrumental to improve the description of water via DFT.

Such connection has already been explored in the literature via e.g. free-volume models [38, 39], relationships between $g(r)$ and glass transition temperature [40], and the proposition of different structural descriptors [41, 42, 43], among which the entropy excess scaling, already employed for AIMD simulations of liquid metals [44] or water [45], stands out [46, 47, 48]. The excess of entropy, which can be decomposed in terms on the $N-$body radial distribution functions [49], has been proven to exhibit an exponential relation with the diffusion coefficient for multiple systems [50, 51]. In particular, for glass forming liquids such as supercooled binary mixtures and water, the approximation of the entropy excess by its two body contribution (related to an integral of a function of $g(r)$) has been shown to work well for a broad range of temperatures [52, 53, 54, 48, 55]. One of the main limitations for AIMD is its great need of resources as compared to their classical counterparts. Nevertheless, if the link between dynamics and structure is established, we would be able to predict the transport coefficients from structural properties, which require shorter simulation times to converge [56]. Aside of this, entropy excess scaling has also been used as a tool to bring insights into the molecular mechanisms underlying the SE relation [55].

In this report, we determine from a selection of density functionals commonly used to characterise water [11, 10, 57], which one better describes the temperature dependence of the water viscosity and self-diffusion coefficient in its liquid and supercooled state, in comparison with FF simulations using the TIP4P/2005 water model [58]. Additionally, we explore the connection between structural properties and the transport coefficients for all the functionals proposed via the two-body entropy parameter, presenting this physical descriptor as a path to develop better functionals and better compare with experimental results. We used AIMD simulations, describing the electronic structure within DFT, in the NVT ensemble to determine hydrodynamic bulk transport coefficients of three different density functionals: PBE [59] functional with Grimme’s D3 corrections [60, 61] (namely PBE-D3), optB88-vdW [62, 63] and SCAN [64]. While describing the electronic structure with the SCAN functional, we also included the role of NQEs by performing PIMD simulations for a limited set of temperatures, by employing a recently introduced machine learning approach to speed-up the calculation of the electronic structure problem [65]. Further simulation details can be found in Materials and Methods.

(a)

(b)

Figure 1: Temperature evolution for different functionals of (a) shear viscosity and (b) diffusion coefficient. A good agreement is found between the hydrodynamic radius measures (dotted lines) and the Green-Kubo ones (dashed lines), implying all the functionals verify Stokes-Einstein relation with the same hydrodynamic radius $R_{\mathrm{h}}=1\,$Å (see text for detail).

Results and Discussion

We display in Fig. 1(a) (dashed lines) the temperature evolution of the shear viscosity, determined from the long-time plateau of the Green-Kubo integral, $\eta_{\mathrm{GK}}$ (Eq. (8) in Materials and methods), for the different functionals. No plateau was observed for PBE-D3 at $T<360~$K and optB88-vdW at $T=260~$K. To benchmark the results, the same procedure was carried for FF simulations with the TIP4P/2005 water model [58], which provides an excellent description of experimental results for both viscosity and diffusion coefficient [66, 67, 37]. In Fig. 1(a) one can see that the viscosity obtained from the SCAN functional is in better agreement with FF at $330~$K and $360~$K, although between $260~$K and $300~$K it overestimates the viscosity by more than a factor of $1.7$. With regard to PBE-D3 and optB88-vdW, one observe from Fig. 1(a) that both functionals overestimate $\eta_{\mathrm{GK}}$ value. Overall, all functionals fail at describing the temperature evolution of the shear viscosity.

The diffusion coefficient $D_{\mathrm{PBC}}$ was determined from the slope of the mean squared displacement in the diffusive regime (see the supporting information, SI). In practice, because of hydrodynamic interactions between the periodic image boxes, a finite size correction for the diffusion coefficient has to be introduced [68, 66, 69]. For a cubic simulation box of size $L_{\mathrm{box}}$ with periodic boundary conditions:

$$D_{\mathrm{GK}}=D_{\mathrm{PBC}}+2.837\frac{k_{\mathrm{B}}T}{6\pi\eta_{\mathrm{GK}}L_{\mathrm{box}}},$$ (1)

with $k_{\mathrm{B}}$ the Boltzmann constant and $T$ the temperature. We denoted $D_{\mathrm{GK}}$ the diffusion coefficient obtained through Eq. (1) because we used $\eta_{\mathrm{GK}}$ in it. $D_{\mathrm{PBC}}$ could not be determined within our simulation times for PBE-D3 at $T<360~$K and optB88-vdW at $T=260~$K because the system did not enter in the diffusive regime. This result is consistent with the absence of a plateau for $\eta_{\mathrm{GK}}$. The corrected diffusion coefficient $D_{\mathrm{GK}}$ results are displayed in Fig. 1(b) (dashed lines). In analogy to $\eta_{\mathrm{GK}}$, one observes in this figure that SCAN is the functional that better describes water diffusion coefficient at high temperatures, although it fails at low $T$.

Generally, viscosity $\eta$ and diffusion $D$ are related through the SE relation:

$$D=\frac{k_{\mathrm{B}}T}{6\pi\eta R_{\mathrm{h}}},$$ (2)

with $R_{\mathrm{h}}$ the effective hydrodynamic radius of the molecules [70]. Even though the failure of this relation is well known at low temperatures [27, 28, 29, 30, 31, 32], it still remains valid for a broad range of temperature. We verified this statement by computing $R_{\mathrm{h}}$ for FF simulations, obtaining constant $R_{\mathrm{h}}\sim 1~$Å for the range of temperatures considered in the present study (see the SI).

Taking into account $D$ size correction Eq. (1) and SE relation Eq. (2), one can relate the viscosity to $D_{\mathrm{PBC}}$ and to the hydrodynamic radius:

$$\eta_{R_{\mathrm{h}}}=\frac{k_{\mathrm{B}}T}{6\pi D_{\mathrm{PBC}}}\quantity( \frac{1}{\Rh} - \frac{2.837}{\Lbox} ).$$ (3)

In the same way, one can also determine a relation for $D_{R_{\mathrm{h}}}$ independent of $\eta$ from Eq. (1) and Eq. (2):

$$D_{R_{\mathrm{h}}}=\frac{D_{\mathrm{PBC}}}{1-\frac{2.837R_{\mathrm{h}}}{L_{\mathrm{box}}}}.$$ (4)

Therefore, viscosity and diffusion can be determined exclusively from the slope of the mean squared displacement at long times by imposing the hydrodynamic radius $R_{\mathrm{h}}$. In order to test the applicability of this prediction, in Fig. 1 we display the results for $\eta_{R_{\mathrm{h}}}$ from Eq. (3) and $D_{R_{\mathrm{h}}}$ from Eq. (4) by imposing $R_{\mathrm{h}}=1\,$Å (value in agreement with the FF measures, see the SI). In Fig. 1 one can see a good match between the Green-Kubo and the hydrodynamic radius measures for both transport coefficients and for all the functionals considered, meaning that, although all the functionals fail in predicting viscosity and diffusion temperature dependence, all of them verify the SE relation with the same constant value of $R_{\mathrm{h}}$.

Figure 2: Dimensionless two-body entropy $s_{2}/k_{\mathrm{B}}$ for different functionals and for FF as a function of the temperature.

(a)

(b)

(c)

(d)

Figure 3: Radial distribution functions (a) and (b) and $-s_{2}/k_{\mathrm{B}}$ running integrals (c) and (c) of water at 300 K and 260 K obtained from the SCAN functional with classical and quantum nuclei (SCAN+NQEs) as well as from the FF simulations.

Having determined the transport properties for the different functionals, we proceed to explore their connection with the structure of water, given by the radial distribution function, $g(r)$. Specifically, we computed the structural descriptor $s_{2}$ (two-body excess entropy, see the SI), given by the integral [49]:

$$\frac{s_{2}}{k_{\mathrm{B}}}=-2\pi n\int_{0}^{\infty}r^{2}\quantity( g(r)\ln g(r) - g(r) +1 )\,\differential r,$$ (5)

with $n$ the number density of the system.

Figure 2 presents the temperature dependence of the dimensionless two-body entropy $s_{2}/k_{\mathrm{B}}$ for the different functionals, compared with FF results. We note that, although the $g(r)$ simulated with the TIP4P/2005 FF exhibits some discrepancies with respect to the experimental $g(r)$ [58], the change of $s_{2}$ with temperature is in qualitative agreement with experiments for a wide range of temperatures, down to the supercooled regime [71].

One can observe that SCAN is the functional that better describes the $s_{2}$ temperature evolution, as compared to FF. Interestingly, accounting for NQEs (see the SCAN+NQE in Fig. 2) does not produce significant changes in the behavior of $s_{2}/k_{\mathrm{B}}$, even at the low temperature of $260$ K. Whereas at the highest temperature of $360\,$K the optB88-vdW functional reproduces the water structure – as represented by the $s_{2}$ order parameter, using this functional in the supercooled regime gives rise to serious over structuring, thus leading to an overestimation of more than two times in the value of $s_{2}$ at $260\,$K. Instead, PBE-D3 fails at recovering the liquid two-body entropy at any temperature, and the integral in Eq. 5 reaches a plateau for the lowest temperatures, i.e. the $g(r)$ oscillations amplitude is not significantly affected by temperature, hinting at a possible glass transition [40].

To elucidate the connection between the water structure and the transport coefficients through the $s_{2}$ order parameter, we compare the radial distribution functions and the $s_{2}/k_{\mathrm{B}}$ running integrals for the FF and the SCAN functional at 260 K and at 300 K, both with classical and with quantum nuclei (see Fig. 3), while the comparison with the PBE-D3 and optB88-vdW functionals can be found in the SI. At 300 K the quantum nuclear $g(r)$ displays a lower first peak compared to the classical one and to the $g(r)$ from the FF, while the rise in the first peak also occurs at a slightly shorter distance. Thus, NQEs give rise to a less structured first coordination shell, in qualitative agreement with experimental results in H₂O and D₂O which display these same characteristics, respectively [72]. Beyond the first peak, an increased structure is observed in the $g(r)$ predicted with SCAN+NQE, compared to the case with classical nuclei and with the FF result, in disagreement with experiments on light and heavy water [72]. Such discrepancy has been observed in previous PIMD simulations obtained with semi-local density functionals [12, 73], where it was pointed out that the origin of the increased structure of the second peak arises from the destabilization of interstitial hydrogen bonded configurations occurring in quantum nuclear simulations, and more accurate descriptions might alter this balance and lead to a less structured second and third solvation shells [16, 74]. Shifting the focus on the $s_{2}/k_{\mathrm{B}}$ running integral (see Fig. 3(c)), it can be noticed that although the largest contribution to the limiting value of $s_{2}/k_{\mathrm{B}}$ arises from the first solvation shell, a non-negligible part is also due to the oscillations beyond the first solvation shell. Thus, as a result of a less structured first solvation shell and of more structured second and third shells, the SCAN+NQE functional predicts a limiting value of $s_{2}/k_{\mathrm{B}}$ that is similar to SCAN with classical nuclei, whereas the $s_{2}/k_{\mathrm{B}}$ value predicted with FF is visibly lower. Interestingly, the differences between the quantum and classical nuclear $g(r)$, observed at $300$ K, are attenuated at $260\,$K, also when comparing with the FF results. This leads to a value of $s_{2}/k_{\mathrm{B}}$ that is remarkably similar for all the three methods at $260\,$K, as seen in Fig. 3(d).

Finally, in order to establish a relation between structure and transport coefficients, we proceed to test the entropy excess scaling laws. With that regard, it has been verified that the entropy excess $s_{\mathrm{ex}}$ can be approximated by the two-body contribution $s_{\mathrm{ex}}\simeq s_{2}$ for water and supercooled binary mixtures [52, 53, 54, 55]. As described in Ref. 75, $s_{2}$ is constructed from the oxygen-oxygen, oxygen-hydrogen, and hydrogen-hydrogen pair distributions. Still, the scaling laws hold well by just computing $s_{2}$ from the oxygen-oxygen $g(r)$, amounting to only considering the translational 2-body entropy, with a difference of a factor 3 between both estimates, so $s_{\mathrm{ex}}\simeq 3s_{2}$ from our results. It can be shown [76] (see the SI) that the dimensionless diffusion coefficient $D/D_{0}$ is expected to scale as:

$$\frac{D}{D_{0}}=A\exp(-B\, s_2/\kB),$$ (6)

with $D_{0}=l_{0}\sqrt{k_{\mathrm{B}}T/m}$ (where $l_{0}=n^{-1/3}$ is the average interparticle distance and $m$ is the fluid mass), and $A$ and $B$ dimensionless constants at a given fluid density. Considering this Eq. (6), together with the SE equation, and assuming $R_{\mathrm{h}}\sim l_{0}$, one can expect a scaling for the dimensionless viscosity $\eta/\eta_{0}$ as:

$$\frac{\eta}{\eta_{0}}=A^{\prime}\exp(-B' \, s_2/\kB),$$ (7)

with $\eta_{0}=\sqrt{mk_{\mathrm{B}}T}/l_{0}^{2}$. From the SE relation one expects $B^{\prime}=-B$.

(a)

(b)

Figure 4: Reduced (a) viscosity $\eta/\eta_{0}$ and (b) diffusion coefficient $D/D_{0}$, defined in Eq. (6) and Eq. (7), as a function of the dimensionless two-body entropy $s_{2}/k_{\mathrm{B}}$ for different functionals and FF simulations. In continuum line are represented the respective exponential fits for each functional. The fit results are detailed in Table 1. The color and marker style representing the different functionals is the same as in Fig. 2.

We tested Eq. (6) and Eq. (7) for the different functionals. Figure 4 shows the results for the dimensionless transport coefficients as a function of the two-body entropy excess for the different functionals. One can see that, although the functionals predict different transport coefficients (Fig. 1) and $s_{2}$ results (Fig. 5), all of them verify an exponential scaling of $\eta_{\mathrm{GK}}/\eta_{0}$ and $D_{\mathrm{GK}}/D_{0}$ with $s_{2}$. Therefore, we fitted the relations in Eq. (6) and Eq. (7) for optB88-vdW, SCAN, and FF (continuous lines in Fig. 4). No fit was performed for PBE-D3 due to the single value measure we could report for this functional. The fit results are indicated in Table 1. One can observe that, although out of the error bars, the fit parameters for SCAN and FF are the closest ones and that, for all functionals, $B^{\prime}=-B$, implying a verification of the SE relation, Eq. (2).

Table 1: Fit parameters of the two-body excess entropy scaling relation for the dimensionless viscosity and diffusion coefficient, for different functionals and FF simulations. The fit corresponds to the function $y=A\exp(-B s_2/\kB)$ with $y$ the dimensionless viscosity $\eta_{\mathrm{GK}}/\eta_{0}$, following Eq. (7), and diffusion coefficient $D_{\mathrm{GK}}/D_{0}$, following Eq. (6).

	$\eta_{\mathrm{GK}}/\eta_{0}$		$D_{\mathrm{GK}}/D_{0}$
	$A^{\prime}$	$B^{\prime}$	$A$	$B$
	$(\times 10^{-1})$		$(\times 10^{-1})$
optB88-vdW	$4.29(1.39)$	$4.11(0.34)$	$3.52(0.29)$	$-3.97(0.09)$
SCAN	$1.79(0.48)$	$5.31(0.31)$	$8.17(2.49)$	$-5.24(0.36)$
FF	$1.92(0.26)$	$4.52(0.18)$	$7.73(1.25)$	$-4.58(0.21)$

One can exploit the exponential relationship between the transport coefficients and $s_{2}$ to predict transport properties from structural ones: once the fitting parameters in Eqs. (6-7) have been extracted by calculating the dependence of $\eta$ and $D$ on $s_{2}$ for a limited set of temperatures, the value of the transport coefficients can be obtained just from the calculation of the $s_{2}$ order parameter via the radial distribution function also for a wider temperature range. Indeed, generally structural properties such as the $g(r)$ require shorter simulations to converge, especially when using force based methods, as the one proposed in Ref. 56, to reduce the variance when compared to the conventional strategies based on particles positions binning. Figure 5 presents the Green-Kubo results and their comparison with the prediction resulting from the fit via $s_{2}$. One can see good agreement between the explicit calculation of transport coefficients and their and predictions via $s_{2}$ for all the data. Also, although the transport coefficients could not be calculated explicitly for the optB88-vdWfunctional at the lowest temperature of $260$ K, they could be determined from the exponential relationship with $s_{2}$, yielding an exceedingly high viscosity and low diffusion, and thus verifying the failure of this functional in reproducing the temperature dependence of both transport coefficients.

(a)

(b)

Figure 5: Temperature evolution for different functionals of (a) shear viscosity from Eq. (7) and (b) diffusion coefficient from Eq. (6) with fit parameters from Table 1. A good agreement is found between the $s_{2}$ based prediction (dotted lines) and the explicit calculation (dashed lines), verifying the link between the structure and the transport coefficients. The color and marker style representing the different functionals is the same as in Fig. 1.

Conclusions

In summary, among the selected DFT functionals, SCAN best captures the temperature evolution of water transport properties – as described by the the accurate TIP4P/2005 FF – despite the disagreement observed at temperatures of 300 K and below. We detected large discrepancies between functionals, with a major failure of PBE-D3, which is far too viscous. Despite these discrepancies, the SE relation was observed to hold in the considered temperature range for all the functionals and, moreover, all of them predicted the same hydrodynamic radius $R_{\mathrm{h}}\sim 1~$Å. This property, together with the finite size correction for the diffusion coefficient, allowed us to propose a measure of viscosity $\eta_{R_{\mathrm{h}}}$ and diffusion coefficient $D_{R_{\mathrm{h}}}$, based only on the slope of the mean squared displacement in the diffusive regime $D_{\mathrm{PBC}}$, for known box size and fixed $R_{\mathrm{h}}$.

Motivated by a possible connection between dynamics and structure, we computed the radial distribution functions for the different functionals. Analogously to the transport coefficient results, we observed that SCAN radial distribution function is the one that better compares to FF, with little differences at the lowest temperatures in the second and third solvation shells, whereas PBE-D3 was far more structured that SCAN and optB88-vdW at high temperatures, in agreement with the high viscosity value measured for this functional. An explicit relationship between dynamic and structure can be established through the two-body entropy excess, which is an integral of a function of the radial distribution function. We verified the exponential relation on $s_{2}$ for both reduced bulk transport coefficients and, although the connection between $s_{2}$ and transport properties is not univocal, as different fitting parameters were used for each functional, the fitting parameters all have the same order of magnitude. We suggest that the non-universality of the exponential relation is due to the use of the translational two-body excess entropy; a more universal relation could be observed by using the full two-body excess entropy, although the related structural features would be less easy to interpret.

Finally, based on the established exponential relation between the bulk transport coefficients, we computed both viscosity and diffusion coefficient from the $s_{2}$ results and the fitting parameters. This allowed us to estimate transport coefficients for functionals strongly structured (for instance optB88-vdW at $260\,$K), which present such a high viscosity value that longer simulations are needed in order to observe a well-defined plateau in the Green-Kubo integral. Therefore we propose here that, once the exponential dependence has been determined for a few temperatures, the viscosity and diffusion coefficient can be determined only from structural properties, which typically require shorter simulation times to converge [56]. This can be a useful technique to apply in order to calculate transport coefficients for very viscous systems, where the associated time-scales are far from the ones computationally reachable with AIMD simulations. The connection between transport coefficients and the radial distribution function via the two-body entropy excess also establishes some guidelines to choose a functional for simulations of nanofluidic systems, where a functional which better reproduces water’s structure will more likely reproduce its dynamical properties. The $s_{2}$ order parameter can be also employed as an instrumental tool to gauge the potential of DFT or of high accuracy methods in describing dynamics without computing transport properties explicitly, where the comparison between different $s_{2}$ values becomes more straightforward than the comparison between two $g(r)$ profiles, or just the value of the $g(r)$ minimum or maximum, which does not ensure a full structure correlation. Indeed, from the $s_{2}(r)$ running integrals, we showed the importance of reproducing not only the first solvation shell of the $g(r)$ but also the long range structure, which is a non-negligible contribution to $s_{2}$ value. This feature, together with the scaling behavior of the bulk transport coefficients as a function of entropy, suggests that it is important that DFT reproduces not only the first peak in the $g(r)$, but also its long range behavior, which is critical to obtain an accurate description of dynamical properties such as viscosity and diffusion coefficient.

It is worth discussing the possible origins of the discrepancies of the temperature evolution of of the viscosity and of the diffusion coefficient, especially striking for the SCAN functional, which shows good agreement with FF and experiments at high temperatures but it overestimates the viscosity and underestimates the diffusivity at low temperatures. Although one might expect that the inclusion of zero-point energy and quantum tunneling, which become increasingly relevant at lower temperatures, would play an important role, we have shown that taking NQEs into account did not improve upon the results obtained with classical nuclei. As such, the most likely source of discrepancy lies in the approximate description of the electronic structure with the chosen functionals. Capturing the delicate balance between van der Waals dispersion and exchange interactions constitutes the main challenge for the description of water [10] and it is also critical in order to predict water transport properties below room temperature. It remains to be seen whether the use of high-accuracy methods such as the random-phase approximation (RPA), Møller–Plesset perturbation theory [11] and quantum Monte Carlo [77] would improve upon the current description of water transport properties also in super-cooled conditions. Recent results on the diffusion coefficient for a wide range of temperatures obtained with RPA, also including the role of NQEs are promising in this regard [78]. As a further interesting perspective of this work, the established connection between the structure and dynamics might reveal what tips the balance between the strengthening and the weakening quantum delocalization effects of the H-bond network in water. Since diffusion is found experimentally to vary significantly between D₂O and H₂O [12], and $s_{2}$ is directly connected to diffusion, one can gauge the impact of the competing H-bond strengthening and weakening NQEs on the structure and the dynamics directly from the experimental measurements of the diffusivity and of the radial distribution function of light and heavy water at different temperatures.

Materials and methods

Simulation Details

We performed AIMD simulations of 32 water molecules in bulk using DFT with the CP2K code (development version) [79], which employs the Gaussian and Plane waves (GPW) method to describe the wave-function and the electron density and to solve the Kohn-Sham equations [80]. Three different density functionals were considered: PBE [59] functional with Grimme’s D3 corrections [60, 61], optB88-vdW [62, 63] and SCAN [64]. The electronic structure problem was solved within the Born-Oppenheimer approximation for 5 different temperatures ($T=\{260,270,300,330,360\}\,$K (the two lowest ones corresponding to the expected supercooled regime) controlled via the Nosé-Hoover thermostat. We worked at constant volume with a box size such that $\rho=1\,$g/cm³ ($L_{\mathrm{box}}=9.85\,$Å for 32 water molecules). The running time was $\simeq 120\,$ps for all functionals and temperatures except optB88-vdW and $T=\{260,270\}\,$K, with running time $\simeq 240\,$ps. The timestep considered was $0.5\,$ps. The initial configuration for all the functionals corresponded to the steady state at the given temperature obtained from force field (FF) MD after a running time of $200\,$ps. The energy cutoff for plane waves was $600\,$Ry for PBE-D3 and optB88-vdW, and $800\,$Ry for SCAN, and the localized Gaussian basis set was short range molecularly optimized double-$\zeta$ valence polarized (DZVP-MOLOPT-SR) [81].

NQEs were modelled using PIMD simulations with a thermostatted ring polymer contraction scheme using 24 replicas [82, 83]. PIMD simulations were performed using the i-PI code [84] together with CP2K, where the former is used to propagate the quantum nuclear dynamics, whereas the latter is used for the optimization of the wave-function and to calculate the forces on each atom. Sampling in the canonical ensemble has been carried out at 260 K and at 300 K for 35 ps. Further, in the case of the PIMD simulations, the electronic structure problem is solved using subsystem density functional theory within the Kim-Gordon (KG) scheme, where the shortcomings of the electronic kinetic energy term of KG-DFT are addressed by correcting this term via a $\Delta$-machine learning approach [85]. Specifically, we use a neural-network potential based on the Behler-Parrinello scheme [86] to learn the difference in the total energy and atomic forces between KS-DFT and KG-DFT (see Ref. 65 and the SI for further details). The resulting $\Delta$-learning method provides the same accuracy of KS-DFT at the lower computational cost of KG-DFT.

We also performed force field (classical MD) simulations via the LAMMPS package [87]. Analogously to AIMD, we worked in the NVT ensemble with the temperature controlled via a Nosé-Hoover thermostat and with a volume such that $\rho=1\,$g/cm³. Three different box sizes were considered: $32$ water molecules ($L_{\mathrm{box}}=9.85\,$Å), $64$ water molecules ($L_{\mathrm{box}}=12.42\,$Å) and $128$ water molecules ($L_{\mathrm{box}}=15.64\,$Å). The water model considered in all the cases was TIP4P/2005 [58].

Shear Viscosity and Diffusion Coefficient

For all the simulations, we determined the shear viscosity from the Green-Kubo relation:

$$\eta_{\mathrm{GK}}=\frac{V}{k_{\mathrm{B}}T}\int_{0}^{\infty}\expectationvalue{p_{ij}(t)p_{ij}(0)}\differential t,$$ (8)

with $V$ the volume, $k_{\mathrm{B}}$ the Boltzmann constant, $T$ the temperature and $p_{ij}=\{p_{xy},p_{xz},p_{yz}\}$ the non-diagonal components of the stress tensor.

The error bars were computed within $60\%$ of confidence level. For viscosity, the total MD stress was divided in three time slots of equal length, each of them containing three independent measures of $\eta$. $\eta_{\mathrm{GK}}$ was measured at the time where the $\eta(t)$ running integral reached a time plateau (see the SI). Therefore 9 independent viscosity values were computed for each functional at a given temperature. For the diffusion coefficient, the first $20\,$ps were removed from the trajectory so the system equilibration from the initial configuration does not affect the mean squared displacement results. From them, 3 independent measures of $D_{\mathrm{PBC}}$ were obtained from the three independent Cartesian components.