A computer‐simulation study of hydrophobic hydration of rare gases and of methane. I. Thermodynamic and structural properties

Crustwater: Modeling Hydrophobic Solvation

We describe Crustwater, a statistical mechanical model of nonpolar solvation in water. It treats bulk water using the Cage Water model and introduces a crust, i.e., a solvation shell of coordinated partially structured waters. Crustwater is analytical and fast to compute. We compute here solvation vs temperature over the liquid range, and vs pressure and solute size. Its thermal predictions are as accurate as much more costly explicit models such as TIP4P/2005. This modeling gives new insights into the hydrophobic effect: (1) that oil-water insolubility in cold water is due to solute-water (SW) translational entropy and not water-water (WW) orientations, even while hot water is dominated by WW cage breaking, and (2) that a size transition at the Angstrom scale, not the nanometer scale, takes place as previously predicted.

On the Solute-Induced Structure-Making/Breaking Effect: Rigorous Links among Microscopic Behavior, Solvation Properties, and Solution Non-Ideality

We studied the solute-induced perturbation of the solvent environment around a solute species from a microscopic viewpoint and propose a novel approach to the understanding of the structure-making/breaking process, regardless of the type and nature of the solute-solvent interactions. Based on the Kirkwood-Buff fluctuation formalism, we present a rigorous statistical mechanics description of the evolution of the solvent structure around the solute, analyze its response to small perturbations of the ( TP) state conditions and composition of the system, and make direct connections between a few equivalent micro- and macroscopic manifestations as probes for, and targets of, experimental measurements. We illustrate the analysis with theoretical results from integral equation calculations of model fluids and experimental evidence from available data for a variety of aqueous electrolyte and nonelectrolyte real fluid solutions. Finally, we provide a critical discussion about the inadequacy underlying a widely used de facto criterion for the classification of structure-making/breaking solutes.

Methane Hydration-Shell Structure and Fragility

The influence of oily molecules on the structure of liquid water is a question of importance to biology and geology and many other fields. Previous experimental, theoretical, and simulation studies of methane in liquid water have reached widely conflicting conclusions regarding the structure of hydrophobic hydration-shells. Herein we address this question by performing Raman hydration-shell vibrational spectroscopic measurements of methane in liquid water from -10 °C to 300 °C (at 30 MPa, along a path that parallels the liquid-vapor coexistence curve). We show that, near ambient temperatures, methane's hydration-shell is slightly more tetrahedral than pure water. Moreover, the hydration-shell undergoes a crossover to a more disordered structure above ca. 85 °C. Comparisons with the crossover temperature of aqueous methanol (and other alcohols) reveal the stabilizing influence of an alcohol OH head-group on hydrophobic hydration-shell fragility.

A closure relation to molecular theory of solvation for macromolecules

We propose a closure to the integral equations of molecular theory of solvation, particularly suitable for polar and charged macromolecules in electrolyte solution. This includes such systems as oligomeric polyelectrolytes at a finite concentration in aqueous and various non-aqueous solutions, as well as drug-like compounds in solution. The new closure by Kobryn, Gusarov, and Kovalenko (KGK closure) imposes the mean spherical approximation (MSA) almost everywhere in the solvation shell but levels out the density distribution function to zero (with the continuity at joint boundaries) inside the repulsive core and in the spatial regions of strong density depletion emerging due to molecular associative interactions. Similarly to MSA, the KGK closure reduces the problem to a linear equation for the direct correlation function which is predefined analytically on most of the solvation shells and has to be determined numerically on a relatively small (three-dimensional) domain of strong depletion, typically within the repulsive core. The KGK closure leads to the solvation free energy in the form of the Gaussian fluctuation (GF) functional. We first test the performance of the KGK closure coupled to the reference interaction site model (RISM) integral equations on the examples of Lennard-Jones liquids, polar and nonpolar molecular solvents, including water, and aqueous solutions of simple ions. The solvation structure, solvation chemical potential, and compressibility obtained from RISM with the KGK closure favorably compare to the results of the hypernetted chain (HNC) and Kovalenko-Hirata (KH) closures, including their combination with the GF solvation free energy. We then use the KGK closure coupled to RISM to obtain the solvation structure and thermodynamics of oligomeric polyelectrolytes and drug-like compounds at a finite concentration in electrolyte solution, for which no convergence is obtained with other closures. For comparison, we calculate their solvation structure from molecular dynamics (MD) simulations. We further couple the 3D-RISM integral equation with the 3D-version of the KGK closure, and solve it for molecular mixtures as well as oligomeric polyelectrolytes and drug-like molecules in electrolyte solutions.

Small molecule hydration energy and entropy from 3D-RISM

Implicit solvent models offer an attractive way to estimate the effects of a solvent environment on the properties of small or large solutes without the complications of explicit simulations. One common test of accuracy is to compute the free energy of transfer from gas to liquid for a variety of small molecules, since many of these values have been measured. Studies of the temperature dependence of these values (i.e. solvation enthalpies and entropies) can provide additional insights into the performance of implicit solvent models. Here, we show how to compute temperature derivatives of hydration free energies for the 3D-RISM integral equation approach. We have computed hydration free energies of 1123 small drug-like molecules (both neutral and charged). Temperature derivatives were also used to calculate hydration energies and entropies of 74 of these molecules (both neutral and charged) for which experimental data is available. While direct results have rather poor agreement with experiment, we have found that several previously proposed linear hydration free energy correction schemes give good agreement with experiment. These corrections also provide good agreement for hydration energies and entropies though simple extensions are required in some cases.

A Theoretical Study of the Hydration of Methane, from the Aqueous Solution to the sI Hydrate-Liquid Water-Gas Coexistence

Monte Carlo and molecular dynamics simulations were done with three recent water models TIP4P/2005 (Transferable Intermolecular Potential with 4 Points/2005), TIP4P/Ice (Transferable Intermolecular Potential with 4 Points/ Ice) and TIP4Q (Transferable Intermolecular Potential with 4 charges) combined with two models for methane: an all-atom one OPLS-AA (Optimal Parametrization for the Liquid State) and a united-atom one (UA); a correction for the C–O interaction was applied to the latter and used in a third set of simulations. The models were validated by comparison to experimental values of the free energy of hydration at 280, 300, 330 and 370 K, all under a pressure of 1 bar, and to the experimental radial distribution functions at 277, 283 and 291 K, under a pressure of 145 bar. Regardless of the combination rules used for σ_C,O, good agreement was found, except when the correction to the UA model was applied. Thus, further simulations of the sI hydrate were performed with the united-atom model to compare the thermal expansivity to the experiment. A final set of simulations was done with the UA methane model and the three water models, to study the sI hydrate-liquid water-gas coexistence at 80, 230 and 400 bar. The melting temperatures were compared to the experimental values. The results show the need to perform simulations with various different models to attain a reliable and robust molecular image of the systems of interest.

Effect of hydrostatic pressure on gas solubilization in micelles

Molecular dynamics simulations of anionic sodium decylsulfate and nonionic pentaethylene glycol monodecyl ether micelles in water have been performed to examine the impact of hydrostatic pressure on argon solubilization as a function of pressure. The potential-of-mean force between the micelles and argon demonstrates that nonpolar gases are attracted to the interiors of both micelles. The affinity of argon for micelle interiors, however, decreases with increasing pressure as a result of the comparatively higher molar volume of argon inside assemblies. We evaluate solubility enhancement coefficients, which describe the drop in the solute chemical potential as a function of the micellized surfactant concentration, to quantify the impact of micellization on gas solubilization. While argon is similarly attracted to the hydrophobic cores of both micelles, the gas is more effectively sequestered within nonionic micelles compared with anionic micelles as a result of salting out by charged head groups and accompanying counterions. The solubility enhancement coefficients of both micelles decrease with increasing pressure, reflecting the changing forces observed in the potentials-of-mean force. An analytical liquid drop model is proposed to describe the pressure dependence of argon solubilization within micelles that captures the simulation solubility enhancement coefficients after fitting an effective micelle radius for each surfactant.

Thermodynamics and the hydrophobic effect in a core-softened model and comparison with experiments

A simple and computationally inexpensive core-softened model, originally proposed by Franzese [G. Franzese, J. Mol. Liq. 136, 267 (2007)], was adopted to show that it exhibits properties of waterlike fluid and hydrophobic effect. The potential used between particles is spherically symmetric with two characteristic lengths. Thermodynamics of nonpolar solvation were modeled as an insertion of a modified Lennard-Jones particle. It was investigated how the anomalous predictions of the model as well as the nonpolar solvation compare with the experimental data for water anomalies and the temperature dependence of noble gases hydration. It was shown that the model qualitatively follows the same trends as water. The model is able to reproduce waterlike anomalous properties (density maximum, heat capacity minimum, isothermal compressibility, etc.) and hydrophobic effect (minimum solubility for nonpolar solutes near ambient conditions, increased solubility of larger noble gases, etc.). It is argued that the model yields similar results as more complex and computationally expensive models.

Thermodynamics of liquids: standard molar entropies and heat capacities of common solvents from 2PT molecular dynamics

We validate here the Two-Phase Thermodynamics (2PT) method for calculating the standard molar entropies and heat capacities of common liquids. In 2PT, the thermodynamics of the system is related to the total density of states (DoS), obtained from the Fourier Transform of the velocity autocorrelation function. For liquids this DoS is partitioned into a diffusional component modeled as diffusion of a hard sphere gas plus a solid component for which the DoS(υ) → 0 as υ→ 0 as for a Debye solid. Thermodynamic observables are obtained by integrating the DoS with the appropriate weighting functions. In the 2PT method, two parameters are extracted from the DoS self-consistently to describe diffusional contributions: the fraction of diffusional modes, f, and DoS(0). This allows 2PT to be applied consistently and without re-parameterization to simulations of arbitrary liquids. We find that the absolute entropy of the liquid can be determined accurately from a single short MD trajectory (20 ps) after the system is equilibrated, making it orders of magnitude more efficient than commonly used perturbation and umbrella sampling methods. Here, we present the predicted standard molar entropies for fifteen common solvents evaluated from molecular dynamics simulations using the AMBER, GAFF, OPLS AA/L and Dreiding II forcefields. Overall, we find that all forcefields lead to good agreement with experimental and previous theoretical values for the entropy and very good agreement in the heat capacities. These results validate 2PT as a robust and efficient method for evaluating the thermodynamics of liquid phase systems. Indeed 2PT might provide a practical scheme to improve the intermolecular terms in forcefields by comparing directly to thermodynamic properties.

Protruding interfacial OH groups and 'on-water' heterogeneous catalysis

The key aspect of the remarkable organic catalysis that is observed to occur at the organic/water phase boundary, the so-called 'on-water' catalysis (Narayan et al 2005 Angew. Chem. 44 3275), was recently proposed to be the protruding OH groups of water molecules at the interface that interact with the transition state (TS) via hydrogen bonding and lower activation barriers (Jung and Marcus 2007 J. Am. Chem. Soc. 129 5492). In particular, the cycloaddition reaction of quadricyclane (Q) with dimethyl azodicarboxylate (DMAD) on-water was calculated to be more than 100,000 times more efficient in terms of rate constant than the neat reaction. In this paper, we review and consider a related reaction of Q with dimethyl acetylenedicarboxylate, where nitrogen, a good H-bond acceptor, in DMAD is replaced by carbon, a poor H-bond acceptor. A very low rate acceleration of acetylenedicarboxylate on-water relative to the neat reaction is obtained theoretically, as compared to DMAD on-water, due to the relatively low H-bonding ability of acetylenedicarboxylate with water at the TS relative to the reactants. We suggest that there may also be an 'intrinsic steric effect' or orientational advantage in the on-water catalysis in general, and both electronic and steric effects may be in operation for the smaller on-water catalysis for the cycloaddition reaction of quadricyclane and acetylenedicarboxylate. A preliminary quantum mechanical/molecular mechanical (QM/MM) simulation including 1264 water molecules for the on-water reaction of DMAD + Q also suggests that there are indeed approximately two-four more H-bonds between the TS and the dangling OH groups than between the reactants and the surface.

Monte Carlo simulations of the hydrophobic effect in aqueous electrolyte solutions

The hydrophobic interaction between two methane molecules in salt-free and high salt-containing aqueous solutions and the structure in such solutions have been investigated using an atomistic model solved by Monte Carlo simulations. Monovalent salt representing NaCl and divalent salt with the same nonelectrostatic properties as the monovalent salt have been used to examine the influence of the valence of the salt species. In salt-free solution the effective interaction between the two methane molecules displayed a global minimum at close contact of the two methane molecules and a solvent-separated secondary minimum. In 3 and 5 M monovalent salt solution the potential of mean force became slightly more attractive, and in a 3 M divalent salt solution the attraction became considerably stronger. The structure of the aqueous solutions was determined by radial distribution functions and angular probability functions. The distortion of the native water structure increased with ion valence. The increase of the hydrophobic attraction was associated with (i) a breakdown of the tetrahedral structure formed by neighboring water molecules and of the hydrogen bonds between them and (i) the concomitant increase of the solution density.

Comparison of Thermodynamic Properties of Coarse-Grained and Atomic-Level Simulation Models

Thermodynamic data are often used to calibrate or test amomic-level (AL) force fields for molecular dynamics (MD) simulations. In contrast, the majority of coarse-grained (CG) force fields do not rely extensively on thermodynamic quantities. Recently, a CG force field for lipids, hydrocarbons, ions, and water, in which approximately four non-hydrogen atoms are mapped onto one interaction site, has been proposed and applied to study various aspects of lipid systems. To date, no extensive investigation of its capability to describe salvation thermodynamics has been undertaken. In the present study, a detailed picture of vaporization, solvation, and phase-partitioning thermodynamics for liquid hydrocarbons and water was obtained at CG and AL resolutions, in order to compare the two types or models and evaluate their ability to describe thermodynamic properties in the temperature range between 263 and 343 K. Both CG and AL models capture the experimental dependence of the thermodynamic properties on the temperature, albeit a systematically weaker dependence is found for the CG model. Moreover, deviations are found for solvation thermodynamics and for the corresponding enthalpy-entropy compensation for the CG model. Particularly water/oil repulsion seems to be overestimated. However, the results suggest that the thermodynamic properties considered should be reproducible by a CG model provided it is reparametrized on the basis of these liquid-phase properties.

Calculation of the aqueous solvation energy and entropy, as well as free energy, of simple polar solutes

With the advent of more powerful computers, the question of calculating thermodynamic quantities, such as the energy and the entropy, in solute-solvent systems is revisited. The calculation of these thermodynamic quantitites was limited in the past by their slow convergence relative to the free energy. Using molecular dynamics simulations, the energy, entropy, and free energy of solvation of NMA and CH(3)NH(2), as well as their relative values, have been determined. Three different methods (the thermodynamic perturbation method, the thermodynamic integration method, and a finite-difference method) are compared. The finite difference method gives the best results and accurate values for the entropy and energy were obtained using a reasonable amount to computer time. The results suggest that a meaningful thermodynamic description of biomolecular processes can be realized with present methods and the available computer time.

Mode-coupling study on the dynamics of hydrophobic hydration

The molecular motion of water in water-hydrophobic solute mixtures was investigated by the mode-coupling theory for molecular liquids based on the interaction-site description. When the model Lennard-Jones solute was mixed with water, both the translational and reorientational motions of solvent water become slower, in harmony with various experiments and molecular dynamics simulations. We compared the mechanism of the slowing down with that of the pressure dependence of the molecular motion of neat water [T. Yamaguchi, S.-H. Chong, and F. Hirata, J. Chem. Phys. 119, 1021 (2003)]. We found that the decrease in the solvent mobility caused by the solute can essentially be elucidated by the same mechanism: That is, the fluctuation of the number density of solvent due to the cavity formation by the solute strengthens the friction on the collective polarization through the dielectric friction mechanism: We also employed the solute molecule that is the same as solvent water except for the amount of partial charges, in order to alter the strength of the solute-solvent interaction continuously. The mobility of the solvent water was reduced both by the hydrophobic and strongly hydrophilic solutes, but it was enhanced in the intermediate case. Such a behavior was discussed in connection with the concept of positive and negative hydrations.

Mode-coupling study on the dynamics of hydrophobic hydration II: Aqueous solutions of benzene and rare gases

The dynamic properties of both the solute and solvent of the aqueous solution of benzene, xenon and neon are calculated by the mode-coupling theory for molecular liquids based on the interaction-site model. The B-coefficients of the reorientational relaxation and the translational diffusion of the solvent are evaluated from their dependence on the concentration of the solute, and the reorientational relaxation time of water within the hydration shell is estimated based on the two-state model. The reorientational relaxation times of water in the bulk and within the hydration shell, that of solute, and the translational diffusion coefficients of solute and solvent, are calculated at 0-30 degrees C. The temperature dependence of these dynamic properties is in qualitative agreement with that of NMR experiment reported by Nakahara et al. (M. Nakahara, C. Wakai, Y. Yoshimoto and N. Matubayasi, J. Phys. Chem., 1996, 100, 1345-1349, ref. 36), although the agreement of the absolute values is not so good. The B-coefficients of the reorientational relaxation times for benzene, xenon and neon solution are correlated with the hydration number and the partial molar volume of the solute. The proportionality with the latter is better than that with the former. These results support the mechanism that the retardation of the mobility of water is caused by the cavity formation of the solute, as previously suggested by us (T. Yamaguchi, T. Matsuoka and S. Koda, J. Chem. Phys., 2004, 120, 7590-7601, ref. 34), rather than the conventional one that the rigid hydration structure formed around the hydrophobic solute reduces the mobility of water.

Direct Measure of the Hydration Number of Aqueous Methane

A hydration number of 20 has been measured for aqueous CH4 in the temperature range 273-298 K on the basis of the simultaneous observation of the 13C isotropic resonance lines of CH4 in both aqueous solution and the dodecahedral cage of CH4 structure I hydrate. Additional experimental evidence and analysis reported in this paper suggest that the dominant aqueous hydration number of methane is 20.

Water's Hydrogen Bonds in the Hydrophobic Effect: A Simple Model

We propose a simple analytical model to account for water's hydrogen bonds in the hydrophobic effect. It is based on computing a mean-field partition function for a water molecule in the first solvation shell around a solute molecule. The model treats the orientational restrictions from hydrogen bonding, and utilizes quantities that can be obtained from bulk water simulations. We illustrate the principles in a 2-dimensional Mercedes-Benz-like model. Our model gives good predictions for the heat capacity of hydrophobic solvation, reproduces the solvation energies and entropies at different temperatures with only one fitting parameter, and accounts for the solute size dependence of the hydrophobic effect. Our model supports the view that water's hydrogen bonding propensity determines the temperature dependence of the hydrophobic effect. It explains the puzzling experimental observation that dissolving a nonpolar solute in hot water has positive entropy.

Fluctuating charge normal modes: An algorithm for implementing molecular dynamics simulations with polarizable potentials

A new method for performing molecular dynamics simulations with fluctuating charge polarizable potentials is introduced. In fluctuating charge models, polarizability is treated by allowing the partial charges to be variables, with values that are coupled to charges on the same molecule as well as those on other molecules. The charges can be efficiently propagated in a molecular dynamics simulation using extended Lagrangian dynamics. By making a coordinate change from the charge variables to a set of normal mode charge coordinates for each molecule, a new method is constructed in which the normal mode charge variables uncouple from those on the same molecule. The method is applied to the TIP4P-FQ model of water and compared to other methods for implementing the dynamics. The methods are compared using different molecular dynamics time steps.

Estimating entropies from molecular dynamics simulations

While the determination of free-energy differences by MD simulation has become a standard procedure for which many techniques have been developed, total entropies and entropy differences are still hardly ever computed. An overview of techniques to determine entropy differences is given, and the accuracy and convergence behavior of five methods based on thermodynamic integration and perturbation techniques was evaluated using liquid water as a test system. Reasonably accurate entropy differences are obtained through thermodynamic integration in which many copies of a solute are desolvated. When only one solute molecule is involved, only two methods seem to yield useful results, the calculation of solute-solvent entropy through thermodynamic integration, and the calculation of solvation entropy through the temperature derivative of the corresponding free-energy difference. One-step perturbation methods seem unsuitable to obtain entropy estimates.

Energy controlled insertion of polar molecules in dense fluids

We present a method to search low energy configurations of polar molecules in the complex potential energy surfaces associated with dense fluids. The search is done in the configurational space of the translational and rotational degrees of freedom of the molecule, combining steepest-descent and Newton-Raphson steps which embed information on the average sizes of the potential energy wells obtained from prior inspection of the liquid structure. We perform a molecular dynamics simulation of a liquid water shell which demonstrates that the method enables fast and energy-controlled water molecule insertion in aqueous environments. The algorithm finds low energy configurations of incoming water molecules around three orders of magnitude faster than direct random insertion. This method represents an important step towards dynamic simulations of open systems and it may also prove useful for energy-biased ensemble average calculations of the chemical potential.

Fast anomalous diffusion of small hydrophobic species in water

Using Car-Parrinello molecular dynamics a structural diffusion mechanism for the simplest hydrophobic species in water, an H atom, is proposed. The hydrophobic solvation cavity is a highly dynamical aggregate that actually drives, by its own hydrogen-bond fluctuations, the diffusion of the enclosed solute. This makes possible an anomalously fast diffusion that falls only short of that of "Grotthuss structural diffusion" of H+ in water. Here, the picture of a static, i.e., "iceberglike," clathrate cage is a misleading concept. The uncovered scenario is similar to the "dynamical hole mechanism" found in a very different context, that is, large molecules moving in hot polymeric melts.