Density and Compressibility of Liquid Water and Ice from First

Jun 7, 2015 - As Figure 1 shows, substantial pressure oscillations in individual runs persist for ∼10 ps but level out afterward. Simulating an ense...
0 downloads 9 Views 705KB Size
Subscriber access provided by UNIV OF CALIFORNIA SAN DIEGO LIBRARIES

Letter

Density and Compressibility of Liquid Water and Ice from First-Principles Simulations with Hybrid Functionals Alex P. Gaiduk, Francois Gygi, and Giulia Galli J. Phys. Chem. Lett., Just Accepted Manuscript • Publication Date (Web): 07 Jun 2015 Downloaded from http://pubs.acs.org on June 7, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry Letters is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

Density and Compressibility of Liquid Water and Ice from First-Principles Simulations with Hybrid Functionals Alex P. Gaiduk,∗,† Fran¸cois Gygi,‡ and Giulia Galli∗,† †Institute for Molecular Engineering, The University of Chicago, Chicago, IL 60637 ‡Department of Computer Science, University of California, Davis, Davis, CA 95616 E-mail: [email protected]; [email protected]

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract We determined the equilibrium density and compressibility of water and ice from first-principles molecular dynamics simulations using gradient-corrected (PBE) and hybrid (PBE0) functionals. Both functionals predicted the density of ice to be larger than that of water, by 15% (PBE) and 35% (PBE0). The PBE0 functional yielded a lower density of both ice and water with respect to PBE, leading to better agreement with experiment for ice but not for liquid water. Approximate inclusion of dispersion interactions on computed molecular-dynamics trajectories led to a substantial improvement of the PBE0 results for the density of liquid water, which, however, resulted to be slightly lower than that of ice.

Graphical TOC Entry

Keywords Ab initio molecular dynamics; hybrid functionals; PBE0; water; ice; equilibrium density; compressibility.

2 ACS Paragon Plus Environment

Page 2 of 25

Page 3 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

Water is one of the most important liquids on our planet, yet several of its properties are not well understood. 1–4 Experimental techniques such as neutron or X-ray scattering 5–9 have provided valuable information on the structure of water; however, experiments give access only to average structural properties of liquids. Atomistic simulations complement measurements yielding a detailed description of the microscopic properties of liquid water and are useful not only as predictive tools but also to interpret a variety of experimental results. Most simulations of water and aqueous solutions appeared in the literature of the last decades have relied on empirical force fields to model interatomic interactions. 10–12 Such models, while being computationally affordable, are sometimes not transferable between thermodynamic conditions different from those they were fitted for. First-principles simulations use approximate solutions of the basic equations of quantum mechanics to determine the forces acting on atoms, 13–16 and provide a consistent description of all interactions in a given system, e.g. water and aqueous solutions. Most first-principles simulations are based on density-functional theory, especially on generalized gradient approximations (GGA) to the exact exchange-correlation functional EXC . 17,18 These approximations describe electron-electron interactions fairly accurately in a variety of systems; 19 however, they suffer from two important drawbacks—presence of spurious delocalization errors 20–22 and lack of dispersion interactions. In the case of water, the delocalization error is likely responsible, at least in part, for excessively strong hydrogen bonds between molecules and, as a result, for an overstructured liquid. The lack of dispersion interaction is partially responsible for an underestimate of the equilibrium density of liquid water 23,24 —common GGAs such as PBE 25,26 and BLYP 27,28 predict an equilibrium density that is 20–30% lower than in the experiment. The first ab initio determination of the equilibrium density of water was reported by McGrath et al., 29–31 who computed a density of 0.8 g/ml in isobaric–isothermal (NPT) Monte–Carlo simulations with the BLYP functional. Schmidt et al. 32 later performed NPT molecular dynamics (MD) simulations with the BLYP and PBE functionals and determined

3 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

similar mass densities: 0.75 (BLYP) and 0.88 g/ml (PBE). A recent series of simulations by Wang et al. 23 and Corsetti et al. 33 focused on the effect of dispersion interactions. 34,35 These authors found that in the absence of dispersion corrections revPBE, 36 BLYP, and PBE functionals predict equilibrium densities of 0.63–0.75, 0.76–0.85, and 0.85–0.90 g/ml, respectively. 23 Accounting for van der Waals interactions with nonlocal vdW-DF 34 and VV10 35 functionals yielded a density of 1.02 g/ml for vdW-DF-corrected revPBE, 23 and 1.10–1.20 g/ml for vdW-DF-corrected PBE and for VV10 functionals. 33,37 The equilibrium density of hexagonal ice (Ih) computed with revPBE, BLYP, and PBE functionals was found to be 0.89, 0.91, 38 and 0.97–0.99 g/ml, 39,40 respectively. Some of these values (revPBE, BLYP) are below the experimental ice density of 0.93 g/ml, 41 others (PBE) are above it. We note that, among these three functionals, the density of ice and liquid water described with revPBE is the lowest, while the density described with PBE is the highest. This probably means that the strength of hydrogen bonds increases in the row revPBE, BLYP, PBE. The deviation of the computed densities from the corresponding experimental values was much smaller for ice than for liquid water, resulting, at least in part, from the larger bulk modulus of the solid (12.1 GPa 42 ) compared to that of the liquid (2.2 GPa 43 ). The vdW-DF dispersion-corrected functional predicted the density of ice Ih to be 0.87 g/ml; 39 more recent version of this approximation (vdW-DF2 44 ) yielded a higher density of 0.89–0.90 g/ml. 39,45 PBE corrected by the dispersion functional of Ref. 46 yielded the equilibrium density of ice in worse agreement with experiment (1.01–1.02 g/ml). 39 Hybrid functionals 47 (e.g. PBE0 48,49 ) include a fraction of Hartree–Fock exchange in EXC , which alleviates the delocalization error present within GGAs. Recent studies 50,51 showed that hybrid functionals yield weaker hydrogen bonds in liquid water compared to GGAs, and provide a more realistic description of its structure. Whether the observed change in the structure also affects its equilibrium mass density is yet unclear: The only work 52,53 appeared so far that reported the equilibrium water density at the PBE0 level of theory, treated the Hartree–Fock energy functional approximately 54 and included Grimme

4 ACS Paragon Plus Environment

Page 4 of 25

Page 5 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

dispersion corrections. 55 That study found a negligible difference in the density of water simulated using PBE and PBE0 functionals, with the PBE0 result slightly (∼ 0.03 g/ml) lower. A similar small decrease in the equilibrium density of ice Ih (0.006 g/ml) was reported when PBE0 was used instead of PBE. 39 The goal of the present work was to determine the density and compressibility of liquid water and ice using the PBE0 functional and compare them with the results obtained at the PBE level of theory. Using ab initio MD simulations with plane-wave basis sets and pseudopotentials, we show that both PBE and PBE0 approximations overestimate the density of the cubic (Ic) and hexagonal (Ih) ice but underestimate that of liquid water. The PBE0 hybrid functional predicts a density lower than the PBE functional, leading to a better agreement with experiment for the solid but worse agreement for the liquid. Approximate inclusion of dispersion interactions in our calculations appeared to improve the agreement with experiment for the liquid when using hybrid functionals, but is likely to be accompanied by a worsening of the results for ice. 39 The equilibrium density is the point at which the internal pressure of the system is equal to the external, atmospheric, pressure. In this work, internal pressure P was computed as the average of the diagonal elements of the stress tensor σ: 56

P =

1 (σxx + σyy + σzz ) . 3

(1)

In the case of ice, the equilibrium density was determined by relaxing the positions of atoms and the periodic cell size; at equilibrium, the forces acting on nuclei and the elements of the stress tensor vanish. In the case of water, the density was derived from Born– Oppenheimer simulations in the canonical (NVT) ensemble. NPT simulations, common for liquids, usually require longer equilibration times than NVT ones, 32 which imposes an additional computational burden on ab initio MD. In addition, when the size of the periodic cell varies in NPT simulations, the plane-wave basis set resolution changes; this may lead

5 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

to an inconsistent description of the system at different densities 29,32 if the chosen cutoff is not large enough to describe with the same accuracy the pressure of the system at different densities. One solution to this problem is to carry out simulations at constant plane-wave resolution, by using a large reference cell and a confinement potential. 16,57–59 However, the use of a large reference cell along with an increased energy cutoff is not optimal from a computational standpoint, especially in the case of hybrid functionals. In this work, the equilibrium density of water was obtained by computing the pressure for samples with different densities. The equilibrium density was then determined by interpolation as the one at which the internal pressure of the liquid equals the atmospheric pressure. In practice, the pressure was interpolated to P = 0 because the typical errors in the computed pressure are on the order of 10–50 MPa, i.e. larger than the value of the atmospheric pressure. All of our calculations were carried out with the plane-wave first-principles molecular dynamics code Qbox. 60,61 Hamann–Schl¨ uter–Chiang–Vanderbilt (HSCV) pseudopotentials 62,63 were used to represent the core region of oxygen and hydrogen atoms. 64 The 1s electron of hydrogen and 2s2 2p4 electrons of oxygen were treated explicitly. To accurately compute the stress tensor, we used a plane-wave cutoff of 160 Ry (see Supporting Information). Most water simulations were performed at the kinetic energy cutoff of 85 Ry; however, several additional sets of PBE runs were carried out with a cutoff of 160 Ry, to check if MD runs conducted at the same cutoff as the stress-tensor calculations yielded different values of the pressure. We employed 64-water-molecule supercells for all liquid water simulations, 54-molecule cells as a model of cubic ice Ic, and 96-molecule supercells as a model of proton-disordered ice Ih (optimized cell parameters for ices Ic and Ih are reported in the Supporting Information). Water simulations at the experimental density of 1.00 g/ml (cell size of 12.414 ˚ A) were carried out for a set of 8 (4) independent samples with the PBE (PBE0) functional. In addition, we performed a set of one-sample simulations at lower densities of 0.83 (cell size of

6 ACS Paragon Plus Environment

Page 6 of 25

Page 7 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

13.203 ˚ A), 0.80 (13.372 ˚ A), 0.74 (13.732 ˚ A), and 0.71 g/ml (13.944 ˚ A), with both PBE and PBE0 functionals. Starting geometries of all samples were taken from a data set of PBE simulations of water at 400 K equilibrated for at least 10 ps. 65 We performed additional equilibration of all samples for at least ∼ 5–10 ps. Starting geometries of ice Ic and Ih were taken from our previous work. 45 We used deuterium instead of hydrogen in all calculations to take advantage of a larger simulation time step of 0.242 fs. (Yet, the densities reported here use the hydrogen atomic mass of 1.00794.) All simulations were performed in the NVT ensemble at the temperature of 400 K using the Bussi–Donadio–Parrinello thermostat. 66 The thermostat time constant was varied between 25 and 240 fs during equilibration, and was set to 240 fs in production NVT runs. An elevated temperature was chosen in order to obtain good agreement with experiment for the pair-correlation functions and the diffusion coefficients in PBE simulations; 67,68 for consistency, we used the same temperature in our PBE0 runs. A set of simulations was performed in the microcanonical (NVE) ensemble to verify that the choice of the statistical ensemble does not alter the computed pressure. All molecular-dynamics simulations with the PBE0 functional were carried out using a bisection algorithm 69,70 to accelerate the computation of the Hartree–Fock exchange integrals. We used a bisection threshold of 0.02 and a two-level recursion, consistent with our previous work. 71 Pressure calculations with the PBE0 functional were carried out with the full exact-exchange functional (no bisection). Bulk properties of cubic (Ic) and proton-disordered (Ih) ice were computed as follows. Equilibrium densities were determined by fully relaxing atomic coordinates and cell parameters in the ice samples. After optimization, the structures were uniformly expanded and contracted by 0.3% with respect to the equilibrium volume, followed by another optimization of the atomic positions. We then computed pressures in the resulting structures and

7 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 25

Table 1: Equilibrium volume V , density ρ, and bulk modulus B0 of cubic (Ic) and proton-disordered hexagonal (Ih) ice computed using generalized-gradient (PBE) and hybrid (PBE0) density functionals. The structures were optimized at 160 Ry using 54-molecule supercells for ice Ic and 96-molecule supercells for ice Ih with the Qbox code 61 by performing Γ-point-only sampling. Ice, Method Ic, PBE Ih, PBE Ic, PBEa Ih, PBEb Ih, PBEc Ih, PBEd Ih, PBE+TSe

V/H2 O (cm3 /mol) ρ (g/cm3 ) 18.37 0.98 18.40 0.98 18.37 0.98 18.20 0.99 18.45 0.98 18.54 0.97 17.88 1.01

B0 (GPa) 14.7 14.8 14.67 15.6 – – –

Ic, PBE0 Ih, PBE0 Ih, PBE0d Ih, PBE0+TSe

18.74 18.86 18.66 17.99

0.96 0.95 0.97 1.00

13.4 13.0 – –

Expt.

19.29f

0.93

8.33f –12.1g

a

Ref. 45. Ref. 40. c Ref. 38. d Ref. 39. e Ref. 39; “TS” denotes Tkatchenko–Scheffler dispersion corrections. 46 f Density and bulk modulus of ice Ih extrapolated to 0 K from Ref. 41. g Ref. 42. b

8 ACS Paragon Plus Environment

Page 9 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

obtained the bulk moduli from the expression:

B = −V

∆P . ∆V

(2)

As Table 1 shows, our calculations predict the same density for both ices—0.98 g/ml at the PBE level, consistent with Ref. 45. Use of the PBE0 functional lowered the density to 0.95–0.96 g/ml, bringing it closer to the experimental value of 0.93 g/ml. These results are in agreement with the recent study of Santra et al., 39 who found that the density of ice Ih decreased from 0.972 to 0.966 g/ml when using the PBE0 functional instead of PBE. Interestingly, adding Tkatchenko–Scheffler dispersion correction 46 to both PBE and PBE0 functionals increased the density of ice to more than 1 g/ml 39 and worsened the agreement with experiment. The measured bulk modulus of ice is yet uncertain, as indicated by the scattered experimental values; 41,72,73 nevertheless, all experimental data (9–12 GPa) are lower than those predicted by the PBE functional (14–15 GPa). The use of the PBE0 functional reduces the bulk modulus of ice by ∼ 1.5 GPa, bringing it closer to experiment. To model liquid water, we performed a series of simulations at the plane-wave cutoff of 160 Ry. As Figure 1 shows, substantial pressure oscillations in individual runs persist for ∼ 10 ps but level out afterwards. Simulating an ensemble of 8 samples for only 5 ps suffices to get a stable average pressure (red curve in Figure 2); yet most of our runs were at least 20-ps-long to ensure good convergence of the total pressure value. To check if the choice of the statistical ensemble affected the computed pressure, 71 we removed the thermostat after 20 ps of simulation in the canonical ensemble and continued simulating 8 samples of water for another 20 ps at the same cutoff of 160 Ry. The resulting running average pressure P (t) is compared with that obtained in the NVT ensemble in Figure 2. Although converging the pressure takes slightly longer in the NVE ensemble, the two P (t) curves fall within each other’s statistical error bars. The final pressure obtained in the NVT ensemble (0.69 ± 0.01

9 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

GPa) is essentially the same as that computed in the NVE ensemble (0.70 ± 0.02 GPa). 0.9

0.8

P [GPa]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 25

0.7

0.6

0.5

0

5

10

15

20

t [ps]

Figure 1: Running average pressure (P ) for 8 independent samples of 64 water molecules simulated using the PBE functional at a constant temperature of 400 K for ∼ 20 ps. The simulation was carried out at a plane-wave energy cutoff of 160 Ry, and the pressure was computed at each simulation step using the same cutoff. Different colors indicate different samples. It was previously determined that structural properties of water described with normconserving Hamann 62 and HSCV pseudopotentials are converged at the energy cutoff of 85 Ry (Fig. S4). 50,67,68 To verify if this cutoff is also sufficient for our study of the equilibrium density, 29 we performed a series of MD runs at 85 Ry and computed the pressure at 160 Ry for select snapshots, every 19.4 fs (see Supporting Information). The average pressure determined from 8 independent trajectories of water at the experimental density at 85 Ry is 0.57 GPa, compared to 0.69 GPa determined from simulations at 160 Ry. Given the small error bars associated with both values (0.01–0.03 GPa), this difference is significant. The positive value of the pressure obtained at the experimental density indicates that the theoretical density is lower than experiment. To determine the equilibrium density, we chose four additional densities—0.71, 0.74, 0.80, and 0.83 g/ml—and simulated these samples using the PBE functional at 85 Ry, and then at 160 Ry. The results are summarized in Table SI 10 ACS Paragon Plus Environment

Page 11 of 25

0.9

NVT NVE

0.8

P [GPa]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

0.7

0.6

0.5

0

5

10

15

20

t [ps]

Figure 2: Running average pressure (P ) for NVT and NVE simulations of liquid water with the PBE functional. Each curve was computed as an average between 8 independent samples simulated in identical conditions. The red curve is an average of the curves from Figure 1. The shaded area represents an error bar computed as two standard deviations of the mean above and below each curve. of the Supporting Information and in Table 2, and are plotted in Figure 3. We fitted the values of the pressure at different densities to the virial equation of state (see Supporting Information): P = Aρ + Bρ2 + Cρ3 .

(3)

For PBE simulations at 160 Ry we obtained the equilibrium density of 0.81 ± 0.01 g/ml; for simulation at 85 Ry, the value of ρ increased to 0.86 ± 0.02 g/ml. This finding is consistent with the result of McGrath et al., 29 who found that water modeled at 70 Ry is denser than water modeled at 300 Ry by 0.1 g/ml. The difference between our 85-Ry and 160-Ry simulations (0.05 g/ml) is smaller than that of Ref. 29, but it is outside the error bars for each individual density. Overall, our result for the density of PBE water is in very good agreement with recent determinations using NPT 32,37 and NVT 23 simulations techniques. We now turn to the analysis of PBE0 simulation results. We simulated water with the PBE0 functional at 85 Ry energy cutoff, as the difference in densities obtained at 85 and 11 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 2: Equilibrium densities and compressibilities of liquid water computed using the PBE and PBE0 functionals at 400 K with (+D) and without dispersioninteraction corrections (see text) and temperature corrections (+T). The errors given in the Table were computed from errors in the pressure (see Supporting Information). Method

Density, g/ml

Compressibility, Mbar−1

PBE (160 Ry)a

0.81±0.01

69±9

PBEb PBE+Dc PBE+Dc +Td

0.86±0.02 1.02 1.04

44±7 21 19

PBE0b PBE0+Dc PBE0+Dc +Td

0.71±0.02 0.96 0.99

108±35 35 30

Expt.e

1.00

45

a

Both trajectories and pressures computed at 160 Ry. Stress tensor computed at 160 Ry for trajectories obtained with 85 Ry. c Dispersion contribution to pressure estimated using Grimme’s functional with zero-damping 55 for trajectories obtained using PBE and PBE0 functionals. d Kinetic contribution to the stress tensor (N/V ) kB T rescaled from T = 400 to 330 K to evaluate the effect of a lower temperature. e Ref. 43. b

12 ACS Paragon Plus Environment

Page 12 of 25

Page 13 of 25

0.8

PBE0 at 85 Ry PBE at 160 Ry PBE at 85 Ry

0.6 0.4 P [GPa]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

0.2 0.0 -0.2 -0.4

0.7

0.8

0.9

1.0

ρ [g/ml]

Figure 3: Computed pressure of water in ab initio simulations with the PBE and PBE0 functionals at different densities. The equilibrium densities were determined by a fit to the virial equation of state, Eq. (3). The bold horizontal lines represent the ranges of possible equilibrium densities for each simulation protocol, determined by adding or subtracting error bars from the pressure values and fitting to the virial equation (see Supporting Information). 160 Ry is likely smaller than that expected from changing functionals. Table SI summarizes the values of pressure obtained in PBE0 simulations at the experimental density and at four lower densities chosen to be the same as in PBE simulations. Using the values reported in Table SI, the equilibrium density of PBE0 water was determined to be 0.71 ± 0.02 g/ml. Taking into account that the density obtained from 160 Ry trajectories is expected to be lower by ∼ 0.05 g/ml, the equilibrium density of PBE0 water could be as low as 0.65 g/ml. Our findings are summarized in Figure 3, which compares the pressure–density curves determined by fitting the simulation data to the virial equation of state, Eq. (3). The point at which each curve crosses the P = 0 line corresponds to the equilibrium density, with the slope at crossing being related to the isothermal compressibility κT = (1/ρ0 ) (∂ρ/∂P )T . The compressibility of water simulated using the PBE functional at 85 Ry is 44 ± 7 Mbar−1 , close to the experimental measurement of 45 Mbar−1 . 43 Our result is consistent with the value of 32 ± 7 Mbar−1 , determined by fitting the PBE simulation data of Wang et at. 23

13 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

to Eq. (3); the discrepancy could be attributed to a lower temperature T = 360 K chosen in Ref. 23 and/or to the use of different basis sets. The good agreement of these values with experiment appears to be fortuitous as the compressibility of water simulated using PBE at 160 Ry cutoff is higher, 69 ± 9 Mbar−1 , overestimating the experimental value by ∼ 50%. The compressibility of water described at the PBE0 level (108 ± 35 Mbar−1 ) is also 3.6

Hydrogen bonds / H2O molecule

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3.2

2.8

PBE PBE0 2.4

0.7

0.8

0.9

1.0

ρ [g/ml]

Figure 4: Number of hydrogen bonds as a function of density (ρ) from PBE and PBE0 simulations carried out in this work. The number of hydrogen bonds was determined using the geometrical criterion of Ref. 74: Water molecules were considered hydrogen-bonded if the distance between oxygen atoms was smaller than 3.35 ˚ A and the angle O· · · O—H was within 30◦ . Error bars were computed as two standard deviations of the mean for ρ = 1.0 g/ml and two standard deviations (taken from simulations at ρ = 1.0 g/ml with the same functional) for all other densities. higher than in the experiment. We note that the compressibility determined using the PBE0 functional has a large statistical uncertainty because the absolute value of the pressure in the low-density region of the P –ρ plot is comparable to the magnitude of error bars. Our analysis confirms previous reports that PBE underestimates the density of liquid water. 23,32,37 The new finding that the PBE0 functional does not improve upon the PBE results might appear surprising, given that hybrid functionals improve the description of the structural 50,51 and electronic properties 71,75 of liquid water. It was previously concluded that 14 ACS Paragon Plus Environment

Page 14 of 25

Page 15 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

one of the reasons why GGAs underestimate the density of water is the lack of dispersion interactions in these functionals. 76 PBE overestimates the hydrogen bond strength, leading to an overstructured liquid and partially compensating for missing dispersion interactions. Compared to the PBE functional, PBE0 has lower self-interaction error, weakening the hydrogen bonds between water molecules (Figure 4). Hence, the use of the PBE0 functional does not lead to cancellation of errors, and water described with this approximation has a lower density. Although our results would seem to indicate that the PBE0 functional provides a less accurate description of the density of liquid water than PBE, we note that adding dispersion corrections 35,44,46,55,77–81 to the two functionals may bring the PBE0 results in better agreement with experiments. As a first-order estimate of the effect of dispersion interactions, we computed the contribution to the stress tensor arising from the Grimme dispersion functional 55 with zero-damping for PBE and PBE0 trajectories. For each trajectory, we used the specific set of parameters tailored to the corresponding functional, as implemented in the dftd3 code. 82 Although such post-MD scheme is not equivalent to simulations with the dispersion-corrected functional, it is nevertheless a reasonable first-order approximation to the effect of dispersion interactions. After adding corrections to pressure values computed along PBE and PBE0 trajectories at 85 Ry, we found an equilibrium density of 1.02 g/ml at the PBE+D level and 0.96 g/ml at the PBE0+D level of theory (Table 2). Compressibilities from dispersion-corrected PBE and PBE0 functionals are 21 and 35 Mbar−1 respectively, with PBE0+D results now in much better agreement with experiment. The density of water at the PBE0+D level is still slightly further away from the experiment than that at the PBE+D level. A possible reason could be the effect of the high simulation temperature.1 As mentioned above, simulations in our work were performed at T = 400 K using both PBE and PBE0 functionals. This temperature, shown to be optimal for the PBE functional in order to obtain good agreement with measured structural 1

Experimental density of water at 373 K and atmospheric pressure is 0.958 g/ml, while that at 400 K and 0.25 MPa is 0.937 g/ml. 43

15 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

properties, 67,68 is probably too high to obtain the same agreement with measured structural properties, in the case of PBE0.2 To estimate the effect of a lower simulation temperature, we rescaled the kinetic contribution to the pressure, (N/V ) kB T . Subtracting the difference of this term between 400 and 330 K from the pressures reported in “DFT+D” columns of Table SI of Supporting Information, we estimated that the density of PBE0+D water would become almost exact (0.99 g/ml), while that of PBE+D water would further deviate from experiment (1.04 g/ml). Compressibilities of liquid water would decrease to 19 and 30 Mbar−1 at the PBE+D+T and PBE0+D+T levels, respectively. Overall, our dispersion- and temperature-corrected compressibilities of water turned out to be close to compressibilities determined at the vdW-DF (14.8 Mbar−1 ) and VV10 (25.8 Mbar−1 ) levels of theory. 33 In summary, we found that the PBE0 hybrid functional lowers the density of ice (by 3%) and water (by 17%) compared to the gradient-corrected PBE functional. This leads to a better agreement with measurements for ice but a worse one for liquid water. Adding Grimme’s dispersion corrections to PBE and PBE0 functionals improves the equilibrium density and compressibility of water; however, as shown by Santra et al., 39 the agreement for ice might get worse. These results clearly call for development of functionals inclusive of van der Waals interactions capable of describing the hydrogen bonded configurations present in both water and ice at the same level of accuracy. Use of hybrid functionals together with dispersion-corrected functionals could lead to improved description of hydrogen bonding, electronic structure, and bulk properties of solid and liquid water, and provide an all-around first-principles model of water. 2

PBE0 simulations in our previous work 50,71 were carried out at 380 K. We believe the choice of higher temperature did not affect the main message of the present study—that the density of water is lower at the PBE0 than at the PBE level of theory—because the PBE and PBE0 simulations were performed in the same conditions.

16 ACS Paragon Plus Environment

Page 16 of 25

Page 17 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

Acknowledgement This work was supported by DOE/BES (Grant No. DE-SC0008938). An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. A.P.G. acknowledges helpful discussions with Dr. Ikutaro Hamada.

Supporting Information Available This material is available free of charge via the Internet at http://pubs.acs.org/.

References (1) Ball, P. Water: Water—an enduring mystery. Nature 2008, 452, 291–292. (2) Ball, P. In Astrochemistry and Astrobiology; Smith, I. W. M., Cockell, C. S., Leach, S., Eds.; Springer-Verlag: Berlin Heidelberg, 2013; pp 169–210. (3) Nilsson, A.; Pettersson, L. G. M. Perspective on the structure of liquid water. Chem. Phys. 2011, 389, 1–34. (4) Brovchenko, I.; Oleinikova, A. Multiple phases of liquid water. ChemPhysChem 2008, 9, 2660–2675. (5) Fischer, H. E.; Barnes, A. C.; Salmon, P. S. Neutron and x-ray diffraction studies of liquids and glasses. Rep. Prog. Phys. 2006, 69, 233–299. (6) Soper, A. K. Joint structure refinement of x-ray and neutron diffraction data on dis-

17 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ordered materials: Application to liquid water. J. Phys.: Condens. Matter 2007, 19, 335206. (7) Soper, A.; Benmore, C. J. Quantum differences between heavy and light water. Phys. Rev. Lett. 2008, 101, 065502. (8) Skinner, L. B.; Huang, C.; Schlesinger, D.; Pettersson, L. G. M.; Nilsson, A.; Benmore, C. J. Benchmark oxygen–oxygen pair-distribution function of ambient water from x-ray diffraction measurements with a wide Q-range. J. Chem. Phys. 2013, 138, 074506. (9) Skinner, L. B.; Benmore, C. J.; Neuefeind, J. C.; Parise, J. B. The structure of water around the compressibility minimum. J. Chem. Phys. 2014, 141, 214507. (10) Ichiye, T. In Advances in Chemical Physics; Rice, S. A., Dinner, A. R., Eds.; John Wiley & Sons, Inc., 2014; Vol. 155; pp 161–199. (11) Vega, C.; Abascal, J. L. F. Simulating water with rigid non-polarizable models: A general perspective. Phys. Chem. Chem. Phys. 2011, 13, 19663–19688. (12) Guillot, B. A reappraisal of what we have learnt during three decades of computer simulations on water. J. Mol. Liq. 2002, 101, 219–260. (13) Car, R.; Parrinello, M. Unified approach for molecular dynamics and density-functional theory. Phys. Rev. Lett. 1985, 55, 2471–2474. (14) Galli, G.; Pasquarello, A. In Computer Simulation in Chemical Physics; Allen, M. P., Tildesley, D. J., Eds.; Springer Netherlands, 1993; pp 261–313. (15) Car, R.; de Angelis, F.; Giannozzi, P.; Marzari, N. In Handbook of Materials Modeling; Yip, S., Ed.; Springer Netherlands, 2005; pp 59–76. (16) Marx, D.; Hutter, J. Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods; Cambridge University Press: New York, 2009. 18 ACS Paragon Plus Environment

Page 18 of 25

Page 19 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

(17) Scuseria, G. E.; Staroverov, V. N. In Theory and Applications of Computational Chemistry: The First Forty Years; Dykstra, C. E., Frenking, G., Kim, K. S., Scuseria, G. E., Eds.; Elsevier: Amsterdam, 2005. (18) Becke, A. D. Perspective: Fifty years of density-functional theory in chemical physics. J. Chem. Phys. 2014, 140, 18A301. (19) Perdew, J. P.; Kurth, S. In A Primer in Density Functional Theory; Fiolhais, C., Nogueira, F., Marques, M., Eds.; Springer: Berlin, 2003. (20) Perdew, J. P.; Zunger, A. Self-interaction correction to density-functional approximations for many-electron systems. Phys. Rev. B 1981, 23, 5048–5079. (21) Cohen, A. J.; Mori-S´anchez, P.; Yang, W. Insights into current limitations of density functional theory. Science 2008, 321, 792–794. (22) Cohen, A. J.; Mori-S´anchez, P.; Yang, W. Challenges for density functional theory. Chem. Rev. 2012, 112, 289–320. (23) Wang, J.; Rom´an-P´erez, G.; Soler, J. M.; Artacho, E.; Fern´andez-Serra, M.-V. Density, structure, and dynamics of water: The effect of van der Waals interactions. J. Chem. Phys. 2011, 134, 024516. (24) Bankura, A.; Karmakar, A.; Carnevale, V.; Chandra, A.; Klein, M. L. Structure, dynamics, and spectral diffusion of water from first-principles molecular dynamics. J. Phys. Chem. C 2014, 118, 29401–29411. (25) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (26) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple [Phys. Rev. Lett. 77, 3865 (1996)]. Phys. Rev. Lett. 1997, 78, 1396.

19 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(27) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 1988, 38, 3098–3100. (28) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle–Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785–789. (29) McGrath, M. J.; Siepmann, J. I.; Kuo, I.-F. W.; Mundy, C. J.; VandeVondele, J.; Hutter, J.; Mohamed, F.; Krack, M. Isobaric–isothermal Monte Carlo simulations from first principles: Application to liquid water at ambient conditions. ChemPhysChem 2005, 6, 1894–1901. (30) McGrath, M. J.; Siepmann, J. I.; Kuo, I.-F. W.; Mundy, C. J.; VandeVondele, J.; Hutter, J.; Mohamed, F.; Krack, M. Simulating fluid-phase equilibria of water from first principles. J. Phys. Chem. A 2005, 110, 640–646. (31) McGrath, M. J.; Siepmann, J. I.; Kuo, I.-F. W.; Mundy, C. J. Vapor–liquid equilibria of water from first principles: Comparison of density functionals and basis sets. Mol. Phys. 2006, 104, 3619–3626. (32) Schmidt, J.; VandeVondele, J.; Kuo, I.-F. W.; Sebastiani, D.; Siepmann, J. I.; Hutter, J.; Mundy, C. J. Isobaric–isothermal molecular dynamics simulations utilizing density functional theory: An assessment of the structure and density of water at near-ambient conditions. J. Phys. Chem. B 2009, 113, 11959–11964. (33) Corsetti, F.; Artacho, E.; Soler, J. M.; Alexandre, S. S.; Fern´andez-Serra, M.-V. Room temperature compressibility and diffusivity of liquid water from first principles. J. Chem. Phys. 2013, 139, 194502. (34) Dion, M.; Rydberg, H.; Schr¨oder, E.; Langreth, D. C.; Lundqvist, B. I. Van der Waals density functional for general geometries. Phys. Rev. Lett. 2004, 92, 246401.

20 ACS Paragon Plus Environment

Page 20 of 25

Page 21 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

(35) Vydrov, O. A.; Van Voorhis, T. Nonlocal van der Waals density functional: The simpler the better. J. Chem. Phys. 2010, 133, 244103. (36) Zhang, Y.; Yang, W. Comment on “Generalized gradient approximation made simple”. Phys. Rev. Lett. 1998, 80, 890. (37) Miceli, G.; de Gironcoli, S.; Pasquarello, A. Isobaric first-principles molecular dynamics of liquid water with nonlocal van der Waals interactions. J. Chem. Phys. 2015, 142, 034501. (38) Feibelman, P. J. Lattice match in density functional calculations: Ice Ih vs. β-AgI. Phys. Chem. Chem. Phys. 2008, 10, 4688–4691. (39) Santra, B.; Klimeˇs, J.; Tkatchenko, A.; Alf`e, D.; Slater, B.; Michaelides, A.; Car, R.; Scheffler, M. On the accuracy of van der Waals inclusive density-functional theory exchange-correlation functionals for ice at ambient and high pressures. J. Chem. Phys. 2013, 139, 154702. (40) Hamada, I. A van der Waals density functional study of ice Ih. J. Chem. Phys. 2010, 133, 214503. (41) Whalley, E. Energies of the phases of ice at zero temperature and pressure. J. Chem. Phys. 1984, 81, 4087–4092. (42) Gammon, P. H.; Kiefte, H.; Clouter, M. J. Elastic constants of ice samples by Brillouin spectroscopy. J. Phys. Chem. 1983, 87, 4025–4029. (43) Wagner, W.; Pruß, A. The IAPWS formulation 1995 for the thermodynamic properties of ordinary water substance for general and scientific use. J. Phys. Chem. Ref. Data 2002, 31, 387–535. ´ D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Higher-accuracy (44) Lee, K.; Murray, E. van der Waals density functional. Phys. Rev. B 2010, 82, 081101. 21 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

´ Galli, G. Dispersion interactions and vibrational effects in ice as a function (45) Murray, E.; of pressure: A first principles study. Phys. Rev. Lett. 2012, 108, 105502. (46) Tkatchenko, A.; Scheffler, M. Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data. Phys. Rev. Lett. 2009, 102, 073005. (47) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648–5652. (48) Perdew, J. P.; Ernzerhof, M.; Burke, K. Rationale for mixing exact exchange with density functional approximations. J. Chem. Phys. 1996, 105, 9982–9985. (49) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158–6170. (50) Zhang, C.; Donadio, D.; Gygi, F.; Galli, G. First principles simulations of the infrared spectrum of liquid water using hybrid density functionals. J. Chem. Theory Comput. 2011, 7, 1443–1449. (51) DiStasio Jr., R. A.; Santra, B.; Li, Z.; Wu, X.; Car, R. The individual and collective effects of exact exchange and dispersion interactions on the ab initio structure of liquid water. J. Chem. Phys. 2014, 141, 084502. (52) Del Ben, M.; Sch¨onherr, M.; Hutter, J.; VandeVondele, J. Bulk liquid water at ambient temperature and pressure from MP2 theory. J. Phys. Chem. Lett. 2013, 4, 3753–3759. (53) Del Ben, M.; Sch¨onherr, M.; Hutter, J.; VandeVondele, J. Correction to “Bulk liquid water at ambient temperature and pressure from MP2 theory”. J. Phys. Chem. Lett. 2014, 5, 3066–3067. (54) Guidon, M.; Hutter, J.; VandeVondele, J. Auxiliary density matrix methods for Hartree–Fock exchange calculations. J. Chem. Theory Comput. 2010, 6, 2348–2364. 22 ACS Paragon Plus Environment

Page 22 of 25

Page 23 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

(55) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H–Pu. J. Chem. Phys. 2010, 132, 154104. (56) Martin, R. M. Electronic Structure: Basic Theory and Practical Methods; Cambridge University Press: Cambridge, 2004. (57) Focher, P.; Chiarotti, G. L.; Bernasconi, M.; Tosatti, E.; Parrinello, M. Structural phase transformations via first-principles simulation. Europhys. Lett. 1994, 26, 345–351. (58) Bernasconi, M.; Chiarotti, G. L.; Focher, P.; Scandolo, S.; Tosatti, E.; Parrinello, M. First-principle-constant pressure molecular dynamics. J. Phys. Chem. Solids 1995, 56, 501–505. (59) Focher, P. First-principle studies of structural phase transformations. Ph.D. thesis, Scuola Internazionale Superiore di Studi Avanzati (SISSA), Trieste (1994). (60) Gygi, F. Architecture of Qbox: A scalable first-principles molecular dynamics code. IBM J. Res. Dev. 2008, 52, 137–144. (61) Qbox code: http://eslab.ucdavis.edu/software/qbox/ (retrieved June 6, 2015). (62) Hamann, D.; Schl¨ uter, M.; Chiang, C. Norm-conserving pseudopotentials. Phys. Rev. Lett. 1979, 43, 1494–1497. (63) Vanderbilt, D. Optimally smooth norm-conserving pseudopotentials. Phys. Rev. B 1985, 32, 8412–8415. (64) HSCV pseudopotential table:

http://fpmd.ucdavis.edu/potentials/index.htm

(retrieved June 6, 2015). (65) Water PBE400 dataset: http://www.quantum-simulation.org (retrieved June 6, 2015). 23 ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(66) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. (67) Grossman, J. C.; Schwegler, E.; Draeger, E. W.; Gygi, F.; Galli, G. Towards an assessment of the accuracy of density functional theory for first principles simulations of water. J. Chem. Phys. 2004, 120, 300–311. (68) Schwegler, E.; Grossman, J. C.; Gygi, F.; Galli, G. Towards an assessment of the accuracy of density functional theory for first principles simulations of water. II. J. Chem. Phys. 2004, 121, 5400–5409. (69) Gygi, F. Compact representations of Kohn–Sham invariant subspaces. Phys. Rev. Lett. 2009, 102, 166406. (70) Gygi, F.; Duchemin, I. Efficient computation of Hartree–Fock exchange using recursive subspace bisection. J. Chem. Theory Comput. 2013, 9, 582–587. (71) Gaiduk, A. P.; Zhang, C.; Gygi, F.; Galli, G. Structural and electronic properties of aqueous NaCl solutions from ab initio molecular dynamics simulations with hybrid density functionals. Chem. Phys. Lett. 2014, 604, 89–96. (72) Gold, L. W. Some observations on the dependence of strain on stress for ice. Can. J. Phys. 1958, 36, 1265–1275. (73) Yoshimura, Y.; Stewart, S. T.; Somayazulu, M.; Mao, H. K.; Hemley, R. J. Highpressure x-ray diffraction and Raman spectroscopy of ice VIII. J. Chem. Phys. 2006, 124, 024502. (74) Luzar, A.; Chandler, D. Structure and hydrogen bond dynamics of water–dimethyl sulfoxide mixtures by computer simulations. J. Chem. Phys. 1993, 98, 8160–8173. (75) Zhang, C.; Pham, T. A.; Gygi, F.; Galli, G. Electronic structure of the solvated chloride anion from first principles molecular dynamics. J. Chem. Phys. 2013, 138, 181102. 24 ACS Paragon Plus Environment

Page 24 of 25

Page 25 of 25

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

(76) Jonchiere, R.; Seitsonen, A. P.; Ferlat, G.; Saitta, A. M.; Vuilleumier, R. Van der Waals effects in ab initio water at ambient and supercritical conditions. J. Chem. Phys. 2011, 135, 154503. (77) Grimme, S. Accurate description of van der Waals complexes by density functional theory including empirical corrections. J. Comput. Chem. 2004, 25, 1463–1473. (78) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787–1799. (79) Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456–1465. (80) Becke, A. D.; Johnson, E. R. A density-functional model of the dispersion interaction. J. Chem. Phys. 2005, 123, 154101. (81) Hamada, I. Van der Waals density functional made accurate. Phys. Rev. B 2014, 89, 121103. (82) dftd3 code, v.3.1 rev. 0: http://www.thch.uni-bonn.de/tc/index.php?section= downloads&subsection=getd3&lang=english (retrieved June 6, 2015).

25 ACS Paragon Plus Environment