Hydrated Excess Proton at WaterâHydrophobic Interfaces - The

J. Phys. Chem. B 2009, 113, 4017–4030

4017

Hydrated Excess Proton at Water-Hydrophobic Interfaces† Satoru Iuchi, Hanning Chen, Francesco Paesani, and Gregory A. Voth* Department of Chemistry and Center for Biophysical Modeling and Simulation, UniVersity of Utah, 315 South 1400 East, Room 2020, Salt Lake City, Utah 84112-0850 ReceiVed: June 16, 2008; ReVised Manuscript ReceiVed: August 9, 2008

The behavior of the hydrated excess proton is investigated at the water-vapor, water-hydrophobic wall, and water-carbon tetrachloride interfaces through molecular dynamics simulations using the third-generation multistate empirical valence bond model (MS-EVB3). The MS-EVB3 simulations show a surface preference of the excess proton at the water-vapor interface, consistent with the discovery of this effect using an earlier version of the MS-EVB model (Petersen et al. J. Phys. Chem. B 2004, 108, 14804) and with the experimental results. The preference of the hydrated excess proton for other water-hydrophobic interfaces is also analyzed for the first time. The hydrated proton structures and charge defect delocalization effects at these interfaces are discussed in detail. By decomposing the free energy profiles into the internal energy and entropic contributions, the thermodynamic (free energy) driving forces for the surface preference of the excess proton are also elaborated. These results indicate that the “rigid” hydrated proton structures at all the interfaces are energetically (as opposed to entropically) stabilized due to the “amphiphilic” nature of the hydrated excess proton, resulting in its overall interfacial concentration enhancement. The effects of acid concentration and nuclear quantization are also explored. 1. Introduction Proton solvation and transport in aqueous environments has been investigated intensively due to its importance in chemical and biological phenomena. On theoretical grounds, molecular dynamics (MD) simulation has been a powerful tool to investigate the molecular level details of such phenomena. To obtain an accurate description of the hydrated excess proton behavior through MD simulations, however, it is important to incorporate the inherent delocalized nature of the excess proton charge defect over the water hydrogen bond (HB) network. The multistate empirical valence bond (MS-EVB) method, which includes the delocalization effects of the excess proton, has successfully provided molecular level mechanisms of proton solvation and transport in numerous systems (see refs 1 and 2 and references therein). Ab initio MD simulations also take into account the delocalized nature of the excess proton within the framework of density functional theory calculations. As a result, they have also provided molecular level insights into these phenomena (see, e.g., ref 3), although for substantially smaller systems of much shorter simulation time scales given their substantial computational cost. As an outgrowth of these efforts, our group has proposed that an “amphiphilic” nature of the hydrated proton may play an important role in various chemical and biological phenomena.4-7 The hydronium ion, formed by a water molecule and an excess proton, has a pyramidal equilibrium structure with three O-H bonds and a lone pair region. The three O-H bonds can form HBs with three water molecules, thus exhibiting a strong hydrophilic character. On the other hand, the lone pair region does not form HBs with waters due to the overall positive charge moiety in the cation, resulting in a “hydrophobic” character. Due to this unusual amphiphilic character of the hydronium ion, † Part of the special section “Aqueous Solutions and Their Interfaces”. * To whom correspondence should be addressed. E-mail: [email protected].

the hydrated excess proton was predicted to show a preference for the wate-vapor interface. This behavior was first demonstrated by Petersen et al. with an MS-EVB simulation,4 but it was also later supported by several types of MD simulations8-12 and experimental evidence.13-16 This result also has important implications for the understanding of numerous phenomena. For example, interesting behavior attributed to the amphiphilic nature of hydrated proton has been observed through MD simulations in large water clusters,7,17,18 methanol-water solutions,5 and concentrated acid solution.6 The prediction of enhanced hydronium concentration at the aqueous liquid-vapor interface has also led to some degree of controversy. For example, on the basis of the earlier prediction of Petersen et al.,4 other authors have proposed that even neutral water at pH 7 has an effectively lower pH at the liquid-vapor interface by more than 2 pH units10,11 (i.e., a more acidic character). This latter proposal has been vigorously challenged on the basis of experimental evidence,19 which in turn has been rebutted20 and further questioned on the basis of additional experimental evidence.21 While this controversy over the effective pH of the neutral water surface has not yet been resolved, it should be understood that the original simulation of Petersen et al. predicting the enhancement of hydronium concentration at the aqueous liquid-vapor interface was carried out for a 0.055 M HCl aqueous solution,4 which is quite different from neutral pH conditions.22 This prediction of an enhanced hydronium concentration at aqueous liquid-vapor interfaces has also been the subject of some confusion. For example, it is sometimes implied via referencing10,11,20 that the interfacial enhancement effect of hydronium was first seen in studies of small protonated water clusters having 20 or so water molecules.23,24 In fact, while the excess proton is indeed found on the surface of such small clusters, these clusters are essentially all surface and have no “bulk” region. The small protonated water cluster results, while very interesting in their own right, are not clearly relevant to

10.1021/jp805304j CCC: $40.75  2009 American Chemical Society Published on Web 09/27/2008

4018 J. Phys. Chem. B, Vol. 113, No. 13, 2009 the enhancement of hydronium at the aqueous liquid-vapor interface and at other condensed-phase mixed dielectric interfaces. Considering this current state of affairs, it is important to re-examine and better understand the behavior of the hydrated excess proton at water-hydrophobic interfaces, which may also provide valuable insight into proton solvation and transport in systems containing hydrophobic regions such as fuel cell membranes and biomolecular systems such as proteins. In fact, understanding hydrophobicity in general is an important subject for various chemical and biological phenomena.25-27 Although the water-vapor interface containing the excess proton has been studied through MD simulations as mentioned above, a detailed free energy analysis to determine the origin of the observed effects has never been performed. Moreover, the hydrated proton at other water-hydrophobic interfaces has not been well studied except for an ab initio MD study of limited size and duration.28 In this paper, therefore, the behavior of the hydrated excess proton is systematically investigated through MS-EVB MD simulations at the liquid-air (water-vapor), liquid-solid (water-hydrophobic wall), and liquid-liquid (water-carbon tetrachloride) interfaces. For this purpose, the recently refined third-generation MS-EVB model (MS-EVB3)29 is utilized. Effects of proton concentration and of nuclear quantization are also briefly discussed. The use of the MS-EVB3 model is necessary due to the inherent delocalized nature of the charge defect associated with the excess proton. The delocalization of the charge defect is an inherent property due to the Grotthuss shuttling of the excess proton in aqueous systems.1-3 Classical hydronium cation models (i.e., H3O+ with no possible Grotthuss shuttling) omit this essential physics from the simulation.8-12 Moreover, unlike ab initio MD simulations, long-time simulations for sufficiently large systems can be performed using the MS-EVB model.1-5 Therefore, free energy profiles can be obtained to analyze the excess proton behavior at the interfaces. The water-vapor interface can be regarded as a hydrophobic interface in the sense that the vapor (air) region has a low dielectric constant (ε ) 1) and obviously does not allow an HB network. For example, an experimental study using the vibrational sum frequency generation (SFG) spectra showed that various hydrophobic interfaces, including the water-air, water-surfactant-coated quartz, and water-hexane surfaces, are all characterized by water dangling bonds.30 Therefore, one might also see a preference of the hydrated excess proton at water-hydrophobic interfaces as is observed for the water-vapor interface. In addition, one might simply expect to see a contact between the hydronium ion and the hydrophobic surface due to the hydrophobic region of the hydronium ion.4 However, a simulation result by Wallqvist and Berne31 demonstrated that a solvent-separated pair between a methane molecule and hydrophobic wall is more likely than the contact pair between them. It is therefore possible that different types of hydrophobic surfaces affect rearrangements of water and hydrated protons at interfaces in different ways, resulting in a different excess proton hydration structure compared to that at the water-vapor interface.Comparisonsofresultsforthreedifferentwater-hydrophobic interfaces can also provide useful insight into the hydrated excess proton behaviors in various systems including membranes and proteins. In this study, it is found that the hydrated proton structures are formed with the lone pair region of the hydronium ion being pointed outward from the bulk water region at all the interfaces examined and that the excess proton charge defect delocalizes over such structures. The latter behavior is not captured by classical hydronium models, and it contributes a

Iuchi et al. non-negligible fraction of the free energy of hydrated proton stabilization at the interface. In addition to the detailed analysis of hydrated proton structures at the three hydrophobic interfaces, the thermodynamic mechanisms to determine the surface preference are investigated in the present work. It has been discussed that a favorable energy contribution coming from the surface stabilization of the excess proton may be a major source of the preference for the water-vapor interface.10,11 However, the free energy components have not been studied consistently using a physically accurate model that includes the charge delocalization. In addition, it is not clear whether the surface preference at other hydrophobic interfaces can be explained in the same fashion. In this paper, therefore, the free energy profiles are decomposed into the internal energy and entropic contributions. It is shown from such decomposition analyses that forming rigid hydrated proton structures is associated with an entropic penalty at all the interfaces, but such hydration structures are more than compensated for by the energetic stabilization, resulting in an overall surface (free energy) preference for the hydrated excess proton. The present analyses on both interface solvation structures and free energy contributions confirm the conclusions of an earlier study with the second-generation MS-EVB (MSEVB2) model4 that the amphiphilic nature of the hydronium ion plays a major role in the preference of the excess proton for the water-vapor interface. In addition, this concept is generalized here to the cases of the water-hydrophobic wall and water-CCl4 interfaces. Interestingly, the amphiphilic effect has an energetic, as opposed to entropic, origin. This paper is organized as follows. In section 2, details of the simulation conditions using the MS-EVB3 model are described. In section 3, the detailed analyses for the water-vapor interfaces are presented, and the effects of proton concentration and nuclear quantum effects are briefly discussed as well. In sections 4 and 5, the results for the water-hydrophobic wall and water-CCl4 interfaces are discussed, respectively. Conclusions are given in section 6. 2. Simulation Details 2.1. MS-EVB3 Simulations. Since the details of the MSEVB3 model are described in the original paper,29 only the simulation details related to this work are described in this section. In the MS-EVB3 framework, as in the earlier versions of the MS-EVB models,32,33 each EVB state represents a possible classical hydronium ion configuration in waters. The EVB Hamiltonian is constructed on the basis of those EVB states, and the lowest energy state obtained from the diagonalization of the EVB Hamiltonian matrix corresponds to the adiabatic potential energy surface of the protonated system that includes the full physics of Grotthuss shuttling and the resulting charge defect delocalization. The eigenvector of the lowest energy state is given by a set of coefficients, c ) (c1, c2,..., cNEVB), where NEVB is the number of EVB states, and the amplitude cI2 represents the weight of the Ith EVB state.1 The center of excess charge (CEC) is defined as NEVB

rCEC )

∑ cI2rCOC I

(1)

I)1

where rICOC represents the center of charges for the Ith EVB state.1,33 The CEC position reflects the spatial delocalization of the excess proton, and thus, it is used to investigate properties of the excess proton in this work. The highest and second highest amplitudes, c12 and c22, reflect the degree of delocalizations

Excess Proton at Water-Hydrophobic Interfaces along the water HB network: (c12, c22) ) (0.5, 0.5) represents the Zundel cation, H5O2+, and a larger value of c12 (∼ 0.7) corresponds to the Eigen cation-like complex, H9O4+.1,32,33 The hydronium ion in the highest amplitude EVB state is denoted as a “pivot” hydronium ion, and the oxygen and hydrogen atoms of the pivot are denoted by O* and H* hereafter. As stated before, three kinds of water-hydrophobic interfaces are investigated in the present study: (i) water-vapor, (ii) water-hydrophobic wall, and (iii) water-CCl4 interfaces. Unless stated otherwise, in all of the interface systems, 500 water molecules, 1 excess proton, and 1 chloride anion were contained in a rectangular box with periodic boundary conditions for all the directions. These simulation conditions correspond to an about 0.1 M HCl solution (i.e., pH 1 conditions). The parameters for the interactions for the MS-EVB3 model were taken from the original paper, and the interactions coming from Cl- were taken into account through the partial charge of -1.0 and the 12-6 Lennard-Jones (LJ) parameters, σCl ) 4.4 Å and εCl ) 0.1 kcal/mol.6,34 In addition, LJ interactions between Cl- and water hydrogen (HW) and between Cl- and hydronium ion hydrogen (HT) were added.6,34 The standard 3D Ewald sum method35 was used for the electrostatic interactions. The real space electrostatic and LJ interactions were cut off at 12.0 Å. Constant NVT simulations at 298.15 K were performed to calculate nondynamical properties using the Nose´-Hoover thermostat36 with a relaxation time of 0.5 ps. A time step of 1 fs29 was used with the leapfrog Verlet algorithm.35 From 10 different initial configurations, at least a 500 ps run was performed for equilibration, followed by a 1 ns run in each simulation (10 ns in total) with configurations saved every 100 fs. For dynamical properties, constant NVE simulations were performed using a time step of 0.5 fs. From 10 configurations obtained from the constant NVT simulations, a ∼100 ps equilibration run was carried out, followed by a 250 ps production run in each simulation (2.5 ns in total) with configurations saved every 10 fs. The average temperatures in all the constant NVE simulations were ∼298 K. To examine the charge defect delocalization effects at the interfaces, simulations for the three interface systems were carried out without allowing the excess proton to delocalize over the water HB network. This was done using a “one-state” MSEVB3 model in which the repulsive interaction terms between waters and the hydronium ion in the EVB diagonal elements (eqs 7 and 8 in ref 29) were set to zero.37 These simulations are denoted later as ones using a “classical” hydronium ion model. The hydronium oxygen and hydrogen atoms are denoted as OT and HT, respectively. The code for the MS-EVB3 model29 modified from DL_POLY version 2.13 software38 was used for all the simulations in this work. 2.2. Water-Vapor Interface. For this system, a “slab” geometry was used.4 The system size was 24.54 × 24.54 × 75.0 Å3, where a length of 24.54 Å was determined from the liquid density of 1.012 g/cm3 for the SPC/Fw water model39 utilized in the MS-EVB3 model. It was shown that the use of the 3D Ewald sum method does not yield artifacts when a box length in the z direction is sufficiently large.40 The origin z ) 0 is set to be the center of mass (COM) of the system. The excess proton behavior was investigated in most systems containing a single excess proton. However, to examine the effects of higher acid concentrations, self-consistent iterative (SCI) MS-EVB3 simulations6,34 were also carried out for a 1.1 M HCl solution liquid-vapor interface. Ten excess protons and Cl- anions were put into the same size box as in the single-

J. Phys. Chem. B, Vol. 113, No. 13, 2009 4019 proton simulations. From a 5 ns simulation using the classical hydronium model, 10 configurations were picked. From these 10 initial configurations, the SCI-MS-EVB3 simulations were performed. After a ∼100 ps run for equilibration, ten 120 ps production runs (1.2 ns in total) were performed from each initial configuration with configurations saved every 10 fs. The effect of nuclear quantization on the surface preference of the excess proton was also investigated using path integral molecular dynamics (PIMD) simulation.41,42 Due to the large computational overhead, a smaller system composed of 255 waters, 1 Cl-, and 1 classical hydronium ion was used. The interactions were based on the qMS-EVB3 model,29 which was designed for quantum MS-EVB3 simulations. The system size was 19.72 × 19.72 × 75.0 Å3, and a cutoff length of 9.85 Å was used for the real space electrostatic and LJ interactions. The time propagation was based on the normal mode representation of path integrals42 using the velocity Verlet algorithm35 with a time step of 0.2 fs. Four Nose´-Hoover chains were attached to each path-integral normal mode to realize the canonical (constant NVT) ensemble efficiently,42 and 24 quasiparticles were used to represent each atom. From five initial configurations, ∼250 ps PIMD equilibration runs were performed, followed by a 600 ps production run in each simulation (3.0 ns in total) with configurations saved every 100 steps. For a direct comparison, the classical MD simulations using the MSEVB3 parameters were also performed with the same conditions as in the larger system except for the cutoff length and system size, which were the same as in the PIMD simulations. 2.3. Water-Hydrophobic Wall Interface. Several MD simulations have been previously performed for investigating the water-hydrophobic wall interface system.31,43-46 In the present work, the hydrated excess proton behavior is investigated at the water-solid hydrophobic interface using the model wall developed by Lee et al.43 In this hydrophobic wall system, it is assumed that molecules interact with uniformly distributed hydrophobic molecules through the usual 12-6 LJ interactions, and the z-dependent interaction between them is derived to be43,47

Uwall(z) ) ≡

4πFεσ12 2πFεσ6 45|z - zw|9 3|z - zw|3 A C 9 |z - zw| |z - zw|3

(2)

where σ and ε are the 12-6 LJ parameters for interactions with hydrophobic molecules. The wall origins, zw, were placed at (14.744 Å so that the “effective width” 43 of the present system was 24.54 Å. The system size was the same as in the water-vapor interface, 24.54 × 24.54 × 75.0 Å3. It has been shown that the use of the 3D Ewald sum method is one of the best choices for the pure water-wall interface system.45 The origin of z ) 0 corresponds to the middle position between the two walls. Values of A and C for the water oxygen (OW)-wall interaction were directly taken from the literature.43 Coefficients A and C for other interactions, summarized in Table 1, were determined through eq 2 on the basis of the usual combination rule for the 12-6 LJ interaction35 so as to be consistent with the OW-wall interaction. It is noted that the HW-wall interactions were not included in some of the simulations, consistent with the other simulation work.31,43-46 In a different simulation study for a water-hydrophobic interface, however, the HW-wall interaction was employed to avoid nonphysical penetration of hydrogen atoms into the wall.48 Considering this, MD simulations including the hydrogen-wall interactions were

4020 J. Phys. Chem. B, Vol. 113, No. 13, 2009

Iuchi et al.

TABLE 1: Parameters for Atom-Wall Interactions atom symbol

atom

OW OT ClHW HT

water oxygen hydronium oxygen chloride anion water hydrogen hydronium hydrogen

A [kcal/mol 4170.0a 3061.0 23970.0 618.8b 458.9b

Å9

]

( ∂∆F(z) ∂T )

-T∆S(z) ) T 3

C [kcal/mol Å ] 18.2a 13.92 39.08 4.957b 3.810b

a Taken from ref 43. b Most of the simulation results presented in this work do not include hydrogen-wall interactions. The details are described in section 2.3.

also performed. The parameters for such interactions, included in Table 1, were determined by following the procedure in ref 48. As discussed later, it was found that these interactions do not significantly affect the preference of the hydrated excess proton for the water-hydrophobic wall interface. Therefore, only the simulation results without the hydrogen-wall interactions will be discussed unless otherwise noted. 2.4. Water-CCl4 Interface. The rectangular system, 24.54 × 24.54 × 51.32 Å3, determined by a density of 1.58 g/cm3 for liquid CCl449 contained 500 waters, 1 excess proton, 1 chloride anion, and 100 CCl4 molecules. The origin of z ) 0 was set to be the middle position of the water density profile. For the CCl4 molecules, a nonpolarizable and flexible five-site model was employed50 so as to be consistent with the MS-EVB3 water model, which is also nonpolarizable. 2.5. Umbrella Simulations. Because the hydrated proton shows an apparent preference for the interfaces, the MS-EVB simulations without any constraints do not provide sufficient sampling of configurations where the CEC resides in the bulk phase. To sample these configurations more efficiently and to obtain reasonably accurate free energy profiles, biased MSEVB3 simulations with an umbrella potential were performed. A harmonic umbrella potential was applied to the CEC position only in the z direction with a force constant of 5 kcal/mol. The technical details of CEC umbrella simulations are described in ref 51. In addition, the umbrella potential is not a pairwise interaction, and thus, it produces a nonzero total sum of forces, resulting in the system COM drift. To avoid this drift, the nonzero total sum of forces was uniformly redistributed onto all atoms and then subtracted from forces acting on all atoms.52 This procedure was not applied to the water-hydrophobic wall system, where the wall potential is not a pairwise interaction. However, the origin in the z direction can be clearly defined as the middle position between the two walls. A 1.5 ns production simulation (3.0 ns for the water-CCl4 interface) after a 200 ps equilibration run was performed in each window at 278.15, 298.15, and 318.15 K. The windows were centered on positions from z ) 1.0 Å to interface regions with an increment of 1 Å. The CEC z positions were saved every five steps and were then used to determine the free energy profiles. Other simulation conditions were the same as in the case of the unbiased simulations. The potentials of mean force (PMF), i.e., the relative Helmholz free energy profiles, ∆F(z), were reconstructed using the weighted histogram analysis method (WHAM).53,54 The bulk regions (average values from z ) 0.0 to z ) 4.0 Å) were taken as the common origin of the free energy profiles for all the interfaces. The standard errors were calculated by dividing trajectories into five blocks in each window. The free energy profile, ∆F(z), consists of two contributions: the internal energy change, ∆E(z), and the entropic contribution, -T∆S(z). The entropic contribution can be determined through the derivative of the free energy profile:

N,V

(3)

The derivative of free energy at 298.15 K with respect to temperature in eq 3 was evaluated by a finite-difference method55,56 on the basis of the free energy profiles at 278.15 and 318.15 K. The errors for the entropic contributions were estimated by using those of ∆F(z) at the two temperatures.56 To discuss the probability of the EVB amplitudes in different slab regions, the umbrella simulations were also utilized. From the biased conditional probability, 〈F(cI2;zCEC)〉biased, the twodimensional probability, 〈F(cI2,zCEC)〉, can be calculated via57,58

〈F(cI2, zCEC)〉 ) 〈F(zCEC)〉〈F(cI2 ;zCEC)〉biased

(4)

More statistically reliable probability distributions for the EVB amplitudes with respect to zCEC can be obtained by this procedure than the unbiased simulations for wider range regions. At the water-hydrophobic wall interface, interactions from the wall contribute to the free energy profile. It is useful to remove the contribution of interactions between the wall and the proton to see only the solvation effects.31 Because the CEC is not an interaction site, it is easier to employ the O* free energy profiles to eliminate the wall contribution. The O* free energy profiles were determined as follows. First, the biased probability distribution, 〈F(zO*,zCEC)〉biased, was determined from the umbrella simulations in each window. Then, the unbiased probability distribution was calculated using the two-dimensional WHAM54 with a zero umbrella potential for the O* position and a umbrella potential of 5 kcal/mol for the CEC position. The final O* unbiased probability distribution can be calculated by

〈F(zO*)〉 )

∫ dzCEC〈F(zO*, zCEC)〉 ∫ ∫ dzO* dzCEC〈F(zO*, zCEC)〉

(5)

Because the O* positions are always located near the CEC positions without an explicit umbrella potential for the O* position, the probability distribution 〈F(zO*,zCEC)〉 has nonzero values only around zCEC, and thus, a reliable O* distribution is obtained except for two edge regions (around 0.0 and 12.5 Å) where nonzero distributions are not included in the integrations. The wall contribution was then eliminated from the O* free energy profiles ∆F(zO*) determined from eq 5. The O* atom can be both a hydronium oxygen and a water oxygen in the MS-EVB method, so the wall contribution in the free energy profile is given by a linear combination of the two wall potentials as wall ∆Ewall(z) ) RUOT (z) + (1 - R)Uwall OW (z)

(6)

where R ) 0.6 was used. The details of this procedure are described in the Appendix. 3. Results: Water-Vapor Interface It was demonstrated from the MS-EVB2 simulation that the excess proton shows a preference for the water-vapor interface due to its amphiphilic character of the hydrated proton.4 In this section, such a preference is discussed in more detail on the basis of the MS-EVB3 simulation. 3.1. Interface Excess Proton Solvation Structures. Figure 1 shows the z-dependent density profiles for OW, CEC, Cl-, and O*, where the origin of z ) 0 corresponds to the COM of the system. Since the two interfaces in the positive and negative z regions are equivalent by symmetry, the density profiles were averaged with respect to z ) 0 to obtain better statistics. The

Excess Proton at Water-Hydrophobic Interfaces

J. Phys. Chem. B, Vol. 113, No. 13, 2009 4021

Figure 1. Normalized density profiles of hydrated excess proton CEC, OW, Cl-, and O* for the water-vapor interface. Dashed lines represent results from the classical hydronium model. The origin of z ) 0 corresponds to the center of mass of the system.

Figure 2. Probability distributions of angle φ between the surface normal and O-H bond of the pivot hydronium ion in different interface layers: |zO*| ) 13.5-15.0 Å (solid line), 12.0-13.5 Å (dashed line), 10.5-12.0 Å (dashed-dotted line), and 7.5-9.0 Å (dotted line).

density profiles at the water-vapor interface are qualitatively similar to those obtained from the earlier MS-EVB2 simulation.4 The density profiles for both the excess charge defect CEC and O* clearly show the preference of the excess proton at the water-vapor interface, consistent with behavior inferred from experimental evidence.13-16 In contrast to the surface preference of the excess proton, the chloride anion has no distinctive surface preference, in agreement with recent experimental observations that show a small effect of the sodium and chloride concentrations on the SFG spectra.14,59 The current MS-EVB3 model does not have explicit electronic polarization effects for all of the molecules in the system. Nevertheless, the simulation results clearly show the apparent surface preference of the excess proton, implying that, although polarizabilities of the hydronium ion and/or of water may enhance the surface preference, the preference is not dominated by this factor. Indeed, the classical hydronium model also shows the same surface preference, though the degree is smaller as observed in Figure 1. The classical hydronium model with the MS-EVB2 parameters also showed this surface preference.4 It should be noted, however, that polarization effects for the hydrated proton complex are effectively included in the MS-EVB framework to some extent1 and also that the dominant electronic polarization effects in water solvent can in principle be incorporated into a pairwise decomposable empirical water potential.60 A comparison between the density profiles of the CEC and O* in Figure 1 provides further insight into the hydrated proton structure at the water-vapor interface. Because the CEC is given by the expectation value of the center of excess charges (eq 1), it can be different from the O* (most hydronium-like oxygen) position, which is in fact observed in Figure 1. At the water-vapor interface, the CEC density profile at the interface is shifted to the bulk phase with respect to that of O* by ∼0.75 Å. This implies that the hydronium cation tends to orient at the interfaces so that the three O-H sides of the ion are toward the bulk phase and the excess proton is delocalized in that direction. This is consistent with an earlier simulation result by Petersen et al. for the water-vapor interface that the lone pair side of the hydronium ion is oriented toward the vapor region4 as well as subsequent simulations using a classical hydronium model.8,9 Interestingly, the ability of the excess proton charge defect to delocalize into the bulk region has the effect of enhancing the effective hydronium oxygen density at the interface over the classical hydronium model (cf. the solid versus dashed purple curves in Figure 1). It is estimated from the magnitudes of the density maxima that the nonclassical charge delocalization effect accounts for about 25% of the hydronium concentration enhancement at the interface.

The solvation structure around the hydronium ion was investigated by dividing the slab systems into several layers with 1.5 Å thickness from |zO*| ) 7.5 Å. The probability distributions of the angle φ between the three O-H bonds of the pivot hydronium ion and surface normal are shown in Figure 2. The surface normal is taken as a vector in the z direction from the bulk to vapor region, and thus, a negative value of cos φ means that the hydrogen atom in the O-H bond is toward the bulk region. It is noted that a very similar analysis has been done based on classical hydronium simulations.8,9 Consistent with the insight from the density profile mentioned above, Figure 2 indicates that most of the O-H bonds are oriented toward the bulk region and the lone pair side is toward the vapor region. A plausible reason is that such an orientation allows H3O+ to form HBs with water molecules without affecting the water HB network due to the hydrophobic part of H3O+.4 Nearly uniform distributions are observed in the layer from |zO*| ) 7.5 Å to |zO*| ) 9.0 Å, indicating that the specific orientation of the pivot hydronium ion already disappears by this region away from the interface. It is noted that similar features were observed in the results from the classical hydronium model (data not shown). The probability distributions of the same angle for the water molecules near the pivot hydronium ion were also investigated. To calculate these properties, water molecules around O* within an O-O* distance of 2.8 Å were sampled. Typically, three water molecules were taken into account by this procedure. Then the angle between the O-H bonds in such waters and the surface normal was calculated. The resultant angle distributions in the topmost interface layer from |zO*| ) 13.5 Å to |zO*| ) 15.0 Å are shown in Figure 3. For comparison, the same distributions calculated from all the water contributions inside the same layer are also displayed. As seen in Figure 3a, the distribution P(cos φ) calculated for all the water contributions has a peak around cos φ ) 1, indicating that some of the O-H bonds are directed toward the vapor region. These dangling bonds are consistent with experimental observation30 and with other simulation results.61 By contrast, there are no distinctive peaks in the positive cos φ region for the waters around the pivot hydronium, indicating that most of the O-H bonds in the waters near the hydronium tend to be directed toward the bulk region. These results show that dangling O-H bonds are suppressed around the hydronium ion. Indeed, Petersen et al.4 showed a negative correlation between the hydronium ion and the dangling O-H bonds at the water-vapor interface by analyzing a radial distribution function obtained from the MS-EVB2 simulations. They argued that such a negative correlation is related to the experimental finding that there is a decrease in the SFG spectral intensity


Figure 3. Probability distributions of two angles for the water-vapor interface in the layer |zO*| ) 13.5-15.0 Å. (a) Distribution of angle φ between the surface normal and O-H bond of waters around the pivot hydronium ion. (b) Distribution of angle θ between the surface normal and dipole moment of the same waters. Solid and dashed curves represent the results from the MS-EVB3 simulation and classical hydronium simulation, respectively. For comparison, distributions including all the water contributions obtained from the MS-EVB3 simulation in the same layer (dashed-dotted line) and in the inner layer, |zO*| ) 12.0-13.5 Å (dotted line), are shown.

assigned to a dangling O-H peak with an increase in acid concentration.15,62,63 The current result is consistent with their finding. Further insight can be obtained from the probability distributions of the angle θ between the water permanent dipole moment and the surface normal. As seen in Figure 3b, the distribution P(cos θ) for all the water contributions at the topmost layer has a peak around cos θ ) 0.25, indicating that the water dipole moment tends to direct slightly outward to the bulk phase, consistent with other simulation results.61 Such a tendency becomes less distinctive in the inner interface layer from |zO*| ) 12.0 Å to |zO*| ) 13.5 Å. By contrast, the distribution for waters around the pivot hydronium ion shows a peak around cos θ of -0.75, indicating that the dipole moments of such waters have the opposite tendency. Because of the preferred orientation of the hydronium ion discussed above, the dipole moment of the hydronium ion tends to orient toward the bulk region. It is thus expected that such inverted water dipoles around the hydronium ion stabilize the hydrated proton structure. It should be noted that the insights into the protonated solvation structures presented here are similar to those in the bulk phase discussed in ref 29 in terms of the amphiphilic nature of the hydrated proton. As discussed later, water molecules near the pivot hydronium have non-negligible probabilities to form the hydronium ion due to the excess proton charge defect delocalization. Thus, the two O-H bonds in such water molecules may still have a tendency to orient toward the bulk, because the lone pair side of the hydronium ion tends to be toward the vapor region. This perspective can also explain the features of the distributions discussed above. To analyze this in more detail, the same distributions were calculated from simulations with the classical hydronium ion model, which are included in Figure 3. Although there are slight quantitative differences, the main features in

Iuchi et al.

Figure 4. (a) Probability distributions of the highest and second highest EVB amplitudes for the water-vapor interface in the layers |zCEC| ) 13.5-15.0, 12.0-13.5, 10.5-12.0, and 0.0-7.5 Å. (b) Free energy profiles with respect to the highest EVB amplitude.

the distributions are the same. Therefore, the water orientations around the hydronium ion seem to mainly reflect the inherent hydrated proton structure due to the amphiphilic nature of the hydronium ion. In fact, the experimental SFG spectra show the increase of the peak assigned to the HB structure with respect to the change from neat water to HCl solution.14,15,64 It was argued from an MD simulation result that the orientational ordering of water molecules within the solvation shell of the hydronium ion is a source of this increase.65 The present results seem to be related to this argument. The delocalization of the excess proton charge defect around the pivot hydronium ion was also investigated. As discussed above, the difference of maximum positions in the CEC and of the O* density profiles implies that the excess proton is delocalized to some degree through the HB network at the interface. Here, the degree of delocalization is investigated in terms of the distributions of the largest and second largest EVB amplitudes, P(c12) and P(c22), in several interface water layers: |zCEC| ) 13.5-15.0, 12.0-13.5, 10.5-12.0, and 0.0-7.5 Å, where z ) 0 is the COM of the system. The resultant distributions are shown in Figure 4, where the distribution in the layer from 0.0 to 7.5 Å was found to be very similar to the distribution obtained from a 1 ns bulk water MS-EVB3 simulation. As seen in Figure 4, while the Eigen-like form is the dominant species at the interface regions, which is manifested in a peak around c12 ≈ 0.6, the probability of the Zundellike form around c12 ≈ c22 ≈ 0.45 is not negligible. The present result thus demonstrates that the excess proton clearly delocalizes through the water HB networks at the water-vapor interface and the excess proton samples a continuous distribution of species between the Eigen-like and Zundel-like limiting forms. Compared to the distribution in the bulk region, the distributions around 0.45 become larger for both c12 and c22, indicating that the probability of the Zundel-like form becomes larger at the interface regions than in the bulk. This insight can also be seen in the free energy profiles ∆F(c12) obtained from P(c12). As shown in Figure 4b, the relative free energy difference between the Zundel-like and the Eigen-like forms becomes smaller in the interface region. It should be noted that Zundel-like forms at the surface of protonated water clusters have been observed



TABLE 2: Lateral Diffusion Constantsa (Å2/ps) CECb vapor

wall

|z| region, Å

EVB

classical

12.0-15.0 9.0-12.0

0.34(0.09)c 0.34(0.06)c

0.16(0.05) 0.16(0.04)

EVB 0.30(0.08)

CCl4 classical

EVB

classical

0.18(0.08)

0.26(0.08)c 0.31(0.07)c

0.13(0.04) 0.14(0.02)

EVB d

0.30,

bulk

classical

0.29e

0.13 Waterf

vapor

wall

CCl4

|z| region, Å

EVB

classical

EVB

classical

EVB

classical

12.0-15.0 9.0-12.0

0.55(0.05) 0.30(0.02)

0.52(0.05) 0.29(0.02)

0.26(0.01)

0.28(0.01)

0.29(0.02) 0.24(0.02)

0.30(0.04) 0.24(0.02)

d

bulk

EVB

classical

0.23

0.23,0.232g

a Values in parentheses represent the standard deviations determined from 10 trajectories. b OT is used in the classical hydronium cases. c To increase statistics, the CEC which goes outside the layer within 2 ps was taken into account (cf. ref 61). d From 250 ps long constant NVE simulations for the bulk systems. EVB system: 500 waters, 1 Cl-, and 1 H+. Classical system: 499 waters, 1 Cl-, and 1 classical H3O+. e Reference 29. f Contributions from water molecules which participate in the hydrated proton complexes were not taken into account for the water diffusion constants. g Reference 39.

in ab initio MD simulations7 and in an MS-EVB simulation18 using the MS-EVB1 model,32 where the slightly enhanced probability of Zundel-like forms was observed in the latter simulation. The existence of both the Eigen and Zundel cations at the water-vapor interface is also suggested by experiments.14,15 3.2. Dynamical Properties. It is also interesting to investigate dynamical properties of the hydrated excess proton at the water-vapor interface. Here, the lateral diffusion constant was calculated on the basis of the mean square displacements35,61 by dividing the slab system into several layers. To be specific, the lateral diffusion constants of water and the excess proton CEC (or OT in the classical hydronium model) were evaluated in two layers, |z| ) 9.0-12.0 Å (region 1) and |z| ) 12.0-15.0 Å (region 2), which are summarized in Table 2. Compared to the MS-EVB3 water diffusion constant of 0.23 Å2/ps in the bulk phase, the lateral diffusion constant of water becomes larger because of the lower water density near the vapor: 0.30 and 0.55 Å2/ps in regions 1 and 2, respectively. This observation is consistent with other simulation results in the pure water-vapor interface system.61 On the other hand, the increase of the CEC diffusion constant from the bulk to interface regions is not large: 0.30, 0.34, and 0.34 Å2/ps in the bulk, region 1, and region 2, respectively. It is interesting to note that while the excess proton CEC does not diffuse faster at the interface than in the bulk in an absolute sense, the ratio of the CEC lateral diffusion constant to the underlying lateral water diffusion constant decreases by a factor of ∼2. This result suggests that the excess proton is also dynamically stabilized at the water-vapor interface. A comparison of the diffusion constants between the EVB and classical hydronium models highlights the role of Grotthuss effects of the excess proton. As tabulated in Table 2, the classical hydronium model gives significantly smaller diffusion constants of OT in all the regions: 0.13, 0.16, and 0.16 Å2/ps in the bulk, region 1, and region 2, respectively. 3.3. Multiple Proton Simulation. Throughout this study, the excess proton behavior is investigated primarily in systems containing one excess proton. On the other hand, some classical MD simulations have been carried out for the water-vapor interface using classical hydronium models with a concentration

Figure 5. Normalized density profiles of the hydrated excess proton CEC, O*, OW, and Cl- for the 1.1 M HCl solution liquid-vapor interface obtained from SCI-MS-EVB3 simulations. For comparison, results of the single-proton case corresponding to the 0.1 M HCl solution are also displayed (dashed lines).

of HCl at ∼1.2 M to compare to experimental SFG spectra.8,9 To examine the effects of higher proton concentrations, density profiles for the 1.1 M HCl solution liquid-vapor interface using the multiproton SCI-MS-EVB method34 are shown in Figure 5 and compared with results for the single-proton case corresponding to the 0.1 M HCl solution. Figure 5 clearly shows the surface preference of the excess proton at this concentration, though the probability of the CEC at the interface relative to the bulk is smaller than that in the single-proton case. This may be because hydronium ions at the interfaces prevent other hydronium ions from going to the interface. Indeed, in the 10 × 120 ps length trajectories, typically a few excess protons were seen in the bulk phase. In addition to the decrease of the CEC density maximum, it is seen in Figure 5 that the chloride anions have a small maximum in the density near the maximum in the CEC distribution. This double-layer structure is thought to be formed due to the hydronium ions at the interface.9 The lack of explicit polarization effects of the chloride anions may underestimate this trend, which may further affect the surface preference of the excess proton. In fact, the double-layer structure seems to be less pronounced compared to that in a classical hydronium MD simulation with a polarizable model.9


Iuchi et al.

Figure 6. Normalized density profiles of OT, HT, OW, and Cl- for the water-vapor interface obtained from quantum PIMD simulations (solid lines) and classical MD simulations (dashed lines). The classical hydronium model is used in both simulations.

3.4. Nuclear Quantum Effects. It is known that some properties of the excess proton such as the diffusion constant are affected by nuclear quantum effects.29 Therefore, the explicit inclusion of nuclear quantum effects might also affect the interface structures and the surface preference of the excess proton. In this section, therefore, nuclear quantum effects on the surface preference of the excess proton are briefly addressed using a PIMD simulation with the classical hydronium model. PIMD simulations with the full MS-EVB3 model were too computationally demanding given the available computational resources for this project. The density profiles obtained from the PIMD simulations are shown in Figure 6, where those obtained from the classical simulations are included for comparison. The surface preference of the hydronium ion is clearly observed, and there is no distinctive difference compared with the classical result. In addition, the shift between OT and HT distributions at the interface region due to the preferred orientation of the hydronium ion is observed in both the PIMD and classical simulations. As discussed above, the excess proton still forms a distinct hydrated proton structure at the interface, and thus, it seems to be reasonable to conclude that the explicit inclusion of nuclear quantum effects does not affect the surface preference to a large degree. It should be noted, however, that the current results do not imply that nuclear quantum effects on other properties such as diffusion and spectroscopy may not be important. 3.5. Free Energy Profiles. It has been demonstrated that the amphiphilic nature of the hydrated proton is manifested in the anisotropic hydronium ion orientation and solvation structures at the water-vapor interface. Such an anisotropic hydration structure allows the delocalization of excess proton charge defect at the interface through the water HB network, which may provide additional stabilization of the proton solvation structure. To fully understand this behavior, since it is necessary to analyze the free energy profile and its contributions consistently, it is not clear whether the entropic contribution is negligible. For example, it was discussed that the entropic contributions are non-negligible for Cl- surface solvation in a water cluster.66 It is also possible that accommodating the hydronium ion in the bulk region is entropically unfavorable due to the hydrophobic character in the lone pair region,7 and thus, an entropic contribution might lead to the surface preference of the excess proton in the traditional sense of a hydrophobic effect. A systematic analysis of the free energy profiles thus provides critical insight into the thermodynamic driving forces for the excess proton surface preference. Figure 7 shows the free energy profiles of the excess proton charge defect CEC at different temperatures. At 298.15 K, the free energy of the CEC in the interface region is lower than

Figure 7. (a) Hydrated excess proton CEC free energy profiles for the water-vapor interface at different temperatures obtained from umbrella sampling simulations. The dashed line represents the free energy profile of classical hydronium oxygen at 298.15 K obtained from the density profile in Figure 1. (b) Decompositions of the CEC free energy profiles. Red and green lines correspond to entropic and internal energy contributions, respectively.

that of the CEC in the bulk region by ∼1.8 kcal/mol, which is about 3 times greater than the thermal fluctuation, kBT ≈ 0.6 kcal/mol. This value can be compared to a value of ∼2.4 kcal/ mol obtained from the MS-EVB1 model32 at 280 K for the water droplet system.18 The delocalization effect of the excess proton is seen by comparing the MS-EVB3 result with the free energy profile of the classical hydronium ion model, where a smaller free energy minimum of ∼0.9 kcal/mol is observed. This result suggests that the excess proton charge delocalization nearly doubles the free energy of surface stabilization and that classical hydronium models are missing this key physics.8-12 Further insight can be obtained from the temperature dependence of the CEC free energy profiles. As seen in Figure 7a, the free energy profiles indicate that the stabilization in the free energy from the bulk to the interface decreases with respect to an increase of temperature. As a result, the entropic contribution, -T∆S(z), increases from the bulk to the interface, corresponding to a decrease of the entropy, as seen in Figure 7b. This entropy decrease suggests a more ordered protonated water structure at the interface than in the bulk, consistent with the preferred orientations of the interface hydrated proton structure as discussed earlier. It should be noted that the CEC diffusion constant at the water-vapor interface was not very different from that in the bulk. Therefore, the lateral diffusive motion may not affect the change of the entropic contribution very much.48 The ordered solvation structures, however, are energetically favorable since the hydronium cation can still form HBs with water molecules at the interface without perturbing the water HB network by the lone pair side4 as discussed in the previous section. Such a favorable change of the internal energy term, ∆E(z), leads to an overall surface stabilization seen in Figure 7. As seen in Figure 7b, the entropy term, -T∆S(z), decreases beyond ∼12 Å, though the hydrated proton structure is more “rigid” as seen in Figures 2 and 3. This is interpreted as an increase of the number of accessible translational states as a result of fewer solvent molecules near the vapor region. Such



Figure 9. Normalized density profiles of hydrated excess proton CEC, OW, Cl-, and O* for the water-hydrophobic wall interface. Dashed and dotted lines represent results from the classical hydronium simulation and from the MS-EVB3 simulation including H-wall interactions, respectively. The origin of z ) 0 corresponds to the middle point between two walls.

Figure 8. (a) Normalized density profiles of the water oxygen for the water-hydrophobic wall interface (solid line), the water-hydrophobic wall including H-wall interactions (dashed line), and the water-vapor interface (dashed-dotted line). (b) Probability distributions of angle φ between the surface normal and O-H bond and of angle θ between the surface normal and water dipole moment in the region from |zOW| ) 12.0 Å to |zOW| ) 13.5 Å for the water-hydrophobic wall interface (solid line) and for the same interface including H-wall interactions (dashed line).

configurations may, however, become energetically unfavorable since the stabilization of the hydronium ion and/or water HB network around the excess proton is likely smaller in a lower water density region, as seen in the free energy profile. A classical hydronium polarizable model simulation demonstrated that adding explicit polarizability into the water and classical hydronium ion models enhances its surface preference for the water-vapor interface.8 The favorable internal energy change presented here seems to be consistent with this argument, since many-body polarization effects would stabilize the surface anisotropic protonated water structure. The present argument is also consistent with the behavior of cations at the water-vapor interface. It has been shown from MD simulations that small cations are repelled from the surface layer of the water-vapor interface due to their smaller polarizabilities.8,40 This is in contrast to the case of the excess proton which forms an energetically stable protonated solvation structure at the interface due in part to the amphiphilic nature of the hydronium ion and in part to the highly effective polarization of the hydrated excess proton structures. 4. Results: Water-Hydrophobic Wall Interface As mentioned in the Introduction, the behavior of the excess proton at water-hydrophobic interfaces is of fundamental importance for understanding proton solvation and dynamics in various systems such as membranes and proteins. In this section, the hydrated proton behavior at a water-hydrophobic wall interface is analyzed and compared to the case of the water-vapor interface. 4.1. Interface Excess Proton Solvation Structures. Figure 8a shows the density profiles of water molecules for the water-hydrophobic wall interface, where the origin of z ) 0 corresponds to the middle point between the two walls. As seen in Figure 8a, the oscillations around the solid wall are distinctive,

whereas such oscillations are not seen in the profile for the water-vapor interface. This is because the solid wall forces the waters to form a quasi-layered structure.43,67 The HW-wall interaction mentioned in section 2.3 gives additional attractive/ repulsive interactions from the wall, resulting in more pronounced oscillations as seen in Figure 8a. The water orientations at these interfaces were investigated in terms of the angle distributions, P(cos φ) and P(cos θ), where φ and θ are angles between the surface normal and the O-H bond and between the surface normal and the water dipole moment, respectively. The resulting distributions at the interface region from |zOW| ) 12.0 Å to |zOW| ) 13.5 Å are shown in Figure 8b. A comparison between the results of the water-vapor (Figure 3) and water-hydrophobic wall interfaces reveals that the main features of the distributions are very similar. In particular, both the P(cos φ) distributions have peaks around cos φ ) 1, corresponding to the O-H dangling bonds. This is consistent with other wall simulations43,45 and with experimental evidence that both the water-vapor and water-solid hydrophobic interfaces are characterized by the dangling O-H bonds.30 As seen in Figure 8b, adding an HW-wall interaction leads to a different feature in the angle distributions. Because of the interaction between the HW and wall, the hydrogen atoms cannot penetrate into the wall side, and consequently, the peak around cos φ ) 1 disappears in the interface region. Consistently, the dipole direction is inverted compared to the result without the HW-wall interaction. This seems to be inconsistent with the experimental observation mentioned above, and the HW-wall interaction employed here is probably too large. However, the surface preference of the excess proton, which will be discussed below, was not affected significantly by this factor. Therefore, the simulation results without the hydrogen-wall interactions will be mainly discussed for the water-hydrophobic wall interface system. The normalized z-dependent density profiles of the hydrated excess proton for the water-hydrophobic wall interface are shown in Figure 9. As in the case of the water-vapor interface, the surface preference of the excess proton and the difference between the CEC and O* density profiles at the interface are observed. It is noted that the preference of the proton for the hydrophobic surface was also observed in an ab initio MD simulation.28 Due to the solid wall, the CEC does not penetrate into the hydrophobic region, resulting in a sharper CEC density profile than that of the vapor interface. It is also observed that the CEC density profile obtained from the simulation including the HW-wall and HT-wall interactions is sharper than that


Figure 10. Normalized probability distributions of two angles for the water-hydrophobic interface in the layers |zO*| ) 12.0-13.5 Å (solid line) and 10.5-12.0 Å (dashed line). (a) Distribution of angle φ between the surface normal and O-H bond. (b) Distribution of angle θ between the surface normal and dipole moment. Blue, red, and green lines represent the results of waters near the pivot, of all waters, and of the pivot hydronium ion, respectively.

without these interactions. However, the difference in the two CEC density profiles is not very large: the difference in the corresponding free energy stabilizations from the bulk to the interface regions was estimated to be ∼0.2 kcal/mol. The hydrated proton structures were investigated in terms of the same angle probability distributions for the water-vapor interface discussed in section 3.1. The resulting distributions at the surface layers from |zO*| ) 12.0 Å to |zO*| ) 13.5 Å and from |zO*| ) 10.5 Å to |zO*| ) 12.0 Å are shown in Figure 10. As in the case of the water-vapor interface, most of the O-H bonds of the pivot hydronium ion are oriented toward the bulk region and the lone pair side is toward the wall, which is especially prominent at the topmost layer. The main features of the probability distributions for the water molecules near the pivot hydronium ion are also similar to those of the water-vapor interface. These results suggest that the overall feature of the proton solvation is governed by the amphiphilic nature of the hydronium ion and is not sensitive to the differences in the two interfaces. Although the main features are similar, the distribution of the hydronium O-H bonds is narrower than that of the water-vapor interface, implying a slightly more Eigen-like ordered solvation structure due to the solid wall. This is manifested in a slight change of delocalization as seen in Figure 11. At the topmost surface layer, from |zCEC| ) 11.3 Å to |zCEC| ) 12.3 Å, the population of the Zundel-like form relative to the Eigen-like form is smaller than that of the topmost layer at the water-vapor interface (Figure 4). From the experimental SFG spectra of pure water-hydrophobic interfaces (without the excess proton), it has been argued that the solid wall does not allow flexibility and the spectrum of the pure water-solid hydrophobic interface shows a more icylike feature than that of the water-vapor interface.30,63,68 Similarly, the present result implies that the solid wall enhances ordering of the hydrated proton structures. It should be noted again, however, that the

Iuchi et al.

Figure 11. (a) Probability distributions of the highest and second highest EVB amplitudes for the water-hydrophobic wall interface in the layers |zCEC| ) 11.3-12.3, 10.3-11.3, 9.3-10.3, and 0.0-8.3 Å. (b) Free energy profiles with respect to the highest EVB amplitude.

degree of this effect is not large and the main features of the solvation structures are similar to those of the water-vapor interface. 4.2. Dynamical Properties. As tabulated in Table 2, the lateral diffusion constant in the layer from |zCEC| ) 9.0 Å to |zCEC| ) 12.0 Å was calculated to be 0.30 Å2/ps, which is the same as in the bulk phase. This is similar to the water-vapor interface case as discussed in section 3.2, where relatively little change in the diffusion constants was observed. This result is, for example, in contrast to that for a phospholipid membrane surface, where the proton diffusion constant becomes smaller than in the bulk due to the polar group region of the membrane.69 4.3. Free Energy Profiles. Compared to the free energy profile for the water-vapor interface (Figure 7), the Helmholz free energy profile, ∆F(z), for the water-hydrophobic wall interface has a narrower minimum region as seen in Figure 12a. The stabilization in the CEC free energy from the bulk to the interface regions is ∼2.2 kcal/mol, which is slightly larger than that of the water-vapor interface, ∼1.8 kcal/mol. To discuss differences between the free energy profiles for these two interface systems, the contribution of wall interactions was eliminated from the O* free energy profile by a procedure described in section 2.5. The resulting O* free energy profile without the wall contribution displayed in Figure 12b shows a further decrease of the free energy beyond ∼12.0 Å, indicating that the sharper and narrower profile for the water-hydrophobic wall interface is due to the solid wall. In addition, the free energy difference between the bulk and interface regions becomes closer to that of the water-vapor interface when the wall contribution is not taken into account. In these respects, it is reasonable to conclude that the difference between the free energy profiles at the water-vapor and the water-hydrophobic interfaces is relatively small. A small hydrophobic molecule is solvated without disrupting the surrounding water HB network very much, but enhancing the order of the surrounding water molecules.25,67 From this viewpoint, it is possible that the hydronium ion in bulk water may be entropically unfavorable,4,7 and thus, the driving force for the surface preference of the excess proton near a hydro-



Figure 13. Normalized density profiles of hydrated excess proton CEC, OW, Cl-, O*, and CCl4 carbon (CT) and Cl for the water-CCl4 interface. Dashed lines represent results from the classical hydronium model. The origin of z ) 0 corresponds to the center of the water density profile.

Figure 12. (a) Hydrated excess proton CEC free energy profiles for the water-hydrophobic wall interface at different temperatures obtained from umbrella sampling simulations. The dashed line represents the free energy profile of classical hydronium oxygen at 298.15 K obtained from the density profile in Figure 9. Error bars are only shown in the profile at 298.15 K, and their magnitudes at other temperatures are similar. (b) O* free energy profiles with and without wall-O* interaction contribution (blue and red lines, respectively). (c) Decompositions of hydrated excess proton CEC free energy profiles. Red and green lines correspond to entropic and internal energy contributions, respectively.

phobic wall might have a favorable entropic contribution. This argument is similar to the fact that a contact between two hydrophobes leads to a favorable entropic contribution as entropically unfavorable water molecules are excluded from the region between two hydrophobes.67 However, the resultant behaviors of ∆E(z) and -T∆S(z) for the water-hydrophobic wall interface in Figure 12c again indicate that this is not the case. The behaviors of two contributions are qualitatively similar to those for the water-vapor interface: the entropic contribution, -T∆S(z), increases from the bulk to the interface regions, while the internal energy decreases, resulting in the overall surface preference, even though the individual magnitudes of these free energy terms in the water-wall case are substantially smaller than in the water-vapor case (Figure 7). As discussed in the water-vapor interface case, the unfavorable change of the entropic contribution is explained by the preferred hydrated proton structure at the interface. The rigid proton solvation structure at the interface is energetically favorable, overcoming such an entropic loss. The smaller increase of -T∆S(z) than in the water-vapor interface may be rationalized considering that the solid wall does not allow water and hydronium ion molecules to move to a large degree near the wall region. As a consequence, the space near the wall becomes a more moleculefree space and the entropy contribution rapidly decreases in that space as seen in Figure 12c. Considering the present results together with those for the water-vapor interface in section 3.5, the excess proton plays a role as a “structure maker” in the sense that a rigid hydrated proton structure is formed at the interfaces with an entropy reduction (i.e., an unfavorable entropic contribution).66 Such a rigid hydrated proton structure leads to a favorable energy change, resulting in the overall surface preference.

The present results also seem to be consistent with the other interface simulation results in some ways. First, the simulation results for a liquid water system between two flat walls demonstrated that the water molecules at the interface are rearranged and have dangling O-H bonds so as to maximize the total interaction energy.43 Similarly, the excess proton is energetically stabilized by forming a preferred hydration structure due to its amphiphilic nature, which overcomes an associated entropic penalty as mentioned above. Second, the simulation result by Wallqvist and Berne31 demonstrated that a solvent-separated pair between a methane molecule and hydrophobic wall is more plausible than a contact pair between them, suggesting that a contact between a small hydrophobic molecule and flat hydrophobic wall is not always achieved. This seems to be consistent with the present insight that the surface preference of the excess proton does not originate from the favorable entropic contribution, but from the energetically stabilized surface hydration structure of the hydrated excess proton. 5. Results: Water-CCl4 Interface In this section, analyses similar to those in sections 3 and 4 are carried out for the water-hydrophobic liquid-liquid interface. To be specific, the water-CCl4 interface system is chosen because this system has been studied experimentally.70,71 5.1. Interface Excess Proton Solvations and Dynamics. Figure 13 shows the z-dependent density profiles for OW, CEC, Cl-, and O*, where the positive and negative z regions were assumed to be equivalent and averaged with respect to z ) 0 (the middle point of the water density profile). As seen in Figure 13, the density profiles for the water-CCl4 interface are very similar to those for the water-vapor interface (Figure 1). Experimental results have suggested that the excess proton exists at the water-vapor interface on the basis of the SFG spectral change from pure water to acidic solutions.8,14,15 On the other hand, the recent SFG experiments for the water-CCl4 interface showed little spectral change in the free water region from low pH (HCl solution) to neutral conditions.70,71 The connection of these experimental results to the preference of the excess proton for the water-CCl4 interface obtained from this work is therefore not completely clear. The same angle distributions as in the previous sections for the water-vapor and water-hydrophobic wall interfaces were also calculated for the water-CCl4 interface. The resultant distributions at the surface layer from |zO*| ) 13.5 Å to |zO*| ) 15.0 Å are displayed in Figure 14. Compared to the distributions for the water-vapor (Figure 3) and for the water-hydrophobic


Figure 14. Probability distributions of two angles for the water-CCl4 interface in the layer |zO*| ) 13.0-15.0 Å. (a) Distributions of angle φ between the surface normal and O-H bonds. (b) Distributions of angle θ between the surface normal and dipole moment. Solid and dashed lines represent the results of the waters near the pivot and of the pivot hydronium ion, respectively. For comparison, distributions including all the water contributions in the same layer (dashed-dotted line) and in the inner layer |zO*| ) 12.0-13.5 Å (dotted line) are shown.

wall (Figure 10) interfaces, the main features of the distributions for the water-CCl4 interface are closer to those of the distributions for the water-vapor interface. This seems to be reasonable, since the water-CCl4 interface is as rough as the water-vapor interface, which is manifested, for example, in the water density profile. The similar EVB amplitude distributions (data not shown here) also support this viewpoint. The lateral diffusion constants were also evaluated for the water-CCl4 interface system and are tabulated in Table 2. The CEC diffusion constants at the interface layers are slightly smaller than that in the bulk, presumably because of the mixing of CCl4 molecules. This trend is opposite that of the other two interface cases. 5.2. Free Energy Profiles. Figure 15a shows the free energy profiles of the CEC at different temperatures for the water-CCl4 interface. At 298.15 K, the CEC free energy stabilization at the interface region is ∼1.6 kcal/mol, which is similar to ∼1.8 kcal/ mol for the water-vapor interface. As seen in the CEC density profile from the unbiased simulations (Figure 13), the CEC in the bulk phase was less probable compared to that of the water-vapor interface case (Figure 1), though the unbiased simulation lengths were the same. In this regard, a slight shoulder of ∼0.3 kcal/mol between the bulk and interface regions appeared in the CEC free energy profile and may not be due to a numerical error. The resultant behaviors of ∆E(z) and -T∆S(z) in Figure 15b show that the internal energy contribution is again responsible for the overall surface preference at z ≈ 12 Å, which is the same as in the other interface cases. Such a favorable internal energy relative to that of the bulk region is again explained by the amphiphilic nature of the hydrated proton: the preferred orientation of the hydrated proton structure can be energetically stabilized. However, at the intermediate region between the bulk and interface, the behaviors of ∆E(z) and -T∆S(z) are different from those of the other two interfaces. Because the density of the CCl4 molecule starts increasing at z ≈ 10 Å as seen in Figure

Iuchi et al.

Figure 15. (a) Hydrated excess proton CEC free energy profiles for the water-CCl4 interface at different temperatures obtained from umbrella sampling simulations. The dashed line represents the free energy profile of the classical hydronium oxygen at 298.15 K obtained from the density profile in Figure 13. Error bars are only shown in the profile at 298.15 K, and their magnitudes at other temperatures are similar. (b) Decompositions of the hydrated excess proton CEC free energy profiles. Red and green lines correspond to entropic and internal energy contributions, respectively.

13, it is likely that the solvation environment around the excess proton becomes more inhomogeneous and more perturbed by the CCl4 molecules when the excess proton moves from the bulk region to such an intermediate region. Therefore, the decrease of -T∆S(z) from z ≈ 6 Å to z ≈ 9 Å may be interpreted as the entropy increase associated with such an environmental change. In this intermediate region, the internal energy is less favorable than that in the bulk possibly due to the onset of mixing with the CCl4 molecules. 6. Conclusions In this study, the behavior of the excess proton was extensively investigated through the MS-EVB3 simulations for three different types of water-hydrophobic interfaces: the water-vapor, water-hydrophobic wall, and water-CCl4 interfaces. It was seen that the hydrated excess proton exhibits a preference for all of the interfaces, with the lone pair side of the hydronium ion being directed toward the hydrophobic regions. The orientations of the water molecules around the hydronium ion were also found to be restricted, indicating rigid hydrated proton structures at the interfaces. It was found that the excess proton charge defect delocalizes at all the interfaces with a slight probability increase of the Zundel-like form compared to that of the bulk region. Such a proton delocalization enhances the surface stabilization of the hydrated proton structures over the classical hydronium limit. The thermodynamic driving force for the surface preference of the excess proton was systematically investigated by decomposing the free energy profiles into the internal energy and entropic contributions. As a result, it was found that the overall surface preference is governed by the favorable internal energy change at all of the interfaces investigated here. Because the hydrated proton structures at the interfaces do not disturb the water hydrogen bond networks very much, the hydrated proton



structures can be energetically stabilized. The surface stabilization overcomes the unfavorable entropic contributions associated with the entropy-lowering structural arrangements of water and hydronium ion molecules. The present results may also provide useful insight into the hydrated excess proton behavior in other important mixed dielectric environments such as membranes and proteins. Acknowledgment. This research was supported by a grant from the National Science Foundation (CHE-0719522). We thank Prof. Feng Wang and Dr. Mark Maupin for valuable discussions and comments.

The internal energy change at z ) zi with respect to the reference point z ) z0 is given by

)

∂(-β∆A(zi)) ∂(-β)

∫ δ(z - zi)Ue-βU dΓ ∫ δ(z - z0)Ue-βU dΓ ∫ δ(z - zi)e-βU dΓ ∫ δ(z - z0)e-βU dΓ

) 〈U(zi)〉zi - 〈U(z0)〉z0

(7)

where Γ and U represent all the coordinates and the potential energy, respectively. The wall contribution in the internal energy change is then given by

∫ δ(z - zi)Uwall(z)e-βU dΓ ∆Ewall(zi) ) ∫ δ(z - zi)e-βU dΓ ∫ δ(z - z0)Uwall(z)e-βU dΓ ∫ δ(z - z0)e-βU dΓ

(8)

In the usual classical simulations, eq 8 becomes ∆Ewall(zi) ) Uwall(zi), with an assumption of Uwall(z0) ) 0, because the wall potential Uwall(z) only depends on z. Here, the contribution from the interaction between the wall and O* is subtracted from the O* free energy profile obtained from the MS-EVB simulation. Each EVB state feels a different wall potential, and thus, eq 8 becomes somewhat complicated as

∫ δ(z - zi)c12(Γ)e-βU dΓ + ∫ δ(z - zi)e-βU dΓ N ∫ δ(z - zi)cI2(Γ)e-βU dΓ wall ∑ UI (zi) I)2 ∫ δ(z - zi)e-βU dΓ

∆Ewall(zi) ) Uwall 1 (zi) EVB

NEVB

≡ Uwall 1 (zi)R1(zi) +

∑ Uwall I (zi) RI(zi)

wall ∆Ewall(zi) ) UOT (zi) R1(zi) + Uwall OW (zi){1 - R1(zi)} (10)

Because the difference between the OT-wall and OW-wall interactions is not large due to the similar LJ parameters for OT and OW, it was found that the wall contribution ∆Ewall(zi) is not sensitive to R1(zi). As a good approximation, therefore, R1 ) 0.6 was simply used for all zi, because the distributions of c12 have a peak around this value in all slab regions as seen in Figure 11. It is certainly possible to calculate R1(zi) rigorously if necessary. References and Notes

Appendix: Eliminating the Wall Contribution from the Free Energy Profile

∆E(zi) )

OW-wall interactions:

(9)

I)2

where O* is an OT in EVB state 1 and an OW in the other NEVB RI ) 1 derived states (I * 1). Considering the condition ∑ I)1 from the normalization condition of the EVB coefficients, eq 9 is rewritten as a linear combination of the OT-wall and

(1) Voth, G. A. Acc. Chem. Res. 2006, 39, 143. (2) Swanson, J. M. J.; Maupin, C. M.; Chen, H.; Petersen, M. K.; Xu, J.; Wu, Y.; Voth, G. A. J. Phys. Chem. B 2007, 111, 4300. (3) Marx, D. ChemPhysChem. 2006, 7, 1848. (4) Petersen, M. K.; Iyengar, S. S.; Day, T. J. F.; Voth, G. A. J. Phys. Chem. B 2004, 108, 14804. (5) Petersen, M. K.; Voth, G. A. J. Phys. Chem. B 2006, 110, 7085. (6) Wang, F.; Izvekov, S.; Voth, G. A. J. Am. Chem. Soc. 2008, 130, 3120. (7) Iyengar, S. S.; Day, T. J. F.; Voth, G. A. Int. J. Mass Spectrom. 2005, 241, 197. (8) Mucha, M.; Frigato, T.; Levering, L. M.; Allen, H. C.; Tobias, D. J.; Dang, L. X.; Jungwirth, P. J. Phys. Chem. B 2005, 109, 7617. (9) Ishiyama, T.; Morita, A. J. Phys. Chem. A 2007, 111, 9277. (10) Buch, V.; Milet, A.; Vaćha, R.; Jungwirth, P.; Devlin, J. P. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 7342. (11) Vaćha, R.; Buch, V.; Milet, A.; Devlin, J. P.; Jungwirth, P. Phys. Chem. Chem. Phys. 2007, 9, 4736. (12) Wick, C. D.; Kuo, I.-F. W.; Mundy, C. J.; Dang, L. X. J. Chem. Theory Comput. 2007, 3, 2002. (13) Petersen, P. B.; Saykally, R. J. J. Phys. Chem. B 2005, 109, 7976. (14) Tarbuck, T.; Ota, S. T.; Richmond, G. L. J. Am. Chem. Soc. 2006, 128, 14519. (15) Levering, L. M.; Sierra-Hernańdez, M. R.; Allen, H. C. J. Phys. Chem. C 2007, 111, 8814. (16) Pegram, L. M.; Record, M. T. J. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 14278. (17) Burnham, C. J.; Petersen, M. K.; Day, T. J. F.; Iyengar, S. S.; Voth, G. A. J. Chem. Phys. 2006, 124, 024327. (18) Ko¨finger, J.; Dellago, C. J. Phys. Chem. B 2008, 112, 2349. (19) Beattie, J. K. Phys. Chem. Chem. Phys. 2008, 10, 330. (20) Vaćha, R.; Buch, V.; Milet, A.; Devlin, J. P.; Jungwirth, P. Phys. Chem. Chem. Phys. 2008, 10, 332. (21) Petersen, P. B.; Saykally, R. J. Chem. Phys. Lett. 2008, 458, 255. (22) It should be noted that the authors of refs 10 and 11 utilize an empirical model of a classical hydronium cation and classical hydroxide anion contained within a water “slab” geometry to mimic the neutral watervapor interface. In their simulations the hydronium and hydroxide do not have the possibility of chemically recombining into a water molecule as they do in nature (with a thermodynamic driving force of about 19 kcal/ mol, to be compared to their predicted ∼3 kcal/mol hydronium surface enhancement). Furthermore, the hydronium cation in their empirical simulation cannot participate in Grotthuss shuttling, and their hydroxide anion does not have the possibility of “inverse” Grotthuss shuttling. Simulation of the thermodynamics of a pure water system of pH 7 would require more than 500 million water molecules in addition to the physical features described above. The empirical simulations reported in refs 10 and 11 are therefore actually of an acidic solution, described by a classical hydronium cation and a classical hydroxide anion, at a concentration of 0.13 mol/L. (23) Shin, J.-W.; Hammer, N. I.; Diken, E. G.; Johnson, M. A.; Walters, R. S.; Jaeger, T. D.; Duncan, M. A.; Christie, R. A.; Jordan, K. D. Science 2004, 304, 1137. (24) Kusaka, I.; Oxtoby, D. W. J. Chem. Phys. 2000, 113, 10100. (25) Chandler, D. Nature 2005, 437, 640. (26) Pratt, L. R.; Pohorille, A. Chem. ReV. 2002, 102, 2671. (27) Dill, K. A. Biochemistry 1990, 29, 7133. (28) Kudin, K. N.; Car, R. J. Am. Chem. Soc. 2008, 130, 3915. (29) Wu, Y.; Chen, H.; Wang, F.; Paesani, F.; Voth, G. A. J. Phys. Chem. B 2008, 112, 467. (30) Du, Q.; Freysz, E.; Shen, Y. R. Science 1994, 264, 826. (31) Wallqvist, A.; Berne, B. J. Chem. Phys. Lett. 1988, 145, 26. (32) Schmitt, U. W.; Voth, G. A. J. Chem. Phys. 1999, 111, 9361. (33) Day, T. J. F.; Soudackov, A. V.; Cuma, M.; Schmitt, U. W.; Voth, G. A. J. Chem. Phys. 2002, 117, 5839. (34) Wang, F.; Voth, G. A. J. Chem. Phys. 2005, 122, 144105.

4030 J. Phys. Chem. B, Vol. 113, No. 13, 2009 (35) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1987. (36) Melchionna, S.; Ciccoti, G.; Holian, B. L. Mol. Phys. 1993, 78, 533. (37) Wu, Y.; Ilan, B.; Voth, G. A. Biophys. J. 2007, 92, 61. (38) Smith, W.; Forester, T. R. DL_POLY, version 2.13; CCLRC Daresbury Laboratory: Daresbury, Warrington, England, 1999. (39) Wu, Y.; Tepper, H. L.; Voth, G. A. J. Chem. Phys. 2006, 124, 024503. (40) Jungwirth, P.; Tobias, D. J. Chem. ReV. 2006, 106, 1259. (41) Berne, B. J.; Thirumalai, D. Annu. ReV. Phys. Chem. 1986, 37, 401. (42) Tuckerman, M. E. In Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms; Grotendorst, J., Marx, D., Muramatsu, A., Eds.; John von Neumann Institute for Computing: Ju¨lich, Germany, 2002; Vol. 10, p 269. (43) Lee, C. Y.; McCammon, J. A.; Rossky, P. J. J. Chem. Phys. 1984, 80, 4448. (44) Lee, C. Y.; Rossky, P. J. J. Chem. Phys. 1994, 100, 3334. (45) Shelley, J. C.; Patey, G. N. Mol. Phys. 1996, 88, 385. (46) Wallqvist, A. Chem. Phys. Lett. 1990, 165, 437. (47) Hill, T. L. An Introduction to Statistical Thermodynamics; Dover: Mineola, NY, 1986. (48) Zangi, R.; Engberts, J. B. F. N. J. Am. Chem. Soc. 2005, 127, 2272. (49) Chang, T.-M.; Dang, L. X. J. Chem. Phys. 1996, 104, 6772. (50) Schweighofer, K. J.; Essmann, U.; Berkowitz, M. J. Phys. Chem. B 1997, 101, 3793. (51) Tepper, H. L.; Voth, G. A. Biophys. J. 2005, 88, 3095. (52) Tepper, H. L.; Voth, G. A. J. Phys. Chem. B 2006, 110, 21327. (53) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. J. Comput. Chem. 1992, 13, 1011.

Iuchi et al. (54) Roux, B. Comput. Phys. Commun. 1995, 91, 275. (55) Smith, D. E.; Haymet, A. D. J. J. Chem. Phys. 1993, 98, 6445. (56) Wan, S.; Stote, R. H.; Karplus, M. J. Chem. Phys. 2004, 121, 9539. (57) Allen, T. W.; Andersen, O. S.; Roux, B. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 117. (58) Maupin, C. M.; Wong, K. F.; Soudackov, A. V.; Kim, S.; Voth, G. A. J. Phys. Chem. A 2006, 110, 631. (59) Liu, D.; Ma, G.; Levering, L. M.; Allen, H. C. J. Phys. Chem. B 2004, 108, 2252. (60) Iuchi, S.; Izvekov, S.; Voth, G. A. J. Chem. Phys. 2007, 126, 124505. (61) Taylor, R. S.; Dang, L. X.; Garrett, B. C. J. Phys. Chem. 1996, 100, 11720. (62) Radu¨ge, C.; Pflumio, V.; Shen, Y. R. Chem. Phys. Lett. 1997, 274, 140. (63) Miranda, P. B.; Shen, Y. R. J. Phys. Chem. B 1999, 103, 3292. (64) Baldelli, S.; Schnitzer, C.; Shultz, M. J. Chem. Phys. Lett. 1999, 302, 157. (65) Buch, V.; Tarbuck, T.; Richmond, G. L.; Groenzin, H.; Li, I.; Shultz, M. J. J. Chem. Phys. 2007, 127, 204710. (66) Herce, D. H.; Perera, L.; Darden, T. A.; Sagui, C. J. Chem. Phys. 2005, 122, 024513. (67) Israelachvili, J. N. Intermolecular and Surface Forces; Academic Press: New York, 1991. (68) Shen, Y. R.; Ostroverkhov, V. Chem. ReV. 2006, 106, 1140. (69) Smondyrev, A. M.; Voth, G. A. Biophys. J. 2002, 82, 1460. (70) Scatena, L. F.; Brown, M. G.; Richmond, G. L. Science 2001, 292, 908. (71) Scatena, L. F.; Richmond, G. L. J. Phys. Chem. B 2001, 105, 11240.

JP805304J

Hydrated Excess Proton at WaterâHydrophobic Interfaces - The

Recommend Documents

Hydrated Excess Proton at WaterâHydrophobic Interfaces - The