Ab initio theory and modeling of water

Ab initio theory and modeling of water Mohan Chen,1 Hsin-Yu Ko,2 Richard C. Remsing,3, 4 Marcos F. Calegari Andrade,2 Biswajit Santra,2 Zhaoru Sun,1 Annabella Selloni,2 Roberto Car,2 Michael L. Klein,1, 3, 4, ∗ John P. Perdew,1, 3 and Xifan Wu1, 4, † 1

Department of Physics, Temple University, Philadelphia, Pennsylvania 19122, USA Department of Chemistry, Princeton University, Princeton, New Jersey 08544, USA 3 Department of Chemistry, Temple University, Philadelphia, Pennsylvania 19122, USA 4 Institute for Computational Molecular Science, Temple University, Philadelphia, Pennsylvania 19122, USA (Dated: October 2, 2017)

arXiv:1709.10493v1 [cond-mat.soft] 29 Sep 2017

2

Water is of the utmost importance for life and technology. However, a genuinely predictive ab initio model of water has eluded scientists. We demonstrate that a fully ab initio approach, relying on the strongly constrained and appropriately normed (SCAN) density functional, provides such a description of water. SCAN accurately describes the balance among covalent bonds, hydrogen bonds, and van der Waals interactions that dictates the structure and dynamics of liquid water. Notably, SCAN captures the density difference between water and ice Ih at ambient conditions, as well as many important structural, electronic, and dynamic properties of liquid water. These successful predictions of the versatile SCAN functional open the gates to study complex processes in aqueous phase chemistry and the interactions of water with other materials in an efficient, accurate, and predictive, ab initio manner. Keywords: water, ab initio theory, hydrogen bonding, density functional theory, molecular dynamics. PACS numbers:

I.

SIGNIFICANCE

Water is vital to our everyday life, but its structure at a molecular level is still not fully understood from either experiment or theory. The latter is hampered by our inability to construct a purely predictive, first principles model. The difficulty in modeling water lies in capturing the delicate interplay among the many strong and weak forces that govern its behavior and phase diagram. Herein, molecular simulations with a recently proposed non-empirical quantum mechanical approach (the SCAN density functional), yield an excellent description of the structural, electronic, and dynamic properties of liquid water. SCAN-based approaches, which describe diverse types of bonds in materials on an equal, accurate footing, will likely enable efficient and reliable modeling of aqueous phase chemistry.

Water is arguably the most important molecule for life and is involved in almost all biological processes. Without water, life, as we know it, would not exist, earning water the pseudonym matrix of life, among others1 . Despite the apparent simplicity of an H2 O molecule, water in the condensed phase displays a variety of anomalous properties which originates from its complex structure. In an ideal arrangement, water molecules form a tetrahedral network of hydrogen (H) bonds with each vertex being occupied by a water molecule. This tetrahedral network is realized in the solid phase ice Ih, but thermal fluctuations disrupt the H-bond network in the liquid state, with the network fluctuating on picosecond to nanosecond timescales. Due to the complexity of the H-bond network and its competition with thermal fluctuations, a precise molecular-level understanding of the structure of liquid water remains elusive. Major challenges lie in unambiguously capturing the atomic-scale fluctuations in water experimentally. Current approaches such as timeresolved spectroscopy2,3 and diffraction measurements4,5 may be able to resolve changes on picosecond timescales, but rely on interpretation through models, which often cannot describe all the details of liquid water with quantitative accuracy. Not surprisingly, the nature of the Hbond network in liquid water continues to be at the center of scientific debate and advances in both experiment and theory are needed, especially with regard to quantitative modeling of aqueous phase chemistry. Ab initio molecular dynamics (AIMD) simulation6 is an ideal approach for modeling the condensed phases of water across the phase diagram and aqueous phase chemistry using quantum mechanical principles7–11 , although for some applications, such as the study of liquid vapor phase equilibria12 , Monte Carlo methods are better

2 suited. In particular, Kohn-Sham density functional theory (DFT)13 — used to model the system in its electronic ground state — provides an efficient framework that enables the simulation of the length and time scales needed to converge many statistical mechanical averages in disordered, liquid state systems. The DFT formalism is exact for the electronic ground-state energy and density, but in practice approximations must be adopted to describe many-body effects, included in the exchange-correlation (XC) functional. XC functionals can be conceptually arranged, by accuracy and computational efficiency, according to Jacob’s ladder14 , with the simplest local density approximation (LDA)15,16 on the bottom rung of the ladder, followed by generalized gradient approximations (GGAs)17–19 , meta-GGAs, hybrid functionals20,21 , and so on. The past three decades have witnessed widespread successes of DFT in elucidating and predicting properties of materials. However, water still presents a major challenge, with many DFT-based simulations yielding results that are not even qualitatively consistent with experimental measurements. The H-bonds formed between gas-phase water clusters were first treated within the LDA22,23 , which overestimates H-bond strengths and yields inter-water distances that are too close. This overbinding is largely corrected by GGA-level functionals, which became a class of popular functionals to study liquid water within the last two decades 10 . Despite the improvements over LDA that are provided by GGAs, H-bond strengths are overestimated and, consequently, the dynamical properties predicted by GGAs are generally much too slow. Worse still, GGAs predict that ice sinks in water, that is, water has a lower density than ice11,11,24–26 . These disagreements remain even after considering hybrid functionals11 and accounting for nuclear quantum effects (NQEs)27 , illustrating that the deficiencies are a manifestation of errors within the underlying GGA to the XC functional. The difficulty in modeling liquid water with DFT arises from the delicate nature of the H-bond network. A Hbond is a directional attractive force between the oxygen of one molecule and the protons of another. While mainly electrostatic in nature, H-bonds also exhibit a non-negligible covalency. Notably, a covalent O-H bond binds one order of magnitude stronger than a H-bond in water. Therefore, a slightly misbalanced covalent bond inevitably incurs a non-negligible error in the predicted H-bond strength. Moreover, water molecules interact with each other through van der Waals (vdW) dispersion forces at larger distances, which are non-directional and in general weaker than H-bonds by roughly an order of magnitude. Thus, one needs to capture the balance among interactions whose magnitudes vary by orders of magnitude in water. The short-ranged portion of the vdW interactions have been captured by local and semi-local XC functionals. In contrast, the intermediate- and long-ranged parts of the vdW interactions have not been captured by any general-

purpose GGA. Recent studies have identified vdW interactions as an important determinant of water structure; vdW interactions often lead to more disordered water structures, more accurate water densities, and improved dynamic properties 11,24–26,28–30 . Thus, the H-bond network of liquid water is produced by a delicate competition among covalent bonds, H-bonds, and vdW interactions, and describing this complex interplay of interactions continues to be a highly challenging task. In this regard, non-empirical, general purpose XC functionals that describe all types of interactions on an equal footing are imperative but still largely absent in the literature. To address the above issues, we performed AIMD simulations of liquid water in the isothermal-isobaric ensemble31 , employing the strongly constrained and appropriately normed (SCAN) meta-GGA functional32 . SCAN is inherently non-empirical, developed by satisfying all 17 known exact constraints on semi-local XC functionals. Thus, the results obtained from SCAN are purely predictive and do not rely on training data. SCAN was shown to predict the energetics of gas-phase water hexamers and ice phases with quantitative accuracy, while other XC functionals, even with vdW corrections, were unable to make even qualitative predictions33 . This suggests that SCAN possesses the ingredients necessary to describe liquid water. Indeed, we demonstrate that SCAN predicts structural, electronic, and dynamic properties of liquid water in excellent agreement with experimental measurements. In particular, due to its ability to describe vdW interactions on intermediate length-scales, SCAN yields the correct density ordering between liquid water and ice, correctly predicting that ice floats on liquid water. The dynamics of liquid water are also improved to near quantitative agreement with experiments. We expect the computationally-efficient and accurate SCAN functional to serve as a major quantum mechanics-based tool for studying chemical processes in aqueous media.

II.

MOLECULAR AND ELECTRONIC STRUCTURE OF LIQUID WATER

The pair structure of liquid water can be measured by X-ray diffraction4,5 and neutron diffraction experiments4 , from which structural information is contained in the resulting radial distribution functions (RDFs). We compare the RDFs obtained from AIMD simulations with SCAN and the Perdew-Burke-Ernzerhof (PBE)19 GGA, as well as the experimental data. Here we compare two fully ab initio density functionals, without an empirical dispersion (D) correction to either. While such a correction improves PBE for solids and liquids34 , it slightly worsens PBE’s unacceptable overbinding of molecules, and thus PBE-D is not recommended for reactions in solvents. Figs. 1(a) and 1(b) show the oxygen-oxygen and oxygen-hydrogen RDFs, gOO (r) and gOH (r), respectively. SCAN dramatically improves almost all features in gOO (r) and gOH (r), producing a pair structure in much

3

Lone pair

Bonding pair

FIG. 1: Radial distribution functions (a) gOO (r) and (b) gOH (r) of liquid water predicted by PBE and SCAN at 330 K, as well as that from X-ray diffraction experiments5 for gOO (r) and joint X-ray/neutron diffraction experiments4 for gOH (r). An elevated temperature of 30 K was utilized in AIMD simulations to mimic NQEs35 .

better agreement with experimental measurements than PBE. The first peak of gOH (r) contains all correlations within the covalent O-H bonds. SCAN enhances the covalency of water molecules, shortening the covalent bond length to 0.977 ˚ A (first maximum in gOH (r)), in comparison to the 0.989 ˚ A from PBE. The shorter O-H bond length indicates that the oxygen and protons bind more strongly. Consequently, the protons of water molecules are less easily donated to form H-bonds. Correlations between H-bonded neighbors are contained in the first peak of gOO (r) and the second peak of gOH (r). As evidenced by Fig. 1, SCAN captures these correlations with high accuracy due to its ability to describe H-bonding. The region between the first and second peaks of gOO (r) predominantly consists of non-Hbonded water molecules that occupy the interstitial space between H-bonded neighbors; the increased number of water molecules in the interstitial regions is due to vdW interactions, as discussed further below. Subsequent coordination shells are also captured by SCAN, evidenced by the good agreement between the second and third peaks in gOO (r). We emphasize that the near perfect agreement between the SCAN gOO (r) and experiment is non-trivial, because the structure of water is a manifestation of the delicate interplay among covalent bonds, H-bonds, and vdW interactions. The strength of directional H-bonds is largely determined by the electronic structure of water molecules. The electronic density of states (DOS) of liquid water, averaged over trajectories, is shown in Fig. 2(a) and compared to the DOS measured by full valence band photoemission

FIG. 2: (a) Density of states (DOS) of liquid water, averaged over SCAN and PBE trajectories, as well as from photoemission spectroscopy43 . The peaks are labeled according to the symmetric orbitals of a water molecule with C2v symmetry. Data are aligned44 to the 1b1 peak of the experimental (EXP) data. (b) Distributions of the centers of maximally localized Wannier functions (MLWFs) with respect to the oxygen position for lone and bonding electron pairs. The inset shows a representative snapshot of the MLWFs of a water molecule; lone and bonding pair MLWFs are colored green and blue, respectively.

spectroscopy43 . The four peaks of the DOS are assigned to the 2a1 , 1b2 , 3a1 , and 1b1 orbitals based on the spatial symmetries of the water molecule. The simulated DOS are aligned at the position of the 1b1 orbital peak44 . The energy difference between the 2a1 peak predicted by PBE and experiment is 2.3 eV. SCAN substantially lowers this energy difference to 0.9 eV, providing a much better description of the strongly bound 2a1 orbital than the GGA-level description provided by PBE. Note that the strongly bound 2a1 orbital is mainly composed of the characteristic 2s orbital and is close to the oxygen atom. The above four orbitals are related to the two lone electron pairs and two bonding electron pairs of a water molecule; the lone electron pairs are closely connected to the 2a1 and 1b1 orbitals while the bonding electron pairs have a strong relation to the 1b2 and 3a1 orbitals. Therefore, the improved DOS by SCAN implies that the lone and bonding electron pairs are better captured than those from PBE. We examine the lone and bonding electron pairs on an equal footing through maximally localized Wannier functions (MLWFs)45 , which are generated from a unitary transformation of the occupied KohnSham eigenstates. Fig. 2(b) shows the distributions of the centers of the MLWFs. The lone electron pairs are closer to the oxygen atom in the SCAN description of wa-

4 TABLE I: Properties of water (330 K) and ice Ih (273 K) predicted by SCAN and PBE functionals in the isobaric-isothermal ensemble: densities of water (ρw ) and ice Ih (ρIh ), density difference (∆ρ), density ratio ρw /ρIh , dipole moments of water (µw ) and ice Ih (µIh ), band gap (Eg ), tetrahedral order parameter (q), diffusion coefficient (D), and rotational correlation time (τ2 ). The temperatures for experimental data (EXP) ρw , ρIh , µw , D, and τ2 are 30036 , 27337 , 29838 , 29839 , and 300 K40 , respectively. The experimental q value4 was obtained by combining X-ray diffraction at 296 K and neutron diffraction data at 298 K in a structural model using empirical potential structural refinement. No experimental data of µIh are found but an induction model gave rise to 3.09 D for µIh 41 . Experimental data for q, D, and τ2 are for D2 O chosen for consistency with the masses used in simulations for the dynamic properties. Error bars correspond to one standard deviation. Method ρw (g/mL) ρIh (g/mL) ∆ρ (g/mL) ρw /ρIh µw (D) µIh (D) SCAN 1.050±0.027 0.964±0.023 0.086±0.035 1.089±0.038 2.97±0.29 3.29±0.21 PBE 0.850±0.016 0.936±0.013 -0.086±0.021 0.908±0.021 3.12±0.28 3.35±0.21 EXP 0.9965636 0.916737 0.080 1.087 2.9±0.638

ter than PBE, while the bonding electron pairs only differ slightly between the two XC functionals. The smaller distance between lone electron pairs and oxygen in SCAN leads to a less negative environment around the lone electron pairs and explains the lower energy of the 2a1 orbital in comparison to that of PBE. Meanwhile, the nearly unchanged description of bonded electron pairs in the two functionals is consistent with the observation that 1b2 and 3a1 states are also similar. Consequently, electrostatic attractions between oxygen nuclei and protons of neighboring water molecules are weaker in SCAN than PBE, weakening the directional H-bond strength. In addition to improving the intermolecular structure, the reduced H-bond strength in SCAN also improves the intramolecular structure of water. The shorter distance between the lone electron pairs and the oxygen nucleus weakens the capability to accept H-bonds and water molecules become less polarizable. The reduction in polarizability is expected to improve other electronic properties of liquid water, moving them in closer agreement with experimental measurements. Indeed, the dipole moment µ of liquid water, computed via MLWFs, is reduced by SCAN. Table 1 shows that µ = 3.12 D with PBE, while µ reduces to 2.97 D with SCAN, in better agreement with experimental measurements of 2.9±0.6 D38 . This improvement indicates that the important dipoledipole interactions in liquid water are better described by SCAN. We also estimate the band gap of water, Eg , by averaging over eight randomly selected configurations from the trajectories. SCAN and PBE predict Eg = 4.92 and 4.43 eV, respectively. While SCAN improves Eg by about 0.5 eV, it differs significantly from the experimental value of 8.7 eV42 . We attribute this discrepancy to the well-known underestimation of band gaps by GGAs and meta-GGAs. The SCAN functional can describe the intermediateranged vdW interactions33 , which shift the first minimum and the second maximum of gOO (r) toward the first peak, with respect to that of PBE (without vdW interactions), as shown in Fig. 1(a). Water molecules beyond the first coordination shell experience non-directional attractions from surrounding water molecules in SCAN and are pulled into the interstitial spaces between H-bonded

Eg (eV) q D (˚ A2 /ps) τ2 (ps) 4.92±0.14 0.68±0.18 0.190±0.025 2.9±0.4 4.43±0.13 0.83±0.11 0.018±0.002 7.1±0.5 8.7±0.642 0.5934 0.18739 2.440

waters by these vdW forces. Consequently, the peak position of the second coordination shell shifts inward towards the central oxygen, and the population of interstitial waters increases, illustrated by the increase in the height of the first minimum in gOO (r). Thus, the inclusion of non-directional vdW interactions on intermediate length-scales leads to a more disordered and highlypacked water structure. From the increased packing, one expects the density of liquid water predicted by SCAN to be larger than that from PBE. Moreover, the dominant effect of vdW interactions is to provide cohesive interactions between molecules in condensed phases. Within the vdW picture of liquids, this leads to a cohesive pressure of magnitude −aρ2w , which “squeezes” water molecules closer together29 ; a is the vdW constant and a measure of the strength of these attractive interactions, and ρw is the density of liquid water. Indeed, the SCAN functional predicts ρw to be significantly higher than that predicted by PBE, as shown in Table 1. Another problem of paramount importance is that solid water, ice Ih, floats on liquid water near ambient conditions. This is probably the most widely known anomalous property of water. Yet, almost all DFT-based approaches, except some of those relying on empirical parameters, predict a solid phase that is denser than the liquid. In this regard, we also carried out AIMD simulations of ice Ih containing 96 water molecules at 273 K. SCAN predicts a ρw that is larger than the density of ice Ih (ρIh ), Table 1, while for PBE, ρIh > ρw . The water density from SCAN is ' 5% larger than that determined experimentally, which is a significant improvement over the 15%, 25%, and 39% underestimation by the PBE, BLYP17,18,24 , and PBE011 functionals, respectively. Compared to PBE, SCAN increases ρw and ρIh by 21% and 3%, respectively. The 21% increase of ρw by SCAN is vital in correcting the density ordering between the two phases by other functionals. Indeed, the experimental density difference between liquid water and ice Ih, ∆ρ = ρw − ρIh , is correctly predicted by SCAN as 0.086 g/mL, while the opposite sign is predicted by PBE. SCAN also correctly predicts ρw /ρIh ' 1.089, in agreement with the experimental value of ' 1.087, and in contrast to the 0.908 ratio via PBE. Note that the wa-

5

(a)

(b)

(d)

PBE

PBE Tetrahedral Structure

(c)

(e) SCAN

F (kBT) 3

0

FIG. 3: (a) Distributions of the number of hydrogen bonds in liquid water from SCAN and PBE. The inset illustrates an ideal tetrahedral H-bonding structure. Oxygen and hydrogen atoms are respectively depicted in red and white; H-bonds are shown with dashed lines. (b,c) Bond angle distributions POOO (θ) from (b) PBE and (c) SCAN. POOO (θ) is decomposed into contributions arising from waters with a fixed number of HBs (2, 1, and 0) between a central oxygen and its two nearest neighbors. The experimental POOO of D2 O is inferred from experiments4 , and the area of POOO is normalized to unity. (d,e) Free energies (F ) as a function of θ and the oxygen-oxygen distance r from (d) PBE and (e) SCAN. The free energy minimum is identified by the red circle and referenced to zero. The direction of change of the free energy minimum with increasing r is shown with a red arrow. The cutoff distance used for computing the free energies is the same as that for POOO and is shown with a dashed red line.

ter density obtained by PBE is slightly lower than the results of previous studies and we discuss the difference in the SI.

III.

TETRAHEDRAL STRUCTURE OF THE H-BOND NETWORK

With the translational order encoded in RDFs well captured by SCAN, we now focus on the tetrahedral orientational ordering of liquid water induced by the Hbond network. An ideal tetrahedral H-bonding structure shown in the inset of Fig. 3(a) is formed because a water molecule can possess four optimal H-bonds: two accepting and two donating. Thermal fluctuations break and reform H-bonds, causing the tetrahedral structures in liquid water to be distorted or broken by entropic effects. This, combined with the increased packing due to vdW interactions, leads to an average number of H-bonds per molecule slightly less than four in liquid water. To illustrate the impact of SCAN on the H-bond network, distributions of the number of H-bonds per molecule are presented in Fig. 3(a). The percentage of water molecules participating in four H-bonds drops from 72% in PBE to 56% in SCAN. This suggests that H-bonds are weaker with SCAN than with PBE. SCAN predicts an average of 3.61 H-bonds per molecule, smaller than the 3.77 obtained from PBE. This reduction in the number of H-bonds is consistent with the influences of the underlying SCAN functional on liquid water: directional H-bonds are weakened and more easily broken by thermal

fluctuations. The increased disorder is further stabilized by the inclusion of the intermediate-ranged vdW interactions naturally arising in SCAN. The reduction of H-bonds produced by SCAN disrupts the tetrahedral structure of liquid water. To quantify the amount of tetrahedral order, we adopt the tetrahedral order parameter q 9 . A perfect tetrahedral local environment corresponds to q = 1, and q decreases as the local structure becomes less tetrahedral. Following experimental work4 , we evaluate q using a cutoff radius that yields an average coordination number of 4. The resulting cutoffs are 3.15 and 3.45 ˚ A for SCAN and PBE, respectively, with SCAN in better agreement with the cutoff of 3.18 ˚ A inferred from experiment4 . Despite the high first peak in the PBE gOO (r), the shorter cutoff from SCAN suggests a more compact first coordination shell, consistent with the higher density of liquid water it predicts. PBE results in an overly tetrahedral liquid (Table 1). SCAN, however, yields q in better agreement with experiments on heavy water4 , suggesting that SCAN provides a more accurate structural description of the fluctuating H-bond network. Three-body correlations in water can be quantified by the bond angle distribution POOO (θ), where θ is the angle formed by an oxygen of a water molecule and two of its oxygen neighbors; neighbors are defined using the same cutoff as above4 . The PBE POOO in Fig. 3(b) displays a high peak centered around the tetrahedral angle, 109.5◦ and is a much narrower distribution than that from experiment. This indicates that PBE overestimates the tetrahedral character of the liquid, consistent with the

6 above-described overstructuring. In stark contrast, the SCAN POOO is in excellent agreement with experiment, with almost exactly the same widths and intensities of the two peaks close to 109.5◦ and 55◦ , Fig. 3(c). The peak located near 109.5◦ arises from tetrahedral structures. The peak at θ ' 55◦ is related to broken H-bonds and interstitial, non-H-bonded water, and major differences between SCAN and PBE are observed in this region of the distribution. We decompose POOO (θ) into three contributions according to the number of Hbonds a water molecule formed within a water triplet. The POOO (θ) of triplets formed with 2, 1, and 0 H-bonds are plotted in Figs. 3(b) and 3(c). Triplets involving 2 H-bonds dominate the PBE POOO (θ), while triplets with broken (0 or 1) H-bonds contribute much less. In contrast, triplets with less than 2 H-bonds contribute significantly to the POOO (θ) predicted by SCAN, especially near θ ≈ 55◦ . The free energy as a function of θ and the distance r between neighboring oxygen atoms in the triplet reveals additional insights, as shown in Figs. 3(d) and 3(e). As expected, the minimum free energy corresponds to tetrahedral-like structures with θ ' 109.5◦ and r ≈ 2.7 ˚ A. In contrast to PBE, SCAN predicts a significant fraction of triplets with θ far from 109.5◦ , indicating that the SCAN liquid is more disordered. The free energies suggest that water molecules in the first coordination shell experience a smaller free energy barrier to adopt a broad range of θ-values with SCAN than with PBE. Importantly, there are substantial differences between the two functionals in describing the dependence of the free energy on r. With PBE, as r is increased away from the free energy minimum, θ hardly moves from 109.5◦ , as depicted by the red arrow in Fig. 3(d). This is consistent with the over-structuring of water by PBE and implies that θ is weakly influenced by fluctuations of the first coordination shell. In contrast, SCAN produces a stronger correlation between r and θ, such that the free energy is lowered at larger r by decreasing θ, illustrated by the red arrow in Fig. 3(e). This is consistent with the higher population of non-H-bonded, interstitial waters in the SCAN prediction. These non-tetrahedrally oriented water molecules contribute significantly to POOO (θ) below 109.5◦ and highlight the reduced tetrahedrality of the SCAN H-bond network.

IV.

DYNAMICS

Changes in the H-bond energy alter the delicate enthalpy-entropy balance in liquid water that dictates its dynamic properties; for example, breakage and formation of H-bonds through thermal fluctuations controls diffusion. Thus, stronger H-bonds tilt the enthalpy-entropy balance toward energetic contributions, reducing the tendency to break H-bonds and consequently lowering the diffusion coefficient D. We estimate D from the longtime limit of the mean squared displacement, averaged

over the oxygen and hydrogen atoms (see SI). Indeed, the D value of PBE is an order of magnitude smaller than that of experiment, while SCAN improves the estimate of D to near agreement with experiment. H-bond dynamics are more directly probed via the second-order rotational correlation function of the O-H bond vector rOH , C2 (t) = hP2 (rOH (t) · rOH (0))i/hP2 (rOH (0)2 )i, where P2 (x) is a second-order Legendre polynomial. The integral of C2 (t) yields the rotational correlation time τ2 of the O-H bond; correlation functions and details surrounding τ2 computation are given in the SI. SCAN predicts a value of τ2 in agreement with nuclear magnetic resonance spectroscopy40 , Table 1, while rotational dynamics are slowed in the PBE system. The mechanism for rotational relaxation of the O-H bond vector is associated with breaking a H-bond. In PBE, H-bonds are too strong, significantly hindering this pathway. SCAN appropriately predicts the weight of these pathways for rotational relaxation due to its accurate description of H-bonding.

V.

CONCLUSIONS AND OUTLOOK

The SCAN density functional provides a genuinely predictive ab initio model of liquid water. Importantly, SCAN is a long-awaited exchange-correlation functional that can correctly predict liquid water that is denser than ice at ambient conditions. SCAN excellently describes covalent and H-bonds due to an improved description of electronic structure, and captures intermediate-ranged vdW interactions that further improve the structure and thermodynamics of liquid water. These vdW forces can play a critical and active role at interfaces, for example, underlying drying transitions28,29,46 , instilling confidence that SCAN will enable predictive modeling of heterogeneous chemical environments. However, there are still improvements to be made regarding the water structure. SCAN predicts a slightly over-structured first peak of gOO (r). Previous studies have attributed the over-structuring to self-interaction errors16 , which can be mitigated by including a fraction of exact exchange in hybrid functionals. Moreover, the first peak in gOH (r) is too narrow, and the error is dominated by the lack of NQEs of hydrogen35 . The widths and intensities of peaks in the computed DOS are also respectively narrower and higher than those in the experimental DOS. In fact, DFT is not rigorous for photoemission spectra, and does not include lifetime broadening; NQEs, however, can additionally broaden the DOS bringing the resulting widths and intensities in closer agreement with experiment44 . NQEs can be accounted for within the Feynman discretized path-integral approach27,35,44 . In conclusion, the SCAN XC functional within DFT shows promising predictive power and will likely enable confident ab initio predictions for complex systems at the forefront of physics, chemistry, biology, and materials science47 .

7 VI.

METHODS

We performed Car-Parrinello molecular dynamics6 in Quantum ESPRESSO48 . We employed the Hamann-Schl¨ uter-Chiang-Vanderbilt pseudopotentials49 generated using PBE. The valence electrons, including the 1s electron of H and the 2s2 p4 electrons of O, were treated explicitly. The energy cutoff was 150 Ry. Simulations were performed in the isothermal-isobaric ensemble (constant NpT) by using the Parrinello-Rahman barostat31 and a single Nos´e-Hoover thermostat50 with a frequency of 60 THz to maintain a constant pressure (p) and temperature (T ), respectively. T = 330 K for liquid water and 273 K for ice Ih; the 30 K increase above ambient conditions in the former mimics NQEs on the liquid structure35 . We adopted a cubic cell with N =64 water molecules. The fictitious mass of the electrons was set to 100 au and the corresponding mass pre-conditioning with a kinetic energy cutoff of 25 Ry was used to all Fourier components of wavefunctions51 . The deuterium mass was used instead of hydrogen to enable the use of a timestep of 2 au; dynamics were compared to D2 O instead of H2 O. SCAN and PBE trajectories for water were 30.0 and 20.0 ps in length, respectively. Corresponding trajectories for ice Ih were 11.1 and 13.8 ps, respectively. The first 5 ps of each trajectory was used for equilibration and the remainder used for analysis. We utilized a standard geometric criterion for hydrogen bonding; covalently bonded O-H are associated with an O-H distance less than 1.24 ˚ A and H-bonds have an O-O distance less than 3.5 ˚ A and a 6 OOH angle less than 30◦52 . Additional details are in the SI.

VII.

AUTHOR CONTRIBUTIONS

X.W. and M.C. designed research. M.C. and Z.S. performed research. H.-Y.K., M.F.C.A., and B.S. contributed new methods. M.C. and R.C.R. analyzed data. All authors contributed to writing the paper.

VIII.

ACKNOWLEDGEMENTS

This work was supported by the U.S. Department of Energy (DOE) SciDac under Grant #de-sc0008726. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. DOE under Contract #DE-AC02-05CH11231. RCR, ZS, and JPP were supported as part of the Center for the Computational Design of Functional Layered Materials, an Energy Frontier Research Center funded by the U.S. DOE, Office of Science, Basic Energy Sciences under Award #de-sc0012575. MFCA is partially supported by the CNPq - Brazil. XW is partially supported by the National Science Foundation (NSF), DMR under Award #DMR-1552287.

8 SUPPORTING INFORMATION A.

Computational Details

The NpT algorithm was implemented in the Quantum ESPRESSO48 package. In our water simulations, all of the plane waves {G} with kinetic energies below 150 Ry were included, and we followed Ref.53 to maintain a constant plane wave kinetic energy cutoff of E0 =130 Ry for a fluctuating cell by adding a smooth step function with height A = 200 Ry and width σ = 15 Ry to the plane 2 0 )], wave kinetic factor as G2 → G2 + A[1 + erf( G /2−E σ where erf is the error function. The reference cells were chosen to be cubes with side lengths of 14.3345 ˚ A and 12.6579 ˚ A for PBE- and SCAN-based AIMD simulations, respectively. We ran parallel AIMD simulations using both strongly constrained and appropriately normed (SCAN)32 and Perdew-Burke-Ernzerhof (PBE)19 exchange-correlation functionals on 216 computer cores and recorded the wall times of 100 MD steps. The system was bulk liquid water consisting of 64 water molecules as utilized in this work. Both simulations were carried out on nodes with 2 x Intel Ivy Bridge @ 2.4 GHz and up to 64 GB RAM. We obtained 6.44 and 3.89 seconds/step for SCAN and PBE functionals, respectively, with SCAN being only 1.66 times more costly than PBE. Therefore, we conclude the cost of using the SCAN functional in studying liquid water is not dramatically more expensive than that of using the PBE functional and can be considered on the same level.

B.

Bulk Densities

Bulk densities as a function of time are shown in Fig. 4 for the SCAN and PBE descriptions of liquid water and ice Ih. The ice Ih phase remained solid throughout and did not transform to the liquid phase in the AIMD trajectories. The averaged densities are listed in Table 1 in the main text with the errors bars corresponding to one standard deviation. With PBE the dynamics of ambient liquid water is very sluggish and we find that at 330 K the mean density of liquid water (0.85 g/mL as shown in Table 1 in the main text) computed with PBE can vary within 0.01 g/mL (as estimated by another independent run of more than 60 ps by members of Roberto Cars group with PBE and all the same parameters) depending on the initial configuration and trajectory length. Since PBE liquid density is well below PBE ice Ih density we did not try to reduce the statistical uncertainty in the PBE liquid density. In addition, a recent study found structural, dynamical, and electronic properties of liquid water as obtained by AIMD simulations utilizing CP2K54 and Quantum ESPRESSO packages55 compared well, with the latter equilibrium density 0.02 g/mL lower, and gOO (r) first-

peak 0.1 higher, than corresponding CP2K values. These differences are comparable to the statistical uncertainties. Therefore the small difference in the PBE liquid density found here: 0.85 versus 0.865-0.887 g/mL by24 , as well as the slightly higher gOO (r) first-peak: 3.61 versus 3.363.5424 is attributed to both statistical uncertainties and numerical differences of these approaches. The densities clearly fluctuate around an average value for each trajectory, illustrating equilibration of the trajectories. Moreover, the fact that water is denser than ice Ih in the SCAN prediction is clearly observed, while the opposite is found for PBE. Finally, we note that fluctuations in the density of water are larger in the SCAN trajectory than in the PBE trajectory. This indicates that water is more compressible in the SCAN description than PBE, which produces a more rigid and ordered liquid structure. 1 .1 5 (a ) 1 .1 0

D e n s ity ( g /m L )

IX.

1 .0 5 1 .0 0 0 .9 5 0 .9 0 W

a te a te Ic e Ih Ic e Ih W

0 .8 5 0 .8 0 0

5

1 0

1 5

2 0

r (S r (P (S (P

2 5

C A B E C A B E

N , 3 3 0 , 3 3 0 K N , 2 7 3 , 2 7 3 K

K ) ) K ) )

3 0

T im e ( p s )

FIG. 4: Density fluctuations of liquid water and ice Ih as obtained from both SCAN and PBE trajectories using the isobaric-isothermal ensemble (NpT). A relatively shorter trajectory was generated for ice Ih because its density of solid phase converges quickly. The PBE functional incorrectly predicts that ice Ih (green line) is denser than water (orange line), while the SCAN functional successfully captures the larger density of water (blue line) than that of ice Ih (pink line). The black dashed and dotted lines represent the approximate experimental values of water density at ambient conditions (300 K) and ice density at 273 K under ambient pressure. An additional 30 K was applied to water to mimic the nuclear quantum effect35 . The averaged densities are listed Table 1 in the main text.

C.

van der Waals Interactions in SCAN and PBE

As discussed extensively in the main text, an accurate description of water and ice from first principles is challenging because the H-bond network of water arises from a delicate balance of strong intra-molecular covalent bonds, weak inter-molecular H-bonds, and even weaker vdW interactions. GGA functionals exhibit delocalization problems and the intermediate- and longranged vdW attraction is strongly underestimated. To be specific, the exchange energy density obtained from GGA is much more negative than LDA in regions with a

9

S C A N P B E 1 0

5

0 0

2

4

6

t (p s )

8

1 0

1 2

FIG. 5: Mean squared displacements for the SCAN and PBE systems consisting of 64 water molecules. Shaded regions indicate one standard error.

E.

Rotational Time Correlation Functions

FIG. 6: Second-order rotational time correlation functions, C2 (t), for the O-H bond vector of water as described by SCAN and PBE. Panel (b) is the same as (a), but zoomed in on the region from 0 to 0.5 ps.

The second-order rotational correlation function for the O-H bond vector rOH was calculated according to C2 (t) = hP2 (rOH (t) · rOH (0))i/hP2 (rOH (0)2 )i, where

3 .0

3 2 H

2 .5

2

O , S C A N -N V T , 3 3 0 K Q u a n tu m V A S P

2 .0

E S P R E S S O

1 .5

g

O O

large reduced density gradient. Therefore, the attractive vdW interactions between water molecules are missing in GGAs, which results in more ordered water molecules and a lower bulk density. The SCAN functional captures this delicate balance between covalent bonds, H-bonds, and vdW interactions in water, which is critical in accurately describing the water density. By including a dimensionless variable α in the kinetic energy density, SCAN can reduce to different GGAs by recognizing covalent single bonds when α = 0, slowly varying densities when α ' 1, and non-covalent bonds when α > 1. It was recently demonstrated that SCAN captures the intermediate-ranged vdW interactions for a variety of materials33 . To further illustrate this point, we applied the Tkatchenko-Scheffler (TS) vdW scheme56 to both the PBE and SCAN functionals. The TS scheme determines the ”turn-on” radius for atom pairs based on the XC functional used, and a larger TS scaling parameter SR implies that a larger radius is adopted to turn on the vdW interactions. The scaling parameters SR were obtained by fitting to the S22 molecular database. We find that in both water and ice structures, the scaling parameters are 0.94 and 1.17 for the PBE and SCAN functionals, respectively. The 24.5% larger scaling parameter in SCAN demonstrates that it captures the vdW interactions out to significantly larger distances than PBE.

gion of the MSDs was performed to obtain the D values. Finally, the obtained D values were averaged and the result is listed in Table 1 of the main text; the slope of the linear region is equal to 6D. Clearly, water diffuses much faster in the SCAN description than PBE, due to the weaker H-bonds predicted by the SCAN functional. The computed D from SCAN and PBE are 0.190 and 0.018 ˚ A2 /ps, respectively. The diffusion coefficient from PBE is close to the 0.020 ˚ A2 /ps reported in a previ9 ous work . However, we note that correlations may exist in divided sections and affect the accuracy of computed diffusion coefficients57,58 . More accurate diffusion coefficients require longer simulation times and we leave this investigation for future work.

(r)

M e a n S q u a r e D is p la c e m e n t ( Å

2

)

1 5

1 .0

D.

Mean Squared Displacements

0 .5 0 .0

The mean squared displacements (MSDs) computed from SCAN and PBE trajectories are shown in Fig. 5, where the center-of-mass positions of water molecules were used to compute MSDs. Each trajectory was divided into sections to compute the MSDs, each 12 ps in length and separated by 3.0 ps. We chose five and three sections for SCAN and PBE trajectories, respectively. Next, a linear fitting of the long-time, linear re-

2 .5

3 .0

3 .5

r (Å )

4 .0

4 .5

FIG. 7: Oxygen-oxygen radial distribution functions gOO (r) for 32 water molecules in condensed phase as obtained from the SCAN functional implemented in VASP and Quantum ESPRESSO electronic structure packages.

10

F.

Validation of SCAN in Different Packages

To further validate our results employing the SCAN functional for liquid water, we ran AIMD simulations on a cell of 32 water molecules for 20 ps by employing both Vienna Ab initio Simulation Package (VASP)59 and Quantum ESPRESSO48 packages. The cell was chosen to be a cube of side length 9.877 ˚ A. We adopted the NVT ensemble with the Nos´e-Hoover thermostat and the temperature was set to 330 K. We used the mass of deuterium instead of hydrogen to speed up the convergence. In VASP, we used projector-augmented-wave (PAW) potentials with configurations of [O]2s2 2p4 and [H]1s1 . In particular, we chose the hard PAW potentials for oxygen and hydrogen atoms and set the energy cutoff to 1200 eV in order to converge our results. The BornOppenheimer molecular dynamics was performed with a time step of 0.5 fs. In Quantum ESPRESSO, we carried out Car-Parrinello molecular dynamics6 and the settings were chosen to be the same as described in the main text. In Fig. 7, the oxygen-oxygen radial distribution functions gOO (r) are shown for the above two calculations. We find both electronic structure packages yield almost the same gOO (r) features. In particular, the first peak from both packages are almost identical. The results suggest that the properties of liquid water as predicted by

the SCAN functional are reliable and are reproducible with converged basis set and electron dynamics.

G.

Radial Distribution Function gOH

We plot in Fig. 8 the zoomed-in first peak of radial distribution function gOH from both PBE- and SCAN-based AIMD simulations. The first peak position represents the length of O-H covalent bond and we can see that SCAN predicts a slightly shorter covalent bond (0.977 ˚ A) than that from PBE (0.989 ˚ A). P B E (3 3 0 K ) S C A N (3 3 0 K )

4 0

O H

(r)

3 0

2 0

g

P2 (x) is a second-order Legendre polynomial. These time correlation functions as predicted for water by PBE and SCAN are shown in Fig. 6. Clearly, SCAN results in significantly faster rotational dynamics, evidenced by the much faster decay of the SCAN C2 (t) than that of PBE. As discussed in the main text, this is because SCAN results in weaker and more physical hydrogen bonding interactions than PBE. The short-time behavior of C2 (t) is shown in Fig. 6(b). We find that the initial short-time decay (