The molecular charge distribution, the hydration shell, and the unique properties of liquid water.

The molecular charge distribution, the hydration shell, and the unique properties of liquid water Ming-Liang Tan, Joseph R. Cendagorta, and Toshiko Ichiye Citation: The Journal of Chemical Physics 141, 244504 (2014); doi: 10.1063/1.4904263 View online: http://dx.doi.org/10.1063/1.4904263 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/24?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Effect of bound state of water on hydronium ion mobility in hydrated Nafion using molecular dynamics simulations J. Chem. Phys. 141, 104904 (2014); 10.1063/1.4894813 A recipe for free-energy functionals of polarizable molecular fluids J. Chem. Phys. 140, 144504 (2014); 10.1063/1.4870653 A first principles molecular dynamics study of lithium atom solvation in binary liquid mixture of water and ammonia: Structural, electronic, and dynamical properties J. Chem. Phys. 134, 024519 (2011); 10.1063/1.3511702 Relaxation of Voronoi shells in hydrated molecular ionic liquids J. Chem. Phys. 131, 174509 (2009); 10.1063/1.3256003 Can the dodecahedral water cluster naturally form in methane aqueous solutions? A molecular dynamics study on the hydrate nucleation mechanisms J. Chem. Phys. 128, 194504 (2008); 10.1063/1.2919558

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.225.1.226 On: Tue, 30 Dec 2014 15:48:24

THE JOURNAL OF CHEMICAL PHYSICS 141, 244504 (2014)

The molecular charge distribution, the hydration shell, and the unique properties of liquid water Ming-Liang Tan, Joseph R. Cendagorta, and Toshiko Ichiyea) Department of Chemistry, Georgetown University, Washington DC 20057, USA

(Received 29 September 2014; accepted 1 December 2014; published online 29 December 2014) The most essential features of a water molecule that give rise to its unique properties are examined using computer simulations of different water models. The charge distribution of a water molecule characterized by molecular multipoles is quantitatively linked to the liquid properties of water via order parameters for the degree (S 2) and symmetry (∆S 2) of the tetrahedral arrangement of the nearest neighbors, or “hydration shell.” ∆S 2 also appears to determine the long-range tetrahedral network and interfacial structure. From the correlations, some models are shown to be unable to reproduce certain properties due to the limitations of the model itself rather than the parameterization, which indicates that they are lacking essential molecular features. Moreover, since these properties depend not only on S 2 but also on ∆S 2, the long-range structure in these models may be incorrect. Based on the molecular features found in the models that are best able to reproduce liquid properties, the most essential features of a water molecule in liquid water appear to be a charge distribution with a large dipole, a large quadrupole, and negative charge out of the molecular plane, as well as a symmetrically ordered tetrahedral hydration shell that results from this charge distribution. The implications for modeling water are also discussed. C 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4904263]

I. INTRODUCTION

The unique properties of liquid water are generally attributed to the hydrogen (H) bonded tetrahedral network formed by water molecules.1,2 These properties include a temperature of maximum density (TMD), anomalous state dependence of the self-diffusion constant (D), and high dielectric permittivity (ε).3 Transport properties of supercooled water have been linked to the TMD,4 and anomalies in transport properties have been quantified using local translational and tetrahedral order parameters.5 Also, the unusual thermodynamic behavior has been linked to competition of bond and density order parameters.6 The molecular dipole and quadrupole have also been connected to ε.7,8 However, the most important molecular features responsible for the network are still not universally accepted despite years of intense study.1,9 Since experimental examinations of a water molecule in liquid water are difficult, computer simulations have played an important role. Moreover, since characteristics of water molecules cannot be varied experimentally, computer simulations can be used to query the molecular origins of liquid properties by determining the essential features in empirical models to reproduce them. Interestingly, relatively recent nonpolarizable models match certain properties surprisingly well with widely different molecular charge distributions10,11 and actually better than many models of the same vintage that include electronic polarizability.12,13 Conversely, the molecular charge distribution apparently affects the surface potential φ of the water-vapor interface,14 hydrophobic hydration,15 and possibly ion hydration.16 Several molecular features of water have been proposed as the minimal features necessary to describe a water a)[email protected]

0021-9606/2014/141(24)/244504/7/$30.00

molecule: the increased dipole in the liquid phase due to electronic polarization, a large quadrupole, and lone pairs (L).13 The larger dipole is generally accepted as important. However, since the four-site models have large quadrupoles but no lone pairs while the five-site models have sites representing tetrahedral lone pairs that cannot contribute to the quadrupole and thus have small quadrupoles, the predominance of nonpolarizable11 and polarizable17 four-site models implicitly assumes the importance of the large quadrupole18 over lone pairs. Also, in quantum mechanical/molecular mechanical calculations19–21 and ab initio molecular dynamics (AIMD) simulations,22 the molecular charge distribution of a water molecule in a liquid environment has the large quadrupole seen in the gas phase.23 On the other hand, comparisons of molecular dynamics (MD) with AIMD simulations also indicate that tetrahedral lone pairs improve over charge only in the molecular plane in the pure liquid24 and around ions,25 although the accuracy of AIMD studies of water has been limited by simple functionals and basis sets.26 Recent examination of MP2/aug-cc-pVQZ calculations shows that the large quadrupole is due not only to the charge on the hydrogens but also to the electron density, which is concentrated near the oxygen center with p orbital-type character perpendicular to the molecular plane (referred to here as p⊥ character), while tetrahedral lone pairs cannot contribute to the quadrupole.19 Thus, the electron density gives rise to a significant out-of-plane charge and the resulting cubic octupole is intermediate in value between the typical four- and five-site models. The success of very recent polarizable large quadrupole models such as the six-site TL6P27,28 model (polarized via a Gaussian inducible dipole) with negatively charged out-ofplane sites at p⊥ locations and the QDO29,30 and BK331 models with Gaussian negative charge density centered close to the

141, 244504-1

© 2014 AIP Publishing LLC


244504-2

Tan, Cendagorta, and Ichiye

oxygen (BK3 also has Gaussian positive charge density on the hydrogens) also indicates the importance of out-of-plane charge. Thus, quantifying the relationship between the molecular charge distribution and liquid properties can lead to a fundamental understanding of the importance of features such as the large quadrupole, out-of-plane charge, and electronic polarizability. In addition, determining the essential features is also important for improving potential energy functions for water, which are still widely used in computer simulations of the pure liquid5 and of solvent around ions,32 hydrophobes,33 and biomolecules34 given the current limitations of AIMD simulations for such systems.13 Our aim is to determine the most essential features of a water molecule in the liquid phase that give rise to its unusual properties, rather than to determine the best method for modeling these features or the best parameters for a method. The most essential features are defined here as the minimal features needed to reproduce liquid state properties, which are thus the features that a model should capture. Given that nonpolarizable models appear to mimic many properties quite well, the assumption is made that the average charge distribution, if modeled properly, is sufficient for mimicking pure liquid properties. As long as the assumption is correct, a model with a good average charge distribution for a given state point should be able to reproduce all non-electronic pure liquid properties near that state point; conversely, if no nonpolarizable model can reproduce them, then the ability to polarize must be important even at the state point where the model was parameterized. Specifically, while matching properties of liquid water can help verify a model, the inability to match certain properties can point to missing features. Of course, a nonpolarizable model will not hold far from where it was parameterized, and so that exact reproduction of the TMD is not to be expected even though the temperature dependence of the density has been used as a target for optimizing parameters for nonpolarizable models. However, since the average charge distribution does not change too much with temperature, a TMD close to experiment is expected. Here, the TMD, D, and ε of liquid water and φ of the watervapor interface are quantitatively linked to the charge distribution of a water molecule via the hydration shell using MD simulations. The charge distribution of a water molecule is described by its molecular multipoles and the hydration shell of a water molecule by two order parameters for the degree (S 2) and symmetry (∆S 2) of the tetrahedral arrangement of its nearest neighbors. The correlation between the order parameters and the molecular multipoles is first demonstrated, followed by a demonstration of the correlations between the order parameters and each of the liquid properties.

II. METHODS A. Water models

Water models were chosen from the literature that exhibit reasonable density and radial distribution functions. The water models examined were TIP3P,35 TIP3P-A,36 SPC/E,37 TIP4P,35 TIP4P-Ew,38 TIP4P/2005,39 TIP5P-0.5,40 TIP5P0.6,40 TIP5P-0.65,40 TIP5P,40 TIP5P-E,7 and SSDQO1.41

J. Chem. Phys. 141, 244504 (2014)

SSDQO2 uses the same van der Waals and electrostatic parameters as SSDQO1 except for Ω0 and Ω2 as noted in Table S1 in supplementary material.54 The three models labeled TIP5Px, where x is a fraction, have bOL = x Å except that x = 0.5 has bOL = 0.4875 Å, where bOL is the distance between the O and the lone pair site. These models all account for the average polarization of a water molecule in the liquid phase by a larger dipole given the robustness of previous results.10,11,13 The models can be divided into four types. First, the baseline behavior of point charges only on the nuclei is considered using three “simple models:” TIP3P-A, TIP3P, and SPC/E. Second, the large quadrupole (Q) is considered using three “large-Q” models with an additional negatively charged massless M site located at distance bOM from the oxygen towards the hydrogens: TIP4P, TIP4P-Ew, and TIP4P/2005. Third, tetrahedral lone pairs (LP) out of the molecular plane are considered using five “LP” models with an additional two negatively charged massless L sites located tetrahedrally at distance bOL from the oxygen: TIP5P-0.5, TIP5P-0.6, TIP5P-0.65, TIP5P, and TIP5P-E. The large-Q and LP models reduce to simple models if bOM or bOL, respectively, is zero. Finally, p⊥-type charge density is considered using two “p⊥” models: SSDQO1 and SSDQO2; a multipole model is used for efficiency since at least six charged sites are necessary for consistency with multipoles from quantum mechanical studies.19,27,42 For brevity, negative charge only in the molecular plane is referred to as two-dimensional (2D) while out-of-plane is referred to as three-dimensional (3D). Thus, simple models have small quadrupoles and 2D charge, large-Q models have large quadrupoles and 2D charge, LP models have small quadrupoles and 3D charge, and p⊥ models have both large quadrupoles and 3D charge. B. Simulation methods

The MD simulations were performed using the molecular mechanics package CHARMM version c36.43 The soft-sticky dipole-quadrupole-octupole (SSDQO) water model44 has been incorporated into CHARMM using MSCALE.45 The electrostatic interactions were calculated via the particle mesh Ewald method,46 with a β-spline coefficient equal to 6 and a k value of 0.34 Å−1. The SHAKE47 algorithm was used to constrain covalent bonds involving hydrogen atoms. All simulations employed the leapfrog Verlet algorithm48 with a time step of 1 fs and cubic periodic boundary conditions. The MD simulations were carried out in the NPT ensemble with 512 (360 for TMD studies) water molecules at a target temperature of 300 K and pressure of 1 atm using the Nosé-Hoover algorithm.49,50 Velocities were assigned according to a Gaussian distribution every 200 fs for 50 ps and then scaled every 1 ps if the temperature exceeded 300 K ± 5 K for 250 ps. The system was then allowed to run unrestrained for 9.5 ns and the last 6 ns of data was analyzed for each system. The coordinates at 1 ps intervals were analyzed for all average properties and 0.1 ps intervals for all time dependent properties. The translational diffusion coefficients D were calculated by a least-square fit of the mean-square displacement from simulation for t = 120-600 ps. The dielectric permittivities ε were calculated from the Kirkwood gK-factors51


244504-3


J. Chem. Phys. 141, 244504 (2014)

using Neumann’s approach.52 The surface potentials φ were from Ref. 53. C. Tetrahedral order parameters

To quantify the tetrahedrality of the hydration shell of a water molecule, two tetrahedral order parameters S 2 = (S 2H +S 2L)/2 and ∆S 2 = (S 2H −S 2L)/2 are proposed in terms of order parameters in the α directions, where α is approximately in either the H or so-called L directions, of a central water, given by S 2α =

  2 1 2

 P2(uα j · dα j ) .

(1)

j=1

P2 is the second order Legendre polynomial, d is a unit vector between the central water and the jth nearest neighbor, and u are unit vectors between the central water and the four corners of a perfect tetrahedron with two corners aligned in the molecular plane symmetrically about the hydrogens and the other two aligned symmetrically out-of-plane (Fig. 1). The four u vectors are exactly tetrahedral with respect to each other, each pointing in a α direction; an α neighbor lies in the α direction with one neighbor per direction. These order parameters are defined with respect to the molecular frame of the central molecule and the terms in the summation are independent of each other. While S 2 gives somewhat similar values to Debenedetti’s tetrahedral order parameter q,5 the main motivation for defining new order parameters is that an equivalent quantity to ∆S 2 cannot be defined for q, which is not defined with respect to the molecular frame. In addition, unlike Debenedetti’s q, the terms in the summation are independent of each other. S 2 = 1 for perfect tetrahedral order and S 2 = 0 for a random distribution, while ∆S 2 = 0 for symmetric order and ∆S 2 = ±1/2 for perfect order only around either the H (upper sign) or L (lower sign) directions with complete disorder in the opposite directions. The values of u·d are used to choose whether the water molecules are H or L neighbors and which of the four tetrahedral directions is used in the calculation of S 2. First, all possible u·d between the central water and the four nearest waters

FIG. 1. Definitions of u and d vectors. uH (blue) is in the hydrogen direction, uL (red) is in the lone pair direction, and d (green) is pointing to a neighbor in the H direction.

(i.e., 16 angles) are calculated. The water with the largest u·d among the four neighbors is assigned as a H or L neighbor based on whether d points to the H or L direction, respectively, and that u is used in the summation for S 2. Next, the procedure is repeated for the remaining three waters excluding the tetrahedral direction assigned in the previous step and again for the remaining two waters excluding the two tetrahedral directions assigned in the previous two steps. Finally, the remaining water is assigned to the remaining direction. D. Linear regressions

Fits of tetrahedral order parameters as a function of multipole moments. The fits of S 2 and ∆S 2 were to the coefficients of the 1/r 5 of the multipole interaction energy expansion, which is the lowest order term with sufficient angular dependence to give rise to tetrahedral order. In addition, the 1/r 5 dominates at contact distances, where higher power terms are still small. These fits are complicated because the moments are correlated for the site models. Therefore, rather than a single fit to all of the moments, the fits were performed in stages. The assumptions on the ratios of moments were based on SPC/E and a hypothetical model similar to SPC/E with the same dipole moment but with the negative charge placed on tetrahedral “lone pair” sites anti-symmetric to the hydrogen atoms (referred to as AST for anti-symmetric tetrahedral as in Ref. 19), for the same reasons as given above. The initial fit of S 2 assumed that asymmetry of the models was the cause of the largest possible variance. Thus, S 2 of the 3D models plus TIP3P (as a 3D model with bOL = 0) were fit for the slope and intercept of the dependence on µ0(Ω2 + 54 Ω0), with the ratio of Ω0 to Ω2 chosen such that this term disappears for SPC/E. Next, the size of the quadrupole was assumed to be the cause of the next largest variance. S 2 of all models were fit for the slope and intercept of the dependence on Θ22, with the assumption that the dependence on µ0(Ω2 + 54 Ω0) had the slope determined based on the fit of the 3D models. Dependence on the densities and Θ20 were also checked and found to be small. The constant S 2 (Ω0 = Ω2 = Θ2 = 0) = 0.52 seems reasonable since it should be somewhat greater than S 2 = 0.445, the value for point particles because of the limitation to the nearest tetrahedral direction, due to the finite size of the neighbors. The initial fit of ∆S 2 assumed that it accounts for all of the asymmetry between the H and L directions so there is no constant term. First, it was fit against only µ0Ω2 since all of the asymmetry between the H and L directions should lie in the energy terms with Ω2. However, it was clear that this did not account for all of the asymmetry since systematic errors could be seen for each type of model. Thus, it was fit against both µ0Ω2 and Θ22 + 98 µ0Ω0, which would still be zero for an AST type model, allowing the slopes for both to vary. Finally, a systematic variation was observed for the SSDQO-type models, which is not surprising since both the dipole-octupole and quadrupole-quadrupole terms are only calculated approximately. The final values of the coefficients for S 2 and ∆S 2 fit were obtained leaving out the SSDQO-type models. Fits of liquid properties as a function of order parameters. The linear regressions (LR) of the TMD and D to S 2 had two


244504-4


J. Chem. Phys. 141, 244504 (2014)

FIG. 2. Examples of the four types of models (charge shift indicated by blue ball, lone pairs, and p ⊥-density by pink balls) and their three-dimensional distribution of neighboring water molecules around a water molecule (contour levels 3.0 times bulk). (Graphics program courtesy of P. E. Mason and J. W. Brady.) See text for definitions.

parameters. The multiple linear regressions of gK and φ to S 2 and ∆S 2 had three parameters. The R2 are Pearson’s correlation coefficients between the predicted and actual values of the dependent variable in linear regression models that include an intercept.

III. RESULTS

The two tetrahedral order parameters S 2 = (S 2H + S 2L)/2 and ∆S 2 = (S 2H − S 2L)/2 are used to quantify the tetrahedrality of the hydration shells of water molecules in the simulations. The order parameters in each tetrahedral direction, S 2H and S 2L, calculated from simulations are consistent with threedimensional distributions of neighbors around a water molecule from the same simulations; the two H neighbors are localized in all of the water models while the two L neighbors show very model-dependent behavior (Fig. 2). However, experimental three-dimensional distributions have not been determined. The charge distribution of a water molecule should affect its hydration shell and thus S 2 and ∆S 2. Generally, the models with 2D charge have large ∆S 2 while the models with 3D charge have small ∆S 2, reflecting the degree of molecular symmetry between the location of charge in the H and L directions. In addition, the correlation of S 2H and S 2L with the charge distribution can be quantified via the 1/r 5 terms (i.e., quadrupole-quadrupole and dipole-octupole) of a multipole expansion of the electrostatic interaction energy, which have the strongest angular dependence at contact (i.e., with the nearest neighbors). Multiple linear regressions of S 2 and ∆S 2 with the spherical multipoles give S 2 ≈ 0.52 + 0.038Θ22 − 0.057µ0(Ω2 + 45 Ω0), ∆S 2 ≈ 0.030µ0Ω2 − 0.025(Θ22 + 98 µ0Ω0),

LP models as bOL approaches bOH, the O-H bond length. Thus, a large quadrupole or 3D charge contributes to more order in S 2, while 2D charge or a small quadrupole contributes to more asymmetry in ∆S 2. The effects of the tetrahedral order of the hydration shell on liquid properties can be quantified by their dependence on S 2 and ∆S 2. The values calculated for each model from the simulations are given in Table S2 of the supplementary material.54 The calculated TMD shows good correlation with a linear increase in S 2, with a linear regression 1 2 ∆TMD = TMD −TMD, exp ≈ −2TMD, exp 1 − 0.7 S , (3) where TMD,exp = 277 K with a correlation coefficient of R2 = 0.998 (Fig. 4(a), bottom). While it is not surprising that the results converge to the correct TMD since it was used in the optimization of some of these models, it is somewhat

(2)

with Pearson’s type coefficients of R2 = 0.9994 and 0.990, respectively, where µ0 is the dipole, Θ2 is the planar quadrupole, Ω0 is the linear octupole, and Ω2 is the cubic octupole (see Table S1 of supplementary material54). Plots of S 2H and S 2L from simulation show good correlation with those calculated by substituting multipoles of the model into the linear regression expressions (Fig. 3). The Θ22 terms increase in the large-Q models as bOM increases, while the µ0Ω2 terms decrease in the

FIG. 3. Correlation of the tetrahedral order parameter S 2 with multipoles. Top right corner is [S 2H]MD with [S 2H]LR (solid symbols) and bottom left corner is [S 2L]MD with [S 2L]LR (open symbols); one-to-one correspondence (dotted line), where the final subscript indicates from the MD or LR. Simple (brown and orange-tone triangles), charge shift (green-tone diamonds), lone pair (blue-tone squares), and p ⊥ (red-tone circles) models are identified in the legend (see supplementary material54). Increasing b OM and b OL is also indicated.


244504-5


J. Chem. Phys. 141, 244504 (2014)

FIG. 4. Simulation results plotted as differences from experiment for the (a) diffusion constant D (solid symbols, left axis), temperature of maximum density 2 TMD (open symbols, right axis), ( ) and linear regressions of D (dashed line) and TMD (dotted line) as functions of S ; (b) Kirkwood g-factor g K (solid symbols) 1 2 2 and (g K − 1) / 1 − 0.18 ∆S (open symbols) with the linear regression (dashed line) as functions of S ; and (c) surface potential φ (solid symbols) and 1 φ + 1.43(1 − 0.74 S 2) (open symbols) with the linear regression (dashed line) as function of ∆S 2. Symbol types are given in Fig. 3. Since only the experimental value of µ 02g K can be determined, µ 0 is assumed to be 2.3 D in the plots based on those models that give the experimental ε of 78. Also, since a range of experimental values (see supplementary material54) has been measured for φ, an average value of ∼0.1 eV is used.

surprising that they converge with the same dependence on S 2. In addition, since D decreases with the size of the diffusing entity and presumably a more tetrahedrally ordered hydration shell increases the effective size, the calculated D shows good correlation with a linear decrease in S 2, with a linear regression 1 2 ∆D = D − Dexp ≈ 9Dexp 1 − 0.7 S , (4) where Dexp = 2.32×10−9 m2/s with R2 = 0.985 (Fig. 4(a), top). Thus, both the TMD and D depend on the average (i.e., of the H and L directions) order of the hydration shell and both converge to the experimental values at S 2 ≈ 0.7 (Fig. 4(a)). The Kirkwood gK-factor,51 which is related to ε by gK = k BT(ε −1)(2ε +1)/(4π ρε µ20), is the correlation of the dipole vector of a central water with that of the rest of the water. The degree of positional order determines the fluctuations in this angle, so gK shows a linear decrease with S 2 (Fig. 4(b), solid symbols), but now the scaling differs by a multiplicative factor depending on whether the model has 2D or 3D charge. The models with 2D charge match the experiment for ε at S 2 ≈ 0.63, which is smaller than that needed to match the experiment for TMD and D, while the models with 3D charge match the experiment at S 2 ≈ 0.71 and thus can match the experimental TMD, D, and ε simultaneously. This indicates 3D charge is necessary to describe all three properties with the same model. In addition, the multiplicative factor scaling the dependence of gK on S 2 is apparently due to increased symmetry of the hydration shell leading to increased overall order, with a multiple linear regression of the gK for all models as a function of S 2 and ∆S 2 1 1 ∆S 2 1 − 0.75 S2 (5) gK ≈ 1 + 48 1 − 0.18 with Pearson’s type correlation coefficient of R2 = 0.996. Plots 1 2 of (gK − 1)/ 1 − 0.18 ∆S , which remove the contribution of ∆S 2 (Fig. 4(b), open symbols), illustrate how the larger ∆S 2 of the models with 2D charge reduces gK and thus ε compared to the models with 3D charge. Moreover, this means that the tetrahedral network predicted by 2D charge is more disordered, as

evidenced by the negligible contributions to gK from beyond the second shell in these models, while models with 3D charge can have significant contributions (Fig. 5). However, the first hydration shell of p⊥ models has smaller dipolar orientation than the LP models, similar to the large-Q models since both the p⊥ and large-Q models have large quadrupoles (Table S1 in supplementary material54). Thus, ∆S 2 also determines the range of the tetrahedral network since more symmetric

FIG. 5. Contribution to g K of four nearest neighbors (symbols to left of dashed line) and of nearest and next nearest neighbors (symbols to right of dashed line). When the next nearest neighbors are included, the simple (orange-tone triangles) and M-site (green-tone diamonds) models are close to the one-to-one correspondance line (dotted) with the total g K while the L-site (blue-tone squares) and p ⊥ (red-tone circles) models still have not reached the total g K, indicating contributions beyond the next nearest neighbors and thus longer range order. Note also that the p ⊥ models have somewhat smaller contributions from the first shell, like the M-site models, because both have large quadrupoles.


244504-6


J. Chem. Phys. 141, 244504 (2014)

hydration shells can form longer range networks, and more symmetric hydration shells indicated by small ∆S 2 give the best agreement with experimental ε. The surface potential of the water-vapor interface measures the orientation of the water dipoles at the interface. Since orientational effects are due to the asymmetry of the water molecule, φ correlates directly with symmetry ∆S 2 of the hydration shell in the pure liquid (Fig. 4(c), solid symbols). However, the order, S 2, of the hydration shell also affects the degree of orientation, with a multiple linear regression φ as a function of ∆S 2 and S 2, 1 S 2) φ ≈ 4.38∆S 2 − 1.43(1 − 0.74

(6)

1 with R2 = 0.983. Plots of φ + 1.43(1 − 0.74 S 2) which remove 2 the contribution of S (Fig. 4(c), open symbols) illustrate the effects of disorder. As for ε, the best agreement with the experiment is for small ∆S 2, i.e., for models with 3D charge.

IV. DISCUSSION

The correlation of the tetrahedral order parameters of the hydration shells from simulation with the multipole moments of the water model used in the simulation shows that the nature of the hydration shell is linked to the charge distribution of a water molecule. In addition, the correlation of the simulated liquid properties studied here with the simulated tetrahedral order parameters shows that these liquid properties are linked to the nature of the hydration shell, and thus to the charge distribution of a water molecule. From these correlations, 2D models are shown to be unable to reproduce certain properties due to the limitations of the model itself rather than the parameterization, which indicates that they are lacking essential molecular features. In addition, the structure of water predicted by models with 2D charge and those with 3D charge is significantly different, because of differences in gK and φ for the two types of models. Specifically, the persistence of contributions to gK beyond the second shell in models with 3D charge indicates that the tetrahedral network extends further than in models with only 2D charge. The structure predicted by the 3D models appears more correct since the 2D charge models are unable to predict the correct ε. Of the models with 3D charge, one LP model and the two p⊥ models appear to give equally good agreement with the experiment. However, in addition to the inconsistency of lone pairs with quantum mechanical results,19–22 the LP model has density ρ = 0.979 gm/cc at standard temperatures and pressures while the two p⊥ models have ρ = 0.993 and 0.997 gm/cc while the experimental value is ρ = 0.997 gm/cc. A re-parameterization of the lone pair model with the same electrostatic parameters (partial charges) but different LennardJones parameters has ρ = 1.003 gm/cc closer to experimental density, but gives poorer agreement with the other liquid properties examined here (Fig. 4, Table S2 of supplementary material54). Thus, both a large quadrupole and 3D charge, which is manifested in p⊥ charge, appear necessary to give good agreement with experimental liquid properties, although further investigation is warranted.

These studies also have implications for modeling water, which should be considered in light of the need for accuracy versus computational efficiency. For instance, the correct average charge distribution due to electronic polarization appears to be more important than including electronic polarizability for pure water at a given state point, although accurate polarizability of a good charge distribution would be important for modeling phase space behavior and in inhomogeneous environments. Specifically, the inability of nonpolarizable models to reproduce certain properties such as ε appears to be due to the lack of a good description of 3D charge rather than polarizability, as often proposed.11 One indicator of a good description of the 3D charge is an octupole intermediate between a five-site model and a four-site model based on quantum mechanical/molecular mechanical studies.19 A further implication is that a nonpolarizable model with a good average distribution may perform better than a polarizable model, if the polarization is compensating for a poor charge distribution. However, for changes in liquid properties as a function of temperature or water around solutes or at interfaces, polarizable models should perform better than the best nonpolarizable models since the charge distribution can be expected to change from its average value where a nonpolarizable model was parameterized. A test for a good polarizable model is the prediction of the correct average charge distribution over a range of states, for instance, between the gas and liquid phase. In particular, since both 3D charge and a large quadrupole are present in both the gas and liquid phase, they should be present in polarizable models in both phases as well, as suggested by the success of recent polarizable models such as QDO,29,30 BK3,31 and TL6P27,28 that have out-of-plane negative charge. V. CONCLUSIONS

Overall, the results show that the charge distribution of a water molecule determines its hydration shell, and that the degree and symmetry of the tetrahedral order of the hydration shell determines its liquid state properties. By quantifying the order of the hydration shell via multipoles, liquid properties (TMD, D, ε, and φ) are shown to be functions of the molecular multipoles, although limitations of the site representation, quantum effects, and charge exchange may affect the exact numbers. The results here indicate that a large dipole, a large quadrupole, and 3D charge consistent with p⊥ density are the most important features of a water molecule in the liquid phase and that symmetric tetrahedral order is an important characteristic of liquid water. These results also have implications for the molecular understanding of water and aqueous solutions that were based on computer simulations that have mainly used models with only 2D charge. For instance, since symmetric order gives rise to longer range order in the tetrahedral network, as indicated by the correlation of gK with ∆S 2 shown here, liquid water near the TMD may be more ice-like than thought based on large-Q models, consistent with the observation that large-Q models have melting temperatures Tm significantly too low with respect to the TMD while an LP model has both Tm and TMD in good agreement with the experiment.11 In addition,


244504-7


since symmetric tetrahedral order affects orientation of water at interfaces, as indicated by the correlation of φ with ∆S 2, water may be more tetrahedral around hydrophobic moieties than generally thought, which is consistent with observations in simulations using simple, large-Q, and p⊥ models for ethanol-water mixtures, in which the p⊥ model gives much better agreement with experimental results for density and neutron diffraction coordination numbers as a function of ethanol concentration, as well as spectroscopic data at low concentrations.15 Therefore, a molecular understanding of the unique properties of water lies not only in the nuclei but also the electron density distribution in a water molecule. ACKNOWLEDGMENTS

The authors are grateful for support from the National Science Foundation through Grant No. CHE-1158267 and from the McGowan Foundation. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. OCI-1053575; the Matrix and Medusa clusters, maintained by University Information Services at Georgetown University; and the Lobos cluster at the Laboratory for Computational Biology, National Heart, Lung, and Blood Institute, National Institutes of Health through the generosity of Bernard R. Brooks. T.I. also thanks Glenn J. Martyna and Thomas L. Beck for helpful comments. 1P.

Ball, Nature 452, 291 (2008). H. Poole, F. Sciortino, T. Grande, H. E. Stanley, and C. A. Angell, Phys. Rev. Lett. 73, 1632 (1994). 3D. Eisenberg and W. Kauzmann, The Structure and Properties of Water, 1st ed. (Oxford University Press, New York, 1969). 4C. A. Angell, Nature 331, 206 (1988). 5J. R. Errington and P. G. Debenedetti, Nature 409, 318 (2001). 6H. Tanaka, Phys. Rev. Lett. 80, 5750 (1998). 7S. W. Rick, J. Chem. Phys. 120, 6085 (2004). 8S. L. Carnie and G. N. Patey, Mol. Phys. 47, 1129 (1982). 9J. D. Bernal and R. H. Fowler, J. Chem. Phys. 1, 515 (1933). 10T. Ichiye, Adv. Chem. Phys. 155, 161 (2014). 11C. Vega and J. L. F. Abascal, Phys. Chem. Chem. Phys. 13, 19663 (2011). 12B. Guillot, J. Mol. Liq. 101, 219 (2002). 13W. L. Jorgensen and J. Tirado-Rives, Proc. Natl. Acad. Sci. USA 102, 6665 (2005). 14S. M. Kathmann, I.-F. W. Kuo, C. J. Mundy, and G. K. Schenter, J. Phys. Chem. B 115, 4369 (2011). 15M.-L. Tan, J. R. Cendagorta, and T. Ichiye, J. Am. Chem. Soc. 135, 4918 (2013). 16T. L. Beck, Chem. Phys. Lett. 561, 1 (2013). 17G. Lamoureux, A. D. Mackerell, Jr., and B. Roux, J. Chem. Phys. 119, 5185 (2003). 2P.

J. Chem. Phys. 141, 244504 (2014) 18C.

Vega and J. L. F. Abascal, J. Chem. Phys. 123, 144504 (2005). Niu, M.-L. Tan, and T. Ichiye, J. Chem. Phys. 134, 134501 (2011). 20K. Coutinho, R. C. Guedes, B. J. C. Cabral, and S. Canuto, Chem. Phys. Lett. 369, 345 (2003). 21A. Osted, J. Kongsted, K. V. Mikkelson, P.-O. Åstrand, and O. Christiansen, J. Chem. Phys. 124, 124503 (2006). 22P. L. Silvestrelli and M. Parrinello, J. Chem. Phys. 111, 3572 (1999). 23J. Verhoevan and A. Dymanus, J. Chem. Phys. 52, 3222 (1970). 24M. V. Fernández-Serra and E. Artacho, Phys. Rev. Lett. 96, 016404 (2006). 25R. C. Remsing, M. D. Baer, G. K. Schenter, C. J. Mundy, and J. D. Weeks, J. Phys. Chem. Lett. 5, 2767 (2014). 26M. D. Ben, M. Schönherr, J. Hutter, and J. VandeVondele, J. Phys. Chem. Lett. 4, 3753 (2013). 27P. Tröster and P. Tavan, J. Phys. Chem. Lett. 5, 138 (2014). 28P. Tröster, K. Lorenzen, and P. Tavan, J. Phys. Chem. B 118, 1589 (2014). 29A. Jones, F. Cipcigan, V. P. Sokhan, J. Crain, and G. J. Martyna, Phys. Rev. Lett. 110, 227801 (2013). 30V. P. Sokhan, A. Jones, F. S. Cipcigan, J. Crain, and G. J. Martyna, “Signature properties of water: Their molecular electronic origins,” Proc. Natl. Acad. Sci. USA (submitted). 31P. T. Kiss and A. Baranyai, J. Chem. Phys. 138, 204507 (2013). 32J. L. Skinner, Science 328, 985 (2010). 33D. Chandler, Nature 437, 640 (2005). 34A. A. Golosov and M. Karplus, J. Phys. Chem. B 111, 1482 (2007). 35W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983). 36D. L. Price and C. L. Brooks III, J. Chem. Phys. 121, 10096 (2004). 37H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987). 38H. W. Horn, W. C. Swope, J. W. Pitera, J. D. Madura, T. J. Dick, G. L. Hura, and T. Head-Gordon, J. Chem. Phys. 120, 9665 (2004). 39J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005). 40M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112, 8910 (2000). 41J. A Te and T. Ichiye, J. Chem. Phys. 132, 114511 (2010). 42W. Yu, P. E. M. Lopes, B. Roux, and A. D. Mackerell, Jr., J. Chem. Phys. 138, 034508 (2013). 43B. R. Brooks, C. L. Brooks III, A. D. MacKerell, Jr., L. Nilsson, R. J. Petrella, B. Roux, Y. Won, G. Archontis, C. Bartels, S. Boresch, A. Caflisch, L. Caves, Q. Cui, A. R. Dinner, M. Feig, S. Fischer, J. Gao, M. Hodoscek, W. Im, K. Kuczera, T. Lazaridis, J. Ma, V. Ovchinnikov, E. Paci, R. W. Pastor, C. B. Post, J. Z. Pu, M. Schaefer, B. Tidor, R. M. Venable, H. L. Woodcock, X. Wu, W. Yang, D. M. York, and M. Karplus, J. Comput. Chem. 30, 1545 (2009). 44T. Ichiye and M.-L. Tan, J. Chem. Phys. 124, 134504 (2006). 45H. L. Woodcock, B. T. Miller, M. Hodoscek, A. Okur, J. D. Larkin, J. W. Ponder, and B. R. Brooks, J. Chem. Theory Comput. 7, 1208 (2011). 46S. E. Feller, R. W. Pastor, A. Rojnuckarin, S. Bogusz, and B. R. Brooks, J. Phys. Chem. 100, 17011 (1996). 47J. P. Rychaert, G. Ciccotti, and H. J. C. Berendsen, J. Comput. Phys. 23, 327 (1977). 48L. Verlet, Phys. Rev. 165, 201 (1968). 49S. Nosé, J. Chem. Phys. 81, 511 (1984). 50W. G. Hoover, Phys. Rev. A: At Mol Opt. Phys. 31, 1695 (1985). 51J. G. Kirkwood, J. Chem. Phys. 4, 592 (1936). 52M. Neumann, J. Chem. Phys. 82, 5663 (1985). 53J. R. Cendagorta and T. Ichiye, “The Surface Potential of the Water-Vapor Interface from Classical Simulations,” J. Phys. Chem. B (to be published). 54See supplementary material at http://dx.doi.org/10.1063/1.4904263 for tables of multipole moments and liquid properties for the models studied here and experiment. 19S.


Hydration properties of lanthanoid(III) carbonate complexes in liquid water determined by polarizable molecular dynamics simulations.

The SAMPL4 hydration challenge: evaluation of partial charge sets with explicit-water molecular dynamics simulations.

The hydration shell of myoglobin.

Spatial and Orientational Structure of the Hydration Shell of Benzene in Sub- and Supercritical Water.

Liquid but durable: molecular dynamics simulations explain the unique properties of archaeal-like membranes.

Non-linear infrared spectroscopy of the water bending mode: direct experimental evidence of hydration shell reorganization?

Hydration shell effects in the relaxation dynamics of photoexcited Fe-II complexes in water.

Hydration of the cyanide ion: an ab initio quantum mechanical charge field molecular dynamics study.

Water dynamics in protein hydration shells: the molecular origins of the dynamical perturbation.

Shape dependence of the radial distribution function of hydration water around proteins.

Influence of the shell thickness and charge distribution on the effective interaction between two like-charged hollow spheres.

A Theoretical Study of the Hydration of Methane, from the Aqueous Solution to the sI Hydrate-Liquid Water-Gas Coexistence.

Introducing Charge Hydration Asymmetry into the Generalized Born Model.

Water is a molecular liquid.

Adaptive resolution simulation of a biomolecule and its hydration shell: Structural and dynamical properties.

vapor interface.

Ion hydration and associated defects in hydrogen bond network of water: observation of reorientationally slow water molecules beyond first hydration shell in aqueous solutions of MgCl2.

Water structure and hydration.

Molecular dynamics simulation of the hydration shell of a B-DNA decamer reveals two main types of minor-groove hydration depending on groove width.

Hydrogen Bonding and Related Properties in Liquid Water: A Car-Parrinello Molecular Dynamics Simulation Study.

Hydration control of the mechanical and dynamical properties of cellulose.

Metastable liquid-liquid transition in a molecular model of water.

The Hydration Equation: Update on Water Balance and Cognitive Performance.

Water Intake and Hydration Indices in Healthy European Adults: The European Hydration Research Study (EHRS).