First-Principles Simulations of Liquid Water Using a Dielectric

geometrical criteria, the number of hydrogen bonds per water molecule in SCAN simulations is ∼ 3.6 Å, based on the data ..... 4, so we did not incl...
0 downloads 2 Views 3MB Size
Subscriber access provided by Kaohsiung Medical University

Chemical and Dynamical Processes in Solution; Polymers, Glasses, and Soft Matter

First-Principles Simulations of Liquid Water Using a Dielectric-Dependent Hybrid Functional Alex P. Gaiduk, Jeffrey A. Gustafson, Francois Gygi, and Giulia Galli J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/acs.jpclett.8b01017 • Publication Date (Web): 16 May 2018 Downloaded from http://pubs.acs.org on May 16, 2018

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 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 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.

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 21 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

First-Principles Simulations of Liquid Water Using a Dielectric-Dependent Hybrid Functional Alex P. Gaiduk,†,‡ Jeffrey Gustafson,¶ Francois Gygi,§ and Giulia Galli∗,†,¶,‡ †Institute for Molecular Engineering, The University of Chicago, Chicago, Illinois 60637, United States ‡Materials Science Division, Argonne National Laboratory, Argonne, Illinois 60439, United States ¶Department of Chemistry, The University of Chicago, Chicago, Illinois 60637, United States §Department of Computer Science, University of California, Davis, California 95616, United States E-mail: [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 carried out first-principles simulations of liquid water at ambient conditions using a dielectric-dependent hybrid functional, where the fraction of exact exchange is set equal to the inverse of the high-frequency dielectric constant of the liquid. We found excellent agreement with experiment for the oxygen-oxygen partial correlation function at the experimental equilibrium density and 311 ± 3 K. Other structural and dynamical properties such as the diffusion coefficient, molecular dipole moments, and vibrational spectra, are also in good agreement with experiment. Our results, together with previous findings on electronic properties of the liquid with the same functional, show that the dielectric-dependent hybrid functional accurately describes both the structural and electronic properties of liquid water.

Graphical TOC Entry

Keywords Water structure; simulations; ab initio molecular dynamics; hybrid functionals; dielectricdependent functionals

2

ACS Paragon Plus Environment

Page 2 of 21

Page 3 of 21 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

First-principles simulations of water are a challenging problem. It has so far been difficult to find a density functional or hybrid functional that can provide an equally accurate description of the different types of bonds present in water: intramolecular covalent bonds, intermolecular hydrogen bonds and van der Waals interactions, as well as the electronic and dielectric properties of the liquid. 1 However, in recent years substantial progress has been reported in the description of the structural properties of liquid water. For example, DiStasio et al. 2 assessed the performance of increasingly complex density functionals to describe water, including hybrid functionals with dispersion corrections, 3 and they concluded that both exact exchange and dispersion corrections are necessary to faithfully reproduce the structure and the dynamical properties of liquid water. Chen et al. 4 used the SCAN meta-GGA functional 5 to simulate water and reported results comparable in accuracy to those obtained with the computational protocol of DiStasio et al. 2 In addition, SCAN correctly predicts that ice floats on water—an important property that generalized-gradient approximations (e.g. PBE 6–8 ) and several hybrid functionals 9,10 cannot reproduce. The simulations of Refs. 2 and 4 were carried out at a temperature of 330 K—closer to ambient conditions than the temperature of 400 K required for a satisfactory agreement of PBE results with experiments. 9,11–15 Several advances have also been reported in the description of the quantum nature of protons in water. 16,17 Recent work from Markland et al., 18–20 for example, explored the explicit inclusion of nuclear quantum effects (NQEs) using a ring-polymer contraction method combined with a multiple time-step scheme. 18 Using this computational protocol, the authors performed simulations of several hundreds of picoseconds with various gradient-corrected and hybrid density functionals, and determined how nuclear quantum effects affect the structure, dynamics, and vibrational properties of water. 20 One of their findings (corroborated by other researchers 21–24 ) is that the inclusion of NQEs red-shifts the O–H stretching mode of the infrared spectrum of water—this result will be important in the context of our work. Recently, we studied the electronic properties of liquid water 25 and solutions of salts 26 and

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

determined that among several generalized-gradient approximations and hybrid functionals, dielectric-dependent hybrids (DDH) 27,28 predict the most accurate photoelectron spectra and energies of the unoccupied states of water. In this work we use first-principles molecular dynamics to simulate water using a DDH functional 27 with a fraction of exact exchange equal to the inverse of the high-frequency dielectric constant of the liquid, and assess its performance for structural, dielectric, and vibrational properties. We find that this DDH functional describes the structure and dielectric properties of the liquid accurately at room temperature, while the position of the O–D stretch peak in the infrared (IR) spectrum appears to be blue-shifted, compared to experiment. We argue that when adding nuclear quantum effects, the O–O partial correlation function would only be negligibly affected, while the position of the O–D stretching band would be in better agreement with experiments. We start by describing the structural properties of the liquid. Figure 1 compares the oxygen–oxygen gO–O and oxygen–hydrogen radial distribution functions computed using a DDH functional at 311 ± 3 K, with experiment and with results obtained using the PBE (400 K) and SCAN (330 K) 4 functionals. As noted by many authors, 1,2,4,11,12,15,29–31 even at 400 K, the PBE gO–O is appreciably overstructured compared to experiment, with the position of the first peak centered at smaller distances. Both DDH and SCAN functionals yield a gO–O in much closer agreement with experiments, although the simulations with DDH were done at a temperature of 311 ± 3 K, closer to ambient conditions. A comparison with the results of MB-pol 32–34 model is given in the Supporting Information. Both DDH (311 K) and SCAN (330 K) functionals also predict the oxygen–hydrogen radial distribution function in closer agreement to experiment than PBE (400 K), as shown in the bottom panel of Fig. 1. The DDH functional predicts an intramolecular O–H bond length of ∼ 0.95 ˚ A, while SCAN and PBE yield 0.98 and 1.0 ˚ A, respectively. The inclusion of nuclear quantum effects is expected to lengthen the water intramolecular bond length by 0.02 ˚ A 25 with respect to that obtained with classical simulations. Hence simulations inclusive

4

ACS Paragon Plus Environment

Page 4 of 21

Page 5 of 21

gO-O(r)

3

PBE SCAN DDH

2

1

0

2

4

6

8

3

gO-H(r)

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

2

1

0

0

2

4

6

8

r [Å]

Figure 1: Oxygen–oxygen (gO–O ) and oxygen–hydrogen (gO–H ) radial distribution functions of liquid water as computed using a dielectric-dependent hybrid functional at 311 ± 3 K, the PBE functional at 400 K, and the SCAN functional at 330 K (from Ref. 4). The experimental O–O radial distribution function, shown by the solid gray shaded area, is from Skinner et al., 35 and the O–H distribution function for light water is from Soper. 36 of NQEs with the DDH-functional would yield a bond length of approximately 0.97 ˚ A, in closer agreement with the experimental value of 0.99 ˚ A inferred from neutron diffraction data. 37 At the same time, inclusion of NQEs is not expected to significantly change the O–O radial distribution. 20 Using a geometrical definition of the hydrogen bond 38,39 with an O–O cutoff distance of 3.5 ˚ A and O· · · OD cutoff angle of 30◦ , we estimated the number of hydrogen bonds in the simulations with the DDH functional to be 3.39 per water molecule. Using the same geometrical criteria, the number of hydrogen bonds per water molecule in SCAN simulations is ∼ 3.6 ˚ A, based on the data of Fig. 3A of Ref. 4. The number of hydrogen bonds per molecule is 3.54/D2 O for the PBE trajectory, lower than ∼ 3.75/D2 O reported in Ref. 4 5

ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

θ

P(q)

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 6 of 21

PBE DDH

0.0

0.2

0.4

0.6

0.8

1.0

q

Figure 2: Distribution of the tetrahedrality index q computed using Eq. (1) for liquid water simulated with PBE (400 K) and DDH (311 ± 3 K) functionals. The inset illustrates the definition of the angle θij in Eq. (1). due to Chen et al. performing simulations at a lower temperature of 330 K. Experimentally, the number of hydrogen bonds per water molecule is not known but estimates based on Xray absorption, neutron scattering, and NMR 40 suggest values between 3 and 4 at ambient conditions, although some studies reported much lower number of hydrogen bonds per H2 O of 2.2. 41 Figure 2 shows the distribution P (q) of the tetrahedrality index defined as 42,43 * q =1−

2 + 3 4  3X X 1 cos θij + , 8 i=1 j=i+1 3

(1)

where θij is the angle formed between the oxygen atom of a given water molecule and the oxygen atoms i and j of its four nearest neighbors (Figure 2 inset). The parameter q = 1 corresponds to a perfect tetrahedral arrangement and q = 0 to a completely disordered system. 42 In PBE simulations, P (q) has a main peak at ∼ 0.8 and a secondary peak at ∼ 0.5. As pointed out by Guardia et al., 43 the skewed shape of this distribution arises from most

6

ACS Paragon Plus Environment

Page 7 of 21

0.020

POOO(θ) sin θ [deg-1]

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

Expt. PBE SCAN DDH

0.015

0.010

0.005

0.000

40

60

80

100

120

140

160

180

θ [deg]

Figure 3: The oxygen–oxygen–oxygen angular distribution functions POOO (θ) computed from our simulations with PBE (400 K), SCAN (330 K), and DDH (311 ± 3 K) functionals, and compared with the empirical potential structural refinement of X-ray/neutron scattering data Rπ for D2 O. 44 Similar to Ref. 2, the POOO (θ) functions were normalized to 0 dθ POOO (θ) sin θ = 1. water molecules participating in 4 hydrogen bonds. The simulations performed with the DDH functional exhibit a less pronounced main peak and a more pronounced secondary peak, consistent with fewer hydrogen bonds per water molecule. The average tetrahedrality index obtained here with the PBE functional at 400 K (0.68 ± 0.02) is similar to the result of Chen et al. obtained with the SCAN functional at 330 K (0.68 ± 0.18). Our DDH result (0.63 ± 0.01) is similar to the estimate of DiStasio et al. obtained using PBE0 and dispersion corrections (0.64); 2 the estimate obtained from empirical fits to the experimental data of Soper and Benmore is 0.576. 44 Additional information on the local structure of water may be obtained by computing the oxygen–oxygen–oxygen triplet angular distribution function POOO (θ). 2 The results, displayed in Figure 3, show that the PBE POOO (θ) function peaks near 105◦ , close to the perfect tetrahedral angle of 109.47◦ . The small shoulder near 55◦ , indicative of distorted arrangement of hydrogen bonds in the first coordination shell, is underestimated compared to the empirical potential structural refinement of X-ray/neutron scattering data. 44 Use of the DDH (311 K) 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

and SCAN (330 K) functionals improves the agreement with experiment; the results of simulations with these functionals agree with each other, up to the statistical uncertainty of these simulations The diffusion coefficient of liquid water described using the DDH functional at 311 K is (2.52 ± 0.44) · 10−5 cm2 /s, which agrees well with the experimental result of 2.58 · 10−5 cm2 /s.1 The diffusion coefficient computed with the SCAN functional (1.90 ± 0.25) · 10−5 cm2 /s at 330 K is slightly lower than the DDH result (although still within wide 13,46 statistical error bars), indicating a slightly less diffusive liquid, consistent with the larger number of hydrogen bonds per water molecule and slightly higher tetrahedral order in water simulated with SCAN. To estimate the equilibrium density of water described with the DDH functional, we computed the stress tensor along our molecular-dynamics trajectories. The internal pressure, determined as an average of the trace of the stress tensor, is equal to 0.3 GPa. Assuming that the compressibility of DDH water is close to the experimental value 47 of 45 Mbar−1 (for example, the compressibility of water in the PBE simulations is 44±7 Mbar−1 ), 9 the equilibrium density of water obtained using the DDH functional is estimated to be ∼ 0.9 g/mL. The density of water described with SCAN (1.05 g/mL) is slightly higher than the experimental result. (Note that accurate estimates of the equilibrium density of water predicted by the DDH functional would require carrying out NPT simulations 4 or computing the pressure for several different water densities; 9 these simulations are beyond the scope of the present work, and here we just report a preliminary estimate of the liquid density.) Next, we discuss the dielectric properties of liquid water. The upper panel of Figure 4 compares distributions of distances between oxygen atoms and the maximally localized Wannier functions (MLWF) centers (dO–MLWF ) predicted by the different functionals. The peak at larger distances corresponds to the MLWF centers located along the O–D bonds (bonding electron pairs) while that at shorter distances corresponds to the oxygens’ lone pairs. The 1 Obtained by linear interpolation of self-diffusion coefficients in heavy water 45 at 278 K (1.015 · 10−5 cm2 /s), 298 K (1.872 · 10−5 cm2 /s), and 318 K (2.979 · 10−5 cm2 /s), to 311 K.

8

ACS Paragon Plus Environment

Page 8 of 21

Page 9 of 21

Probability [arb. units]

PBE SCAN DDH

0.3

0.4

0.5

dO-MLWF [Å]

PBE SCAN DDH

Probability [arb. units]

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

1.5

2.0

2.5

3.0

3.5

4.0

4.5

Dipole moment [D]

Figure 4: (Top panel) Distribution of distances between oxygen atom and centers of maximally localized Wannier functions (MLWF) in water simulated with PBE, SCAN (ref. 4), and DDH functionals. (Bottom panel) Distribution of water dipole moments in the same simulations. Vertical bars mark average dipole moments of water molecules obtained with the DDH (2.75 D), SCAN (2.97 D), and PBE (3.05 D) functionals. The experimental dipole moment of water, 2.9 D with the statistical uncertainty of 0.6 D, 48 is shown by the gray shaded region. bond-pair peak of DDH is located at ∼ 0.49 ˚ A, while that of PBE is at 0.5 ˚ A (SCAN results 4 are between DDH and PBE), consistent with the difference of OD bond lengths found in the three trajectories (Figure 1). The lone-pair peak predicted by the DDH functional is at a shorter distance from the

9

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

oxygen atom than of PBE. Also, the PBE lone-pair peak has a very pronounced shoulder, while the corresponding peak in the DDH simulation is symmetric. As pointed out in Ref. 49, the shoulder on the lone-pair peak is a fingerprint of broken hydrogen bonds. Absence of this shoulder indicates that the lone pairs participating in hydrogen bonding are more localized with hybrid functionals than with PBE, and less available for hydrogen bonding. Again, the position and shape of the lone-pair peak predicted by the SCAN functional 4 is intermediate between the PBE and DDH results. The bottom panel of Figure 4 shows the distributions of the dipole moments in the liquid peaked at 3.05 D (PBE) and 2.75 D (DDH). The average dipole moment of water described with the SCAN functional is 2.97 D. 4 The results from these three functionals are consistent with the X-ray and neutron scattering results of 2.9 D, within the large experimental error of 0.6 D. 48 Figure 5 compares the infrared spectra computed with the PBE and DDH functionals. In the spectrum, the peak near 2500 cm−1 corresponds to intramolecular O–D bond vibrations. Its position is significantly affected by the electronic-structure method used—for the PBE functional, the peak is red-shifted by 150 cm−1 with respect to the experiment, while for the DDH functional, it is blue-shifted by about 250 cm−1 . However, as pointed out by Marsalek and Markland 20 (as well as several other workers 21–23 ) this agreement with experiment would likely change after taking into account nuclear quantum effects. Ring-polymer and centroid molecular-dynamics simulations of water yield an O–H stretch peak red-shifted by about 175 to 200 cm−1 compared to classical molecular dynamics with ab initio potentials, 20–24 at least in the case of light water. (Note that these methods of accounting for nuclear quantum effects may produce spurious frequency shifts at low temperatures 52,53 but not at 311 K used in our work.) The magnitude of the red shift of the vibrational frequency of O–D in heavy water is approximately half of that in light water. 22,50 As such, upon inclusion of NQEs, the O–D peak would be red-shifted by ∼ 80–90 cm−1 —this change would improve the agreement of DDH results with experiments.

10

ACS Paragon Plus Environment

Page 10 of 21

Page 11 of 21

Expt. PBE DDH

Intensity [arb. units]

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

500

1000

1500

2000

2500

3000

Wavenumber [cm-1]

Figure 5: Computed and experimental infrared spectra (IR) of liquid water. The calculations were carried out with the PBE functional at 400 K, and with a DDH functional at 311 K. Red dashed curve shows the portion of the DDH curve between 2100 and 3200 cm−1 redshifted by 85 cm−1 to indicate a possible correction to the position the O–D stretch peak due to proton quantum effects (the magnitude of the shift for O–D stretching frequency is unknown, but it is estimated to be approximately half 50 of the O–H stretching frequency shift of 175 23 –200 20 cm−1 ). The experimental spectrum at room temperature is from Ref. 51. No analysis of the vibrational properties of water was reported in Ref. 4, so we did not include the SCAN functional into our comparison. In conclusion, in this work we assessed the performance of a dielectric-dependent hybrid functional for the description of the structure, as well as the dielectric and vibrational properties of liquid water. The fraction of exact exchange chosen here is equal to the inverse of the high-frequency dielectric constant of the liquid. We found that the DDH functional predicts a realistic water structure at room temperature (311 K), similar to that obtained with the SCAN functional 4 but at a temperature closer to ambient conditions. Although the IR spectrum from DDH simulations appears to show a stretching band blue-shifted with respect to experiment, we expect the agreement with measured spectra to substantially improve if nuclear quantum effects are taken into account. Our results, combined with previous reports that global 27 and range-separated 28 dielectric-dependent density functionals predict

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

accurate photoelectron spectra of water 25 and solutions of salts, 26 indicate that these functionals may serve as good theoretical models of both structural and electronic properties of water and aqueous solutions. Our work indicates that several properties of water, including its structure, dielectric and electronic properties, as well as vibrational spectra, can be obtained without explicit inclusion of dispersion interactions beyond what is contained in hybrid functionals. However a full assessment of the role of dispersion forces will require the inclusion of NQEs in the simulations, and work is in progress along these directions.

Computational Methods We modeled water using periodic 64-water-molecule cubic supercells with a side length of 12.414 ˚ A. Deuterium was used instead of hydrogen to allow for a longer simulation time step of 0.242 fs. First-principles simulations were carried out using the Qbox code, 54 employing the Bussi–Donadio–Parrinello thermostat 55 with the time constant of 242 fs. Calculations with the DDH functional used the recursive subspace bisection algorithm 56,57 to accelerate calculation of the exact-exchange integrals. Consistent with our previous work, 9,39 we adopted a bisection threshold of 0.02 and a two-level recursion. Pressure calculations with the DDH functional were carried out with the full exact-exchange functional (no bisection). Optimized norm-conserving Vanderbilt (ONCV) pseudopotentials 58–60 were used to represent the core regions of oxygen and hydrogen. The 1s electron of H and 2s2 2p4 electrons of O were treated explicitly. The simulations were conducted at an energy cutoff of 60 Ry. This computational protocol is at variance with our previous simulations 9,26,39,49 carried out with Hamann–Schl¨ uter–Chiang–Vanderbilt pseudopotentials 61,62 at 85 Ry, but we verified that these protocols provide nearly identical structural properties of the liquid, as shown in the Supporting Information. We used dielectric-dependent density functional of Skone et al. 27 with the fraction of the exact exchange (0.565) equal to the inverse of the experimental high-frequency dielec-

12

ACS Paragon Plus Environment

Page 12 of 21

Page 13 of 21 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

tric constant of water (1.77) derived from the refraction index 63 of water for visible light (1.33 64 ). We chose to use the experimental high-frequency dielectric constant instead of the computed one 27 self-consistently at each step because the experimental and DDH results for representative configurations were found to be close to each other. 65 The choice of the experimental high-frequency dielectric constant is also consistent with the protocol used to compute electronic properties of water. 25 We performed simulations using 8 uncorrelated trajectories running in parallel with both PBE and DDH functionals. Previously, we have successfully used the same computational setup to collect large molecular simulations datasets. 9,13,39,49 Performing several runs permits to determine statistical error bars on computed quantities, and was reported to be a preferable way of estimating dynamical quantities such as diffusion coefficients, compared to a single run of the same length. 46 Overall, we prepared 8 trajectories, 10 ps each, using both PBE and DDH functionals, accumulating 80 ps of total simulation data for each functional. When analyzing the data, we estimated the error as a standard deviation of the mean σmean across the 8 samples. To account for the finite sample size, we multiplied errors by a factor obtained from the Student’s t-test, yielding a two-sided interval at the 98% level of confidence. For 8 independent measurements, the number of degrees of freedom is 7, and σmean was multiplied by 2.998 to yield a ∼ 6σmean interval for the true mean. We adopted this value in all of the analyses with error bars presented in this work.

Acknowledgement The authors gratefully acknowledge helpful discussions with Federico Giberti and Viktor Rozsa, and thank Francesco Paesani for the estimate of the red shift for the O–D stretch mode in the infrared spectrum of heavy water due to nuclear quantum effects. The work was supported by MICCoM as part of the Computational Materials Sciences Program funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES),

13

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

Materials Sciences and Engineering Division (5J-30161-0010A). Work by A.P.G. was also supported by the postdoctoral fellowship from the Natural Sciences and Engineering Research Council of Canada. An award of computer time was provided by the INCITE program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under contract DEAC02-06CH11357.

Supporting Information Available Comparisons of the structure of water simulated with optimized norm-conserving pseudopotentials at 60 Ry and 85 Ry, and with Hamann–Schl¨ uter–Chiang–Vanderbilt pseudopotentials at 85 Ry.

References (1) Gillan, M. J.; Alf`e, D.; Michaelides, A. Perspective: How good is DFT for water? J. Chem. Phys. 2016, 144, 130901. (2) 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. (3) 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. (4) Chen, M.; Ko, H.-Y.; Remsing, R. C.; Andrade, M. F. C.; Santra, B.; Sun, Z.; Selloni, A.; Car, R.; Klein, M. L.; Perdew, J. P.; Wu, X. Ab initio theory and modeling of water. Proc. Natl. Acad. Sci. U. S. A. 2017, 201712499.

14

ACS Paragon Plus Environment

Page 14 of 21

Page 15 of 21 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

(5) Sun, J.; Ruzsinszky, A.; Perdew, J. Strongly Constrained and Appropriately Normed Semilocal Density Functional. Phys. Rev. Lett. 2015, 115, 036402. (6) 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. (7) 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. (8) 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. (9) Gaiduk, A. P.; Gygi, F.; Galli, G. Density and compressibility of liquid water and ice from first-principles simulations with hybrid functionals. J. Phys. Chem. Lett. 2015, 6, 2902–2908. (10) Ben, M. D.; Hutter, J.; VandeVondele, J. Probing the structural and dynamical properties of liquid water with models including non-local electron correlation. J. Chem. Phys. 2015, 143, 054506. (11) 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. (12) 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.

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

(13) Dawson, W.; Gygi, F. Equilibration and analysis of first-principles molecular dynamics simulations of water. J. Chem. Phys. 2018, 148, 124501. (14) 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. (15) 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. (16) Ceriotti, M.; Fang, W.; Kusalik, P. G.; McKenzie, R. H.; Michaelides, A.; Morales, M. A.; Markland, T. E. Nuclear Quantum Effects in Water and Aqueous Systems: Experiment, Theory, and Current Challenges. Chem. Rev. 2016, 116, 7529– 7550. (17) Markland, T. E.; Ceriotti, M. Nuclear quantum effects enter the mainstream. Nat. Rev. Chem. 2018, 2, 0109. (18) Marsalek, O.; Markland, T. E. Ab initio molecular dynamics with nuclear quantum effects at classical cost: Ring polymer contraction for density functional theory. J. Chem. Phys. 2016, 144, 054112. (19) Schran, C.; Marsalek, O.; Markland, T. E. Unravelling the influence of quantum proton delocalization on electronic charge transfer through the hydrogen bond. Chem. Phys. Lett. 2017, 678, 289–295. (20) Marsalek, O.; Markland, T. E. Quantum Dynamics and Spectroscopy of Ab Initio Liquid Water: TheInterplay of Nuclear and Electronic Quantum Effects. J. Phys. Chem. Lett. 2017, 8, 1545–1551. (21) Habershon, S.; Fanourgakis, G. S.; Manolopoulos, D. E. Comparison of path integral

16

ACS Paragon Plus Environment

Page 16 of 21

Page 17 of 21 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

molecular dynamics methods for the infrared absorption spectrum of liquid water. J. Chem. Phys. 2008, 129, 074501. (22) Paesani, F.; Xantheas, S. S.; Voth, G. A. Infrared Spectroscopy and Hydrogen-Bond Dynamics of Liquid Water from Centroid Molecular Dynamics with an Ab Initio-Based Force Field. J. Phys. Chem. B 2009, 113, 13118–13130. (23) Paesani, F. Getting the Right Answers for the Right Reasons: Toward Predictive Molecular Simulations of Water with Many-Body Potential Energy Functions. Acc. Chem. Res. 2016, 49, 1844–1851. (24) Ojha, D.; Henao, A.; K¨ uhne, T. D. Nuclear quantum effects on the vibrational dynamics of liquid water. J. Chem. Phys. 2018, 148, 102328. (25) Gaiduk, A. P.; Pham, T. A.; Govoni, M.; Paesani, F.; Galli, G. Electron affinity of liquid water. Nat. Commun. 2018, 9, 247. (26) Gaiduk, A. P.; Govoni, M.; Seidel, R.; Skone, J. H.; Winter, B.; Galli, G. Photoelectron spectra of aqueous solutions from first principles. J. Am. Chem. Soc. 2016, 138, 6912– 6915. (27) Skone, J. H.; Govoni, M.; Galli, G. Self-consistent hybrid functional for condensed systems. Phys. Rev. B 2014, 89, 195112. (28) Skone, J. H.; Govoni, M.; Galli, G. Nonempirical range-separated hybrid functionals for solids and molecules. Phys. Rev. B 2016, 93, 235106. (29) Laasonen, K.; Sprik, M.; Parrinello, M.; Car, R. “Ab initio” liquid water. J. Chem. Phys. 1993, 99, 9080–9089. (30) Mattsson, A. E.; Mattsson, T. R. AM05 Density Functional Applied to the Water Molecule, Dimer, and Bulk Liquid. J. Chem. Theory Comput. 2009, 5, 887–894.

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

(31) Zhang, C.; Wu, J.; Galli, G.; Gygi, F. Structural and vibrational properties of liquid water from van der Waals density functionals. J. Chem. Theory Comput. 2011, 7, 3054–3061. (32) Babin, V.; Leforestier, C.; Paesani, F. Development of a “first principles” water potential with flexible monomers: Dimer potential energy surface, VRT spectrum, and second virial coefficient. J. Chem. Theory Comput. 2013, 9, 5395–5403. (33) Babin, V.; Medders, G. R.; Paesani, F. Development of a “first principles” water potential with flexible monomers. II: Trimer potential energy surface, third virial coefficient, and small clusters. J. Chem. Theory Comput. 2014, 10, 1599–1607. (34) Medders, G. R.; Babin, V.; Paesani, F. Development of a “first-principles” water potential with flexible monomers. III. Liquid phase properties. J. Chem. Theory Comput. 2014, 10, 2906–2910. (35) 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. (36) Soper, A. K. The Radial Distribution Functions of Water as Derived from Radiation Total Scattering Experiments: Is There Anything We Can Say for Sure? ISRN Phys. Chem. 2013, 2013, e279463. (37) Zeidler, A.; Salmon, P. S.; Fischer, H. E.; Neuefeind, J. C.; Simonson, J. M.; Markland, T. E. Isotope Effects in Water as Investigated by Neutron Diffraction and Path Integral Molecular Dynamics. J. Phys.: Cond. Matter 2012, 24, 284126. (38) Luzar, A.; Chandler, D. Effect of Environment on Hydrogen Bond Dynamics in Liquid Water. Phys. Rev. Lett. 1996, 76, 928–931.

18

ACS Paragon Plus Environment

Page 18 of 21

Page 19 of 21 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

(39) 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. (40) Rastogi, A.; Ghosh, A. K.; Suresh, S. In Thermodynamics - Physical Chemistry of Aqueous Systems; Moreno Piraj´an, J. C., Ed.; InTech, 2011; pp 351–365. (41) Wernet, P.; Nordlund, D.; Bergmann, U.; Cavalleri, M.; Odelius, M.; Ogasawara, H.; N¨aslund, L. A.; Hirsch, T. K.; Ojam¨ae, L.; Glatzel, P.; Pettersson, L. G. M.; Nilsson, A. The Structure of the First Coordination Shell in Liquid Water. Science 2004, 304, 995– 999. (42) Errington, J. R.; Debenedetti, P. G. Relationship between structural order and the anomalies of liquid water. Nature 2001, 409, 318–321. (43) Guardia, E.; Skarmoutsos, I.; Masia, M. Hydrogen Bonding and Related Properties in Liquid Water: A CarParrinello Molecular Dynamics Simulation Study. J. Phys. Chem. B 2015, 119, 8926–8938. (44) Soper, A. K.; Benmore, C. J. Quantum differences between heavy and light water. Phys. Rev. Lett. 2008, 101, 065502. (45) Mills, R. Self-diffusion in normal and heavy water in the range 1–45◦ . J. Phys. Chem. 1973, 77, 685–688. (46) Pranami, G.; Lamm, M. H. Estimating error in diffusion coefficients derived from molecular dynamics simulations. J. Chem. Theory Comput. 2015, 11, 4586–4592. (47) 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.

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

(48) Badyal, Y. S.; Saboungi, M.-L.; Price, D. L.; Shastri, S. D.; Haeffner, D. R.; Soper, A. K. Electron distribution in water. J. Chem. Phys. 2000, 112, 9206–9208. (49) Gaiduk, A. P.; Galli, G. Local and Global Effects of Dissolved Sodium Chloride on the Structure of Water. J. Phys. Chem. Lett. 2017, 8, 1496–1502. (50) Ivanov, S. D.; Witt, A.; Shiga, M.; Marx, D. Communications: On artificial frequency shifts in infrared spectra obtained from centroid molecular dynamics: Quantum liquid water. J. Chem. Phys. 2010, 132, 031101. (51) Max, J.-J.; Chapados, C. Isotope effects in liquid water by infrared spectroscopy. III. H2 O and D2 O spectra from 6000 to 0 cm−1 . J. Chem. Phys. 2009, 131, 184505. (52) Rossi, M.; Liu, H.; Paesani, F.; Bowman, J.; Ceriotti, M. Communication: On the consistency of approximate quantum dynamics simulation methods for vibrational spectra in the condensed phase. J. Chem. Phys. 2014, 141, 181101. (53) Witt, A.; Ivanov, S. D.; Shiga, M.; Forbert, H.; Marx, D. On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy. J. Chem. Phys. 2009, 130, 194510. (54) Gygi, F. Architecture of Qbox: A scalable first-principles molecular dynamics code. IBM J. Res. Dev. 2008, 52, 137–144. (55) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. (56) Gygi, F. Compact representations of Kohn–Sham invariant subspaces. Phys. Rev. Lett. 2009, 102, 166406. (57) Gygi, F.; Duchemin, I. Efficient computation of Hartree–Fock exchange using recursive subspace bisection. J. Chem. Theory Comput. 2013, 9, 582–587.

20

ACS Paragon Plus Environment

Page 20 of 21

Page 21 of 21 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

(58) Hamann, D. R. Optimized norm-conserving Vanderbilt pseudopotentials. Phys. Rev. B 2013, 88, 085117. (59) Schlipf, M.; Gygi, F. Optimization algorithm for the generation of ONCV pseudopotentials. Comput. Phys. Commun. 2015, 196, 36–44. (60) SG15 ONCV pseudopotentials: http://www.quantum-simulation.org/potentials/ sg15_oncv/ (retrieved May 15, 2018). (61) Hamann, D. R.; Schl¨ uter, M.; Chiang, C. Norm-conserving pseudopotentials. Phys. Rev. Lett. 1979, 43, 1494–1497. (62) Vanderbilt, D. Optimally smooth norm-conserving pseudopotentials. Phys. Rev. B 1985, 32, 8412–8415. (63) McQuarrie, D. A. Statistical mechanics; Harper & Row: New York, 1976. (64) Thorm¨ahlen, I.; Straub, J.; Grigull, U. Refractive Index of Water and Its Dependence on Wavelength, Temperature, and Density. J. Phys. Chem. Ref. Data 1985, 14, 933–945. (65) Pan, D.; Wan, Q.; Galli, G. The refractive index and electronic gap of water and ice increase with increasing pressure. Nat. Commun. 2014, 5, 3919.

21

ACS Paragon Plus Environment