Heat capacity effects associated with the hydrophobic hydration and interaction of simple solutes: a detailed structural and energetical analysis based on molecular dynamics simulations

We examine the SPCE [H. J. C. Berendsen et al., J. Chem. Phys. 91, 6269 (1987)] and TIP5P [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys 112, 8910 (2000)] water models using a temperature series of molecular-dynamics simulations in order to study heat-capacity effects associated with the hydrophobic hydration and interaction of xenon particles. The temperature interval between 275 and 375 K along the 0.1-MPa isobar is studied. For all investigated models and state points we calculate the excess chemical potential for xenon employing the Widom particle insertion technique. The solvation enthalpy and excess heat capacity is obtained from the temperature dependence of the chemical potentials and, alternatively, directly by Ewald summation, as well as a reaction field based method. All three methods provide consistent results. In addition, the reaction field technique allows a separation of the solvation enthalpy into solute/solvent and solvent/solvent parts. We find that the solvent/solvent contribution to the excess heat capacity is dominating, being about one order of magnitude larger than the solute/solvent part. This observation is attributed to the enlarged heat capacity of the water molecules in the hydration shell. A detailed spatial analysis of the heat capacity of the water molecules around a pair of xenon particles at different separations reveals that even more enhanced heat capacity of the water located in the bisector plane between two adjacent xenon atoms is responsible for the maximum of the heat capacity found for the desolvation barrier distance, recently reported by Shimizu and Chan [J. Am. Chem. Soc. 123, 2083 (2001)]. The about 60% enlarged heat capacity of water in the concave part of the joint xenon-xenon hydration shell is the result of a counterplay of strengthened hydrogen bonds and an enhanced breaking of hydrogen bonds with increasing temperature. Differences between the two models with respect to the heat capacity in the xenon-xenon contact state are attributed to the different water model bulk heat capacities, and to the different spatial extension of the structure effect introduced by the hydrophobic particles. Similarities between the different states of water in the joint xenon-xenon hydration shell and the properties of stretched water are discussed.     

Structural properties of water: comparison of the SPC, SPCE, TIP4P, and TIP5P models of water

Molecular-dynamics simulations were carried out for the SPC, SPCE, TIP4P, and TIP5P models of water at 298 K. From these results we determine the following quantities: the absolute entropy using the two-particle approximation, the mean lifetime of the hydrogen bond, the mean number of hydrogen bonds per molecule, and the mean energy of the hydrogen bond. From the entropy calculations we find that nearly all contributions to the total entropy originates from the orientation effects. Moreover, we determine the contributions to the total entropy which originate from the first, second, and higher solvation shells. It is interesting that the limits between solvation shells are clearly visible. The first solvation shell (0.22 < r < 0.36 nm) contributes approximately 43 J mol K to the total entropy; the second solvation shell (0.36 < r < 0.60 nm) contributes approximately 12 J mol K, while contributions from the third and other solvation shells are very small, approximately 2 J mol K in summary. This indicates that water molecules are strongly ordered up to 0.55-0.6 nm around the central water molecule, and beyond this limit the ordering diminishes. The results of calculations (entropy and hydrogen bonds) are compared with the experimental data for the choosing of the best water model. We find that the SPC and TIP4P models reproduce the best experimental values, and we recommend these models for computer simulations of the aqueous solution of biomolecules.     

Relation between the melting temperature and the temperature of maximum density for the most common models of water

Water exhibits a maximum in density at normal pressure at 4 degrees above its melting point. The reproduction of this maximum is a stringent test for potential models used commonly in simulations of water. The relation between the melting temperature and the temperature of maximum density for these potential models is unknown mainly due to our ignorance about the melting temperature of these models. Recently we have determined the melting temperature of ice I(h) for several commonly used models of water (SPC, SPC/E, TIP3P, TIP4P, TIP4P/Ew, and TIP5P). In this work we locate the temperature of maximum density for these models. In this way the relative location of the temperature of maximum density with respect to the melting temperature is established. For SPC, SPC/E, TIP3P, TIP4P, and TIP4P/Ew the maximum in density occurs at about 21-37 K above the melting temperature. In all these models the negative charge is located either on the oxygen itself or on a point along the H-O-H bisector. For the TIP5P and TIP5P-E models the maximum in density occurs at about 11 K above the melting temperature. The location of the negative charge appears as a geometrical crucial factor to the relative position of the temperature of maximum density with respect to the melting temperature.     

Hydration thermodynamic properties of amino acid analogues: a systematic comparison of biomolecular force fields and water models

We present an extensive study on hydration thermodynamic properties of analogues of 13 amino acid side chains at 298 K and 1 atm. The hydration free energies DeltaG, entropies DeltaS, enthalpies DeltaH, and heat capacities Deltac(P)() were determined for 10 combinations of force fields and water models. The statistical sampling was extended such that precisions of 0.3, 0.8, 0.8 kJ/mol and 25 J/(mol K) were reached for DeltaG, TDeltaS, DeltaH, and Deltac(P)(), respectively. The three force fields used in this study are AMBER99, GROMOS 53A6, and OPLS-AA; the five water models are SPC, SPC/E, TIP3P, TIP4P, and TIP4P-Ew. We found that the choice of water model strongly influences the accuracy of the calculated hydration entropies, enthalpies, and heat capacities, while differences in accuracy between the force fields are small. On the basis of an analysis of the hydrophobic analogues of the amino acid side chains, we discuss what properties of the water models are responsible for the observed discrepancies between computed and experimental values. The SPC/E water model performs best with all three biomolecular force fields.     

Calculation of the water-cyclohexane transfer free energies of neutral amino acid side-chain analogs using the OPLS all-atom force field

We calculated the free energy of solvation of the neutral analogs of 18 amino acid side-chains (not including glycine and proline) using the OPLS all-atom force field in TIP4P water, SPC water, and cyclohexane by molecular dynamics simulation and thermodynamic integration. The average unsigned errors in the free energies of solvation in TIP4P, SPC, and cyclohexane are 4.4, 4.9, and 2.1 kJ/mol respectively. Most of the calculated hydration free energies are not favorable enough compared to experiment. The largest errors are found for tryptophan, histidine, glutamic acid, and glutamine. The average unsigned errors in the free energy of transfer from TIP4P to cyclohexane and from SPC to cyclohexane are 4.0 and 4.1 kJ/mol, respectively. The largest errors, of more than 7.5 kJ/mol, are found for histidine, glutamine, and glutamatic acid.     

Computer simulation of two new solid phases of water: Ice XIII and ice XIV

NpT Monte Carlo simulations have been performed for two recently discovered solid phases of water which have been denoted as ice XIII and ice XIV C. G. Salzmann et al. [Science311, 1758 (2006)]. Several potential models of water were considered, namely, the traditional SPC/E, TIP4P, and TIP5P and the more recent TIP5P-E, TIP4P-Ew, TIP4P/Ice, and TIP4P/2005 models. Significant differences in density and oxygen-oxygen radial distribution functions are found between the predictions of the SPC/E, TIP5P, and the models of the TIP4P family. The models TIP4P/Ice and TIP4P/2005 provide the best estimates of the density.     

Solvation free energies of amino acid side chain analogs for common molecular mechanics water models

Quantitative free energy computation involves both using a model that is sufficiently faithful to the experimental system under study (accuracy) and establishing statistically meaningful measures of the uncertainties resulting from finite sampling (precision). In order to examine the accuracy of a range of common water models used for protein simulation for their solute/solvent properties, we calculate the free energy of hydration of 15 amino acid side chain analogs derived from the OPLS-AA parameter set with the TIP3P, TIP4P, SPC, SPC/E, TIP3P-MOD, and TIP4P-Ew water models. We achieve a high degree of statistical precision in our simulations, obtaining uncertainties for the free energy of hydration of 0.02-0.06 kcal/mol, equivalent to that obtained in experimental hydration free energy measurements of the same molecules. We find that TIP3P-MOD, a model designed to give improved free energy of hydration for methane, gives uniformly the closest match to experiment; we also find that the ability to accurately model pure water properties does not necessarily predict ability to predict solute/solvent behavior. We also evaluate the free energies of a number of novel modifications of TIP3P designed as a proof of concept that it is possible to obtain much better solute/solvent free energetic behavior without substantially negatively affecting pure water properties. We decrease the average error to zero while reducing the root mean square error below that of any of the published water models, with measured liquid water properties remaining almost constant with respect to our perturbations. This demonstrates there is still both room for improvement within current fixed-charge biomolecular force fields and significant parameter flexibility to make these improvements. Recent research in computational efficiency of free energy methods allows us to perform simulations on a local cluster that previously required large scale distributed computing, performing four times as much computational work in approximately a tenth of the computer time as a similar study a year ago.     

Liquid-liquid phase transitions in supercooled water studied by computer simulations of various water models

Liquid-liquid and liquid-vapor coexistence regions of various water models were determined by Monte Carlo (MC) simulations of isotherms of density fluctuation-restricted systems and by Gibbs ensemble MC simulations. All studied water models show multiple liquid-liquid phase transitions in the supercooled region: we observe two transitions of the TIP4P, TIP5P, and SPCE models and three transitions of the ST2 model. The location of these phase transitions with respect to the liquid-vapor coexistence curve and the glass temperature is highly sensitive to the water model and its implementation. We suggest that the apparent thermodynamic singularity of real liquid water in the supercooled region at about 228 K is caused by an approach to the spinodal of the first (lowest density) liquid-liquid phase transition. The well-known density maximum of liquid water at 277 K is related to the second liquid-liquid phase transition, which is located at positive pressures with a critical point close to the maximum. A possible order parameter and the universality class of liquid-liquid phase transitions in one-component fluids are discussed.     

The melting temperature of the most common models of water

The melting temperature of ice I(h) for several commonly used models of water (SPC, SPC/E,TIP3P,TIP4P, TIP4P/Ew, and TIP5P) is obtained from computer simulations at p = 1 bar. Since the melting temperature of ice I(h) for the TIP4P model is now known [E. Sanz, C. Vega, J. L. F. Abascal, and L. G. MacDowell, Phys. Rev. Lett. 92, 255701 (2004)], it is possible to use the Gibbs-Duhem methodology [D. Kofke, J. Chem. Phys. 98, 4149 (1993)] to evaluate the melting temperature of ice I(h) for other potential models of water. We have found that the melting temperatures of ice I(h) for SPC, SPC/E, TIP3P, TIP4P, TIP4P/Ew, and TIP5P models are T = 190 K, 215 K, 146 K, 232 K, 245 K, and 274 K, respectively. The relative stability of ice I(h) with respect to ice II for these models has also been considered. It turns out that for SPC, SPC/E, TIP3P, and TIP5P the stable phase at the normal melting point is ice II (so that ice I(h) is not a thermodynamically stable phase for these models). For TIP4P and TIP4P/Ew, ice I(h) is the stable solid phase at the standard melting point. The location of the negative charge along the H-O-H bisector appears as a critical factor in the determination of the relative stability between the I(h) and II ice forms. The methodology proposed in this paper can be used to investigate the effect upon a coexistence line due to a change in the potential parameters.     

Capillary waves at the liquid-vapor interface and the surface tension of water

Capillary waves occurring at the liquid-vapor interface of water are studied using molecular dynamics simulations. In addition, the surface tension, determined thermodynamically from the difference in the normal and tangential pressure at the liquid-vapor interface, is compared for a number of standard three- and four-point water models. We study four three-point models (SPC/E, TIP3P, TIP3P-CHARMM, and TIP3P-Ew) and two four-point models (TIP4P and TIP4P-Ew). All of the models examined underestimate the surface tension; the TIP4P-Ew model comes closest to reproducing the experimental data. The surface tension can also be determined from the amplitude of capillary waves at the liquid-vapor interface by varying the surface area of the interface. The surface tensions determined from the amplitude of the logarithmic divergence of the capillary interfacial width and from the traditional thermodynamic method agree only if the density profile is fitted to an error function instead of a hyperbolic tangent function.     

Higher-order virial coefficients of water models

We use the Mayer sampling method, with both direct and overlap sampling, to calculate and compare classical virial coefficients up to B6 for various water models (SPC, SPC/E, MSPC/E, TIP3P, and TIP4P). The precision of the computed values ranges from 0.1% for B2 to an average of 25% for B6. When expressed in a form scaled by the critical properties, the values of the coefficients for SPC water are observed to greatly exceed the magnitude of corresponding coefficients for the simple Lennard-Jones model. We examine the coefficients in the context of the equation of state and the Joule-Thomson coefficient. Comparisons of these properties are made both to established molecular simulation data for each respective model and to real water. For all models, the virial series up to B5 describes the equation of state along the saturated vapor line better than the series that includes B6. At supercritical temperatures, however, the sixth-order series often describes pressure-volume-temperature behavior better than the fifth-order series. For example, the sixth-order virial equation of state for SPC/E water predicts the 673 K isotherm within 8% of published molecular simulation values up to a density of 9 mol/L (roughly half the critical density of SPC/E water).     

Common features of orientational order at the temperature of maximum density for various water models: molecular dynamics study

Canonical ensembles for liquid water were obtained from molecular dynamics simulations at various temperatures using the TIP5P, TIP4P-FQ, TIP4P, and SPC/E water models at a fixed density of 1 g/cm3. From these ensembles, it was found that the distributions of the orientational order parameter q of these models showed similar patterns as temperature changed except that the distributions were shifted relative to each other by the difference of their temperature of maximum density (TMD). The four models exhibited similar distributions and average values of orientational order around their respective TMDs, and these common features were investigated in detail especially. The current study suggests that the unique microscopic configuration of water molecules cause TMD phenomenon in any reasonable water model. This finding provides a useful tool in the development of new water potentials by offering guidelines to predict the TMD, avoiding troublesome isothermal-isobaric ensemble simulations.     

Study of multipole contributions to the structure of water around ions in solution using the soft sticky dipole-quadrupole-octupole (SSDQO) model of water

The solvation of ions in the soft sticky dipole-quadrupole-octupole (SSDQO) model for liquid water is presented here. This new potential energy function for liquid water describes water-water interactions by a Lennard-Jones term plus a sticky potential consisting of an approximate moment expansion with point dipole, quadrupole, and octupole moments. The SSDQO potential energy function using the moments from extended simple point charge (SPC/E), TIP3P, or TIP5P reproduces the pair potential energy functions and radial distribution functions of the respective multipoint model but it is much faster than even the three-point models. Here, the solvation of ions in SSDQO water is studied using ion-water potential energy functions consisting of moment expansions up to the charge-quadrupole term, up to the charge-octupole term, and up to an approximate charge-hexadecapole term using the moments of SPC/E water. The radial distributions from Monte Carlo simulations show the best agreement with the results for ions in SPC/E water for the expansion up to the charge-hexadecapole term. Thus, the best results are obtained when the water-water and ion-water potentials are exact up to the 1r(4) term and also contain an approximate 1r(5) term. Overall, the simplicity, efficiency, and accuracy of the SSDQO potential make it potentially very useful for computer simulations of aqueous solvation.     

Simulated surface tensions of common water models

Initial simulated values of the surface tension for the SPC/E water model have indicated excellent agreement with experiment. More recently, differing values have been obtained which are significantly lower than previous estimates. Here, we attempt to explain the differences between the previous studies and show that a variety of simulation conditions can affect the final surface tension values. Consistent values for the surface tensions of six common fixed charge water models (TIP3P, SPC, SPC/E, TIP4P, TIP5P, and TIP6P) are then determined for four temperatures between 275 and 350 K. The SPC/E and TIP6P models provide the best agreement with experiment.     

Sources of the deficiencies in the popular SPC/E and TIP3P models of water

Motivated by the results of Vega et al. [J. Phys. Condens. Matter 20, 153101 (2008)] about the phase diagram of water, and by the results of Kiss and Baranyai [J. Chem. Phys. 131, 204310 (2009)] about the properties of gas-phase clusters, we carried out a comparative study of the structure modeled by SPC∕E and TIP3P interactions in ambient liquid water. The gas-phase clusters of SPC∕E and TIP3P models show erroneous structures, while TIP4P-type models, either polarizable or not, provide qualitatively correct results. The trimers of SPC∕E and TIP3P are planar in gas phase, contrary to experimental and TIP4P-type models. The aim of this study was to see whether traces of these false geometries characteristic to SPC∕E and TIP3P in gas phase can also be found in the liquid phase. For this purpose we selected trimers formed by adjacent neighbors of water molecules in the liquid and calculated their geometrical features. We determined angles formed by the HO bonds of the molecules with OO vectors and with the normal vector of the OOO plane in the selected trimers. Our results showed that, despite high temperature, the SPC∕E and TIP3P water contains larger number of planar arrangements than other TIP4P-type models. Although structural differences presented in this study are small, they are accurately detectable. These results weaken the reliability of studies obtained by the SPC∕E or TIP3P models even in the liquid phase.     

Comparison of solvation-effect methods for the simulation of peptide interactions with a hydrophobic surface

In this study we investigated the interaction behavior between thirteen different small peptides and a hydrophobic surface using three progressively more complex methods of representing solvation effects: a united-atom implicit solvation method [CHARMM 19 force field (C19) with Analytical Continuum Electrostatics (ACE)], an all-atom implicit solvation method (C22 with GBMV), and an all-atom explicit solvation method (C22 with TIP3P). The adsorption behavior of each peptide was characterized by the calculation of the potential of mean force as a function of peptide-surface separation distance. The results from the C22/TIP3P model suggest that hydrophobic peptides exhibit relatively strong adsorption behavior, polar and positively-charged peptides exhibit negligible to relatively weak favorable interactions with the surface, and negatively-charged peptides strongly resist adsorption. Compared to the TIP3P model, the ACE and GBMV implicit solvent models predict much stronger attractions for the hydrophobic peptides as well as stronger repulsions for the negatively-charged peptides on the CH(3)-SAM surface. These comparisons provide a basis from which each of these implicit solvation methods may be reparameterized to provide closer agreement with explicitly represented solvation in simulations of peptide and protein adsorption to functionalized surfaces.     

Structural properties and thermodynamics of water clusters: a Wang-Landau study

The temperature dependence of structural properties and thermodynamic behavior of water clusters has been studied using Wang-Landau sampling. Four potential models, simple point charge/extended (SPC/E), transferable intermolecular potential 3 point (TIP3P), transferable intermolecular potential 4 point (TIP4P), and Gaussian charge polarizable (GCP), are compared for ground states and properties at finite temperatures. Although the hydrogen bond energy and the distance of the nearest-neighbor oxygen pair are significantly different for TIP4P and GCP models, they approach to similar ground state structures and melting transition temperatures in cluster sizes we considered. Comparing with TIP3P, SPC/E model provides properties closer to that of TIP4P and GCP.     

Prediction of self-diffusion coefficient and shear viscosity of water and its binary mixtures with methanol and ethanol by molecular simulation

Density, self-diffusion coefficient, and shear viscosity of pure liquid water are predicted for temperatures between 280 and 373 K by molecular dynamics simulation and the Green-Kubo method. Four different rigid nonpolarizable water models are assessed: SPC, SPC/E, TIP4P, and TIP4P/2005. The pressure dependence of the self-diffusion coefficient and the shear viscosity for pure liquid water is also calculated and the anomalous behavior of these properties is qualitatively well predicted. Furthermore, transport properties as well as excess volume and excess enthalpy of aqueous binary mixtures containing methanol or ethanol, based on the SPC/E and TIP4P/2005 water models, are calculated. Under the tested conditions, the TIP4P/2005 model gives the best quantitative and qualitative agreement with experiments for the regarded transport properties. The deviations from experimental data are of 5% to 15% for pure liquid water and 5% to 20% for the water + alcohol mixtures. Moreover, the center of mass power spectrum of water as well as the investigated mixtures are analyzed and the hydrogen-bonding structure is discussed for different states.