In this work, we report the results of a series of density functional theory (DFT) based *ab initio* molecular dynamics (AIMD) simulations of ambient liquid water using a hierarchy of exchange-correlation (XC) functionals to investigate the individual and collective effects of exact exchange (Exx), via the PBE0 hybrid functional, non-local van der Waals/dispersion (vdW) interactions, via a fully self-consistent density-dependent dispersion correction, and an approximate treatment of nuclear quantum effects, via a 30 K increase in the simulation temperature, on the microscopic structure of liquid water. Based on these AIMD simulations, we found that the collective inclusion of Exx and vdW as resulting from a large-scale AIMD simulation of (H_{2}O)_{128} significantly softens the structure of ambient liquid water and yields an oxygen-oxygen structure factor, *S*_{OO}(*Q*), and corresponding oxygen-oxygen radial distribution function, *g*_{OO}(*r*), that are now in quantitative agreement with the best available experimental data. This level of agreement between simulation and experiment demonstrated herein originates from an increase in the relative population of water molecules in the interstitial region between the first and second coordination shells, a collective reorganization in the liquid phase which is facilitated by a weakening of the hydrogen bond strength by the use of a hybrid XC functional, coupled with a relative stabilization of the resultant disordered liquid water configurations by the inclusion of non-local vdW/dispersion interactions. This increasingly more accurate description of the underlying hydrogen bond network in liquid water also yields higher-order correlation functions, such as the oxygen-oxygen-oxygen triplet angular distribution, *P*_{OOO}(θ), and therefore the degree of local tetrahedrality, as well as electrostatic properties, such as the effective molecular dipole moment, that are in much better agreement with experiment.

## I. INTRODUCTION

Water is arguably the most important molecule on Earth. Without water, proteins would not fold, salts would not dissolve, and there would certainly not be life, at least as we know it. While a single water molecule has a simple and well-known structure, consisting of an oxygen atom covalently bound to two hydrogen atoms, the unique physical and chemical properties of water in the solid and liquid phases are largely due to the subtle interplay between the underlying hydrogen bond network and entropic effects, which collectively govern the interactions between water molecules. As such, an accurate understanding of the intricate details of the microscopic structure of liquid water—the net result of this competition between hydrogen bond energetics (which favor ordered structures) and temperature/entropic effects (which favor disordered structures)—requires an accurate and balanced theoretical treatment of both the nuclear and electronic degrees of freedom, and has therefore posed a substantial challenge for modern state-of-the-art *ab initio* quantum mechanical methodologies.

A highly accurate and detailed understanding of the microscopic structure of liquid water is of great importance to a number of fields, ranging from biology/biochemistry (i.e., selective ion channel functioning, solvent-mediated protein folding, etc.) to environmental science (i.e., “green” solvents, cloud formation, etc.) to energy storage/electrochemistry (i.e., aqueous electrolytes, charge-transfer interactions, catalyst design, etc.).^{1} While there is no experimental methodology currently available to directly obtain the real-space microscopic structure of liquid water,^{2,3} *ab initio* molecular dynamics (AIMD)^{4} simulations employing Density Functional Theory (DFT)^{5} as the source of the underlying quantum mechanical potential can furnish such structural information. With the AIMD technique, the nuclear potential energy surface is generated “on the fly” from the electronic ground state^{6} without the need for empirical input, thereby allowing for a quantum mechanical treatment of not only the structure and dynamics of a given molecular system of interest, but also its electronic and dielectric properties, as well as potential chemical reactions (i.e., the breaking and forming of chemical bonds). Since the initial pioneering simulations of liquid water,^{7,8} AIMD has been applied to many complex problems in biology, chemistry, and energy research, e.g., the designing of efficient catalysts for hydrogen production^{9} and the modeling of the autoionization^{10} and electrocatalytic splitting of water,^{11} to name a few.

In general, the predictive power of DFT-based AIMD simulations depends crucially upon the accuracy of the underlying exchange-correlation (XC) functional utilized in the quantum mechanical treatment of the electronic degrees of freedom. In current AIMD simulations of liquid water, the most widely used XC functionals involve the generalized gradient approximation (GGA), and it is now evident that the use of GGA-DFT has severe limitations when applied to liquid water^{12–41} as well as the crystalline phases of ice.^{42–49} For one, GGA-DFT suffers from the well-known self-interaction error^{50} which leads to excessive proton delocalization in liquid water, yielding red shifts in the predicted O–H stretching frequencies as compared to experiment.^{29,35} Second, GGA-DFT ignores non-local electron correlation effects that are responsible for van der Waals (vdW) or dispersion interactions, which are thought to be crucial, for instance, in correctly obtaining the equilibrium density of liquid water.^{34,39,40} Beyond the choice of the XC functional, most AIMD simulations of liquid water performed to date have adopted classical mechanics for the nuclear equations of motion and therefore completely neglect nuclear quantum effects (NQE), and this approximation has also been deemed insufficient for a highly accurate quantitative description of the microscopic structure of liquid water. In the case of liquid water, light atoms, such as hydrogen in particular, deviate significantly from classical behavior even at room temperature,^{51–53} as evidenced by experimental isotope effect studies which demonstrated softening of the liquid structure (when comparing H_{2}O to D_{2}O),^{54} hence neglecting the quantum nature of the nuclei in liquid water leads to overstructuring in the predicted radial distribution functions.^{52,55–57}

A commonly adopted method to alleviate the self-interaction error in GGA-DFT is the use of hybrid XC functionals, wherein a fraction of the exact exchange (Exx) energy is included in the density functional approximation. Due to the relatively high computational cost associated with these XC functionals, applications of hybrid-based DFT over the past several years have been mostly restricted to small gas-phase clusters of water,^{58–62} although recently hybrid functionals have been applied in the study of several crystalline phases of ice.^{43–46,63} In comparison to GGAs, these studies demonstrated that the energetic, structural, and vibrational properties of these systems, as predicted by hybrid DFT calculations, are generally in closer agreement with the available experimental data.^{43–46,60–63} Although applications of hybrid functionals to liquid water have been relatively rare in the literature, it was found that PBE0^{64,65} greatly improves upon the accuracy of the vibrational spectrum of liquid water.^{29,35} In the same breath, however, ambiguities exist and still remain concerning the effects of Exx on the structure of liquid water; while some studies have inferred a softening of the structure with hybrid functionals over GGAs,^{23,29,35} others have found the effects of Exx to be negligible in this regard.^{25,41} As a further complication, the increased computational cost associated with hybrid XC functionals has restricted most of the AIMD simulations performed at this level of theory to small system sizes and/or relatively short simulation times, both of which have the potential to prohibit an accurate assessment of the effects of Exx on the structural properties of liquid water. With that being said, it is generally accepted that more in-depth studies will be required in order to definitively ascertain the effects of Exx on the microscopic structure of liquid water—and this represents one of the main goals of the research reported in this paper.

In addition, both hybrid and GGA XC functionals lack the ability to describe vdW/dispersion interactions, which arise from non-local dynamical electron correlation and have a substantial effect on the microscopic structure of water in both the solid and liquid phases. In this regard, the explicit inclusion of pairwise-additive vdW interactions in DFT has been shown to significantly improve upon the theoretical description of the transition pressures among the high-pressure phases of ice^{43,46,47} and the predicted equilibrium density of liquid water.^{39–41} Many recent studies have concluded that the structure of liquid water significantly softens when vdW interactions are accounted for in the underlying XC potential; however, the extent to which these non-local vdW forces affect the structure of liquid water is largely dependent upon the given approach utilized to facilitate vdW-inclusive DFT.^{32–41} For instance, the pairwise-additive *C*_{6}/*R*^{6} based vdW correction approaches of Grimme,^{66,67} when used in conjunction with the popular PBE^{68} and BLYP^{69,70} GGA functionals, tend to reduce the intensity of the first maximum in the oxygen-oxygen radial distribution function (*g*_{OO}(*r*)) of liquid water over a fairly wide range, i.e., by approximately 5% for PBE-D^{39} and 10%–17% for BLYP-D.^{32,33,39} Other popular vdW-inclusive DFT approaches include the so-called vdW-DF functionals,^{71} which incorporate vdW interactions via the use of explicit non-local correlation functionals, have demonstrated severe shortcomings in reproducing the second coordination shell in ambient liquid water.^{34,36} Again, an accurate assessment of the effects of non-local vdW/dispersion interactions on the microscopic structure of liquid water, in particular when these interactions are treated in conjunction with a hybrid XC functional, represents an important step in increasing our understanding of this fundamental aqueous system, and this indeed is the main focus of the research reported herein.

In order to understand and quantify the individual and collective effects of Exx and vdW on the microscopic structure of liquid water, we have systematically performed a series of AIMD simulations of liquid water at ambient conditions using a hierarchy of XC functionals with increasing accuracy. In particular, the sequence of XC functionals employed in this work includes the standard semi-local GGA of Perdew, Burke, and Ernzerhof^{68} (PBE), the corresponding PBE0^{64,65} hybrid, and self-consistent (SC) dispersion-corrected analogs^{72} of each, i.e., PBE+TS-vdW(SC) and PBE0+TS-vdW(SC), based on the Tkatchenko-Scheffler^{73} density-dependent vdW/dispersion functional (TS-vdW). To successfully perform these AIMD simulations, we have employed a linear scaling O(*N*) algorithm to compute the exact exchange energy,^{74} which allows for relatively long (*t*_{sim} > 25 ps) hybrid functional based simulations using large system sizes (e.g., the explicit treatment of simulation cells including 128 water molecules). In addition, we have also developed and utilized a linear scaling O(*N*) self-consistent implementation of the TS-vdW dispersion correction, denoted by TS-vdW(SC) herein, which provides a framework for computing atomic *C*_{6} dispersion coefficients as an explicit functional of the charge density, thereby accounting for the local chemical environment surrounding each atom.^{73} In this paper, we demonstrate that the individual and collective effects of Exx and vdW interactions significantly soften the structure of liquid water and predict structural properties that are in closer agreement with the available experimental data as compared to the now well-established predictions of the PBE GGA functional.

Since it is also now known that a classical treatment of the nuclear degrees of freedom has been deemed insufficient for obtaining a quantitatively accurate description of the microscopic structure of liquid water, we have also performed an additional AIMD simulation at the slightly elevated temperature of 330 K, a technique suggested by the lowest-order perturbative expansion in ℏ of the free energy, which leads to a momentum distribution of the particles that is still Gaussian (as in classical physics) but with an increased (mass-dependent) temperature for each particle.^{75} Water molecules are composed of oxygen and hydrogen nuclei, and a uniform temperature increase of ≈30 K is in fact a crude approximation for the NQE in liquid water; however, this approximate treatment has been found to mimic the effect of quantum nuclei on the *g*_{OO}(*r*) of liquid water, both in *ab initio*^{52} and empirical potential simulations.^{56,57} By approximately accounting for the NQE as described above in AIMD simulations at the PBE0+TS-vdW(SC) level of theory, which is the most accurate XC functional considered in this work, we have obtained quantitative agreement with scattering experiments in the *g*_{OO}(*r*) and significant systematic improvements in many other structural properties (i.e., the oxygen-hydrogen radial distribution function, local tetrahedrality and the oxygen-oxygen-oxygen angular distribution function, characterization of the underlying hydrogen bonding network) as well as electrostatic properties (such as the molecular dipole moment) of liquid water.

The remainder of the paper is organized as follows. In Sec. II, we describe the computational details of the AIMD simulations performed herein. In Sec. III, we present the theoretical methodologies developed and utilized in this work to study the structural properties of liquid water at ambient conditions. Section IV contains an in-depth discussion of the results of these simulations and a comparative analysis with the currently available theoretical and experimental literature on ambient liquid water. The paper is then completed with Sec. V, which provides some brief conclusions as well as the future outlook of AIMD simulations of liquid water.

## II. SIMULATION DETAILS

In this work, we have systematically performed DFT-based AIMD simulations of ambient liquid water using a hierarchy of different XC functionals. The sequence of XC functionals employed herein includes the standard semi-local PBE-GGA,^{68} the corresponding hybrid PBE0^{64,65} which includes 25% exact exchange,^{76} and the self-consistent (SC) dispersion-corrected analogs^{72} thereof, i.e., PBE+TS-vdW(SC) and PBE0+TS-vdW(SC), based on the Tkatchenko-Scheffler^{73} density-dependent vdW/dispersion functional (TS-vdW). Further details of the theoretical methods employed in this work are provided directly below in Sec. III.

Each of these AIMD simulations of liquid water was performed in the canonical (*NVT*) ensemble using periodic simple cubic simulation cells with lattice parameters set to reproduce the experimental density of liquid water at ambient conditions. The Car-Parrinello (CP)^{6} equations of motion for the nuclear and electronic degrees of freedom were integrated using the standard Verlet algorithm and a time step of 4.0 a.u. (≈0.1 fs). The ionic temperatures were controlled with Nosé-Hoover chain thermostats, each with a chain length of 4.^{77} To achieve rapid equipartition of the thermal energy, we employed one Nosé-Hoover chain thermostat per atom (i.e., the so-called “massive” Nosé-Hoover thermostat) and also rescaled the fictitious thermostat masses by the atomic masses, so that the relative rates of the thermostat fluctuations were inversely proportional to the masses of the atoms to which they were coupled.^{78} The interactions between the valence electrons and the ions (consisting of the nuclei and their corresponding frozen-core electrons) were treated with Troullier-Martins type norm-conserving pseudopotentials.^{79} The electronic wavefunctions were expanded using a plane wave basis set with a kinetic energy cutoff of 71 Ry, a basis set expansion which was deemed sufficiently converged for the accurate determination of structural properties.^{80} To ensure an adiabatic separation between the electronic and nuclear degrees of freedom in the CP dynamics, we used a fictitious electronic mass of 300 a.u., which was found to be a reasonable choice for the simulation of water,^{13} and the nuclear mass of deuterium for each hydrogen atom. Mass preconditioning was applied to all Fourier components of the electronic wavefunctions having a kinetic energy greater than 3 Ry.^{81}

Using a simulation cell with *L* = 12.4 Å, four AIMD simulations were performed on (H_{2}O)_{64} at 300 K using the PBE, PBE0, PBE+TS-vdW(SC), and PBE0+TS-vdW(SC) XC functionals. Since it is now well-known that a classical treatment of the nuclear degrees of freedom is insufficient for a quantitatively accurate description of the microscopic structure of liquid water at room temperature, we have also performed an additional AIMD simulation on (H_{2}O)_{128} with a temperature of 330 K at the PBE0+TS-vdW(SC) level, which is the most accurate XC functional considered herein. As mentioned above, an increase of ≈30 K in the simulation temperature has been found to mimic the NQE in structural quantities such as the oxygen-oxygen radial distribution function (*g*_{OO}(*r*)) in both DFT^{52} and force-field^{56,57} based MD simulations.

All calculations reported in this work utilized a modified development version of the Quantum ESPRESSO (QE) software package.^{82} All of the simulations were initially equilibrated for approximately 2 ps and then continued for at least an additional 20 ps for data collection. The largest AIMD simulation performed herein (i.e., (H_{2}O)_{128} at the PBE0+TS-vdW(SC) level of theory) employed 1024 cores on the Cray XE6 platform at the National Energy Research Scientific Computing Center (NERSC) for ≈5 days to carry out 1 ps of simulation.

## III. THEORETICAL METHODS

### A. Hybrid XC functionals

The utilization of hybrid XC functionals requires significant additional computational cost relative to GGA-DFT-based AIMD simulations. To meet these additional computational demands and thereby allow for large-scale AIMD simulations based on hybrid XC functionals, we have employed a linear scaling O(*N*) exact exchange algorithm, which has been developed^{74} and extensively optimized in our research group. The computational savings afforded by this algorithm originate from the maximally localized Wannier function (MLWF)^{83} representation of the occupied subspace,

In this approach, the (closed-shell) exact exchange energy,

and orbital-dependent “force” acting on the electronic wavefunctions in the CP dynamics,^{6}

can both be written in terms of *v*_{ij}(**r**), the electrostatic potential generated by the pair density,

The potential *v*_{ij}(**r**) is the solution of the Poisson equation,

subject to the boundary condition that *v*_{ij}(**r**) → 0 at infinity. By virtue of the localized character of

*v*

_{ij}(

**r**) at sufficient distances from the center of the pair density is exactly represented by a multipolar expansion,

in which

*R*

_{PE}), i.e., for |

**r**| <

*R*

_{PE}. Eq. (4) is then solved numerically in real-space (i.e., on the dense real-space grid) within the sphere defined by

*R*

_{PE}, with the boundary condition provided by Eq. (5) for |

**r**| =

*R*

_{PE}.

To efficiently compute the exact exchange energy at every AIMD step, the required MLWFs were evaluated by minimizing the total spread functional,

^{−6}in the optimal spread

^{2}during the PBE0 and PBE0+TS-vdW(SC) simulations performed in this work. For a given pair of occupied Kohn-Sham electronic states, the exact exchange energy was computed when the corresponding MLWF centers were within 3.70 Å. The Poisson equation cutoff radius (

*R*

_{PE}) above was defined with respect to the center of a given pair density

*l*

_{max}= 6, which has previously been found to be sufficiently accurate in computing exact exchange energies in a variety of systems.

^{74}This procedure provides the exact exchange energy with an accuracy of at least 10

^{−4}hartree/molecule, as was tested by systematically increasing the

*R*

_{PE}cutoff radius and including pairs of MLWFs up to distances of half the simulation cell length.

### B. Self-consistent dispersion-corrected XC functionals

To account for the non-local vdW/dispersion interactions in our AIMD simulations of liquid water, we developed a fully self-consistent (SC) O(*N*) implementation^{72} of the density-dependent dispersion correction of Tkatchenko and Scheffler^{73} (TS-vdW) into a development version of QE. The TS-vdW energy is constructed as an explicit functional of the charge density and is written as a sum over atomic pair energies for atoms *A* and *B* as

in which *f*_{AB}[ρ(**r**)] is a Fermi-type damping function (used to ensure that the vdW correction goes continuously to zero at short interatomic distances) that explicitly depends on the charge density through the vdW radii of atoms *A* and *B*,

*C*

_{6, AB}[ρ(

**r**)] is the effective heteronuclear dispersion coefficient that is also dependent on the charge density via the static dipole polarizabilities and homonuclear dispersion coefficients of atoms

*A*and

*B*,

*R*

_{AB}is the scalar distance between atoms

*A*and

*B*,

*R*

_{AB}= |

**R**

_{A}−

**R**

_{B}|. As such, the TS-vdW correction provides a theoretical framework which accounts for local chemical environment effects, such as hybridization, exchange-correlation, and Pauli repulsion, in the calculation of the vdW energy and forces, and is therefore not strictly pairwise-additive

^{85}as, e.g., the DFT+D approaches of Grimme.

^{66,67}

The implementation utilized herein in conjunction with the PBE and PBE0 XC functionals, namely PBE+TS-vdW(SC) and PBE0+TS-vdW(SC), respectively, allows us to self-consistently account for the changes in the total energy and charge density (*E*[ρ(**r**)] = *E*_{DFT}[ρ(**r**)] + *E*_{TS-vdW}[ρ(**r**)]) that arise from the fact that the vdW/dispersion correction itself is an explicit correlation functional of the charge density. The fully self-consistent implementation of this energy expression for use in AIMD simulations therefore requires a corresponding electronic potential and “force” acting on the electronic wavefunctions in the CP dynamics,^{6}

in which *V*_{TS-vdW}(**r**) is the local vdW/dispersion potential derived via

which is added to the local XC potential, *V*_{XC}(**r**), at each MD step.

During the AIMD simulations at the PBE+TS-vdW(SC) and PBE0+TS-vdW(SC) levels of theory, the *s*_{R} scaling parameter in the damping function above, which determines the onset of the TS-vdW correction for a given XC functional, was set to 0.94 and 0.96, respectively, as recommended in Ref. 73. The free atom reference quantities for the static dipole polarizabilities (

^{3}and

^{6}for oxygen and

^{3}and

^{6}for hydrogen, and were taken from the database of Chu and Dalgarno

^{86}(as recommended in Ref. 73). These quantities have an estimated error of approximately 1%–3% and were computed using self-interaction corrected time-dependent DFT (SIC-TDDFT) and then scaled to reproduce highly accurate many-body quantum mechanical references for the alkali, alkaline earth, and noble gas atoms. The free atom reference quantities for the vdW radii (

^{87}The TS-vdW energy (and resultant ionic forces and electronic potential) were computed in real-space by explicitly summing over the atoms in the simulation and neighboring periodic cells to ensure convergence in

*E*

_{TS-vdW}to at least 10

^{−5}hartree/molecule.

## IV. RESULTS AND DISCUSSION

### A. Comparative analysis of the oxygen-oxygen structure factor

We begin our analysis by comparing the oxygen-oxygen structure factor, *S*_{OO}(*Q*), obtained from the highest quality AIMD simulations performed herein at the PBE0+TS-vdW(SC) (330 K) level of theory and various X-ray scattering measurements (see Fig. 1(a)). Such X-ray scattering experiments directly probe the electron charge density (which is centered on the oxygen atoms in the case of water), thereby allowing for *S*_{OO}(*Q*) to be accurately obtained from the directly measured scattering intensity, *I*(*Q*), and knowledge of the corresponding form factor, *F*^{2}(*Q*), via *S*(*Q*) = *I*(*Q*)/*F*^{2}(*Q*), in which *Q* represents the change in the wavevector between the incident and scattered radiation. Once *S*(*Q*) has been experimentally determined, the so-called pair correlation or radial distribution function, *g*(*r*), which is the most commonly utilized measure of the microscopic structure of a liquid in real space, can then be obtained via the inverse Fourier transform, i.e.,

wherein *Q*_{max} is the maximum value of *Q* in the given scattering experiment and ρ is the atomic number density. Hence, the process of obtaining an accurate *g*(*r*) is highly sensitive to the *Q*_{max} utilized in a given X-ray scattering experiment, which still represents the main technical challenge for scattering experiments to date, and the *a priori* knowledge of the corresponding form factor, *F*^{2}(*Q*), required for deriving *S*(*Q*) from the experimentally measured scattering intensity.

As shown in Fig. 1(a), the *S*_{OO}(*Q*) taken from various X-ray scattering experiments^{88–90} with increasing values of *Q*_{max} are very similar across the range of *Q* values accessible in a given scattering experiment; despite this fact, each of these X-ray scattering experiments yields significantly different *g*_{OO}(*r*), especially with respect to the intensity and position of the first peak.^{88–90} In this regard, we note that these differences in the experimentally derived *g*_{OO}(*r*) are almost entirely due to the aforementioned *Q*_{max} limitation inherent to the X-ray scattering experiments and not due to errors in the corresponding *F*^{2}(*Q*). To justify this claim, we first computed the difference between the *F*^{2}(*Q*) most commonly utilized in these X-ray scattering experiments, which are obtained from high-level quantum chemistry calculations on a single gas-phase water molecule,^{91,92} and the corresponding quantity at the DFT level of theory. Finding nearly quantitative agreement in this quantity at both levels of theory for a single gas-phase water molecule, as also demonstrated by Krack *et al.*,^{93} we then computed the difference between the DFT-based *F*^{2}(*Q*) for the single gas-phase water molecule and the *F*^{2}(*Q*) corresponding to condensed-phase liquid water obtained from the DFT-based AIMD simulations performed herein.^{94} In performing this analysis, we found that the differences between the form factors computed with respect to gas- and condensed-phase water are extremely small (i.e., to within 1%–2%) and are mainly concentrated in the low-*Q* region (i.e., for *Q* < 2 Å^{−1}); hence, such differences would be essentially negligible when computing the corresponding *g*_{OO}(*r*).

Using an AIMD trajectory obtained at the PBE0+TS-vdW(SC) (330 K) level of theory, we have computed *S*_{OO}(*Q*) as

in which *N* is the number of oxygen atoms in the simulation cell, *N*_{Q} is the degeneracy of the wavevector shell of length *Q* = |**Q**|, and **R**_{i} is the real-space position vector of *i*th oxygen atom. Unlike the experimentally derived *S*_{OO}(*Q*), the accuracy of the computed *S*_{OO}(*Q*) in the low-*Q* region is subject to a large degree of numerical noise, a consequence of the fact that there is a relatively small number of **Q** vectors present in this sector. To combat this fact (and to minimize finite size effects), we performed the highest quality AIMD simulations in this work (at the PBE0+TS-vdW(SC) (330 K) level) using a periodic simple cubic box containing 128 water molecules (*L* = 15.7 Å), wherein the smallest accessible *Q* value is 0.4 Å^{−1}, which helps in improving the accuracy of the theoretically determined *S*_{OO}(*Q*) in the low-*Q* region.

From Fig. 1(a), it is clear that the theoretically determined *S*_{OO}(*Q*) is in nearly quantitative agreement with the experimental results across the entire *Q* region accessible to the X-ray scattering experiments, with only a slightly noticeable shift towards higher *Q* values. This level of agreement between theory and experiment has been achieved via the utilization of the PBE0+TS-vdW(SC) XC functional in AIMD simulations performed at 330 K—a level of theory which (i) reduces the level of self-interaction with respect to standard GGA-DFT by including 25% exact exchange, (ii) accounts for non-local vdW/dispersion interactions by a self-consistent density-dependent *C*_{6}/*R*^{6} correction, and (iii) approximately corrects for nuclear quantum effects with a 30 K increase in the simulation temperature (see Sec. III for more details on the theoretical methodologies employed in this work). In what follows, we will focus on the corresponding *g*_{OO}(*r*), a real-space quantity which allows for a more straightforward comparative analysis and delineation of the individual and collective effects stemming from each of these improvements in the underlying XC functional.

### B. Comparative analysis of the oxygen-oxygen radial distribution function

Within the last year, it has been shown that the accuracy of the oxygen-oxygen radial distribution function, *g*_{OO}(*r*), directly obtained from X-ray scattering experiments can be systematically improved upon via the use of larger *Q*_{max} values, which has yielded essentially converged results (i.e., to within 1%) for the intensity and position of the first peak in the *g*_{OO}(*r*) of liquid water with *Q*_{max} > 20.0 Å^{−1}.^{90} In Fig. 1(b), this latest X-ray scattering data,^{90} which is now in quantitative agreement with the *g*_{OO}(*r*) obtained from empirical potential structural refinement (EPSR) based on joint X-ray/neutron data by Soper and Benmore,^{54} are both compared against the highest quality AIMD simulations performed herein at the PBE0+TS-vdW(SC) (330 K) level of theory. From this figure, one immediately notices that the *g*_{OO}(*r*) obtained from PBE0+TS-vdW(SC) (330 K) AIMD simulations is in excellent agreement with both of the experimentally derived results, a key finding of this work which will be discussed in greater detail below.

As mentioned above, this level of nearly quantitative agreement with experiment has been achieved via the collective treatment of exact exchange (via the PBE0 hybrid functional), non-local vdW/dispersion interactions (via a self-consistent density-dependent *C*_{6}/*R*^{6} correction) and approximate nuclear quantum effects (via a 30 K increase in the simulation temperature). To illustrate this point, we will now compare the individual and collective effects of Exx and vdW on the predicted *g*_{OO}(*r*) of liquid water by considering Fig. 2, which provides the *g*_{OO}(*r*) results obtained from AIMD simulations based on each of the different XC functionals utilized in this work and Table I, which summarizes the positions and intensities at various key points in each *g*_{OO}(*r*), the oxygen-oxygen coordination numbers, *n*_{OO}, associated with the first shell in liquid water, as well as other structural properties which will be discussed in detail in Secs. IV D–IV F.

Method . | T
. | $r^{\rm max}_{1}$ $r1 max $
. | $g^{\rm max}_{1}$ $g1 max $
. | $r^{\rm min}_{1}$ $r1 min $
. | $g^{\rm min}_{1}$ $g1 min $
. | $r^{\rm max}_{2}$ $r2 max $
. | $g^{\rm max}_{2}$ $g2 max $
. | n_{OO}
. | q
. | n_{HB}
. | μ . | D
. |
---|---|---|---|---|---|---|---|---|---|---|---|---|

PBE | 300 | 2.69 | 3.28 | 3.28 | 0.37 | 4.44 | 1.44 | 4.03 | 0.78 | 3.83 | 3.19 ± 0.29 | 0.020 |

PBE+TS-vdW(SC) | 300 | 2.71 | 2.99 | 3.27 | 0.54 | 4.40 | 1.33 | 4.01 | 0.74 | 3.74 | 3.12 ± 0.31 | 0.044 |

PBE0 | 300 | 2.71 | 2.96 | 3.30 | 0.53 | 4.39 | 1.35 | 4.10 | 0.73 | 3.71 | 3.00 ± 0.27 | 0.067 |

PBE0+TS-vdW(SC) | 300 | 2.72 | 2.76 | 3.31 | 0.70 | 4.47 | 1.20 | 4.22 | 0.69 | 3.62 | 2.96 ± 0.27 | 0.098 |

PBE0+TS-vdW(SC) | 330 | 2.74 | 2.51 | 3.33 | 0.84 | 4.44 | 1.22 | 4.32 | 0.64 | 3.48 | 2.91 ± 0.28 | 0.146 |

Expt.^{90} | 295 | 2.80 | 2.57 | 3.45 | 0.84 | 4.50 | 1.12 | 4.3(1) | … | … | … | … |

Expt.^{54} | 298 | 2.76 | 2.62 | 3.42 | 0.82 | 4.43 | 1.14 | 4.67(5) | 0.58 | … | … | … |

Expt.^{95} | … | … | … | … | … | … | … | … | … | … | 2.9 ± 0.6 | … |

Expt.^{96–98} | … | … | … | … | … | … | … | … | … | … | … | 0.23–0.24 |

Method . | T
. | $r^{\rm max}_{1}$ $r1 max $
. | $g^{\rm max}_{1}$ $g1 max $
. | $r^{\rm min}_{1}$ $r1 min $
. | $g^{\rm min}_{1}$ $g1 min $
. | $r^{\rm max}_{2}$ $r2 max $
. | $g^{\rm max}_{2}$ $g2 max $
. | n_{OO}
. | q
. | n_{HB}
. | μ . | D
. |
---|---|---|---|---|---|---|---|---|---|---|---|---|

PBE | 300 | 2.69 | 3.28 | 3.28 | 0.37 | 4.44 | 1.44 | 4.03 | 0.78 | 3.83 | 3.19 ± 0.29 | 0.020 |

PBE+TS-vdW(SC) | 300 | 2.71 | 2.99 | 3.27 | 0.54 | 4.40 | 1.33 | 4.01 | 0.74 | 3.74 | 3.12 ± 0.31 | 0.044 |

PBE0 | 300 | 2.71 | 2.96 | 3.30 | 0.53 | 4.39 | 1.35 | 4.10 | 0.73 | 3.71 | 3.00 ± 0.27 | 0.067 |

PBE0+TS-vdW(SC) | 300 | 2.72 | 2.76 | 3.31 | 0.70 | 4.47 | 1.20 | 4.22 | 0.69 | 3.62 | 2.96 ± 0.27 | 0.098 |

PBE0+TS-vdW(SC) | 330 | 2.74 | 2.51 | 3.33 | 0.84 | 4.44 | 1.22 | 4.32 | 0.64 | 3.48 | 2.91 ± 0.28 | 0.146 |

Expt.^{90} | 295 | 2.80 | 2.57 | 3.45 | 0.84 | 4.50 | 1.12 | 4.3(1) | … | … | … | … |

Expt.^{54} | 298 | 2.76 | 2.62 | 3.42 | 0.82 | 4.43 | 1.14 | 4.67(5) | 0.58 | … | … | … |

Expt.^{95} | … | … | … | … | … | … | … | … | … | … | 2.9 ± 0.6 | … |

Expt.^{96–98} | … | … | … | … | … | … | … | … | … | … | … | 0.23–0.24 |

The first point to note here is the fact that the *g*_{OO}(*r*) obtained from the PBE based AIMD simulations performed herein are in good agreement with the previous reports in the literature, wherein the intensity of the first maximum,

^{13,14,26}followed by a first minimum that is too deep with respect to experiment. This level of overstructuring in the

*g*

_{OO}(

*r*) of liquid water generated at the PBE level of theory is fairly well-known at this point and is a manifestation of the key limitations of GGA-DFT in the theoretical description of liquid water, namely, the presence of self-interaction error and the lack of non-local vdW/dispersion interactions in the underlying XC potential. In other words, these inherent limitations in the underlying XC potential essentially yield a system that is more representative of deeply supercooled liquid water when the corresponding PBE-based AIMD simulations are performed at ambient conditions.

By including a fraction (i.e., 25%) of exact exchange via the use of the PBE0 hybrid functional, an immediate softening of the *g*_{OO}(*r*) is observed in which the intensity of the first maximum,

*g*

_{OO}(

*r*) generated at the PBE level of theory (see Fig. 2 and Table I). These noticeable differences in the

*g*

_{OO}(

*r*) are primarily due to the fact that the PBE0 hybrid functional, which reduces the amount of self-interaction error in the underlying XC potential and thereby corrects for the overly delocalized electron density at the GGA-DFT level of theory, weakens the hydrogen bond strength in liquid water relative to PBE, a finding that is also in agreement with the previously observed improvements in the infrared spectrum of PBE0-based liquid water.

^{29}In fact, this effect has also been demonstrated in previous studies on both gas-phase water clusters

^{60–62}and ice,

^{43,45,46}in which a significant reduction in the strength of the hydrogen bonds was observed in each of these aqueous systems with the use of hybrid XC functionals. By weakening the individual hydrogen bonds, simulations at the PBE0 level of theory yield a marked increase in the population of water molecules located in the interstitial region, i.e., the region between the first and second coordination shells. However, this trend of destructuring or softening in the

*g*

_{OO}(

*r*) of liquid water observed here is in contrast with a few previous studies at the PBE0 level of theory, in which no effect was observed on the

*g*

_{OO}(

*r*) when exact exchange was accounted for over the standard PBE functional.

^{25,41}In this regard, we note the existence of a theoretical difference in the treatment of the exact exchange contributions between this work, in which the exact exchange contributions were exactly computed numerically

^{74}(see Sec. III), and that of Guidon

*et al.*,

^{25}in which the Coulomb potential in the exact exchange integrals was replaced by a finite-range approximant in terms of the complementary error function, erfc(

*r*).

By including a self-consistent treatment of the non-local vdW/dispersion interactions via the density-dependent Tkatchenko-Scheffler dispersion correction,^{73} we observed a very similar effect on the *g*_{OO}(*r*) of liquid water; in fact, the radial distribution functions obtained at the PBE0 and PBE+TS-vdW(SC) levels of theory are almost identical over the entire distance range considered (see Fig. 2 and Table I). Unlike exact exchange, however, vdW forces are non-directional and strengthen the interactions between a given water molecule and the water molecules in its first and second coordination shells, which actually causes the water molecules in the second coordination shell to move inward and populate the interstitial region. As a result, the hydrogen bonding (see Sec. IV E) in the first coordination shell is weakened, causing the water molecules in the first coordination shell to be pushed outward as well. Here vdW forces tend to stabilize such disordered configurations—structures that would be energetically disfavored at the PBE level of theory—and this effect leads to the same net result in the *g*_{OO}(*r*) as that observed in PBE0-based liquid water. In this regard, there is only a subtle difference in the *g*_{OO}(*r*) of liquid water between the PBE+TS-vdW(SC) and PBE0 levels of theory: for PBE+TS-vdW(SC), the position of the first minimum,

^{99}

*n*

_{OO}, actually decreases slightly when compared to PBE (4.01 vs. 4.03), whereas at the PBE0 level,

*n*

_{OO}increases to a value of 4.10 since

Several studies performed to date have concluded that the structure of liquid water significantly softens when vdW interactions are accounted for in AIMD simulations; however, the extent to which non-local vdW/dispersion forces affect the structure of liquid water is largely dependent on the given approach utilized to facilitate vdW-inclusive DFT.^{32–41} For instance, the most commonly utilized approach involves the pairwise-additive *C*_{6}/*R*^{6} vdW/dispersion correction of Grimme^{66,67} in conjunction with the popular PBE and BLYP GGA functionals, which tends to reduce the intensity of the first maximum in the *g*_{OO}(*r*) of liquid water over a fairly wide range, i.e., by approximately 5% for PBE-D^{39} and 10%–17% for BLYP-D.^{32,33,39} Other popular vdW-inclusive DFT approaches for studying liquid water include the so-called vdW-DF functionals,^{71} which incorporate vdW/dispersion interactions via the use of explicit non-local correlation functionals. This class of vdW-inclusive functionals have severe shortcomings in reproducing the second coordination shell in the *g*_{OO}(*r*), especially when employed in conjunction with the revPBE exchange functional;^{34,36} the use of the PBE exchange functional instead recovers this limitation to some extent^{34,35,100} and reduces the intensity of the first maximum by ≈25% with respect to PBE.^{35} For comparison, the self-consistent implementation of the density-dependent Tkatchenko-Scheffler^{73} vdW/dispersion correction employed in this work, PBE+TS-vdW(SC), yields an 8.8% reduction in the intensity of the first maximum with respect to PBE.

The collective effect of treating exact exchange and non-local vdW/dispersion interactions, as depicted by the AIMD simulations of liquid water at the PBE0+TS-vdW(SC) level of theory, yields a *g*_{OO}(*r*) that is even closer to the experimental findings. In this case, both of these improvements to the underlying XC functional work together to reduce the intensity of the first maximum,

*E*

_{vdW}=

*E*

_{vdw}[

*n*(

**r**)]). In the same breath, the collective increase in the first minimum of 0.33 at the PBE0+TS-vdW(SC) level of theory is exactly equal to the sum of the individual increases observed when utilizing PBE0 (0.16) and PBE+TS-vdW(SC) (0.17) alone, reflecting the fact that both of these effects definitively work in unison to weaken the hydrogen bond network and increase the relative population of water molecules in the interstitial region. We note in passing that similar conclusions can be drawn by considering the resultant shifts in the positions of the first maximum and minimum of the PBE0+TS-vdW(SC)

*g*

_{OO}(

*r*): the collective shift in the position of the first maximum of 0.03 Å is only slightly smaller than the sum of the individual shifts of 0.02 Å (PBE0) and 0.02 Å (PBE+TS-vdW(SC)), while the increase in the position of the first minimum of 0.03 Å represents an interesting collective effect between PBE0 (which increased

Although a more detailed analysis of the nature and extent of the hydrogen bonding network generated from the AIMD simulations performed in this work will be discussed in Sec. IV E, another interesting point can be made here concerning the contributions arising from select neighboring water molecules to the overall *g*(*r*) as a function of the underlying XC potential. In Fig. 3, the contributions to the *g*_{OO}(*r*) of liquid water arising from the fourth (4th) and fifth (5th) neighboring water molecules are plotted for each of the XC functionals employed in this work. From this figure, one can immediately notice an apparent shift in the *g*_{OO}(*r*) contributions occurring between the 4th and 5th neighbors as the underlying XC functional is systematically improved by accounting for exact exchange and/or non-local vdW/dispersion interactions. For the first four neighboring water molecules (with only the 4th neighbor shown in Fig. 3 for clarity), we observe that the individual and collective effects of Exx and vdW tend to simultaneously increase the position and decrease the intensity of this contribution to the overall *g*_{OO}(*r*); on the other hand, we observe the exact opposite for the fifth and sixth neighboring water molecules (with only the 5th neighbor shown in Fig. 3 for clarity), which tend to simultaneously decrease the position and increase the intensity of their contribution to the overall *g*_{OO}(*r*). Both of these findings are indicative of a deviation from the perfect tetrahedral hydrogen bonding network observed in ice (and to a large extent in PBE liquid water) and a relative increase in the population of water molecules in the interstitial region. As such, this clearly illustrates how these systematically improved XC functionals soften the structure of liquid water with respect to GGA-DFT and thereby generate oxygen-oxygen radial distribution functions that are increasingly in better agreement with the experimental data.

Even further improvements with respect to the experimental findings is possible when nuclear quantum effects are mimicked in the *g*_{OO}(*r*) of liquid water by performing the simulation at the PBE0+TS-vdW(SC) level of theory with a 30 K increase in the simulation temperature (see Fig. 2 and Table I). At 330 K, we observe further softening in the *g*_{OO}(*r*) of liquid water, with peak positions and intensities that are now in quantitative agreement with the *g*_{OO}(*r*) directly obtained from X-ray^{90} scattering experiments. In this regard, approximately accounting for NQE leads to a further decrease in the intensity of the first maximum of 0.25, which is now approximately 2.3%–4.2% lower than the experimental data, and a further increase in the intensity of the first minimum of 0.14, which is in exact agreement with Skinner *et al.*^{90} and only 2.4% higher than the value provided by Soper and Benmore.^{54} Similar agreement can be found when analyzing the positions of the first maximum (with an error of 0.7%–2.1%) and the first minimum (with an error of 2.6%–3.5%) with respect to the experimental findings, as well as the resulting coordination number. This level of agreement with experiment in the *g*_{OO}(*r*) of liquid water has been quite difficult to achieve to date and clearly requires an increase in the accuracy of the underlying XC functional, and at the same time, suggests the need for a quantum mechanical treatment of the nuclear degrees of freedom.

### C. Comparative analysis of the oxygen-hydrogen radial distribution function

A similar trend of systematic improvements in the theoretical description of the oxygen-hydrogen radial distribution function, *g*_{OH}(*r*), of liquid water (see Fig. 4) can be observed with the inclusion of exact exchange (via the PBE0 hybrid functional), non-local vdW/dispersion interactions (via a self-consistent density-dependent *C*_{6}/*R*^{6} correction), and approximate nuclear quantum effects (via a 30 K increase in the simulation temperature). In this case, however, the width and intensity of the first peak of the *g*_{OH}(*r*) of liquid water, which describes the covalent O–H bond, can only be described by properly accounting for nuclear quantum effects in the simulation of liquid water, i.e., the approximate treatment of nuclear quantum effects via a 30 K increase in the simulation temperature is not sufficient to capture the quantum nature of the lighter hydrogen atoms, which deviate significantly from classical behavior at room temperature.^{52,53} As such, the explicit inclusion of NQE in the *ab initio* simulations of liquid water via the Feynman discretized path integral (PI) scheme is the focus of intense ongoing research in our group and will be addressed in future work.

The largest effects arising from the individual and collective treatment of exact exchange and non-local vdW/dispersion interactions can be found in the position and intensity of the second maximum in the *g*_{OH}(*r*), which is the signature of the hydrogen bond network in liquid water. In this regard, we observe very similar trends to those reported above for the *g*_{OO}(*r*) of liquid water, in that the AIMD simulations at the PBE0 and PBE+TS-vdW(SC) levels of theory predict nearly identical *g*_{OH}(*r*) over the entire distance range considered (see Fig. 4). Considering for instance the intensity of the second maximum in the *g*_{OH}(*r*), there is an overall reduction or destructuring of approximately 0.22 (PBE0) and 0.15 (PBE+TS-vdW(SC)), which is consistent with the fact that both of these XC functionals tend to disrupt the hydrogen bond network and increase the relative population of water molecules in the interstitial region. In this case, the effect of including a fraction of exact exchange on the *g*_{OH}(*r*) is only slightly larger than the observed effect resulting from an explicit treatment of vdW/dispersion interactions. When used in combination, the PBE0+TS-vdW(SC) XC functional again displays a nearly additive effect on the *g*_{OH}(*r*) by reducing the intensity of the second maximum by 0.33, which is to be contrasted against the sum of the aforementioned individual contributions at 0.37.

An even further reduction in the intensity of the second maximum in the *g*_{OH}(*r*) is provided by PBE0+TS-vdW(SC) AIMD simulations performed at 330 K, which yield an overestimation of approximately 11.6% with respect to the joint X-ray/neutron scattering results of Soper and Benmore.^{54} We again stress here that structural quantities that directly involve the hydrogen atoms in liquid water, such as the *g*_{OH}(*r*) and *g*_{HH}(*r*), require an explicit treatment of nuclear quantum effects; hence, errors in excess of 10% as reported above are completely reasonable and expected in such quantities when generated from simulations employing classical mechanics for the nuclear equations of motion. In the same breath, structural quantities like the *g*_{OO}(*r*), which was considered in detail in Sec. IV B, were less sensitive to an approximate treatment of NQE, and nearly quantitative accuracy can be achieved for this property with a 30 K increase in the simulation temperature, provided a sufficiently accurate underlying XC potential is utilized.

### D. Tetrahedrality and the oxygen-oxygen-oxygen triplet distribution function

In order to further analyze the local arrangement of water molecules in condensed-phase liquid water, we have computed the distribution of triplet oxygen-oxygen-oxygen angles, *P*_{OOO}(θ), within the first coordination shell (see Fig. 5). In general, three-body correlation functions such as *P*_{OOO}(θ) are not directly accessible from scattering experiments; however, an experimental distribution of θ was extracted via EPSR based on joint X-ray/neutron scattering data.^{54} In this regard, we note that the aforementioned quantitative agreement between the oxygen-oxygen radial distribution functions (see Fig. 2) directly obtained via X-ray scattering experiments^{90} and EPSR based on joint X-ray/neutron scattering data^{54} supports the effectiveness of this procedure. To extract the corresponding *P*_{OOO}(θ) quantity from the AIMD simulations performed in this work, three oxygen atoms were considered as part of a given triplet if two of the oxygen atoms were within a prescribed cutoff distance from the third. Following Ref. 54, this cutoff distance was chosen to yield an average oxygen-oxygen coordination number of ≈4, and fell within the range of 3.25–3.27 Å for all of the AIMD simulations considered herein.

From Fig. 5, the EPSR-experimental *P*_{OOO}(θ) shows a broad peak around 100°, which is indicative of the fact that the local tetrahedral network in liquid water is substantially more disordered than in crystalline ice. At the PBE-GGA level of theory, we find that the peak of *P*_{OOO}(θ) is very close to the perfect tetrahedral angle of 109.5° and the overall distribution is much too narrow as compared to the EPSR-experimental results—another manifestation of the fact that liquid water generated using PBE is overstructured and the degree of local tetrahedral order is much higher than that observed in experiment. In fact, by computing the tetrahedral order parameter, *q*, as^{101}

in which θ_{ij} is the angle formed by a central given water molecule and its nearest neighbors *i* and *j*, we find that PBE yields an average *q* value of 0.78, which is indeed much higher than the estimated EPSR-experimental value of *q* = 0.576^{54} (see Table I). For reference, a system with perfect tetrahedral order would yield *q* = 1, whereas *q* = 0 for the ideal gas (i.e., a system with random positions).

As seen above, the degree of local tetrahedral order in liquid water can be significantly reduced when utilizing an XC potential that accounts for exact exchange and non-local vdW/dispersion interactions. As shown in Fig. 5, the individual and collective effects of Exx and vdW shift the peak of *P*_{OOO}(θ) towards lower θ values and make the overall distribution much broader, both of which result in better agreement with the EPSR-experimental triplet distribution functions. Here we note that a reduction in the tetrahedral order parameter with respect to the GGA-DFT values has been previously reported using other vdW-inclusive functionals^{33,102} with values of *q* ≈ 0.7, which is in reasonable agreement with the value of *q* = 0.74 computed herein at the PBE+TS-vdW(SC) level of theory; however, capturing this trend using several vdW-inclusive functionals has come with the deleterious cost of a nearly vanishing second coordination shell in the *g*_{OO}(*r*).^{36} At the PBE0 level, we computed a value of *q* = 0.73 for the tetrahedral order parameter, which is extremely close to the value obtained above in the PBE+TS-vdW(SC) case and in good agreement with the previous study of Zhang *et al.*^{102} When used in conjunction, the self-consistent PBE0+TS-vdW(SC) exchange-correlation functional yields a value of *q* = 0.69 for the tetrahedral order parameter, which is exactly additive in terms of the aforementioned individual effects of Exx and vdW and therefore closer to the EPSR-experimental measure of the local tetrahedrality.

As seen above in Secs. IV A–IV C, even further improvement comes with the approximate inclusion of NQE, wherein we computed a value of *q* = 0.64, which is now approximately 10.3% larger than the EPSR-experimental value. Interestingly, the shoulder present in the EPSR-experimental *P*_{OOO}(θ) distribution at approximately 60°, which is indicative of a highly distorted hydrogen bond network in the first coordination shell, only becomes noticeably visible at the PBE0+TS-vdW(SC) (330 K) level of theory. Nevertheless, this quantity is sensitive to the positions of the hydrogen atoms, and as such, will require a proper treatment of nuclear quantum effects in order to accurately capture these finer details.

### E. Hydrogen bond analysis in liquid water

The fluidity in liquid water comes from the fact that the hydrogen bonds, which create the underlying random tetrahedral network, are continuously breaking and forming. Unlike ice, in which each water molecule is involved in four hydrogen bonds with its nearest neighbors (i.e., each water molecule is simultaneously donating and accepting two hydrogen bonds), liquid water contains some fraction of broken hydrogen bonds on average.^{103} In this regard, any quantitative measure of the number of either intact or broken hydrogen bonds in liquid water is somewhat ambiguous, since the notion of a hydrogen bond itself is not uniquely defined. However, qualitative agreement between many proposed definitions for intact hydrogen bonds (i.e., those based on geometry, topology, and even electronic structure) have been deemed satisfactory over a wide range of thermodynamic conditions.^{104} Since we are interested in the qualitative change in the statistics of the number and types of hydrogen bonds in liquid water as a function of the individual and collective treatment of exact exchange, non-local vdW/dispersion interactions, and approximate nuclear quantum effects, we have chosen to utilize the popular hydrogen bond definition proposed by Luzar and Chandler,^{105} wherein the defining parameters for an intact hydrogen bond include both a distance and angular criterion, namely, *R*_{OO} < 3.5 Å and β < 30°, with β ≡ ∠O_{A}⋅⋅⋅O_{D}−H_{D}.

Using this definition, we computed the average number of intact hydrogen bonds per water molecule^{106} for each of the AIMD simulations performed in this work and reported these numbers in Table I. From this data, we first observed that the average number of intact hydrogen bonds per water molecule decreases in the following order as the underlying XC potential is systematically improved: PBE (3.83), PBE+TS-vdW(SC) (3.74), PBE0 (3.71), PBE0+TS-vdW(SC) (3.67), and PBE0+TS-vdW(SC) at 330 K (3.62). This decrease in the number of intact hydrogen bonds per water molecule (or corresponding increase in the number of broken hydrogen bonds per water molecule) is representative of a larger degree of disorder in the local tetrahedral network in liquid water, and allows for an increase in the relative population of water molecules in the interstitial region, a feature that has already manifested itself in several of the structural quantities considered above, in particular in the reduction of the intensity of the first minimum in the *g*_{OO}(*r*) in Fig. 2.

As mentioned above, this trend also reflects an increase in the fluidity of *ab initio* liquid water generated with XC functionals that account for exact exchange and dispersion interactions; as the underlying XC functional is systematically improved, the computed diffusion coefficients (*D*) continue to increase towards the experimental values^{96–98} in conjunction with a simultaneous growth in the percentage of broken hydrogen bonds (see Table I). Extracted from the slope of the molecular mean squared displacement (MSD) vs. simulation time, these diffusion coefficients are given for illustrative purposes only, as the AIMD simulations reported herein utilized massive Nosé-Hoover thermostatting, which can affect the computation of dynamical properties in an uncontrolled manner. Moreover, diffusion coefficients computed from periodic simulation cells that include 64, or even 128, water molecules are subject to non-negligible finite size effects.^{26} In this regard, we should also note that the explicit treatment of NQE should also increase the self-diffusion coefficient of liquid water at ambient conditions. According to the estimates provided in Ref. 107, this increase can be as large as a factor of 1.5 relative to the classical result.

To provide a more detailed look into the hydrogen bonding in liquid water, we also performed a decomposition of the intact hydrogen bonds per water molecule into acceptor-(A) and donor-(D) types, i.e., O_{A}⋅⋅⋅H_{D}−O_{D}, in Fig. 6. From this figure, one immediately notices that the largest percentage of the intact hydrogen bonds for each of the AIMD simulations considered herein are of the A_{2}D_{2} type—the hydrogen bonding motif which represents the traditional picture of a water molecule accepting and donating two hydrogen bonds. At the PBE level of theory, the percentage of A_{2}D_{2} type hydrogen bonds is a majority at 75.4%, a quantity which is reduced to just less than 50% when Exx and vdW effects are accounted for in the underlying XC potential at the PBE0+TS-vdW(SC) (330 K) level, a trend that again reflects the increasing level of disorder in the local tetrahedral network in liquid water.

Hydrogen bond types in which a given water molecule is involved in three hydrogen bonds, as in the A_{1}D_{2} and A_{2}D_{1} types, comprise the next largest percentage of the intact hydrogen bonds in liquid water. Here we find an accompanying increase in the relative population of these hydrogen bond types as the underlying XC potential is improved from PBE to PBE0+TS-vdW(SC) (330 K); in this case, the percentage of the A_{1}D_{2} and A_{2}D_{1} hydrogen bond types were observed to approximately double from 10.7% to 19.5% and 5.9% to 12.7%, respectively. Expectedly, the number of intact hydrogen bonds in which a given water molecule is involved in less than three hydrogen bonds, as in the A_{1}D_{1} type, or more than four hydrogen bonds, as in the A_{3}D_{2} type, are less likely to occur and this trend is also reflected in Fig. 6. In this regard, we observed that the percentage of the A_{1}D_{1} hydrogen bond type approximately increased by a factor of four from 2.4% (PBE) to 9.3% (PBE0+TS-vdW(SC) (330 K)), while the percentage of the A_{3}D_{2} type remained essentially constant and therefore independent of the underlying XC potential.

This decomposition analysis of the intact hydrogen bonds as a function of the underlying XC potential again reflects the individual and collective effects of Exx, vdW, and an approximate treatment of NQE in the microscopic structure of liquid water; with Exx and the increased simulation temperature weakening the hydrogen bond strength and non-local vdW/dispersion interactions stabilizing disordered water configurations, there is a relative increase in the population of interstitial water molecules which serves to introduce an increasing degree of disorder into the underlying tetrahedral network.

### F. Dipole moment analysis in liquid water

With the inclusion of exact exchange, non-local vdW/dispersion interactions, and approximate nuclear quantum effects, the modifications to the local microscopic structure of liquid water observed in this work can also be correlated with molecular electrostatic properties which govern the strength of the hydrogen bonding network and influence solvation behavior. In this section, we investigate one such electrostatic property, namely the distribution of molecular dipole moments—the first non-zero multipole moments in the individual water molecules comprising the condensed phase. For reference, the dipole moment of a single water molecule in the gas phase is accurately known to be 1.855 D from spectroscopic measurements,^{108,109} and several DFT functionals can reproduce this experimental value to within 3%.^{35} However, there still remains a large uncertainty in the value of the molecular dipole moment in liquid water as extracted from X-ray scattering form factors (2.9 ± 0.6 D),^{95} and in this regard DFT-based AIMD has the potential to provide molecular dipole moments in the condensed phase with a larger degree of certainty.

The calculation of molecular dipole moments in condensed-phase systems requires partitioning of the electron density, which again cannot be accomplished in a uniquely defined manner; as a result, the magnitude of the molecular dipole moment in ice Ih, for instance, can vary between 2.3 and 3.1 D depending on the partitioning scheme employed.^{110} In this regard, maximally localized Wannier functions (MLWFs), which are obtained via a unitary transformation of the occupied Kohn-Sham electronic states, have been shown to be an extremely useful tool in computing molecular dipole moments in condensed-phase environments.^{35,111–114} For the case of liquid water, the MLWFs associated with different molecules have only a very minor overlap,^{111} and this fact allows for a nearly unique definition of the charges belonging to a given water molecule, thereby eliminating to a large extent the aforementioned partitioning issues when computing molecular dipole moments in the condensed-phase.^{113,114}

Using the four MLWFs corresponding to each water molecule (which represent the eight associated valence electrons), the molecular dipole moment (** μ**) for a water molecule in the condensed-phase can then be computed as

in which

**R**

_{O}are the Cartesian coordinates of the hydrogen and oxygen atoms comprising the water molecule, respectively, and

Using this prescription, we have computed the distribution of molecular dipole moments in liquid water for each of the DFT-based AIMD simulations performed in this work and plotted the resulting data in Fig. 7(a). From this figure, it is evident that the peak positions of the dipole moment distributions are shifted towards smaller values as Exx and vdW effects are accounted for in the underlying XC potential. In fact, the average magnitude of the molecular dipole moment in liquid water decreases in the following order as the underlying XC potential is systematically improved: PBE (3.19 ± 0.29 D), PBE+TS-vdW(SC) (3.12 ± 0.31 D), PBE0 (3.00 ± 0.27 D), PBE0+TS-vdW(SC) (2.96 ± 0.27 D), and PBE0+TS-vdW(SC) at 330 K (2.91 ± 0.28 D). Unlike the structural properties considered above, we note that the inclusion of exact exchange in the XC potential clearly has a larger effect on the magnitude of the molecular dipole moment in liquid water than the albeit self-consistent treatment of non-local vdW/dispersion interactions employed herein. Since the dipole moment is an electrostatic property, this result is not surprising as the hybrid PBE0 functional induces changes in the underlying local electronic structure of the individual water molecules that comprise the condensed-phase.^{35} Although each of the aforementioned functionals yields a molecular dipole moment magnitude that is within the (relatively large) error bar of the experimental data, the value of 2.91 ± 0.28 D obtained at the PBE0+TS-vdW(SC) (330 K) level of theory is very close to the mean experimental value and is likely to be the most accurate estimate of this quantity among the series of XC functionals considered herein.

To further investigate the variation in the molecular dipole moment in liquid water as a function of the underlying XC potential, we also computed the distribution of distances between the oxygen atoms in a given water molecule and the centers of the corresponding MLWFs (*R*_{O-MLWF} = |**R**_{O} − **R**_{W}|) for each AIMD simulation performed in this work (see Fig. 7(b)). In a given water molecule, the corresponding set of four associated MLWFs can straightforwardly be identified as (i) two bonding electron pairs, centered along the O–H covalent bonds with *R*_{O-MLWF} ≈ 0.5 Å from the central oxygen atom and (ii) two lone electron pairs, oriented in an essentially tetrahedral fashion with respect to the bonding electron pairs, but centered at significantly shorter distances from the central oxygen atom (*R*_{O-MLWF} = 0.30–0.35 Å). From Fig. 7(b), it is clear that the distribution of bonding pair O-MLWF distances is essentially independent of the underlying XC potential, whereas the peak positions of the lone pair O-MLWF distance distributions systematically shift toward shorter distances (by approximately 0.01–0.02 Å) as Exx, vdW, and approximate NQE are accounted for in the AIMD simulations. This net reduction in the lone pair O-MLWF distance is indicative of a decrease in the magnitude of the local molecular dipole moment and therefore a weakened hydrogen bond network in liquid water, since a given water molecule accepts a hydrogen bond via the electrostatic interaction between one of its lone electron pairs (represented here by the MLWF) and a hydrogen atom situated on the corresponding donor water molecule. In fact, this correlation is immediately visible from the plot of the variation in the magnitude of the molecular dipole moments with the number of intact hydrogen bonds (*N*_{HB}) per water molecule as a function of the underlying XC potential in Fig. 7(c). From this data, one immediately notices an appreciably direct (and nearly linear) correlation between the number of intact hydrogen bonds per water molecule and the strength of the molecular dipole moment. Hence the systematic reduction in the strength and extent of the hydrogen bonding network as a function of systematically improving the underlying XC potential as discussed above is intimately related to the reduction in the magnitude of the local molecular dipole moment in liquid water. This trend was observed with all of the functionals considered herein and is consistent with the systematic decrease in the population of water molecules having four intact hydrogen bonds (and the corresponding increase in the population of water molecules with two and three intact hydrogen bonds) as illustrated in Sec. IV E and previous studies on gas-phase water clusters.^{110,115,116}

One can take this analysis one step further by investigating the correlation between the small shoulder situated at ≈0.31 Å in the lone pair O-MLWF distance distribution—a feature which becomes increasingly more pronounced with the collective inclusion of Exx and vdW in the underlying XC potential (see Fig. 7(b))—and the number of intact acceptor-type hydrogen bonds per water molecule. In Fig. 7(d), we performed such an analysis by assigning^{117} the intact acceptor-type hydrogen bonds to either of the two MLWFs associated with the lone electron pairs of the corresponding acceptor water molecule for the AIMD simulation at the PBE0+TS-vdW(SC) (330 K) level of theory. In this regard, we observed that the number of acceptor-type hydrogen bonds decreases as the lone pair O-MLWF distance becomes shorter, a result that is again consistent with our above findings and indicative of a strong correlation between the underlying electronic structure of a given water molecule and the overall strength and integrity of the hydrogen bond network in liquid water.

## V. CONCLUSIONS

In this work, we have performed a series of DFT-based AIMD simulations of ambient liquid water using a hierarchy of XC functionals to investigate the individual and collective effects of Exx, non-local vdW/dispersion interactions, and an approximate treatment of nuclear quantum effects on the microscopic structure of liquid water. Utilizing highly efficient linear scaling algorithms as developed and extensively optimized in our group,^{74} we have employed the PBE0 hybrid functional to account for Exx effects in conjunction with a fully self-consistent implementation of the charge density-dependent dispersion correction of Tkatchenko and Scheffler^{73} to account for vdW interactions. Based on these AIMD simulations, we found that the inclusion of Exx and vdW systematically reduces the overstructuring in structural quantities such as the oxygen-oxygen (*g*_{OO}(*r*)) and oxygen-hydrogen (*g*_{OH}(*r*)) radial distribution functions that is commonly observed in liquid water generated by GGA-DFT based AIMD simulations. Moreover, the collective effects of Exx and vdW, as resulting from a large-scale AIMD simulation of (H_{2}O)_{128} at the PBE0+TS-vdW(SC) level of theory (coupled with a 30 K increase in the simulation temperature to mimic the NQE) yield an oxygen-oxygen structure factor, *S*_{OO}(*Q*), and corresponding *g*_{OO}(*r*) that are in quantitative agreement with the best available experimental data. The molecular conformations obtained from the AIMD simulations performed in this work depend on the microscopic interactions and the sampling methodology utilized herein. Thus, the details of the simulated structures may change reflecting improvements in the underlying XC functional, stricter criteria for adiabatic separation of the nuclear and electronic coordinates, as well as an explicit quantum mechanical treatment of the nuclear degrees of freedom. While these factors may affect the outcome of future simulations of liquid water, the qualitative effects found here should be robust and include: an increase in the relative population of water molecules in the interstitial region between the first and second coordination shells, a collective reorganization in the liquid phase which is facilitated by a weakening of the hydrogen bond strength by the use of the PBE0 hybrid XC functional, coupled with a relative stabilization of the resultant disordered liquid water configurations (i.e., configurations which significantly deviate from the perfect tetrahedral network in crystalline ice) by the inclusion of non-local vdW/dispersion interactions. This increasingly more accurate description of the underlying hydrogen bond network in liquid water obtained at the PBE0+TS-vdW(SC) (330 K) level of theory also yields higher-order correlation functions, such as the oxygen-oxygen-oxygen triplet angular distribution, *P*_{OOO}(θ), and therefore the degree of local tetrahedrality, as well as electrostatic properties, such as the effective molecular dipole moment, that are in much better agreement with experiment. As such, future research directions along this line include an investigation of the individual and collective effects of Exx, vdW, and NQE on the local environment and equilibrium density of ambient liquid water—additional structural properties that should be heavily influenced by an improved underlying description of the microscopic structure of liquid water. Furthermore, the accurate description of the underlying hydrogen bond network in liquid water—as provided by the theoretical methodologies outlined in this work—provides a firm basis for spectroscopic studies of liquid water (such as X-ray absorption and photo-electron spectroscopy),^{118,119} as well as the study of solvation structure in aqueous ionic solutions, which play a key role in biology and energy research.

Although the simulations performed in this work provide us with a more accurate description of the microscopic structure of liquid water, an explicit treatment of the NQE (i.e., via the Feynman discretized PI approach) is still lacking in our current approach and will be required to quantitatively describe certain structural properties that directly involve the lighter hydrogen atoms, which can significantly deviate from classical behavior even at room temperature. The explicit inclusion of NQE via PI-AIMD simulations coupled with an underlying XC functional that accounts for both Exx and vdW interactions provides an accurate and balanced theoretical treatment of both the nuclear and electronic degrees of freedom in liquid water; as such, this is the focus of intense ongoing research in our group and will be addressed in future work.

Additional future work should also be devoted to addressing remaining deficiencies in the underlying XC functional by further reducing the deleterious effects of electron self-interaction. In addition, the theoretical description of the vdW/dispersion interactions could also be improved via the inclusion of beyond-pairwise interactions, e.g., as provided by the many-body dispersion (MBD) framework.^{120–124} Moving beyond the realm of DFT, the theoretical description of the electronic degrees of freedom could also be addressed using more accurate quantum chemistry approaches and/or quantum Monte Carlo; however, the high computational cost associated with these methodologies has restricted their use to date in large-scale AIMD simulations of condensed-phase systems.

## ACKNOWLEDGMENTS

R.A.D., B.S., and R.C. acknowledge support from the Department of Energy (DOE) under Grant Nos. DE-SC0005180 and DE-SC000826. R.A.D., Z.L., and R.C. also acknowledge support from the National Science Foundation (NSF) under Grant No. CHE-0956500 during the early stages of this work. X.W. acknowledges support from the American Chemical Society Petroleum Research Fund (ACS PRF) under Grant No. 53482-DNI6 and the DOE under Grant No. DE-SC0008726. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This research used resources of the Argonne Leadership Computing Facility at Argonne National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-06CH11357. Additional computational resources were provided by the Terascale Infrastructure for Groundbreaking Research in Science and Engineering (TIGRESS) High Performance Computing Center and Visualization Laboratory at Princeton University.

## REFERENCES

*n*

_{OO}is very sensitive to

*g*

_{OO}(

*r*), since

*n*

_{OO}is computed as

*W*

_{c}) and the donating hydrogen atom is less than 2.7 Å and the angle ∠W

_{c}− O⋅⋅⋅H is less than 30°. For a given hydrogen bond, if both of the MLWFs satisfy these two geometric criteria (either partially or fully), then the two MLWFs share that particular hydrogen bond with a distance-weighted fraction.