Common Water Models Research Articles

Entropy of Liquid Water as Predicted by the Two-Phase Thermodynamic Model and Data-Driven Many-Body Potentials.

We investigate the entropy of liquid water at ambient conditions using the two-phase thermodynamic (2PT) model, applied to both common pairwise-additive water models and the MB-pol and MB-pol(2023) data-driven many-body potentials. Our simulations demonstrate that the 2PT model yields entropy values in semiquantitative agreement with experimental data when using MB-pol and MB-pol(2023). Additionally, our analyses indicate that the entropy values predicted by pairwise-additive water models may benefit from error compensation between the inherent approximations of the 2PT model and the known limitations of these water models in describing many-body interactions. Despite its approximate nature, the simplicity of the 2PT model makes it a valuable tool for estimating relative entropy changes of liquid water across various environments, especially when combined with water models that provide a consistently robust representation of the structural, thermodynamic, and dynamical properties of liquid water.

Molecular dynamics calculations of the enthalpy of vaporization for different water models

Phase-change heat transfer processes involving water play a key role in numerous societal applications, such as electrical power generation, microelectronics cooling, or water desalination. Advanced computational tools, such as molecular dynamics (MD), are increasingly being used to simulate these processes, especially on nanostructured surfaces and in nanoparticle suspensions, to gain insight into the unique nanoscale physics at play. In these studies, it is critical to use models of the water molecule that yield accurate predictions of water’s thermophysical properties, particularly its enthalpy of vaporization. Here, we use MD to determine the enthalpy of vaporization (hfg) for six common water models and four different values of the cutoff radius (i.e., the intermolecular distance beyond which interactions are simplified or ignored). The most accurate prediction of hfg (within 2 % of the value tabulated by ASHRAE) comes from the TIP3P water model with a cutoff radius of 14 Å. We also conducted a preliminary study of the variation of hfg with pressure and found reasonable agreement between simulations and tabulated ASHRAE data (<7% error). This work provides physical insights into the sensitivity of the predicted macroscopic thermophysical properties of water to minute changes in the atomic-level properties and is expected to inform future MD-based studies of phase change processes of water, especially as they occur at elevated pressures and on nanostructured surfaces and in nanoconfined geometries.

Dependence of amino-acid dielectric relaxation on solute-water interaction: Molecular dynamics study

Computational molecular methods represent an ultimate microscope to visualize and quantify dynamics of biomolecules and nanostructures on the atomic scale and also a tool to link this dynamics to macroscopic observables.However, these potentially powerful computational methods have been so far poorly explored to gain a deeper knowledge of broadband biomolecular dielectric properties which are essential for novel non-invasive electromagnetic methods in biomedical diagnostics and treatment.To fill in this gap, we explore here a broad range of parameter space of molecular dynamics simulation (two common water models, four force fields, eight types of amino acids) and demonstrate their effects on dielectric increment and relaxation time. We found that the force field and water parameters influence the extracted dielectric increment and relaxation time of the biomolecule-water solution dielectric spectra significantly. To understand this observation, we dissected the effect of individual parameters of molecular force field on the dielectric spectra. For the first time, we showed the charges on the atoms, and the bond length of amino acids are the determining factors of both their dielectric increment and relaxation time in molecular dynamics simulations.Hence, our results leverage atomic resolution-based computational methods for a deeper understanding of biomolecular microwave dielectric properties with potential biological, medical and nanotechnological applications.

Water Models in Molecular Dynamics Simulation Prediction of Dielectric Properties of Biomaterials

To develop and reliably use diagnostic and therapeutic methods employing microwaves, we need to have an accurate knowledge of biological dielectric properties. Traditionally, dielectric properties of biosamples are determined experimentally. However, such measurements require dedicated hardware and physical availability of sufficient volume of biological samples of interest. Instead, here, we demonstrate the prediction of complex permittivity of a simple biomolecular sample using computational molecular dynamics simulations. We focus here on the role of a molecular model of water since it is the major compound determining microwave dielectric properties of biological tissues and wet samples. Here, for the first time, we analyze how the common molecular water models (SPC/E, TIP3P, and TIP4P) affect complex permittivity of biomolecular solutions predicted by molecular dynamics simulations. We found that the type of the molecular water model used in the simulation affects not only water contribution but also biomolecule contribution to the permittivity spectra. Our results contribute to in silico prediction and understanding of dielectric properties of biomaterials.

A reliable approach for calculating thermophysical properties of liquid using molecular dynamics simulations

This paper presents a new data processing method for calculating more reliable thermophysical properties of liquid using Green-Kubo Method based molecular dynamics (MD) simulation. In this study, MD simulations of water are first performed using three common water models. A new approach for analysing the simulation data is then developed to obtain statistically meaningful thermophysical properties such as thermal conductivity and viscosity. It is demonstrated that for a given desired standard deviation for a specific simulation, a suitable simulation box size can be determined. The proposed approach can assure the repeatability and reliability of the calculated thermophysical properties of liquid.

A corresponding-states analysis of the liquid-vapor equilibrium properties of common water models.

Many atomistic potential models have been proposed to reproduce the properties of real water and to capture as many of its anomalies as possible. The large number of different models indicates that this task is by no means an easy one. Some models are reasonably successful for various properties, while others are designed to account for only a very few specific features of water accurately. Among the most popular models are SPC/E, TIP4P, TIP4P/2005, TIP4P/Ice, and TIP5P-E. Here, we report the equilibrium properties of the liquid-vapor coexistence, such as the densities of the liquid phase and the vapor phase, the interfacial tension between them, and the vapor pressure at saturation. From these data, the critical parameters are determined and subsequently used to cast the liquid-vapor coexistence properties into a corresponding-states form following Guggenheim's suggestions. Doing so reveals that the three TIP4P-based models display the same corresponding-states behavior and that the SPC/E model behaves quite similarly. Only the TIP5P-E model shows clear deviations from the corresponding-states properties of the other models. A comparison with data for real water shows that the reduced surface tension is well described, while the reduced coexistence curve is too wide. The models underestimate the critical compressibility factor and overestimate Guggenheim's ratio as well as the reduced boiling temperature (Guldberg's ratio). As demonstrated by the collapse of the data for the TIP4P-based models, these deviations are inherent to the specific model and cannot be corrected by a simple reparametrization. For comparison, the results for two recent polarizable models, HBP and BK3, are shown, and both models are seen to perform well in terms of absolute numbers and in a corresponding-states framework. The kind of analysis applied here can therefore be used as a guideline in the design of more accurate and yet simple multi-purpose models of water.

Protein-Ligand Electrostatic Binding Free Energies from Explicit and Implicit Solvation.

Accurate yet efficient computational models of solvent environment are central for most calculations that rely on atomistic modeling, such as prediction of protein-ligand binding affinities. In this study, we evaluate the accuracy of a recently developed generalized Born implicit solvent model, GBNSR6 (Aguilar et al. J. Chem. Theory Comput. 2010, 6, 3613-3639), in estimating the electrostatic solvation free energies (ΔG(pol)) and binding free energies (ΔΔG(pol)) for small protein-ligand complexes. We also compare estimates based on three different explicit solvent models (TIP3P, TIP4PEw, and OPC). The two main findings are as follows. First, the deviation (RMSD = 7.04 kcal/mol) of GBNSR6 binding affinities from commonly used TIP3P reference values is comparable to the deviations between explicit models themselves, e.g. TIP4PEw vs TIP3P (RMSD = 5.30 kcal/mol). A simple uniform adjustment of the atomic radii by a single scaling factor reduces the RMS deviation of GBNSR6 from TIP3P to within the above "error margin" - differences between ΔΔG(pol) estimated by different common explicit solvent models. The simple radii scaling virtually eliminates the systematic deviation (ΔΔG(pol)) between GBNSR6 and two out of the three explicit water models and significantly reduces the deviation from the third explicit model. Second, the differences between electrostatic binding energy estimates from different explicit models is disturbingly large; for example, the deviation between TIP4PEw and TIP3P estimates of ΔΔG(pol) values can be up to ∼50% or ∼9 kcal/mol, which is significantly larger than the "chemical accuracy" goal of ∼1 kcal/mol. The absolute ΔG(pol) calculated with different explicit models could differ by tens of kcal/mol. These discrepancies point to unacceptably high sensitivity of binding affinity estimates to the choice of common explicit water models. The absence of a clear "gold standard" among these models strengthens the case for the use of accurate implicit solvation models for binding energetics, which may be orders of magnitude faster.

Mean ionic activity coefficients in aqueous NaCl solutions from molecular dynamics simulations.

The mean ionic activity coefficients of aqueous NaCl solutions of varying concentrations at 298.15 K and 1 bar have been obtained from molecular dynamics simulations by gradually turning on the interactions of an ion pair inserted into the solution. Several common non-polarizable water and ion models have been used in the simulations. Gibbs-Duhem equation calculations of the thermodynamic activity of water are used to confirm the thermodynamic consistency of the mean ionic activity coefficients. While the majority of model combinations predict the correct trends in mean ionic activity coefficients, they overestimate their values at high salt concentrations. The solubility predictions also suffer from inaccuracies, with all models underpredicting the experimental values, some by large factors. These results point to the need for further ion and water model development.

The water supercooled regime as described by four common water models

The temperature scale of simple water models in general does not coincide with the natural one. Therefore, in order to make a meaningful evaluation of different water models, a temperature rescaling is necessary. In this paper, we introduce a rescaling using the melting temperature and the temperature corresponding to the maximum of the heat capacity to evaluate four common water models (TIP4P-Ew, TIP4P-2005, TIP5P-Ew and Six-Sites) in the supercooled regime. Although all the models show the same general qualitative behavior, the TIP5P-Ew appears as the best representation of the supercooled regime when the rescaled temperature is used. We also analyze, using thermodynamic arguments, the critical nucleus size for ice growth. Finally, we speculate on the possible reasons why atomistic models do not usually crystalize while the coarse grained mW model do crystallize.

Characteristics of thermal conductivity in classical water models

The thermal conductivities of common water models are compared using equilibrium (EMD) and non-equilibrium molecular dynamics (NEMD) simulation. A complete accounting for electrostatic contributions to the heat flux was found to resolve the previously reported differing results of NEMD and EMD Green-Kubo measurements for the extended simple point-charge (SPC/E) model. Accordingly, we demonstrate the influence of long-range electrostatics on the thermal conductivity with a simple coulomb cutoff, Ewald summation, and by an extended particle-particle particle-mesh method. For each water model, the thermal conductivity is computed and decomposed in terms of frequency-dependent thermodynamic and topological contributions. The rigid, three-site SPC, SPC/E, and transferable intermolecular potential (TIP3P-Ew) water models are shown to have similar thermal conductivity values at standard conditions, whereas models that include bond stretching and angle bending have higher thermal conductivities.

Dynamics and energetics of hydrophobically confined water

The effects of water confined in regions between self-assembling entities is relevant to numerous contexts such as macromolecular association, protein folding, protein-ligand association, and nanomaterials self-assembly. Thus assessing the impact of confined water, and the ability of current modeling techniques to capture the salient features of confined water is important and timely. We present molecular dynamics simulation results investigating the effect of confined water on qualitative features of potentials of mean force describing the free energetics of self-assembly of large planar hydrophobic plates. We consider several common explicit water models including the TIP3P, TIP4P, SPC/E, TIP4P-FQ, and SWM4-NDP, the latter two being polarizable models. Examination of the free energies for filling and unfilling the volume confined between the two plates (both in the context of average number of confined water molecules and "depth" of occupancy) suggests TIP4P-FQ water molecules generally occupy the confined volume at separation distances larger than observed for other models under the same conditions. The connection between this tendency of TIP4P-FQ water and the lack of a pronounced barrier in the potential of mean force for plate-plate association in TIP4P-FQ water is explored by artificially, but systematically, populating the confined volume with TIP4P-FQ water at low plate-plate separation distances. When the critical separation distance [denoting the crossover from an unoccupied (dry) confined interior to a filled (wet) interior] for TIP4P-FQ is reduced by 0.5 Å using this approach, a barrier is observed; we rationalize this effect based on increased resistant forces introduced by confined water molecules at these low separations. We also consider the dynamics of water molecules in the confined region between the hydrophobes. We find that the TIP4P-FQ water model exhibits nonbulklike dynamics, with enhanced lateral diffusion relative to bulk. This is consistent with the reduced intermolecular water-water interaction indicated by a decreased molecular dipole moment in the interplate region. Analysis of velocity autocorrelation functions and associated power spectra indicate that the interplate region for TIP4P-FQ at a plate separation of 14.4 Å approaches characteristics of the pure water liquid-vapor interface. This is in stark contrast to the other water models (including the polarizable SWM4-NDP model).

A Monte Carlo Simulation Study of Methane Clathrate Hydrates Confined in Slit-Shaped Pores

Monte Carlo simulations are used to study the structure, stability, and dissociation mechanisms of methane hydrate crystals inside carbon-like slit-shaped pores. The simulation conditions used mimic experimental studies of the dissociation of methane and propane hydrates in mesoporous silica gels (Handa, Y. P.; Stupin, D. J. Phys. Chem. 1992, 96, 8599). Simulations are performed under conditions of fixed methane pressure and fixed water loading, with the temperature increased in steps, with long equilibrations at each temperature. The initial structures of the confined hydrates are taken to be bulk-like, and pore widths chosen to accommodate integer or half-integer numbers of hydrate unit cells. Density profiles and orientational order parameter profiles are obtained and used to understand the structural changes associated with hydrate dissociation. Three different common water models, SPC/E, TIP4P, and TIP4P/2005, are used and the results compared. For water modeled using either the TIP4P or TIP4P/2005 potentials, dissociation temperatures are depressed proportionally to the inverse pore width, as predicted by the macroscopic Gibbs-Thomson equation. This behavior is observed for pores small enough that only half-cages of the clathrate structure are present. Experimental work has verified Gibbs-Thomson behavior for pores as small as 2 nm (Seshadri, K.; Wilder, J. W.; Smith, D. H. J. Phys. Chem. B 2001, 105, 2627); micropores of the size studied here have not yet been studied by experiment. Interestingly, the dissociation of hydrates modeled using the SPC/E water potential does not display the predicted pore-size dependence, and the dissociation mechanisms in this model seem to be quite different than those in the TIP4P-type models. In the SPC/E hydrates, with increasing temperature, cage dissocation occurs before methane desorption. In TIP4P-type hydrates, these processes occur either at the same temperature (to within the resolution of this study) or with dissociation occurring at higher temperatures than desorption. These simulations show that a variety of interesting clathrate structures and phase behaviors may be accessed in suitably designed microporous materials, with potentially useful applications in gas storage or separations.

Molecular dynamics simulations of vapor/liquid coexistence using the nonpolarizable water models

The surface tension, vapor-liquid equilibrium densities, and equilibrium pressure for common water models were calculated using molecular dynamics simulations over temperatures ranging from the melting to the critical points. The TIP4P/2005 and TIP4P-i models produced better values for the surface tension than the other water models. We also examined the correlation of the data to scaling temperatures based on the critical and melting temperatures. The reduced temperature (T/T(c)) gives consistent equilibrium densities and pressure, and the shifted temperature T + (T(c, exp) - T(c, sim)) gives consistent surface tension among all models considered in this study. The modified fixed charge model which has the same Lennard-Jones parameters as the TIP4P-FQ model but uses an adjustable molecular dipole moment is also simulated to find the differences in the vapor-liquid coexistence properties between fixed and fluctuating charge models. The TIP4P-FQ model (2.72 Debye) gives the best estimate of the experimental surface tension. The equilibrium vapor density and pressure are unaffected by changes in the dipole moment as well as the surface tension and liquid density.

Protein Simulations with an Optimized Water Model: Cooperative Helix Formation and Temperature-Induced Unfolded State Collapse

A recognized shortcoming in current protein simulations is that most force fields are parametrized with relatively primitive three-site water models. Since the deficiencies of the common three-site water models in reproducing the phase diagram of water are well-known, an improved description of the solvent will be required, for example, to study proteins in molecular simulations at thermodynamic conditions other than standard temperature and pressure. Here, we combine a protein force field derived from Amber ff03 together with the highly optimized TIP4P/2005 water model, with a small backbone modification to match the population of helical states obtained with the new water model to experiment. Remarkably, we find that the resulting force field, Amber ff03w, produces a more cooperative helix-coil transition, compared with the similarly "backbone-corrected" Amber ff03* model with TIP3P water, with calculated helix propagation parameters in good agreement with the experiment. The radius of gyration for nonhelical conformations is significantly larger for Amber ff03w than for Amber ff03* and shows a collapse with increasing temperature as found in single-molecule experiments on longer proteins. The origin of the collapse appears to be a more favorable enthalpic component of the peptide-solvent interaction and is correlated with increasing turn formation, in accord with the experiment. In addition to this enhanced cooperativity, we verify that, with the new force field, replica exchange folding simulations of the GB1 hairpin and Trp cage result in folded structures, starting from completely unfolded initial conditions; simulations of folded proteins are also stable. These results together suggest that Amber ff03w (with TIP4P/2005) will be well suited for studying protein folding and properties of unfolded state and intrinsically disordered proteins over a wide range of thermodynamic conditions.

The shear viscosities of common water models by non-equilibrium molecular dynamics simulations

We have performed non-equilibrium molecular dynamics simulations to simulate the shear viscosity of five commonly used water models, SPC, SPC/E, TIP3P, TIP4P and TIP5P, using the periodic perturbation method. By comparing the computed viscosities of the water models, we find that the viscosities of SPC/E and TIP5P are in better agreement with the experimental value than the others. Furthermore, we have systematically investigated the effects of system size, acceleration amplitude, electrostatic treatment and cut-off length on the viscosities of the SPC/E and TIP5P models. Finally, by simulating the viscosity of water as a function of temperature, our work suggests that the SPC/E and TIP5P models reasonably reproduce the water viscosity vs. temperature experimental profile from 283 to 373 K at 1 bar, although SPC/E gives a better viscosity at high temperatures and TIP5P reproduces closer viscosity values at low temperatures.

Assessing the thermodynamic signatures of hydrophobic hydration for several common water models

Following the conclusions of an information theory analysis that hydrophobic hydration is dictated by the equation of state of liquid water, we perform simulations of ten different water models to examine the correlation between the fidelity of each model to the experimental density of liquid water and the accuracy of its description of methane hydration. We find that the three- and five-point water models provide an inferior description of both the liquid density and methane solubility compared to the four-point water models. Of the four-point water models, TIP4P/2005 provides the best description of both the aqueous equation-of-state and methane hydration thermodynamics. When the optimized potentials for liquid simulation united-atom description for methane is used, we find that while the entropy and heat capacity of methane hydration are in excellent agreement with experiment, the chemical potential and enthalpy are systematically shifted upwards. We subsequently reoptimize the methane interaction to accurately reproduce the experimental solubilities as a function of temperature by accounting for missing attractive interactions.

Problems with Some Current Water Models for Close Pair Interactions That Are Not Near the Minimum Energy.

The ability of an empirical, polarizable model of water to predict a thermal ensemble of molecular configurations at ambient conditions was examined using first-principle quantum mechanics. The empirical model of water selected for this evaluation was the TTM2-F model. The quantum mechanical methodology selected was the second-order Møller-Plesset model (MP2). Only pairwise interaction energies were considered. Significant deviations from the empirical model were found. Similar results were found for ad-hoc comparisons with several other common water models including the TIP3P, TIP4P, TIP4P-FQ, TIP5P, TTM2.1-F, TTM2.2-F, TTM3-F, and POL5/QZ potential models. Our results show that spatially close dimer configurations with interaction energies notably above the potential well minimum (but are still thermally accessible at ambient conditions) are the source of the largest deviations. To assist others in future water model parametrizations we report the MP2 near complete basis set limit energies for 840 water dimer configurations sampled from an approximate thermal ensemble at ambient conditions.

Charge Asymmetries in Hydration of Polar Solutes

We study the solvation of polar molecules in water. The center of water's dipole moment is offset from its steric center. In common water models, the Lennard-Jones center is closer to the negatively charged oxygen than to the positively charged hydrogens. This asymmetry of water's charge sites leads to different hydration free energies of positive versus negative ions of the same size. Here, we explore these hydration effects for some hypothetical neutral solutes, and two real solutes, with molecular dynamics simulations using several different water models. We find that, like ions, polar solutes are solvated differently in water depending on the sign of the partial charges. Solutes having a large negative charge balancing diffuse positive charges are preferentially solvated relative to those having a large positive charge balancing diffuse negative charges. Asymmetries in hydration free energies can be as large as 10 kcal/mol for neutral benzene-sized solutes. These asymmetries are mainly enthalpic, arising primarily from the first solvation shell water structure. Such effects are not readily captured by implicit solvent models, which respond symmetrically with respect to charge.

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.

“Tetrahedrality” and the Relationship between Collective Structure and Radial Distribution Functions in Liquid Water

Molecular dynamics simulations of pure liquid water under ambient conditions using four common empirical water models have been analyzed to determine how well the oxygen-oxygen radial distribution function, g(OO)(r), used as the sole criterion of congruence with experiment, captures variations in the actual anisotropic collective structuring for these models. The largest systematic deviations from tetrahedrality were found to be due to deformations of the angle between the two closet hydrogen bond donor neighbors, but for intrinsic geometric reasons, these were found to contribute less to g(OO)(r) than deformations of the angles between one hydrogen bond donor neighbor and one hydrogen bond acceptor neighbor. Relying exclusively on a qualitative characterization of the second peak in g(OO)(r) seems to overemphasize the differences between the structuring in some of these models.