We perform ab initio molecular dynamics (AIMD) simulation of liquid water in the canonical ensemble at ambient conditions using the strongly constrained and appropriately normed (SCAN) meta-generalized-gradient approximation (GGA) functional approximation and carry out systematic comparisons with the results obtained from the GGA-level Perdew-Burke-Ernzerhof (PBE) functional and Tkatchenko-Scheffler van der Waals (vdW) dispersion correction inclusive PBE functional. We analyze various properties of liquid water including radial distribution functions, oxygen-oxygen-oxygen triplet angular distribution, tetrahedrality, hydrogen bonds, diffusion coefficients, ring statistics, density of states, band gaps, and dipole moments. We find that the SCAN functional is generally more accurate than the other two functionals for liquid water by not only capturing the intermediate-range vdW interactions but also mitigating the overly strong hydrogen bonds prescribed in PBE simulations. We also compare the results of SCAN-based AIMD simulations in the canonical and isothermal-isobaric ensembles. Our results suggest that SCAN provides a reliable description for most structural, electronic, and dynamical properties in liquid water.

Liquid water is one of the most important solvents for chemical and biological processes.1 Bonded by relatively strong covalent bonds, the structure of a single water molecule is now well understood. However, in its condensed phases, water molecules form tetrahedral or near-tetrahedral structures through hydrogen bonds (H-bonds). Within such a network, water molecules are bonded to each other through directional H-bonds resulting from weak electrostatic attractions between a lone pair electron of oxygen (O) and hydrogen (H) from the neighboring water molecule. In contrast to covalent bonds, the H-bonds are energetically much weaker. Therefore, the structure of the H-bond network is not only easily disturbed by thermal fluctuations and applied pressures as evidenced by the richness of ice phases and liquid structures2–4 but also largely affected by relatively weak physical interactions such as van der Waals (vdW) dispersion forces.5–7 Furthermore, the above picture becomes even more complex when the quantum nature of the nuclei is taken into consideration.8–10 Because of the delicacy of the H-bond network, the structure of liquid water has neither been concluded from experiment nor from theory. Therefore, the nature of the H-bond structure in liquid water continues to be at the center of scientific debate. Not surprisingly, research in this field is progressing by continuous joint efforts from both experiment and theory.

Several experimental techniques were used to probe the microscopic properties of liquid water, including X-ray and neutron scattering,11 X-ray absorption spectroscopy (XAS),12 infrared spectroscopy,13 Raman spectroscopy,14 and photoemission spectroscopy (PES),15,16 etc. X-ray and neutron scattering experiments directly measure the statistically averaged structural information. On the one hand, X-ray is more sensitive to O than H atoms since it is scattered by electrons and can only capture O–O correlations with high clarity. On the other hand, neutrons scattered from nuclei are ideal for obtaining correlations involving H atoms. The scattered intensities are used to extract the O–O, O–H, and H–H radial distribution functions (RDFs). However, experimental determination of the RDFs is very challenging and typically requires scattering experiments with isotope substitutions and theoretical input such as empirical potential structural refinement (EPSR).10–12,17 In addition, XAS measurements are able to probe the electronic structure in different phases of water, from which the tetrahedral H-bonding and shell structure information are indirectly obtained.12 Furthermore, the infrared and Raman spectroscopies allow the extraction of the O–H vibrational dynamics of water.13,14 Finally, the valence band photoemission experiments measure the occupied electronic states, which reveal the local bonding environment in liquid water.15,16

Nevertheless, due to breaking and reforming of H-bonds on the time scale of picoseconds (ps),18 an experimental technique that could directly detect the real-space microscopic structure of liquid water is still lacking. Importantly, structural properties of liquid water could be obtained via computer simulations. For example, ab initio molecular dynamics (AIMD)19 simulations employing density functional theory (DFT)20 model a given system based on quantum mechanical principles and in general need no empirical inputs. Specifically, AIMD simulations generate a nuclear potential energy surface “on the fly” from the electronic ground state and yield a variety of properties including structural, electronic, and dynamical properties of the system. Therefore, AIMD is an ideal approach for modeling important properties of water across the phase diagram.21–24 Indeed, during the past two decades, much progress has been made on the AIMD simulations of liquid water by adopting different levels of functional approximations.25–33 These studies have advanced our understandings on various properties of liquid water and the ice phases. Other methods have been also applied to liquid water in recent years, e.g., Monte Carlo simulations,34,35 Møller-Plesset perturbation theory (MP2),36 and random phase approximation.37 However, DFT is still the most popular ab initio method to study liquid water due to the balance between accuracy and computational costs.

However, challenges still remain in providing a more accurate water model by DFT-based AIMD simulations, which depend crucially upon the underlying exchange-correlation (XC) functional. For instance, gas-phase water clusters were first modeled by the local density approximation (LDA)38 with apparently over-estimated H-bond strengths and distances between water molecules that are too close. Currently, one of the most widely used category of XC functionals for AIMD simulations of liquid water is the generalized-gradient approximation (GGA), which largely corrects the over-binding provided by the LDA functional but still largely deviate from the experiments. Among various GGA-level XC functionals, the Perdew-Burke-Ernzerhof (PBE)39 and Becke-Lee-Yang-Parr (BLYP)40,41 functionals are two of the most popular forms.25–27,29,42 However, the modeling of liquid water by GGA still exhibits a rigid solvation shell structure and much smaller diffusivity when compared to experimental data. Furthermore, the inherent self-interaction error (SIE) existing in the GGA-level functionals unavoidably leads to excessive proton delocalization in liquid water and overestimation of the molecular polarizability.43,44 What is worse, GGA functionals fail to describe the correct density ordering between liquid water and hexagonal ice, leading to the consequence that the GGA ice would sink in GGA water.5,27,42 The above is due to the fact that GGA ignores the non-local electron correlation effects that are responsible for vdW or dispersion interactions, which are crucial for improving water density properties.5,27,45,46 Although the short-range portion of the vdW interactions is described by the local and semi-local XC functionals, the intermediate- and long-range parts are not captured by any general-purpose GGA functionals. Besides the choice of XC functionals, the nuclear quantum effects (NQEs) are essential for light atoms such as hydrogen and deuterium, which are not considered in the framework of classical molecular dynamics. Studies8,9 have shown that the inclusion of NQEs in AIMD simulations further softens the shell structure of water molecules in the liquid phase. Therefore, an accurate description of the balance among covalent bonds, H-bonds, and vdW interactions in liquid water that could eventually improve its H-bond network remains a challenge.

Several studies have tried to address the issues in different ways. For instance, to reconcile the absence of intermediate- and long-range vdW interactions, several empirical47 and non-empirical48 vdW corrections have been proposed. Particularly, the vdW density functional5,49 and density-dependent Tkatchenko-Scheffler dispersion correction6,27,28,50 were found to soften the solvation shell structures of water molecules and largely correct the density of liquid water. Moreover, the hybrid XC functional such as PBE0, wherein a fraction of the exact exchange energy is included in the density functional approximation, was used for the simulation of liquid water in order to mitigate the SIE.27,28,51 Properties of liquid water with the use of the PBE0 and Tkatchenko-Scheffler dispersion correction were thoroughly studied in Ref. 28, and it was found that the substitution of PBE by PBE0 functional in the canonical ensemble has similar effects as vdW interactions which soften the shell structures of water molecules. However, by only including PBE0, the density of liquid water is further lower compared to that from PBE.27 In summary, although simulations of liquid water utilizing GGA XC functionals exhibited several critical shortages on the description of the H-bond network in liquid water, other non-empirical, general-purpose XC functionals that could describe all types of interactions in liquid water on an equal footing are still largely absent. Instead, one could only combine different XC functionals, most of which include empirical parameters, to improve the description of liquid water at the molecular level.

The recently developed strongly constrained and appropriately normed (SCAN) meta-GGA functional52 can be used within the AIMD approach to address the above issues. SCAN satisfies all 17 exact constraints appropriate for semilocal XC functionals and has been shown to predict accurate geometries and energies of diversely bonded molecules and condensed matter materials while maintaining the computational cost similar to GGA.52,53 Furthermore, Ref. 53 reports that the SCAN functional is highly accurate for covalent bonds in water monomers; it also captures the intermediate-range vdW interactions in ice and gas-phase water hexamers. This study implies that SCAN possesses the ingredients necessary to describe liquid water. In a recent work,54 AIMD simulations of liquid water in the isothermal-isobaric (NpT) ensemble employing the SCAN functional were carried out, reporting that SCAN correctly predicts higher density of liquid water than ice and gives a delicately balanced description of covalent bonds, H-bonds, and vdW interactions in liquid water.

In this work, in order to compare the different vdW effects in the SCAN and PBE+vdW functionals, as well as the effects of different ensembles on properties of liquid water, we performed AIMD simulations of liquid water in the canonical (NVT) ensemble with PBE, PBE+vdW, and SCAN functionals. We demonstrate in detail that a variety of properties of liquid water by SCAN are systematically improved toward the experimental direction when compared to the semi-local PBE and PBE+vdW functionals. Essentially, SCAN improves the description of the electronic structure of liquid water and captures the intermediate-range vdW interactions. As a result, structural properties including radial distribution functions, local tetrahedrality, oxygen-oxygen-oxygen angular distribution function, and statistics of H-bond network properties, as well as electrostatic properties such as molecular dipole moment and band gap in liquid water have been systematically improved. We find that both SCAN and PBE+vdW functionals outperform the PBE functional. Moreover, although vdW interactions can be described by both SCAN and PBE+vdW functionals, the two functionals perform differently; SCAN has a more significant improvement over the structural properties of liquid water than PBE+vdW does. This is understood by the fact that the vdW functional does not affect the directional H-bond strengths, while SCAN improves the covalent bonds in water molecules and reduces the overly strong directional H-bond strength as described by the PBE functional. Besides, we also compare the simulation results between the NVT and the NpT ensembles, which illuminate how ensembles and slightly different water densities influence the structural properties of liquid water.

The paper is organized as follows. First, we describe the simulation methods and list the computational details in Sec. II. We then show the AIMD simulation results of structural, electronic, and dynamical properties of liquid water as obtained by employing different XC functionals in Sec. III. Finally, we conclude the results in Sec. IV.

All simulations were performed with the Car-Parrinello (CP) molecular dynamics19 implemented in the Quantum ESPRESSO package.55,56 We set up a cubic cell with L = 12.445 Å and periodic boundary conditions which consists of 64 H2O molecules. The electronic wavefunctions were expanded using a plane wave basis set with an energy cutoff of 85 Ry. The norm-conserving pseudopotentials57 in the form of Hamann-Schlüter-Chiang-Vanderbilt (HSCV) were adopted for the O and H atoms. The functionals employed herein include the semi-local PBE GGA functional, PBE with the non-local vdW/dispersion interactions in the form of the density-dependent Tkatchenko-Scheffler50 dispersion correction, i.e., PBE+vdW, and the SCAN meta-GGA functional. The parameters used for the Tkatchenko-Scheffler vdW were adopted from Ref. 28.

All AIMD simulations were performed in the NVT ensemble at 330 K by following Ref. 28. The PBE, PBE+vdW, and SCAN systems were initially equilibrated for about 8 ps and then continued for 55, 54, and 50 ps, respectively, for data collection. The nuclear mass of deuterium, 2.014 amu, was used for each hydrogen atom, and the mass of oxygen is 15.999 amu. The ionic temperature was controlled by a single Nosé-Hoover thermostat58–60 with a frequency of 60 THz. The CP equations of motion were integrated using the standard Verlet algorithm with a time step of 2.0 a.u. (∼0.0484 fs). To ensure an adiabatic separation between the nuclear and electronic degrees of freedom in the CP dynamics, an appropriate electronic fictitious mass μ should be chosen in CP molecular dynamics. It was discussed in Ref. 26 that compared to Born-Oppenheimer (BO) dynamics, simulations of liquid water with CP dynamics yield lower maximums and higher minimums in the radial distribution functions. In particular, a larger μ used in CP dynamics further softens the structure of liquid water. In this work, we set the electronic fictitious mass μ to be 100 a.u. with an electronic mass cutoff of 25 a.u.61 Additionally, we compare the NVT data to those obtained with CP molecular dynamics simulations employing PBE and SCAN functionals in the NpT ensemble, and the parameters of NpT simulations are described in Ref. 54.

We compare the computed results of radial distribution functions with the experimental results. In each of the O–O, O–H, and H–H radial distribution functions, we utilize the experimental data obtained from the EPSR based on joint X-ray/neutron data by Soper and Benmore10 and compare them with our simulation results performed with PBE, PBE+vdW, and SCAN functionals. In general, we find that while both PBE+vdW and SCAN functionals soften the structure of liquid water compared to that obtained from the PBE functional, SCAN has a larger correction than PBE+vdW does, and the reasons are explicated in the following paragraphs.

Figure 1 shows the O–O radial distribution function, gOO(r), of liquid water. The intensities and positions of maximums and minimums of gOO(r) found from our simulations, together with those measured by the X-ray diffraction experiments with high Q-range and low statistical noise,17 are shown in Table I. At the PBE level of theory, due to the lack of intermediate- and long-range vdW interactions and the SIE that lead to overly strong H-bonds, the liquid water system has a relatively more rigid shell structure compared to experiment. By including vdW correction in the PBE simulation, the structure of liquid water is improved, which is consistent with previous reports.5,27,28,45,51 The gOO(r) is softened because of the non-directional vdW interactions which provide attractive forces between a given molecule and the molecules around it. Therefore, the water molecules originally located at the second coordination shells move inwards to the interstitial region between the first and second shells. Meanwhile, some water molecules in the first shell are pushed outwards due to the breaking of H-bonds. As a result, the locally disordered configurations with O–O distance in the interstitial region are energetically more favorable. In other words, the inclusion of vdW interactions corrects the overly strong H-bonding network by populating the interstitial region and softening the over-structured configurations. However, since vdW interactions alone do not change the description of electronic structure of liquid water, the overly strong directional H-bonding strength prescribed in the PBE density functional remains unchanged.

FIG. 1.

O–O radial distribution functions gOO(r) of liquid water as obtained from AIMD simulations in the canonical ensemble employing PBE (red solid line), PBE+vdW (black dashed line), and SCAN exchange-correlation functionals (blue dashed line). The experimental data are from joint X-ray/neutron experiments.10 

FIG. 1.

O–O radial distribution functions gOO(r) of liquid water as obtained from AIMD simulations in the canonical ensemble employing PBE (red solid line), PBE+vdW (black dashed line), and SCAN exchange-correlation functionals (blue dashed line). The experimental data are from joint X-ray/neutron experiments.10 

Close modal
TABLE I.

Tabulated summaries of the structural properties of liquid water from the O–O radial distribution functions gOO(r) obtained from AIMD simulations in the canonical ensemble and experimental data. Three different XC functionals, namely, PBE, PBE+vdW, and SCAN are used in AIMD simulations. From left to right, positions (in Å) and intensities of the first maximum (r1max and g1max), first minimum (r1min and g1min), and second maximum (r2max and g2max) of the corresponding gOO(r).

gOO(r)r1maxg1maxr1ming1minr2maxg2max
PBE 2.71 3.59 3.28 0.31 4.44 1.55 
PBE+vdW 2.72 3.24 3.25 0.51 4.38 1.38 
SCAN 2.75 2.96 3.36 0.76 4.47 1.15 
Expt.17  2.80 2.57 3.45 0.84 4.50 1.12 
Expt.10  2.75 2.62 3.45 0.84 4.43 1.13 
gOO(r)r1maxg1maxr1ming1minr2maxg2max
PBE 2.71 3.59 3.28 0.31 4.44 1.55 
PBE+vdW 2.72 3.24 3.25 0.51 4.38 1.38 
SCAN 2.75 2.96 3.36 0.76 4.47 1.15 
Expt.17  2.80 2.57 3.45 0.84 4.50 1.12 
Expt.10  2.75 2.62 3.45 0.84 4.43 1.13 

By using the SCAN XC functional, we observe that the gOO(r) of liquid water is further softened as compared to the gOO(r) from the PBE+vdW AIMD simulation. Specifically, as shown in both Fig. 1 and Table I, the intensities of the first maximum (g1max), first minimum (g1min), and second maximum (g2max) are all improved; additionally, SCAN corrects the position of the first minimum, r1max, toward the experimental value. The excellent performance by SCAN to soften the structure of liquid water is partly attributed to its ability to capture the intermediate-range vdW interactions. Moreover, SCAN provides improvement over the short-range interactions in liquid water, i.e., the O–H covalent bond and the overly strong directional H-bonds as described by PBE. As will be discussed with more details in the next paragraph and Sec. III E, we find that SCAN improves the length of O–H covalent bond and electronic structure of liquid water. Consequently, the weakened H-bonds lead to more broken H-bonds and more disordered local tetrahedral structure in liquid water, which contribute to the additional softening of the shell structure observed in gOO(r). Although the gOO(r) obtained from SCAN and the experimental result are almost identical beyond 3.5 Å, we notice that within the first coordination shell SCAN still gives slightly over-structured liquid water as compared to the experimental result. The still over-structured gOO(r) described by SCAN in the first coordination shell is possibly due to the lack of NQEs and the existence of SIE in our current modeling8,28,43 and can be taken into account with more advanced simulation methods such as hybrid SCAN functionals and will be addressed in future work.

We now turn our attention to the O–H radial distribution function, gOH(r), as shown in Fig. 2, where a similar trend of systematic softening of structure can be observed with the use of PBE+vdW and SCAN functionals as compared to the PBE results. Again, we find that SCAN outperforms PBE+vdW for the description of the gOH(r) in several ways. First, the position of the first maximum of the gOH(r), r1max, which corresponds to the O–H covalent bond, is shortened from 0.994 Å in PBE to 0.978 Å in SCAN AIMD simulations; however, there is only a minor change in the PBE+vdW simulation, which gives the O–H bond length of 0.990 Å. This shorter length of covalent bond as described by SCAN confirms the enhancement of the covalency of water molecules. This result is also consistent with the previous work,53 where it was reported that SCAN is highly accurate for properties of gas-phase water monomers including O–H covalent bond length and H–O–H covalent bond angle. As the O and H atoms in a molecule bind more strongly, it is more difficult for an H atom to form an H-bond with its neighboring O atoms, resulting in more broken H-bonds in a softer structure of liquid water.

FIG. 2.

O–H radial distribution functions gOH(r) of liquid water as obtained from AIMD simulations in the canonical ensemble employing PBE (red solid line), PBE+vdW (black dashed line), and SCAN exchange-correlation functionals (blue dashed line). The experimental data are from joint X-ray/neutron experiments.10 

FIG. 2.

O–H radial distribution functions gOH(r) of liquid water as obtained from AIMD simulations in the canonical ensemble employing PBE (red solid line), PBE+vdW (black dashed line), and SCAN exchange-correlation functionals (blue dashed line). The experimental data are from joint X-ray/neutron experiments.10 

Close modal

The second maximum of the gOH(r) helps to understand better the H-bond network in liquid water, as it corresponds to the intermolecular O–H distances. Table II shows that both PBE+vdW and SCAN yield intensities of the second maximum g2max of the gOH(r) closer to the experimental value as compared to PBE. More interestingly, we observe that SCAN outperforms PBE+vdW in correcting both g2max and g2min (the second minimum) of the gOH(r) because the two functionals disrupt the H-bond network of liquid water in slightly different ways. On the one hand, in the PBE+vdW simulation, the interstitial region between the first and second solvation shells is populated due to the attractive forces among water molecules provided by vdW interactions, which further disrupts the first shell; however, strengths of directional H-bonds in liquid water are barely changed. On the other hand, SCAN not only captures the intermediate-range vdW interactions but also corrects the unphysical overly strong directional H-bonding strength prescribed in PBE; both corrections by SCAN contribute to the improved description of the H-bond network in liquid water. We notice that both PBE+vdW and SCAN functionals increase the position of the second maximum, r2max, from 1.72 Å in PBE to 1.75 Å in PBE+vdW, and 1.79 Å in SCAN, while the experimental value is 1.73 Å. Finally, for the third maximum of the gOH(r), we observe that the intensity of g3max from PBE+vdW only reduces by 0.01 with a 0.02 Å shift of the position r3max outwards as compared to PBE, while SCAN reduces g3max by 0.08 and a 0.04 Å shift of r3max outwards as compared to PBE. As a result, the difference between the intensities of the third maximum and the second minimum Δg2,3max=g3maxg2max improves from −0.20 in PBE to +0.18 in SCAN, leading to a significant improvement toward the experimental value of +0.4. However, SCAN still gives an over-structured gOH(r), and the discrepancy is larger than that of gOO(r). We rationalize this discrepancy based on the fact that NQEs affect H atoms more severely than O atoms;8,9 an explicit treatment of NQEs would be necessary to further examine the gOH(r).

TABLE II.

Tabulated summaries of the structural properties of liquid water from the O–H radial distribution functions gOH(r) obtained from AIMD simulations in the canonical ensemble and experimental data. From left to right, positions (in Å) and intensities of the first maximum (r1max and g1max), second maximum (r2max and g2max), second minimum (r2min and g2min), and third maximum (r3max and g3max) of the corresponding gOH(r).

gOH(r)r1maxg1maxr2maxg2maxr2ming2minr3maxg3max
PBE 0.994 32.1 1.72 1.88 2.49 0.06 3.21 1.68 
PBE+vdW 0.990 32.6 1.75 1.66 2.45 0.09 3.23 1.67 
SCAN 0.978 36.4 1.79 1.42 2.40 0.18 3.25 1.60 
Expt.17  … … 1.73 1.01 2.39 0.28 3.30 1.41 
gOH(r)r1maxg1maxr2maxg2maxr2ming2minr3maxg3max
PBE 0.994 32.1 1.72 1.88 2.49 0.06 3.21 1.68 
PBE+vdW 0.990 32.6 1.75 1.66 2.45 0.09 3.23 1.67 
SCAN 0.978 36.4 1.79 1.42 2.40 0.18 3.25 1.60 
Expt.17  … … 1.73 1.01 2.39 0.28 3.30 1.41 

After discussions of the gOO(r) and gOH(r), we discuss the H–H radial distribution function, gHH(r), as illustrated in Fig. 3. As listed in Table III, for the first maximum of the gHH(r), we observe that while PBE and PBE+vdW give nearly identical results, SCAN yields a higher first maximum g1max and a decreased position r1max. As the first peak corresponds to the distance between the two H atoms within a water molecule, these changes support that SCAN alters the intramolecular short-range interactions and enhances the covalency of water molecules. Meanwhile, the improvements over the second maximum and the second minimum by PBE+vdW and SCAN compared to PBE are similar to the previously discussed cases of the first maximum of the gOO(r) and the second maximum of the gOH(r). These improvements indicate that SCAN corrects the water structure more significantly than PBE+vdW does because the former functional not only captures the intermediate-range vdW interactions but also reduces the directional H-bond strength.

FIG. 3.

H–H radial distribution functions gHH(r) of liquid water as obtained from AIMD simulations in the canonical ensemble employing PBE (red solid line), PBE+vdW (black dashed line), and SCAN XC functionals (blue dashed line). The experimental data are from joint X-ray/neutron experiments.10 

FIG. 3.

H–H radial distribution functions gHH(r) of liquid water as obtained from AIMD simulations in the canonical ensemble employing PBE (red solid line), PBE+vdW (black dashed line), and SCAN XC functionals (blue dashed line). The experimental data are from joint X-ray/neutron experiments.10 

Close modal
TABLE III.

Tabulated summaries of the structural properties of liquid water from the H–H radial distribution functions gHH(r) obtained from AIMD simulations in the canonical ensemble and experimental data. From left to right, positions (in Å) and intensities of the first maximum (r1max and g1max), second maximum (r2max and g2max), second minimum (r2min and g2min), and third maximum (r3max and g3max) of the corresponding gHH(r).

gHH(r)r1maxg1maxr1ming1minr2maxg2maxr2ming2min
PBE 1.59 2.98 1.79 0.08 2.25 1.76 2.90 0.44 
PBE+vdW 1.58 3.05 1.79 0.07 2.29 1.63 2.87 0.55 
SCAN 1.56 3.31 1.79 0.05 2.31 1.48 2.93 0.69 
Expt.10  … … … … 2.38 1.29 2.98 0.82 
gHH(r)r1maxg1maxr1ming1minr2maxg2maxr2ming2min
PBE 1.59 2.98 1.79 0.08 2.25 1.76 2.90 0.44 
PBE+vdW 1.58 3.05 1.79 0.07 2.29 1.63 2.87 0.55 
SCAN 1.56 3.31 1.79 0.05 2.31 1.48 2.93 0.69 
Expt.10  … … … … 2.38 1.29 2.98 0.82 

In order to further analyze the local arrangement of water molecules in liquid phase, we have computed the distribution of triplet oxygen-oxygen-oxygen angles, POOO(θ), within the first coordination shell (see Fig. 4) and the tetrahedral order parameter q (see Table IV). To compute the POOO(θ), three O atoms were considered as part of a given triplet if two of the O atoms were within a cutoff distance from the third one so that the average O–O coordination number reaches 4.0 at the cutoff distance.10,28 The cutoff distances were selected as 3.26, 3.24, and 3.20 Å from PBE, PBE+vdW, and SCAN AIMD simulations, respectively. The tetrahedral order parameter q is determined by

q=138i=13j=i+14(cos(θij)+13)2
(1)

and varies between 0 and 1, where 0 represents an ideal gas, while 1 represents a perfect tetrahedral structure. The computed results of POOO(θ) and q from our AIMD simulations are also compared with results of EPSR based on the joint X-ray/neutron data.10 

FIG. 4.

O–O–O triplet angular distribution functions POOO(θ) of liquid water obtained from AIMD simulations and EPSR based on joint X-ray/neutron scattering data.10 The triplet angular distribution functions shown here were normalized to 0πdθPOOOθsinθ=1.

FIG. 4.

O–O–O triplet angular distribution functions POOO(θ) of liquid water obtained from AIMD simulations and EPSR based on joint X-ray/neutron scattering data.10 The triplet angular distribution functions shown here were normalized to 0πdθPOOOθsinθ=1.

Close modal
TABLE IV.

Tabulated summary of the properties of liquid water obtained via the AIMD simulations performed in this work and various experiments. From left to right, the average number of H-bonds (nHB) per water molecule; the average tetrahedrality parameter (q) per water molecule; the average size of the rings (Sring) threaded through one water molecule; the average number of the rings (Nring) threaded through one water molecule; the average diffusion coefficient (D in Å2/ps) per AIMD simulation; the average dipole moment (μ in D) per water molecule; the energy difference between the peaks of 2a1 and 1b1 valence molecular orbits (ΔE=E1b1E2a1 in eV) in DOS; the electronic band gaps (Eg) from each AIMD simulations (in eV). For reference, the available corresponding experimental data are provided in the last row.

MethodnHBqSringNringDμΔEEg
PBE 3.85 ± 0.52 0.895 ± 0.459 6.64 ± 1.69 9.63 ± 2.84 0.015 ± 0.012 3.19 ± 0.29 17.63 4.51 ± 0.10 
PBE+vdW 3.77 ± 0.64 0.839 ± 0.478 7.01 ± 2.04 9.81 ± 3.39 0.041 ± 0.023 3.13 ± 0.30 15.63 4.32 ± 0.14 
SCAN 3.60 ± 0.77 0.709 ± 0.473 7.75 ± 2.58 10.01 ± 4.22 0.129 ± 0.046 2.94 ± 0.28 18.98 4.90 ± 0.08 
Expt. 3.5883  0.57610  … … 0.18684  2.976  19.7415  8.7 ± 0.672  
MethodnHBqSringNringDμΔEEg
PBE 3.85 ± 0.52 0.895 ± 0.459 6.64 ± 1.69 9.63 ± 2.84 0.015 ± 0.012 3.19 ± 0.29 17.63 4.51 ± 0.10 
PBE+vdW 3.77 ± 0.64 0.839 ± 0.478 7.01 ± 2.04 9.81 ± 3.39 0.041 ± 0.023 3.13 ± 0.30 15.63 4.32 ± 0.14 
SCAN 3.60 ± 0.77 0.709 ± 0.473 7.75 ± 2.58 10.01 ± 4.22 0.129 ± 0.046 2.94 ± 0.28 18.98 4.90 ± 0.08 
Expt. 3.5883  0.57610  … … 0.18684  2.976  19.7415  8.7 ± 0.672  

It can be inferred from Fig. 4 that the EPSR-experimental POOO(θ) exhibits a broad distribution with its maximum at around 100.5°, suggesting that the local tetrahedral network in liquid water is considerably more disordered than that in ice. At the PBE level of theory, the overall distribution of POOO(θ) is significantly narrower compared to the EPSR-experimental result. Similarly, the tetrahedral order parameter q in our PBE simulation has a value of 0.895, which is much higher than the EPSR-experimental value of 0.576.10 Both POOO(θ) and q computed from the PBE-based AIMD simulation corroborates with the aforementioned fact that liquid water modeled by the PBE functional tends to be over-structured with overly strong H-bonding. In Ref. 28, the simulation of liquid water by PBE at 300 K gives a POOO(θ) distribution with the peak at ∼109°, which is very close to the perfect tetrahedral angle of 109.5°. In our simulation, the peak position is shifted to 104° at an elevated temperature of 330 K. We consider that the ∼5° difference of the peak position originates from the additional 30 K in our simulation, as when temperature increases, the topological structure close to the perfect tetrahedral network is less favored. Noticeably, two simulations utilizing the PBE0+vdW level of theory in Ref. 28 at 300 and 330 K gave rise to peak positions of POOO(θ) at around 106° and 102°, respectively; the ∼4° difference caused by elevated temperature is in consistency with our simulation results.

The degree of tetrahedral order in liquid water is slightly reduced when vdW interactions are included in the simulation. It is observed that the distribution of POOO(θ) as calculated from the PBE+vdW trajectory is widened with a noticeable reduction of peak intensity when compared to that obtained from the PBE trajectory, while the tetrahedral order parameter q reduces from 0.895 (PBE) to 0.839 (PBE+vdW). SCAN further improves the tetrahedral order as evidenced by the fact that the peak intensity is lowered, the distribution of POOO(θ) is widened, and the q value is reduced to 0.709. Also, we find that the peak position of POOO(θ) keeps the same value of 104° by both PBE and PBE+vdW but reduces to 101° by SCAN, which is very close to the experimental value of 100.5°. This is rationalized by the fact that the peak position is correlated with the most perfect local tetrahedral structure that is formed by directional H-bonds, which are weakened by SCAN as previously mentioned. The results shown here suggest that SCAN also accurately captures the local tetrahedral structure of liquid water.

In an ice system, generally, each water molecule stably donates two and accepts two H-bonds, resulting in four H-bonds (intact H-bonds) on average. The fluidity of liquid water comes from the fact that H-bonds, which form the tetrahedral-like network, are continuously breaking and forming on a time scale of ps. By comparing the statistics of H-bonds in different AIMD trajectories, we obtain information on the H-bond strength and the fluidity of the system. The notion of H-bond itself, however, is not uniquely defined; although any quantitative description of the statistics for intact H-bonds yields certain degree of ambiguity, a general agreement among many proposed definitions for intact H-bonds still exist.62 We follow the H-bond definition proposed by Luzar and Chandler63 and define an intact H-bond when O–O distance ROO < 3.5 Å and β < 30°, where β ≡ ∠OA…OD − HD denotes the bond angle formed by an oxygen atom OA that accepts a H-bond (HD) from a nearby oxygen atom OD.

Using this definition of H-bonds, we compute the average numbers of H-bonds per water molecule in different AIMD trajectories (see Table IV), which are 3.85, 3.77, and 3.60 for PBE, PBE+vdW, and SCAN functionals, respectively. We observe that while both PBE+vdW and SCAN functionals increase the proportion of broken H-bonds and the degree of disorder in the local tetrahedral network of liquid water, the former functional has a smaller effect than the latter one. As discussed in Secs. III A and III B, vdW interactions disturb the H-bond network by populating the interstitial region and disrupting the local tetrahedral structure. SCAN also captures the intermediate-range vdW interactions but additionally improves the covalency of water molecules and electronic structures (see Sec. III E). Moreover, SCAN weakens the overly strong directional H-bonding strength described by PBE, which vdW interactions fail to correct. The improved description of directional H-bonding by SCAN contributes to the further reduction of the average number of H-bonds when compared to PBE+vdW.

To further illustrate the H-bond statistics in liquid water, we compute the time-averaged distribution of H-bonds in terms of their numbers (see Fig. 5). In the PBE trajectory, the proportion of water molecules that have four H-bonds is 79.5%; this number reduces to 69.7% and 56.4% in the PBE+vdW and SCAN trajectories, respectively. Therefore, we find that both PBE+vdW and SCAN functionals reduce the proportion of water molecules with four H-bonds and increase the proportions of water molecules with 1, 2, 3, and 5 H-bonds. Once again, SCAN functionals affect the H-bond statistics more than PBE+vdW does. This result is consistent with our previous analyses that the PBE+vdW and SCAN functionals influence the H-bond network of liquid water differently.

FIG. 5.

Distribution of number of H-bonds in liquid water obtained from AIMD simulations employing PBE, PBE+vdW, and SCAN XC functionals.

FIG. 5.

Distribution of number of H-bonds in liquid water obtained from AIMD simulations employing PBE, PBE+vdW, and SCAN XC functionals.

Close modal

The improvement over the average number of H-bonds from our AIMD simulations also implies an enhancement of the fluidity of liquid water since the breakage and formation of H-bonds through thermal fluctuations affect the diffusivity of water molecules: as there are less H-bonds with the H-bonding strength reduced, the diffusion coefficient D increases (see Table IV). These D values are extracted from the slope of the molecular mean square displacement (MSD) versus the simulation time. However, these D values can be further improved based on the following considerations. First, small periodic cells have a non-negligible finite-size effect on D,26,64 in which diffusion coefficients converge slowly with the increased sizes of the simulation cell. Second, a much longer trajectory or multiple runs can reduce the uncertainty of D.65 

We can gain more understanding of the liquid water structure by analyzing the ring statistics. The distribution of closed rings was not only used to study topological networks in silicate network structures66,67 but also in liquid water systems,68–70 where rings are threaded by O atoms that are connected with H-bonds. In this work, rings were defined adopting the shortest-path definition from Ref. 68, in which only the shortest circuit with two continuous H-bonds from one water molecule is recorded. The analyses of rings include average size of rings, distribution of ring sizes, average number of rings, and distribution of ring numbers (see Fig. 6 and Table IV). Note that a ring size is defined as the number of O atoms in a closed ring.

FIG. 6.

(a) Probability distribution of rings as a function of ring sizes obtained from AIMD simulations by using PBE, PBE+vdW, and SCAN functionals. The integration of the area for each curve is the average number of rings threading one water molecule. (b) Distribution of number of rings threaded through one water molecule obtained from AIMD simulations by employing PBE, PBE+vdW, and SCAN functionals.

FIG. 6.

(a) Probability distribution of rings as a function of ring sizes obtained from AIMD simulations by using PBE, PBE+vdW, and SCAN functionals. The integration of the area for each curve is the average number of rings threading one water molecule. (b) Distribution of number of rings threaded through one water molecule obtained from AIMD simulations by employing PBE, PBE+vdW, and SCAN functionals.

Close modal

In an ice Ih system, there are in general 12 6-membered rings threading through each water molecule. In other words, a water molecule is involved in 12 rings, with each ring threading 6 water molecules. However, in a liquid water system, one would expect different ring sizes and ring numbers for different sizes, as the broken H-bonds and disordered H-bonding network give rise to rings with irregular shapes and sizes other than 6. We define small rings as those threading 4 molecules and less, medium rings as those threading 5–8 molecules, and large rings as those threading more than 8 molecules, by loosely choosing the cross points in the probability distributions of ring sizes by different XC functionals in Fig. 6(a). The existence of small and medium rings is mainly caused by the disruption of local tetrahedrality in liquid water, while the large rings are mainly formed by the combination of two or more small and/or medium rings when shared H-bond(s) between adjacent rings break by thermal fluctuations. As shown in Fig. 6(a), both PBE+vdW and SCAN reduce the numbers of small and medium rings and increase the numbers of large rings, with SCAN outperforming PBE+vdW. This result is consistent with the aforementioned structural analyses and further confirms our observation that SCAN describes a more disordered H-bonding network in liquid water.

We further examine the statistics for the sizes of rings as shown in Fig. 6(a). Both PBE+vdW and SCAN functionals reduce the portion of medium rings and increase the portions of small and large rings as compared to the PBE results. Particularly, the proportion of 6-membered rings, which is the dominating size of rings found in liquid water in all the three functionals tested, has been reduced by PBE+vdW as compared to the PBE result; the proportion of 6-membered ring is further reduced by SCAN. Notably, SCAN significantly increases the portion of rings threading more than 8 water molecules. For instance, the average numbers of 10-membered rings obtained by PBE, PBE+vdW, and SCAN functionals are 0.12, 0.21, and 0.60, respectively. Accordingly, the average ring sizes for PBE, PBE+vdW, and SCAN functionals are 6.64, 7.01, and 7.75, respectively. The reduction of medium rings and increase of small and large rings in both PBE+vdW and SCAN functionals are associated with the fact that both functionals disrupt the local tetrahedrality, soften the coordination shells, and increase the portion of broken H-bonds (see Secs. III B and III C). However, the influence on the distribution of ring sizes by SCAN is more significant than PBE+vdW because SCAN not only captures the intermediate-range vdW interactions but also corrects the overly strong directional H-bond strength in the PBE simulation.

Figure 6(b) shows the distributions of number of rings threading one water molecule, in which we find that the distribution is broadened with more small and large numbers of rings by PBE+vdW as compared to the one by PBE and is even more broadened when SCAN is considered. A broader distribution reflects an H-bond network that is softer and more distorted, which is consistent with the fact that SCAN has a larger effect on disrupting the H-bond network.

In this section, we discuss electronic properties of liquid water which govern the strength of H-bonds and influence solvation behaviors. We first analyze the electronic density of states (DOS) and band gaps of liquid water obtained by AIMD simulations with different XC functionals and then discuss molecular dipole moments. As we will see, the analyses on electronic properties show that the SCAN functional improves the electronic structure of molecules in liquid water, which leads to the correction of the overly strong H-bonds by PBE, while PBE+vdW does not provide such improvement.

The improvement in the electronic properties of liquid water achieved by the SCAN functional can be seen through the comparison of the computed electronic DOS with the full valence band photoemission spectroscopy15 (see Fig. 7). The DOS of liquid water in each AIMD simulation was computed through a Gaussian broadening of 0.5 eV71 for eigenstates of the system and averaged over 10 snapshots through the trajectory. The four peaks of the DOS are assigned to the 2a1, 1b2, 3a1, and 1b1 valence molecular orbitals based on the spatial symmetries of a water molecule, with the lone electron pairs closely connected to the 2a1 and 1b1 orbitals and the bonding electron pairs related to the 1b2 and 3a1 orbitals. The simulated DOS and the photoemission spectra are aligned following Ref. 72 in such a way that the onset position, which was taken as the x intercept of the linear fitting for the tail of the 1b1 distribution, is set to 0 eV.

FIG. 7.

Density of states (DOS) of water obtained from AIMD simulations employing PBE, PBE+vdW, and SCAN XC functionals. The experimental data of DOS from photoemission spectra15 are also shown.

FIG. 7.

Density of states (DOS) of water obtained from AIMD simulations employing PBE, PBE+vdW, and SCAN XC functionals. The experimental data of DOS from photoemission spectra15 are also shown.

Close modal

From Fig. 7, we notice that PBE+vdW yields an almost identical DOS compared to PBE, which is in accordance with the previous report that the inclusion of vdW interactions do not lead to any significant effect on either DOS or the band gap,73 whereas the electronic structure is improved by SCAN. We find that in SCAN, the 2a1 peak shifts 1.38 eV to a smaller energy compared to PBE, which is consistent with previous study54 and agrees better with the experimental value. As a result, the energy differences ΔE between the 2a1 and 1b1 peaks predicted by PBE, PBE+vdW, and SCAN are 17.63, 17.63, and 18.98 eV, respectively, as compared to the experimental value of 19.74 eV.15 This suggests that compared to the value obtained by PBE, SCAN significantly improves the energy difference between the two states and better describes the lone pair electrons. We have also computed the electronic band gaps as obtained through an average of 10 configurations for each functional. The band gaps are 4.51 ± 0.10, 4.32 ± 0.14, and 4.90 ± 0.08 eV in AIMD simulations utilizing PBE, PBE+vdW, and SCAN functionals, respectively. The band gaps obtained with the GGA level of theory are consistent with previous reports,6 while the band gap calculated with SCAN is corrected toward, although still much lower than, the experimental value of 8.7 eV.74 

Next, we examine the distribution of the molecular dipole moments in liquid water. For reference, the dipole moment of a single water molecule in the gas phase is accurately known to be 1.855 D from the electroscopic experiment,75 and several DFT functionals are capable of reproducing this experimental value to within 3%.6 However, the molecular dipole moment of liquid water extracted from X-ray scattering form factors yields a large uncertainty (2.9 ± 0.6 D).76 In the meantime, for the calculation of dipole moments from AIMD simulations of condensed-phased systems, the partitioning of the electron density would be required and different partitioning schemes would lead to widely different values of dipole moments.77 In this regard, maximally localized Wannier functions (MLWFs), which are obtained through a unitary transformation of occupied Kohn-Sham orbitals, have been shown to be very useful in computing dipole moments in condensed-phase environments.6,78–80 It has been shown that in bulk water the amount of overlap between the charge distribution of Wannier functions associated with different molecules is less than 1% of the norm,81 which allows for a nearly unique definition of the charges belonging to a given water molecule and eliminates to a large extent the partitioning issues when computing dipole moments in liquid water.

In our simulation, four doubly occupied valence MLWFs represent centers of eight valence electrons associated with each water molecule. Two of the MLWFs centered on the O–H covalent bonds are defined as bonding pair electrons, while the other two being approximately centered on the remaining tetrahedral directions are labeled as lone pair electrons. The dipole moment (μ) for a water molecule can be computed as follows:

μ=RH1+RH2+6RO2i=14RWi,
(2)

in which RH1, RH2, and RO are the Cartesian coordinates of the two H atoms and one O atom of a water molecule, respectively, and RWi are the coordinates of the four corresponding MLWF centers. With this formula, we have computed the distribution of dipole moments in liquid water from AIMD trajectories with PBE, PBE+vdW, and SCAN XC functionals, and the results are shown in Fig. 8(a). It is evident that the peak positions of the dipole distributions in Fig. 8(a) obtained from PBE+vdW and SCAN functionals are shifted toward smaller values as compared to that from PBE. In particular, the shift in the dipole maximum is much larger with SCAN than PBE+vdW. The average dipole moments are 3.19, 3.13, and 2.94 D for simulations with PBE, PBE+vdW, and SCAN functionals, respectively. Additionally, we note that the height of the peak for the distribution of dipole moments from the PBE+vdW simulation is the lowest among the three AIMD simulations, which is associated with the widening of distribution of lone pair O-MLWF distances, as will be explained next.

FIG. 8.

(a) Probability density plots of the distributions of molecular dipole moment magnitudes (μ in D) in liquid water obtained from the AIMD simulations employing PBE, PBE+vdW, and SCAN XC functionals. (b) Probability density plots of the distributions of distance between the oxygen atoms in a given water molecule and the centers of the four nearest maximally localized Wannier functions (RO-MLWF in Å).

FIG. 8.

(a) Probability density plots of the distributions of molecular dipole moment magnitudes (μ in D) in liquid water obtained from the AIMD simulations employing PBE, PBE+vdW, and SCAN XC functionals. (b) Probability density plots of the distributions of distance between the oxygen atoms in a given water molecule and the centers of the four nearest maximally localized Wannier functions (RO-MLWF in Å).

Close modal

To further understand the variation of dipole moments in liquid water given by different XC functionals, we also illustrate in Fig. 8(b) the electronic properties of liquid water by computing the distribution of distance between an O atom of water molecule and the centers of its four nearest MLWFs (ROMLWF=RORW). In general, Fig. 8(b) shows that the two bonding pair electrons have a distribution of distances longer than ROMLWF= 0.45 Å, while the two lone pair electrons are well separated from the bonding pair electrons and have a distribution of distances shorter than ROMLWF=0.40 Å. We also list the O-MLWF distances for a single water molecule for reference. Specifically, the PBE, PBE+vdW, and SCAN functionals, respectively, yield 0.523, 0.523, and 0.518 Å as the distances between bonding pair electrons and the oxygen. Moreover, the distances between lone pair electrons and the oxygen are 0.289, 0.289, and 0.290 Å as obtained from the PBE, PBE+vdW, and SCAN functionals, respectively. In the liquid phase, a water molecule accepts an H-bond through the electrostatic interaction between one of its lone pair electrons and an H atom from one of its neighboring water molecules. If both the H atom and the lone pair electron are closer to the O atom, then the accepting side of the H-bond, i.e., the O atom, is in a less negative electric environment. Meanwhile, the donating side of the H-bond, i.e., the H atom is in a less positive electric environment.

Notably, the distributions of bonding pair electrons by the three XC functionals only exhibit small differences, while the PBE+vdW and SCAN functionals change more significantly the distribution of the lone pair electrons. Specifically, the PBE+vdW functional decreases the peak position of the lone pair electrons from 0.338 Å, as obtained by the PBE functional, to 0.334 Å and lowers the peak intensity and broadens the distribution. As discussed before, smaller distances between O atoms and their corresponding lone pair electrons lead to less negative environments for the O atoms, which in turn accept less H-bonds. This is consistent with the fact that single molecules in gas phase have no H-bonds and relatively smaller electric dipole moments as compared to those in condensed phases. Therefore, these changes of dipole moments and centers of MLWFs can be rationalized by the inclusion of vdW interactions which populate the interstitial region between the first and second solvation shells and reduce the average number of H-bonds in liquid water.

The SCAN functional yields a distribution of lone pair electrons with almost unchanged broadness, but the peak position shifts 0.012 Å inwards to 0.326 Å as compared to that from PBE. The appreciable net reduction of the ROMLWF distance and covalent bond length (see Sec. III A) by SCAN contribute as the main sources to the decrease of the average dipole moment for liquid water. Importantly, SCAN changes covalent bond strengths in liquid water, which results in the reduced polarizability and the weakened H-bond strength; this effect is absent in the modeling of liquid water using PBE+vdW. In conclusion, this unique change of electrostatic interactions by SCAN functionals provides a more accurate description of H-bond network in liquid water and is consistent with previous analyses on the structural and dynamical properties of liquid water by the SCAN functional that differ from the simulation with the PBE+vdW functional.

In this section, we compare results obtained from simulations using different ensembles (NVT and NpT) with PBE and SCAN functionals for ambient liquid water. We have taken the NpT ensemble results from the previous work by Chen et al.54 and compared with our current NVT simulations. In principle, these two ensembles should yield the same structure of liquid water if the underlying potential is sufficiently accurate to predict the correct equilibrium volume (V) at a given external pressure (p) and the simulation time is long enough for sampling the phase space. However, the equilibrium volumes/densities at 1 bar predicted with both PBE and SCAN functionals in the NpT ensemble deviate from the experimental value (∼1 g/mL). As a result, we find in Fig. 9 that the NVT ensemble yields a slightly more structured liquid water [higher first and second peaks in gOO(r)] when compared to that from the NpT ensemble.

FIG. 9.

(a) O–O radial distribution functions gOO(r) of liquid water obtained from PBE-based AIMD simulations in the NVT (red solid line) and the NpT (purple dashed line) ensembles. (b) O–O radial distribution functions gOO(r) of liquid water obtained from SCAN-based AIMD simulations in the NVT (blue solid line) and the NpT (green dashed line) ensembles.

FIG. 9.

(a) O–O radial distribution functions gOO(r) of liquid water obtained from PBE-based AIMD simulations in the NVT (red solid line) and the NpT (purple dashed line) ensembles. (b) O–O radial distribution functions gOO(r) of liquid water obtained from SCAN-based AIMD simulations in the NVT (blue solid line) and the NpT (green dashed line) ensembles.

Close modal

According to the recent study by Chen et al., the average densities of liquid water obtained from AIMD simulations with PBE and SCAN functionals in the NpT ensemble were 0.850 and 1.050 g/mL, respectively;54 another work on the simulation of liquid water also shows that utilizing SCAN overcorrects the density of water compared to experiment.82 On the other hand, the density of liquid water is fixed as 1.0 g/mL in our NVT simulation. Therefore, the smaller density as predicted by the PBE AIMD trajectory in the NpT ensemble implies that water molecules are more separated than those in the SCAN AIMD trajectory. Indeed, we observe that with the PBE functional, the positions of maximums and minimums of gOO(r) in the NpT ensemble shift outwards compared to those in the NVT ensemble. Specifically, the position of the first maximum r1max increases from 2.71 Å in the NVT ensemble to 2.73 Å in the NpT ensemble. Accordingly, the averaged pressures of liquid water, calculated from 50 snapshots from the NVT ensemble with the kinetic energy cutoff set as 150 Ry, are 0.9 ± 3.7, −7.0 ± 4.1, and −6.5 ± 3.4 kbar in the trajectories adopting the PBE, PBE+vdW, and SCAN functionals, respectively. By softening the structure of liquid water, SCAN corrects the large density error in PBE in the NpT ensemble. However, the use of the experimental density in the NVT ensemble (1.0 g/mL) than the one predicted by NpT (1.05 g/mL) mitigates the correction effect by SCAN. As a result, water molecules in the NVT ensemble are more apart with their shell structures become slightly over-structured; we find that the positions of maximums and minimums in gOO(r) as obtained from the SCAN trajectory in the NpT ensemble are squeezed slightly inwards compared to those from the NVT ensemble, which is consistent with the slightly higher density obtained in the NpT ensemble. To be specific, the position of the first maximum r1max decreases from 2.75 Å in the NVT to 2.74 Å in the NpT. Consistently, the self-diffusion coefficient of 0.129 Å2/ps in liquid water obtained from the SCAN trajectory in the NVT ensemble is smaller than the value of 0.190 Å2/ps from the NpT ensemble.

We have performed AIMD simulations of liquid water at ambient conditions in the canonical ensemble with the focus on comparing the results as obtained from the newly developed meta-GGA SCAN XC functional to those calculated with the PBE and PBE+vdW functionals. Importantly, SCAN captures the intermediate-range vdW interaction that is crucial in improving the structure of liquid water. Besides, SCAN provides a better description of electronic structures that lead to a more delicate balance between covalent bonds and H-bonds in liquid water, largely curing the overly strong H-bonds and softening the rigid solvation shell structure as prescribed by the PBE functional. We find that SCAN yields more accurate structural properties like radial distribution functions and O–O–O triplet angular distribution functions, improves the H-bond network as can be shown in the H-bond statistics and ring statistics, and predicts electronic properties including DOS, band gap, and dipole moment that are closer to experimental data. Furthermore, we discussed properties of liquid water by SCAN in both NVT and NpT ensembles. In general, SCAN performs better than both PBE and PBE+vdW functionals for describing liquid water. However, there are still improvements to be made regarding the structure of liquid water. For instance, the current modeling by SCAN still prescribes a slightly over-structured solvation shell structures. The correction for SIE by the hybrid functional was reported to further mitigate the over-structuring,28,43 but the performance of combining SCAN with the hybrid functional is still unknown and can be investigated in future studies. Additionally, the explicit treatment of NQEs through the Feynman discretized path-integral approach8,9,72 can be used to correct the discrepancies regarding the widths and intensities of radial distribution functions, especially of gOH(r). Nevertheless, the SCAN-based AIMD simulation in the canonical ensemble already describes a delicate balance among covalent bonds, H-bonds, and vdW interactions that lead to reasonably accurate structures and dynamics of liquid water. Therefore, the SCAN XC functional provides a promising ab initio model of liquid water and should be considered for future work on modeling of aqueous solutions.

This work was supported by National Science Foundation through Award No. DMR-1552287 (X.W.). This research used resources of the National Energy Research Scientific Computing Center, which is supported by the U.S. Department of Energy (DOE), Office of Science under Contract No. DEAC0205CH11231. The work of Z.S. was supported as part of the Center for the Computational Design of Functional Layered Materials, an Energy Frontier Research Center funded by the U.S. DOE, Office of Science, Basic Energy Sciences under Award No. DESC0012575. The work of L.Z. was partially supported by the American Chemical Society Petroleum Research Fund (ACS PRF) under Grant No. 53482-DNI6.

2.
F.
Labat
,
C.
Pouchan
,
C.
Adamo
, and
G. E.
Scuseria
,
J. Comput. Chem.
32
,
2177
(
2011
).
3.
B.
Santra
,
J.
Klimeš
,
A.
Tkatchenko
,
D.
Alfè
,
B.
Slater
,
A.
Michaelides
,
R.
Car
, and
M.
Scheffler
,
J. Chem. Phys.
139
,
154702
(
2013
).
4.
O.
Kambara
,
K.
Takahashi
,
M.
Hayashi
, and
J.-L.
Kuo
,
Phys. Chem. Chem. Phys.
14
,
11484
(
2012
).
5.
J.
Wang
,
G.
Romn-Prez
,
J. M.
Soler
,
E.
Artacho
, and
M. V.
Fernndez-Serra
,
J. Chem. Phys.
134
,
024516
(
2011
).
6.
C.
Zhang
,
J.
Wu
,
G.
Galli
, and
F.
Gygi
,
J. Chem. Theory Comput.
7
,
3054
(
2011
).
7.
B.
Santra
,
J.
Klimeš
,
D.
Alfè
,
A.
Tkatchenko
,
B.
Slater
,
A.
Michaelides
,
R.
Car
, and
M.
Scheffler
,
Phys. Rev. Lett.
107
,
185701
(
2011
).
8.
J. A.
Morrone
and
R.
Car
,
Phys. Rev. Lett.
101
,
017801
(
2008
).
9.
M.
Ceriotti
,
W.
Fang
,
P. G.
Kusalik
,
R. H.
McKenzie
,
A.
Michaelides
,
M. A.
Morales
, and
T. E.
Markland
,
Chem. Rev.
116
,
7529
(
2016
).
10.
A. K.
Soper
and
C. J.
Benmore
,
Phys. Rev. Lett.
101
,
065502
(
2008
).
11.
T.
Head-Gordon
and
H.
Greg
,
Chem. Rev.
102
,
2651
(
2002
).
12.
A.
Nilsson
and
L. G. M.
Pettersson
,
Chem. Phys.
389
,
1
(
2011
).
13.
C. J.
Fecko
,
J. D.
Eaves
,
J. J.
Loparo
,
A.
Tokmakoff
, and
P. L.
Geissler
,
Science
301
,
1698
(
2003
).
14.
E. T. J.
Nibbering
and
E.
Thomas
,
Chem. Rev.
104
,
1887
(
2004
).
15.
B.
Winter
,
R.
Weber
,
W.
Widdra
,
M.
Dittmar
,
M.
Faubel
, and
I. V.
Hertel
,
J. Phys. Chem. A
108
,
2625
(
2004
).
16.
B.
Winter
and
F.
Manfred
,
Chem. Rev.
106
,
1176
(
2006
).
17.
L. B.
Skinner
,
C.
Huang
,
D.
Schlesinger
,
L. G. M.
Pettersson
,
A.
Nilsson
, and
C. J.
Benmore
,
J. Chem. Phys.
138
,
074506
(
2013
).
18.
A. K.
Soper
,
ISRN Phys. Chem.
2013
,
279463
(
2013
), 67 pages.
19.
R.
Car
and
M.
Parrinello
,
Phys. Rev. Lett.
55
,
2471
(
1985
).
20.
K.
Burke
,
J. Chem. Phys.
136
,
150901
(
2012
).
21.
M.
Tuckerman
,
L.
Kari
,
S.
Michiel
, and
P.
Michele
,
J. Chem. Phys.
103
,
150
(
1995
).
22.
D.
Marx
,
Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods
(
Cambridge University Press
,
New York
,
2009
).
23.
K.
Laasonen
,
M.
Parrinello
,
R.
Car
, and
C.
Lee
,
J. Chem. Phys.
99
,
9080
(
1993
).
24.
B.
Santra
,
R. A.
DiStasio
,
F.
Martelli
, and
R.
Car
,
Mol. Phys.
113
,
2829
(
2015
).
25.
D.
Asthagiri
,
L. R.
Pratt
, and
J. D.
Kress
,
Phys. Rev. E
68
,
041505
(
2003
).
26.
J. C.
Grossman
,
E.
Schwegler
,
E. W.
Draeger
,
F.
Gygi
, and
G.
Galli
,
J. Chem. Phys.
120
,
300
(
2004
).
27.
A. P.
Gaiduk
,
F.
Gygi
, and
G.
Galli
,
J. Phys. Chem. Lett.
6
,
2902
(
2015
).
28.
R. A.
DiStasio
, Jr.
,
B.
Santra
,
Z.
Li
,
X.
Wu
, and
R.
Car
,
J. Chem. Phys.
141
,
084502
(
2014
).
29.
E.
Schwegler
,
C. G.
Jeffrey
,
F.
Gygi
, and
G.
Galli
,
J. Chem. Phys.
121
,
5400
(
2004
).
30.
K.
Laasonen
,
F.
Csajka
, and
M.
Parrinello
,
Chem. Phys. Lett.
194
,
172
(
1992
).
31.
M.
Chen
,
L.
Zheng
,
B.
Santra
,
H.-Y.
Ko
,
R. A.
DiStasio
, Jr.
,
M. L.
Klein
,
R.
Car
, and
X.
Wu
,
Nat. Chem.
10
,
413
(
2018
).
32.
C. W.
Swartz
and
X.
Wu
,
Phys. Rev. Lett.
111
,
087801
(
2013
).
33.
Y.
Yao
and
Y.
Kanai
,
J. Chem. Theory Comput.
14
,
884
(
2018
).
34.
M. J.
McGrath
,
J. I.
Siepmann
,
I. F. W.
Kuo
,
C. J.
Mundy
,
J.
Vandevondele
,
J.
Hutter
,
F.
Mohamed
, and
M.
Krack
,
J. Phys. Chem. A
110
,
640
(
2006
).
35.
M. J.
McGrath
,
J. I.
Siepmann
,
I. F. W.
Kuo
,
C. J.
Mundy
,
J.
Vandevondele
,
J.
Hutter
,
F.
Mohamed
, and
M.
Krack
,
ChemPhysChem
6
,
1894
(
2005
).
36.
M.
Del Ben
,
M.
Schönherr
,
J.
Hutter
, and
J.
Vandevondele
,
J. Phys. Chem. Lett.
4
,
3753
(
2013
).
37.
M.
Del Ben
,
J.
Hutter
, and
J.
VandeVondele
,
J. Chem. Phys.
143
,
054506
(
2015
).
38.
K.
Laasonen
,
M.
Parrinello
,
R.
Car
,
C.
Lee
, and
D.
Vanderbilt
,
Chem. Phys. Lett.
207
,
208
(
1993
).
39.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
40.
A. D.
Becke
,
Phys. Rev. A
38
,
3098
(
1988
).
41.
C.
Lee
,
W.
Yang
, and
R. G.
Parr
,
Phys. Rev. B
37
,
785
(
1988
).
42.
J.
Schmidt
,
J.
VandeVondele
,
I. F. W.
Kuo
,
D.
Sebastiani
,
J. I.
Siepmann
,
J.
Hutter
, and
C. J.
Mundy
,
J. Phys. Chem. B
113
,
11959
(
2009
).
43.
J. P.
Perdew
and
A.
Zunger
,
Phys. Rev. B
23
,
5048
(
1981
).
44.
M. J.
Gillan
,
D.
Alfe
, and
A.
Michaelides
,
J. Chem. Phys.
144
,
130901
(
2016
).
45.
R. C.
Remsing
,
J. M.
Rodgers
, and
J. D.
Weeks
,
J. Stat. Phys.
145
,
313
(
2011
).
46.
A.
Bankura
,
B.
Santra
,
R. A.
DiStasio
,
C. W.
Swartz
,
M. L.
Klein
, and
X.
Wu
,
Mol. Phys.
113
,
2842
(
2015
).
47.
S.
Grimme
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
1
,
211
(
2011
).
48.
J.
Hermann
,
R. A.
Distasio
, and
A.
Tkatchenko
,
Chem. Rev.
117
,
4714
(
2017
).
49.
M.
Dion
,
H.
Rydberg
,
E.
Schrder
,
D. C.
Langreth
, and
B. I.
Lundqvist
,
Phys. Rev. Lett.
92
,
246401
(
2004
).
50.
A.
Tkatchenko
and
M.
Scheffler
,
Phys. Rev. Lett.
102
,
073005
(
2009
).
51.
C.
Zhang
,
D.
Donadio
,
F.
Gygi
, and
G.
Galli
,
J. Chem. Theory Comput.
7
,
1443
(
2011
).
52.
J.
Sun
,
A.
Ruzsinszky
, and
J. P.
Perdew
,
Phys. Rev. Lett.
115
,
036402
(
2015
).
53.
54.
M.
Chen
 et al,
Proc. Natl. Acad. Sci. U. S. A.
114
,
10846
(
2017
).
55.
P.
Giannozzi
 et al,
J. Phys.: Condens. Matter
21
,
395502
(
2009
).
56.
P.
Giannozzi
 et al,
J. Phys.: Condens. Matter
29
,
465901
(
2017
).
57.
D.
Vanderbilt
,
Phys. Rev. B
32
,
8412
(
1985
).
58.
G. J.
Martyna
,
M. L.
Klein
, and
M.
Tuckerman
,
J. Chem. Phys.
97
,
2635
(
1992
).
59.
S.
Nosé
,
J. Chem. Phys.
81
,
511
(
1984
).
60.
W. G.
Hoover
,
Phys. Rev. A
31
,
1695
(
1985
).
61.
F.
Tassone
,
F.
Mauri
, and
R.
Car
,
Phys. Rev. B
50
,
10561
(
1994
).
62.
D.
Prada-Gracia
,
R.
Shevchuk
, and
F.
Rao
,
J. Chem. Phys.
139
,
084501
(
2013
).
63.
A.
Luzar
and
D.
Chandler
,
Phys. Rev. Lett.
76
,
928
(
1996
).
64.
T. D.
Kuhne
,
M.
Krack
, and
M.
Parrinello
,
J. Chem. Theory Comput.
5
,
235
(
2009
).
65.
W.
Dawson
and
F.
Gygi
,
J. Chem. Phys.
148
,
124501
(
2018
).
66.
67.
X.
Yuan
and
A. N.
Cormack
,
Comput. Mater. Sci.
24
,
343
(
2002
).
68.
A. C.
Belch
and
S. A.
Rice
,
J. Chem. Phys.
86
,
5676
(
1987
).
69.
M.
Matsumoto
,
A.
Baba
, and
I.
Ohmine
,
J. Chem. Phys.
127
,
134504
(
2007
).
70.
A.
Hassanali
,
F.
Giberti
,
J.
Cuny
,
T. D.
Kuhne
, and
M.
Parrinello
,
Proc. Natl. Acad. Sci. U. S. A.
110
,
13723
(
2013
).
71.
D.
Prendergast
,
J. C.
Grossman
, and
G.
Galli
,
J. Chem. Phys.
123
,
014501
(
2005
).
72.
W.
Chen
,
F.
Ambrosio
,
G.
Miceli
, and
A.
Pasquarello
,
Phys. Rev. Lett.
117
,
186401
(
2016
).
73.
F.
Ambrosio
,
G.
Miceli
, and
A.
Pasquarello
,
J. Chem. Phys.
143
,
244508
(
2015
).
74.
A.
Bernas
,
C.
Ferradini
, and
J. P.
Jay-Gerin
,
Chem. Phys.
222
,
151
(
1997
).
75.
S. L.
Shostak
,
W. L.
Ebenstein
, and
J. S.
Muenter
,
J. Chem. Phys.
94
,
5875
(
1991
).
76.
Y. S.
Badyal
,
M. L.
Saboungi
,
D. L.
Price
,
S. D.
Shastri
,
D. R.
Haeffner
, and
A. K.
Soper
,
J. Chem. Phys.
112
,
9206
(
2000
).
77.
E. R.
Batista
,
S. S.
Xantheas
, and
H.
Jnsson
,
J. Chem. Phys.
111
,
6011
(
1999
).
78.
M.
Sharma
,
R.
Resta
, and
R.
Car
,
Phys. Rev. Lett.
98
,
247401
(
2007
).
79.
N.
Marzari
,
A. A.
Mostofi
,
J. R.
Yates
,
I.
Souza
, and
D.
Vanderbilt
,
Rev. Mod. Phys.
84
,
1419
(
2012
).
80.
N.
Marzari
and
D.
Vanderbilt
,
Phys. Rev. B
56
,
12847
(
1997
).
81.
P. L.
Silvestrelli
and
M.
Parrinello
,
Phys. Rev. Lett.
82
,
3308
(
1999
).
82.
J.
Wiktor
,
F.
Ambrosio
, and
A.
Pasquarello
,
J. Chem. Phys.
147
,
216101
(
2017
).
83.
A. K.
Soper
,
F.
Bruni
, and
M. A.
Ricci
,
J. Chem. Phys.
106
,
247
(
1997
).
84.
R.
Mills
,
J. Phys. Chem.
77
,
685
(
1973
).