Hydrophobic interaction and hydrogenbond network for a methane pair in liquid water
 *Department of Chemistry, Princeton University, Princeton, NJ 08544;
 ^{†}NEC Laboratories America, Inc., 4 Independence Way, Princeton, NJ 08540;
 ^{§}California Institute for Quantitative Biomedical Research, Departments of Biopharmaceutical Sciences and Biochemistry and Biophysics, University of California, San Francisco, CA 94143; and
 ^{¶}Department of Molecular Biology, Princeton University, Princeton, NJ 08544
See allHide authors and affiliations

Communicated by Morrel H. Cohen, Rutgers, The State University of New Jersey, Bridgewater Township, NJ, December 19, 2006 (received for review February 3, 2006)
Abstract
We employ fully quantummechanical molecular dynamics simulations to evaluate the force between two methanes dissolved in water, as a model for hydrophobic association. A stable configuration is found near the methane–methane contact separation, while a shallow second potential minimum occurs for the solventseparated configuration. The strength and shape of the potential of mean force are in conflict with earlier classical forcefield simulations but agree well with a simple hydrophobic burial model which is based on solubility experiments. Examination of solvent dynamics reveals stable water cages at several specific methane–methane separations.
Hydrophobicity is the molecular driving force behind numerous important biological processes, including protein folding and the formation of biological membranes (1–3). A quantitative understanding of hydrophobic interactions is crucial for modeling protein structures, protein functions, or manipulation of hydrophobic nanoparticles in aqueous solutions (4).
Experimentally, the strength of the hydrophobic effect (hydration potential) can be measured by the solubility of hydrocarbons (5, 6). However, the detailed shape of the potential of mean force (PMF) between two hydrocarbon molecules has only been probed indirectly (7, 8). Considerable effort has been expended in studying hydrophobic interactions and hydration by using classical Lennard–Jones potentials and various water models (9–11). The model parameters were typically chosen for consistency with bulk thermodynamic quantities. However, the hydrophobic effect for dissolved molecules originates largely from the hydrogenbond network in the first solvation shell (12), and the properties of interfacial water differ substantially from those of bulk water. Indeed, hydrogen bonding remains quite difficult to represent effectively with simple (atom–atom) molecularmechanics force fields (13).
In firstprinciples molecular dynamics (FPMD) (14), interatomic forces are derived directly from quantummechanical calculations. FPMD has been successfully applied to ice (15), water clusters (16), bulk liquid water (17), and water in the solvation shell of a dissolved ion (18) or methane (19). Here, we report determination of the PMF between a pair of methane molecules in water by FPMD. In classical simulations, the general features of the PMF are a stable freeenergy minimum at contact separation, with a second but pronounced freeenergy minimum at a distance where the two methanes are separated by a single layer of solvent (20). However, rather small changes in the classical methane–water interaction parameters can lead to reordering of the stability of these two minima (21). Our quantummechanical simulation reveals an effective hydrophobic surface tension. The result is a stable configuration of two methanes near contact separation with only a shallow potential minimum for the solventseparated configuration. The depth of the stable potential minimum is roughly in accord with solubility measurement (5) but is much deeper than the potential minimum found in previous classical forcefield MD simulations.
Our simulations also allowed us to examine the dynamics of water surrounding a hydrophobic solute. We found stable water cages at several specific methane–methane separations. These clathratelike cages illustrate the effects of solute size on the local water structure.
Overview of Approaches
If one simulates Brownian motion of two methanes in water, each methane is in constant collision with water molecules, so naively the two methanes should drift away from each other. Instead, even though there is hardly any direct attractive force between two methane molecules at distances >4.5 Å, they are bound together by the surrounding water during the simulation time. To quantify the effect of solvent, we can define an effective potential W(r) between the two solute molecules. For a given separation r between the two molecules of interest, the (thermal ensemble) mean effective force acting on them is where T and V are the temperature and volume of the system, and N is the number of particles. We shall call W(r) the potential of mean effective force (PMEF), and the corresponding Helmholtz free energy is denoted by F(r; T, V, N). A related and more widely used quantity to describe interactions in fluids is the PMF. The PMF is defined through the equilibrium probability of finding two molecules a certain distance apart in the solvent. Specifically, the PMF is given by where g(r) is the radial distribution function of the solute molecules. The normalization is chosen such that g(r) tends to unity as r → ∞. The average effective force is related to the PMF by We shall refer to the extra term, 2k_{B}T/r, as the volumeentropy force. It is due to changes in the free volume of the spherical shell available to the solutes. A derivation of Eq. 2 is provided in the supporting information (SI).
Various methods have been used to compute the PMF between two methane molecules in water: freeenergy perturbation (22), thermodynamic integration (23), and umbrella sampling (11, 24). In the present study, we used a constrainedmoleculardynamics method (25, 26) in which a holonomic constraint fixes the methanepair separation. To this end, we introduced a Lagrange multiplier λ to fix the distance between the two methane molecules. The new Lagrangian ′ and the constraint force f are where is the original Lagrangian of the Nparticle system, and r_{1} and r_{2} are the positions of the two methanes. The time average of the constraint force, Eq. 3, is equivalent to the average effective force f̄(r). A simple integration of f̄(r) over constrained variable r gives the freeenergy difference ΔW(r). Finally, the PMF (up to a constant) is (see the SI) The distinction between w(r) and W(r) is critical when comparing with experiments and relating various approaches.
Results
Potential of Mean Force.
For each methane–methane separation, the constraint force is recorded at each MD time step, and the average is taken over the simulation time (6–10 ps). An example of data collection is shown in Fig. 1. Because of the methane C–H stretch modes, the constraint force between two methane molecules switches directions from repulsive to attractive roughly every 27 fs. As shown in Fig. 1, the time average of the constraint force converges in a few picoseconds. In previous, classical mechanics studies, using a Lennard–Jones potential and water models, the simulation time was 500 ps for each 0.5Å window to converge the PMF within 0.1 kcal·mol^{−1} (11). Runtime considerations for our fully quantummechanical MD method induced us to use shorter simulation times, typically 6 ps for each methane–methane separation. We tested the convergence of the average constraint force for several separations between 3.8 and 5.8 Å using a longer MD simulation time (10 ps). We also checked the difference between using a smaller unit cell with 63 waters, versus a larger unit cell with 139 waters, for methane–methane separations of r = 4.4 and 6.0 Å. The resulting differences were typically of the order of 10^{−4} a.u. in the constraint force (see the SI for details of data collection and error analysis). Forces are expressed in atomic units: 1 a.u. = 8.2353 × 10^{−8} J/m.
Fig. 2a shows the results of meaneffectiveforce calculations from r = 2.8 to 7.8 Å. When the methane–methane separation is less than the contact distance (≈3.9 Å), the overlap of methane electronic orbitals results in a strong repulsive force, which dominates effects from the surrounding water. To save computing time, we used the repulsive force between two methane molecules in vacuum when r is <3.6 Å. For methane–methane separations between 4.2 and 7.8 Å, we found a net attractive force, except near 6.5 Å. The PMEF, obtained by integrating the mean effective force, is shown in Fig. 2b. A small potential barrier around 6.2 Å separates the contact potential minimum from the second potential minimum near 6.7 Å. The PMF in Fig. 2b is obtained by subtracting the volumeentropy term, 2k_{B}T ln r, from the PMEF.
Hydration Structure.
Our quantummechanical simulations allow us to study in detail the hydration structure around a methane pair. Near the pair, water forms a clathratelike cage structure with a thickness of ≈5.5 Å (27). In Fig. 3, we plot the number of water molecules within the first solvation shell (coordination number) versus the methane–methane separation. Interestingly, the coordination number shows plateaus at several methane–methane separations.
The presence of plateaus in Fig. 3 suggests the existence of stable well structured hydration shells at particular methane–methane separations. The dynamic behavior of water in the solvation shell can be visualized by plotting the rms displacement of each water molecule. For each MD step, we transform (translation plus rotation) the coordinate frame to the one where the two carbon atoms are stationary. Fig. 4 shows the rms displacement over simulation times of 6 ps for several methane–methane separations. The rms displacements are smallest at r = 4.8 Å, which is in the center of the first plateau in the coordination number plot.
Discussion
Experimental Determination of the PMF.
An important motivation for computing PMFs is to obtain accurate solvation parameters for molecular modeling studies. Solvation free energies have been obtained from measurements of partition coefficients of hydrocarbons between aqueous and nonpolar phases (5, 28). The transfer of a hydrophobic solute into water is accompanied by an increase in free energy, part of which results from structural changes in the water around each solute molecule. At equilibrium, this transfer free energy is balanced by an increase in volume entropy. Experimentally, the transfer free energy for liquid nalkanes to water is highly linear with respect to the solventaccessible area (29). So to a good approximation, when the solventaccessible area changes, the corresponding freeenergy change ΔG can be described via a surfacetension parameter σ, where ΔA is the change in the solventaccessible area.** Before the transfer of the solute, the solventaccessible area A is buried within the nonpolar solute liquid; after the transfer, the area A is immersed in liquid water. The surfacetension parameter, σ, is estimated from experiment to be 47 cal·mol^{−1}·Å^{−2}, valid for all nalkanes in the series from methane up to decane (5, 30).^{††}
The surfacetension model can be used to estimate the PMF between two methanes (31).^{‡‡} The solventaccessible area of a methane molecule is defined as a sphere with radius ρ_{0} = (1.95 + 1.40) Å, the sum of methane and water radii. As illustrated schematically in Fig. 5, the solventaccessible area of a methane pair is smaller when the two are in contact than when they are separated. In the contact configuration, the “buried” surface area is 20.8955% of the total available solventaccessible area, 8πρ_{0}^{2}. Thus, the depth of the PMF should be roughly ΔG = σ·0.2·8πρ_{02} ≈ 2.7 kcal·mol^{−1}. Within the surfacetension model, the hydrophobic force between two methanes is where r is the separation between the two methane molecules. So surface tension gives rise to a constant attractive force (0.00083 a.u., or 0.98 kcal·mol^{−1}·Å^{−1}) for r between 3.9 Å (contact separation) and 6.7 Å (separated by a water diameter).
To compare the estimated PMF from the surface tension model with our calculated PMF, we need to correct for the temperature difference between the experiment (room temperature) and the simulation (T = 343 K). The freeenergy change in Eq. 4 is determined by the molar volume partition coefficient K_{M} between liquid hydrocarbons and water, where R is the gas constant and V_{s} (V_{w}) is the molar volume of solute (water) (32). The hydrophobic free energy increases rapidly above room temperature. In addition to the explicit temperature dependence in RT, the solubility (measured by K_{M}) of hydrocarbons in water decreases with increasing temperature. Because of the low boiling temperatures of small nalkanes (from methane to butane), temperaturedependent solubility data are only available for partitioning between gas phase and water (33, 34). From these gas/water partition coefficients, we obtained a linear relation between the transfer free energy and the solventaccessible area which shows the transfer free energy per solventaccessible area increases by ≈40% between room temperature and T = 343 K (details of this analysis can be found in the SI). Therefore, we estimate that the adjusted potential depth of the PMF at room temperature is around ≈2.8 kcal·mol^{−1} from our simulation, in good agreement with the estimate from the surfacetension model. Indeed, the shape of our calculated PMF is also in good agreement with the surfacetension model [using σ(T = 343 K) = 65 cal·mol^{−1}·Å^{−2} in Fig. 2b].
Comparison of QuantumMechanical Versus Classical PMF.
The PMF we obtained from a fully quantummechanical simulation reasonably resembles the surfacetension potential but is distinctly different from earlier classical results. The PMF in Fig. 2b has a characteristic contact minimum at 3.9 Å with potential depth of ≈3.9 kcal·mol^{−1}, which is much deeper than all previous studies [usually between 0.5 kcal·mol^{−1} (10) and 0.9 kcal·mol^{−1} (11)]. Even after adjusting to room temperature, we predict the potential depth of the PMF to be around −2.8 kcal·mol^{−1}.^{§§}
The question is whether and why classical forcefield simulations underestimate the hydrophobic effect. Not unexpectedly, there are many differences between quantum and classical MD simulations in the methane/water system, such as the angular distribution of waters and hydrogenbondring statistics (19). More studies comparing the two approaches remain to be done. However, we expect that the qualitative feature of the PMF in Fig. 2b, a small barrier separating the first and second potential minima, should be a rather robust feature of the densityfunctional theory calculations. This is in contrast to the occurrence of a pronounced solventseparated minimum in classical simulations, which resembles the feature of hardsphere liquids and may result from overly repulsive shortdistance forces of Lennard–Jones potentials.
Stability of the Hydration Shell.
For a more stable aqueous solvation shell (like the one at r = 4.8 Å in Fig. 4), one would expect a fewer exchanges of water between the solvation shell and bulk water, less variation of the coordination number, and a higher ratio of pentagon to hexagon rings in the nonshortcircuited hydrogenbonded network distribution (35). Indeed, the meansquared deviation of the coordination number recorded during a 6ps MD run is 3.7 at r = 4.4 Å but decreases to 1.6 at r = 4.8 Å.
To make sure that the results in Fig. 4 were reproducible, we computed rms displacements at methane–methane separations of r = 4.4 and 4.8 Å starting from new initial configurations. Each frame in Fig. 6 represents the rms displacements over a 4ps MD simulation. The rms displacements at r = 4.8 Å are consistently smaller than the rms displacements at r = 4.4 Å. For both r = 4.4 and 4.8 Å, the water in the cap regions is typically more fluid (larger rms displacements) than the water in the equatorial plane between the methanes.
Given the small rms displacements of the waters in the solvation shell at r = 4.8 Å in both Figs. 4 and 6, one might conjecture that the waters assume a unique stable structure near the methanes. To test this idea, we determined the average water positions for two independent MD simulations at r = 4.8 Å. Even after applying all possible symmetry operations with respect to the methane–methane axis, we found that the two water structures could not be overlaid on each other. The hydrogenbonded rings near the equator of the hydration shell coincided well between the two simulations, but there remained considerable mismatch between water locations in the cap regions.
The stability of the solvation shell at r = 4.8 Å, even in the absence of a unique water structure, can be understood as follows. For a particular methane–methane separation, the surrounding water may be able to form a well packed cage, as shown schematically in Fig. 7a. For a larger or smaller methane separation, the surrounding water will not pack well, as indicated in Fig. 7b. Poor water packing will allow competing configurations with different numbers of waters in the solvation shell, leading to large rms displacements and rapid exchange of waters with the bulk. In contrast, for the well packed cage, rms displacements of water will be small and there will be relatively little exchange of waters with the bulk, even though the detailed packing in the solvation shell may not be unique.
Conclusion
The PMF for two methane molecules in water was calculated by using a constrained FPMD method, where the effective force between two methanes is directly computed by a Lagrange multiplier. We highlighted the importance of a volume entropy term which leads to a distinction between the PMF (free energy due to the intermolecular interaction) and the PMEF (total free energy).
The PMF has a deep stable minimum near contact and a shallow solventseparated minimum. The magnitude of the hydrophobic force agrees well with measurements of the solubility of hydrocarbons in water. Our results therefore provide a good starting point to parameterize hydrophobic association in largescale simulations. There is a considerable difference between the PMF in the current study (based on quantummechanical densityfunctional theory) and previous results (based on classical Lennard–Jones potentials and water models). The stability of local water structures near a methane pair depends on the methane–methane separation. We analyzed the dynamic behavior of water in the solvation shell and found well structured solvation shells for particular methane–methane separations.
Methods
In our FPMD simulations, the interatomic forces were calculated in the Born–Oppenheimer approximation within densityfunctional theory. We used the PBEGGA exchangecorrelation functional (36) and ultrasoft pseudopotentials (37). The electronic Kohn–Sham orbitals were expanded in a planewave basis set with a cutoff energy of 25 Ryd. The details of the densityfunctional theory in the simulation of water/methane systems (cutoffs, exchangecorrelation functionals, pseudopotentials, etc.) were considered in several earlier studies (38–40). The time step for MD was 0.24 fs. Two methane molecules and a number of waters were put in a cubic unit cell with periodic boundary conditions. For methane–methane separation smaller than 6.0 Å, a unit cell of side length 12.5 Å containing 63 waters was used. For methane–methane separation larger than 6.0 Å, a unit cell of side length 16.2 Å with 139 waters was used. We started data collection after an equilibration time of 3 ps. Throughout the simulation, the temperature was thermostated at 343 K. The water structure has been checked carefully at this temperature (41, 42). The simulation temperature is chosen to ensure that the system is in a liquid phase; although a complete phase diagram of the ab initio water (17) has not been mapped out, the freezing temperature is likely to be >273 K (40).
Acknowledgments
We are grateful to Frank Stillinger, and to David Chandler for pointing out to us the role of hydrophobicity at extended surfaces. This work was supported by National Science Foundation Grants CHE0121432 and DMR0313129. J.L.L. is a Kenda Foundation Golden Jaden Fellow.
Footnotes
 ↵^{‖}To whom correspondence should be addressed. Email: wingreen{at}princeton.edu

Author contributions: R.C. and N.S.W. designed research; J.L.L. performed research; and C.T. contributed new reagents/analytic tools.

The authors declare no conflict of interest.

See Commentary on page 2557.

This article contains supporting information online at www.pnas.org/cgi/content/full/0610945104/DC1.

↵** For a linear hydrocarbon chain such as an alkane, its solventaccessible area scales linearly with volume. Indeed, the free energy of the creation of a small cavity in water can be shown to be approximately linear in excluded volume. Therefore, the “cavity volume” freeenergy contribution can be absorbed into the surface tension parameter σ in Eq. 4.

↵†† The value σ = 47 cal·mol^{−1}·Å^{−2} is obtained by adopting Flory–Huggins theory (32), whereas Chan and Dill (6) deduced the value 34 cal·mol^{−1}·Å^{−2} from “classical theory” and cyclohexane–water transfer data (see discussion in ref. 30).

↵‡‡ Different surface measures may be more adequate for different thermodynamic properties (see, for example, discussions in ref. 31).

↵§§ Had we used σ = 34 cal·mol^{−1}·Å^{−2} for the surface tension parameter (30), ΔG is ≈1.9 kcal·mol^{−1}, still considerably larger than classical forcefield MD results. The difference between the two estimated ΔG values, 1.9 and 2.7 cal·mol^{−1}·Å^{−2}, may well be within the uncertainties of current densityfunctional theory approximations.
Abbreviations
 FPMD,
 firstprinciples molecular dynamics;
 PMEF,
 potential of mean effective force;
 PMF,
 potential of mean force.
 Received February 3, 2006.
 © 2007 by The National Academy of Sciences of the USA
References
 1.↵
 2.↵
 3.↵
 4.↵
 5.↵
 Sharp KA,
 Nicholls A,
 Fine RF,
 Honig B
 6.↵
 7.↵
 Thompson PT,
 Davis CB,
 Wood RH
 8.↵
 Wood RH,
 Thompson PT
 9.↵
 10.↵
 11.↵
 11.↵
 Bromberg S,
 Dill KA
 13.↵
 Morozov AV,
 Kortemme T,
 Tsemekhman K,
 Baker D
 14.↵
 15.↵
 16.↵
 17.↵
 18.↵
 19.↵
 20.↵
 21.↵
 22.↵
 23.↵
 Pearlman DA
 24.↵
 25.↵
 26.↵
 27.↵
 HeadGordon T
 28.↵
 McAuliffe C
 29.↵
 30.↵
 31.↵
 32.↵
 33.↵
 Clever HL,
 Young CL
 Battino R
 34.↵
 Lide DR
 Gevantman LH
 35.↵
 36.↵
 37.↵
 38.↵
 Ikeda T,
 Terakura K
 39.↵
 Xu X,
 Goddard WA III.
 40.↵
 Sit PHL,
 Marzari N
 41.↵
 Pasquarello A,
 Petri I,
 Salmon PS,
 Parisel O,
 Car R,
 Tóth É,
 Powell DH,
 Fischer HE,
 Helm L,
 Merbach AE
 42.↵
Citation Manager Formats
Related Article
 Dissecting hydrophobicity Feb 13, 2007