Molecular Polarizability of Water from the Local Dielectric Response Theory

TL;DR

This paper introduces a fully ab initio local dielectric response theory to accurately compute water's molecular polarizability, revealing anisotropic effects and charge transfer contributions that enhance understanding of water's dielectric properties.

Contribution

It develops a novel ab initio method based on local dielectric response theory to analyze electronic excitations in water beyond the dipole approximation.

Findings

Hydrogen-bond network causes strong anisotropic effects on polarizability.

Charge transfer increases isotropic polarizability by 5-6%.

Out-of-plane component of polarizability is significantly enhanced.

Abstract

We propose a fully ab initio theory to compute the electron density response under the perturbation in the local field. This method is based on our recently developed local dielectric response theory [Phys. Rev. B 92, 241107(R), 2015], which provides a rigorous theoretical framework to treat local electronic excitations in extended systems beyond the commonly employed dipole approximation. We have applied this method to study the electronic part of the molecular polarizability of water in ice Ih and liquid water. Our results reveal that the crystal field of the hydrogen-bond network has strong anisotropic effects, which significantly enhance the out-of-plane component and suppress the in-plane component perpendicular to the bisector direction. The contribution from the charge transfer is equally important, which increases the isotropic molecular polarizability by 5-6%. Our study provides new insights into the dielectric properties of water, which form the basis to understand the electronic excitations in water and to develop accurate polarizable force fields of water.

Tables

Table 1. Table 1: Molecular polarizability of water (in Å 3 ) of the gas phase monomer, ice I h and liquid water. Results in liquid water are the statistical average ( ⟨ α ⟩ delimited-⟨⟩ 𝛼 \langle\alpha\rangle ) and the standard deviation ( σ 𝜎 \sigma ). The x 𝑥 x axis is the dipole axis, and the z 𝑧 z axis is the normal direction of the molecular plane. Gas phase molecules in the rows of ice I h and liquid water have the same geometry as in the condensed phase. CF and CT denote the crystal field and charge transfer, respectively.

		${\alpha }_{xx}$	${\sigma }_{xx}$	${\alpha }_{yy}$	${\sigma }_{yy}$	${\alpha }_{zz}$	${\sigma }_{zz}$	${\overline{\alpha }}_{H_{2}O}$	$\sigma$
monomer	PBEa	1.58		1.60		1.58		1.59
	BLYPa	1.60		1.63		1.59		1.61
	BLYPb,c	1.47		1.53		1.42
	expd	$1.468\pm 0.003$		$1.528\pm 0.013$		$1.415\pm 0.013$		$1.470\pm 0.003$
	gas phase	1.65		1.71		1.61		1.66
ice Ih	CF	1.58		1.44		1.76		1.60
	CF+CT	1.69		1.54		1.88		1.70
	gas phase	1.64	0.05	1.69	0.08	1.61	0.03	1.65	0.05
water	CF	1.66	0.10	1.55	0.09	1.78	0.14	1.66	0.08
	CF+CT	1.76	0.11	1.63	0.10	1.86	0.16	1.75	0.09
	BLYPb	1.45		1.42		1.48
	BLYPc	1.44	0.03	1.41	0.03	1.49	0.03

Equations

$$|w_{\mathbf{R}n}⟩=\frac{\Omega }{(2\pi {)}^3}{\int }_{BZ}d\mathbf{k}{e}^{-i\mathbf{k}\cdot \mathbf{R}}\sum _{m=1}^{n_v}{U}_{mn}^{(\mathbf{k})}|{\psi }_m^{\mathbf{k}}⟩,$$

$$|w_n⟩=\sum _{m=1}^{n_v}{U}_{mn}|{\psi }_m⟩.$$

$${\mathbf{\mu }}_M=\sum _{i\in M}{Z}_i{\mathbf{R}}_i^{ion}-2\sum _{n\in M}{\mathbf{r}}_n,$$

$$x_n=-\frac{L}{2\pi }\Im \ln ⟨w_n|e^{-i\frac{2\pi }{L}x}|w_n⟩,$$

$$\Delta {\mathbf{\mu }}_M={\alpha }_M^{gas}\mathbf{E},$$

$$\Delta {\mathbf{\mu }}_M={\alpha }_M{\mathbf{E}}_M^{loc},$$

$${\mathbf{E}}_M^{ind}\approx \sum _{N\ne M}{T}_{MN}\Delta {\mathbf{\mu }}_N,$$

$$\Delta {\mathbf{\mu }}_M={\alpha }_M(\mathbf{E}+\sum _{N\ne M}{T}_{MN}\Delta {\mathbf{\mu }}_N).$$

$${\chi }^0(\omega ;\mathbf{r},{\mathbf{r}}^{\prime })$$

$$\chi (\omega ;\mathbf{r},{\mathbf{r}}^{\prime })$$

$$\Delta \rho$$

$$\Delta {\rho }_v$$

$$({\epsilon }_v-\hat{H}-\alpha \hat{P}\pm \omega )|\Delta {\psi }_v^{\pm }⟩=\hat{Q}\Delta V_{scf}|{\psi }_v⟩.$$

$$\Delta \rho =\sum _{n}2(|\Delta W_n^{+}⟩+|\Delta W_n^{-}⟩)⟨W_n|,$$

$$\sum _{n^{\prime }}({\tilde{\epsilon }}_{n{n}^{\prime }}-\hat{H}-\alpha \hat{P}\pm \omega )|\Delta W_{n^{\prime }}^{\pm }⟩=\hat{Q}\Delta V_{scf}|W_n⟩.$$

$$|\Delta W_n^{\pm }⟩$$

$$\equiv \sum _{n^{\prime }}[\tilde{\epsilon }-(\hat{H}+\alpha \hat{P}\mp \omega )\hat{I}{]}_{n{n}^{\prime }}^{-1}\hat{Q}\Delta V_{scf}|W_{n^{\prime }}⟩,$$

$$\Delta \rho$$

$$\Delta {\rho }_{n{n}^{\prime }}$$

$${\chi }_{n{n}^{\prime }}=\sum _{n^{\prime \prime }}\sum _{m=0}^{\infty }{\chi }_{n{n}^{\prime \prime }}^0{\left(K{\chi }^0\right)}_{n^{\prime \prime }{n}^{\prime }}^m.$$

$${\chi }_n={\chi }_n^0(1-K{\chi }^0{)}^{-1}.$$

$$\Delta {\rho }_M={\chi }_M^{us}\Delta V_{loc},$$

$${\chi }_M^{us}={\chi }_M^0+{\chi }_M^{0}K{\chi }_M^0+{\chi }_M^0(K{\chi }_M^0{)}^2+\cdots ={\chi }_M^0+{\chi }_M^{0}K{\chi }_M^{us}.$$

$${\alpha }_{M,ij}=-\int d{r}_{i}d{r}_j^{\prime }{r}_i{\chi }_M^{us}(\omega ;\mathbf{r},{\mathbf{r}}^{\prime })r_j^{\prime }\equiv -Tr\left[r_i({\chi }_M^{us}{r}_j^{\prime })\right],$$

$$\Delta {\rho }_E={\chi }_E^{us}\Delta V_{loc}.$$

$${\chi }_S^{us}$$

$${\chi }_S$$

$${\chi }_M$$

$${\chi }_E$$

$$\Delta {\mu }_{M,i}={\alpha }_{M,ij}{E}_j-Tr\left[r_i({\chi }_M^{us}K{\chi }_E{r}_j^{\prime })\right]E_j.$$

$$\Delta {\mu }_{M,i}={\alpha }_{M,ij}{E}_j+\sum _{N\in E}\sum _{kl}{\alpha }_{M,ik}{T}_{MN,kl}\Delta {\mu }_{N,l},$$

$${\alpha }_M^{eff}={\alpha }_M+\sum _{N\in E}{\alpha }_M{T}_{MN}{\alpha }_N^{eff}.$$

$$\left(\begin{array}{c}{\chi }_M\\ {\chi }_E\end{array}\right)={\left(\begin{array}{cc}1& -{\chi }_M^{us}K\\ -{\chi }_M^{us}K& 1\end{array}\right)}^{-1}\left(\begin{array}{c}{\chi }_M^{us}\\ {\chi }_E^{us}\end{array}\right).$$

Full text

I Introduction

Water is one of the most important substances for life, many fields of science, and numerous technological applications. Despite its simple molecular structure, water has many anomalous behaviors, e.g., its density reaching the maximum at 4°C, that have attracted intensive research for decades. One intriguing aspect of the fundamental properties of water is the electronic excitation, which is essential to the understanding of a broad range of problems, such as solvation, water/solid interfaces, and electrochemical reactions in liquid solution. The quantum theory of electronic excitation has been widely used to interpret, e.g., the x-ray absorption Cavalleri et al. (2002); Hetényi et al. (2004); Prendergast and Galli (2006); Chen et al. (2010) and optical absorption spectra of water Hahn et al. (2005); Garbuio et al. (2006), which in turn provide important physical insights into the atomic structures and the electronic structure of water. Besides, it has been well established that van der Waals dispersion forces, arising from the coupling of instantaneous and induced excitations, play a critical role in accurately describing the structure of water in ab initio molecular dynamics simulations Santra et al. (2008); Zhang et al. (2011); Santra et al. (2011); Wang et al. (2011); Gillan et al. (2016).

In the linear response regime, how to describe the electronic part of the molecular polarizability of water ($\alpha_{\rm{H_{2}O}}$) remains a subject of debate, despite the rather simple textbook picture of dielectrics Kittel (1966). Experimentally, $\alpha_{\rm{H_{2}O}}$ is known in the gas phase only Murphy (1977). In the condensed phase, the average value of the molecular polarizability is typically estimated from the refraction index measurement using the Lorentz-Lorenz relation Debye (1912). Although this approach is suitable for gases and nonpolar liquids, it fails for polar liquids, such as water, due to the inaccurate description of the local field experienced by the polar molecule. Many competing models have been proposed to improve Debye’s dipole theory Debye (1912), including models from Onsager Onsager (1936) and Kirkwood Kirkwood (1939). Due to the deficiency of a classical treatment of the local field in liquid water, this problem has not been solved yet. Besides the actual value of $\alpha_{\rm{H_{2}O}}$, it is more important to understand how $\alpha_{\rm{H_{2}O}}$ changes in different chemical environments (e.g. from the gas phase to the condensed phase) or under different boundary conditions (e.g., bulk water or confined water, such as water at the solid / water interface).

From the computational point of view, although the macroscopic average of the dielectric response of water is very well captured by the density functional theory (DFT) or beyond-DFT methods, such as the many-body perturbation theory Hahn et al. (2005); Garbuio et al. (2006), the chemical nature of the microscopic counterpart, e.g., the microscopic electric susceptibility ($\chi$), remains poorly understood. Currently, a rigorous, fully ab initio theory to compute $\alpha_{\rm{H_{2}O}}$ is lacking, and existing models predict contradicting trends in $\alpha_{\rm{H_{2}O}}$. Using point charge models to represent the solvent molecules, Gubskaya and Kusalik Gubskaya and Kusalik (2001) found an increase of $\alpha_{\rm{H_{2}O}}$ in liquid. Morita Morita (2002) proposed a cluster model, where the solute polarizability in the cluster is approximated by the difference of the polarizability of the cluster and that of the solvent only, $\Delta\alpha=\alpha^{\rm{tot}}-\alpha^{\rm{solv}}$. The non-additive correction was taken into account by introducing a dielectric continuum model of the solvation shell Morita (2002). In contrast to the point charge model, the cluster model predicts the isotropic molecular polarizability, $\bar{\alpha}_{\rm{H_{2}O}}=\frac{1}{3}\rm{Tr}\,(\alpha_{\mathrm{H_{2}O}})$, in the liquid is reduced from that in the gas phase by $7-9$%, and the reduction was attributed to the electron repulsion of the ambient solvent molecules that perturb and confine the spatially diffuse tail of the electron cloud of the solute Morita (2002).

Most of other computational studies employed an extension Heaton et al. (2006) of the interactive dipole model (IDM) (Silberstein, 1917a; *SILB17B; *SILB17C; Applequist et al., 1972; Thole, 1981) to the condensed phase, where the electrostatics of the electron density response is approximated at the dipole-dipole interaction level. IDM calculations found that $\bar{\alpha}_{\rm{H_{2}O}}$ in the liquid water is reduced by less than $2$% Salanne et al. (2008); Buin and Iftimie (2009) or the same as Wan et al. (2013) that of the gas phase. These results also showed an anisotropic effect, with the in-plane components reduced, while the out-of-plane component enhanced Salanne et al. (2008); Buin and Iftimie (2009); Wan et al. (2013). The origin of these changes is unclear, as the effects of the crystal field are intertwined with structural changes of water monomer geometries under the thermal fluctuation.

Because of the uncertainties in determining $\alpha_{\rm{H_{2}O}}$ from both experiment and theory, there is no clear recipe on how to choose the right $\alpha_{\rm{H_{2}O}}$ in a polarizable force field, given the knowledge of the gas phase value Ponder et al. (2010). To address this open question, in this work we proposed a fully ab initio method to compute $\alpha_{\rm{H_{2}O}}$ based on the our recently developed local dielectric response theory (LDRT) Ge and Lu (2015), which is among the techniques to compute the dielectric response function without explicitly referring to the empty states Baroni et al. (1987, 2001); Wilson et al. (2008, 2009); Walker et al. (2006); Rocca et al. (2008, 2010). The new method is general and can be used to study $\alpha_{\rm{H_{2}O}}$ in different chemical environments. We proved that the widely used IDM is the dipole limit of our theory. In the numerical study of $\alpha_{\rm{H_{2}O}}$ in liquid water, we paid special attention to separate the environmental effects due to the crystal field and charge transfer from those caused by intramolecular thermal fluctuations.

II Method

II.1 Molecular polarizability in extended systems under the dipole approximation

In order to quantify molecular polarizabilities in an extended system, it is necessary to express the dielectric response of a many-electron system in a local representation, instead of the Bloch representation. A formal way to proceed is to use the Wannier function (WF) formalism, such as the maximally localized Wannier function (MLWF) Marzari and Vanderbilt (1997); Souza et al. (2001); Marzari et al. (2012), where the WFs ($|w_{{\bf R}n}\rangle$) are constructed from unitary transformations of the Bloch orbitals ($|\psi_{m}^{\bf k}\rangle$). In systems with a finite band gap, we consider only the occupied bands:

[TABLE]

where $n_{v}$ is the total number of the occupied bands, and $\Omega$ is the real space primitive cell volume. Unitary matrices ($U^{({\bf k})}$) in the MLWF formalism minimize the spatial spreads of the WFs labeled by the lattice vector (${\bf R}$) and the Wannier center index ($n$). In this study, we focus on isolated systems or periodic systems with a large supercell, where a single $\Gamma$-point sampling in the reciprocal space is used. The extension to general cases of the $k$-point sampling, although in principle feasible, is beyond the scope of this work. In the $\Gamma$-point formulation, the above expression is simplified to

[TABLE]

Within the modern theory of the electric polarization in crystalline dielectrics King-Smith and Vanderbilt (1993); Vanderbilt and King-Smith (1993); Souza et al. (2000); Resta (1994), the dipole moment of a sub-system (e.g., a ion or molecule labeled by $M$) is defined in atomic units as

[TABLE]

where $Z_{i}$ and $\bm{R}_{i}^{ion}$ are the charges and positions of ions. Positions of the Wannier centers that belong to $M$, $\bm{r}_{n}=x_{n}\hat{i}+y_{n}\hat{i}+z_{n}\hat{k}$, can be computed from Marzari and Vanderbilt (1997)

[TABLE]

where $L$ is the size of the supercell.

Upon a small perturbation from an external electric field ($\bm{E}$) at a given frequency, the electronic contribution to the net induced dipole of $M$ in the gas phase is given by

[TABLE]

where $\alpha_{M}^{gas}$ is the frequency-dependent gas phase molecular polarizability tensor. Here we dropped the explicit frequency dependence to simplify the notation. Once the molecule is embedded in a chemical environment, e.g. a cation in a solid or a solute molecule surrounded by solvent molecules, the induced dipole arises from the response to the local field that is the sum of the external field and the induced field from the environment, $\bm{E}_{M}^{loc}=\bm{E}+\bm{E}_{M}^{ind}$. Unlike $\bm{E}$, $\bm{E}_{M}^{ind}$ is a microscopic quantity, whose direction and amplitude vary within the size of the molecule. In practice, this finite size effect is often ignored, which leads to the approximate expression,

[TABLE]

where $\alpha_{M}$ is the environment-dependent molecular polarizability tensor.

Under the dipole approximation, $\bm{E}_{M}^{ind}$ in the IDM is approximated by the dipole field of environment molecules ($N\neq M$),

[TABLE]

where $T_{MN}=\nabla\nabla r_{MN}^{-1}$ is the dipole-dipole interaction tensor with $r_{MN}$ the inter-molecular distances. Substituting Eq. 7 into Eq. 6, one can derive an equation of interacting dipoles (Silberstein, 1917a, b, c; Applequist et al., 1972; Thole, 1981; Heaton et al., 2006),

[TABLE]

$\alpha_{M}$ in Eq. 8 can be solved conveniently, once $\Delta\bm{\mu}$ are obtained from finite external field calculations Heaton et al. (2006); Buin and Iftimie (2009) or linear response theory Wan et al. (2013).

II.2 Local dielectric response theory

Although the IDM combined with DFT has been widely used to calculate molecular polarizabilities Heaton et al. (2006); Salanne et al. (2008); Buin and Iftimie (2009); Rotenberg et al. (2010); Wan et al. (2013), it has several drawbacks. First of all, as mentioned above, Eq. 6 neglects the finite size effect in the induced field of the environment. Secondly, the dipole approximation neglects all the higher order terms. It might be justifiable at the far field, but when molecules get close to each other, the validity of the dipole approximation is questionable. More severely, the use of $T_{MN}$ does not satisfy the positive definite requirement of Eq. 8 Thole (1981). As a consequence, molecular polarizabilities from the IDM can diverge or become negative in simple systems, e.g., the cooperative (head to tail) induced dipoles in the direction of the line connecting the two Thole (1981); Van Duijnen and Swart (1998); Allen (2004). Finally, at the interface between $M$ and the environment (see Fig. 1b), wavefunctions of excited environmental electrons can have finite overlap with the occupied valence electron orbitals of $M$, leading to a charge transfer (CT) type of density response on $M$. The CT contribution is a quantum effect, and can not be captured at the level of dipole-dipole interaction.

In order to develop a rigorous, fully quantum mechanical treatment of $\alpha_{M}$, let us first consider a slightly broader scenario, the local excitations of a molecule embedded a chemical environment as shown in Fig. 1a. The whole system can thus be divided into three regions (R${}_{1}\in$ R${}_{2}\in$ R3) as shown in Fig. 1b, R1: the region of the molecule of interest ($M$), R3: the whole system providing the crystal field (CF), and R2: a subset of R3 coupled to R1, such that CT can happen between R2 and R1 upon local excitations. We will introduce a formal theory based on the microscopic susceptibility, $\chi$. We note that Allen Allen (2004) developed a model of charge-dipole interaction along the same line, where charges are treated by quantum mechanics, and charge interactions are solved in the random phase approximation (RPA). The remaining obstacle in our approach is to reformulate $\chi$, a non-local quantity by definition, in the local representation.

An early proposal was made by Hanke to address this issue by localizing electrons and holes separately when building an explicit electron-hole pair basis Hanke (1973). In practice, this method is inconvenient, because of a) the difficulty to numerically converge the number of unoccupied bands, and b) even if the convergence can be reached, the poor locality of the high energy unoccupied bands. Both limitations of Hanke’s method Hanke (1973) can be avoided, because it is unnecessary to localize the unoccupied bands. In other words, we formulate the theory on the occupied manifold only, without using the explicit electron-hole pair basis. In the following, we first summarize the recently developed local dielectric response theory (LDRT) Ge and Lu (2015) that provides the theoretical framework to study local excitations. Then we use the LDRT to develop a fully ab initio theory to calculate the dielectric response of the perturbation in the local field.

The central quantities in the linear response theory are the bare, $\chi^{0}(\omega;\bm{r},\bm{r}^{\prime})$, and the screened susceptibility, $\chi(\omega;\bm{r},\bm{r}^{\prime})$, which are the functional derivatives of the charge density response with respect to perturbations in the self-consistent field potential ($\delta V_{scf}$) and the external potential ($\delta V_{ext}$):

[TABLE]

In the following, we adopt the shorthand notation: $\Delta\rho=\chi^{0}\,\Delta V_{scf}=\chi\,\Delta V_{ext}$, where integration on common variables is implied. $\chi$ can be solved from $\chi^{0}$ through Dyson’s equation Fetter and Walecka (2003), $\chi=\chi^{0}+\chi^{0}\,K\,\chi$, where $K=v_{H}+K_{xc}$ with $v_{H}$ and $K_{xc}$ being the Coulomb and exchange-correlation kernel, respectively. In the language of the linear response theory, $\chi^{0}$ is the building block.

Under the Bloch representation, both $\Delta\rho$ and $\chi^{0}$ can be partitioned according to occupied bands,

[TABLE]

where $|\Delta\psi_{v}^{\pm}\rangle$ are the solution of the Sternheimer equation at $+$ and $-$ frequencies Baroni et al. (2001),

[TABLE]

Here $\varepsilon_{v}$ are energy levels of the occupied states of the Kohn-Sham (KS) Hamiltonian ($\hat{H}$); $\hat{P}=\sum_{v}\left|\psi_{v}\right\rangle\langle\psi_{v}|$ and $\hat{Q}=\hat{I}-\hat{P}$ are projectors onto the occupied and unoccupied state manifolds, which are introduced to avoid the explicit reference to the unoccupied states Baroni et al. (2001). The term $\alpha\hat{P}$ is introduced to remove the singularity of Eq. 11.

The essence of the LDRT is to recognize that $\Delta\rho$ is invariant under the unitary transformation of the occupied orbitals, which allows $\Delta\rho$ and $\chi^{0}$ to be formulated under the Wannier representation,

[TABLE]

where $|\Delta W_{n}^{\pm}\rangle$ are solutions of the generalized Sternheimer equation Ge and Lu (2015),

[TABLE]

Since $|W_{n}\rangle$ are not the eigenstates of the KS Hamiltonian, $\varepsilon_{v}$ are replaced by the coupling matrix elements, $\tilde{\varepsilon}_{n,n^{\prime}}=\langle W_{n}|H|W_{n^{\prime}}\rangle$. In contrast to Eq. 11, linear response equations in Eq. 13 are entangled due to $\tilde{\varepsilon}_{n,n^{\prime}}$. The variation of $|W_{n}\rangle$ caused by the perturbation at $|W_{n^{\prime}}\rangle$ can be obtained from

[TABLE]

where $\hat{I}$ is an $N_{v}\times N_{v}$ identity matrix. The indices of $\Delta W_{nn^{\prime}}$ denote the perturbation site (right) and the response site (left), respectively. For systems with a finite band gap, Ge and Lu Ge and Lu (2015) proved that the spatial distribution of $|\Delta W_{nn^{\prime}}^{\pm}\rangle$ decays exponentially in real space, and its magnitude decays exponentially as the distance between sites $n$ and $n^{\prime}$.

Combining Eqs. 13 and 14, one can formally construct the partial response densities and partial microscopic susceptibilities (PMSs) on Wannier centers according to

[TABLE]

Similarly, $\chi$ can also be partitioned into PMSs through the Dyson’s equation,

[TABLE]

The partition of $\chi^{0}$ ($\chi$) into PMSs is illustrated in Fig. 2 for a system with two occupied Wannier centers. On each site, there are four possible diagrams: local perturbation (LP) only, local response (LR) only, local perturbation plus local response (LPR), and an empty site (ES). Therefore, $\chi^{0}$ ($\chi$) of the whole system can be partitioned into four terms: two diagonal terms (LPR on site 1 and 2) and two off-diagonal terms (LP on one site and LR on the other site).

PMSs like $\chi_{nn^{\prime}}$ are important quantities to study the non-locality of the response functions. In particular, $\chi_{nn^{\prime}}$ associates the response density at site $n$ to the external perturbation, $\Delta V_{ext}$, at site $n^{\prime}$. In the quantum chemistry language, similar notions have been introduced based on the local basis set formalism in the conceptual density functional theory Geerlings et al. (2003), to define the local reactivity index Torrent-Sucarrat et al. (2008), the atom-condensed linear response matrix Sablon et al. (2010), and a measure of aromaticity Sablon et al. (2012); Fias et al. (2013). It is worth noting that PMSs defined through the LDRT are formally additive, and do not suffer basis set-dependent errors in the local basis set approach. By summing over the perturbation index, one can construct local response density, $\Delta\rho_{n}=\sum_{n^{\prime}}\Delta\rho_{nn^{\prime}}$, and local susceptibilities, $\chi^{0}_{n}=\sum_{n^{\prime}}\chi^{0}_{nn^{\prime}}$ and $\chi_{n}=\sum_{n^{\prime}}\chi_{nn^{\prime}}$, or their symmetric form, $\tilde{\chi}^{0}_{n}=\frac{1}{2}\sum_{n^{\prime}}(\chi^{0}_{nn^{\prime}}+\chi^{0}_{n^{\prime}n})$ and $\tilde{\chi}_{n}=\frac{1}{2}\sum_{n^{\prime}}(\chi_{nn^{\prime}}+\chi_{n^{\prime}n})$. A Dyson-like relation exists between $\chi_{n}$ and $\chi_{n}^{0}$,

[TABLE]

$\chi_{n}$ provides a local measure of excited state properties projected onto a Wannier orbital, such as the bond polarizability Ge and Lu (2015).

II.3 Electronic density response of the local perturbation

We start by dividing the Wannier centers of the whole system into two sub-systems: the molecule of interest ($M$) and the environment ($E$) as shown in Fig. 1b. It is convenient to define local quantities associated with each sub-system, $S$ ($S$ = $M$ or $E$): $\Delta\rho_{S}=\sum_{n\in S}\Delta\rho_{n}$, $\chi^{0}_{S}=\sum_{n\in S}\chi^{0}_{n}$, and $\chi_{S}=\sum_{n\in S}\chi_{n}$. A quantum mechanical description of the molecular polarizability embedded in a medium relies on the density response to the perturbation in the local field, or in short, the local perturbation. Conceptually, within the linear response theory it implies a) to freeze both the external perturbation potential and the induced potential due to the polarization of the environment, and b) to obtain the self-consistent solution of the response density of the molecule. In the same spirit, Buin and Iftimie Buin and Iftimie (2009) proposed the frozen orbitals polarizability model within the IDM based on the MLWFs, as an alternative to Heaton et al.’s formulation Heaton et al. (2006). In practice, one has to “unscreen” (US) the PMSs on the molecule ($\chi_{S}$) to remove the screening effects from the environment. For this purpose, we define the unscreened molecular susceptibility ($\chi_{M}^{us}$) as

[TABLE]

where $\Delta V_{loc}=\Delta V_{ext}+K\,\Delta\rho_{M}$ is the perturbation in the local potential. Clearly, this quantity is different from $\chi$ that is related to external perturbation, $\Delta V_{ext}$, and $\chi^{0}$ that is related to the self-consistent field perturbation, $\Delta V_{scf}$. The unscreening procedure implies that

[TABLE]

The molecular polarizability is given by

[TABLE]

where $i$ and $j$ denote Cartesian axes. Because $\chi_{M}^{0}$ and $\chi_{M}$ are additive, so is the effective molecular polarizability, $\alpha_{M,ij}^{eff}=-Tr\left[r_{i}(\chi_{M}r^{\prime}_{j})\right]$. On the other hand, because $\chi_{M}^{us}$ is not additive, neither is $\alpha_{M}$.

In order to reveal the relation between the above quantum mechanical description and the IDM of Eq. 8, we similarly define $\chi_{E}^{us}$ as the unscreened susceptibility of the sub-system $E$ under the perturbation of the local field,

[TABLE]

It allows us to express the unscreened and screened susceptibilities of each sub-system,

[TABLE]

Employing the additivity relation, $\chi^{0}=\chi^{0}_{M}+\chi^{0}_{E}$, we have

[TABLE]

Consider $\Delta V_{ext}(\bm{r})=E_{j}\,r_{j}$. Multiply $-r_{i}$ from left and $E_{j}r^{\prime}_{j}$ from right on both sides of Eq. 23 and integrate. It follows that

[TABLE]

Ignore $K_{xc}$ in $K$, and take the dipole approximation to expand $v_{H}=\sum_{N\in E}|(\bm{r}_{M}+\bm{r}^{\prime\prime})-(\bm{r}_{N}+\bm{r}^{\prime\prime\prime})|^{-1}$ to the second order of $r_{MN}$. The only non-vanishing terms arise from $\bm{r}^{\prime\prime}\,T_{MN}\,\bm{r}^{\prime\prime\prime}$ due to the charge neutrality condition, $\int d\bm{r}\Delta\rho(\bm{r})=0$. Consequently, the second term on the right hand side of Eq. 25 becomes $\sum_{N\in E}\sum_{kl}\langle r_{i}|\chi_{M}^{us}|r^{\prime\prime}_{k}\rangle\,T_{MN,kl}\,\langle r^{\prime\prime\prime}_{l}|\chi_{E}|r^{\prime}_{j}\rangle\,E_{j}$, where $T_{MN,kl}=\partial_{k}\partial_{l}r_{MN}^{-1}$. It is straightforward to show that

[TABLE]

which reproduces Eq. 8. Alternatively, by dividing $E_{j}$ on both sides of Eq. 26, one obtains the IDM without explicit dependence on the external field,

[TABLE]

We have proved that the IDM is the classical limit of Eq. 23 under the dipole approximation.

One may also arrange the screened and unscreened PMSs in the matrix form through

[TABLE]

By expanding the matrix inversion in Taylor series, the first and second order corrections to $\chi_{M}$ are $\chi_{M}^{us}\,K\,\chi_{E}^{us}$ and $\chi_{M}^{us}\,K\,\chi_{E}^{us}\,K\,\chi_{M}^{us}$, respectively. In the limit that $M$ and $E$ are spatially fully separated, Dobson Dobson (1994) derived the same low order correction terms to $\chi_{M}$, which have been used by one of us to derive three-body terms in the RPA correlation energy Lu et al. (2010). Because Eq. 28 includes contributions at infinite orders, it is also valid in the regime where the electron densities of $M$ and $E$ overlap. We also note the similarity between Eq. 23 in this work and Eq. 41 in the linear response theory of subsystem TDDFT Pavanello (2013); Krishtal et al. (2015).

III Computational Details

In this study, ground state and linear response calculations were performed using the Perdew-Burke-Ernzerhof (PBE) Perdew et al. (1996); *PERD97 exchange-correlation functional with the SG15 optimized norm-conserving Vanderbilt (ONCV) pseudopotentials Hamann (2013); Schlipf and Gygi (2015); onc together with the $\Gamma$-point sampling. The kinetic energy cutoff of the planewave basis set was chosen at $65$ Rydberg. All the gas phase calculations were performed using a simple cubic supercell of $20$ Å. MLWFs were constructed with Wannier90 Mostofi et al. (2014). $\alpha_{M}$ is calculated from the self-consistent solution of $\chi_{M}^{us}$ according to Eq. 19 through the DFPT implemented in a customized version of Quantum ESPRESSO Giannozzi et al. (2009). At each self-consistent step, the generalized Sternheimer equations are solved, and the density response is projected onto $M$ as $\Delta\rho_{M}=\chi^{0}_{M}\Delta V_{loc}$. The converged density response is used to compute $\alpha_{M}$.

In practice, we first compute the ground state of the system, which can be either a finite system in a supercell or an extended system described under the periodic boundary condition. Since we are using a $\Gamma$-point formalism for extended systems, the unit cell size has to be sufficiently large. Next we construct Wannier orbitals for occupied states, and associate Wannier orbitals to their corresponding water molecules. Different environmental effects can be quantified separately based on the choices of R1, R2 and R3 denoted by $(n_{R_{1}}:n_{R_{2}}:n_{R_{3}})$ in the subsequent linear response calculations. Each of R1, R2 and R3 denotes a subset of water molecules. Specifically in the study of $\alpha_{\mathrm{H_{2}O}}$, R1 is restricted to one water molecule, i.e. $n_{R_{1}}=1$. For the CF effect, R3 contains the whole system, which means that $V_{scf}$ of the whole system is used in the Sternheimer equations. We restrict the charge transfer effects inside the subsystem defined by R2. If $n_{R_{2}}=n_{R_{3}}$, no truncation is applied, and the generalized Sternheimer equations are solved for all the occupied Wannier orbitals. If $n_{R_{2}}\in n_{R_{3}}$, the generalized Sternheimer equations are truncated so that only the occupied Wannier orbitals corresponding to R2 are included.

The effect of the crystal field is investigated by comparing water clusters with different sizes with extended systems. To this end, we consider two extended systems (ice Ih and liquid water) and three types of water clusters: water monomers ($n_{R_{3}}=1$), water clusters including the first solvation shell ($n_{R_{3}}=5$), and water clusters including the first and second solvation shells ($n_{R_{3}}=17$). The structures of gas phase water clusters are fixed at the same geometry as those in ice Ih or liquid water in order to eliminate the ambiguity due to the structural changes. To study the distance dependence of the CT effect, the region of R2 is varied in size, with the lower limit being one water molecule (R1) and the upper limit being the whole system (R3).

The ice Ih structure is modeled by an orthorhombic supercell ($n_{R_{3}}=96$; $a=13.30$ Å, $b=15.36$ Å, and $c=14.47$ Å) constructed from the hexagonal unit cell ($a=b=7.78$ Å, $c=7.33$ Å, and $\gamma=60$°) of Ref. Santra et al., 2013. Although the effect of proton disorder can in principle be studied by using different proton-ordered, energetically quasi-degenerate ice Ih structures, it is beyond the scope of the current work. The structures of liquid water were taken from the trajectories of the ab initio molecular dynamics simulation of $64$-molecule water samples at 400 K, the water PBE400 dataset pbe ($n_{R_{3}}=64$). The supercell size is $12.41$ Å. $\alpha_{\mathrm{H_{2}O}}$ are averaged over the first $20$ snapshots of the PBE$400$_$64$ subset, resulting in $1280$ water molecule geometries.

IV Results and Discussions

IV.1 Water monomer in the gas phase

$\alpha_{\mathrm{H_{2}O}}$ of the gas phase water monomer optimized with the PBE exchange-correlation functional is listed in Table 1. $\bar{\alpha}_{\mathrm{H_{2}O}}=1.59$ Å3, and the three principal components are $\alpha_{xx}=1.58$ Å3, $\alpha_{yy}=1.60$ Å3, and $\alpha_{zz}=1.58$ Å3, respectively. The principal axes are define as the following: $x$ along the bisector direction, $y$ along the in-plane perpendicular direction, and $z$ along the normal direction of the molecular plane. These values are in close agreement with the previous PBE results Wan et al. (2013). $\bar{\alpha}_{\mathrm{H_{2}O}}$ overestimates the experimental value by $8$%, in line with the error expected from PBE. This overestimation is associated with the underestimation of band gap and band width well-known for semi-local approximations to the exchange and correlation potential and in line with Penn’s model Penn (1962) of the dielectric constant. Nevertheless, our main focus is to understand the general trend of the environmental effects on $\alpha_{\rm{H_{2}O}}$ in the condensed phase, i.e., changes in $\alpha_{\rm{H_{2}O}}$, not the absolute value.

We also checked our gas phase results using the BLYP functional Becke (1988); Lee et al. (1988) in order to compare with the BLYP results in the literature Salanne et al. (2008); Buin and Iftimie (2009). As shown in Table 1, BLYP results in Refs. Salanne et al., 2008; Buin and Iftimie, 2009 are systematically smaller than our results, with $\bar{\alpha}_{\mathrm{H_{2}O}}$ is about $8$% smaller than our result. This difference is caused by the supercell size convergence issue in the finite field calculations Umari and Pasquarello (2003) in previous studies Salanne et al. (2008); Buin and Iftimie (2009). As pointed out by Buin and Iftimie Buin and Iftimie (2009), $\alpha_{xx}$ obtained from a $10$ Å cubic box is known to underestimate the fully converged value by $8$%, and similar underestimations are likely to occur also in their liquid water calculations. Their observation is consistent with the fact that $\alpha_{xx}$ in Refs. Salanne et al., 2008; Buin and Iftimie, 2009 is about $8$% smaller than our result. On the other hand, our linear response method does not suffer this convergence issue. For example, gas phase $\alpha_{ii}$ calculated with 20 Å and 30 Å cubic supercells are within 0.04%. If we compare PBE and BLYP results in this work, $\bar{\alpha}_{\mathrm{H_{2}O}}^{\mathrm{BLYP}}$ is slightly larger than $\bar{\alpha}_{\mathrm{H_{2}O}}^{\mathrm{PBE}}$ by $1$%.

IV.2 Isotropic molecular polarizability of water

Environmental effects on $\bar{\alpha}_{\mathrm{H_{2}O}}$ in ice Ih are shown in Fig. 3. Because the geometries of each water molecule in the ground state ice Ih model are almost identical, the computed $\bar{\alpha}_{\mathrm{H_{2}O}}$ differ by less than $0.15\%$. Therefore we present the results from a single water molecule in ice Ih. We consider four systems: three gas phase water clusters with $n_{R_{3}}=1$ (monomer), $5$ (first shell), and $17$ (second shell) and ice Ih with $n_{R_{3}}=96$. Here gas phase clusters refer to water molecules with the same geometry as in ice Ih to enable a straight comparison. First, we focus on the CF effect by fixing $n_{R_{2}}=1$ and varying $n=n_{R_{3}}$. The effect of the CF is to reduce $\bar{\alpha}_{\mathrm{H_{2}O}}$. Compared to $\bar{\alpha}_{\mathrm{H_{2}O}}=1.66$ Å3 of the gas phase monomer, $\bar{\alpha}_{\mathrm{H_{2}O}}$ decreases by $3.1$% caused by the CF of the first solvation shell. When the full CF of ice Ih is included, $\bar{\alpha}_{\mathrm{H_{2}O}}$ decreases by $3.8$%. Gubskaya and Kusalik Gubskaya and Kusalik (2001) reported an opposite trend, i.e., an increase of $\bar{\alpha}_{\mathrm{H_{2}O}}$ by $1.2\sim 5.3$%, which is likely due to the approximate nature of their local field models represented by the electrostatic field of a few point charges.

Next, we analyze the CT effect in ice Ih by fixing $n_{R_{3}}=96$ and varying $n=n_{R_{2}}$. In general, we expect the CT effect to enhance $\alpha_{M}$, as the excitation of neighboring water molecules has a constructive contribution. Indeed after we include the CT effect from the first solvation shell, $\bar{\alpha}_{\mathrm{H_{2}O}}$ increases from $1.60$ to $1.70$ Å3, which is more significant than the total CF effect. Since the CT effect is short-ranged, it saturates readily at $n_{R_{2}}=5$. At $n_{R_{2}}=96$, $\bar{\alpha}_{\mathrm{H_{2}O}}$ further increases by only $0.003$ Å3. Finally, we combine the CF and CT effects by simultaneously varying $n=n_{R_{2}}=n_{R_{3}}$. As shown in Fig. 3, the CT effect dominates within the first solvation shell, and $\bar{\alpha}_{\mathrm{H_{2}O}}$ reaches the maximum at $1.72$ Å3. Beyond that, the CT effect fades away, and the tail of the CF effect takes over. $\bar{\alpha}_{\mathrm{H_{2}O}}$ decreases slowly by $0.01$ Å3 at $n=96$. Overall, $\bar{\alpha}_{\mathrm{H_{2}O}}$ in ice Ih is larger than the gas phase value by $2.8$%.

The inset of Fig. 3 shows the isosurface of the electron density response in ice Ih due to the charge transfer under a uniform electric field along the $x$ axis (the red arrow), i.e., $\Delta\Delta\rho^{\rm{CT}}=\Delta\rho(1,96,96)-\Delta\rho(1,1,96)$. Clearly one can see a positive induced dipole moment emerging from the four nearest neighbor water molecules (highlighted in purple) that are H-bonded to the central water molecule. Previously, Lu et al. Lu et al. (2008) applied the MLWF procedure to localize the eigenvectors of the static dielectric matrices in ice Ih and liquid water, and identified dominant screening modes that are either localized on individual water molecules or involving H-bonded water dimers. Based on our study, the physical origin of these modes becomes clear, which are intramolecular excitations and intermolecular charge transfer excitations.

Next we compute $\bar{\alpha}_{\mathrm{H_{2}O}}$ in liquid water, where the thermal disorder effects play a key role in contrast to the results of ice Ih at zero temperature. Here we divided the thermal disorder effects into intramolecular contributions on individual water molecules and intermolecular contributions from the crystal field including H-bonds and long-range electrostatic effects. To isolate the intramolecular component, we sampled $\bar{\alpha}_{\mathrm{H_{2}O}}$ of gas phase water monomers with the same geometries as those in liquid water. The intermolecular contributions were extracted by comparing $\bar{\alpha}_{\mathrm{H_{2}O}}$ of the gas phase and liquid water that includes environmental effects.

As shown in Table 1 and Fig. 4, $\langle\bar{\alpha}_{\mathrm{H_{2}O}}\rangle$ is $1.65$ Å3 in the gas phase, which is averaged over $1280$ monomer configurations. Here $\langle\cdots\rangle$ denotes the statistical average. This mean value is nearly the same as that in ice Ih ($1.66$ Å3), and the intramolecular thermal disorder in the PBE water leads to a standard deviation of $\sigma=0.05$ Å3. However, once the CF effect is included, $\langle\bar{\alpha}_{\mathrm{H_{2}O}}\rangle$ slightly increases to $1.66$ Å3, which is substantially larger than that in ice Ih ($1.60$ Å3) by $4$%. At the same time, $\sigma$ increases to $0.08$ Å3, indicating a notable intermolecular contribution. The effect of the CT in liquid water is similar to that in ice Ih, which further increases $\bar{\alpha}_{\mathrm{H_{2}O}}$ by $0.09$ Å3 with $\sigma$ almost unaffected. This suggests that the crystal field is the primary source of the intermolecular thermal disorder effects.

IV.3 Anisotropic effects of the molecular polarizability of water

The crystal field in the condensed phases can modify the the electronic properties of individual water molecules significantly. For example, solvated water molecules form hydrogen-bonds (H-bonds) with their neighbors in a tetrahedral geometry, which has a strong influence on the Wannier function centers (WFCs: two from the O-H covalent bonds and two from the oxygen lone pairs). As a consequence, the O-H bond WFCs are pulled in and the lone pair WFCs are pulled out as compared to the gas phase monomer Silvestrelli and Parrinello (1999). In order to gain insights into the chemical origin of the CF and CT effects, we examine the anisotropic effects in ice Ih and liquid water by decomposing the $\alpha_{\mathrm{H_{2}O}}$ tensor onto three internal principal axes, and correlating the changes in the molecular polarizability ($\Delta\alpha_{ii},\,i=x,y,z$) with the changes in the O-WFC pair correlation function, $g(r)$, in terms of $\Delta r_{\mathrm{O-WFC}}$.

As shown in Table 1, the CF of ice Ih causes $\alpha_{xx}$ and $\alpha_{yy}$ to decrease by $4$% and $16$%, and $\alpha_{zz}$ to increase by $9$%. It appears that the in-plane suppression effect is more pronounced in the $y$ direction. Interestingly, this anisotropic CF effect can be qualitatively captured using the point charge model Gubskaya and Kusalik (2001) except for the $x$ component, which yields $\Delta\alpha_{xx}=4\%$, $\Delta\alpha_{yy}=-4\%$, and $\Delta\alpha_{zz}=16\%$ from their model II. This strong anisotropic CF effect is thus characteristic of the H-bond network in the condensed phase. On the other hand, the effect of the CT is mostly isotropic, increasing each component by about $0.1$ Å3.

Similar to the isotropic molecular polarizability, $\langle\alpha_{ii}\rangle$ of gas phase monomers in liquid water ($1.64$, $1.69$ and $1.61$ Å3) are very close to those in ice Ih ($1.65$, $1.71$ and $1.61$ Å3) as shown in Table 1 and Fig. 5, with the largest deviation of $1$% from the $y$ direction. The largest spread is also found in the $y$ direction ($\sigma_{yy}=0.08$ Å3), followed by $x$ ($\sigma_{xx}=0.05$ Å3) and $z$ ($\sigma_{zz}=0.03$ Å3). The CF in liquid water also exhibits a significant anisotropic effect, which leads to $\Delta\langle\alpha_{ii}\rangle=+1$%, $-9$%, and $+11$% as compared to the gas phase values. Another important feature is that the $z$ component acquires the largest spread ($\sigma_{zz}=0.14$ Å3), and the spreads in $x$ and $y$ directions are much smaller ($\sigma_{xx}=0.10$ Å3 and $\sigma_{xx}=0.09$ Å3). The CT effect in liquid water increases $\langle\alpha_{ii}\rangle$ almost uniformly by $0.1$ Å3, and the spreads are only slightly increased.

The origin of the CF effects in ice Ih and liquid water was investigated using $\Delta r_{\mathrm{O-WFC}}$, which are the differences of O-WFC distances between the condensed phases and the gas phase. As shown in Fig. 6, in ice Ih $\Delta r_{\mathrm{O-WFC}}=0.05$ and $-0.04$ Å for the oxygen lone pair and the O-H bond, respectively. In other words, O-H bond WFCs are pulled in from $0.53$ Å by $0.04$ Å, while lone pair Wannier centers are pulled out from $0.30$ Å by $0.05$ Å as compared to the gas phase monomer, in line with the established trend in liquid water Silvestrelli and Parrinello (1999). The formation of the H-bonds therefore weakens the intramolecular O-H bond, making it more ionic, and loosens the oxygen lone pairs. Since the O-H WFs become more tightly bound to the oxygen atom, the in-plane components of $\alpha_{\mathrm{H_{2}O}}$ are suppressed. In contrast, since the lone pair WFs become more loosely bound to the oxygen atom, the out-of-plane component is enhanced. In liquid water, both changes get smaller, i.e. $\Delta\langle r_{\mathrm{O-WFC}}\rangle=0.03$ and $-0.03$ Å, indicating overall softer H-bonds than ice Ih. Consequently, the CF effects have a smaller impact on the in-plane components of liquid water ($\Delta\langle\alpha_{xx}\rangle=1$% and $\Delta\langle\alpha_{yy}\rangle=-9$%) than ice Ih ($\Delta\alpha_{xx}=-4$% and $\Delta\alpha_{yy}=-16$%), while $\Delta\langle\alpha_{zz}\rangle$ in both systems are comparable (11% and 9%). The oxygen lone pair has a larger spread ($\sigma_{\rm{O-WFC}}=0.015$ Å) than the O-H bond ($\sigma_{\rm{O-WFC}}=0.012$ Å), which is consistent with the largest spread in $\alpha_{zz}$ (see Fig. 5 and Table 1).

To investigate any possible correlation between the CF and CT effects in liquid water, we track $\Delta\alpha_{ii}$ from CF and CT effects for each water molecule. As shown in Fig. 7, while the distribution of $\Delta\alpha_{xx}^{\mathrm{CF}}$ is nearly symmetric around zero, the signs of $\Delta\alpha_{yy}^{\mathrm{CF}}$ and $\Delta\alpha_{zz}^{\mathrm{CF}}$ are opposite, with one being predominantly negative and the other positive. No apparent correlation was found for the in-plane components, while there is a weak correlation between $\Delta\alpha_{zz}^{\mathrm{CF}}$ and $\Delta\alpha_{zz}^{\mathrm{CT}}$. Therefore our procedure of treating CF and CT effects separately is justified.

The main differences between our model and previous models can be understood as the following. In the point charge model Gubskaya and Kusalik (2001), the CF is over-simplified. In the IDM, although the full CF is included at the ground state level, the Coulomb kernel is approximated by the dipole interaction. It is not straightforward to compare the absolute values of $\alpha_{\rm{H_{2}O}}$ of this work with the IDM results in the literature Salanne et al. (2008); Buin and Iftimie (2009), because different exchange-correlation functionals (PBE and BLYP) were used, and results in previous studies are likely not fully converged. A likely meaningful comparison is the trend of the environmental effects on $\alpha_{\rm{H_{2}O}}$ with respect to the gas phase values computed from the optimized geometry. In this work, the relative changes of $\langle\alpha_{ii}\rangle$ are 11%, 6% and 18% in $x$, $y$, and $z$, respectively. The trends are qualitatively different from the IDM that yields -2$\sim$-1%, -8$\sim$-7% and 4$\sim$5%, respectively Salanne et al. (2008); Buin and Iftimie (2009). Besides the effects from different functionals, these differences may be attributed to several limitations in the IDM, such as the lack of the finite size effect and the use of the dipole-dipole approximation.

In summary, we have developed a fully ab initio method to compute the electron density response to the perturbation in the local field based on the local dielectric response theory. We applied this method to compute the molecular polarizability of water in condensed phases. Using the same molecular geometries as in the condensed phases, we found that the effects of the crystal field is to reduce $\alpha_{yy}$ and enhance $\alpha_{zz}$, and that the charge transfer effect increases all the principal components uniformly. Our study provides a rigorous theoretical framework to determine $\alpha_{\rm{H_{2}O}}$, essential to both the physical understanding of water and computer simulations using polarizable force field models. As the electron-based spectroscopy techniques, e.g., electron energy loss spectroscopy (EELS), has reached subangstrom spatial resolution Kapetanakis et al. (2015), we expect these experimental techniques can provide more details about the microscopic dielectric response of water in the future.

We thank Mark Hybertsen, Xifan Wu, Roberto Car, Annabella Selloni, Philip B. Allen and Michele Pavanello for their helpful discussions. This research used resources of the Center for Functional Nanomaterials, which is a U.S. DOE Office of Science Facility, at Brookhaven National Laboratory under Contract No. DE-SC0012704. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.