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.     

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.     

Entropy, diffusivity, and structural order in liquids with waterlike anomalies

The excess entropy, defined as the difference between the entropies of the liquid and the ideal gas under identical density and temperature conditions, is studied as a function of density and temperature for liquid silica and a two-scale ramp potential, both of which are known to possess waterlike liquid state anomalies. The excess entropy for both systems is evaluated using a fairly accurate pair correlation approximation. The connection between the excess entropy and the density and diffusional anomalies is demonstrated. Using the pair correlation approximation to the excess entropy, it can be shown that if the energetically favorable local geometries in the low and high density limits have different symmetries, then a structurally anomalous regime can be defined in terms of orientational and translational order parameters, as in the case of silica and the two-scale ramp system but not for the one-scale ramp liquid. Within the category of liquids with waterlike anomalies, we show that the relationship between the macroscopic entropy and internal energy is sufficient to distinguish between those with local anisotropy and consequent open packings at low densities and those with isotropic interactions but multiple length scales. Since it is straightforward to evaluate the pair correlation entropy and internal energy from simulations or experimental data, such plots should provide a convenient means to diagnose the existence as well as type of anomalous behavior in a range of liquids, including ionic and intermetallic melts and complex fluids with ultrasoft repulsions.     

Molecular dynamics simulations of water droplets on polymer surfaces

Molecular dynamics simulations were used to study the wetting of polymer surfaces with water. Contact angles of water droplets on crystalline and two amorphous polyethylene (PE) and poly(vinyl chloride) (PVC) surfaces were extracted from atomistic simulations. Crystalline surfaces were produced by duplicating the unit cell of an experimental crystal structure, and amorphous surfaces by pressing the bulk polymer step by step at elevated temperature between two repulsive grid surfaces to a target density. Different-sized water droplets on the crystalline PE surface revealed a slightly positive line tension on the order of 10(-12)-10(-11) N, whereas droplets on crystalline PVC did not yield a definite line tension. Microscopic contact angles produced by the simple point charge (SPC) water model were mostly a few degrees smaller than those produced by the extended SPC model, which, as the model with lowest bulk energy, presents an upper boundary for contact angles. The macroscopic contact angle for the SPC model was 94 degrees on crystalline PVC and 113 degrees on crystalline PE. Amorphicity of the surface increased the water contact angle on PE but decreased it on PVC, for both water models. If the simulated contact angles on crystalline and amorphous surfaces are combined in proportion to the crystallinity of the polymer in question, simulated values in relatively good agreement with measured values are obtained.     

A first principles study on the structure of ice-VI: static distortion, molecular geometry, and proton ordering

We have studied the structure of ice-VI by examining all ice-rule-allowed structures in its primary unit cell of 10 water molecules with first principles methods. A significant amount of static distortions in the oxygen positions away from their crystallographic positions are found, which is in good agreements with significant higher-order terms in the atomic displacement parameters obtained from X-ray and neutron diffraction data. Structural anomalies (such as exceptionally short OH bonds and small H-O-H angles) noted in conventional crystal structure refinements were not seen in our ab initio calculations, and it is evident that these structural anomalies arose from oversimplified models in which static distortions are not properly accounted for. Our results also show that the molecular geometry of water in ice-VI is similar to but richer than those in ice-Ih and ice-VII. Larger distortions in bond lengths/angles and correlation between the molecular geometry and the neighboring environments were found. Different proton-ordering schemes proposed in the literature were examined, and our calculations provide evidence in favor of a ferroelectric phase of the proton-ordered counterpart of ice-VI at about 80 K.     

Computational investigation of order, structure, and dynamics in modified water models

Model liquids have been constructed to study the role of local structure in the anomalous properties of liquid water. The intermolecular potentials were modified by increasing the weight of the Lennard-Jones term relative to the electrostatic term in the SPC/E model for water. The resulting family of liquids varies from SPC/E water to a Lennard-Jones-like liquid. Properties were measured as a function of density and temperature. The local structure was described by two order parameters, one measuring the tetrahedral order and the other measuring the translational order. The translational order parameter was found to be large for both tetrahedral and Lennard-Jones liquids, but to go through a minimum as the potentials were modified, demonstrating that the two types of structure are incompatible. Just as in water several properties (e.g., the translational diffusion coefficient, entropy) exhibit anomalous density dependence as a result of the breakdown of local tetrahedrality, we observed nonmonotonic behavior of the translational diffusion constant and reorientational relaxation rate as the fluids were transformed from tetrahedral to Lennard-Jones-like. This is also an indication of the incompatibility between Lennard-Jones and water-like structure.     

The thickness of a liquid layer on the free surface of ice as obtained from computer simulation

Molecular dynamic simulations were performed for ice I(h) with a free surface by using four water models, SPC/E, TIP4P, TIP4P/Ice, and TIP4P/2005. The behavior of the basal plane, the primary prismatic plane, and of the secondary prismatic plane when exposed to vacuum was analyzed. We observe the formation of a thin liquid layer at the ice surface at temperatures below the melting point for all models and the three planes considered. For a given plane it was found that the thickness of a liquid layer was similar for different water models, when the comparison is made at the same undercooling with respect to the melting point of the model. The liquid layer thickness is found to increase with temperature. For a fixed temperature it was found that the thickness of the liquid layer decreases in the following order: the basal plane, the primary prismatic plane, and the secondary prismatic plane. For the TIP4P/Ice model, a model reproducing the experimental value of the melting temperature of ice, the first clear indication of the formation of a liquid layer, appears at about -100 degrees C for the basal plane, at about -80 degrees C for the primary prismatic plane, and at about -70 degrees C for the secondary prismatic plane.     

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).     

Effects of water model and simulation box size on protein diffusional motions

We performed molecular dynamics simulations of ubiquitin with the distinct water models, TIP3P, SPC, SPC/E, and SPC/Fw, in different system sizes with different box shapes. The translational diffusion constants of pure water linearly depend on the effective box length, which is known as finite size effect, whereas the first and second rotational times of pure water are nearly constant. We then observed that both the overall translation and rotational motions of the protein are linearly correlated to the viscosity of pure water. As expected from the finite size effect in pure water, the translational diffusion of the protein is significantly affected by the system size, and rotational diffusion is nearly size-independent. After correction for the finite size effect, the SPC/E and SPC/Fw models reproduce both the translational and rotational motion of the protein relatively well. Thus, water models that reproduce the experimentally derived diffusional properties of pure water accurately are expected to also be suitable for simulating protein diffusion quantitatively.     

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.     

Dynamical properties of the soft sticky dipole-quadrupole-octupole water model: a molecular dynamics study

The dynamical properties of the soft sticky dipole-quadrupole-octupole (SSDQO) water model using SPC/E moments are calculated utilizing molecular dynamics simulations. This new potential for liquid water describes the water-water interactions by a Lennard-Jones term and a sticky potential, which is an approximate moment expansion with point dipole, quadrupole, and octupole moments, and reproduces radial distribution functions of pure liquid water using the moments of SPC/E [Ichiye and Tan, J. Chem. Phys. 124, 134504 (2006)]. The forces and torques of SSDQO water for the dipole-quadrupole, quadrupole-quadrupole, and dipole-octupole interactions are derived here. The simulations are carried out at 298 K in the microcanonical ensemble employing the Ewald method for the long-range dipole-dipole interactions. Here, various dynamical properties associated with translational and rotational motions of SSDQO water using the moments of SPC/E (SSDQO:SPC/E) water are compared with the results from SPC/E and also experiment. The self-diffusion coefficient of SSDQO:SPC/E water is found to be in excellent agreement with both SPC/E and experiment whereas the single particle orientational relaxation time for dipole vector is better than SPC/E water but it is somewhat smaller than experiment. The dielectric constant of SSDQO:SPC/E is essentially identical to SPC/E, and both are slightly lower than experiment. Also, molecular dynamics simulations of the SSDQO water model are found to be about twice as fast as three-site models such as SPC/E.     

Self-diffusion of supercooled o-terphenyl near the glass transition temperature

Self-diffusion coefficients for the low molecular weight glass former o-terphenyl have been measured near Tg by isothermally desorbing thin film bilayers of deuterio and protio o-terphenyl in a vacuum chamber. We observe translational diffusion that is about 100 times faster at Tg + 3 K than the Stokes-Einstein prediction. Predictions from random first order transition theory and a dynamic facilitation approach are in reasonable agreement with our results; in these approaches, enhanced translational diffusion is associated with spatially heterogeneous dynamics. Self-diffusion controls crystallization in o-terphenyl for most of the supercooled liquid regime, but at temperatures below Tg + 10 K, the reported crystallization rate increases suddenly while the self-diffusion coefficient does not. This work and previous work on trisnaphthylbenzene both find a self-diffusion-controlled crystal growth regime and an enhancement in self-diffusion near Tg, suggesting that these phenomena are general characteristics of fragile low molecular weight glass formers. We discuss the width of the relaxation time distributions of o-terphenyl and trisnaphthylbenzene as they relate to the observation of enhanced translational diffusion.     

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.     

Analysis of the subcritical carbon dioxide-water interface

We follow the evolution of the H(2)O/CO(2) interface at 300 K from the low pressure limit to near-critical pressures in molecular dynamics simulations using the SPC water and EPM2 carbon dioxide models. The intrinsic structure of the interface is elucidated by accumulating density profiles relative to the fluctuating capillary wave surface. Our main finding is that a carbon dioxide film of increasing density and thickness grows in two stages at the interface while the structure of the water surface barely changes. At low density, the entire film density profile grows linearly with the bulk CO(2) density. This regime continues up to a bulk CO(2) density of roughly 0.00095 ?(-3). At pressures above this point, we observe a distinct second peak in the CO(2) density, along with a tail of excess density that decays exponentially with distance from the interface. The decay length of the exponential tail diverges with increasing CO(2) pressure according to an inverse power law decay. Over the entire range of pressures, the CO(2) film had no detectable effect on the orientational order of the water surface. As expected, when the film of excess CO(2) at the interface grows, we find that the surface tension drops with increasing pressure. This is in qualitative accord with existing measurements, although the rate at which the surface tension falls with increasing pressure according to the SPC and EPM2 models is too small, indicating that the surface excess of CO(2) is underestimated by these models.     

Waterlike structural and excess entropy anomalies in liquid beryllium fluoride

The relationship between structural order metrics and the excess entropy is studied using the transferable rigid ion model (TRIM) of beryllium fluoride melt, which is known to display waterlike thermodynamic anomalies. The order map for liquid BeF2, plotted between translational and tetrahedral order metrics, shows a structurally anomalous regime, similar to that seen in water and silica melt, corresponding to a band of state points for which average tetrahedral (q(tet)) and translational (tau) order are strongly correlated. The tetrahedral order parameter distributions further substantiate the analogous structural properties of BeF2, SiO2, and H2O. A region of excess entropy anomaly can be defined within which the pair correlation contribution to the excess entropy (S2) shows an anomalous rise with isothermal compression. Within this region of anomalous entropy behavior, q(tet) and S2 display a strong negative correlation, indicating the connection between the thermodynamic and the structural anomalies. The existence of this region of excess entropy anomaly must play an important role in determining the existence of diffusional and mobility anomalies, given the excess entropy scaling of transport properties observed in many liquids.     

Ionic melts with waterlike anomalies: thermodynamic properties of liquid BeF2

Thermodynamic properties of liquid beryllium difluoride (BeF(2)) are studied using canonical ensemble molecular dynamics simulations of the transferable rigid ion model potential. The negative slope of the locus of points of maximum density in the temperature-pressure plane is mapped out. The excess entropy, computed within the pair correlation approximation, is found to show an anomalous increase with isothermal compression at low temperatures which will lead to diffusional as well as structural anomalies resembling those in water. The anomalous behavior of the entropy is largely connected with the behavior of the Be-F pair correlation function. The internal energy shows a T(35) temperature dependence. The pair correlation entropy shows a T(-25) temperature dependence only at high densities and temperatures. The correlation plots between internal energy and the pair correlation entropy for isothermal compression show the characteristic features expected of network-forming liquids with waterlike anomalies. The tagged particle potential energy distributions are shown to have a multimodal form at low temperatures and densities similar to those seen in other liquids with three-dimensional tetrahedral networks, such as water and silica.     

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.     

Influence of bond flexibility on the vapor-liquid phase equilibria of water

The authors performed Gibbs ensemble simulations on the vapor-liquid equilibrium of water to investigate the influence of incorporating intramolecular degrees of freedom in the simple point charge (SPC) water model. Results for vapor pressures, saturation densities, heats of vaporization, and the critical point for two different flexible models are compared with data for the corresponding rigid SPC and SPC/E models. They found that the introduction of internal vibrations, and also their parametrization, has an observable effect on the prediction of the vapor-liquid coexistence curve. The flexible SPC/Fw model, although optimized to describe bulk diffusion and dielectric constants at ambient conditions, gives the best prediction of saturation densities and the critical point of the examined models.     

Thermodynamic and dynamic anomalies in a one-dimensional lattice model of liquid water

We investigate the occurrence of waterlike thermodynamic and dynamic anomalous behavior in a one dimensional lattice gas model. The system thermodynamics is obtained using the transfer matrix technique and anomalies on density and thermodynamic response functions are found. When the hydrogen bond (molecules separated by holes) is more attractive than the van der Waals interaction (molecules in contact) a transition between two fluid structures is found at null temperature and high pressure. This transition is analogous to a 'critical point' and intimately connects the anomalies in density and in thermodynamic response functions. Monte Carlo simulations were performed in the neighborhood of this transition and used to calculate the self diffusion constant, which increases with density as in liquid water.