Molecular Multipole Potential Energy Functions for Water - American

Nov 12, 2015 - Recently, models that utilize a multipole expansion located on a single site in the ... by its electrostatic multipole moments, which d...
0 downloads 0 Views 943KB Size
Download PDF
Article pubs.acs.org/JPCB

Molecular Multipole Potential Energy Functions for Water Ming-Liang Tan,†,‡ Kelly N. Tran,† Frank C. Pickard, IV,‡ Andrew C. Simmonett,‡ Bernard R. Brooks,‡ and Toshiko Ichiye*,†,‡ †

Department of Chemistry, Georgetown University, Washington, DC 20057, United States Laboratory of Computational Biology, National Heart, Lung, and Blood Institute, National Institutes of Health, Rockville, Maryland 20892, United States

‡

ABSTRACT: Water is the most common liquid on this planet, with many unique properties that make it essential for life as we know it. These properties must arise from features in the charge distribution of a water molecule, so it is essential to capture these features in potential energy functions for water to reproduce its liquid state properties in computer simulations. Recently, models that utilize a multipole expansion located on a single site in the water molecule, or “molecular multipole models”, have been shown to rival and even surpass site models with up to five sites in reproducing both the electrostatic potential around a molecule and a variety of liquid state properties in simulations. However, despite decades of work using multipoles, confusion still remains about how to truncate the multipole expansions efficiently and accurately. This is particularly important when using molecular multipole expansions to describe water molecules in the liquid state, where the short-range interactions must be accurate, because the higher order multipoles of a water molecule are large. Here, truncation schemes designed for a recent efficient algorithm for multipoles in molecular dynamics simulations are assessed for how well they reproduce results for a simple three-site model of water when the multipole moments and Lennard-Jones parameters of that model are used. In addition, the multipole analysis indicates that site models that do not account for out-of-plane electron density overestimate the stability of a non-hydrogen-bonded conformation, leading to serious consequences for the simulated liquid.

■

octupole.7 Furthermore, a recent comparison of 13 different nonpolarizable water models12 shows that a charge distribution with accurate multipoles in comparison to quantum mechanical calculations gives rise to a symmetric, ordered hydration shell and the most accurate properties in liquid state simulations, indicating that both the hydrogens and out-of-plane electron density are important molecular features. Thus, empirical potential energy functions for water must accurately represent the charge distribution due to the nuclei and the electron density of a water molecule. However, most models of water are based on the interaction energies between water molecules. To improve the models, a common approach has been to match various thermodynamic properties, which has led to models such as TIP4P-Ew,13 TIP4P/2005,14 and TIP5P.15 Another approach for improving interaction energy based models has been to model the functional form of the interaction energy based on ab initio potential energies, which has led to models such as TTM3-F,16 DPP2,17 and MB-Pol.18 This approach has led to very accurate results especially for the MB-Pol model; however, the potential energy form is complex and so far impractical for biological simulations. Only recently have a few models been proposed on the basis of the

INTRODUCTION Water has many unique properties as a pure liquid and as a solvent, which make it necessary for life on earth. However, although water is the most common liquid on this planet, the molecular bases for its properties as a liquid are still not understood well. The tetrahedral network of hydrogen bonds between molecules is thought to be responsible,1,2 but the essential features of a water molecule that give rise to this network and the actual structure of the network are still unclear.1 Therefore, water is unsurprisingly difficult to model using empirical potential energy functions in computer simulations. Furthermore, although ab initio molecular dynamics (AIMD) simulations hold promise, they have reached neither the accuracy nor the computational efficiency necessary for most problems involving water. Thus, empirical potential energy functions remain the most practical method for typical simulations.3 The molecular features of a water molecule can be quantified by its electrostatic multipole moments, which describe the charge distribution. A water molecule in the gas phase has been shown in experiments4,5 and quantum mechanical calculations6 to have relatively large dipole and quadrupole moments. In addition, quantum mechanical/molecular mechanical (QM/ MM)7−9 and AIMD10,11 simulations indicate that a molecule in the liquid phase has an even larger dipole and a somewhat larger quadrupole. A detailed study shows that the large quadrupole is due to both the hydrogens and the electron density from the p-type orbital on the oxygen perpendicular to the molecular plane, which gives rise to a medium-size © XXXX American Chemical Society

Special Issue: Bruce C. Garrett Festschrift Received: September 30, 2015 Revised: November 11, 2015

A

DOI: 10.1021/acs.jpcb.5b09565 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

represent the electrostatic potential of water molecules even at distances corresponding to the contact distance between water in the liquid phase. In particular, the octupole has long been understood to be necessary for water around ions,28 because it describes the asymmetry between the positive charge in the plane of the water molecule due to the hydrogens and the negative charge out of the molecular plane due to the electron density distribution.7 On the other hand, including higher order multipoles that only refine the energy in regions of phase space that are not often explored may reduce the efficiency unnecessarily. In the liquid phase, these regions could be high energy, for instance, unfavorable in the dimer interaction energy,6 or less sampled because other conformations are more favored in the hydrogen-bonded network of the liquid phase. Unfortunately, the correct distribution for liquid water is not known and the favorable conformations found in the first hydration shell were found to be different in simulations of SPC/E, TIP3P, TIP4P, and ST2.29 Since the interaction energy between multipoles is also an expansion, it is also necessary to truncate this expansion accurately and efficiently. Since higher order multipoles become more important at small distances from the multipole center r because of the inverse dependence on r, this is especially important at contact distances, which is essential in simulations of the liquid state. In particular, as r decreases, higher powers of 1/r become important, and for a given power, it is essential to include all large contributions due to the multipoles. Of course, as in any discrete representation of diffuse charge including partial charge descriptions, charge penetration effects at small distances cannot be accounted for30 except effectively by modification of the parameters. For instance, the expansion used in “soft−sticky dipole−quadrupole−octupole” (SSDQO)type models has terms up to 1/r5 term so that dipole−dipole, dipole−quadrupole, quadrupole−quadrupole, and dipole−octupole (where 2l1-pole-2l2-pole refers to both 2l1-pole-2l2-pole and 2l2-pole-2l1-pole throughout this work) interactions are included.31 Since the 1/r5 terms are approximated, the expansion is referred to as an approximate multipole expansion (AME). However, the recent implementation of the computationally efficient MPOLE module32 for multipole expansions with Ewald methods into the CHARMM molecular mechanics/ dynamics program33 uses an algorithm that truncates the interaction energy at order l of the multipole rather than power n of 1/r. Specifically, for a single site on the oxygen of a water molecule, an expansion up to the octupole is necessary and sufficient for an accurate representation of the electrostatic potential of a water molecule, but the interaction energy up to the octupole in the current MPOLE algorithm includes quadrupole−quadrupole and dipole−octupole interactions in the 1/r5 term, quadrupole−octupole interactions in the 1/r6 term, and octupole−octupole interactions in the 1/r7 term. However, the exact multipole expansion also includes the dipole−hexadecapole contribution in the 1/r6 term and the quadrupole− hexadecapole and dipole−25-pole contributions in the 1/r7 term. If the hexadecapole and 25-pole moments are large, this can lead to serious errors. Since MPOLE greatly enhances the utility of the molecular multipole approach for biomolecular simulations because of its efficiency, the accuracy of the truncation schemes should be explored. It is also important to emphasize that the differences between the truncation schemes become moot if the higher order multipole moments are small, in addition to when they are separated by large distances. For instance, when multipoles are

electrostatic potential of a water molecule, such as the SSDQO119 and SCME20 models. Another question in the accuracy of empirical potential energy functions for water is the method used to represent the charge distribution of a water molecule. The most common representations use nonpolarizable “partial charges” placed on the nuclei and sometimes on additional dummy sites.3 Models with sites only on the nuclei (i.e., three-site models) do not represent out-of-plane electron density and are deficient in reproducing many liquid properties. Four-site models have a dummy site in the molecular plane to increase the quadrupole, while five-site models have two dummy sites out of the molecular plane to represent “lone pair” character. Studies show that models with six sites are necessary to reproduce the quantum mechanical multipoles;7,21 however, the additional sites greatly decrease computational efficiency. More fundamentally, the use of dummy sites implies electron density is being represented by point charges. In addition, partial charges on even nuclear sites are not measurable quantities. Thus, partial charges must be obtained by fitting the electrostatic potential due to the electron density predicted by a quantum mechanical calculation on selected sites or by using schemes to partition that electron density. Other approaches for improvement over three-site models instead of (or in addition to) dummy sites include adding polarizability,22 higher multipoles,23,24 and nonpolarizable25 or polarizable26 Gaussian charge density onto the nuclear sites, all at an added cost over three-site models. A particular concern with adding polarizability to three- and four-site models is the lack of out-ofplane negative charge in the gas phase, which is seen in quantum mechanical calculations regardless of phase.7 Since computational efficiency is critical for many problems such as solvation of biological macromolecules, another approach has been to reduce the number of sites to a single site on the oxygen but to include orientational dependence of the electrostatics via a multipole expansion of the charge distribution of a water molecule. Thus, this approach gives models that are based on the electrostatic potential of a water molecule rather than the interaction energy between water molecules. Multipole expansions of the electrostatic potential are a formally exact description outside the charge distribution in the limit of infinite multipoles as well as the limit of infinite distances.27 In other words, multipole expansions are just another way of representing the electrostatics of a molecule. Although often thought of as less accurate than site models, multipoles actually form a complete orthogonal basis set. In addition, molecular multipoles can be measured in experiments up to the quadrupole in the gas phase or calculated directly from electronic structure calculations to any order. Finally, the multipole approach gains in computational efficiency over multisite models because only a single distance is needed between pairs of water molecules but can lose in computational efficiency as higher order multipoles are included. The loss in efficiency is because a multipole of order l, or 2l-pole, is represented in Cartesian space by a tensor of order l and length 3 (i.e., the dimensionality in Cartesian space) so the efficiency of the algorithms used for tensor operations becomes increasingly important for higher order multipoles. Thus, a careful balance is needed between accuracy obtained by including additional higher order multipoles and efficiency obtained by including the minimal number of multipoles. Other work indicates that molecular multipoles centered on the oxygen up to the octupole are necessary and sufficient to B

DOI: 10.1021/acs.jpcb.5b09565 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B placed on each nucleus, or “atomic multipoles”, as in the AMOEBA model, 34 they are representing the charge distribution of a single nucleus and surrounding electron density so the higher order multipoles may be relatively small at intermolecular distances, although recent studies of such a model indicate truncation at the 1/r5 terms rather than order l of the multipoles is also needed to model liquid properties.35 However, molecular multipoles are due to the charge distribution of multiple nuclei (i.e., the oxygen and two hydrogens in a water molecule) and associated electron density so that the charges are larger and more widely separated and the higher order moments are larger. The object here is to examine how well a molecular multipole expansion on a single site can reproduce the electrostatic potential due to the charge distribution of a small molecule with only one heavy atom. Specifically, the truncation of the multipole expansion is explored for molecular multipole moments centered on the water oxygen. Truncation at an order of multipole, i.e., either at the quadrupole− quadrupole term or at the octupole−octupole term, is compared to truncating at 1/r5, which includes both the quadrupole−quadrupole and dipole−octupole terms. As a test, comparisons are made between a site model for water and multipole models with the different truncation schemes that use multipoles calculated from the site model, similar to previous studies of the SSDQO expansion.31 The SPC/E model for water, which is a simple three-site model with an O−H bond length of 1 Å and a tetrahedral H−O−H bond angle,36 is used because the exact tetrahedral angle simplifies the analysis. The Lennard-Jones parameters from SPC/E are also used for the multipole models. First, the different truncated multipole expansions are compared to the site electrostatics in the interaction potential energies between two water molecules in certain orientations; special attention is paid to how well the truncated multipole expansions reproduce the potential energy near conformations known to be minima in the gas or liquid phase. Second, the truncated multipole expansions are compared to the site electrostatics for the structure and other properties of the liquid in MD simulations; this can be viewed as using a sampling procedure to assess how well the electrostatic interaction energy is reproduced for the most important contact conformations in the liquid state as well as the indirect effects of molecules beyond the first hydration shell.

■

Uij(r) =

1 [μ ·μ − 3(n ·μi )(μj ·n)] r3 i j

1 [−2μi ·Θj ·n + 2n ·Θi ·μj + 5(n ·μi )(Θj:n(2)) r4 1 − 5(n(2):Θi)(μj ·n)] + 5 [2Θi:Θj − 20(n ·Θi) ·(Θj ·n) 3r +

+ 35(n(2):Θi)(Θj:n(2)) + 9μi ·(Ωj:n(2)) + 9(n(2):Ωi) ·μj − 21(n ·μi )(Ωj ∴ n(3)) − 21(n(3) ∴ Ωi)(μj ·n)] +

1 [−2Θi:(Ωj ·n) + 2(n ·Ωi):Θj + 14(n ·Θi) ·(Ωj:n(2)) r6

− 14(n(2):Ωi) ·(Θj ·n) − 21(n(2):Θi)(Ωj ∴ n(3)) − 21(n(3) ∴ Ωi)(Θj:n(2)) + dipole−hexadecapole term] +

1 [2Ωi ∴ Ωj − 42(n ·Ωi):(Ωj ·n) 5r 7

+ 126(n(2):Ωi)· (Ωj:n(2)) − 231(n(3) ∴ Ωi)(Ωj ∴ n(3)) + quadrupole−hexadecapole term + dipole−25 decapole term] + terms of higher power of 1/ r

(1)

where r = rn is the internuclear vector from water a to b and μ, Θ, and Ω are the traceless dipole, quadrupole, and octupole moment tensors, respectively. All terms including multipoles up to the octupole are shown explicitly. Also, the dyadic products are denoted by [n(2)]ij = ninj and [n(3)]ijk = ninjnk, and the matrix contractions are denoted by A·B = ∑i AiBi, A:B = ∑ij AijBij, and A∴B = ∑ijk AijkBijk. From eq 1, it is apparent that the 1/r5 term contains contributions from both the quadrupole−quadrupole and dipole−octupole interactions. The geometry of a water molecule is defined by its OH bond length bOH and its HOH bond angle θHOH. To define the multipoles, a molecular coordinate system is located with the center at the oxygen, the z-axis along the θHOH bisector with the positive z-direction pointing toward the hydrogens, and the yz plane defined by the oxygen and hydrogens. Also, two sites L are located at a distance bOL in the negative z-direction in the xz plane with the bisector of the LOL angle θLOL along the z-axis. In this coordinate system, the elements of the μ, Θ, and Ω are μ = (0, 0, μ)

(2)

⎛ ⎞ 1 0 ⎟ ⎜−Θ2 − Θ0 0 2 ⎜ ⎟ ⎜ ⎟ 1 Θ= Θ2 − Θ0 0 ⎟ ⎜0 2 ⎜⎜ ⎟⎟ 0 Θ0 ⎠ ⎝0

(3)

⎛ ⎞ 1 0 ⎟ ⎜−Ω 2 − Ω 0 0 2 ⎜ ⎟ ⎟ 1 [Ω]ijz = ⎜ Ω2 − Ω0 0 ⎟ ⎜0 2 ⎜⎜ ⎟⎟ 0 Ω0 ⎠ ⎝0

(4)

The other elements of the Ω matrix can be obtained by symmetry. More details can be found in Niu et al.7 The lowest four multipole moments of the SPC/E model, which has partial charges qH on each H and −2qH on the O

THEORY

The multipole expansion for the electrostatic interaction energy between two water molecules can be written as C

DOI: 10.1021/acs.jpcb.5b09565 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

The NVT MD simulations were at ρ = 0.995 g/cm3, the density of SPC/E at 1 atm,3 and were maintained at T = 298 K using the Nosé−Hoover39,40 thermostat. All simulations employed the leapfrog Verlet algorithm41 with a time step of 1 fs and cubic periodic boundary conditions. The simulations were initially prepared using coordinates from a preequilibrated box of 512 water molecules at 0.995 g/cm3. Velocities were assigned according to a Gaussian distribution every 200 fs for 50 ps and then scaled every 1 ps if the temperature goes outside the range 298 ± 5 K for 250 ps. The system was then equilibrated without restraints for 1.2 ns; the trajectory for analysis was the subsequent 12 ns of simulation. The average properties were calculated using coordinates at 1 ps intervals. To estimate error bars, the standard deviations were calculated by dividing the 12 ns trajectories into three 4 ns segments. The diffusion constants were not corrected for box size dependence;42 they are estimated to be ∼0.25 × 10−5 g/ cm3 too small. Quantum mechanical calculations were performed using Gaussian 09.43 The calculations were at the second order Møller−Plesset (MP2) perturbation theory level using Dunning’s correlation consistent quadruple-ζ basis set with valence polarization and diffuse functions (aug-cc-pVQZ)44 as in previous work.6,7 For each conformation, three calculations were carried out: (a) a single point energy using the TIP3P45 molecular geometry; (b) a constrained geometry optimization with TIP3P internal coordinates fixed; and (c) a full geometry optimization.

(where qH = 0.4238 e), bOH = 1 Å and θHOH = 109.47°, are summarized in Table 1. For contrast, the multipoles of an “antiTable 1. Multipole Moments moment

SPC/E

AST

μ0 (D) Θ0 (D Å) Θ2 (D Å) Ω0 (D Å2) Ω2 (D Å2) H0 (D Å3) H2 (D Å3) H4 (D Å3)

2.350 0.000 2.035 −1.567 1.959 −1.583 1.131 0.989

2.350 0.000 2.035 −1.567 0.000 0.000 1.131 0.000

symmetric tetrahedral”, or AST, model, which has partial charges 1/2qH on each H and −1/2qH on two “lone pair” sites L, the same bOH and θHOH as SPC/E, bOL = 1 Å, and θLOL = 109.47°, are also given. While AST has the exact same dipole and quadrupole moments as SPC/E, the importance of the octupole can be seen because it is the lowest moment that differentiates the AST charge distribution from SPC/E. In addition, it is apparent that the hexadecapole moments are large for a water molecule.

■

METHODS MD simulations were performed using the molecular mechanics package CHARMM version c40a2.33 The multipole expansions up to the quadrupole−quadrupole or up to the octupole−octopole (i.e., truncated in order l of the multipole) are evaluated using the MPOLE module.32 The multipole expansion up to all 1/r5 terms is evaluated using MPOLE up to the quadrupole−quadrupole term and adding a separate dipole−octupole term. MPOLE uses the particle mesh Ewald (PME) method for the long-range interactions,37 with a Bspline of order 6, a k value of 0.55 Å−1, and a real space truncation at 7 Å; in addition, the long-range interactions of the explicitly added exact dipole−octupole term used the CHARMM switching function from 5 to 7 Å. For the simulations of the SPC/E site model, the Coulomb interactions between partial charges also used PME, with a B-spline of order 4, a k value of 0.4 Å−1, and a real space truncation at 10 Å. The long-range interactions of the Lennard-Jones potential for both the site and multipole models were handled using the CHARMM switching function from 8 to 10 Å. The cutoff for the pair and image lists was 12 Å, which was updated heuristically. The SHAKE38 algorithm was used to maintain rigid bodies.

■

RESULTS Different truncation schemes for multipole expansions for water centered on the oxygen are assessed by comparing results from truncated expansions using multipole moments calculated from SPC/E partial charges in the SPC/E geometry (Table 1) with results from the site SPC/E model. The truncation schemes for the expansion (eq 1) are after the quadrupole−quadrupole term, referred to as the Θ−Θ truncation, after the quadrupole− quadrupole and dipole−octupole terms, referred to as the μ−Ω truncation, after the quadrupole−octupole term, referred to as the Θ−Ω truncation, and after the octupole−octupole term, referred to as the Ω−Ω truncation, where only multipoles up to the octupole are considered. The truncations in multipole order, i.e., the Θ−Θ and Ω−Ω truncations, and in power of 1/r, i.e., the μ−Ω truncation, are focused on. In all simulations, the Lennard-Jones parameters for the 12−6 potential centered on the oxygen are σ = 3.166 Å and ε = 0.155 kcal/mol as in the SPC/E model. Four conformations of a pair of water molecules are focused on because they were found to be important in simulations of

Figure 1. Schematics of water dimer conformations with angle definitions (a) Hbcis, (b) Hbnlp, (c) Hbtrn, and (d) SbS. The center of the coordinate system is defined on the oxygen of the left-hand water molecule. D

DOI: 10.1021/acs.jpcb.5b09565 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B liquid water using TIP3P, SPC/E, TIP4P, and ST2.29 The conformations are defined with respect to a central water located with its dipole vector in the positive z-direction by first rotating the second water molecule around its dipole by ϕμμ, next rotating its dipole vector by θμμ with respect to the zdirection, and finally shifting it by a oxygen−oxygen distance rOO at an angle of θOO with respect to the z-direction (Figure 1). The first three conformations are hydrogen-bonded, in which the OH bond of the central molecule points to different parts of the second molecule. First, the central OH points to a “lone pair” of the second molecule oriented with its dipole also in the z-direction (Figure 1a), which is referred to as Hbcis, since the hydrogens of the two molecules are cis to each other. This is the most probable conformation in the liquid phase.29 Second, the OH points to the dipole vector of the second molecule (Figure 1b), which is referred to as Hnlp for “no lone pair”, since it should be the preferred direction for a dimer if there are no lone pair effects. This is the minimum energy conformation for an SPC/E dimer46 but is not the preferred conformation in the hydrogen-bonded tetrahedral network of the liquid. Third, the OH points to the other “lone pair” (Figure 1c), which is referred to as Hbtrn, since the hydrogens of the two molecules are trans to each other. This is the minimum energy conformation of a water dimer according to quantum mechanical calculations.6 The last is a side-by-side conformation (Figure 1d), which is referred to as SbS. This is a non-hydrogen-bonded conformation found in the liquid for SPC/E and TIP4P but much less so for TIP3P and especially ST2.29 Interaction Energies between a Pair of Water Molecules. The interaction energies between two water molecules were calculated to see how well the truncation schemes reproduced the energies of SPC/E. First, the dependences of the interaction energy on θOO for the different truncation schemes are examined (Figure 2), especially the behavior near the minima. Three different orientations of the water molecule, in which the orientation of the two molecules with respect to each other are defined by ϕμμ and θμμ of the Hbcis (Figure 2a), Hbtrn (Figure 2b), and SbS (Figure 2c) conformations, are shown. The plots show that the Θ−Θ truncation is not able to reproduce the site model. In addition, careful examination shows the μ−Ω truncation is preferable to the Ω−Ω truncation. The μ−Ω truncation reproduces the energy dependence of SPC/E for the Hbcis orientation, the most common conformation in the liquid, and the minimum energy of the Hbtrn orientation, although it somewhat overestimates the exact value of θOO where this minimum occurs. In contrast, the Ω−Ω truncation predicts minima for the hydrogen-bonded conformations that are ∼2 kcal/mol too high in energy compared to SPC/E. For the SbS orientation, the μ−Ω truncation underestimates the depth of the SPC/E minimum by ∼2 kcal/mol, while the Ω−Ω truncation overestimates it by ∼2 kcal/mol, which is examined further below. The dependences of the interaction energy on θμμ when ϕμμ = 90° and θOO = 1/2θtet, as in the three hydrogen-bonded configurations of Figure 1, for the different types of truncation are also shown (Figure 3). Again, it is clear that adding the octupole contributions is an improvement over the Θ−Θ truncation. Interestingly, the μ−Ω truncation gives good agreement with SPC/E for θμμ < ∼100°, which includes the Hbcis and Hbnlp conformations, while the Ω−Ω truncation is good for θμμ > ∼130°. However, since the conformations with

Figure 2. Electrostatic potential in kcal/mol as a function of θOO at r = 2.8 Å for (a) ϕμμ = 90°, θμμ = 0° with Hbcis (∗), (b) ϕμμ = 90°, θμμ = θtet with Hbtrn (●), and (c) ϕμμ = 0°, θμμ = 180° with SbS (×). SPC/E (black) and multipole expansion up to quadrupole−quadrupole terms (blue), up to dipole−octupole terms (green), up to quadrupole− octupole terms (orange), up to octupole−octupole terms (red), and up to octupole−octupole terms but with Ω2 at half of the SPC/E value (magenta dashed).

Figure 3. Electrostatic potential in kcal/mol as a function of θμμ at r = 2.8 Å for (a) ϕμμ = 90°, θOO = 1/2θtet with Hbcis (∗), Hbnlp (▲), and Hbtrn (●) indicated. SPC/E (black) and multipole expansion up to quadrupole−quadrupole terms (blue), up to dipole−octupole terms (green), up to quadrupole−octupole terms (orange), up to octupole− octupole terms (red), and up to octupole−octupole terms but with Ω2 at half of the SPC/E value (magenta dashed).

E

DOI: 10.1021/acs.jpcb.5b09565 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B Table 3. Water Dimer Conformationsa

θμμ > ∼130° are highest in energy, the conformations whose energy are well represented by the Ω−Ω truncation would be least sampled in a simulation. In addition, the relative stabilities of the four conformations for the different truncation schemes are examined (Table 2). Table 2. Water Dimer Energies (kcal/mol)a calculation

Hbcis

Hbnlp

Hbtrn

SbS

SPC/E Θ−Θ μ−Ω Ω−Ω Ω−Ω, with 1/2Ω2 of SPC/E SSDQO single point energy constrained geometry optimization geometry optimization

−6.04 −5.00 −5.95 −5.04 −5.82 −6.10 −3.98 N/A N/A

−8.24 −6.71 −8.26 −6.39 −5.93 −7.19 −4.51 N/A N/A

−8.17 −7.23 −7.40 −6.94 −7.79 −7.51 −4.85 −5.04 −5.09

−6.68 −3.16 −5.58 −8.54 −6.33 −1.73 −3.93 −4.01 −4.04

calculation

conformation

rOO (Å)

θOO (deg)

ϕμμ (deg)

single point energy and potential energy functions

Hbcis

2.80

1

/2θtet

90

0

Hbnlp Hbtrn SbS Hbtrnb

2.80 2.80 2.80 2.91

1

/2θtet /2θtet 100 58

90 90 0 90

/2θtet θtet 180 113

SbS Hbtrnb SbS

2.76 2.90 2.75

103 58 103

0 90 0

180 114 180

constrained geometry optimization geometry optimization

1

θμμ (deg)

1

a

See text for definitions. bBoth the NLP and Hbcis conformations end up in the Hbtrn conformation upon geometry optimization.

Structure of Liquid Water. To compare the interaction energies in a larger region of coordinate space, the different truncation schemes were also compared to SPC/E in MD simulations. The simulations weight the probability of conformations in coordinate space so a good truncation scheme should give good agreement with the structure and properties of SPC/E. First, radial distribution functions, or g(r), for SPC/E and different truncation schemes were examined (Figure 4). For the

a

For site and multipole models, the electrostatic contribution at 2.8 Å is given; the Lennard-Jones interaction at this distance is 1.412 kcal/ mol.

The site model SPC/E predicts Hbnlp is the lowest energy conformation and that SbS is less stable than Hbtrn by ∼1.5 kcal/mol but is actually more stable than Hbcis by ∼0.6 kcal/ mol. The relative stabilities for the μ−Ω truncation are more like SPC/E than the Ω−Ω truncation, consistent with the above discussion. However, a more drastic error becomes apparent in the Ω−Ω truncation because it predicts that the SbS conformation is much more stable than the Hbcis and Hbtrn conformations by 3.5 and 1.5 kcal/mol, respectively, while the μ−Ω truncation predicts the SbS conformation is slightly less stable than the Hbcis conformation by ∼0.4 kcal and even less stable relative to the Hbtrn conformation by ∼2 kcal/mol. Not only are the deviations of the μ−Ω truncation from SPC/E for the relative energy of the SbS conformation smaller than those from the Ω−Ω truncation, but they could be considered even less severe if they would lead to a better representation of real water−water interactions than site SPC/E. To investigate the relative stability of the SbS conformation for realistic systems, quantum mechanical calculations of water dimers in the gas phase were carried out. Although these are in the gas phase and SPC/E accounts for the average polarization in the liquid phase by somewhat larger partial charges, these calculations give an idea of what the relative stability should be. Single point energies were carried out for the Hbcis, Hbnlp, Hbtrn, and SbS conformations (Figure 1). The single point energies indicate that Hbtrn is more stable than SbS by ∼1 kcal/ mol, and the geometry optimizations show that these two conformations are not far from the real minima (Table 3) with similar energy differences (Table 2). The single point energies also indicate that the SbS and Hbcis conformations are approximately equally stable, although the geometry optimizations show that the Hbcis and Hbnlp conformations are not significant local minima because optimizations beginning from these conformations go to the Hbtrn conformation. Thus, the Ω−Ω truncation appears to give a completely incorrect picture of the interaction energy surface because it predicts that the SbS conformation is more stable than any hydrogen bonded conformation by at least ∼1.5 kcal/mol while the μ−Ω truncation appears to be the best for describing the SPC/E interaction energy.

Figure 4. Radial distribution function for SPC/E (black), up to quadrupole−quadrupole terms (blue), up to dipole−octupole terms (green), up to the octupole−octupole (red), and up to octupole− octupole terms but with Ω2 at half of the SPC/E value (magenta dashed). From top to bottom, the OO, HH, and OH are shown.

simulation using the Θ−Θ truncation, gOH(r) indicates too little hydrogen bonding and gOO(r) indicates the liquid is too unstructured compared to SPC/E. For the simulation using the μ−Ω truncation, gOH(r) indicates much better hydrogen bonding, gHH(r) indicates the molecules are reasonably oriented with respect to each other, and gOO(r) indicates reasonable tetrahedral structure, although it is still slightly understructured. However, for the simulation using the Ω−Ω truncation, gOH(r) actually indicates worse hydrogen bonding than the μ−Ω truncation, and gHH(r) indicates that the peak in gOH(r) is actually due to SbS conformations, which have similar OH distances as the hydrogen-bonded conformations, although the donor OH is not in line with the OO internuclear vector and thus not hydrogen-bonded. This is consistent with the results for the pair interaction energies in the previous section, in which the Ω−Ω truncation predicts the global minimum is the F

DOI: 10.1021/acs.jpcb.5b09565 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

around the maxima when the multipole interaction energy is used in MD simulations, since the maxima are rarely sampled. Here, molecular multipole expansions are used to describe the electrostatic interactions even at contact between water molecules, so that higher order multipoles and higher powers of 1/r are needed. The results here further demonstrate the importance of both the quadrupole and the octupole for a molecular multipole description of water and, more importantly, show the multipole expansion for the interaction energy should be truncated according to powers of 1/r. The exact μ−Ω truncation gives a better description of the potential energy of SPC/E near the most probable conformations found in the liquid than the Ω−Ω truncation and is the most efficient way of including octupole effects without including hexadecapole effects. Of course, this does not address atomic multipoles placed at each nuclear center such as in the AMOEBA model and care must be taken in making generalizations from one type of multipole expansion to the other because the charge distributions being represented are very different. Considering a balance of accuracy and CPU speed, a truncation of the multipole expansion at 1/r5 for the electrostatic interaction energy between water molecules works well for water, which is approximately as efficient as three-site models when MPOLE is used. Small hexadecapole effects are apparent but can be corrected for by slight changes. For instance, the van der Waals parameters can be modified slightly to account for the differences in potential energy of the lowest energy conformation in the liquid. However, care must be made in arbitrary adjustment of the multipole moments, since this can result in very different liquid structure not apparent in g(r) but more or less apparent in other properties, as illustrated here by reducing the cubic octupole in a Ω−Ω truncation and many other tests not presented here. In addition, it appears that truncation using the approximate 1/r5 term in the AME of SSDQO is more accurate in representing SPC/E31 than truncation at the exact 1/r5 term presumably because it accounts for some cancelation by higher multipoles. However, only the dipole−dipole interaction in SSDQO is treated using Ewald summations so far. The stability of the different conformations using different truncation schemes can be understood by examining the higher power of 1/r contributions to the interaction energy, which govern the short-range interactions. For instance, the Hbcis conformation is stabilized by the favorable interaction of the cubic quadrupoles (Figure 5a) and between the dipoles and linear octupoles (Figure 5b,c), while the cubic octupoles interact unfavorably, so that including the Ω−Ω term without the other 1/r6 and 1/r7 terms involving higher multipoles in the exact expansion leads to destabilization of this conformation. In contrast, the SbS conformation is stabilized by the favorable

SbS conformation rather than a hydrogen-bonded conformation. Thus, gOO(r) for this truncation actually has a minimum at 4.5 Å where the tetrahedral maximum should be and an unphysical maximum at 5.5 Å. In sum, the μ−Ω truncation shows the most accurate liquid structure of the three schemes. Next, the heat of vaporization ΔHvap, diffusion coefficient D, and dielectric permittivity ε were calculated (Table 4) along Table 4. Properties of Liquids Using Different Truncation Schemes versus the SPC/E Site Modela ΔHvap (kcal/mol) Exp SPC/E site Θ−Θ μ−Ω Ω−Ω Ω−Ω, 1/2Ω2 of SPC/E a b

b

10.51 11.71 ± 0.01 7.69 ± 0.01 9.87 ± 0.01 8.91 ± 0.01 8.81 ± 0.01

D (10−5 cm2/s) 2.30 2.43 7.87 2.73 5.14 4.17

c

± ± ± ± ±

0.04 0.10 0.02 0.03 0.03

ε 78 71 63 63 54 56

d

± ± ± ± ±

2 1 3 2 1

The NVT simulations were at ρ = 0.995 g/cm3 and T = 298 K. Reference 47. cReference 48. dReference 49.

with experimental values;47−49 D and ε probe different parts of water structure.12 For the Θ−Θ and Ω−Ω truncations, ΔHvap is much too small and D is much too large, which indicates that the hydrogen bonding to the nearest neighbors forming the tetrahedral hydration shell is too weak. Also, ε is too low for the Ω−Ω truncation due to the overestimation of SbS conformations, which lead to the largest possible decrease in the Kirkwood g-factor.50 Conversely, the μ−Ω truncation is able to mimic the liquid state properties of SPC/E. Finally, the Ω−Ω truncation but with Ω2 reduced to half of the SPC/E value was also examined to see if the Ω−Ω truncation could be improved by altering the moments. However, arbitrarily adjusting the multipole moments in this truncation scheme is not a good option because it may lead to a very different liquid structure. In particular, although this scheme appears to lead to better agreement with the g(r) for SPC/E than the μ−Ω truncation (Figure 4), the resulting liquid has a very different structure than SPC/E, since the excessively large D indicates that the hydrogen bonding is too weak and the low ε indicates that the amount of Hbcis is too low (Table 4).

■

DISCUSSION The multipole expansion is simply an expansion of the electrostatic potential in an orthogonal basis set, but the expansion weights both maxima and minima in the interaction energy between charge distributions equally. However, accurate descriptions around the minima are more important than

Figure 5. Charge distributions of moments (a) planar quadrupole−planar quadrupole, (b) linear octupole−dipole, (c) dipole−linear octupole, and (d) cubic octupole−cubic octupole. G

DOI: 10.1021/acs.jpcb.5b09565 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

molecular dynamics simulations. For simulations of a condensed phase, it is important to include enough powers of 1/r for the short-range interactions, and therefore, the importance of higher order multipole moments in describing the charge distribution must be assessed. The work presented here uses a site model for water, namely, SPC/E, for assessment of how well truncated multipole expansions reproduce a known electrostatic potential. The results here show that multipole moments up to the octupole are necessary for polar molecules with tetrahedral symmetry like water, since truncating the expansion at the quadrupole−quadrupole term leads to large deviations with SPC/E in the electrostatic potential between water molecules and the properties of liquid water in simulation. More importantly, the results show that truncating the expansion at 1/r5 is a better way of including octupole effects than truncating at the octupole−octupole term without higher order multipole moments, since truncation at 1/r5 leads to results that are closest to SPC/E in the electrostatic potential between water molecules and the properties of liquid water in simulation. In particular, truncation at 1/r5 is important, since the octupole−octupole (specifically, the Ω2−Ω2) term is large for water but canceled by terms involving multipoles of higher order than the octupole. However, the disagreement of the thermodynamic, dynamic, and dielectric properties of the SPC/ E model with experimental values demonstrates that SPC/E is lacking in describing real water. Thus, future research will use multipoles obtained from quantum mechanical electronic structure calculations.

interaction of the quadrupole and cubic octupole (not shown) and the cubic octupoles (Figure 5d), so that including the Ω−Ω terms without the other relevant terms involving higher multipoles leads to overstabilization of the SbS conformation. These studies also provide insights into what the actual potential energy between water molecules looks like and how the molecules are oriented with respect to each other in real liquid water. From the analysis of orientations of waters in the first hydration shell of a central water for SPC/E, TIP3P, TIP4P, and ST2 by Mason and Brady,29 a large population of Hbcis is common to all four models. However, they differ in populations of the other conformations in Figure 1, which will affect different properties. For instance, SPC/E and TIP4P, which have large cubic octupoles like all three- and four-site models, give especially large populations of SbS in the fifth nearest neighbor in the liquid state.29 Since the Ω−Ω terms favor the SbS conformation in the liquid, which tends to reduce the dielectric permittivity because the dipoles of the two water molecules are antiparallel, this indicates that the reason SPC/E and TIP4P predict too low dielectric permittivities is because of too much SbS conformation due to too high cubic octupoles. Conversely, TIP3P and five-site models appear to have high dielectric permittivity because of large populations of Hbcis, which may be overestimated because of the small quadrupole in both TIP3P and five-site models, and at least in part to small cubic octupoles in the five-site models. In addition, the “planar” models SPC/E, TIP3P, and TIP4P all have significant populations of Hbnlp in the four nearest neighbors.29 The electrostatic potentials for “site” SPC/E here and other planar models in other works46 have minima at Hbnlp, consistent with the large population found in simulations of the liquid. However, the QM studies here of the water dimer indicate a flat potential between the three hydrogen-bonded conformations, with a minimum at Hbtrn. In the liquid, Hbtrn apparently is no longer the most favored conformation because the hydrogen-bonded network formed in liquid water favors the Hb cis conformation, but a range of hydrogen-bonded conformations with different θμμ is favored over the SbS conformation. Finally, the disagreement of the thermodynamic, dynamic, and dielectric properties of the SPC/E model with experimental values points out that better multipole moments need to be utilized. Ultimately, rather than using the SPC/E multipoles, multipole moments for water in the liquid phase will be determined. While determining the multipoles in the gas phase is straightforward, it is less straightforward in the liquid phase, both because quantum mechanical simulations of water in the liquid phase are still of limited accuracy and because the best method for partitioning the electron density between molecules in a liquid is not clear.3 Thus, optimization of quantum mechanical moments may be necessary, but the literature values provide a reasonable starting point and any large deviations from these values must be accounted for. The 1/r5 truncation scheme will be made available in the CHARMM developmental code, and optimized multipole moments will be made available via the CHARMMing51 web site.

■

AUTHOR INFORMATION

Corresponding Author

*Phone: 202-687-3724. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors are grateful for support from the National Science Foundation through Grant No. CHE-1464766 and from the McGowan Foundation. This research was also supported in part by the Intramural Research Program of the NIH, National Heart, Lung, and Blood Institute. This work used computer time on the Extreme Science and Engineering Discovery Environment (XSEDE) granted via MCB990010, which is supported by National Science Foundation Grant No. OCI1053575; the Medusa cluster, which is 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.

■

REFERENCES

(1) Ball, P. Water  an Enduring Mystery. Nature 2008, 452, 291− 292. (2) Poole, P. H.; Sciortino, F.; Grande, T.; Stanley, H. E.; Angell, C. A. Effect of Hydrogen Bonds on the Thermodynamic Behavior of Liquid Water. Phys. Rev. Lett. 1994, 73, 1632−1635. (3) Ichiye, T. Water in the Liquid State: A Computational Viewpoint. Adv. Chem. Phys. 2014, 155, 161−200. (4) Verhoeven, J.; Dymanus, A. Magnetic Properties and Molecular Quadrupole Tensor of the Water Molecule by Beam-Maser Zeeman Spectroscopy. J. Chem. Phys. 1970, 52, 3222−3233.

■

CONCLUSIONS The results here demonstrate that truncated molecular multipoles can be used to give an atomistic representation of the electrostatic potential due to a molecule with a complex charge distribution, namely, water. Thus, molecular multipole expansions can be used for potential energy functions for H

DOI: 10.1021/acs.jpcb.5b09565 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (5) Clough, S. A.; Beers, Y.; Klein, G. P.; Rothman, L. S. Dipole Moment of Water from Stark Measurements of H2O, HDO, and D2O. J. Chem. Phys. 1973, 59, 2254−2259. (6) Xantheas, S. S.; Dunning, T. H., Jr. Ab Initio Studies of Cyclic Water Clusters (H2O)N, N = 1−6. I. Optimal Structures and Vibrational Spectra. J. Chem. Phys. 1993, 99, 8774−8782. (7) Niu, S.; Tan, M.-L.; Ichiye, T. The Large Quadrupole of Water Molecules. J. Chem. Phys. 2011, 134, 134501. (8) Coutinho, K.; Guedes, R. C.; Costa Cabral, B. J.; Canuto, S. Electronic Polarization of Liquid Water: Converged Monte-CarloQuantum Mechanics Results for the Multipole Moments. Chem. Phys. Lett. 2003, 369, 345−353. (9) Osted, A.; Kongsted, J.; Mikkelson, K. V.; Åstrand, P.-O.; Christiansen, O. Statistical Mechanically Averaged Molecular Properties of Liquid Water Calculated Using the Combined Coupled Cluster/Molecular Dynamics Method. J. Chem. Phys. 2006, 124, 124503−16. (10) Silvestrelli, P. L.; Parrinello, M. Structural, Electronic, and Bonding Properties of Liquid Water from First Principles. J. Chem. Phys. 1999, 111, 3572−3580. (11) Delle Site, L.; Alavi, A.; Lynden-Bell, R. M. The Electrostatic Properties of Water Molecules in Condensed Phases: An Ab Initio Study. Mol. Phys. 1999, 96, 1683−1693. (12) Tan, M.-L.; Cendagorta, J. R.; Ichiye, T. The Molecular Charge Distribution, the Hydration Shell, and the Unique Properties of Liquid Water. J. Chem. Phys. 2014, 141, 244504. (13) Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. Development of an Improved Four-Site Water Model for Biomolecular Simulations: TIP4P-Ew. J. Chem. Phys. 2004, 120, 9665−9678. (14) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505-1−12. (15) Mahoney, M. W.; Jorgensen, W. L. A Five-Site Model for Liquid Water and the Reproduction of the Density Anomaly by Rigid, Nonpolarizable Potential Functions. J. Chem. Phys. 2000, 112, 8910− 8922. (16) Fanourgakis, G. S.; Xantheas, S. S. Development of Transferable Interaction Potentials for Water. V. Extension of the Flexible, Polarizable, Thole-Type Model Potential (TTM3-F, V. 3.0) to Describe the Vibrational Spectra of Water Clusters and Liquid Water. J. Chem. Phys. 2008, 128, 074506. (17) Kumar, R.; Wang, F.-F.; Jenness, G. R.; Jordan, K. D. A Second Generation Distributed Point Polarizable Water Model. J. Chem. Phys. 2010, 132, 014309. (18) Medders, G. R.; Babin, V.; Paesani, F. Development of a “FirstPrinciples” Water Potential with Flexible Monomers. Iii. Liquid Phase Properties. J. Chem. Theory Comput. 2014, 10, 2906−2910. (19) Te, J. A.; Ichiye, T. Temperature and Pressure Dependence of the Optimized Soft Sticky Dipole-Quadrupole-Octupole Water Model. J. Chem. Phys. 2010, 132, 114511. (20) Wikfeldt, K. T.; Batista, E. R.; Vila, F. D.; Jonsson, H. A Transferable H2O Interaction Potential Based on a Single Center Multipole Expansion: Scme. Phys. Chem. Chem. Phys. 2013, 15, 16542−16446. (21) Yu, W.; Lopes, P. E. M.; Roux, B.; MacKerell, A. D., Jr. Six-Site Polarizable Model of Water Based on the Classical Drude Oscillator. J. Chem. Phys. 2013, 138, 034508. (22) Rick, S. W.; Stuart, S. J. Potentials and Algorithms for Incorporating Polarizability in Computer Simulations. Rev. Comput. Chem. 2002, 18, 89−146. (23) Ren, P.; Ponder, J. W. Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation. J. Phys. Chem. B 2003, 107, 5933−5947. (24) Dykstra, C. E. Electrostatic Interaction Potentials in Molecular Force Fields. Chem. Rev. 1993, 93, 2339−2353. (25) Guillot, B. A Reappraisal of What We Have Learnt During Three Decades of Computer Simulations of Water. J. Mol. Liq. 2002, 101, 219−260.

(26) Jones, A.; Cipcigan, F.; Sokan, V. P.; Crain, J.; Martyna, G. J. Electronically Coarse-Grained Model for Water. Phys. Rev. Lett. 2013, 110, 227801. (27) Jackson, J. D. Classical Electrodynamics, 3rd ed.; Wiley: New York, 1975. (28) Kusalik, P. G.; Patey, G. N. On the Molecular Theory of Aqueous Electrolyte Solutions. Ii. Structural and Thermodynamic Properties of Different Models at Infinite Dilution. J. Chem. Phys. 1988, 89, 5843−5851. (29) Mason, P. E.; Brady, J. W. Tetrahedrality” and the Relationship between Collective Structure and Radial Distribution Functions in Liquid Water. J. Phys. Chem. B 2007, 111, 5669−5679. (30) Stone, A. J. The Theory of Intermolecular Forces; Clarendon: Oxford, U.K., 1996. (31) Ichiye, T.; Tan, M.-L. Soft Sticky Dipole-Quadrupole-Octupole Potential Energy Function for Liquid Model: An Approximate Moment Expansion. J. Chem. Phys. 2006, 124, 134504. (32) Simmonett, A. C.; Pickard, F. C., IV; Schaefer, H. F., III; Brooks, B. R. An Efficient Algorithm for Multipole Energies and Derivatives Based on Spherical Harmonics and Extensions to Particle Mesh Ewald. J. Chem. Phys. 2014, 140, 184101. (33) Brooks, B. R.; Brooks, C. L., III; MacKerell, A. D., Jr.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al. CHARMM: The Biomolecular Simulation Program. J. Comput. Chem. 2009, 30, 1545−1614. (34) Ren, P.; Ponder, J. W. Temperature and Pressure Dependence of the AMOEBA Water Model. J. Phys. Chem. B 2004, 108, 13427− 13437. (35) Shaik, M. S.; Liem, S. Y.; Popelier, P. L. A. Properties of Liquid Water from a Systematic Refinement of a High-Rank Multipolar Electrostatic Potential. J. Chem. Phys. 2010, 132, 174504. (36) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction Models of Water in Relation to Protein Hydration. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, Holland, 1981; pp 331−341. (37) Darden, T. A.; York, D. M.; Pedersen, L. G. Particle Mesh Ewald: An Nlog(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (38) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equation of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (39) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (40) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (41) Verlet, L. Computer ″Experiments″ on Classical Fluids. II. Equilibrium Correlation Functions. Phys. Rev. 1968, 165, 201−214. (42) Yeh, I.-C.; Hummer, G. System-Size Dependence of Diffusion Coefficients and Viscosities from Molecular Dynamics Simulations with Periodic Boundary Conditions. J. Phys. Chem. B 2004, 108, 15873−15879. (43) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, revision A.1; Gaussian, Inc.: Wallingford, CT, 2009. (44) Dunning, T. H., Jr. Gaussian Basis Sets for Use in Correlated Molecular Calculations. I. The Atoms Boron through Neon and Hydrogen. J. Chem. Phys. 1989, 90, 1007−1023. (45) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (46) Kiss, P. T.; Baranyi, A. Clusters of Classical Water Models. J. Chem. Phys. 2009, 131, 204310. (47) Wagner, W.; Pruss, A. International Equations for the Saturation Properties of Ordinary Water Substance. Revised According to the International Temperature Scale of 1990. Addendum to J. Phys. Chem. Ref. Data 16, 893 (1987). J. Phys. Chem. Ref. Data 1993, 22, 783−787. I

DOI: 10.1021/acs.jpcb.5b09565 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B (48) Mills, R. Self-Diffusion in Normal and Heavy Water in the Range 1−45°. J. Phys. Chem. 1973, 77, 685−688. (49) Lide, D. R. CRC Handbook of Chemistry and Physics; CRC Press: Boca Raton, FL, 1991; Vol. 72. (50) Kirkwood, J. G. The Dielectric Polarization of Polar Liquids. J. Chem. Phys. 1939, 7, 911−919. (51) Miller, B. T.; Singh, R. P.; Klauda, J. B.; Hodoscek, M.; Brooks, B. R.; Woodcock, H. L. CHARMMing: A New, Flexible Web Portal for CHARMM. J. Chem. Inf. Model. 2008, 48, 1920−1929.

J

DOI: 10.1021/acs.jpcb.5b09565 J. Phys. Chem. B XXXX, XXX, XXX−XXX