Isotope Effects in Liquid Water via Deep Potential Molecular Dynamics

TL;DR

This study combines advanced quantum mechanical simulations, path-integral methods, and machine learning to accurately predict isotope effects on the structure of liquid water, addressing a key challenge in ab initio modeling.

Contribution

It introduces a novel framework integrating deep potential molecular dynamics with path-integral techniques for detailed isotope effect predictions in water.

Findings

Semi-quantitative agreement with experimental isotope effects

Effective sampling of thermal and nuclear quantum fluctuations

Accurate structural property predictions for H2O and D2O

Abstract

A comprehensive microscopic understanding of ambient liquid water is a major challenge for $ab$ $initio$ simulations as it simultaneously requires an accurate quantum mechanical description of the underlying potential energy surface (PES) as well as extensive sampling of configuration space. Due to the presence of light atoms (e.g., H or D), nuclear quantum fluctuations lead to observable changes in the structural properties of liquid water (e.g., isotope effects), and therefore provide yet another challenge for $ab$ $initio$ approaches. In this work, we demonstrate that the combination of dispersion-inclusive hybrid density functional theory (DFT), the Feynman discretized path-integral (PI) approach, and machine learning (ML) constitutes a versatile $ab$ $initio$ based framework that enables extensive sampling of both thermal and nuclear quantum fluctuations on a quite accurate underlying PES. In particular, we employ the recently developed deep potential molecular dynamics (DPMD) model---a neural-network representation of the $ab$ $initio$ PES---in conjunction with a PI approach based on the generalized Langevin equation (PIGLET) to investigate how isotope effects influence the structural properties of ambient liquid H$_2$O and D$_2$O. Through a detailed analysis of the interference differential cross sections as well as several radial and angular distribution functions, we demonstrate that this approach can furnish a semi-quantitative prediction of these subtle isotope effects.

Figures

Tables

Table 1. Table 1: Comparison of DPMD, PI-DPMD, and experimental/EPSR assigned  [ 3 ] average structural properties (in Å) in liquid \ce H2O and \ce D2O at ambient conditions ( 300 300 300  K,  1 1 1  bar). d \ce ​ O ​ H subscript 𝑑 \ce 𝑂 𝐻 d_{\ce{OH}} ( d \ce ​ O ​ D subscript 𝑑 \ce 𝑂 𝐷 d_{\ce{OD}} ): \ce OH ( \ce OD) covalent bond lengths; d \ce ​ O ​ ⋯ ​ \ce ​ H subscript 𝑑 \ce 𝑂 ⋯ \ce 𝐻 d_{\ce{O}\cdots\ce{H}} ( d \ce ​ O ​ ⋯ ​ \ce ​ D subscript 𝑑 \ce 𝑂 ⋯ \ce 𝐷 d_{\ce{O}\cdots\ce{D}} ): \ce O ⋯ ⋯ \cdots \ce H ( \ce O ⋯ ⋯ \cdots \ce D) HB lengths; d \ce ​ H ​ H subscript 𝑑 \ce 𝐻 𝐻 d_{\ce{HH}} ( d \ce ​ D ​ D subscript 𝑑 \ce 𝐷 𝐷 d_{\ce{DD}} ): (shortest) intermolecular \ce HH ( \ce DD) distances. All standard errors are smaller than 0.01 0.01 0.01  Å and are therefore omitted for clarity.

	$d_{\text{ce}OH}$	$d_{\text{ce}OD}$	$d_{\text{ce}O\cdots \text{ce}H}$	$d_{\text{ce}O\cdots \text{ce}D}$	$d_{\text{ce}HH}$	$d_{\text{ce}DD}$
DPMD	$0.98$	$0.98$	1.76	1.76	2.30	2.30
PI-DPMD	$1.00$	$0.99$	$1.76$	$1.76$	$2.32$	$2.30$
Expt./EPSR (Soper 2008)	$1.01$	$0.98$	$1.74$	$1.81$	$2.42$	$2.37$

Equations

$$G^2\to G^2+A\left[1+\text{erf}\left(\frac{\frac{G^2}{2}-E_0}{\sigma }\right)\right].$$

$$\mathcal{L}(p_{ϵ},p_f,p_{\xi })=p_{ϵ}\Delta {ϵ}^2+\frac{p_f}{3N}\sum _i|\Delta {\mathbf{F}}_i{|}^2+\frac{p_{\xi }}{9}\sum _{ij}|\Delta {\xi }_{ij}{|}^2.$$

$$F_{\mathrm{int}}^{(\mathrm{n})}(Q)=\{\begin{array}{l}{c}_{\text{ce}O}^2{b}_{\text{ce}O}^2{S}_{\text{ce}O\text{ce}O}(Q)+2c_{\text{ce}O}{c}_{\text{ce}H}{b}_{\text{ce}O}{b}_{\text{ce}H}{S}_{\text{ce}O\text{ce}H}(Q)+c_{\text{ce}H}^2{b}_{\text{ce}H}^2{S}_{\text{ce}H\text{ce}H}(Q)\text{for ceH2O}\\ c_{\text{ce}O}^2{b}_{\text{ce}O}^2{S}_{\text{ce}O\text{ce}O}(Q)+2c_{\text{ce}O}{c}_{\text{ce}D}{b}_{\text{ce}O}{b}_{\text{ce}D}{S}_{\text{ce}O\text{ce}D}(Q)+c_{\text{ce}D}^2{b}_{\text{ce}D}^2{S}_{\text{ce}D\text{ce}D}(Q)\text{for ceD2O}\end{array}$$

$$S_{\alpha \beta }(Q)=4\pi \rho \int \text{d}r{r}^2\left[g_{\alpha \beta }(r)-1\right]\frac{\sin Qr}{Qr}.$$

$$\begin{array}{rl}\Delta G_{\text{ce}H}(r)=& {\text{FT}}^{-1}[F_{\mathrm{int}}^{(\mathrm{n},\text{ce}{H2}^{18}\text{ce}\mathrm{O})}(Q)-F_{\mathrm{int}}^{(\mathrm{n},\text{ce}H2O)}(Q)]\\ =& c_{\text{ce}O}^2\left(b_{{}^{18}\text{ce}O}^2-b_{{}^{\mathrm{nat}}\text{ce}O}^2\right)\left[g_{\text{ce}OO}^{\text{ce}(H2O)}(r)-1\right]\\ & +2c_{\text{ce}O}{c}_{\text{ce}H}{b}_{\text{ce}H}\left(b_{{}^{18}\text{ce}O}-b_{{}^{\mathrm{nat}}\text{ce}O}\right)\left[g_{\text{ce}OH}^{\text{ce}(H2O)}(r)-1\right],\end{array}$$

$$\begin{array}{rl}\Delta G_{\text{ce}D}(r)=& c_{\text{ce}O}^2\left(b_{{}^{18}\text{ce}O}^2-b_{{}^{\mathrm{nat}}\text{ce}O}^2\right)\left[g_{\text{ce}OO}^{\text{ce}(D2O)}(r)-1\right]\\ & +2c_{\text{ce}O}{c}_{\text{ce}D}{b}_{\text{ce}D}\left(b_{{}^{18}\text{ce}O}-b_{{}^{\mathrm{nat}}\text{ce}O}\right)\left[g_{\text{ce}OD}^{\text{ce}(D2O)}(r)-1\right],\end{array}$$

Full text

††articletype: INVITED ARTICLE

1 Introduction

Despite the relatively simple structure of a single water molecule, water has an unmatched complexity in the condensed (liquid) phase. This complexity mainly originates from the delicate balance between weak non-covalent intermolecular interactions, e.g., the hydrogen bond (HB) network and van der Waals (vdW) dispersion, and thermal (entropic) effects [1]. On top of this delicate balance, nuclear quantum effects (NQEs)—such as zero-point motion—also affect the structural and dynamical properties of liquid water [2]. This is primarily due to the presence of light atoms (such as \ceH) and the central role played by the tetrahedral HB network in determining aqueous properties in the condensed phase.

A salient example of how NQEs affect the structure of liquid water is the experimentally observed isotope effect upon substitution of \ceH to \ceD. First demonstrated by Soper and Benmore using joint X-ray/neutron scattering experiments [3], isotope effects were shown to manifest as significant covalent-bond contraction (by $\approx 3$%) and HB elongation (by $\approx 4$%) when comparing \ceH2O to \ceD2O. In this study, the positions of the individual atoms were assigned based on the so-called empirical potential structure refinement (EPSR) method [4, 5], in which a physically motivated (but empirical) interatomic potential was utilized in conjunction with reverse Monte Carlo (MC) simulations to generate microscopic snapshots of liquid water (\ceH2O or \ceD2O) that reproduce the experimentally observed interference differential cross section, $F^{(\rm n)}_{\rm int}(Q)$ (vide infra). These EPSR-based atomic positions were then used to generate structural properties such as radial distribution functions (RDFs) and the \ceOOO angular distribution function (ADF), from which one can quantify how the \ceH $\rightarrow$ \ceD isotopic substitution affects the structure of liquid water.

The microscopic structure of liquid water can also be modeled using ab initio molecular dynamics (AIMD) simulations, in which the nuclear potential energy surface (PES) is generated “on the fly” from the electronic ground state without the need for empirical input [6, 7]. By generating a trajectory of configurations based on interatomic forces derived from a first-principles electronic structure theory, this technique allows for a quantum mechanical treatment of the structural and dynamical properties of complex condensed-phase systems (such as liquid water) at a given set of thermodynamic conditions. In addition, AIMD simulations can also describe the electronic and dielectric properties of the system, as well as any chemical reactions that may occur (i.e., bond cleavage/formation).

Due to its quite favorable balance between accuracy and computational cost, density functional theory (DFT) [8, 9, 10] has emerged as the most commonly used electronic structure theory method during AIMD simulations of condensed-phase systems. While DFT is (in principle) an exact theory, the functional form of the exchange-correlation energy still remains unknown; as such, DFT (in practice) relies on an established hierarchy of approximations which allows for (semi-)systematic improvements in accuracy with a corresponding increase in the computational cost [11]. Previous studies have demonstrated that generalized-gradient approximation (GGA) based DFT functionals [12, 13, 14]—due to their propensity for self-interaction error [15] and lack of non-local correlation effects such as vdW dispersion interactions [16]—are inadequate for providing an accurate and reliable description of liquid water. Instead, it is more appropriate to use constraint-based meta-GGA functionals (such as SCAN [17, 18, 19, 20]) or the class of more accurate vdW-inclusive hybrid functionals (such as PBE0-TS [21, 22, 23, 1, 24, 25, 26] or revPBE0-D3 [27, 21, 22, 28, 29]) when describing condensed-phase aqueous systems.

During AIMD simulations, the nuclei are often treated as classical (point) particles moving on the underlying PES generated from an “on-the-fly” and fully quantum mechanical description of the electrons. To rigorously include NQEs while sampling equilibrium structural properties, the Feynman discretized path-integral (PI) [30, 31, 32, 33] approach is typically used. In this approach, each nucleus is mapped onto a classical ring polymer comprised of $P$ beads (coupled via harmonic springs), allowing for a so-called PI-AIMD simulation of the system in which the electronic and nuclear degrees of freedom are both treated quantum mechanically. At finite (but non-zero) temperatures, sampling of the exact quantum distribution of a system (in configuration space) can be accomplished to a desired level of accuracy with PI-AIMD simulations employing a finite number of beads. However, standard PI-AIMD simulations of liquid water at ambient conditions ($300$ K, $1$ bar) are still computationally demanding for the following reasons: (i) they require many beads ($P\geq 32$) in each classical ring polymer to provide a sufficiently converged description of equilibrium structural properties [34, 35, 36, 37]; (ii) they also require relatively long trajectories ($t>20$ ps) to obtain sufficient statistical convergence for many physical quantities, e.g., RDFs and the equilibrium density. The difficulty originating from the required number of beads has been the subject of intense research and can now be alleviated (to varying extents) using a number of different methods, e.g., generalized Langevin equation (GLE) based colored-noise thermostats [38, 39, 40], ring polymer contraction [41, 42, 43, 44, 45], and imaginary-time perturbed path integral approach [46, 47], to name a few. The sampling time issue is an intrinsic limitation of PI-AIMD (and AIMD) simulations that can only be addressed by more efficient algorithms and/or approximations for evaluating the underlying ab initio PES.

Recent progress in the application of machine-learning (ML) techniques to represent PES for complex systems containing many atoms provides a way to overcome this last computational hurdle [48, 49, 50, 51, 52, 53, 54, 55, 56, 57]. After training on high-quality ab initio PES data, ML-based methods—with a linear-scaling associated computational cost that is many orders of magnitude cheaper than AIMD—allow one to perform extensive MD simulations without any significant loss of accuracy. In fact, such methods have already been used to study the complex behavior of liquid water [58, 59, 56, 29]. For example, Cheng and coworkers [29] recently used a neural network (NN) based potential to study the thermodynamics of liquid water and ice I (I$h$ and I$c$) under ambient conditions; in doing so, they accurately predicted the RDFs of liquid water, the melting point of ice I$h$, and the density maximum in liquid water. Moreover, suitable extensions of these methods also allow one to further coarse-grain the system into beads representing each water molecule; this has the potential to significantly reduce the computational cost without losing the main structural information in the original system [60]. Among the suite of existing ML methods, the deep potential molecular dynamics (DPMD) approach [55, 56, 57] utilizes a deep NN to represent the many-body potential energy as a sum of auxiliary “atomic” energies associated with individual atoms in the system. For a given atom, this “atomic” energy depends differentiably on its local environment, i.e., the relative coordinates of its vicinal surrounding atoms within a smooth radial distance cutoff. A novel symmetry-preserving atomic local coordinate system is then introduced as input to the deep NN and an adaptive training procedure on high-quality ab initio PES data is performed to optimize the NN parameters. By reproducing the total potential energy, atomic forces, and stress tensor to well within the accuracy of the ab initio data [55, 56, 57], the resulting DPMD model can recover structural properties that are essentially indistinguishable when compared to the AIMD results for a wide variety of systems (including liquid water).

In this work, we utilize DPMD and PI-DPMD simulations to investigate the experimentally observed structural changes [3] in liquid water upon isotopic substitution of \ceH2O with \ceD2O. In particular, we demonstrate that the DPMD model—when trained on a single, relatively short ($\approx 8$ ps) PI-AIMD simulation of the lighter isotope (liquid \ceH2O)—can be used to provide a (semi-)quantitative description of the isotope effects found in $F^{(\rm n)}_{\rm int}(Q)$, the oxygen–oxygen (\ceOO) RDF, and the oxygen–oxygen–oxygen (\ceOOO) ADF. We rationalize this choice of training data by the fact that \ceH2O and \ceD2O share the same underlying PES, but are subject to different quantum fluctuations; hence, a PI-AIMD simulation of the single lighter isotope (liquid \ceH2O) covers the relevant sector of configuration space for both isotopes. This conclusion is supported by independently testing that the DPMD model trained on quantum \ceH2O describes with similar accuracy classical \ceH2O. In doing so, we show that the combination of vdW-inclusive hybrid DFT as the source of the underlying PES, the Feynman discretized PI approach for treating NQEs, and the highly efficient DPMD model constitutes a versatile ab initio based framework that enables one to attack problems in chemistry, physics, and materials science that require extensive and simultaneous sampling of both thermal and nuclear quantum fluctuations on a quite accurate PES.

The remainder of the manuscript is organized as follows. In Sec. 2, we describe the computational details required for obtaining the initial PES using vdW-inclusive hybrid DFT and the subsequent generation of the DPMD model parameters needed for performing DPMD and PI-DPMD simulations of liquid \ceH2O and \ceD2O. In Sec. 3, we perform a detailed case study of how isotope effects manifest in the equilibrium structural properties of liquid water, including a comparative analysis of our theoretical findings with the experimental work of the Soper and Benmore [3]. In Sec. 4, we provide a brief summary of this work accompanied by some remarks about potential improvements to this approach for the study of liquid water and its many unusual properties.

2 Computational methods

Our study of the isotope effects in liquid \ceH2O and \ceD2O via DPMD and PI-DPMD simulations was performed in three stages. The first stage is the generation of training data for the DPMD model based on a PI-AIMD simulation of liquid \ceH2O with vdW-inclusive hybrid DFT (Sec. 2.1). The second stage is the training of the DPMD model parameters based on this PI-AIMD trajectory (Sec. 2.2). The third and final stage is the utilization of this model in DPMD and PI-DPMD simulations to equilibrate and sample the statistical ensembles of liquid \ceH2O and \ceD2O (Sec. 2.3).

2.1 Generation of DPMD training data via PI-AIMD simulation of liquid \ceH2O

To generate the ab initio PES data for liquid water that is required for training the DPMD model, we performed a single isobaric-isothermal ($NpT$) Born-Oppenheimer PI-AIMD simulation (with a time step of $\Delta t\approx 0.5$ fs) of a cubic periodic cell containing $64$ \ceH2O molecules at ambient ($300$ K and $1.0$ bar) conditions. We started from a previously equilibrated $NpT$ PI-AIMD trajectory at the PBE0-TS level carried out at $330$ K and $1.0$ bar [61], and performed an initial $\approx 1$ ps $NpT$ PI-AIMD simulation in which the temperature was gradually lowered from $330$ K to $300$ K. A final $NpT$ PI-AIMD production run was then performed for $\approx 8$ ps under ambient conditions of $300$ K and $1$ bar; this ab initio PES data will be used to train the DPMD model as described in Sec. 2.2.

In the PI-AIMD simulation, the nuclei were modeled with an $8$-bead ring polymer supplemented with a colored-noise generalized Langevin equation thermostat (i.e., PIGLET) [39, 40] to accelerate convergence with respect to the Trotter dimension. The cell was treated classically using the Raiteri-Gale-Bussi approach [62] and was thermostatted with an additional generalized Langevin equation thermostat [63]. The cell mass was chosen to be consistent with a $200$ fs characteristic timescale, and the corresponding equations of motion for the ionic and cell degrees of freedom were integrated using the i-PI package [64]. To prevent shearing of the liquid, cubic cell symmetry was enforced during all $NpT$ simulations.

The ionic forces and stress tensor (which are required by i-PI) were computed at the PBE0-TS level [21, 22, 23] using a linear-scaling exact-exchange algorithm [65, 1, 66, 67] and a self-consistent implementation of the Tkatchenko-Scheffler (TS) vdW correction [23, 24], as provided in the CP package of Quantum ESPRESSO [68, 69]. The core electrons were treated with Hamann-Schlüter-Chiang-Vanderbilt (HSCV) norm-conserving pseudopotentials [70, 71], as distributed with the Qbox suite of programs [72], while the valence (pseudo-)wavefunctions were represented explicitly with a planewave basis set. To maintain a constant planewave kinetic energy cutoff of $115$ Ry during the $NpT$ simulation, we followed the procedure of Bernasconi et al. [73] by choosing: (i) a cubic reference cell (with $L=24.7$ Bohr) that is large enough to cover the fluctuations along each lattice vector of the simulation cell throughout the $NpT$ trajectory, and (ii) a corresponding planewave basis set with a larger kinetic energy cutoff of $125$ Ry. During the $NpT$ simulation, planewaves with a kinetic energy beyond the desired cutoff of $115$ Ry were smoothly penalized by changing

[TABLE]

With a judicious choice of parameters ($A=200$ Ry, $\sigma=15$ Ry, $E_{0}=115$ Ry), this modification to $G^{2}$ causes the higher-energy ($>115$ Ry) planewaves to become essentially inactive basis functions in the description of the valence (pseudo-)wavefunctions, and thereby leads to $NpT$ dynamics which mimic a constant planewave cutoff of $\approx 115$ Ry.

The electronic ground state was obtained using nested second-order damped Car-Parrinello dynamics [74] in which the ionic and cell degrees of freedom were kept fixed. During the electron minimization, we employed a time step of $\approx 0.125$ fs, a fictitious electron mass of $200$ au, and a mass preconditioning cutoff of $6$ Ry. Additional computational details regarding this PI-AIMD simulation will be provided in a forthcoming paper [61]. The electronic structure calculation is considered converged when the successive second-order-damped CP dynamics has a change in total potential energy below $1.1\times 10^{-6}$ au and a change of mean absolute ionic forces within $3.0\times 10^{-4}$ au.

2.2 Training of the DPMD model

To train the DPMD model, we follow Ref [56] to determine a deep NN potential based on the total potential energy $E$, the ionic forces $\{\bm{F}_{i}\}$ (for all $N$ ions in the system), and the virial $\bm{\Xi}$ (or equivalently, the stress tensor) at each step in the $\approx 8$ ps production PI-AIMD simulation described above. During the training of the DPMD model, we minimize the following loss function $\mathcal{L}$, which depends on the tunable parameters $p_{\epsilon}$, $p_{f}$, and $p_{\xi}$, and is defined as follows:

[TABLE]

In this expression, $\Delta\epsilon$, $\Delta\bm{F}_{i}$, and $\Delta\bm{\xi}$ are the differences between the current DPMD prediction and the training data for the quantities $\epsilon\equiv E/N$, $\bm{F}_{i}$, and $\bm{\xi}\equiv\bm{\Xi}/N$, respectively. To optimize the parameters in the DPMD model, we used the Adam method [75] with a learning rate that decays exponentially with the number of training steps. In practice, we observed that more efficient training can be achieved by linearly varying the parameters with respect to the learning rate (i.e., increasing $p_{\epsilon}$ and $p_{\xi}$ while decreasing $p_{f}$); in doing so, one can achieve a well-balanced training procedure in which the energy, ionic forces, and virial are mutually consistent. All DPMD model training was performed using the DeePMD-kit package [76] interfaced with the TensorFlow library [77].

2.3 DPMD and PI-DPMD simulations of liquid \ceH2O and \ceD2O

To reduce finite-size effects and the statistical errors associated with relatively short PI-AIMD simulations, we extended the system size by one order of magnitude to include $512$ \ceH2O molecules and the simulation timescale by two orders of magnitude to span $\approx 1$ ns during the DPMD and PI-DPMD simulations described below. A single $NpT$ DPMD simulation (in which the nuclei are treated classically) of liquid water at ambient ($300$ K and $1$ bar) conditions was performed using the LAMMPS package [78]. The temperature and pressure were controlled by Nosé-Hoover thermostat and barostat chains [79] (with characteristic timescales of $0.1$ ps and $0.5$ ps, respectively), and an integration time step of $\approx 0.5$ fs. To perform PI-DPMD simulations of liquid \ceH2O and \ceD2O, we integrated the DeePMD-kit and i-PI packages. All PI-DPMD simulations were performed using exactly the same settings as the original PI-AIMD simulation (see Sec. 2.1), including identical thermodynamic conditions, sampling techniques, and time step. We then made the substitution of \ceH to \ceD, and repeated the PI-DPMD simulation for heavy water. For both the DPMD and PI-DPMD simulations, we enforced cubic cell symmetry as described in Sec. 2.1.

3 Results and discussion

By performing DPMD and PI-DPMD simulations trained on ab initio PES data with extended time ($\approx 1$ ns) and length (\ce(H2O)512) scales, this approach mitigates several of the computational challenges associated with modeling complex condensed-phase systems like liquid water. As an application, we use this ab initio based framework to investigate how isotope effects influence the structural properties of ambient liquid \ceH2O and \ceD2O. Based on extensive $NpT$ PI-DPMD simulations of ambient liquid \ceH2O and \ceD2O, we obtained essentially identical equilibrium atomic number densities of $\rho=0.1022\pm 0.0001$ atoms/Å3 for these systems. This value for $\rho$ is slightly higher than the experimental observations of $\rho=0.10007$ atoms/Å3 (\ceH2O) and $0.10000$ atoms/Å3 (\ceD2O) at $298$ K. While the isotope effects on macroscopic properties such as $\rho$ are negligible, the isotope effects on the microscopic structural properties of liquid water are still small (on the order of $3\mathrm{-}4\%$ from experiment [3]), but well above the statistical uncertainty of our simulations. In this section, we investigate the isotope effects on the microscopic structural properties of liquid water predicted by the PI-DPMD simulations. In particular, we first compare against the interference differential cross section $F^{(\rm n)}_{\rm int}(Q)$—a reciprocal-space quantity that is directly measured by experiment—to initially assess the accuracy of our simulations (Sec. 3.1). This is followed by a detailed comparative analysis with the following experimentally (EPSR) determined real-space structural quantities: the oxygen–oxygen (\ceOO), oxygen–hydrogen (\ceOH/\ceOD), and hydrogen–hydrogen (\ceHH/\ceDD) radial distribution functions (RDFs) in Sec. 3.2, and the oxygen–oxygen–oxygen (\ceOOO) angular distribution function (ADF) in Sec. 3.3. We complete this section by describing a potential route towards improving the predictions of our approach using other electronic structure theory models (Sec. 3.4).

3.1 Isotope effects on the microscopic structure of liquid water: Interference differential cross section

To provide a comparison with a directly measurable experimental quantity (and sidestep any potential bias from the EPSR post-processing of the raw experimental data), we first computed $F^{(\rm n)}_{\rm int}(Q)$ from our PI-DPMD trajectories as follows:

[TABLE]

in which $c_{\alpha}$ is the atom fraction for the $\alpha$-th species ($c_{\ce{O}}=1/3$, $c_{\ce{H}}=2/3$, and $c_{\ce{D}}=2/3$), $b_{\alpha}$ is the neutron scattering length of the $\alpha$-th species ($b_{\ce{O}}=5.80$ fm, $b_{\ce{H}}=-3.74$ fm, and $b_{\ce{D}}=6.67$ fm),111Data source: https://www.ncnr.nist.gov/resources/n-lengths/elements/ and $S_{\alpha\beta}(Q)$ is the partial structure factor for the $\alpha\beta$ species pair. In this work, the reciprocal-space $S_{\alpha\beta}(Q)$ quantity was computed based on the corresponding RDF for the $\alpha\beta$ species pair ($g_{\alpha\beta}(r)$) as follows:

[TABLE]

In Fig. 1, we compare the $F^{(\rm n)}_{\rm int}(Q)$ computed from DPMD simulations (in which the nuclei are treated classically) and PI-DPMD simulations (in which the nuclei are treated quantum mechanically) against the experimental results of Soper and Benmore [3]. For $Q<15$ Å*-1*, the line shape of $F^{(\rm n)}_{\rm int}(Q)$ originates from the combined shapes of $S_{\ce{OO}}$, $S_{\ce{OH}}$ ($S_{\ce{OD}}$), and $S_{\ce{HH}}$ ($S_{\ce{DD}}$). For higher $Q$, the \ceOH (\ceOD) contribution dominates $F^{(\rm n)}_{\rm int}(Q)$. As a consequence, the relative amplitude and (approximately) opposite phase in $F^{(\rm n)}_{\rm int}(Q)$ for $Q>15$ Å*-1* between \ceH2O and \ceD2O mainly arise from the ratio in the scattering length $b_{\ce{H}}/b_{\ce{D}}\approx-0.56$. From the apparent differences in $F^{(\rm n)}_{\rm int}(Q)$ between the DPMD and PI-DPMD simulations (for both \ceH2O and \ceD2O), it is clear that NQEs play an important role in determining the structural properties of ambient liquid water. By accounting for NQEs in the PI-DPMD simulations, the theoretical predictions for $F^{(\rm n)}_{\rm int}(Q)$ are in significantly improved agreement with experiment. Even with extensive sampling of thermal and nuclear quantum fluctuations, however, small differences still remain between the theoretical and experimental results. In this regard, more pronounced discrepancies were observed in the $F^{(\rm n)}_{\rm int}(Q)$ for \ceH2O; since this quantity is a linear combination of $S_{\ce{OO}}$, $S_{\ce{OH}}$, and $S_{\ce{HH}}$ via Eqs. (3) and (4), one can argue that such differences primarily arise from deficiencies in the theoretical determinations of these structure factors. Assuming that our sampling of the thermal and nuclear quantum fluctuations in this system is indeed exhaustive and the DPMD model reproduces the underlying ab initio PES [56, 57], this finding suggests that further improvement may be needed in the model used to describe the electronic structure. Although the nuclei in \ceH2O and \ceD2O move on the same PES in the Born-Oppenheimer approximation, the lighter \ceH atoms are more quantum mechanical in nature and will explore more of the PES than the heavier (and more classical) \ceD atoms at a given temperature. As such, PI-DPMD simulations of liquid \ceH2O are more sensitive to the finer details of the underlying PES (e.g., surface curvature, barrier heights, and number/location/topology of local minima). Since small deviations from the true PES (as derived from the exact solution to the time-independent Schrödinger equation) will lead to incorrect sampling of this surface during PI simulations (e.g., overestimation/underestimation of anharmonicity, incorrect relative populations of different minima/basins), further improvements in the electronic structure model are expected to improve our theoretical prediction of $F^{(\rm n)}_{\rm int}(Q)$ for liquid water. Additional support for this conjecture is provided in Sec. 3.2.

3.2 Isotope effects on the microscopic structure of liquid water: Radial distribution functions

To perform a more intuitive analysis of the isotope effects on the microscopic structure of liquid water, we now consider how the isotopic substitution of \ceH to \ceD influences the oxygen–oxygen (\ceOO), oxygen–hydrogen (\ceOH/\ceOD), and hydrogen–hydrogen (\ceHH/\ceDD) radial distribution functions (RDFs) of \ceH2O and \ceD2O. RDFs are real-space quantities that measure the probability of finding a pair of atoms (of a given type) as a function of their radial distance; these quantities are related to the structure factors used in determining $F^{(\rm n)}_{\rm int}(Q)$ above (cf. Eqs. (3) and (4)), and are commonly used in studying the structure of liquids. Based on the DPMD and PI-DPMD simulations described above, we computed the $g_{\ce{OO}}(\bm{r})$, $g_{\ce{OH}}(\bm{r})$, and $g_{\ce{HH}}(\bm{r})$ RDFs for liquid \ceH2O ($g_{\ce{OO}}(\bm{r})$, $g_{\ce{OD}}(\bm{r})$, and $g_{\ce{DD}}(\bm{r})$ for liquid \ceD2O) and plotted them in Fig. 2.

From this figure, one can see that the isotopic substitution \ceH $\rightarrow$ \ceD slightly increases the structure of each RDF. This difference primarily originates from the reduced zero-point motion in \ceD2O, and is consistent with the changes observed in the RDFs of liquid \ceH2O and \ceD2O obtained from simulations which account for NQEs [2]. Before entering into a detailed discussion on the $g_{\ce{OO}}(\bm{r})$—a quantity that can also be determined via X-ray scattering experiments [80]—we first analyze the \ceOH/\ceOD and \ceHH/\ceDD RDFs.

Following the procedure outlined by Soper and Benmore [3], we computed the following structural properties based on our PI-DPMD simulations of liquid \ceH2O and \ceD2O: (i) the \ceOH and \ceOD covalent bond lengths ($d_{\ce{OH}}$ and $d_{\ce{OD}}$), (ii) the \ceO$\cdots$\ceH and \ceO$\cdots$\ceD HB lengths ($d_{\ce{O}\cdots\ce{H}}$ and $d_{\ce{O}\cdots\ce{D}}$), and (iii) the (shortest) intermolecular \ceHH and \ceDD distances ($d_{\ce{HH}}$ and $d_{\ce{DD}}$). This data is summarized in Table 1, where we find that PI-DPMD simulations predict that isotopic substitution from \ceH to \ceD leads to an $\approx 1\%$ decrease in the covalent bond, a negligible ($<1\%$) change in the HB length, and an $\approx 1\%$ decrease in the (shortest) intermolecular \ceHH/\ceDD distance. These findings are consistent with but smaller than the experimental/EPSR assignment of an $\approx 3\%$ contraction in $d_{\ce{OD}}$, an $\approx 4\%$ elongation in $d_{\ce{O}\cdots\ce{D}}$, and an $\approx 2\%$ contraction in $d_{\ce{DD}}$. The experimental/EPSR assignment of these structural quantities requires post-processing of the directly observable experimental data (i.e., $F^{(\rm n)}_{\rm int}(Q)$). A direct experimental probe of the intramolecular $g_{\ce{OH}}(r)$ ($g_{\ce{OD}}(r)$) would be required to better assess the intramolecular isotope effect. This was provided by an elegant oxygen istope substitution experiment by Zeidler et al. [81, 82], in which they measured the difference between the neutron scattering cross section of \ceH218\ceO (i.e., $F_{\rm int}^{\rm(n,\ce{H2}^{18}\ce{O})}(Q)$) and of naturally occurring \ceH2O (i.e., $F_{\rm int}^{\rm(n,\ce{H2O})}(Q)$). If the oxygen isotope effect on $g^{\ce{(H2O)}}_{\ce{OO}}(r)$ and $g^{\ce{(H2O)}}_{\ce{OH}}(r)$ (as well as on $g^{\ce{(D2O)}}_{\ce{OO}}(r)$ and $g^{\ce{(D2O)}}_{\ce{OD}}(r)$) is negligible, as one would expect, the inverse Fourier transform (FT*-1*) of the measured cross sections would give for light and heavy water, respectively:

[TABLE]

[TABLE]

in which $b_{{}^{18}\ce{O}}=6.01$ fm and $b_{{}^{\rm nat}\ce{O}}=5.80$ fm are the relevant neutron scattering lengths. In the $r$ range corresponding intramolecular \ceOH (\ceOD) bonds, $g^{\ce{(H2O)}}_{\ce{OO}}(r)$ and $g^{\ce{(D2O)}}_{\ce{OO}}(r)$ are both equal to zero. Thus, the intramolecular $g^{\ce{(H2O)}}_{\ce{OH}}(r)$ and $g^{\ce{(D2O)}}_{\ce{OD}}(r)$ can be extracted from Eqs. (5) and (6), respectively. The results are displayed in the upper panel of Fig. 3, which also shows in the middle panel the results of our simulation, and in the lower panel the corresponding RDFs from EPSR [3] as reported in Ref. [82].

The data points from the oxygen substitution experiment are few and affected by uncertainty, reflecting the difficulty of the experiment. However, the \ceOH to \ceOD bond contraction is clear and in qualitative agreement with both EPSR and PI-DPMD data. To quantify the effect, Zeidler et al. estimated the peak position of the RDFs by spline fitting obtaining a contraction of $\approx 0.5$%, smaller than the $\approx 3$% effect estimated from EPSR but closer to our PI-DPMD estimate of $\approx 1$%. Interestingly, the PI-DPMD RDFs are asymmetric around the peak position, reflecting the anharmonicity of the potential experienced by \ceH (\ceD), which is steeper when \ceH (\ceD) approaches oxygen. In spite of the limited resolution, such anharmonicity seems to be present in the oxygen substitution experiment. It is absent, however, in the EPSR data in which a harmonic approximation was used for the \ceOH (\ceOD) potential. Further experiments would seem necessary to quantify accurately the \ceOH to \ceOD bond contraction. However, a $3$% contraction seems excessive in view of the fact that a contraction of only $\approx 0.4$% is measured in the vapor phase for the molecule in the roto-vibrational ground state [83], and a similar contraction seems to occur also in ice I$h$, albeit with a larger uncertainty [84].

In spite of the uncertainties that affect both theory and experiment significant physical insight can be gleaned from the comparison of PI-DPMD simulations and EPSR structures. Both approaches characterize the isotope effects as small and subtle changes to the microscopic structure of liquid water. This finding is also consistent with the differences observed between the predicted structural properties from DPMD simulations (in which the nuclei are treated classically) and PI-DPMD simulations of heavy (and hence more classical) liquid water (see Fig. 2).

Now we turn our discussion to $g_{\ce{OO}}(\bm{r})$ in liquid \ceH2O and \ceD2O. In Fig. 4, we plot this quantity as computed from PI-DPMD simulations, the experimental/EPSR procedure [3], and a benchmark X-ray diffraction experiment [80]. Focusing on liquid \ceH2O, PI-DPMD underestimates the \ceOO distance, as its $g_{\ce{OO}}(r)$ is characterized by a first peak that is shorter by $\approx 0.07$ Å (around $2.5$%) than the benchmark X-ray experiment value of $2.80$ Å. Interestingly, the benchmark X-ray assignment of the first-peak position also agrees very well with a more sophisticated quantum Monte Carlo prediction at $2.80$ Å [85], reflecting the limitation of the PBE0-TS model. Such underestimation in the shortest \ceOO distance is consistent with previous simulations [1, 85]. This underestimation in the shortest \ceOO distance is likely a consequence of the residual self-interaction error that is still present in the PBE0 functional and has also been reported for water clusters [86, 87]. When compared with the benchmark X-ray result, the experimental/EPSR assignment also slightly underestimates the position of the first peak in the $g_{\ce{OO}}(\bm{r})$ at $2.77$ Å. While $g_{\ce{OO}}(\bm{r})$ in both assignments are based on X-ray scattering data, the benchmark result seems more accurate due to its higher resolution—a resolution up to $Q_{\rm max}\approx 24$ Å*-1* [80], which is higher than the measurement with $Q_{\rm max}\approx 18$ Å*-1* [88] used in the experimental/EPSR procedure [3]. Besides the $Q_{\rm max}$ probed in a given X-ray experiment, the assignment of $g_{\ce{OO}}(\bm{r})$ should also be sensitive to the form factors used in deriving $S_{\ce{OO}}(\bm{Q})$ from the scattering intensity. In the benchmark experiment [80], the form factors are derived from the electron density of a gas phase \ceH2O molecule [89] instead of one in the liquid phase. However, such condensed-phase effect on the form factors has been demonstrated negligible ($\approx 1$–$2$%) and concentrated in the low-$Q$ region ($Q<2$ Å*-1*) [1], thereby reinforcing the central role played by probing sufficiently high $Q_{\rm max}$ in the accurate experimental determination of the first peak in $g_{\ce{OO}}(\bm{r})$. Apart from the slightly overstructured $g_{\ce{OO}}(\bm{r})$ within the second (coordination) shell (i.e., $r<5.5$ Å), the PI-DPMD prediction agrees quite well with both experimental assignments. Interestingly, the PI-DPMD prediction yields a long-range structural correlation beyond the third shell ($r>8$ Å) as assigned in the benchmark experiment to be observable up to the sixth shell ($r\approx 15$ Å) [80].

When compared with the experimental/EPSR assignment, we observed similar isotope effect in $g_{\ce{OO}}(\bm{r})$ due to \ceH$\rightarrow$\ceD beyond the first peak ($r>3$ Å) as shown in the lower panel of Fig. 4. However, the experimental/EPSR approach assigns an isotope increase in the height of the first peak from $2.59$ (\ceH2O) to $2.84$ (\ceD2O) that is significantly larger than PI-DPMD prediction ($2.58\rightarrow 2.66$). This difference could arise from the limitation of the PBE0-TS model in the PI-DPMD simulation or the challenge to probe the high-$Q$ region in the X-ray scattering experiment. To resolve this issue, it seems necessary to resort to either a high-resolution X-ray experiment on liquid \ceD2O as performed in Ref. [80] or a PI-DPMD simulation trained on realistic electronic structure theory beyond the PBE0-TS level (see Sec. 3.4).

3.3 Isotope effects on the microscopic structure of liquid water: Angular distribution function and tetrahedrality

In this section, we consider how the isotopic substitution of \ceH to \ceD influences the oxygen–oxygen–oxygen (\ceOOO) angular distribution function (ADF) of \ceH2O and \ceD2O. In this analysis, an \ceOOO triplet is defined as any three \ceO atoms ($\ce{O}_{A}$, $\ce{O}_{B}$, and $\ce{O}_{C}$) in which $\ce{O}_{A}$ and $\ce{O}_{B}$ are located within a prescribed cutoff distance $d_{3}$ from $\ce{O}_{C}$. Following the procedure outlined by Soper and Benmore [3], $d_{3}$ was chosen to yield an average \ceO–\ceO coordination number of $4.0$; in this work, we used $d_{3}=3.24$ Å which is consistent with previous AIMD simulations at the PBE0-TS level [1] but larger than the experimental/EPSR value of $3.18$ Å [3]. In Fig. 5, we plot the $P_{\ce{OOO}}(\theta)$ ADF—the probability of finding the triplet angle ($\theta=\angle\ce{O}_{A}\ce{O}_{C}\ce{O}_{B}$)—from the PI-DPMD predictions and experimental/EPSR procedure [3] for liquid \ceH2O and \ceD2O.

In the ADF, the PI-DPMD prediction shows consistent isotope effects in which \ceH2O is less structured than \ceD2O; this manifests as a reduced peak at $\approx 100^{\circ}$ (due to the tetrahedral HB network) and an enhanced peak at $\approx 50^{\circ}$ (due to highly distorted broken HB configurations). Compared with the experiment/EPSR assignment, the isotope effect predicted by the PI-DPMD simulations is once again smaller in extent (particularly for the peak at $\approx 50^{\circ}$). While this discrepancy may due in part to the limitation of the PBE0-TS model, it could also reflect the remaining uncertainty from the experimental/EPSR assignment, which seems to overestimate the isotope effect (e.g., in the \ceOH to \ceOD bond contraction as discussed in Sec. 3.2). On an absolute scale, our predicted ADFs are slightly more structured with a higher main peak at $\approx 100^{\circ}$ and a lower secondary peak at $\approx 50^{\circ}$. This phenomenon is consistent with our observations of the slightly overstructured $g_{\ce{OO}}(\bm{r})$ within the second shell (Sec. 3.2). However, unlike the $g_{\ce{OO}}(\bm{r})$, the ADF is not directly probed from the experiment and relies on the EPSR procedure. As such, more direct experiment or accurate simulation seem necessary to quantify the ADFs and the associated isotope effects in ambient liquid water.

3.4 Potential improvement and future outlook

It is well known that AIMD simulations which employ common GGA functionals (e.g., PBE [14] and BLYP [12, 13]) for the underlying PES yield a liquid water that is substantially overstructured [90, 91, 92, 1, 93, 94, 95, 18, 19]. This effect is reduced but not eliminated with the inclusion of NQEs via PI-AIMD simulations [96]. In the present work, further improvement is achieved by modeling the electronic structure at the PBE0-TS level, which includes a fraction of exact exchange to mitigate the self-interaction error and an effective-pairwise treatment of the non-local vdW interactions. As described throughout this work, this improved electronic structure is not sufficient to provide a quantitative description of the structural properties of liquid water, and there is still room for further improvement. At the level of vdW-inclusive hybrid DFT, the revPBE0-D3 functional [27, 21, 22, 28] appears to predict a $g_{\ce{OO}}(r)$ in better overall agreement with experiment [29]. However, that functional predicts smaller NQEs than PBE0-TS; for example, the differences between the classical and quantum RDFs for \ceH2O reported in Ref. [29] are smaller than that reported in this paper. This observation suggests that revPBE0-D3 should also underestimate the isotope effects in the microscopic structural properties of liquid water. These general difficulties of vdW-inclusive hybrid DFT most likely originate from the imperfect cancellation of the self-interaction error inherent from their approximated functional forms. To overcome these difficulties, it seems necessary to apply a more sophisticated electronic structure theory such as self-interaction-free DFT methods or highly accurate quantum chemical (wavefunction theory) approaches. These approaches are substantially more expensive than hybrid DFT, and can only be used to obtain accurate interatomic interactions (including the total potential energy, atomic forces, and stress tensor) for a small subset of liquid water configurations. As such, definite progress beyond vdW-inclusive hybrid DFT is possible if a ML-based model (like DPMD) could be refined with a select set of accurate quantum chemical data or direct observations from experiments (with approaches such as EPSR). Provided with such an accurate and efficient interatomic potential, we expect that DPMD and PI-DPMD simulations can also be used to study structural and thermal properties of molecular crystals bound by HB interactions and dispersion forces [97, 98, 99], which require the use of hybrid DFT corrected with a beyond effective-pairwise many-body dispersion (MBD) model [100, 101, 102].

4 Conclusions

In this work, we demonstrate that the deep neural network based DPMD model—when trained on a single, relatively short ($\approx 8$ ps) PI-AIMD simulation of liquid \ceH2O—can provide a (semi-)quantitative description of the isotope effects found in $F^{(\rm n)}_{\rm int}(Q)$, the $g_{\ce{OO}}(\bm{r})$, $g_{\ce{OH}}(\bm{r})$/$g_{\ce{OD}}(\bm{r})$, and $g_{\ce{HH}}(\bm{r})$/$g_{\ce{DD}}(\bm{r})$ RDFs, and the $P_{\ce{OOO}}(\theta)$ ADF. Like the experimental/EPSR procedure of Soper and Benmore [3] and the oxygen isotope substitution experiment of Zeidler et al. [81, 82], our PI-DPMD simulations characterize these isotope effects as small and subtle changes to the microscopic structure of liquid water. In doing so, we find that our approach predicts systematically smaller isotope effects in liquid water compared with the experimental/EPSR procedure. This discrepancy can be attributed to the imperfect description of the underlying PES at the vdW-inclusive hybrid DFT (i.e., PBE0-TS) level of theory, or to the remaining uncertainty due to the experimental challenge. Looking ahead, this work suggests that the combination of ML-based techniques (such as the DPMD model) and a sufficient amount of highly accurate ab initio data (beyond the vdW-inclusive hybrid DFT level) could pave the way towards studying subtle physical effects that were computationally prohibited by the length and time scales accessible to standard AIMD and PI-AIMD simulations.

Acknowledgements

HYK, LZ, WE, and RC gratefully acknowledge support from the U.S. Department of Energy (DOE) under Grant No. DE-SC0019394. HW gratefully acknowledge support from the National Natural Science Foundation of China (NSFC) under Grant Nos. 11501039, 11871110, and 91530322, as well as the National Key Research and Development Program of China under Grant Nos. 2016YFB0201200 and 2016YFB0201203. LZ and WE also acknowledge support from the Major Program of the National Natural Science Foundation of China (NNSFC) under Grant Nos. 91130005 and U1430237, as well as the Office of Naval Research (ONR) Grant No. N00014-13-1-0338. RD acknowledges the support of the Center for Alkaline Based Energy Solutions (CABES), an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Award No. DE-SC0019445. This work used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.