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.     

Influence of water-protein hydrogen bonding on the stability of Trp-cage miniprotein. A comparison between the TIP3P and TIP4P-Ew water models

We report extensive replica exchange molecular dynamics (REMD) simulations on the folding/unfolding equilibrium of Trp-cage miniprotein using the Amber ff99SB all atom forcefield and TIP3P and TIP4P-Ew explicit water solvent models. REMD simulation-lengths in the 500 ns to the microsecond regime per replica are required to adequately sample the folding/unfolding equilibrium. We observe that this equilibrium is significantly affected by the choice of the water model. Compared with experimental data, simulations using the TIP3P solvent describe the stability of the Trp-cage quite realistically, providing a melting point which is just a few Kelvins above the experimental transition temperature of 317 K. The TIP4P-Ew model shifts the equilibrium towards the unfolded state and lowers the free energy of unfolding by about 3 kJ mol(-1) at 280 K, demonstrating the need to fine-tune the protein-forcefield depending on the chosen water model. We report evidence that the main difference between the two water models is mostly due to the different solvation of polar groups of the peptide. The unfolded state of the Trp-cage is stabilized by an increasing number of hydrogen bonds, destabilizing the α-helical part of the molecule and opening the R-D salt bridge. By reweighting the strength of solvent-peptide hydrogen bonds by adding a hydrogen bond square well potential, we can fully recover the effect of the different water models and estimate the shift in population as due to a difference in hydrogen bond-strength of about 0.4 kJ mol(-1) per hydrogen bond.     

Density dependence of the entropy and the solvation shell structure in supercritical water via molecular dynamics simulation

We perform molecular dynamics simulations of supercritical water (SCW) with a wide range of densities along a near critical isotherm using the simple point charge extended (SPC/E) pair potential in order to study the entropy and the solvation shell structure around a central water molecule. It is shown that both the translational and orientational two-particle correlation entropy terms can serve as the metrics of the translational-orientational structural orders in water and it is revealed that the translational structural order is very sensitive to the density variation in the gas-like and liquid-like region, while the orientational structural order is much more dependent upon compression in the medium-density SCW region. The comparison of the magnitudes of the full thermodynamic excess entropy and two-particle correlation entropy confirms the recent findings that the many-body terms other than two-body ones also present significant and non-neglectable contributions to the full excess entropy for the highly anomalous fluids like water. The analysis of entropy terms as a function of intermolecular distance and the orientational distribution functions as well as the three-dimensional spatial distribution functions indicate that the structural order occurs only in a much more diffused first solvation shell due to the elongated hydrogen bonds under supercritical conditions. It is revealed that no obvious second or higher neighbor shells occur in SCW, in contrast with the feature of normal liquid water that the anomalous decrease of translational order upon compression occurs mainly in the second shell.     

Free energy of liquid water from a computer simulation via cell theory

A method to calculate the free energy of water from computer simulation is presented. Based on cell theory, it approximates the potential energy surface sampled in the simulation by an anisotropic six-dimensional harmonic potential to model the three hindered translations and three hindered rotations of a single rigid water molecule. The potential is parametrized from the magnitude of the forces and torques measured in the simulation. The entropy of these six harmonic oscillators is calculated and summed with a conformational term to give the total entropy. Combining this with the simulation enthalpy yields the free energy. The six water models examined are TIP3P, SPC, TIP4P, SPC/E, TIP5P, and TIP4P-Ew. The results reproduce experiment well: free energies for all models are within 1.6 kJ mol(-1) and entropies are within 3.6 J K(-1) mol(-1). Approximately two-thirds of the entropy comes from translation, a third from rotation, and 5% from conformation. Vibrational frequencies match those in the experimental infrared spectrum and assist in their assignment. Intermolecular quantum effects are found to be small, with free energies for the classical oscillator lying 0.5-0.7 kJ mol(-1) higher than in the quantum case. Molecular displacements and vibrational and zero point energies are also calculated. Altogether, these results validate the harmonic oscillator as a quantitative model for the liquid state.     

Assessment of biomolecular force fields for molecular dynamics simulations in a protein crystal

Different biomolecular force fields (OPLS‐AA, AMBER03, and GROMOS96) in conjunction with SPC, SPC/E and TIP3P water models are assessed for molecular dynamics simulations in a tetragonal lysozyme crystal. The root mean square deviations for the C_a atoms of lysozymes are about 0.1 to 0.2 nm from OPLS‐AA and AMBER03, smaller than 0.4 nm from GROMOS96. All force fields exhibit similar pattern in B‐factors, whereas OPLS‐AA and AMBER03 accurately reproduce experimental measurements. Despite slight variations, the primary secondary structures are well conserved using different force fields. Water diffusion in the crystal is approximately ten‐fold slower than in bulk phase. The directional and average water diffusivities from OPLS‐AA and AMBER03 along with SPC/E model match fairly well with experimental data. Compared to GROMOS96, OPLS‐AA and AMBER03 predict larger hydrophilic solvent‐accessible surface area of lysozyme, more hydrogen bonds between lysozyme and water, and higher percentage of water in hydration shell. SPC, SPC/E and TIP3P water models have similar performance in most energetic and structural properties, but SPC/E outperforms in water diffusion. While all force fields overestimate the mobility and electrical conductivity of NaCl, a combination of OPLS‐AA for lysozyme and the Kirkwood‐Buff model for ions is superior to others. As attributed to the steric restraints and surface interactions, the mobility and conductivity in the crystal are reduced by one to two orders of magnitude from aqueous solution. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Correlations in liquid water for the TIP3P-Ewald, TIP4P-2005, TIP5P-Ewald, and SWM4-NDP models

Water is one of the simplest molecules in existence, but also one of the most important in biological and engineered systems. However, understanding the structure and dynamics of liquid water remains a major scientific challenge. Molecular dynamics simulations of liquid water were performed using the water models TIP3P-Ewald, TIP4P-2005, TIP5P-Ewald, and SWM4-NDP to calculate the radial distribution functions (RDFs), the relative angular distributions, and the excess enthalpies, entropies, and free energies. In addition, lower-order approximations to the entropy were considered, identifying the fourth-order approximation as an excellent estimate of the full entropy. The second-order and third-order approximations are ~20% larger and smaller than the true entropy, respectively. All four models perform very well in predicting the radial distribution functions, with the TIP5P-Ewald model providing the best match to the experimental data. The models also perform well in predicting the excess entropy, enthalpy, and free energy of liquid water. The TIP4P-2005 and SWM4-NDP models are more accurate than the TIP3P-Ewald and TIP5P-Ewald models in this respect. However, the relative angular distribution functions of the four water models reveal notable differences. The TIP5P-Ewald model demonstrates an increased preference for water molecules to act both as tetrahedral hydrogen bond donors and acceptors, whereas the SWM4-NDP model demonstrates an increased preference for water molecules to act as planar hydrogen bond acceptors. These differences are not uncovered by analysis of the RDFs or the commonly employed tetrahedral order parameter. However, they are expected to be very important when considering water molecules around solutes and are thus a key consideration in modelling solvent entropy.     

Temperature and structural changes of water clusters in vacuum due to evaporation

This paper presents a study on evaporation of pure water clusters. Molecular dynamics simulations between 20 ns and 3 micros of clusters ranging from 125 to 4096 molecules in vacuum were performed. Three different models (SPC, TIP4P, and TIP5P) were used to simulate water, starting at temperatures of 250, 275, and 300 K. We monitored the temperature, the number of hydrogen bonds, the tetrahedral order, the evaporation, the radial distribution functions, and the diffusion coefficients. The three models behave very similarly as far as temperature and evaporation are concerned. Clusters starting at a higher temperature show a higher initial evaporation rate and therefore reach the point where evaporation stop (around 240 K) sooner. The radius of the clusters is decreased by 0.16-0.22 nm after 0.5 micros (larger clusters tend to decrease their radius slightly more), which corresponds to around one evaporated molecule per nm(2). The cluster temperature seems to converge towards 215 K independent of cluster size, when starting at 275 K. We observe only small structural changes, but the clusters modeled by TIP5P show a larger percentage of molecules with low diffusion coefficient as t-->infinity, than those using the two other water models. TIP4P seems to be more structured and more hydrogen bonds are formed than in the other models as the temperature falls. The cooling rates are in good agreement with experimental results, and evaporation rates agree well with a phenomenological expression based on experimental observations.     

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.     

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.     

A molecular dynamics study of water nucleation using the TIP4P/2005 model

Extensive molecular dynamics simulations were conducted using the TIP4P/2005 water model of Abascal and Vega [J. Chem. Phys. 123, 234505 (2005)] to investigate its condensation from supersaturated vapor to liquid at 330 K. The mean first passage time method [J. Wedekind, R. Strey, and D. Reguera, J. Chem. Phys. 126, 134103 (2007); L. S. Bartell and D. T. Wu, 125, 194503 (2006)] was used to analyze the influence of finite size effects, thermostats, and charged species on the nucleation dynamics. We find that the Nose?-Hoover thermostat and the one proposed by Bussi et al. [J. Chem. Phys. 126, 014101 (2007)] give essentially the same averages. We identify the maximum thermostat coupling time to guarantee proper thermostating for these simulations. The presence of charged species has a dramatic impact on the dynamics, inducing a marked change towards a pure growth regime, which highlights the importance of ions in the formation of liquid droplets in the atmosphere. It was found a small but noticeable sign preference at intermediate cluster sizes (between 5 and 30 water molecules) corresponding mostly to the formation of the second solvation shell around the ion. The TIP4P/2005 water model predicts that anions induce faster formation of water clusters than cations of the same magnitude of charge.     

Computational study of structural and dynamical properties of formamide-water mixtures

A molecular dynamics simulation study of structural and dynamical properties in liquid mixtures of formamide and water is presented. Site-site radial pair distribution functions, local mole fractions, pair energy distributions, and tetrahedral orientational order are the quantities analyzed to investigate the local structure in the simulated mixtures, along with a review of the intermolecular structure in terms of the distribution of hydrogen bonds. Our results indicate that there is a substitution of formamide molecules by water in the hydrogen bonds and a formation of a common hydrogen bond network. By analyzing the extent of tetrahedral order in the liquid as a function of composition, it is observed that whereas the tetrahedral network of liquid water is progressively lost by increasing the formamide concentration, the water structure within the first coordination shell is preserved and somewhat enhanced. The hydrogen-bond mean lifetimes were estimated by performing a time integration of the autocorrelation functions of bond occupation numbers. The lifetimes associated with hydrogen bonds between water, formamide, and interspecies pairs are found to increase with increasing formamide concentration. The lifetimes of the water hydrogen bonds show the largest variations, supporting the picture of an enhancement of the water structure among the nearest neighbors within the first coordination shell. We have used two different force field models for water, SPC/E [J. C. Berendsen et al., J. Phys. Chem. 91, 6269 (1987)] and TIP4P/2005 [J. L. F. Abascal and C. Vega, J. Chem. Phys. 123, 234505 (2005)]. Our results for structural and dynamical properties yield very small differences between those models, the TIP4P/2005 predicting a slightly more structured liquid and, consequently, exhibiting a slightly slower translational and librational dynamics.     

Effects of ion concentration on the hydrogen bonded structure of water in the vicinity of ions in aqueous NaCl solutions

Molecular dynamics simulations of dilute and concentrated aqueous NaCl solutions are carried out to investigate the changes of the hydrogen bonded structures in the vicinity of ions for different ion concentrations. An analysis of the hydrogen bond population in the first and second solvation shells of the ions and in the bulk water is done. Although essentially no effect of ions on the hydrogen bonding is observed beyond the first solvation shell of the ions for the dilute solutions, for the concentrated solutions a noticeable change in the average number of water-water hydrogen bonds is observed in the second solvation shells of the ions and even beyond. However, the changes in the average number of hydrogen bonds are found to be relatively less when both water-water and ion-water hydrogen bonds are counted. Thus, the changes in the total number of hydrogen bonds per water are not very dramatic beyond the first solvation shell even for concentrated solutions.     

Two-particle entropy and structural ordering in liquid water

Entropies of simple point charge (SPC) water were calculated over the temperature range 278-363 K using the two-particle correlation function approximation. Then, the total two-particle contribution to the entropy of the system was divided into three parts, which we call translational, configurational, and orientational. The configurational term describes the contribution to entropy, which originates from spatial distribution of surrounding water molecules (treated as points, represented by the center of mass) around the central one. It has been shown that this term can serve as the metric of the overall orientational ordering in liquid water. Analyzing each of these three terms as a function of intermolecular distance, r, we also find a rational definition of the hydration shell around the water molecule; the estimated radii of the first and second hydration shells are 0.35 nm and 0.58 nm, respectively. We find, moreover, that the first hydration shell around the water molecule participates roughly in 70% of the total orientational entropy of water, and this rate is roughly temperature independent.     

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.     

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.     

The importance of polarizability in the modeling of solubility: quantifying the effect of solute polarizability on the solubility of small nonpolar solutes in popular models of water

In recent work by Paschek [J. Chem. Phys. 120, 6674 (2004)] and others [see H. Docherty et al., J. Chem. Phys. 125, 074510 (2006) for a review] it has been suggested that, when coupled to a simple Lennard-Jones model for various small nonpolar solute molecules, the most common models of water (e.g., SPC/E and TIP4P) fail to reproduce quantitatively the solubility of small nonpolar solute molecules in water due in part to failing to account for polarization of the solute molecule. Given the importance of such systems as test-case prototype models of the solubility of proteins and biomolecules, in this work, we investigate the impact of using a polarizable solute model with the SPC/E, TIP3P, TIP4P, TIP4P-Ew, and TIP4P/2005 rigid water models. Specifically we consider Ne, Ar, Kr, Xe, and methane as solutes. In all cases we observe that the use of a polarizable solute improves agreement between experiment and simulations, with the best agreement seen for the largest solutes, Kr, CH(4), and Xe and the modern reparametrizations of the TIP4P model, i.e., the TIP4P-Ew and TIP4P/2005 models.     

Interpretation of the water surface vibrational sum-frequency spectrum

We propose a novel interpretation of the water liquid-vapor interface vibrational sum-frequency (VSF) spectrum in terms of hydrogen-bonding classes. Unlike an absorption spectrum, the VSF signal can be considered as a sum of signed contributions from different hydrogen-bonded species in the sample. We show that the recently observed positive feature at low frequency, in the imaginary part of the signal, is a result of cancellation between the positive contributions from four-hydrogen-bonded molecules and negative contributions from those molecules with one or two broken hydrogen bonds. Spectral densities for each of these subgroups span the entire relevant spectral range. Three-body interactions within our newly developed E3B water simulation model prove to be critical in describing the proper balance between different hydrogen-bonded species, as (two-body) SPC/E, TIP4P, and TIP4P/2005 models fail to reproduce the positive feature. The results clarify the molecular origin of the VSF signal, and highlight the importance of many-body interactions for water in heterogeneous situations.     

Structural and dynamic properties of water within the solvation layer around various conformations of the glycine-based polypeptide

Several conformations of the solvated glycine-based polypeptides were investigated using molecular dynamics simulations. Some properties of water in the neighboring space around these molecules were investigated. It was found that water forms a well-defined layer-the first solvation shell-around the peptide molecule, and thickness of this layer is independent of the peptide structure and is equal to approximately 0.28 nm. Within this layer, water molecules show marked orientations relative to a peptide surface. Using the two-particle contribution to entropy as a measure of structural ordering of water, we found that the first solvation shell contributes 95% or more to the total water ordering around the peptide molecule. In investigating the dynamic properties of water, diffusion coefficients and lifetime of the hydrogen bond, clear differences between solvation layer and the bulk water were observed. It was found that the translational diffusion coefficient, D(T), decreases by 30% or more compared to bulk water; also, the lifetime of the water-water hydrogen bond clearly increases. The rotational diffusion coefficient, however, decreases only slightly, no more than approximately 10%. These differences correspond to the slightly higher energy of the hydrogen bond, and to its slightly distorted geometry. Analyzing the translational dynamics of water in the vicinity of the peptide molecule, it was deduced that the structure of the first solvation shell becomes more rigid than the structure of the bulk water. Investigation of a "pure hydrophobic" form of the polypeptide shows that the structure and the properties of water within the solvation shell are predominantly determined by the hydrophobic effect. The specific interactions between water molecules and various charge groups of the peptide molecule modifies this effect only slightly.