J. Phys. Chem. B 2009, 113, 767–778
767
Electrostatic Polarization Effects and Hydrophobic Hydration in Ethanol-Water Solutions from Molecular Dynamics Simulations Yang Zhong and Sandeep Patel* UniVersity of Delaware, Department of Chemistry and Biochemistry, 238 Brown Laboratory, Newark, Delaware 19716 ReceiVed: August 7, 2008; ReVised Manuscript ReceiVed: October 8, 2008
Nonadditive electrostatic force fields based on the charge equilibration formalism coupled with long timescale molecular dynamics simulations are used to investigate the microscopic structural aspects of hydrophobic hydration in ethanol-water solutions. Employing a combination of polarizable ethanol and water force fields (developed independently), we find that solution properties are satisfactorily reproduced across the ethanol mole fraction range between 0.1 and 0.9. Solution densities are predicted within 3.6% of experimental measurements, while excess mixing enthalpies are overestimated as in earlier studies. The solvation free energy of ethanol in infinite dilution is determined via thermodynamic integration to be 5.70 ( 0.23 kcal/ mol, overestimating the free energetics of solvation relative to experiment (5.01 kcal/mol). Bulk solution dielectric constants and diffusion constants reproduce experimental trends and are in reasonable agreement across the ethanol concentration range studied. Because of explicit accounting of induction effects, ethanol and water exhibit varying molecular dipole moment distributions with concentration. The polarizable ethanol model, possessing higher condensed-phase polarizability relative to the TIP4P-FQ water model (4.54 Å3 versus 1.1 Å3, respectively), displays greater variation upon perturbation by the electric field of water. With regard to hydrophobic hydration, the current force fields indicate positive hydrogen bonding excess for water in the dilute ethanol concentration range, consistent with previous theoretical and experimental studies. Strikingly, we find that there are both positive and negative hydrogen bond excess contributions within the first hydration shell of both the ethanol hydroxyl oxygen and ethylene carbon atoms. The larger positive contributions dominate the overall hydrogen bonding patterns to yield overall net positive excesses. Moreover, we do not find evidence of excess hydrogen bonding vicinal to the nonpolar moieties as has been suggested based on the “iceberg”like models proposed by Frank and Evans. The present results suggest negative excess in the regions surrounding the alkyl groups that Vis-a`-Vis corresponds to a reduction in the average molecular dipole moment of resident water molecules due to smaller dipole induction in the weaker electrostatic fields of the nonpolar groups. I. Introduction Hydrophobic hydration is a key element in myriad physical and physiological processes, ranging from solvation of small molecule solutes (noble gases, alcohols, amines, etc), protein folding and stability, micellization, solute partitioning near lipid/ water interfaces, and wetting/dewetting phenomena.1-7 Despite the decades of work attempting to arrive at a unified microscopic picture and molecular-scale understanding of the effect that can bridge to macroscopic/continuum system properties, there continues to be no single consensus description. Specifically addressing the hydrophobic hydration effect, there are numerous theories regarding the nature of water structure, that is, water hydrogen bonding patterns in solutions of amphiphiles (and the perturbations of this structure from the pure water system). We do not consider all of the vast literature on the subject here, but consider some general pertinent ideas. Viewing the hydrophobic effect as some way being defined by water structure, the classical “iceberg” or “clathrate” model of Frank and Evans8 offers a plausible description of enhanced water structure (hydrogen bonding) vicinal to nonpolar solutes (or nonpolar moieties in mixed systems). The model accounts for greater water-water hydrogen bonded interactions that can compensate for entropic * To whom correspondence should be addressed. E-mail: sapatel@ udel.edu.
losses due to increased structure. Though appealing, there is no direct experimental evidence.3,4,9 Moreover, based on arguments relating differences in freezing entropies for pure water to entropy changes accompanying hydration of small solutes of the dimensions of amino acids, Kauzmann argues that first coordination shell waters around nonpolar solutes are less ordered than in ice.10,11 Other structural explanations of the hydrophobic effect invoke the necessity of forming large “cavities” in the network of small water molecules, thus leading to “disaffinities” of oil for water12-15 Further detailed discussions are found in refs 2, 5, 6, and 16 Alcohol-water systems are considered ideal test cases for studying hydrophobic hydration; we thus choose to explore the ethanol-water system for the present work, particularly since we would like to compare/contrast the results obtained using next-generation force fields for molecular dynamics based on charge equilibration models to alternative models, thus allowing an assessment of the physical descriptions afforded by these differing approaches. Moreover, there exists significant experimental data on alcohol-water mixtures to which simulation results can be referenced.3,4,9,17 Complementing experimental efforts, classical statistical mechanical (empirical force field) approaches to modeling alcohol-water solutions have contributed much to our understanding of these systems. Early Monte Carlo studies investigating small, essentially infinitely dilute
10.1021/jp807053p CCC: $40.75 2009 American Chemical Society Published on Web 12/30/2008
768 J. Phys. Chem. B, Vol. 113, No. 3, 2009 alcohol systems18-20 found support for the Frank and Evans model8 of an enhanced “icelike” water cage structure surrounding the methanol methyl moiety. More recent studies using molecular dynamics simulations and differing interaction models have presented evidence suggesting that the water structure in the vicinity of alcohol hydrophobic groups is not significantly perturbed.21-24 With improved computational hardware resources, ab initio molecular dynamics studies25,26 have begun to shed new light on the microstructure in these systems, again indicating that the presence of methanol does not perturb the local water structure in the infinite dilution limit while also suggesting the formation of a bipercolating network.25,27,28 Moreover, ab initio molecular dynamics methods provide insight into the mechanisms of proton transfer in such systems and have revealed subtle differences in the proton transfer mechanism in concentrated methanol solutions27,28 relative to pure bulk water. Until recently, classical molecular dynamics simulations have employed fixed-charge models to incorporate electrostatic effects, generally with a mean-field approximation of the fixed charges to account for nonadditive electronic contributions. Despite the wide success of these potentials, the lack of explicit treatment of electronic polarization does indeed limit the physics one can study with such potentials. Classical alcohol force fields have only recently incorporated explicit electronic induction effects. Some of the earliest polarizable alcohol models have addressed bulk and interfacial thermodynamic properties. Dang et al. developed a dipole polarizable force field of a six-site methanol molecule with fixed charges and a molecular polarizability placed on the oxygen atom to allow for the induction of dipoles via an iterative self-consistent procedure performed during each dynamics step; this model captures bulk and interfacial energetics, as well as bulk structure, in agreement with experimental data.29 More recently, such models have been employed to study liquid-vapor phenomena and interfacial properties of methanol-water solutions.30 Gao et al. presented a dipole polarizable methanol force field (simultaneously parametrized with other alcohols) in much the same spirit as the Dang model wherein induced atomic dipole moments are determined self-consistently at each Monte Carlo step.31 Caldwell et al. contributed an additional dipole polarizable model of methanol representing accurately the density, vaporization energetics at 298 K, and diffusion coefficient.32 Finally, Patel et al.33,34 recently reported two new alcohol potentials for methanol and ethanol based on the charge equilibration formalism; however, the authors did not extend their analysis to include aqueous solutions. Unfortunately, studies of aqueous ethanol solutions (and aqueous alcohol solutions in general) using polarizable, or nonadditive, force fields are rare due to the lack of systematically developed models available. Recently, Yu et al.35 reported a polarizable force field for methanol-water solutions; the authors discussed bulk and solution properties, emphasizing the ability of their polarizable model to reproduce experimental dielectric constants for the mixtures consistently across the concentration range studied. However, the authors did not discuss structural properties in terms of microstructure or network formation. Noskov et al.36 studied ethanol-water solutions using a Drudeoscillator model for ethanol and the SWM4-DP37 water model. Hernandez-Cobos et al. studied methanol-water solutions using polarizable models.38 The Drude model suggests that a depletion of water-water hydrogen bonds in the first solvation shell of the ethanol methyl group is compensated by a higher solvent density (compaction) in the second shell; this behavior holds for ethanol mole fractions up to 0.4, after which the second
Zhong and Patel shell cannot accommodate sufficient water molecules to compensate for the net depletion in the first shell. This picture suggests that the Frank and Evans iceberg model8 is still valid and in agreement with neutron diffraction data. HernandezCobos,38 on the other hand, suggest from their Monte Carlo simulations, that it is the methanol hydroxyl moiety that is hydrated; the authors further draw attention to the similarities of their work to the Monte Carlo simulations (using additive force fields) of Jorgensen and Madura39 showing a depletion of water-water hydrogen bonds in the methyl cage region relative to average bulk hydrogen bonding. Other models based on the SIBFA40 and TCPE41 methods, parametrized to high-level ab initio calculations on model cluster geometries of methanol and water, have also been introduced; these models, however, have not been extensively applied to long-time molecular dynamics simulations of alcohol-water systems to explore the condensed phase properties predicted by these potentials. In the following, we adopt the philosophy of applying force fields for water and ethanol developed independently to study mixture properties (as has been done in the past)36,42,43 and in particular investigate the nature of hydrophobic hydration as manifested in hydrogen bonding patterns evolving from the application of charge equilibration models of associating (hydrogen-bonding) fluids. In so doing, we hope to contrast the current microscopic picture to that obtained from other physical models of liquids that explicitly consider polarization, or nonadditive effects. Section II describes the force fields and simulation methods/protocols. Sections IVA-IVD discuss bulk solution properties for ethanol mole fractions 0.1 to 0.9. Sections IVC and IVD discuss the nature of hydrogen bonding and its connection to hydrophobic hydration. Section V presents an analysis of clustering in the various solutions. Section VI presents conclusions. II. Force Fields The current simulations employ water-water, alcohol-alcohol, and water-alcohol force fields based on the CHARMM potentials44 for intra- and intermolecular interactions, but extended to include nonadditive electrostatic effects using a charge equilibration (fluctuating charge) formalism. The fluctuating charge model has been applied to various systems over the past decade.33,34,45-50 The method derives from the density functional theory of atoms in molecules as formulated by Yang and Parr.51 More fundamentally, the method is founded on Sanderson’s idea of electronegativity equalization.52,53 In the density functional sense, electronegativity equalization amounts to the equalization of the chemical potential in space. In a molecule, this translates to the redistribution of charge among constituent atoms so as to equalize the electronegativity (chemical potential) at each point. The electrostatic energy of an N-atom molecule in vacuum is N
Eelectrostatic )
N
N
N
∑ χi0Qi + 21 ∑ ∑ ηijQiQj+ ∑ ΦiQi i)1
i)1 j)1
(1)
i)1
The χ’s are atom electronegativities (formally the Mulliken electronegativities) related to the chemical potential of an electron gas surrounding a nucleus by
µi )
∂E ∂E ) -χi0 ) -e ∂N ∂qi
(2)
Hydrophobic Hydration in Ethanol-Water Solutions
J. Phys. Chem. B, Vol. 113, No. 3, 2009 769
The η’s are the atom hardnesses. These values represent a “resistance” to electron flow to/from an atom. The last term in equation 1 represents contributions from an external potential. Heterogeneous elements are derived from the atom type values based on the combining rule:45,54
ηij(Rij, ηi, ηj) )
1 (η + ηj) 2 i
(3)
1 1.0 + (ηi + ηj)2Rij2 4
where Rij is the separation between atoms (or more generally sites) i and j. This local screened coulomb potential has the correct 1/r limiting behavior for separations greater than about 2.5 Å. This interaction is computed for 1-2, 1-3, and 1-4 sites (sites included in bonds, angles, and dihedrals); sites separated by 5 or more sites interact via a standard Coulomb interaction. In the case of interacting molecules, the interaction between sites on different molecules is also of the Coulomb form. We note here that an alternate model for the shielding interaction is the Coulomb overlap integral between atom-centered Slater orbitals. This approximation is computationally tractable for rigid systems where the interaction matrix elements can be computed once at the start of a simulation and then used for the duration; however, for flexible systems, calculation of the overlap integral leads to nontrivial cost. Finally, any external potential is included in the last term of equation 1. In contrast to dipole polarizable models, the fluctuating charge model, by including explicit charge-charge interactions, incorporates higher-order multipoles that are neglected in the former case. With respect to charge dynamics, an extended Lagrangian formalism is used to propagate the charges in time with some general charge constraint, thus strictly providing for electronegativity equalization at each dynamics step. For the present work, molecular charge neutrality is imposed; however, this can be modified to permit charge transfer between molecular entities. The system Lagrangian and associated constraint on total molecular charge for each molecule are then M
L)
Ni
M
∑∑ i)1
Ni
∑ ∑ 21 mQ,iRQ˙iR2 - E(Q, r) -
1 2 + miRr˙iR 2 R)1 i)1
R)1
M
Ni
i)1
R)1
∑ λi ∑ QiR (4) where the first two terms represent the nuclear and charge kinetic energies, the third term is the total potential energy and the fourth term is the molecular charge neutrality constraint with λi the Lagrange multiplier for each molecule, indexed by i. The fictitious charge dynamics, analogous to wave function dynamics in Car-Parrinello (CP) type methods,55 are determined by a charge “mass” (adiabaticity parameter in CP dynamics). The units for this mass are (energy time2/charge2). The charges are propagated based on forces arising from the difference between the average electronegativity of the molecule and the instantaneous electronegativity at an atomic site. The nonbond interaction is of the Lennard-Jones type (standard in the CHARMM56 implementation)
E(R) )
(
∑ ij εij
12 Rmin,ij
rij12
-2
6 Rmin,ij
rij6
)
(5)
TABLE 1: Electrostatic and Nonbond Parameters for Ethanol and TIP4P-FQa atom Rmin/2 molecule type εi (kcal/mol) (Å) χ (kcal/mol)/e η (kcal/mol)/e2 ethanol
CT2 OH1 H HA2 CT3 HA3 TIP4P-FQ OT OM HT
0.070 0.1751 0.046 0.020 0.07 0.02 0.2862
2.015 1.75 0.224 1.32 2.015 1.32 1.7729
0.046
0.224
78.892 145.836 0.0 77.961 83.874 73.576 0.0 68.49 0.0
114.917 169.069 284.497 275.785 132.191 275.785 0.0 371.6 353.0
a
Ethanol parameters taken from Patel and Brooks.34 Water parameters for TIP4P-FQ are taken from Rick et al.50 (Cross terms are not shown and the reader is referred to the work of Rick et al.50 for further details).
where the summation runs over all nonbonded atom pairs. Table 1 shows the electrostatic and van der Waals parameters for the current models. The water-water interactions are treated using the TIP4P-FQ water model as developed by Rick et al.50 The ethanol-ethanol interactions are modeled using the fluctuating charge models developed by Patel and Brooks.33,34 Cross interactions between water and ethanol are treated using the homogeneous interactions parameters with a geometric mean combining rule for the Lennard-Jones energy well parameter and an arithmetic mean combination rule for the separations, Rmin. This is the commonly applied Lorentz-Berthelot combination rule, applied here in order to maintain consistency with the CHARMM force field. We also permit the ethanol to interact with water hydrogen atoms as in the original CHARMM fluctuating charge force field as developed by Patel and Brooks;45 we feel that this allows for a more consistent representation in terms of the global development of a CHARMM polarizable force field. Intramolecular interaction potentials (bond stretching, angle bending, dihedral rotation) are retained from the nonpolarizable CHARMM27 (C27) force field. We note that the Drude polarizable models for alcohols introduce off-diagonal interactions between the oxygen atoms of the water and alcohol to allow greater flexibility in terms of representing the solute-water cross-interactions. For the present study, we use only the combining rules just mentioned. III. Simulation Methods Bulk solution simulations were performed within the isothermal-isobaric ensemble (NPT) at 298K and 1 bar. A constant temperature of 298 K was maintained with a Hoover thermostat,57 and constant pressure was maintained via a Langevin piston method.58 Systems of water and ethanol molecules in cubic boxes of varying dimensions were generated from pre-equilibrated pure solvent systems. Table 2 lists the molecule numbers used for each concentration ranging from 0.1-0.9 ethanol mole fraction. The systems were equilibrated for 1 ns (ns) under constant temperature and pressure at 298 K and 1 bar, followed by 4 ns of molecular dynamics using a Verlet leapfrog integrator with time step of 0.0005 ps (0.5 fs). Nonbonded interactions were switched to zero via a switching function from 9 to 10 Å. We choose this cutoff scheme to be consistent with the water model used, instead of using a value for the pure ethanol state point (the pure ethanol model was parametrized with a switching function between 11.5 and 12.5 Å) particularly as we are concerned with the dilute concentration ranges for our analysis of hydrophobic solvation and making
770 J. Phys. Chem. B, Vol. 113, No. 3, 2009
Zhong and Patel
TABLE 2: Densities, Enthalpies of Mixing, Dielectric Constants and Self-Diffusion Coefficients of Ethanol and Water in the Mixtures, Computed from Simulations of 5ns (statistical uncertainties reported in parentheses) XEtOH 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
F (g/L)
Hmix (kcal/mol)
998.24 (15.62) 974.46 (12.07) 950.02 (10.97) 921.77 (10.43) 894.27 (9.69) 868.79 (9.41) 847.62 (9.91) 826.92 (11.47) 803.80 (12.14) 780.69 (13.06) 756.34 (14.39)
0.0 (0.1970) -0.2305 (0.1861) -0.3286 (0.1773) -0.3100 (0.1703) -0.2762 (0.1596) -0.2509 (0.1632) -0.2358 (0.1768) -0.2215 (0.1959) -0.1551 (0.2146) -0.0775 (0.2354) 0 (0.2554)
ε
72.63 (0.44) 58.93 (0.88) 47.23 (0.67) 45.04 (0.50) 36.58 (0.48) 30.36 (0.22) 24.63 (0.28) 23.91 (0.18) 18.34 (0.09)
Dw (10-5cm2/sec)
1.72 (0.04) 1.30 (0.02) 1.06 (0.03) 0.94 (0.06) 0.89 (0.05) 0.90 (0.09) 0.88 (0.05) 0.94 (0.11) 0.93 (0.25)
comparisons to previous experimental and computational studies. Moreover, it is deemed more consistent to apply the same cutoff for both dilute and concentrated alcohol regimes. Conditionally convergent long-range electrostatic interactions were treated using the smooth particle mesh Ewald59 method with 24 grid points in each dimension and screening parameter of κ ) 0.33. Charge degrees of freedom were maintained at 1 K using a Nose-Hoover bath; all charges were assigned to the same bath and no local “hotspots” in charge were observed. Ethanol and water charge degrees of freedom were assigned masses of 0.00007 and 0.000069 (kcal/mol ps2)/e2, respectively. Charge normalization was maintained over individual molecules with no intermolecular charge transfer; the size of the molecular units is of adequate spatial extent to limit molecular polarization to that dictated by Pauli exclusion in the condensed phase.60 To generate the starting configurations of bulk systems, all coordinates of an equilibrated pure water configuration were shifted up by half of a water cubic box length and then combined with equilibrated coordinates of pure ethanol that were shifted down by half of one ethanol cubic box length. IV. Results IVA. Thermodynamic Properties. To provide a basic validation of the combination of polarizable force fields for water (TIP4P-FQ) and ethanol (charge equilibration), we present select properties based on 5 ns bulk solution simulations. We first consider bulk solution densities and excess mixing enthalpies as a function of ethanol concentration. The predicted densities are in acceptable agreement (relative error
