The vapor-liquid interface potential of (multi)polar fluids and its influence on ion solvation
TThe vapor-liquid interface potential of (multi)polar fluidsand its influence on ion solvation
Lorand Horv´ath
Universit´e de Toulouse; UPS; Laboratoire de Physique Th´eorique (IRSAMC); F-31062 Toulouse,France and CNRS; LPT (IRSAMC); F-31062 Toulouse, France andFaculty of Physics, University “Babe¸s-Bolyai”, 400084 Cluj-Napoca, Romania
Titus Beu
Faculty of Physics, University “Babe¸s-Bolyai”, 400084 Cluj-Napoca, Romania
Manoel Manghi
Universit´e de Toulouse; UPS; Laboratoire de Physique Th´eorique (IRSAMC); F-31062 Toulouse,France and CNRS; LPT (IRSAMC); F-31062 Toulouse, France
John Palmeri
Universit´e de Toulouse; UPS; Laboratoire de Physique Th´eorique (IRSAMC); F-31062 Toulouse,France and CNRS; LPT (IRSAMC); F-31062 Toulouse, France andLaboratoire Charles Coulomb UMR 5221 CNRS-UM2,D´epartement Physique Th´eorique, Universit´e Montpellier 2,Place Eug`ene Bataillon - CC070 F-34095 Montpellier Cedex 5 France ∗ (Dated: November 4, 2018)The interface between the vapor and liquid phase of quadrupolar-dipolar fluids is the seat of anelectric interfacial potential whose influence on ion solvation and distribution is not yet fully under-stood. To obtain further microscopic insight into water specificity we first present extensive classicalmolecular dynamics simulations of a series of model liquids with variable molecular quadrupole mo-ments that interpolates between SPC/E water and a purely dipolar liquid. We then pinpoint theessential role played by the competing multipolar contributions to the vapor-liquid and the solute-liquid interface potentials in determining an important ion-specific direct electrostatic contributionto the ionic solvation free energy for SPC/E water—dominated by the quadrupolar and dipolarcontributions—beyond the dominant polarization one. Our results show that the influence ofthe vapor-liquid interfacial potential on ion solvation is strongly reduced due to the strong partialcancellation brought about by the competing solute-liquid interface potential. I. INTRODUCTION
At the macroscopic interface between a liquid ( l ) andits vapor ( v ) phase there is a spatial inhomogeneity thatinduces a charge imbalance, producing an electric fieldand consequently a potential difference across the inter-face, φ lv = φ l − φ v . Despite extensive molecular sim-ulation studies at both the classical and quantum me-chanical levels over the past few decades [1–10, 12–15],a complete understanding of this potential, how it de-pends on the characteristics of the fluid studied, and itsrole in the solvation of ions is not yet at hand [16, 17].Such an understanding has become a pressing matter,because there is currently much interest in constructingmesoscopic models of electrolytes near vapor-liquid inter-faces and solid membrane surfaces and in nanopores [16–25, 27]. It is also now clear that the value of the interfacepotential observed experimentally depends on the probeused: although electron diffraction and holography tech- ∗ Current address; Electronic address: [email protected] niques may measure the full interface potential, electro-chemical techniques involved in ion solvation seeminglydo not [13–15, 17].Until now mesoscopic approaches to ion distributionhave either completely neglected the contribution of theinterfacial potential [20, 23–25] or, as already demon-strated in [21] and discussed in detail here, severely over-estimated its importance by treating finite size solutes aspoint test charges [22]. For finite size ions a second micro-scopic solute-liquid interface potential, φ ls = φ l − φ s , ex-ists, defined as the potential difference between the bulkliquid ( l ) and the center of a neutral solute ( s ) (fig. 1)[whose size is determined in classical Molecular Dynam-ics (MD) simulations by the short range repulsion of theLennard-Jones potential]. This second solute-liquid con-tribution is missed if the ions are approximated as pointtest charges. The potentially important role played bythis microscopic potential in determining ion distributionnear inhomogeneities needs to be clarified [10, 11] in orderto provide deeper theoretical insight into both molecularsimulations and experimental results. Furthermore, thecoupling between the quadrupolar and dipolar contribu-tions to the interface potential and their respective rolesin governing ion distribution need to be reconsidered. To a r X i v : . [ c ond - m a t . s o f t ] M a r do so we present extensive classical molecular dynam-ics simulations of model liquids that interpolate betweenSPC/E water (a classical three site partial charge model[28]) and a purely dipolar liquid.In physical terms our study can be viewed as part ofthe quest, still far from complete, for the physical compo-nents of the position dependent ionic Potential of MeanForce (PMF), Φ( r ), near dielectric interfaces and sur-faces arising from solvent-ion and ion-ion interactions af-ter the solvent degrees of freedom have been integratedout [16, 20–25]. We focus uniquely on the poorly under-stood role played by the interfacial potential in determin-ing the electrostatic contribution to ion solvation. In-deed, the extremely large discrepancy between the dilutelimit ionic PMF obtained from MD simulations and thosepredicted using an approximate mesoscopic approach in-corporating the contribution of the interfacial potentialin the point ion approximation led the authors of [21]to completely abandon their dilute limit mesoscopic ap-proach; they opted rather for extracting the dilute limitionic PMF directly from MD simulations and then in-jecting it into a generalized Poisson-Boltzmann equationto study salt concentration effects. In more recent workconcerning the optimization of the MD parameters of anon-polarizable model by fitting to experiment, the sameauthors attempted to get around the ambiguities plagu-ing the contribution of the vapor-liquid interfacial po-tential in determining the electrostatic contribution toion solvation by fitting only quantities independent ofthis contribution [26]. Although other contributions tothe ionic PMF and solvation-free energy, such as the hy-drophobic and dispersion ones, may play non-negligibleroles and therefore be important for interpreting molec-ular simulations and understanding experiments, thesecontributions will not be considered here (as they are al-ready fairly well understood thanks to recent progress inthis area) [20–26].One major impediment to obtaining the physicallyidentifiable mesoscopic contributions to the ionic PMF,Φ, arises from questions concerning the amplitude andsign of the vapor-liquid interface potential and the roleit plays in determining the ionic solvation free energy. Inorder to address these questions we compare a direct eval-uation from MD simulations of the two relevant electro-static contributions to the ionic solvation free energy forSPC/E water—a direct interfacial one that does not ac-count for solvent polarization due to the ionic charge anda polarization one that does—with simplified approachespreviously adopted in the literature (namely, a direct oneapproximating ions as point test charges and a simpleBorn-type polarization approximation, defined below). II. VAPOR-LIQUID INTERFACE POTENTIALAND IONIC PMF: STATE-OF-THE-ART
Near a planar vapor-liquid interface the local ion con-centration can be expressed in terms of the PMF, Φ( z ), as ρ i ( z ) = ρ il exp[ − Φ( z ) /k B T ], where z is the normal co-ordinate and ρ il is the ionic concentration in the bulk liq-uid (where Φ is taken to vanish). The total ion solvationfree energy can then be expressed as ∆ G ion = − Φ( z v ),where z v is in the vapor phase (see fig. 1). In theo-retical studies of both vapor-liquid water interfaces andmembrane-liquid surfaces, it has sometimes been hy-pothesized [21–24] that the bare interface potential en-ter the PMF via a simple direct electrostatic contribu-tion Φ (cid:48) pot ( z ) = q [ δφ ( z ) − φ lv ], where q is the ion charge, δφ ( z ) = φ ( z ) − φ v is the local value of the potential differ-ence and φ lv = δφ ( z l ) ( z l is at the center of the liquid slabfar from the interface). This approach, which amounts totreating a finite size ion as a point test charge q [21–24],is critically examined here for SPC/E water.Classical molecular dynamics (MD) simulations pre-dict potentials on the order of − . (cid:48) pot , would seem-ingly yield the dominant contribution ( ∼ k B T formonovalent ions) to the PMF over a substantial part ofthe interfacial region [21]. This approximation, however,predicts incorrect results for the PMF and correspond-ing ion density, when compared with MD simulations,both in the infinitely dilute limit (as already pointed outin [21]) and when incorporated into a modified Poisson-Boltzmann approach (to study higher electrolyte concen-trations [22]): neither the strong build-up of anions neara strongly hydrophobic uncharged surface ([22], Fig. 4a),nor the variations in the dilute limit of the PMF neara membrane surface ([21], Fig. 3) predicted by this ap-proach are in agreement with MD simulations. This ap-proximation also yields a very substantial, albeit seem-ingly undetected, direct contribution to the ionic free en-ergy of solvation,∆ G (cid:48) = − Φ (cid:48) pot ( z v ) = qφ lv . (1)on the order of 25 k B T . Disturbingly, the reasonableagreement between the experimental results for the sur-face tension of electrolyte solutions and certain promisingmesoscopic theoretical approaches that neglect the in-terface potential completely would be severely disruptedif such large interfacial potentials were taken into ac-count [23, 24]. This situation becomes even more com-plicated if one considers that more “realistic” quantummechanical calculations can lead to positive interface po-tentials of much higher amplitude (+3 eV) [12–14], butno signature of such a potential is seen in recent ab initio simulations of ion solvation [15]. III. IONIC FREE-ENERGY OF SOLVATION
The electrostatic (ES) contributions to the total ionsolvation free energy for the models studied here can beextracted directly from the MD simulations via∆ G ionES = ∆ G + ∆ G pol (2) FIG. 1: Electric potential variation, δφ ( z ), across the vapor-liquid and the liquid-solute interfaces for the neutral LJ soluteI immersed in liquid water (SPC/E), at z = 50 ˚A. The Gibbsdividing surface (GDS) of the macroscopic interface and thesolute position are indicated by dashed vertical lines. with ∆ G = qφ sv = q ( φ lv − φ ls ) (3)and ∆ G pol (cid:39) q ( φ ion sv − φ sv ) / φ sv = φ s − φ v = φ lv − φ ls and φ ion sv = φ lv − φ ion ls are,respectively, the total vapor-liquid-solute potential vari-ations for neutral and charged solutes. The ion solute-liquid interface potential, φ ion ls = φ ion l − φ ion s is the po-tential difference between the bulk liquid and the cen-ter of the charged ion (with the bare Coulomb potentialdue to the ion itself subtracted out). The first, direct ,term, ∆ G = qφ sv , may be regarded as the electrostaticfree energy of solvation of an ion placed in the poten-tial φ sv created around its corresponding neutral counter-part; the second polarization term, ∆ G pol , obtained froman approximate generalized “charging” method [29, 30],arises from the response of the solvent to the ion (whichgenerates an overall potential variation φ ion sv much largerin amplitude than φ sv ) [30]. Because φ lv is ion inde-pendent, the polarization contribution can be simplifiedto ∆ G pol (cid:39) q ( φ ls − φ ion ls ) / φ ls − φ ion ls ) ∝ q vanishes in the limit q → G pol ∝ q , as required for a polarization contribution(and also seen in the usual mesoscopic Born term below). A. Mesoscopic Born model
Within the mesoscopic Born model, an ion is modeledas a point charge q sitting in its spherical cavity of ef-fective radius R i in bulk water treated as a continuum of dielectric constant ε w . The radial electric potentialaround the central ionic charge, ϕ ion ( r ), determines theBorn approximation for the polarization contribution via∆ G Bpol = q r → [ ϕ ion ( r ) − ϕ ( r )] = q πε R i (cid:18) ε w − (cid:19) (5)where ε is the vacuum permittivity, ϕ ( r ) is the bareCoulomb potential, and ε w (cid:39)
78 is the dielectric con-stant of bulk water at room temperature. Although thistype of polarization contribution ( ∼ − k B T ) typi-cally dominates the ionic solvation free energy for the ionsstudied here, it is neither clear how accurate the simpleBorn approximation is (due to the neglect of potentiallyimportant ion-solvent correlations near the ion), nor howto choose best R i .Furthermore, despite its dominant role in the globalionic free energy of solvation, the polarization contribu-tion to the local ionic PMF, Φ( z ), is seemingly not thedominant contribution over a significant part of the inter-facial region, which means that the role of other contri-butions must be clarified [21]. Although a Born-type po-larization term is commonly incorporated in mesoscopicapproaches to the free energy of solvation (or PMF), thedirect term, arising from the bare interfacial potential, iseither completely neglected without justification [20, 23–25] or strongly overestimated by incorrectly assuming φ ls = 0 in ∆ G (eq. 3), which leads to the approxima-tion ∆ G (cid:48) = q φ lv (eq. 1) for the direct contribution [orΦ (cid:48) pot ( z ) in the PMF, which includes only the contributionof the vapor-liquid interface [21, 22] and neglects entirelythat of the solute-liquid one] (see fig. 1). IV. MODELS AND METHODS
We compare a direct evaluation from MD simulationsof the two terms contributing to ∆ G ionES , namely ∆ G and∆ G pol , for SPC/E water with the simplified approxima-tions presented above, respectively, ∆ G (cid:48) and ∆ G Bpol . Inorder to shed further light on the interplay between thesolvent molecular dipole and quadrupole moments (andthus water specificity), a series of molecular models hav-ing the same permanent dipole moment as SPC/E, butdifferent quadrupolar ones, were first generated by reduc-ing the H-O-H angle γ , while keeping fixed both the orig-inal partial charges on each site and the distance betweenthe oxygen and the midpoint between hydrogen atoms.The choice of including the variable quadrupole momentmodels was dictated by the need to find a smooth link viaMD simulations between a realistic water model (SPC/E)and the simplified purely dipolar models often studied(due to the inherent difficulty of the problem) using ap-proximate theoretical methods [31, 32]. Due to symme-try the interfacial contributions under scrutiny here forliquids possessing molecular dipole and quadrupole mo-ments must vanish for symmetric purely dipolar mod-els. We also would like to test the approximate formu-lae for the quadrupolar contribution in a more generalsetting (from SPC/E to a purely dipolar model) andto shed light on the coupling between the dipolar andquadrupolar contributions. For all but one molecularmodel the SPC/E parameters were maintained for theLennard-Jones (LJ) interaction centered on the oxygenatom. For the n th molecule of each liquid model we de-fine the molecular dipole moment p n = (cid:80) j q j r jn andquadrupole moment tensor ( Q n ) αβ = 3 (cid:80) j q j x jn,α x jn,β ,with x jn,α the α th Cartesian component of the positionvector r jn of partial charge q j ( j = 1 , ,
3) with respectto the center of charges within molecule n . We can thencompute the macroscopic polarization P ( r ) = (cid:68)(cid:88) n p n δ ( r − r n ) (cid:69) (6)and the macroscopic quadrupole moment density, Q αβ ( r ) = 16 (cid:68)(cid:88) n ( Q n ) αβ δ ( r − r n ) (cid:69) (7)directly from the simulations as ensemble averages [2].The local electric charge density, ρ ( r ), can be eval-uated directly by extracting the partial charge densityassociated with the particular molecular model; and theassociated electric field E and potential φ can then beobtained by integration of the Poisson equation in ap-propriate coordinates: ∇ φ = −∇ · E = − ρε . (8)The mean electric field along the direction normal to theplanar vapor-liquid interface, at distance z is written inCartesian coordinates as: E z ( z ) = z (cid:90) z v ρ ( z (cid:48) ) ε dz (cid:48) , (9)where ρ ( z (cid:48) ) represents the average electric charge den-sity (evaluated within the scope of MD simulations as thevolume density of the sum of partial charges associatedwith a particular molecular model found in the bin cor-responding to the position z (cid:48) ). The coordinate z v is theorigin of integration in the vapor phase (far from the in-terface). Similarly, in spherical coordinates (appropriatefor the curved solute-liquid interface), the mean electricfield is given by E r ( r ) = 1 r r (cid:90) ρ ( r (cid:48) ) r (cid:48) ε dr (cid:48) . (10)The total charge density can also be expressed as amultipole expansion [33], ρ = −∇ · P + (cid:88) α,β ∇ α ∇ β Q αβ + . . . , (11) here truncated after the quadrupolar term. The dipo-lar and quadrupolar electric fields and potentials canthen be obtained from the dipole moment density P andquadrupole moment density Q , respectively: ∇ φ D = −∇ · E D = 1 ε ∇ · P = − ε ρ D (12) ∇ φ Q = −∇ · E Q = − ε (cid:88) α,β ∇ α ∇ β Q αβ = − ε ρ Q (13)After computing the full interface potential from eq. 8, wecompare it with the sum of the dipolar and quadrupolarcontributions computed from eqs. 12 and 13 to assess theaccuracy of the truncated multipole expansion.The dipolar contribution to the total electric field canbe obtained from ρ D = −∇ · P , the mean dipolar chargedensity created by the distribution of the macroscopicdipole moment density P , yielding for the z -component: E Dz ( z ) = − P z ( z ) ε , (14)since P z vanishes at z v in the vapor phase, far from theinterface region. Similarly, the dipolar component of theradial field in spherical coordinates, appropriate for thesolute-liquid interface, reads: E Dr ( r ) = − P r ( r ) ε , (15)where P r ( r ) represents the radial distribution of the den-sity of the dipole moment. The dipole moment density P obtained from the MD simulations as ensemble averagesof molecular dipole moments is presented in detail belowfor the planar l-v interface (Section V A).The determination of the mean electric field (or chargedensity) permits the calculation of δφ ( z ) = φ ( z ) − φ ( z v ),the local electric potential difference evaluated at posi-tion z in the vicinity of the interface: δφ ( z ) = − z (cid:90) z v E z ( z (cid:48) ) dz (cid:48) = 1 ε z (cid:90) z v ρ ( z (cid:48) ) ( z (cid:48) − z ) dz (cid:48) (16)Similarly, the dipolar local potential profiles, δφ D ( z ), canbe obtained from the corresponding electric field, E Dz ( z )(or charge density, ρ D ).The quadrupole moments of the models SPC, SP9-SP5, and S2N range from the SPC/E value down to zero(table I). Because both dipolar S2N and S2L models areasymmetric, a symmetric dumbbell-like model (S2D) wasalso investigated (with two LJ spheres on both ends ofthe dipole). Ions are modeled as simple point chargescarrying an LJ sphere. Simulations of the vapor-liquid in-terface were carried out using a modified parallel versionof the molecular dynamics package Amber 9 [34] and aslab geometry methodology similar to the one often used TABLE I: The system parameters [ γ is the H-O-H angle; µ and Q xx are the permanent molecular dipole and quadrupolemoments; ρ l is the liquid density at the center of the slab (estimated error of ± .
005 g/cm )], total interfacial potential φ lv , thecorresponding quadru- ( φ Qlv ) and dipolar ( φ Dlv ) contributions, and the “isotropic” quadrupolar approximation ( φ Qlv, est , eq. 19);Estimated standard error of ± . γ ( ◦ ) µ (D) Q xx (D˚A) ρ l (g/cm ) φ lv (mV) φ lv − φ Dlv (mV) φ Qlv (mV) φ Qlv, est (mV) φ Qlv /φ lv SPC 109.5 2.347 8.131 0.981 -600.3 -558.8 -559.2 -558.6 0.932SP9 100.0 2.347 5.775 0.892 -445.8 -361.4 -360.1 -360.2 0.808SP8 87.6 2.347 3.736 0.787 -254.8 -206.5 -206.8 -205.9 0.811SP7 75.0 2.347 2.394 0.698 -123.6 -116.1 -116.9 -117.1 0.946SP5 54.7 2.347 1.089 0.611 -21.4 -46.3 -46.2 -46.6 –S2N 0 2.347 0 0.556 56.8 0.2 0 0 0S2L 0 4.065 0 0.696 284.6 -0.4 0 0 0S2D 0 2.347 0 0.658 -0.8 -0.4 0 0 0 in the literature: 1000 liquid molecules placed in a rect-angular unit cell of dimensions 31.04 × × occupying roughly the middle one-third of the availablespace and generating two vapor-liquid interfaces [35, 36].A Lennard-Jones interaction potential is centered on thesolvant oxygen atom, characterized by the SPC/E pa-rameters σ = 3 . ε = 0 . σ and ε and, when polarizable modelsare analyzed, a polarizability. The cross parameters forthe ion-water Lennard-Jones interaction are determinedvia Lorentz-Berthelot mixing rules.Periodic boundary conditions were applied in allthree directions. The long-range charge-charge, charge-dipole and dipole-dipole interactions were treated by theparticle-mesh Ewald summation method for both thecharge and dipole moments [41]. For computational ef-ficiency, in the polarizable simulations, an extended La-grangian method was utilized to compute the induceddipole moments, regarded as additional dynamic vari-ables [42].A cutoff radius of 10 ˚A was used for the short-rangednon-bonded LJ interactions and for the real space compo-nent of the Ewald summation. The geometries of the liq-uid molecules were constrained by applying the SHAKEalgorithm with a relative geometric tolerance of 10 − .The equations of motion were integrated using the veloc-ity Verlet algorithm with a default time step of 1 fs [43].In order to avoid occasional drifts of the slab along the z -axis normal to the interface, the center-of-mass (COM)velocity was removed every 1000 steps. Configurationswere saved every 100 fs in the output trajectories andeach such frame was readjusted with respect to the z -axis, to keep the COM of the electrolyte fixed relative to the simulation cell.Starting from the initial configuration of each simu-lated system, an energy minimization was performed,followed by a 1 ns NVT equilibration at 300 K for theslab systems. For simulations of bulk water ion solvationthe equilibration process was performed in the NPT en-semble, using a weak-coupling pressure regulation with atarget pressure of 1 bar. Subsequently, in both cases, atleast 5 ns of measurements in the NVT ensemble werecarried out, using the Berendsen thermostat with theconfigurational degrees of freedom coupled to a heat bathwith coupling constant τ = 1 ps [44]. In the special caseof polarizability-enabled simulations, the degrees of free-dom related to the induced dipole moment of the ionwere independently coupled to a 1 K heat bath (relax-ation time τ dip = 10 ps), ensuring a proper handling ofthe electronic degrees of freedom [45]. All computed pro-files spanning the vapor-liquid interface were obtained asensemble averages of the instantaneous profiles evaluatedin thin slabs (bins) of thickness 0.2 ˚A parallel to the inter-face. For the radial quantities measured in the bulk sim-ulations of ion solvation, equally distanced, 0.1 ˚A thick,spherical shell bins have been employed. Due to the cu-bic dimensions of the simulation box the radial profilesare relevant up to approximately 15 ˚A. V. MD SIMULATION RESULTS
The MD simulation results show that the bulk regiondensity decreases with decreasing molecular quadrupolemoment (at constant dipole moment), varying by nearlya factor of two in going from SPC/E to the lowest den-sity model, S2N (table I). For the models possessingquadrupole moments, SPC/E – SP5, the density de-creases by less than 40%, despite a decrease by a fac-tor of 8 in the molecular quadrupole moment. Thisdensity variation is an expected physical consequence
TABLE II: Lennard-Jones parameters σ and ε , and polarizability α for ionsIon q ( e ) σ (˚A) ε (kcal/mol) α ref.Na + − -1 3.168 0.2 0.974 [38]Cl − -1 4.339 0.1 3.25 [39]Br − -1 4.700 0.1 4.53 [40]I − -1 5.150 0.1 6.9 [40]FIG. 2: (Color on-line) The dipole moment density, P z ( z ),for all studied liquids at the vapor-liquid interface. Dashedvertical lines represent the GDS of both interfaces for SPC/E. of the reduction in water coordination as the molecu-lar quadrupole moment decreases. The S2N and S2Lmodels form purely dipolar liquids with a characteristicchain-like structure arising from the head-to-tail align-ment of the dipole moments and a lower bulk density dueto the decrease in hydrogen bond coordination from fourto two. The vapor-water interfacial thickness is foundto be approximately 3.8 ˚A at 300 K for SPC/E and in-creases along with the slab thickness as the molecularquadrupole moment decreases. A. Dipolar ordering
Orientational (dipolar) ordering of water takes placenear the l-v interface, which can be seen by plotting P z ( z ) for the series of model liquids (fig. 2). Orienta-tional double layers were found only for SPC/E and SP9with the outer layer dipoles pointing preferentially to-wards the vapor phase and in the opposite direction inthe inner layers (closest to the slab center). For modelswith lower molecular quadrupole moments the molecu-lar dipoles point towards the liquid phase. Because ofthe asymmetry created by the oxygen LJ sphere, asym- metric purely dipolar liquids (S2N and S2L) still pos-sess orientational ordering in the interface region due tothe hydrophobic forces tending to exclude the oxygen LJsphere from the liquid slab. B. Vapor-liquid interface potential: multipolecontributions
In order to illustrate how various multipole momentscontribute to the vapor-liquid interface potential, weobtained the electric potential difference by integrationof the Poisson equation from the charge density ob-tained from the first two terms of the multipole expan-sion [2, 10], ρ ( z ) ≈ ρ D ( z )+ ρ Q ( z ) = − ddz (cid:20) P z ( z ) − ddz Q zz ( z ) (cid:21) (17)leading to the first two (dipolar and quadrupolar) contri-butions to the interface potential: δφ DQ ( z ) ≡ δφ D ( z ) + δφ Q ( z ) = (cid:90) zz v P z ( z ) ε dz − Q zz ( z ) ε (18)since Q zz is taken to vanish in the vapor phase and P z vanishes in both bulk (vapor and liquid) phases. For theplanar interface higher order moments do not contribute.The “exact” model interface potential profile andthe corresponding dipolar and quadrupolar contributions(eq. 18) obtained directly from the simulation data (us-ing, respectively, the “exact” partial charge density, ρ ,and the multipole contributions of eq. 17) are illus-trated in fig. 3. The total vapor-liquid potential reaches − FIG. 3: (Color on-line) (a) Total, (b) dipolar and (c)quadrupolar potential profiles, δφ ( z ) (in volts) for the stud-ied liquids characterized by positive molecular quadrupolarmoments. moment and becomes comparable in absolute value withthe (positive) dipolar contribution for SP5, the near can-cellation in this case leading to a very low (negative)total value. For SPC/E, SP9-SP7 the quadrupole contri-bution provides the major contribution to the interfacepotential (figs. 3 and 4, table I) and thus to the largeinterfacial electric fields ( ∼ Q evaluated in a local referenceframe with the y -axis along the dipole vector and the z -axis out of the molecular plane: δφ Q est ( z ) = − c ( z )6 ε Tr Q , (19)where c ( z ) is the local liquid number density taken fromthe simulations. In this reference frame the only non-zerocomponent is Q xx and therefore in this case Tr Q = Q xx .The estimate φ Qlv, est = δφ Q est ( z l ) for the quadrupolar con-tribution is in excellent agreement with direct determi-nations (table I). We see from eq. 19 that the variationsin vapor-liquid interface potentials for the models withnon-zero quadrupole moments are mainly determined bythe quadrupolar contribution and therefore dominatedby the variations in the molecular quadrupole moments(with the physically relevant density variations playingonly a secondary role). For this reason and because thedipolar contribution is not strictly proportional to the liq-uid density, we do not attempt to normalize the vapor-liquid interface potentials in table I to correct for thevariations in liquid density. We have also checked thatthis simple quadrupolar estimation 19 provides a verygood approximation to the full oscillatory membrane-water surface potential [21], confirming that the air-waterinterfacial and membrane-water surface potentials aremainly determined by the local water density and molec-ular quadrupole moment. Furthermore, we propose thatthis method can be used to estimate the quadrupolar po-tential contribution of any liquid, irrespective of its bulkmolecular quadrupole moment, even those obtained fromab initio quantum mechanical calculations of liquid wa-ter. As an illustration, we have checked that when ab ini-tio values for Tr Q [46, 47], are injected into the simpleapproximation Eq. 19, we find quadrupolar potentials inreasonable agreement with the (quadrupole dominated)total interface potentials (+3 eV) extracted directly fromthe quantum mechanical calculations [13, 14]. C. Solute-Liquid interface potential
In the presence of a solute the planar vapor-liquid in-terface potential has as counterpart the microscopic po-tential between the solute center and the surroundingliquid (of which the first two multipole terms can be ob-tained from the microscopic analog of eq. 18 [30]).The local radial solute-liquid potential profile, definedas δφ r ( r ) = φ r ( r ) − φ r ( r s ) (20)[with φ r ( r s ) = 0 at the center of the solute r s = 0], is ob-tained from the radial charge density ρ r ( r ) by integrating FIG. 4: (Color on-line) Vapor-liquid interface potential dropsfor all studied models: standard (blue) and flipped (F)-charge(black) with respect to | Q xx | (molecular quadrupole moment). the Poisson equation (8) in spherical coordinates: δφ r ( r ) = − r (cid:90) E r ( r (cid:48) ) dr (cid:48) = − r (cid:90) ρ r ( r (cid:48) ) ( r (cid:48) ) ε (cid:18) r (cid:48) − r (cid:19) dr (cid:48) . (21)The dipolar component, δφ Dr ( r ), is determined from thedipole moment density P r ( r ) as: δφ Dr ( r ) = r (cid:90) P r ( r (cid:48) ) ε dr (cid:48) (22)The corresponding radial electric fields, E r ( r ) and E Dr ( r ), are determined from equations (10) and (15).The quadrupolar contribution to the total radial in-terface potential, δφ Qr ( r ) is accessible from the simula-tion data via the radial dependence of the quadrupolemoment density written in spherical coordinates, Q (cid:48) ( r ).We begin by writing the quadrupolar Poisson equation(13) in Cartesian coordinates (centered at the solute po-sition r s = 0), with the tensor elements of the Cartesianquadrupole moment density, Q αβ : ε ∇ · E Q = (cid:88) α,β ∇ α ∇ β Q αβ (23)Using the divergence theorem, we find: ε (cid:90)(cid:90) (cid:13) S E Q · dS = (cid:90)(cid:90) (cid:13) S (cid:88) α,β ∇ α Q αβ · n β dS, (24)where ˆ n is the normal to the interface. Letting F β = (cid:80) α ∇ α Q αβ and S = 4 πr , we obtain the radial depen-dence of the quadrupolar electric field E Qr ( r ), by inte-grating over the angular degrees of freedom:4 πε E Qr ( r ) = (cid:90)(cid:90) (cid:13) S F r sin θdθdφ, (25) with F r = e r · F (where e r is the unit radial vector). Afterperforming the angular integrals, we obtain the radialdependence of the quadrupolar electric field: ε E Qr ( r ) = ∂Q (cid:48) rr ∂r + 1 r (3 Q (cid:48) rr − T rQ ) . (26)Since E Qr ( r ) = − ∂φ Qr ( r ) ∂r and φ Qr (0) = 0, we arrive atthe final solution for the local quadrupolar potential atposition r , with respect to the center of the solute, δφ Qr = φ Qr ( r ) − φ Qr ( r s ), in terms of two components: δφ Qr ( r ) = δφ Q, r ( r ) + δφ Q, r ( r ) (27)with δφ Q, r ( r ) = − ε Q (cid:48) rr ( r ) (28) δφ Q, r ( r ) = 1 ε r (cid:90) dr (cid:48) r (cid:48) [Tr Q ( r (cid:48) ) − Q (cid:48) rr ( r (cid:48) )] . (29)The second contribution is generated by the symmetrybreaking of the diagonal components of Q (cid:48) in the solute-liquid interfacial region. The detailed calculations lead-ing from Q ( r ) to δφ Qr ( r ) will be presented elsewhere [30].Far from the solute center the radial solute-liquid inter-face potential and the various multipole components tendto their respective asymptotic values: φ ls = δφ r ( ∞ ), φ Dls = δφ Dr ( ∞ ), and φ Qls = δφ Qr ( ∞ ).Because the solute-liquid interface is curved, the as-sociated potential depends on the size of the cavity andis therefore ion specific. Thus, there is a first potentialdrop when going from the vapor into the liquid phase,followed by an overall increase near the solute, yieldinga smaller, overall negative, vapor-liquid-solute potentialdrop, φ sv (fig. 1). To obtain a fuller picture, we haveextracted the solute-liquid interfacial potentials (and thedipolar/quadrupolar components) from MD simulationsof neutral ion-like solutes, fixed at the center of a cu-bic box of bulk SPC/E water (using eqs. 21, 22,and 27).The solute-liquid potential φ ls and the correspondingdipolar contribution φ Dls are obtained, respectively, fromthe radial charge density and the dipole moment densitydistributions. The higher order multipole contributions, φ ls − φ Dls , are dominant and, due to ls interface curvatureeffects, not only is the amplitude of the ls quadrupo-lar contribution different from the lv one, but multipoleterms beyond the quadrupolar one play a non-negligiblerole (tables III and IV) [30]. For the halide-like neutralsolutes φ ls decreases in amplitude with increasing solutesize and is smaller than for the neutral sodium like solute,Na . For I , φ ls increases in amplitude by about 6% whenthe polarizability is turned on. Although the φ Q, ls con-tribution is dominant in φ ls , it simply serves, as we shallsee below, to cancel the large vapor-liquid quadrupolarcontribution to the direct contribution to the free energyof ion solvation, ∆ G . The dipole contribution φ Dls goesfrom negative values for small halide-like solutes (F ) topositive values for larger ones. D. Free energy of ion solvation
The solute-liquid potential was then used to evaluatethe direct contribution to the free energy of ion solvation,∆ G = q ( φ lv − φ ls ) (eq. 3), dominated by the differencebetween the quadrupolar lv and ls contributions, whichdo not cancel due to strong curvature effects for the smallions under study. Because the φ Q, ls component dependsonly on the bulk solvent properties, it is in principle equalto φ Qlv and therefore the two terms should cancel in ∆ G ,leaving the dominant ion specific quadrupolar contribu-tion, φ Q, ls (tables III and IV). The dipolar contributionand multipolar contributions higher than quadrupolarplay non-negligible, but secondary roles, in ∆ G . Ourresults for polarizable I − also reveal that ionic polariz-ability plays a minor but non-negligible role in determin-ing ∆ G (tables III and IV). The key point is that ∆ G is much smaller in amplitude than the simple direct es-timate ( | ∆ G (cid:48) | = 0 .
600 eV) because of the strong partialcancellation of the two interface potentials, φ lv and φ ls .Our results show that the direct contribution, ∆ G is re-duced to 30-40% of the simple direct estimate ∆ G (cid:48) anddepends on the sign of the ion charge in such a way asto favor the solvation of cations Na + with respect to theanions (halogens). For the halogens ∆ G is positive andincreases with increasing ion size, thus augmenting theinterface propensity of large anions like I − . Because theamplitude of ∆ G ( ∼ . − .
26 eV) is still between5% (for F − ) and 10% (for I − ) of the dominant polariza-tion contribution, ∆ G pol (see table V), it is on par withother important contributions (such as the hydrophobicone [22–24, 26, 48]) and therefore must be taken intoaccount correctly when considering ion-specific effects.Because of charge-dipole and charge-quadrupole cou-pling, when the ion charge is “turned on” the solventmolecules reorient themselves in order to accommodatethe solute, optimizing as much as possible the orienta-tion and number of hydrogen bonds: the ion polarizesthe medium around it and induces a radial potentialdifference, φ ion ls , dominated by the “long range” dipo-lar contribution. The difference between φ ion ls and φ ls extracted from the simulations then determines ∆ G pol (eq. 4). With our choice of effective ion radius, the simplemesoscopic Born approximation ∆ G Bpol is in reasonablygood agreement with the microscopic polarization contri-bution, ∆ G pol (table V). The polarization contribution tothe solvation free energy, which is the main effect favoringhigh ion solvation, decreases in amplitude with increasinghalogen size. We note also that the simulation results forthe electrostatic contribution to the ionic free energy ofsolvation, ∆ G ionES , follow the experimental trend, ∆ G ionexp ,reasonably well (table V), despite the neglect of certainentropic (hydrophobic) and enthalpic contributions (aris-ing, for example, from the short range repulsion and thelong range van der Waals—dispersion—attraction). VI. CONCLUSION
We have studied a series of liquid models that interpo-late between SPC/E water and pure dipolar liquids andshown that the quadrupolar component of the vapor-liquid interfacial potential typically dominates for thestudied liquids possessing a non-zero quadrupolar mo-ment. In an effort to elucidate the different ion-specificcontributions to the free energy of solvation, we haveshed light on the key role played by the solute-liquidinterface potential and demonstrated that it leads to astrong reduction in the direct electrostatic contributionwith respect to previous estimates based solely on thevapor-liquid potential.We propose that the same mechanism would be atplay if the point partial charge distribution of the sol-vent extracted from classical MD simulations were re-placed by the more realistic extended charge distribu-tions found in ab initio calculations. Indeed, a coarsegraining procedure for the electric potential proposed re-cently [14], which corrects for regions inaccessible to ionicprobes, shows that, encouragingly, both ab initio andpoint charge coarse grained potentials converge to val-ues that are compatible with the results for ∆ G = qφ sv presented in table V. Finally the dominant electrostaticpolarization contribution to the free energy of solvationwas found to agree reasonably well with a Born-type ap-proximation. We conclude that the direct interface po-tential contribution to the ionic free energy of solvation(or PMF) can neither be estimated using the point ionapproximation (leading to a gross overestimate), nor beneglected entirely – the two approximations commonlyadopted in the current literature. An important corol-lary that can be drawn from our study is that, in con-tradistinction to what is sometimes suggested [10], evenpurely quadrupolar liquids should give rise to an inter-face potential contribution to the ionic solvation free en-ergy because of the incomplete cancellation of the lv and ls components (due to solute curvature effects). Themechanism investigated here leading to a strongly re-duced vapor-liquid interfacial potential contribution tothe electrostatic part of the ionic PMF is quite generaland should be applicable not only to membrane-liquidsurfaces, but also other types of solvents and solutes.More complicated, possibly non-spherically symmetric,ions—like large organic ones—can be built for MD sim-ulations from several charged LJ particles and thereforean ionic cavity devoid of solvent will form and give riseto a solute-liquid interface potential.It would also be interesting, albeit difficult, to gener-alize the approximate theoretical statistical mechanicalapproaches developed previously for dipolar liquid mod-els near interfaces [31, 32] to quadrupolar liquids in or-der to capture the effects studied here via MD simula-tions. An important outcome of a reliable mesoscopictheoretical approach to the problem investigated herewould be a greatly enhanced comprehension of the un-derlying physics of ion distributions in inhomogeneous di-0 TABLE III: Solute-liquid interface potentials with multipolar contributions for neutral solutes solvated in SPC/E water (forwhich the vapor-liquid interface potential is φ lv = − . φ Qlv = − . φ Dlv = − . φ Q, ls = φ Qlv . The less than 2% difference between the values obtained from the water slab and bulk simulations can be attributedto the differences in densities generated by the finite slab thickness and the barostat used in the NPT bulk simulations (cf.eq 19). To directly compare solute-liquid and vapor-liquid interface potentials, we correct for this small systematic liquiddensity difference by defining a solute-dependent rescaled vapor-liquid interface potential, φ ∗ lv = Cφ lv , where C ≡ φ Q, ls /φ Qlv (rescaled multipole vapor-liquid interface components are defined similarly). The estimated standard error for the interfacepotentials is 0.15 mV and pol denotes polarizable.Solute φ ls /φ ∗ lv (%) φ ls (mV) φ Dls (mV) φ Qls (mV) φ Q, ls (mV) φ Q, ls (mV) ) φ ls − φ Dls − φ Qls (mV)Na (pol) 62.1 -380.1 14.4 -440.0 -570.5 130.5 45.5TABLE IV: Comparison of the different multipole contributions to the direct electrostatic solvation-free energies, ∆ G ∗ = q ( φ ∗ lv − φ ls ), obtained using the rescaled vapor-liquid interface potentials (*), of various ions as obtained directly from thesimulations (see table II) (∆ G (cid:48) (cid:39) ± . G ∗ (eV) Dipole ∗ (%) Quadrupole ∗ (%) Higher order ∗ (%)Na + -0.2030 -18.1 70.5 47.6F − − − − − (pol) 0.2327 24.4 56.0 19.6TABLE V: Comparison of electrostatic solvation-free energies for various ions as obtained directly from the simulations, ∆ G ionES =∆ G + ∆ G pol , with the estimate ∆ G ionES , est = ∆ G (cid:48) + ∆ G Bpol obtained from the truncated direct contribution [∆ G (cid:48) (cid:39) ± . G Bpol . R i is chosen as the distance from the ion centerwhere the ion-water radial distribution functions first reach 1 [30]. To correct for the small systematic differences in liquiddensity between water slab and bulk simulations, the rescaled vapor-liquid interface potentials, φ ∗ lv , are used in the evaluationof the direct contribution: ∆ G ∗ = q ( φ ∗ lv − φ ls ) (see table III). Estimated error of 0.025˚A for R i ; Estimated standard errors:0.002 eV for ∆ G ; 0.05 eV for the other solvation-free energies. (1 eV= 23 .
06 kcal/mol (cid:39) k B T )Solute R i (˚A) ∆ G ∗ (eV) ∆ G Bpol (eV) ∆ G pol (eV) ∆ G ionES , est (eV) ∆ G ionES (eV) ∆ G ionexp (eV) [46]Na − − − − − + + electric settings with important applications in colloidalscience, nanotechnology (ion transport in artificial nano-pores, or nano-filtration), and biophysics (ion channels,biological membranes, DNA) [16–19]. We also expect that the results presented here transcend the particularchosen models and thus qualitatively illustrate importantphysicochemical mechanisms at play in ion partitioningwithin inhomogeneous dielectric media.1After this work was completed we became aware ofother interesting very recent work covering similar topicsand in which some of the same conclusions were reached[50–53]. In [50] the same problem was studied and thesame interface potential reduction mechanism proposedfrom a different perspective: instead of investigating thevarious multipole contributions, as we do, a novel methodof partitioning the ionic solvation free energy was usedto extract from MD simulations of SPC/E water differ-ent physically identifiable (cavity formation, attractivevan der Waals, local and far-field electrostatic) contribu-tions (using recently optimized MD parameters [26, 54]).Interestingly, for I − , the loss of the first water hydra-tion layer (local electrostatic contribution), favoring ionsolvation, is nearly counterbalanced by the hydropho-bic (cavity formation) contribution favoring desolvation.A remaining net interface potential contribution, favor-ing anion desolvation, due to the competing vapor-liquidand solute-liquid interfaces, was obtained that is consis-tent with the results obtained here. The authors of [50]also proposed, as we did above, that this same interfacepotential reduction mechanism could be used to resolvethe apparent huge discrepancy between the classical andquantum predictions for the interface potential contribu-tion to the ionic solvation free energy. This scenario wasshown to be viable in [53], where both a detailed crit-ical comparison with other quantum simulations and afavorable experimental assessment were carried out. Theresults obtained in [53] are consistent with those of [14]showing that a coarse graining procedure that effectivelyomits certain regions of space in computing average in-terface potentials leads to a closer agreement betweenthe classical and quantum results. Despite these recentadvances concerning the interface potential contributionto the ionic solvation free energy, a definite comparisonwith experiment is for the moment complicated by whatappears to be some model and/or sampling dependency.In [51] the free energy of a single ion close to hydrophobicand hydrophilic surfaces was investigated using a noveltheoretical framework to obtain the position dependentdielectric response for interfacial water using moleculardynamics simulations. A multipole analysis for the pla-nar surface was then carried out, underlining the impor-tance of the quadrupolar contribution to this surface po-tential (consistent with the results presented here). Therole of the solute-liquid interface contribution was not,however, evoked in [51]. In [52] the driving forces for an-ion adsorption to the water vapor-liquid interface werestudied by comparing the results of MD simulations withthose of a simplified mesoscopic theoretical approach.Two different solvent models were used, SPC/E water and a symmetric purely dipolar liquid (the Stockmayermodel, similar to our S2D model, presented above). Onthe one hand, an extra electrochemical surface poten-tial contribution was needed in the mesoscopic theoreti-cal approach to explain the results of the SPC/E simula-tions and the magnitude and sign of this contribution areconsistent with the reduced electrostatic interfacial oneobtained here. On the other hand, no such extra elec-trochemical contribution was needed in the case of theStockmayer model, which is consistent with our resultthat there should be no interfacial electrostatic contribu-tion for symmetric purely dipolar liquids.The simulation results presented here should be notonly a useful building block in the quest for construct-ing the physically relevant mesoscopic components of theion free energies of solvation (or PMF), but also providea challenge for statistical mechanical approaches used toanalyze the subtle interplay between short range steric,dipolar, and quadrupolar interactions in determining theproperties of polar fluids (in particular the dipole distri-bution near interfaces) and their influence on ions.We are currently extending the approach developedhere for the ionic solvation free-energy to the study ofthe local ionic PMF near aqueous interfaces and surfaceswith and without the potentially important effect of wa-ter and ion polarizability [21, 23, 24, 35, 48, 49]. Wenote that unlike certain non-polarizable ion-water mod-els [26, 54], polarizable ones have not yet been properlyoptimized and therefore the relative weight of polariz-ability in driving large anions to interfaces and surfaces(compared with the interfacial electrostatic contributioninvestigated here) remains to be determined. Acknowledgments
We acknowledge financial support from the FrenchAgence universitaire de la francophonie (AUF) andANR (grant ANR-07-NANO-055-01, SIMONANOMEMproject) and the Romanian CNCSIS PN-II 502 and 506grants. This work is partly based on an ANR SIMO-NANOMEM project report posted on the web in 2010 and on the Ph.D. thesis of Lorand Horvath. We wouldalso like to thank B. Coasne for helpful discussionsand T. Beck, S. Kathmann, K. Leung, and R. Netz forcommunicating their work to us. We are also gratefulto R. Netz for helpful comments on an earlier version ofthe manuscript. [1] F. H. Stillinger and A. Ben-Naim J. Chem. Phys. , 4431(1967).[2] M. A. Wilson, A. Pohorille, and L. R. Pratt, J. Phys.Chem. , 4873 (1987). [3] M. A. Wilson, A. Pohorille, and L. R. Pratt, J. Chem.Phys. , 3281 (1988).[4] M. A. Wilson, A. Pohorille, and L. R. Pratt, J. Chem.Phys. , 5211 (1989). [5] V. P. Sokhan and D. J. Tildesley, Mol. Phys. , 625(1997).[6] L. X. Dang and T.-M. Chang, J. Phys. Chem. B , 235(2002).[7] E. N.Brodskaya and V. V. Zakharov, J. Chem. Phys. ,4595 (1995).[8] M. Paluch, Adv. Colloid Interface Sci. , 27 (2000).[9] M. Matsumoto and Y. Kataoka , J. Chem. Phys. , 3233(1988).[10] E. Harder and B. Roux, J. Chem. Phys. , 234706(2008).[11] I. Vorobyov and T. W. Allen, J. Chem. Phys. , 185101(2010).[12] S. M.Kathmann, I.-F. W. Kuo, and C. J. Mundy, J. Am.Chem. Soc. , 17522 (2008).[13] K. Leung, J. Phys. Chem. Lett. , 496 (2010).[14] S. M. Kathmann . W. Kuo, C. J. Mundy, and G. K.Schenter, J. Phys. Chem. B. , 4369 (2011).[15] M. D. Baer and C.J. Mundy, J. Phys. Chem. Lett. ,1088 (2011).[16] W. Kunz, Specific ion effects (World Scientific Publish-ing, Singapore, 2010).[17] P. H¨unenberger, M. Reif, Single-Ion Solvation: Experi-mental and Theoretical Approaches to Elusive Thermo-dynamic Quantities (Royal Society of Chemistry , Cam-bridge, UK, 2011).[18] J. N. Israelachvili, Intermolecular and Surface Forces,2nd edition, (Academic Press, London, 1992).[19] R. J. Hunter Foundations of Colloid Science (Oxford Uni-verstiy Press, Oxford, 2002).[20] M. Bostr¨om, W. Kunz, and B. W. Ninham, Langmuir , 2619 (2005).[21] D. Horinek and R. R. Netz, Phys. Rev. Lett. , 226104(2007).[22] D. M.Huang, C. Cottin-Bizonne , C. Ybert, and L. Boc-quet, Langmuir , 1442 (2008).[23] Y. Levin, Phys. Rev. Lett. , 147803 (2009).[24] Y. Levin, A.P. dos Santos, and A. Diehl, Phys. Rev.Lett. , 257802 (2009).[25] M. Manciu and E. Ruckenstein, Langmuir , 11312(2005).[26] D. Horinek, A. Herz, L. Vrbka, F. Sedlmeier, S. I. Ma-matkulov, and R. R. Netz, Chemical Physics Letters
173 (2009).[27] S. Buyukdagli, M. Manghi, and J. Palmeri, Phys. Rev.Lett. , 158103 (2010).[28] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma,J. Chem. Phys. , 6269 (1987).[29] D. A. McQuarrie, Statistical Mechanics, 2 nd edition(University Science Books, Sausalito, California, 2000), Chap.15.[30] L. Horv´ath, T. A. Beu, M. Manghi, and J. Palmeri, inpreparation.[31] Peter Frodl and S. Dietrich, Phys. Rev. E , 3741(1993).[32] A. Abrashkin, D. Andelman, and H. Orland, Phys. Rev.Lett. , 077801 (2007).[33] J. D. Jackson, Classical Electrodynamics, 3 rd edition(John Wiley & Sons, New York, 1999).[34] D. A. Case et al., AMBER 9 (University of California,San Francisco, 2006).[35] P. Jungwirth and D. J. Tobias, Chem. Rev. , 1259(2006).[36] T.-M. Chang and L. X. Dang, Chem. Rev. , 1305(2006).[37] D. E. Smith and L. X. Dang, J. Chem. Phys. , 3757(1994).[38] L. Perera and M. L. Berkowitz, J. Chem. Phys. , 3085(1994).[39] L. X. Dang, J. Phys. Chem. B , 10388 (2002).[40] G. Markovich, L. Perera, M. L. Berkowitz, and O. Chesh-novsky, J. Chem. Phys. , 2675 (1996).[41] T. Darden, D. York, and L. Pedersen, J. Chem. Phys. , 10089 (1993).[42] A. Toukmaji, C. Sagui, J. Board, and T. Darden, J.Chem. Phys. , 10913 (2000).[43] M. P. Allen and D. J. Tildesley, Computer Simulation ofLiquids (Oxford University Press, New York, 1987).[44] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gun-steren, A. DiNola, and J. R. Haak, J. Chem. Phys. ,3684 (1984).[45] M. Sprik, J. Phys. Chem. , 2283 (1991).[46] A. Kumar, J. Phys. Soc. Jap. , 4 (1992).[47] K. Leung, private commuincation (2011).[48] G. Archontis and E. Leontidis, Chem. Phys. Lett. ,199 (2006).[49] P.-A. Cazade, J. Dweik, B. Coasne, F. Henn, and J.Palmeri, J. Phy. Chem. C , 12245 (2010).[50] Ayse Arslanargin and Thomas L. Beck, J. Chem. Phys. , 104503 (2012).[51] Douwe Jan Bonthuis, Stephan Gekle, and Roland R.Netz, Langmuir , 7679 (2012).[52] Marcel D. Baer, Abraham C. Stern, Yan Levin, DouglasJ. Tobias, and Christopher J. Mundy, J. Phys. Chem.Lett. , 1565 (2012).[53] Thomas L. Beck, Chemical Physics Letters , 1(2013).[54] D. Horinek, S. Mamatkulov, and R. Netz, J. Chem. Phys.130