Water is arguably the most important substance on this planet. The remarkable functionalities of liquid water arise from a series of its unique and mysterious properties. The density anomaly enables water to freeze from the top down, which is essential to sustain the life and natural evolution on earth; 1,2 the frequent autoionization processes and the subsequent fast proton transport determine its amphoteric nature that underpins all the acid-base chemical reactions; 3,4 the hydrophobic effect in an aqueous environment is believed to be the driving force behind protein folding. 5,6 To unravel all these biological and chemical processes, a genuinely predictive model of water holds the key. Not surprisingly, understanding the hydrogen (H)-bond structure of water continues to command intense attention from both experimentalists and theoreticians. 4,[7][8][9][10][11][12][13][14][15][16][17][18] Advanced experimental techniques have been extensively applied to the study of water.
The arrangement of water molecules can be probed by neutron scattering or X-ray diffraction experiments; 9,10 the H-bond dynamics are detectable by infrared (IR), 19,20 sum-frequency vibrational, [21][22][23] and pump-probe 24 spectroscopy measurements; the photoemission 8 and X-ray absorption 7 spectra provide the electronic structure information of water. Nevertheless, the available experimental techniques only provide ensemble averaged information on water. For insight into the microscopic details, scientists in the field mainly rely on the theoretical modeling of molecules to further understand the structure, dynamics, and electronic properties of water. In order to rationalize the experimental data, an accurate first-principles theory is crucial because its predictions are based on quantum mechanical principles without any empirical input. In this regard, ab initio molecular dynamics (AIMD) 25 provides an ideal framework to model liquid water at the microscopic level. In AIMD, the driving force at each timestep is determined by density functional theory (DFT). 26,27 In DFT, the interacting many-electron problem in any condensed system can be exactly mapped onto a single-body problem with the inclusion of an exchange-correlation (XC) term. However, the predictive power of DFT crucially depends on the practical form of the XC functional approximation. The accuracy can be systematically improved by climbing up Perdew's metaphorical Jacob's ladder 28 of XC functionals, with each rung considering additional delicate physical effects. Climbing up Jacob's ladder, which consists of five rungs, to study water is a nontrivial task as more accurate XC functionals require significantly more computational resources. Water is a collection of molecules that are loosely bonded by positively charged H atoms pointing towards negatively charged oxygen (O) atoms, thereby forming the H-bonds that govern many peculiar properties of water. Due to the delicate nature of the H-bond network, it has been widely accepted that the accurate prediction of water properties demands higher rungs of Jacob's ladder than ordinary materials. 13 The first rung of Jacob's ladder is the local density approximation (LDA), 29,30 which is parameterized based on the uniform electron gas and tends to minimize the inhomogeneity of electron density.
While the overestimation of chemical bonds by LDA was observed as a general trend, 31,32 the significantly overestimated H-bond strength means that LDA fails in even qualitatively predicting water properties. Also, first-principles MD calculations of water by LDA only lasted for a short time period in the 1990s. In the following twenty years, simulations of water were dominated by utilizing the second rung of XC functionals with the generalized gradient approximation (GGA). 33,34 By introducing the electron density gradient, the GGA functional allows more constraints to be satisfied than LDA. In particular, the nonempirical GGA functional by Perdew, Burke, and Ernzerhof (PBE) 35 satisfies both the constraints on the electron-correlation hole and electron-correlation energy. PBE-GGA weakens the Hbond strengths towards the experimental direction. 13 Revisions of PBE such as revPBE 36 are constructed by satisfying fewer physical constraints than PBE. However, GGA also has limitations. PBE-GGA intrinsically lacks the intermediate-and long-range van der Waals (vdW) interactions and largely underestimates the liquid water density; notably, incorrectly yielding denser ice than water. 37 Adding explicit dispersion corrections, such as Grimme (D3) 38,39 and Tkatchenko-Sheffler (TS) 40 vdW corrections, to the GGA functional can significantly soften the structure of water and improve the density towards the experimental direction. 11,13,41 Recently, the non-empirical meta-GGA strongly constrained and appropriately normed (SCAN) functional 42,43 was proposed as a general-purpose ab initio functional that belongs to the third rung of Jacob's ladder. With the introduction of kinetic energy density together with a dimensionless variable, SCAN satisfies 17 known constraints on XC functionals as well as the tight lower bound on the exchange energy which has been elusive in PBE-GGA on the second rung of Jacob's ladder. SCAN was found to provide a substantially more accurate description of water as compared to GGAs, including structural, electronic, and dynamical properties. 14,44,45 Most importantly, SCAN includes intermediate-range vdW interactions inherently and the resulting density of water is 1.050 g/cm 3 , which is much closer to the experimental value of 0.997 g/cm 3 , as compared to the density of ∼ 0.850 g/cm 3 from PBE-GGA functionals. 14 Besides the accurate prediction of the density difference between water and hexagonal ice at ambient conditions, 14 the recent combination of the SCAN functional with a machine learning potential found that the SCAN functional can also correctly predict the relative stability between hexagonal ice and cubic ice. 46 In the past few years, SCAN has been one of the most commonly used non-empirical XC functionals in the studies of water.
Notwithstanding the above significant progress, challenges remain. Both GGA and meta-GGA functionals suffer from the self-interaction error. 47,48 The spurious interaction of the electron with itself leads to excessive proton delocalization and an artificially strengthened H-bond. 49 The overestimated H-bond strength is evidenced by the red-shift of the IR stretching band predicted by SCAN in comparison with experiment. 50 The above is also consistent with the slightly overstructured oxygen-oxygen pair distribution predicted by SCAN as compared to the neutron scattering data. 14 Recently, Lambros et al. assessed the accuracy of the SCAN functional by investigating the effects of the exact exchange using the MB-pol approach and found that adding different fractions of exact exchange to the SCAN functional leads to substantial changes in the structure of liquid water. 51 Furthermore, most AIMD simulations of water so far have assumed classical nuclei. However, the light hydrogen atom deviates from classical behavior significantly, even at room temperature. [52][53][54] The Feynman discretized path-integral (PI) method combined with AIMD can treat the nuclear quantum effect (NQEs). 55,56 But for reasons of computational cost, the PI-AIMD simulations have not yet been widely applied to water. More disturbingly, many thermodynamic properties of water are difficult to predict by DFT based on currently available computational resources.
Standard AIMD simulations are typically carried out for a time scale of tens of picoseconds, and the system size is limited to a simulation box containing a few tens to a hundred or so water molecules. In water, the time for the H-bond to break and reform is typically at the picosecond scale. 57 Therefore, the H-bond structure in water can be faithfully predicted within tens of picoseconds by AIMD simulation. However, for many macroscopic properties, such as diffusivity and dielectric constant, the converged predictions require much longer simulation time and larger system size. In order to approach the thermodynamic limit, at least thousands of water molecules are necessary to be considered in the simulation box, which is simply impossible for AIMD.
To address the above issues, we ascend Jacob's ladder of DFT to the fourth rung, using the SCAN0 hybrid functional to study liquid water with NQEs considered simultaneously by the PI method. The predictions of converged thermodynamic properties are largely facilitated by employing state-of-the-art deep potential molecular dynamics (DPMD). [58][59][60] The many-body potential and interatomic forces needed for DPMD are obtained from a deep neural network that is trained with SCAN0 DFT data. DPMD has recently been demonstrated to be able to accurately model the structure of pure water with an accuracy at the DFT level and with the low computational cost of a classical force field. 58,61 For the electronic properties, we adopt the Deep Wannier 62 neural network model, which is trained with the positions of centers of the maximally localized Wannier functions 63,64 that are obtained from SCAN0 DFT calculations. We focus on the SCAN0 results at 330 K so these results can be compared directly with previous simulations of water with the SCAN functional, 14,50 which were also carried out at this temperature. These results reveal the physical corrections caused by including exact exchange and nuclear quantum effects. At the same time, results for SCAN0 at 300 K are also presented for a direct comparison with 298/300 K experimental results. Compared to the third rung of Jacob's ladder of SCAN, SCAN0 presents a better description of liquid water in terms of structure, dynamics, as well as electronic properties. By mixing 10% exact change as suggested for an earlier meta-GGA functional, 65 SCAN0 mitigates the self-interaction error. The less polarizable water together with the NQEs on the protons yield a liquid structure softened towards the experimental direction compared with SCAN. The predicted dynamic properties such as IR spectra and diffusivities and the electric properties such as dielectric constants are in good agreement with available experiments. Finally, we compute the quasiparticle band structure of liquid water based on Hedin's GW approximation 66 as implemented in the BerkeleyGW package. 67,68 The predicted density of states (DOS) and bandgap are in reasonable agreement with recent photoemission spectroscopy (PES) experiments.
To construct an accurate and transferable SCAN0 path-integral deep potential model with a minimal number of the expensive SCAN0 DFT data, an active machine learning procedure called deep potential generator (DP-GEN) 59 was adopted. The procedure is summarized as follows:
(1) Considering that the molecular configurations predicted by SCAN and SCAN0 are close to each other, 900 configurations were uniformly extracted from the 64-molecule 11 ps SCAN PI-AIMD trajectory reported in Ref. 69. The total potential energy E and ionic forces F i of each atom i of these configurations were calculated using Quantum ESPRESSO (QE) 70 with the SCAN0 XC functional. The Hamann-Schlüter-Chiang-Vanderbilt (HSCV) pseudopotentials 71,72 with an energy cutoff of 150 Ry were employed. The obtained E and F i together with atomic positions were adopted as the initial training data.
(2) The DeePMD-kit package 73 was used to train four deep potential models independently for 10 6 steps. The differences between the training processes of the four deep potential models were the initialization random parameters and the activation function. Two models were trained using the Tanh activation function, while the other two models were trained using the GELU activation function. The training followed the procedure described in Refs.
58,60. First, the input data were transformed to local coordinate frames for every atom and its neighbors inside a cutoff distance of 6 Å to preserve the translational, rotational, and permutational symmetries. Afterward, the deep neural network parameters were optimized by the Adam method 74 with the loss function:
where ∆ϵ and ∆F i represent the differences between the training data and current deep potential prediction for the quantities ϵ ≡ E/n and F i , respectively, n is the number of atoms, p ϵ and p f are tunable prefactors. In the training process, p ϵ progressively increases from 0.02 to 1, while p f progressively decreases from 1000 to 1.
(3) To explore the potential energy surface, ten independent PI-DPMD simulations were carried out in the NpT ensemble at 1 bar and 330/300 K for 100 ps with a timestep of 0.5 fs using the i-PI code 56 in conjunction with the DeePMD-kit package. 73 A periodically replicated cubic cell containing 64 water molecules was used. The nuclear degrees of freedom were sampled using eight beads with a colored-noise generalized Langevin equation thermostat 75,76 to accelerate the convergence of the quantum distribution. The pressure was controlled using an isotropic barostat as implemented in the i-PI code with a time constant associated with the dynamics of the piston of 200 fs. The cell was thermostatted by another generalized Langevin equation thermostat 77,78 to enhance sampling efficiency. For each PI-DPMD simulation, all of the four deep potential models were used to predict the force acting on each atom at each timestep but only one model randomly selected from the four models was used to propagate the trajectory. One converged deep potential model was employed to conduct the production run, which was trained using 7349 SCAN0 DFT configurations. Note that the number of training configurations adopted here is not the minimum number of configurations needed for convergence because, in the first few iterations, configurations that have ζ > 0.05 eV/Å were also added to the training data set. These "extra" configurations are not essential for the convergence and will not affect the simulation results. To verify the reliability of the deep potential model, we also conducted a SCAN0 AIMD simulation of liquid water in the N pT ensemble at 1 bar and 330 K for 20 ps with a 64-molecule cell using QE. HSCV pseudopotentials 71,72 with an energy cutoff of 130 Ry were adopted. As shown in Fig. 1, the RDFs predicted by AIMD and DPMD simulations agree well with each other, which proved the convergence of the SCAN0 deep potential model.
To predict the electronic properties, the deep Wannier 62 method was applied to build another deep potential model to map input atomic coordinates into output Wannier centers.
Firstly, the Wannier centers of the 7349 SCAN0 DFT configurations were calculated using QE. Each water molecule i has four valence Wannier centers. The centroid of the four Wannier centers r i w = 1 4 ∑ 4 c=1 r i wc and the atomic positions were adopted as the input training data and trained for 10 6 steps using the DeePMD-kit package. 73
The SCAN0 deep potential model was applied to conduct both DPMD and PI-DPMD simulations at 1 bar and 330/300 K with a simulation cell containing 512 water molecules.
Specially, the dielectric constant and diffusion coefficient were calculated using cells containing 64∼4096 molecules. All static properties were calculated using the N pT ensemble. To avoid the influence of the thermostat and barostat on the dynamics, all dynamical properties were calculated using the microcanonical (NVE ) ensemble with the cell size fixed at the value obtained from the N pT simulation. The DPMD simulation was carried out using LAMMPS 79 in conjunction with the DeePMD-kit package 73 for 2 ns with the first 100 ps discarded for equilibrium. The PI-DPMD simulation was conducted using the i-PI code 56 in connection with the DeePMD-kit package 73 for 500 ps with the first 100 ps discarded for equilibrium. The PI-DPMD simulation parameters related to nuclear and cell dynamics were the same as those adopted in step three of the deep potential training process.
The electronic DOS and bandgap of liquid water were calculated at the level of DFT and G 0 W 0 . All the DFT calculations reported in this work are carried out at the level of SCAN0.
A 0.1 eV Gaussian broadening factor was adopted for the computed eigenstates. The G 0 W 0 calculation was carried out by including dynamical dielectric screening effects using a generalized plasmon-pole model with the BerkeleyGW package. 67,68 Due to the large computational cost of the G 0 W 0 method, the SCAN0 quantum result was calculated using 8 snapshots uniformly extracted from the 500 ps 32-molecule SCAN0 quantum trajectory, with the PBE functional as a starting point. For all the configurations, the dielectric function ϵ and the self-energy Σ were calculated using 20000 bands to ensure convergence. The plane wave energy cutoff for ϵ was 30 Ry. A single gamma point was adopted. Overall, 128 occupied eigenstates and 5 empty eigenstates were calculated.
In this section, the structural, dynamical, and electronic properties of liquid water predicted by the SCAN0 functional at 330 K are presented and compared with previous SCAN AIMD simulations at 330 K to demonstrate the physical corrections caused by adding exact exchange and nuclear quantum effects to results for the SCAN functional. To show a direct comparison with the 298/300 K experimental results, the SCAN0 results at 300 K are also
Table 1: The structural, dynamical, and electronic properties of liquid water predicted by SCAN0 DPMD and PI-DPMD simulations at 330 and 300 K in comparison with the SCAN results at 330 K and the experimental results at 298 and 300 K (298 K for Refs. 9,19,20,80-83, 300 K for Refs. 84-86). The structural properties include the average O-H covalent bond length (d OH ), average number of H-bonds per water molecule (N HB ), and density (ρ). The dynamical properties include the peak positions of the four IR bands (ν T , ν L , ν B , ν S ), the diffusion coefficient (D), and the rotational correlation time (τ 2
). The electronic properties include the dipole moment (µ), the binding energies of 2a 1 , 1b 2 , 3a 1 peaks with respect to the 1b 1 peak in the electronic DOS, the bandgap, and the dielectric constant (ϵ). For the method column in the electronic part, the texts outside/inside the parentheses represent the calculation method of the molecular/electronic structure.
3 ) Structural Exp. 298/300 1.01 9 0.997 84 SCAN AIMD 330 0.984 14 3.61 14 1.050 14 SCAN0 DPMD 330 0.980 3.58 1.030 (H 2 O) SCAN0 PI-DPMD 330 0.999 3.49 1.041 SCAN0 DPMD 300 0.981 3.71 1.032 SCAN0 PI-DPMD 300 1.000 3.61 1.046 Method
ps) τ 2 (ps) Dynamical Exp. 298/300 186 19 486 19 1209 19 2498 19 0.20 85 2.4 86 SCAN AIMD 330 172 50 483 50 1207 50 2448 50 (D 2 O) SCAN0 DPMD 330 227 486 1219 2507 0.26 1.9 SCAN0 DPMD 300 213 514 1222 2479 0.12 4.6 Dynamical Exp. 298/300 176 20 615 20 1646 20 3406 20 0.24 85 1.7∼2.6 80,81 SCAN 330 0.26 87 3.3 87 (H 2 O) SCAN0 DPMD 330 208 675 1667 3475 0.29 1.7 SCAN0 DPMD 300 207 729 1669 3408 0.13 4.1 Method T (K) µ (Debye) ϵ 2a 1 -1b 1 (eV) ϵ 1b 2 -1b 1 (eV) ϵ 3a 1 -1b 1 (eV) bandgap (eV) ϵ Electronic Exp. 298/300 2.9±0.6 82 -19.7 88 -6.2 88 -2.3 88 8.7±0.5 89 78.39 83 SCAN AIMD (DFT) 330 2.97 14 -18.9 14 -5.7 14 -2.0 14 SCAN0 DPMD (DFT) 330 2.96 -19.4 -5.8 -2.0 5.1 76.83 SCAN0 PI-DPMD (DFT) 330 3.08 -19.4 -5.6 -2.0 4.6 84.34 (H 2 O) SCAN0 PI-DPMD (G 0 W 0 ) 330 -19.4 -6.4 -2.6 7.5 SCAN0 DPMD (DFT) 300 3.01 -19.4 -5.8 -2.0 5.1 95.33 SCAN0 PI-DPMD (DFT) 300 3.13 -19.4 -5.6 -2.0 4.6 97.57 SCAN0 PI-DPMD (G 0 W 0 ) 300 -19.4 -6.4 -2.6 7.5 presented in the table and the figures, and are discussed in the last subsection.
Water structure and density predicted at 330 K
We first inspect the structure of water predicted by various DFT functional approximations.
To this end, we analyse the radial distribution functions (RDFs), which show the probability to find a given pair of atoms as a function of distance in real space. The resulting oxygenoxygen and oxygen-hydrogen RDFs, g OO (r) and g OH (r), are shown in Figs. presented in the same figures with SCAN AIMD 14 and experimental results 9 for comparison. Like all local and semi-local DFT functionals, the SCAN functional inherits the spurious self-interaction error, which yields an artificially strengthened tetrahedral structure and delocalized protons. By mixing a fraction (10%) 65 of exact exchange in SCAN0, the above self-interaction error is mitigated which softens the liquid structure towards the experimental direction. The weakened H-bond strength by SCAN0 DPMD relative to SCAN AIMD is evidenced by the right-shift of the first peak position of g OO (r) in Fig. 2 (a), from 2.748 to 2.754 Å and the second peak position of g OH (r) in Fig. 2 (b), from 1.785 to 1.793 Å, both of which are closely associated with the H-bonds. The weaker H-bonds, in turn, enhance the covalency of water molecules in the liquid, as indicated by the shortened OH covalent bond length from 0.984 Å (SCAN AIMD) to 0.980 Å (SCAN0 DPMD). Not surprisingly, the weakened H-bonding also promotes the population of broken H-bonds. Adopting the 69,91 as shown in Figs. 2 (c) and (d), respectively. It can be seen in Fig. 2 (c) that the broken H-bond by the SCAN0 functional originates from the reduced proton displacement along the H-bonding direction. In contrast, from SCAN to SCAN0 functional, the proton displacement along the libration direction barely changes as displayed in Fig. 2 (d). This is consistent with the fact that the hybrid functional reduces the electronegativity of the lone pair electrons, which weakens the directional H-bonding strength. In accordance with the reduced directional H-bond strength, the water molecules are more loosely bonded with an increased O-O distance. Therefore, the predicted water density by SCAN0 DPMD is 1.030 g/cm 3 , which is noticeably smaller than the density of 1.050 g/cm 3 predicted by SCAN AIMD simulation. 14 In SCAN0 PI-DPMD simulations, the NQEs are incorporated in the modeling of water.
The resulting g OO (r) and g OH (r) are presented in Figs. 2 (a) and (b), respectively. In quantum simulations, the configuration space explored by protons is largely extended, in particular for the regions that are inaccessible to classical nuclei. Compared to classical simulations, PI-DPMD simulation results in a significantly broadened first peak of g OH (r), which is a typical effect reported in all path-integral simulations of water. 53,54,69,[91][92][93] Besides the broadening effect, NQEs affect the H-bond network as well. Under the influence of NQEs, the protons are more delocalized along the direction of the stretching mode as well as the direction of the libration mode as shown in Figs. 2 (c 91 Because of the more polarizable nature of water molecules under NQEs, the H-bonding force along the stretching mode becomes more attractive, which prefers a slightly more compactly packed water structure. Indeed, the predicted water density by SCAN0 PI-DPMD is 1.041 g/cm 3 , which is slightly larger than 1.030 g/cm 3 as predicted by SCAN0 DPMD. The larger density under NQEs has also been widely reported in previous studies. 69,92,93 Dynamic properties predicted at 330 K The prediction of the quantum dynamic properties of water is still an open problem due to the application of the complex form of the real-time propagator in the current implementations of Feynman's path integral approach. 94 Therefore, the dynamical properties, including the IR spectra and the diffusivity presented in this work, are generated from classical simulations only.
IR spectroscopy is a powerful experimental technique to probe the vibrational properties of liquid water. The IR spectrum is sensitive to the underlying H-bond network in liquid water and is, therefore, widely adopted to study the microscopic structure of liquid water. 49,50,95 Within linear response theory, the IR absorption coefficient per unit length is given by the Fourier transform of the time correlation function of the overall dipole moment as following: 96
where n(ω) is the refractive index, V is the volume, and β = (k B T ) -1 , with k B and T being Boltzmann's constant and temperature, respectively. The dipole moment µ i of molecule i is calculated as
where r i O , r i H 1 , and r i H 2 is the nuclear coordinates of the ith molecule and r i wc (c = 1∼4) is the four Wannier centers assigned to
the ith water molecule in the system, which is obtained from the Deep Wannier method in this work. Therefore, a precise prediction of the IR spectra demands an accurate description of both the molecular configuration and its dynamical electronic response.
We present the IR spectrum of heavy water predicted by SCAN0 DPMD simulation in improves all the four spectral features towards the experimental direction, in terms of both the spectral energies and their intensities. In order to elucidate the physical origin of the improved spectrum, we further decompose the IR spectrum into its intra-and inter-molecular contributions. 95 The above can be rigorously implemented by separating the time-correlation function in Eq. 1 into contributions from i = j and i̸ =j terms utilizing the additive nature of the two-body correlation function. The resulting decompositions are shown in Fig. 3 (b).
We first inspect the stretching band centered around 2500 cm -1 . According to previous studies, it was found that SCAN functional improves the underestimated stretching frequency compared to the PBE functional. 50 However, the predicted center of the stretching band at ν S =2448 cm -1 50 is still 50 cm -1 smaller than the experimental value of 2498 cm -1 , 19 which indicates that the H-bond strength is still overestimated by SCAN functional. Here, the above overestimated H-bond strength is largely corrected by the SCAN0 functional. The SCAN0 functional yields a stretching band centered at ν S =2507 cm -1 , which is in quantitative agreement with the experimental measurement. We attribute the improvement of the IR stretching band to the reduced self-interaction error by mixing a fraction of exact exchange in SCAN0, which mitigates the overestimated directional H-bonding strength. The softened water structure reduces the tendency of a proton within a water molecule to be donated to the neighboring water molecules, and therefore hinders the vibration of protons as described by the stretching mode. Consistently, the decompositions of the IR spectra show that the improved stretching band originates from both the inter-and intra-molecular contributions, which show a comparable blue shift compared to that from the SCAN functional. This is because the less polarized H-bond network by SCAN0 makes a proton less likely to be donated to a neighboring water molecule and therefore increases the covalent bond strength as evidenced by the shortened OH bond length as shown in Table 1 or Fig. 2. Therefore, it is not surprising that the intra-and inter-molecular contributions to the stretching mode both become stiffened by the exact exchange in SCAN0.
We next draw our attention to the bending band of the IR spectrum, which is centered
Compared with the prediction by SCAN, the SCAN0 functional predicts an improved bending band. In particular, a significant enhancement can be identified in the intensity of the bending motion, which should be attributed to the improved H-bond network as well. Unlike the stretching band, the intra-and inter-molecular contributions to the bending mode are competing with each other. The intra-and inter-molecular dipoles have positive and negative correlations, respectively, which results in the opposite signs of the above two decompositions at roughly the same frequencies as shown in Fig. 3 (b). In crystalline ice, the intensity of the bending band is rather weak due to the strong intermolecular dipole-dipole correlation while the contrary is true in liquid water. 97,98 Due to the artificially strengthened H-bond network described by the SCAN functional, the intermolecular negative dipolar correlation is overestimated, which results in a weaker spectral intensity than that of experiment as shown in Fig. 3 (a). Whereas, the softened structure of liquid water predicted by SCAN0 suppresses the negative correlation of inter-molecular dipole moments and instead promotes the positive intra-molecular dipole correlation. The above effect gives an important correction to the bending band of the IR spectrum as seen by the increased overall intensity towards the experimental direction.
Finally, we analyze spectral features in the low-frequency range, which are composed of the hindered translation band centered at ∼ 200 cm -1 and the libration band center at ∼ 500 cm -1 . It can be seen that the IR spectrum predicted by SCAN0 also shows better agreement with experiment compared to SCAN. In particular, the improvement in the main peak of the libration band is more noticeable. The decompositions of the IR spectrum into its inter-and intra-molecular contributions, shows that the improvement originates from the intra-molecular part. The positive intra-molecular dipolar correlation is slightly enhanced under the weakened H-bond network by SCAN0, which hampers the water libration motion by shifting its band to a slightly higher frequency range.
Besides the improvement to the IR spectrum compared to SCAN, our SCAN0 DPMD simulation also captures the spectral difference between heavy and light water accurately as shown in Figs. 3 (a) and (c). Because of the lighter mass of H than D, all the four spectral features in light water are systematically higher than those in heavy water. A difference can be identified in the stretching/bending band between light water at 3475/1667 cm -1
and heavy water at 2507/1219 cm -1 , respectively, which originates from molecular vibration in the vapor phase. In the low frequencies range, the spectral differences are relatively smaller. The translation/libration band is located at 208/675 cm -1 for light water and 227/486 cm -1 for heavy water. Not surprisingly, the reduced isotope effects are due to the fact that collective vibrations become more important in the lower frequency region, which is less affected by the mass difference between H and D. The self-diffusion coefficient is an important physical quantity that differentiates liquid water from crystalline ice. It has been well recognized that the calculation of the diffusion coefficient of water is a challenge for AIMD. The difficulty lies in the fact that the predicted diffusivity of water, which is modeled by finite box size under periodic boundary conditions, must be corrected before a meaningful comparison can be made with the experimental value D ∞ , which corresponds to the diffusivity of water in its thermodynamic limit. Based on the Kirkwood-Riseman theory of polymer diffusion, the correction term adopts the form of 99 where D(L) is the diffusivity obtained in a system with cell length L, ζ is a numerical coefficient of 2.837, and η is the shear viscosity. In the above, D(L) approaches D ∞ rather slowly with a 1/L scaling. In practice, D ∞ is obtained by linear interpolation. Therefore, the accurate prediction demands water models of large box size containing thousands of molecules, which is computational challenging for DFT-based AIMD.
In this work, utilizing the SCAN0 deep potential model, a series of simulations with cell sizes ranging from 64 to 4096 water molecules were conducted for both light and heavy water.
For each cell size, four independent simulations were conducted for 300 ps. The diffusion coefficients were calculated by linear fittings of the mean square displacement of the molecular center of mass in the time interval of 1∼10 ps. The error bars were estimated as the root mean square error of the four diffusion coefficients obtained from the four independent trajectories. As shown in Fig. 4, the extrapolated values of D to infinite system size are 0.26 and 0.29 Å 2 /ps for D 2 O and H 2 O, respectively, which is in qualitative agreement with the 300 K experimental result of 0.20 and 0.24 Å 2 /ps for D 2 O and H 2 O. 85 It can be noticed that the diffusion coefficient 0.29 Å 2 /ps of light water predicted by the SCAN0 DPMD simulation at 330 K is slightly larger than the 0.26 Å 2 /ps predicted by the SCAN neural network 87 at 330 K. This is because the more weakly H-bonded water molecules described by the SCAN0 than that by the SCAN functional promotes the diffusion of water.
The diffusion coefficient displays the translational dynamics of the H-bonds, while the rotational correlation time reflects the rotational dynamics which can be calculated by the integration of the second order orientational correlation function as:
where P 2 is the second Legendre polynomial and u is the unit vector along the molecular vector of interest. In this work, u is chosen to be the unit vector along the OH covalent bonds of the water molecules. The values of τ 2 for heavy and light water obtained from the SCAN0 DPMD simulation are 1.9 and 1.7 ps, respectively, which is in qualitative agreement with the nuclear magnetic resonance experimental values of 2.4 ps 86 for heavy water and 1.7∼2.6 ps 80,81 for light water. Compared to the τ 2 of 3.3 ps for light water predicted by the SCAN neural network 87 at 330 K, the τ 2 of 1.7 ps predicted by our SCAN0 DPMD simulation indicates a faster rotational dynamics, which is in accordance with the faster translations dynamics by SCAN0 and can also be explained by the weaker hydrogen bonds due to the mitigation of self-interaction error.
Electronic properties predicted at 330 K Electronic DOS and bandgaps
) Electronic DOS of liquid water calculated using the DFT and G 0 W 0 methods based on the configurations predicted by SCAN AIMD, 14 SCAN0 DPMD, and SCAN0 PI-DPMD at 330 K as well as the experimental photoemission spectra. 88 The 1b 1 peak position is used to align all calculations. The inset shows the linearly extrapolated threshold of the valence band maximum. (b) Comparisons between the electronic DOS calculated using the configurations predicted at 330 and 300 K.
The electronic structure of liquid water can be efficiently probed by spectroscopy experiments that involving single-particle excitations. The valence band and conduction band can be extracted from the experimental data obtained by direct and indirect PES measurements, 89 respectively. Moreover, the energy difference between the conduction band minimum and valence band maximum gives the quasiparticle (QP) bandgap. Physically, PES measures the energy of adding or removing a single-particle excitation from the system, which can be understood as electron or hole quasiparticles dressed by interactions with the system. In practice, these quasiparticles can be described by Hedin's GW approximation 66 within many-body perturbation theory. Within the GW approximation, the QP energies are typically calculated as one-shot correction without self-consistency (i.e., G 0 W 0 ) on top of a DFT mean field, making electronic structure calculations within GW an important benchmark for the accuracy of the DFT electronic ground state.
In Fig. 5, we present the DOS predicted by both DFT and G 0 W 0 calculations, with the molecular configurations obtained from SCAN AIMD, 14 SCAN0 DPMD and SCAN0 PI-DPMD simulations. The experimental PES result 88 is also presented for comparison. In Fig. 5, the four peaks with energies smaller than zero belong to the valence band and can be assigned to the 2a 1 , 1b 2 , 3a 1 , and 1b 1 orbitals based on the spatial symmetries of a water molecule, while the peaks having energies above the Fermi level belong to the conduction bands. The DFT results in Fig. 5 show that, as the functional changes from SCAN to SCAN0, the energy of the 2a 1 orbital with respect to the 1b 1 orbital changes from -18.9 to -19.4 eV, the latter of which is closer to the experimental result of -19.7 eV. Therefore, SCAN0 provides a more precise description of the strongly bound 2a 1 peak than SCAN.
The above can be explained by the fact that the 2a 1 orbital is primarily composed of the 2s electrons of the oxygen atom, which are corrected significantly when the self-interaction error is mitigated by the exact exchange in SCAN0. The decrease of the 2a 1 energy explains the increased covalency by SCAN0, as the bonding pairs of electrons in a water molecule are primarily come from a combination of the 1b 2 and 2a 1 orbitals with a small proportion of the 3a 1 orbital. 93 When NQEs are incorporated by the PI method, the distribution of the DOS is broadened and closer to the experimental result because NQEs enable the protons to explore more configuration space that is inaccessible in classical simulations. The broadening of the 2a 1 and 1b 2 orbitals are more obvious than the broadening of the 3a 1 and 1b 1 orbitals because the 2a 1 and 1b 2 orbitals are highly correlated with the bonding pairs of electrons and the hydrogen atoms in water molecules are delocalized significantly under NQEs. It needs to be noted that even with the broadening by NQEs, the calculated 2a 1 peak is still narrower than the experimental result as shown in Fig. 5 (a). This is because the energy loss of the initial 1b 1 electrons results in a broad background of the experimental 2a 1 feature, 12 which is not accounted for in our calculations.
For the SCAN0 PI-DPMD configurations, when the calculation method of the electronic structure changes from DFT to G 0 W 0 , the DOS associated with the valence states barely changes (the small differences between theses two results are within the margin of error associated with the statistical fluctuations, as G 0 W 0 results are only averaged on eight snapshots), but the bandgap at the G 0 W 0 level is significantly improves compared to DFT. While Kohn-Sham DFT is a formally exact theory for the ground-state total energy and density of a material, a Kohn-Sham electronic structure calculation at fixed electron number fails to capture the contribution to the fundamental band gap from the discontinuity of the exact exchange-correlation potential under a change of electron number. 100,101 As a result, it is well-known that the bandgaps in insulating systems are underestimated. 102,103 The underestimations are noticeable in DFT calculations by LDA and GGA functionals and become less severe when the nonlocal XC effects are considered in generalized Kohn-Sham such as hybrid functionals with exact exchange. 104 The above is well illustrated by the comparison of the bandgap of 4.9 eV predicted by the meta-GGA SCAN functional 14 with the bandgap of 5.7 eV predicted by the hybrid SCAN0 functional, the latter of which is closer to the experimental result of 8.7±0.5 eV. 89 It needs to be mentioned that the above bandgaps of 4.9 14 and 5.7 eV are obtained by simply taking the energy difference between the highest occupied orbital and lowest unoccupied orbital. In practice, a better way to calculate the bandgap is determining the position of the valence band maximum and conduction band minimum as the linearly extrapolated threshold as shown in the inset of Fig. 5, which leads to faster convergence of the valence band maximum and conduction band minimum values with respect to cell sizes. 12,105 In this way, the bandgaps predicted by DFT based on the SCAN0 DPMD and PI-DPMD configurations are 5.1 and 4.6 eV, respectively. The quantum simulation reduceds the bandgap, in agreement with previous findings. 12 The G 0 W 0 calculation corrects the bandgap from 4.6 to 7.5 eV, which is close to the experimental result of 8.7±0.5 eV. 89 Further improvement in the bandgap can be achieved by the quasiparticle self-consistent GW method with vertex corrections. 12 However, that is beyond the scope of this work. The large dielectric constant is another anomalous property of liquid water, whose value is significantly larger than other materials with comparable electric dipoles. 83 A microscopic insight into the dielectric constant of liquid water is essential for understanding its properties as a solvent. However, the ab initio calculation of the dielectric constant, ϵ, of liquid water has long been a challenging prediction. Under the periodic boundary condition, the dielectric constant is connected with the fluctuations of electric dipole moments of water molecules and their correlations. Based on linear response theory, the dielectric constant of water can be calculated by the following equation:
where ϵ ∞ = 1.77, µ is the molecular dipole moment, ρ is the molecular number density, k B is the Boltzmann's constant, and T is the temperature. G k is the large-r limit of the correlation factor G k (r), which contains the integration of correlations among the dipole moments of different molecules and can be calculated as: The SCAN0 classical and quantum simulation yield a dielectric constant of 76.83 and 84.34, respectively, which is close to the experimental result of 78.39. 83 As shown by the factor ⟨µ 2 ⟩G k (r) in Eq. 2, the large dielectric constant of liquid water is contributed by both the molecular dipole moment and G k (r). The value of G k (r) predicted by SCAN0 classical and quantum simulations are 2.68 and 2.67, respectively. While the value of ⟨µ 2 ⟩ obtained from SCAN0 classical and quantum simulations are 8.83 and 9.66 Deybe 2 , respectively.
The structure of liquid water predicted by the SCAN0 functional at 300 K is more overstructured compared to that at 330 K as shown by the RDFs in Fig. 2. This is because the weaker thermal fluctuations at a lower temperature break fewer hydrogen bonds, which results in a more tetrahedral-like liquid structure. The differences between the RDFs predicted by the SCAN0 functional at 300 K with the experimental results at 298 K indicate that the hydrogen bond strength is still overestimated by the SCAN0 functional. The stronger hydrogen bonds at 300 K slow down the diffusion of the water molecules as illustrated by the smaller diffusion coefficient and larger rotational correlation time in Table 1. Compared to the IR spectra predicted at 330 K, the 300 K results agree better with the 298 K experimental spectra as displayed in Fig. 3. Specifically, the vibrational frequency of the O-H/O-D stretching mode decreases due to the strengthening of the hydrogen bond; the intensity of the IR spectra increases because the intensity is calculated using dipolar correlation dividing the temperature as shown by Eq. 1. The dielectric constant increases when the temperature decreases as shown in Table 1, which is caused by the increased molecular dipole moment, the stronger correlations between the dipoles of different molecules, and the decreased temperature (see Eq. 2).
The differences between the 330 K classical results and the 300 K quantum results predicted by the SCAN0 functional suggest that increasing the simulation temperature by 30 K does not have the equivalent effect of including nuclear quantum effects, which agrees with the findings in Refs. 87,109. The overstructured H-bond network of liquid water predicted by the SCAN0 functional at 300 K reveals clearly the deficiency of the SCAN0 functional. Compared to DFT functionals at the lower rung of Jacob's ladder of the GGA or meta-GGA level, the improvement of predicted water properties by the SCAN0 functional are qualitatively similar to those by revPBE0 functional with a D3 dispersion correction 92 and the PBE0 functional with a TS vdW correction. 61 However, the construction of SCAN0 functional is motivated by satisfying as many physical constraints as possible in one single functional and simultaneously avoiding many empirical parameters, which increases its transferability within the framework of DFT.
In conclusion, we have investigated liquid water based on DFT, in which the hybrid meta-GGA functional, SCAN0, was employed to approximate the exchange-correlation effect among electrons. The molecular structure of the liquid, the dynamics of the H-bond network, and the electronic properties of water were predicted and compared to available experimental measurements. In particular, we applied deep neural networks as implemented in the DeePMD and Deep Wannier methods, which were trained on the potential energy surface and electronic structure obtained for the electronic ground state described by SCAN0-DFT.
The applications of the deep neural network potentials enable us to model water in much larger simulation boxes with significantly longer simulation time than those in regular DFT, which thereby enabled prediction of the macroscopic water properties by approaching the thermodynamic limit.
In addition to the existing 17 constraints satisfied by the SCAN functional, the nonlocal exact exchange mixed in SCAN0 further alleviates the self-interaction error. The reduced molecular polarizability using SCAN0 softens the H-bond network of water towards the experimental direction as evidenced by the comparison of predicted RDFs with those obtained in scattering experiments. Since SCAN0 water is more loosely bound by the H-bond network, the water molecules are more diffusive and less correlated in spatial and temporal dimensions. The improvements in water dynamics have been demonstrated in IR spectra and diffusivity by comparison between theory and experiments. The spatial correlations, as represented by the dipole-dipole correlation in the dielectric response, were also accurately determined. Moreover, the predicted static dielectric constant is in good agreement with experiments. The electronic band structure and band gap of water were determined by the GW self-energy approach as implemented in the Berkeley GW code package, which agree well with the PES experimental measurements. In all predicted static properties, the effects of including NQEs were found to give rise to mainly a broadening effect due to the delocalized protons as widely reported in the literature.
Our studies present the state-of-the-art theoretical modeling of water on the fourth rung of Jacob's ladder. Although there is improved accuracy in predicted water properties, some important discrepancies remain. The differences between the simulation results at 300 K with the 298/300 K experimental results indicate that the predicted H-bond network by SCAN0 is still slightly overstructured. The water density predicted by SCAN0 is improved over that by SCAN, but it is still a little overestimated. To address these remaining issues, more sophisticated treatment of the self-interaction correction and long-range van der Waals interactions should be included in DFT along with nonlocal exchange-correlation functional constructions at higher rungs of Jacob's ladder.
the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S.