The Individual and Collective Effects of Exact Exchange and Dispersion Interactions on the Ab Initio Structure of Liquid Water

The Individual and Collective Effects of Exact Exchange and Dispersion Interactions on the Ab Initio Structure of Liquid Water

Abstract

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 vdW/dispersion interactions, via a fully self-consistent density-dependent dispersion correction, and approximate nuclear quantum effects (aNQE), 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, vdW, and aNQE as resulting from a large-scale AIMD simulation of (HO) at the PBE0+vdW level of theory, significantly softens the structure of ambient liquid water and yields an oxygen-oxygen structure factor, , and corresponding oxygen-oxygen radial distribution function, , that are now in quantitative agreement with the best available experimental data. This level of agreement between simulation and experiment as 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 the PBE0 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 obtained at the PBE0+vdW+aNQE level of theory yields higher-order correlation functions, such as the oxygen-oxygen-oxygen triplet angular distribution, , and therefore the degree of local tetrahedrality, as well as electrostatic properties, such as the effective molecular dipole moment, that are also 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.). Franks (2000) While there is no experimental methodology currently available to directly obtain the real-space microscopic structure of liquid water, Soper (2013); Head-Gordon and Hura (2002) ab initio molecular dynamics (AIMD) Marx and Hutter (2009) simulations employing Density Functional Theory (DFT) bur () 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 Car and Parrinello (1985) 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, Laasonen et al. (1993); Tuckerman et al. (1995) AIMD has been applied to many complex problems in biology, chemistry, and energy research, e.g., the designing of efficient catalysts for hydrogen production Zipoli et al. (2010) and the modeling of the autoionization Geissler et al. (2001) and electrocatalytic splitting of water, Ikeshoji et al. (2011) 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 Asthagiri et al. (2003); Grossman et al. (2004); Schwegler et al. (2004); Fernández-Serra and Artacho (2004); Kuo et al. (2004); McGrath et al. (2005); VandeVondele et al. (2005); Sit and Marzari (2005); Fernández-Serra et al. (2005); McGrath et al. (2006); Lee and Tuckerman (2006); Todorova et al. (2006); Lee and Tuckerman (2007); Guidon et al. (2008); Kühne et al. (2009); Mattson and Mattson (2009); Yoo et al. (2009); Zhang et al. (2011a); Bartók et al. (2013); Alfè et al. (2013); Lin et al. (2009); Jonchiere et al. (2011); Wang et al. (2011); Zhang et al. (2011b); Møgelhøj et al. (2011); Lin et al. (2012); Yoo and Xantheas (2012); Schmidt et al. (2009); Ma et al. (2012); Ben et al. (2013) as well as crystalline phases of ice. Feibelman (2008); Santra et al. (2011, 2013); Labat et al. (2011); Kambara et al. (2012); Murray and Galli (2012); Fang et al. (2013); Macher et al. (2014) For one, GGA-DFT suffers from the well-known self-interaction error Perdew and Zunger (1981) which leads to excessive proton delocalization in liquid water, yielding red shifts in the predicted O-H stretching frequencies as compared to experiment. Zhang et al. (2011b, a) Secondly, 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. Schmidt et al. (2009); Wang et al. (2011); Ma et al. (2012) 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, Chen et al. (2003); Morrone and Car (2008); Ceriotti et al. (2013) and this leads to overstructuring in the predicted radial distribution functions. Morrone and Car (2008); Ceriotti and Manolopoulos (2012); Paesani et al. (2007); Fanourgakis et al. (2006)

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, Santra et al. (2007, 2008, 2009); Wang et al. (2010); Gillan et al. (2012) although recently hybrid functionals have been applied in the study of several crystalline phases of ice. Santra et al. (2011, 2013); Labat et al. (2011); Kambara et al. (2012); Erba et al. (2009) 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. Santra et al. (2009); Wang et al. (2010); Gillan et al. (2012); Santra et al. (2011); Labat et al. (2011); Kambara et al. (2012); Erba et al. (2009); Santra et al. (2013) Although applications of hybrid functionals to liquid water have been relatively rare in the literature, it was found that PBE0 Perdew et al. (1996a); Adamo and Barone (1999) greatly improves upon the accuracy of the vibrational spectrum of liquid water. Zhang et al. (2011b, a) 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, Zhang et al. (2011b, a); Todorova et al. (2006) others have found the effects of Exx to be negligible in this regard. Guidon et al. (2008); Ben et al. (2013) 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 manuscript.

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 Santra et al. (2011); Kambara et al. (2012); Murray and Galli (2012) and the predicted equilibrium density of liquid water. Schmidt et al. (2009); Ma et al. (2012); Ben et al. (2013) 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. Lin et al. (2009); Jonchiere et al. (2011); Wang et al. (2011); Zhang et al. (2011b); Møgelhøj et al. (2011); Lin et al. (2012); Yoo and Xantheas (2012); Schmidt et al. (2009); Ma et al. (2012); Ben et al. (2013) For instance, the pairwise-additive based vdW correction approaches of Grimme, Grimme (2004); Grimme et al. (2010) when used in conjunction with the popular PBE Perdew et al. (1996b) and BLYP Becke (1988); Lee et al. (1988) GGA functionals, tend to reduce the intensity of the first maximum in the oxygen-oxygen radial distribution function () of liquid water over a fairly wide range, i.e., by approximately 5% for PBE-D Schmidt et al. (2009) and 10-17% for BLYP-D. Lin et al. (2009); Schmidt et al. (2009); Jonchiere et al. (2011) Other popular vdW-inclusive DFT approaches include the so-called vdW-DF functionals, Dion et al. (2004) 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. Wang et al. (2011); Møgelhøj et al. (2011) 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 Perdew et al. (1996b) (PBE), the corresponding PBE0 Perdew et al. (1996a); Adamo and Barone (1999) hybrid, and self-consistent dispersion-corrected analogs dis () of each, i.e., PBE+vdW and PBE0+vdW, based on the Tkatchenko-Scheffler Tkatchenko and Scheffler (2009) density-dependent vdW/dispersion functional (TS-vdW). To successfully perform these AIMD simulations, we have employed a linear scaling O() algorithm to compute the exact exchange energy, Wu et al. (2009) which allows for relatively long ( ps) hybrid functional based simulations using large system sizes (e.g., the explicit treatment of 128 water molecules). In addition, we have also developed and utilized a linear scaling O() self-consistent implementation of the TS-vdW dispersion correction, which provides a framework for computing atomic dispersion coefficients as an explicit functional of the charge density, thereby accounting for the local chemical environment surrounding each atom. Tkatchenko and Scheffler (2009) In this manuscript, 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 which has been shown to mimic the NQE in the Morrone and Car (2008); Paesani et al. (2007) By approximately accounting for NQE in AIMD simulations at the PBE0+vdW level of theory, which is the most accurate XC functional considered in this work, we have obtained quantitative agreement with scattering experiments in the 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 manuscript 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. Sec. 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, Perdew et al. (1996b) the corresponding hybrid PBE0 Perdew et al. (1996a); Adamo and Barone (1999) which includes 25% exact exchange, and the self-consistent dispersion-corrected analogs dis () thereof, i.e., PBE+vdW and PBE0+vdW, based on the Tkatchenko-Scheffler Tkatchenko and Scheffler (2009) density-dependent vdW/dispersion functional. 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 () 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) Car and Parrinello (1985) 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. Martyna et al. (1992) 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. Tobias et al. (1993) The electronic wavefunctions were expanded using a plane wave basis set with a kinetic energy cutoff of 71 Ry. 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. Troullier and Martins (1991) 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, Grossman et al. (2004) 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. Tassone et al. (1994)

Using a simulation cell with Å, four AIMD simulations were performed on (HO) at 300 K using the PBE, PBE0, PBE+vdW, and PBE0+vdW 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 (HO) with a temperature of 330 K at the PBE0+vdW level, which is the most accurate XC functional considered herein. In this regard, an increase of approximately 30 K in the simulation temperature has been shown to mimic the nuclear quantum effects (NQE) in structural quantities such as the oxygen-oxygen radial distribution function () in both DFT Morrone and Car (2008) and force-field Paesani et al. (2007); Fanourgakis et al. (2006) based MD simulations.

All calculations reported in this work utilized a modified development version of the Quantum ESPRESSO (QE) software package. Giannozzi et al. (2009) 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., (HO) at the PBE0+vdW 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

iii.1 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() exact exchange algorithm, which has been developed Wu et al. (2009) and extensively optimized in our research group. The computational savings afforded by this algorithm originate from the maximally localized Wannier function (MLWF) Marzari and Vanderbilt (1997) representation of the occupied subspace, , obtained via a unitary transformation of the occupied Kohn-Sham electronic states, , i.e., , and the subsequent efficient use of sparsity in the numerical evaluation of all required quantities in real-space.

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

(1)

and orbital-dependent “force” acting on the electronic wavefunctions in the CP dynamics, Car and Parrinello (1985)

(2)

can both be written in terms of , the electrostatic potential generated by the pair density, , i.e.,

(3)

The potential is the solution of the Poisson equation,

(4)

subject to the boundary condition that at infinity. By virtue of the localized character of , the potential at sufficient distances from the center of the pair density is exactly represented by a multipolar expansion,

(5)

in which are the multipoles of the charge distribution . The integrals required to evaluate the multipoles are restricted to a volume that fully contains , as delimited by the Poisson equation cutoff radius, i.e., for . Eq. (4) is then solved numerically in real-space (i.e., on the dense real-space grid) within the sphere defined by , with the boundary condition provided by Eq. (5) for .

To efficiently compute the exact exchange energy at every AIMD step, the required MLWFs were evaluated by minimizing the total spread functional, , using second-order damped CP dynamics as described in Ref. [Sharma et al., 2003]. During this procedure, we used a fictitious mass of 500 a.u., a damping coefficient of 0.3, and a time step of 4 a.u. On average, less than 5 iterations per AIMD step were sufficient to reach a convergence of 10 in the optimal spread . The resulting average spread for a single MLWF was 0.52 Å during the PBE0 and PBE0+vdW 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 above was defined with respect to the center of a given pair density and set to a value of 2.65 Å, to ensure an accurate computation of the exact exchange energy. The maximum angular momentum component retained in the multipolar expansion given in Eq. (5) was , which has previously been found to be sufficiently accurate in computing exact exchange energies in a variety of systems. Wu et al. (2009) This procedure provides the exact exchange energy with an accuracy of at least Hartrees/molecule, as was tested by systematically increasing the cutoff radius and including pairs of MLWFs up to distances of half the simulation cell length.

iii.2 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 O() implementation dis () of the density-dependent dispersion correction of Tkatchenko and Scheffler Tkatchenko and Scheffler (2009) (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 and as

(6)

in which 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 and , , 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 and , , and is the scalar distance between atoms and , . 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 Dobson (2014) as e.g., the DFT+D approaches of Grimme. Grimme (2004); Grimme et al. (2010)

The implementation utilized herein in conjunction with the PBE and PBE0 XC functionals, namely PBE+vdW and PBE0+vdW, respectively, allows us to self-consistently account for the changes in the total energy and charge density () 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, Car and Parrinello (1985)

(7)

in which is the local vdW/dispersion potential derived via as

(8)

which is added to the local XC potential, , at each MD step.

During the AIMD simulations at the PBE+vdW and PBE0+vdW levels of theory, the scaling parameter in the damping function above, which determines the onset of the vdW correction for a given XC functional, was set to 0.94 and 0.96, respectively, as recommended in Ref. [Tkatchenko and Scheffler, 2009]. The free atom reference quantities utilized herein were bohr, Hartree bohr, and bohr for oxygen and bohr, Hartree bohr, and bohr for hydrogen. The 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 to at least Hartrees/molecule.

Iv Results and Discussion

iv.1 Comparative Analysis of the Oxygen-Oxygen Structure Factor

We begin our analysis by comparing the oxygen-oxygen structure factor, , obtained from the highest quality AIMD simulations performed herein at the PBE0+vdW (330 K) level of theory and various X-ray scattering measurements (See Fig. 1(a)).

Figure 1: Panel (a): The oxygen-oxygen structure factors, , of liquid water obtained from theory (via the DFT-based AIMD simulations performed in this work) and various X-ray scattering experiments. Sorenson et al. (2000); Huang et al. (2011); Skinner et al. (2013) Panel (b): The oxygen-oxygen radial distribution functions, , of liquid water obtained from theory (via the DFT-based AIMD simulations performed in this work) and various scattering experiments. Skinner et al. (2013); Soper and Benmore (2008)

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 to be accurately obtained from the directly measured scattering intensity, , and knowledge of the corresponding form factor, , via , in which represents the change in the wavevector between the incident and scattered radiation. Once has been experimentally determined, the so-called pair correlation or radial distribution function, , 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.,

(9)

wherein is the maximum value of in the given scattering experiment and is the atomic number density. Hence, the process of obtaining an accurate is highly sensitive to the 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, , required for deriving from the experimentally measured scattering intensity.

As shown in Fig. 1(a), the taken from various X-ray scattering experiments Sorenson et al. (2000); Huang et al. (2011); Skinner et al. (2013) with increasing values of are very similar across the range of values accessible in a given scattering experiment; despite this fact, each of these X-ray scattering experiments yields significantly different , especially with respect to the intensity and position of the first peak. Sorenson et al. (2000); Huang et al. (2011); Skinner et al. (2013) In this regard, we note that these differences in the experimentally derived are almost entirely due to the aforementioned limitation inherent to the X-ray scattering experiments and not due to errors in the corresponding . To justify this claim, we first computed the difference between the 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, Wang et al. (1994); Watanabe et al. (1999) 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.Krack et al. (2002) we then computed the difference between the DFT-based for the single gas-phase water molecule and the corresponding to condensed-phase liquid water obtained from the DFT-based AIMD simulations performed herein. zha () 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- region (i.e., for Å); hence, such differences would be essentially negligible when computing the corresponding .

Using an AIMD trajectory obtained at the PBE0+vdW (330 K) level of theory, we have computed as:

(10)

in which is the number of oxygen atoms in the simulation cell, is the degeneracy of the wavevector shell of length , and is the real-space position vector of -th oxygen atom. Unlike the experimentally derived , the accuracy of the computed in the low- region is subject to a large degree of numerical noise, a consequence of the fact that there is a relatively small number of 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+vdW (330 K) level) using a periodic simple cubic box containing 128 water molecules ( Å), wherein the smallest accessible value is 0.4 Å, which helps in improving the accuracy of the theoretically determined in the low- region.

From Fig. 1(a), it is clear that the theoretically determined is in nearly quantitative agreement with the experimental results across the entire region accessible to the X-ray scattering experiments, with only a slightly noticeable shift towards higher values. This level of agreement between theory and experiment has been achieved via the utilization of the PBE0+vdW 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 (Exx), (ii) accounts for non-local vdW/dispersion interactions by a self-consistent density-dependent correction, and (iii) approximately corrects for nuclear quantum effects (aNQE) 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 , 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.

iv.2 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, , directly obtained from X-ray scattering experiments can be systematically improved upon via the use of larger values, which has yielded essentially converged results (i.e., to within 1%) for the intensity and position of the first peak in the of liquid water with ÅSkinner et al. (2013) In Fig. 1(b), this latest X-ray scattering data, Skinner et al. (2013) which is now in quantitative agreement with the obtained from empirical potential structural refinement (EPSR) based on joint X-ray/neutron data by Soper and Benmore, Soper and Benmore (2008) are both compared against the highest quality AIMD simulations performed herein at the PBE0+vdW (330 K) level of theory. From this figure, one immediately notices that the obtained from PBE0+vdW (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 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, vdW, and aNQE on the predicted of liquid water by considering Fig. 2, which provides the results obtained from AIMD simulations based on each of the different XC functionals utilized in this work and Table 1, which summarizes the positions and intensities at various key points in each , the oxygen-oxygen coordination numbers, , associated with the first shell in liquid water, as well as other structural properties which will be discussed in detail in Sec. IV.5 - IV.6.

Figure 2: The oxygen-oxygen radial distribution functions of liquid water obtained from theory (via the DFT-based AIMD simulations performed in this work). The arrows indicate systematic shifts in the main peak positions and intensities of the computed distributions with improvement of the underlying exchange-correlation functional.
Method
PBE 300 2.69 3.28 3.28 0.37 4.44 1.44 4.03 3.83 0.78 3.19
PBE+vdW 300 2.71 2.99 3.27 0.54 4.40 1.33 4.01 3.74 0.74 3.12
PBE0 300 2.71 2.96 3.30 0.53 4.39 1.35 4.10 3.71 0.73 3.00
PBE0+vdW 300 2.72 2.76 3.31 0.70 4.47 1.20 4.22 3.62 0.69 2.96
PBE0+vdW 330 2.74 2.51 3.33 0.84 4.44 1.22 4.32 3.48 0.64 2.91
Expt. 1 Skinner et al. (2013) 295 2.80 2.57 3.45 0.84 4.50 1.12 4.3(1) - - -
Expt. 2 Soper and Benmore (2008) 298 2.76 2.62 3.42 0.82 4.43 1.14 4.67(5) - 0.58 -
Table 1: Tabulated summary of the structural properties of liquid water obtained via the DFT-based AIMD simulations performed in this work and various scattering experiments. From left to right, the specific XC functional and temperature ( in K) utilized in a given AIMD simulation; positions (in Å) and intensities of the first maximum ( and ), the first minimum ( and ), and the second maximum ( and ) of the corresponding oxygen-oxygen radial distribution functions (See Sec. IV.2); oxygen-oxygen coordination numbers () associated with the first shell, computed by integrating each respective up to the first minimum (See Sec. IV.2); the average tetrahedrality parameter () per water molecule (See Sec. IV.4); the average number of hydrogen bonds () per water molecule (See Sec. IV.5); the average local dipole moment ( in Debye) per water molecule (See Sec. IV.6). For reference, the available corresponding experimental data is provided in the last two rows.

The first point to note here is the fact that the 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, is typically found to be greater than 3.2, Grossman et al. (2004); Schwegler et al. (2004); Kühne et al. (2009) followed by a first minimum that is too deep with respect to experiment. This level of overstructuring in the 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 is observed in which the intensity of the first maximum, , is reduced by approximately 0.32 and the intensity of the first minimum, , is increased by approximately 0.16, coupled with an increase in the positions of the first maximum, , and first minimum, , by 0.02 Å  with respect to the generated at the PBE level of theory (See Fig. 2 and Table 1). These noticeable differences in the 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. Zhang et al. (2011a) In fact, this effect has also been demonstrated in previous studies on both gas-phase water clusters Santra et al. (2009); Gillan et al. (2012); Wang et al. (2010) and ice, Santra et al. (2011); Labat et al. (2011); Kambara et al. (2012) 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 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 when exact exchange was accounted for over the standard PBE functional. Guidon et al. (2008); Ben et al. (2013) 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 Wu et al. (2009) (See Sec. III), and that of Guidon et al.Guidon et al. (2008) in which the Coulomb potential in the exact exchange integrals was replaced by a finite-range approximant in terms of the complementary error function, .

By including a self-consistent treatment of the non-local vdW/dispersion interactions via the density-dependent Tkatchenko-Scheffler dispersion correction, Tkatchenko and Scheffler (2009) we observed a very similar effect on the of liquid water; in fact, the radial distribution functions obtained at the PBE0 and PBE+vdW levels of theory are almost identical over the entire distance range considered (See Fig. 2 and Table 1). 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 to move inward and populate the interstitial region. As a result, the hydrogen bonding (See Sec. IV.5) 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 as that observed in PBE0-based liquid water. In this regard, there is only a subtle difference in the of liquid water between the PBE+vdW and PBE0 levels of theory: for PBE+vdW, the position of the first minimum, is slightly decreased with respect to PBE by approximately 0.01 Å  as opposed to the increase of 0.02 Å  observed above for PBE0. Due to this decrease in at the PBE+vdW level of theory, the integrated number of water molecules in the first coordination shell, coo () , actually decreases slightly when compared to PBE (4.01 vs. 4.03), whereas at the PBE0 level, increases to a value of 4.10 since increases (See Table 1).

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. Lin et al. (2009); Jonchiere et al. (2011); Wang et al. (2011); Zhang et al. (2011b); Møgelhøj et al. (2011); Lin et al. (2012); Yoo and Xantheas (2012); Schmidt et al. (2009); Ma et al. (2012); Ben et al. (2013). For instance, the most commonly utilized approach involves the pairwise-additive vdW/dispersion correction of Grimme Grimme (2004); Grimme et al. (2010) in conjunction with the popular PBE and BLYP GGA functionals, which tends to reduce the intensity of the first maximum in the of liquid water over a fairly wide range, i.e., by approximately 5% for PBE-D Schmidt et al. (2009) and 10-17% for BLYP-D. Lin et al. (2009); Schmidt et al. (2009); Jonchiere et al. (2011) Other popular vdW-inclusive DFT approaches for studying liquid water include the so-called vdW-DF functionals, Dion et al. (2004) 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 , especially when employed in conjunction with the revPBE exchange functional; Wang et al. (2011); Møgelhøj et al. (2011) the use of the PBE exchange functional instead recovers this limitation to some extent Wang et al. (2011); Zhang et al. (2011b); Corsetti et al. (2013) and reduces the intensity of the first maximum by 25% with respect to PBE. Zhang et al. (2011b) For comparison, the self-consistent implementation of the density-dependent Tkatchenko-Scheffler Tkatchenko and Scheffler (2009) vdW/dispersion correction employed in this work, PBE+vdW, 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+vdW level of theory, yields a 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, , by 0.52 and to increase the intensity of the first minimum, , by 0.33, both of which represent larger effects than observed when the PBE0 and PBE+vdW XC functionals were utilized independently. The collective reduction here in the first maximum of 0.52 at the PBE0+vdW level of theory is slightly less than the sum of the individual reductions observed when utilizing PBE0 (0.32) and PBE+vdW (0.29) independently, which is indicative of the fact that there is a compensating effect present here in the short-range (i.e., at distances which include hydrogen-bonded oxygens) when utilizing a self-consistent vdW/dispersion correction scheme that depends on the underlying charge density (i.e., the vdW/dispersion correction is a functional of the electron density, ). In the same breath, the collective increase in the first minimum of 0.33 at the PBE0+vdW level of theory is exactly equal to the sum of the individual increases observed when utilizing PBE0 (0.16) and PBE+vdW (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+vdW : 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+vdW), while the increase in the position of the first minimum of 0.03 Å  represents an interesting collective effect between PBE0 (which increased by 0.02 Å) and PBE+vdW (which decreased by 0.01 Å). Furthermore, this outward shift in at the PBE0+vdW level, while in the right direction, is still not large enough to exactly match the experimental findings, which is suggestive that a further increase in the amount of exact exchange may still be needed in the underlying XC potential.

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.5, another interesting point can be made here concerning the contributions arising from select neighboring water molecules to the overall as a function of the underlying XC potential.

Figure 3: Contributions from the fourth (4th) and fifth (5th) neighbors (neighboring water molecules) to the oxygen-oxygen radial distribution functions of liquid water obtained from theory (via the DFT-based AIMD simulations performed in this work). The arrows indicate systematic shifts in the main peak positions and intensities of the computed distributions with improvement of the underlying exchange-correlation functional.

In Fig. 3, the contributions to the 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 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 ; 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 . 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 of liquid water by performing the simulation at the PBE0+vdW level of theory with a 30 K increase in the simulation temperature (See Fig. 2 and Table 1). At 330 K, we observe further softening in the of liquid water, with peak positions and intensities that are now in quantitative agreement with the directly obtained from X-ray Skinner et al. (2013) 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. Skinner et al. (2013) and only 2.4% higher than the value provided by Soper and Benmore. Soper and Benmore (2008) 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 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.

iv.3 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, , 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 correction), and approximate nuclear quantum effects (via a 30 K increase in the simulation temperature).

Figure 4: The oxygen-hydrogen radial distribution functions of liquid water obtained from theory (via the DFT-based AIMD simulations performed in this work) and scattering experiments. Soper and Benmore (2008) The arrow indicates systematic shifts in the main peak positions and intensities of the computed distributions with improvement of the underlying exchange-correlation functional.

In this case, however, the width and intensity of the first peak of the 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. Morrone and Car (2008); Ceriotti et al. (2013) 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 , 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 of liquid water, in that the AIMD simulations at the PBE0 and PBE+vdW levels of theory predict nearly identical over the entire distance range considered (See Fig. 4). Considering for instance the intensity of the second maximum in the , there is an overall reduction or destructuring of approximately 0.22 (PBE0) and 0.15 (PBE+vdW), 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 is only slightly larger than the observed effect resulting from an explicit treatment of vdW/dispersion interactions. When used in combination, the PBE0+vdW XC functional again displays a nearly additive effect on the 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 is provided by PBE0+vdW 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. Soper and Benmore (2008) We again stress here that structural quantities that directly involve the hydrogen atoms in liquid water, such as the and , 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 , which was considered in detail in Sec. IV.2, 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.

iv.4 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, , within the first coordination shell (See Fig. 5).

Figure 5: The oxygen-oxygen-oxygen triplet angular distribution functions of liquid water obtained from theory (via the DFT-based AIMD simulations performed in this work) and empirical potential structural refinement (EPSR) based on joint X-ray/neutron scattering data. Soper and Benmore (2008) The triplet angular distribution functions show here were normalized to . The arrow indicates systematic shifts in the main peak positions and intensities of the computed distributions with improvement of the underlying exchange-correlation functional.

In general, three-body correlation functions such as are not directly accessible from scattering experiments; however, an experimental distribution of was recently extracted via empirical potential structural refinement (EPSR) based on joint X-ray/neutron scattering data. Soper and Benmore (2008) 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 Skinner et al. (2013) and EPSR based on joint X-ray/neutron scattering data Soper and Benmore (2008) supports the effectiveness of this procedure. To extract the corresponding 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. [Soper and Benmore, 2008], 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 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 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, , as Errington and Debenedetti (2001)

(11)

in which is the angle formed by a central given water molecule and its nearest neighbors and , we find that PBE yields an average value of 0.78, which is indeed much higher than the estimated EPSR-experimental value of  Soper and Benmore (2008) (See Table 1). For reference, a system with perfect tetrahedral order would yield , whereas 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 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 Jonchiere et al. (2011); Zhang et al. (2011c) with values of , which is in reasonable agreement with the value of computed herein at the PBE+vdW 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 Møgelhøj et al. (2011) At the PBE0 level, we computed a value of for the tetrahedral order parameter, which is extremely close to the value obtained above in the PBE+vdW case and in good agreement with the previous study of Zhang et al. Zhang et al. (2011c) When used in conjunction, the self-consistent PBE0+vdW exchange-correlation functional yields a value of 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 Sec. IV.1IV.3, even further improvement comes with the approximate inclusion of NQE, wherein we computed a value of , which is now approximately 10.3% larger than the EPSR-experimental value. Interestingly, the shoulder present in the EPSR-experimental distribution at approximately , which is indicative of a highly distorted hydrogen bond network in the first coordination shell, only becomes noticeably visible at the PBE0+vdW (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.

iv.5 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. Eaves et al. (2005) 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. Prada-Gracia et al. (2013) 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, Luzar and Chandler (1996) wherein the defining parameters for an intact hydrogen bond include both a distance and angular criterion, namely, Å  and , with .

Using this definition, we computed the average number of intact hydrogen bonds per water molecule hbo () for each of the AIMD simulation performed in this work and reported these numbers in Table 1. 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+vdW (3.74), PBE0 (3.71), PBE0+vdW (3.67), and PBE0+vdW+aNQE (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 in Fig. 2.

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., , in Fig. 6.

Figure 6: Percentage-wise decomposition into acceptor- (A) and donor- (D) types of the intact hydrogen bonds per water molecule in liquid water obtained from theory (via the DFT-based AIMD simulations performed in this work). The -axis labels indicate the number of acceptor-type () and donor-type () hydrogen bonds, respectively. All other acceptor–donor combinations, e.g., , , , etc., had contributions of % and were therefore omitted from this figure for clarity. The inset depicts donor (O) and acceptor (O) oxygen atoms participating in a hydrogen bond via the hydrogen (H) atom.

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 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 type hydrogen bonds is a majority at 75.4%, a quantity which is reduced to just less than 50% when Exx, vdW, and aNQE effects are accounted for in the underlying XC potential at the PBE0+vdW (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 and 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+vdW (330 K); in this case, the percentage of the and 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 type, or more than four hydrogen bonds, as in the 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 hydrogen bond type approximately increased by a factor of four from 2.4% (PBE) to 9.3% (PBE0+vdW (330 K)), while the percentage of the 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 aNQE in the microscopic structure of liquid water; with Exx and aNQE 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.

iv.6 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 Debye (D) from spectroscopic measurements, Shostak et al. (1991); Clough et al. (1973) and several DFT functionals can reproduce this experimental value to within 3%. Zhang et al. (2011b) 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 ( D), Badyal et al. (2000) 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–3.1 D depending on the partitioning scheme employed. Batista et al. (1999) 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. Silvestrelli and Parrinello (1999a, b); Sharma et al. (2005); Zhang et al. (2011b); Sharma et al. (2007) For the case of liquid water, the MLWFs associated with different molecules have only a very minor overlap, Silvestrelli and Parrinello (1999a) 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. Sharma et al. (2005, 2007)

Figure 7: Panel (a): Probability density plots of the distributions of molecular dipole moment magnitudes ( in Debye) in liquid water obtained from theory (via the DFT-based AIMD simulations performed in this work). Panel (b): Probability density plots of the distributions of distances between the oxygen atoms in a given water molecule and the centers of the corresponding maximally localized Wannier functions ( in Å). Panel (c): Plots of the variation in the magnitude of the molecular dipole moments ( in Debye) with the number of intact hydrogen bonds () per water molecule. Panel (d): Probability density plot of the O–MLWF distance distribution (red line, left axis) reproduced from Panel (b) and the average number of acceptor-type () hydrogen bonds per water molecule (black circles, right axis) as a function of in liquid water obtained from an AIMD simulation at the PBE0+vdW (330 K) level of theory.

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

(12)

in which , , and are the Cartesian coordinates of the hydrogen and oxygen atoms comprising the water molecule, respectively, and are the coordinates of the four corresponding MLWF centers.

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, vdW, and aNQE 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 ( D), PBE+vdW ( D), PBE0 ( D), PBE0+vdW ( D), and PBE0+vdW+aNQE ( 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. Zhang et al. (2011b) 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 D obtained at the PBE0+vdW (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 () 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 Å  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 ( Å). 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 Å) as Exx, vdW, and aNQE effects 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 () 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.5 and previous studies on gas-phase water clusters. Batista et al. (1999); Kemp and Gordon (2008); Tu and Laaksonen (2000)

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, vdW, and aNQE 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 mlw () 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+vdW (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 exact exchange (Exx), non-local vdW/dispersion interactions, and an approximate treatment of nuclear quantum effects (aNQE) on the microscopic structure of liquid water. Utilizing highly efficient linear scaling algorithms as developed and extensively optimized in our group, Wu et al. (2009) 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 Tkatchenko and Scheffler (2009) 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 () and oxygen-hydrogen () radial distribution functions that is commonly observed in liquid water generated by GGA-DFT based AIMD simulations. Moreover, the collective effects of Exx, vdW, and aNQE, as resulting from a large-scale AIMD simulation of (HO) at the PBE0+vdW level of theory (coupled with a 30 K increase in the simulation temperature to mimic the NQE) yield an oxygen-oxygen structure factor, , and corresponding 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+vdW+aNQE level of theory also yields higher-order correlation functions, such as the oxygen-oxygen-oxygen triplet angular distribution, , 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), Kong et al. (2012); Swartz and Wu (2013) 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 path integral (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. Tkatchenko et al. (2012); R. A. DiStasio Jr. et al. (2012); Tkatchenko et al. (2013); Ambrosetti et al. (2014); R. A. DiStasio Jr. et al. (2014) 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.

Acknowledgements.
R.D., B.S., and R.C. acknowledge support from the Department of Energy (DOE) under Grant Nos. DE-SC0005180 and DE-SC000826. R.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 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

  1. F. Franks, Water: A Matrix of Life, 2nd ed. (The Royal Society of Chemistry, Cambridge, 2000).
  2. A. K. Soper, ISRN Physical Chemistry 2013, 279463 (2013).
  3. T. Head-Gordon and G. Hura, Chem. Rev. 102, 2651 (2002).
  4. D. Marx and J. Hutter, Ab Inito Molecular Dynamics: Basic Theory and Advanced Methods (Cambridge University Press, Cambridge, 2009).
  5. For a recent review of DFT, see e.g., K. Burke, J. Chem. Phys. 136, 150901 (2012).
  6. R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
  7. K. Laasonen, M. Sprik, M. Parrinello,  and R. Car, J. Chem. Phys. 99, 9080 (1993).
  8. M. Tuckerman, K. Laasonen, M. Sprik,  and M. Parrinello, J. Chem. Phys. 103, 150 (1995).
  9. F. Zipoli, R. Car, M. H. Cohen,  and A. Selloni, J. Am. Chem. Soc. 132, 8593 (2010).
  10. P. L. Geissler, C. Dellago, D. Chandler, J. Hutter,  and M. Parrinello, Science 291, 2121 (2001).
  11. T. Ikeshoji, M. Otani, I. Hamada,  and Y. Okamoto, Phys. Chem. Chem. Phys. 13, 20223 (2011).
  12. D. Asthagiri, L. R. Pratt,  and J. D. Kress, Phys. Rev. E 68, 041505 (2003).
  13. J. C. Grossman, E. Schwegler, E. W. Draeger, F. Gygi,  and G. Galli, J. Chem. Phys. 120, 300 (2004).
  14. E. Schwegler, J. C. Grossman, F. Gygi,  and G. Galli, J. Chem. Phys. 121, 5400 (2004).
  15. M. V. Fernández-Serra and E. Artacho, J. Chem. Phys. 121, 11136 (2004).
  16. I.-F. W. Kuo, C. J. Mundy, M. J. McGrath, J. I. Siepmann, J. VandeVondele, M. Sprik, J. Hutter, B. Chen, M. L. Klein, F. Mohamed, M. Krack,  and M. Parrinello, J. Phys. Chem. B 108, 12990 (2004).
  17. 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).
  18. J. VandeVondele, F. Mohamed, M. Krack, J. Hutter, M. Sprik,  and M. Parrinello, J. Chem. Phys. 122, 014515 (2005).
  19. P. H.-L. Sit and N. Marzari, J. Chem. Phys. 122, 204510 (2005).
  20. M. V. Fernández-Serra, G. Ferlat,  and E. Artacho, Molecular Simulation 31, 361 (2005).
  21. M. J. McGrath, J. I. Siepmann, I.-F. W. Kuo,  and C. J. Mundy, Mol. Phys. 104, 3619 (2006).
  22. H.-S. Lee and M. E. Tuckerman, J. Chem. Phys. 125, 154507 (2006).
  23. T. Todorova, A. P. Seitsonen, J. Hutter, I.-F. W. Kuo,  and C. J. Mundy, J. Phys. Chem. B 110, 3685 (2006).
  24. H.-S. Lee and M. E. Tuckerman, J. Chem. Phys. 126, 164501 (2007).
  25. M. Guidon, F. Schiffmann, J. Hutter,  and J. VandeVondele, J. Chem. Phys. 128, 214104 (2008).
  26. T. D. Kühne, M. Krack,  and M. Parrinello, J. Chem. Theory Comput. 5, 235 (2009).
  27. A. E. Mattson and T. R. Mattson, J. Chem. Theory Comput. 5, 887 (2009).
  28. S. Yoo, X. C. Zeng,  and S. S. Xantheas, J. Chem. Phys. 130, 221102 (2009).
  29. C. Zhang, D. Donadio, F. Gygi,  and G. Galli, J. Chem. Theory Comput. 7, 1443 (2011a).
  30. A. P. Bartók, M. J. Gillan, F. R. Manby,  and G. Csányi, Phys. Rev. B 88, 054104 (2013).
  31. D. Alfè, A. P. Bartók, G. Csányi,  and M. J. Gillan, J. Chem. Phys. 138, 221102 (2013).
  32. I.-C. Lin, A. P. Seitsonen, M. D. Coutinho-Neto, I. Tavernelli,  and U. Rothlisberger, J. Phys. Chem. B 113, 1127 (2009).
  33. R. Jonchiere, A. P. Seitsonen, G. Ferlat, A. M. Saitta,  and R. Vuilleumier, J. Chem. Phys. 135, 154503 (2011).
  34. J. Wang, G. Román-Pérez, J. M. Soler, E. Artacho,  and M.-V. Fernández-Serra, J. Chem. Phys. 134, 024516 (2011).
  35. C. Zhang, J. Wu, G. Galli,  and F. Gygi, J. Chem. Theory Comput. 7, 3054 (2011b).
  36. A. Møgelhøj, A. K. Kelkkanen, K. T. Wikfeldt, J. Schiøtz, J. J. Mortensen, L. G. M. Pettersson, B. I. Lundqvist, K. W. Jacobsen, A. Nilsson,  and J. K. Nørskov, J. Phys. Chem. B 115, 14149 (2011).
  37. I.-C. Lin, A. P. Seitsonen, I. Tavernelli,  and U. Rothlisberger, J. Chem. Theory Comput. 8, 3902 (2012).
  38. S. Yoo and S. S. Xantheas, J. Chem. Phys. 134, 121105 (2012).
  39. 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).
  40. Z. Ma, Y. Zhang,  and M. E. Tuckerman, J. Chem. Phys. 137, 044506 (2012).
  41. M. D. Ben, M. Schöner, J. Hutter,  and J. VandeVondele, J. Phys. Chem. Lett. 4, 3753 (2013).
  42. P. J. Feibelman, Phys. Chem. Chem. Phys. 10, 4688 (2008).
  43. B. Santra, J. Klimeš, D. Alfè, A. Tkatchenko, B. Slater, A. Michaelides, R. Car,  and M. Scheffler, Phys. Rev. Lett. 107, 185701 (2011).
  44. B. Santra, J. Klimeš, A. Tkatchenko, D. Alfè, B. Slater, A. Michaelides, R. Car,  and M. Scheffler, J. Chem. Phys. 139, 154702 (2013).
  45. F. Labat, C. Pouchan, C. Adamo,  and G. E. Scuseria, J. Comput. Chem. 32, 2177 (2011).
  46. O. Kambara, K. Takahashi, M. Hayashi,  and J.-L. Kuo, Phys. Chem. Chem. Phys. 14, 11484 (2012).
  47. E. D. Murray and G. Galli, Phys. Rev. Lett. 108, 105502 (2012).
  48. Y. Fang, B. Xiao, J. Tao, J. Sun,  and J. P. Perdew, Phys. Rev. B 87, 214101 (2013).
  49. M. Macher, J. Klimeš, C. Franchini,  and G. Kresse, J. Chem. Phys. 140, 084502 (2014).
  50. J. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
  51. B. Chen, I. Ivanov, M. L. Klein,  and M. Parrinello, Phys. Rev. Lett. 91, 215503 (2003).
  52. J. A. Morrone and R. Car, Phys. Rev. Lett. 101, 017801 (2008).
  53. M. Ceriotti, J. Cuny, M. Parrinello,  and D. E. Manolopoulos, Proc. Natl. Acad. Sci. USA 110, 15591 (2013).
  54. M. Ceriotti and D. E. Manolopoulos, Phys. Rev. Lett. 109, 100604 (2012).
  55. F. Paesani, S. Iuchi,  and G. A. Voth, J. Chem. Phys. 127, 074506 (2007).
  56. G. S. Fanourgakis, G. K. Schenter,  and S. S. Xantheas, J. Chem. Phys. 125, 141102 (2006).
  57. B. Santra, A. Michaelides,  and M. Scheffler, J. Chem. Phys. 127, 184104 (2007).
  58. B. Santra, A. Michaelides, M. Fuchs, A. Tkatchenko, C. Filippi,  and M. Scheffler, J. Chem. Phys. 129, 194111 (2008).
  59. B. Santra, A. Michaelides,  and M. Scheffler, J. Chem. Phys. 131, 124509 (2009).
  60. F.-F. Wang, G. Jenness, W. A. Al-Saidi,  and K. D. Jordan, J. Chem. Phys. 132, 134303 (2010).
  61. M. J. Gillan, F. R. Manby, M. D. Towler,  and D. Alfè, J. Chem. Phys. 136, 244105 (2012).
  62. A. Erba, S. Casassa, L. Maschio,  and C. Pisani, J. Phys. Chem. B 113, 2347 (2009).
  63. J. P. Perdew, M. Ernzerhof,  and K. Burke, J. Chem. Phys. 105, 9982 (1996a).
  64. C. Adamo and V. Barone, J. Chem. Phys. 110, 6158 (1999).
  65. S. Grimme, J. Comput. Chem. 25, 1463 (2004).
  66. S. Grimme, J. Antony, S. Ehrlich,  and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
  67. J. P. Perdew, K. Burke,  and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996b).
  68. A. D. Becke, Phys. Rev. A 38, 3098 (1988).
  69. W. Lee, C. Yang,  and R. G. Parr, Phys. Rev. B 37, 785 (1988).
  70. M. Dion, H. Rydberg, E. Schröder, D. C. Langreth,  and B. I. Lundqvist, Phys. Rev. Lett. 92, 246401 (2004).
  71. R. A. DiStasio Jr., B. Santra, and R. Car, to be published.
  72. A. Tkatchenko and M. Scheffler, Phys. Rev. Lett. 102, 073005 (2009).
  73. X. Wu, A. Selloni,  and R. Car, Phys. Rev. B 79, 085102 (2009).
  74. G. J. Martyna, M. L. Klein,  and M. Tuckerman, J. Chem. Phys. 97, 2635 (1992).
  75. D. J. Tobias, G. J. Martyna,  and M. L. Klein, J. Phys. Chem. 97, 12959 (1993).
  76. N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
  77. F. Tassone, F. Mauri,  and R. Car, Phys. Rev. B 50, 10561 (1994).
  78. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari,  and R. M. Wentzcovitch, J. Phys.: Condens. Matter 21, 395502 (2009).
  79. N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
  80. M. Sharma, Y. Wu,  and R. Car, Int. J. Quant. Chem. 95, 821 (2003).
  81. J. F. Dobson, Int. J. Quantum Chem. XXX, YYY (2014).
  82. J. M. Sorenson, G. Hura, R. M. Glaeser,  and T. Head-Gordon, J. Chem. Phys. 113, 9149 (2000).
  83. C. Huang, K. T. Wikfeldt, D. Nordlund, U. Bergmann, T. McQueen, J. Sellberg, L. G. M. Pettersson,  and A. Nilsson, Phys. Chem. Chem. Phys. 13, 19997 (2011).
  84. L. B. Skinner, C. Huang, D. Schlesinger, L. G. M. Pettersson, A. Nilsson,  and C. J. Benmore, J. Chem. Phys. 138, 074506 (2013).
  85. A. Soper and C. J. Benmore, Phys. Rev. Lett. 101, 065502 (2008).
  86. J. Wang, A. N. Tripathi,  and V. H. Smith Jr., J. Chem. Phys. 101, 4842 (1994).
  87. N. Watanabe, S. Ten-no, S. Pal, S. Iwata,  and Y. Udagawa, J. Chem. Phys. 111, 827 (1999).
  88. M. Krack, A. Gambirasio,  and M. Parrinello, J. Chem. Phys. 117, 9409 (2002).
  89. Z. Li, Improving Ab Initio Molecular Dynamics of Liquid Water, Ph.D. Thesis, Princeton University, 2012.
  90. The quantity is very sensitive to , the position of the first minimum in the , since is computed as , with as the corresponding atomic number density.
  91. F. Corsetti, E. Artacho, J. M. Soler, S. S. Alexandre,  and M.-V. Fernández-Serra, J. Chem. Phys. 139, 194502 (2013).
  92. J. R. Errington and P. G. Debenedetti, Nature 409, 318 (2001).
  93. C. Zhang, L. Spanu,  and G. Galli, J. Phys. Chem. B 115, 14190 (2011c).
  94. J. D. Eaves, J. J. Loparo, C. J. Fecko, S. T. Roberts, A. Tokmakoff,  and P. L. Geissler, Proc. Natl. Acad. Sci. (USA) 102, 13019 (2005).
  95. D. Prada-Gracia, R. Shevchuk,  and F. Rao, J. Chem. Phys. 139, 084501 (2013).
  96. A. Luzar and D. Chandler, Phys. Rev. Lett. 76, 928 (1996).
  97. To obtain the number of distinct hydrogen bonds per water molecule, the numbers provided in Table 1 should be divided by two.
  98. S. L. Shostak, W. L. Ebenstein,  and J. S. Muenter, J. Chem. Phys. 94, 5875 (1991).
  99. S. A. Clough, Y. Beers, G. P. Klein,  and L. S. Rothman, J. Chem. Phys. 59, 2254 (1973).
  100. 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).
  101. E. R. Batista, S. S. Xantheas,  and H. Jónsson, J. Chem. Phys. 111, 6011 (1999).
  102. P. L. Silvestrelli and M. Parrinello, J. Chem. Phys. 111, 3572 (1999a).
  103. P. L. Silvestrelli and M. Parrinello, Phys. Rev. Lett. 82, 3308 (1999b).
  104. M. Sharma, R. Resta,  and R. Car, Phys. Rev. Lett. 95, 187401 (2005).
  105. M. Sharma, R. Resta,  and R. Car, Phys. Rev. Lett. 98, 247401 (2007).
  106. D. D. Kemp and M. S. Gordon, J. Phys. Chem. A 112, 4885 (2008).
  107. Y. Tu and A. Laaksonen, Chem. Phys. Lett. 329, 283 (2000).
  108. A hydrogen bond is assigned to an MLWF if the distance between the center of the MLWF () and the donating hydrogen atom is less than 2.7 Å and the angle 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.
  109. L. Kong, X. Wu,  and R. Car, Phys. Rev. B 86, 134203 (2012).
  110. C. Swartz and X. Wu, Phys. Rev. Lett. 111, 087801 (2013).
  111. A. Tkatchenko, R. A. DiStasio Jr., R. Car,  and M. Scheffler, Phys. Rev. Lett. 108, 236402 (2012).
  112. R. A. DiStasio Jr., O. A. von Lilienfeld,  and A. Tkatchenko, Proc. Natl. Acad. USA 109, 14791 (2012).
  113. A. Tkatchenko, A. Ambrosetti,  and R. A. DiStasio Jr., J. Chem. Phys. 138, 074106 (2013).
  114. A. Ambrosetti, A. M. Reilly, R. A. DiStasio Jr.,  and A. Tkatchenko, J. Chem. Phys. 140, 18A508 (2014).
  115. R. A. DiStasio Jr., V. V. Gobre,  and A. Tkatchenko, J. Phys.: Condens. Matter 26, 213202 (2014).
103710
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
Edit
-  
Unpublish
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel
Comments 0
""
The feedback must be of minumum 40 characters
Add comment
Cancel
Loading ...

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description