Charge-on-spring polarizable water models revisited: from water clusters to liquid water to ice

The properties of two improved versions of charge-on-spring (COS) polarizable water models (COS/G2 and COS/G3) that explicitly include nonadditive polarization effects are reported. In COS models, the polarization is represented via a self-consistently induced dipole moment consisting of a pair of separated charges. A previous polarizable water model (COS/B2), upon which the improved versions are based, was developed by Yu, Hansson, and van Gunsteren. To improve the COS/B2 model, which overestimated the dielectric permittivity, one additional virtual atomic site was used to reproduce the water monomer quadrupole moments besides the water monomer dipole moment in the gas phase. The molecular polarizability, residing on the virtual atomic site, and Lennard-Jones parameters for oxygen-oxygen interactions were varied to reproduce the experimental values for the heat of vaporization and the density of liquid water at room temperature and pressure. The improved models were used to study the properties of liquid water at various thermodynamic states as well as gaseous water clusters and ice. Overall, good agreement is obtained between simulated properties and those derived from experiments and ab initio calculations. The COS/G2 and COS/G3 models may serve as simple, classical, rigid, polarizable water models for the study of organic solutes and biopolymers. Due to its simplicity, COS type of polarization can straightforwardly be used to introduce explicit polarization into (bio)molecular force fields.     

The bend angle of water in ice Ih and liquid water: The significance of implementing the nonlinear monomer dipole moment surface in classical interaction potentials

The implementation of the physically accurate nonlinear dipole moment surface of the water monomer in the context of the Thole-type, polarizable, flexible interaction potential results in the only classical potential, which, starting from the gas phase value for the bend angle (104.52 degrees), reproduces its experimentally observed increase in the ice Ih lattice and in liquid water. This is in contrast to all other classical potentials to date, which predict a decrease of the monomer bend angle in ice Ih and in liquid water with respect to the gas phase monomer value. Simulations under periodic boundary conditions of several supercells consisting of up to 288 molecules of water used to sample the proton disorder in the ice Ih lattice yield an average value of vartheta(HOH)(I(h))=108.4 degrees +/-0.2 degrees for the minimized structures (T=0 K) and 108.1 degrees +/-2.8 degrees at T=100 K. Analogous simulations for liquid water predict an average value of vartheta(HOH)(liquid)=106.3 degrees +/-4.9 degrees at T=300 K. The increase of the monomer bend angle of water in condensed environments is attributed to the use of geometry-dependent charges that are used to describe the nonlinear character of the monomer's dipole moment surface. Our results suggest a new paradigm in the development of classical interaction potential models of water that can be used to describe condensed aqueous environments.     

Efficient sampling of ice structures by electrostatic switching

An electrostatic switching procedure is introduced that enables the systematic generation of high-quality ice configurations at various temperatures. Proton disordered ice Ih configurations were generated for the TIP4P water model at temperatures from 50 to 240 K, for the SPC/Fw water model from 100 to 240 K, and for the DC97 water model at 240 K. The resulting configurations were found to properly sample the canonical ensemble. The dielectric constant of ice Ih was determined from the net dipole fluctuation of the ice configurations. The calculated dielectric constant compares favorably with the study by Rick and Haymet [J. Chem. Phys. 2003, 118, 9291]. However, our method gives smaller error bars, especially at lower temperatures. At temperatures above 200 K, a type of hydrogen-bond defect is identified that cannot be categorized as a D or L type defect.     

The ice-vapor interface and the melting point of ice I(h) for the polarizable POL3 water model

We use molecular dynamics simulations to determine the melting point of ice I(h) for the polarizable POL3 water force field (Dang, L. X. J. Chem. Phys.1992, 97, 2659). Simulations are performed on a slab of ice I(h) with two free surfaces at several different temperatures. The analysis of the time evolution of the total energy in the course of the simulations at the set of temperatures yields the melting point of the POL3 model to be T(m) = 180 ± 10 K. Moreover, the results of the simulations show that the degree of hydrogen-bond disorder occurring in the bulk of POL3 ice is larger (at the corresponding degree of undercooling) than in ice modeled by nonpolarizable water models. These results demonstrate that the POL3 water force field is rather a poor model for studying ice and ice-liquid or ice-vapor interfaces. While a number of polarizable water models have been developed over the past years, little is known about their performance in simulations of supercooled water and ice. This study thus highlights the need for testing of the existing polarizable water models over a broad range of temperatures, pressures, and phases, and developing a new polarizable water force field, reliable over larger areas of the phase diagram.     

A general purpose model for the condensed phases of water: TIP4P/2005

A potential model intended to be a general purpose model for the condensed phases of water is presented. TIP4P/2005 is a rigid four site model which consists of three fixed point charges and one Lennard-Jones center. The parametrization has been based on a fit of the temperature of maximum density (indirectly estimated from the melting point of hexagonal ice), the stability of several ice polymorphs and other commonly used target quantities. The calculated properties include a variety of thermodynamic properties of the liquid and solid phases, the phase diagram involving condensed phases, properties at melting and vaporization, dielectric constant, pair distribution function, and self-diffusion coefficient. These properties cover a temperature range from 123 to 573 K and pressures up to 40,000 bar. The model gives an impressive performance for this variety of properties and thermodynamic conditions. For example, it gives excellent predictions for the densities at 1 bar with a maximum density at 278 K and an averaged difference with experiment of 7 x 10(-4) g/cm3.     

Melting points and thermal expansivities of proton-disordered hexagonal ice with several model potentials

A method of free energy calculation is proposed, which enables to cover a wide range of pressure and temperature. The free energies of proton-disordered hexagonal ice (ice Ih) and liquid water are calculated for the TIP4P [J. Chem. Phys. 79, 926 (1983)] model and the TIP5P [J. Chem. Phys. 112, 8910 (2000)] model. From the calculated free energy curves, we determine the melting point of the proton-disordered hexagonal ice at 0.1 MPa (atmospheric pressure), 50 MPa, 100 MPa, and 200 MPa. The melting temperatures at atmospheric pressure for the TIP4P ice and the TIP5P ice are found to be about T(m)=229 K and T(m)=268 K, respectively. The melting temperatures decrease as the pressure is increased, a feature consistent with the pressure dependence of the melting point for realistic proton-disordered hexagonal ice. We also calculate the thermal expansivity of the model ices. Negative thermal expansivity is observed at the low temperature region for the TIP4P ice, but not for the TIP5P ice at the ambient pressure.     

Comparison of selected polarizable and nonpolarizable water models in molecular dynamics simulations of ice I(h)

We present a molecular dynamics simulation study in which we determined the melting point of ice I(h) for the polarizable SWM4-NDP water model (Lamoureux et al., Chem. Phys. Lett., 2006, 418, 245-249) and compared the performance of several popular water force fields, both polarizable and nonpolarizable, in terms of melting temperature, stability and orientational structuring of ice. The simulations yield the melting temperature of SWM4-NDP ice as low as T(m) = 185 ± 10 K, despite the quadrupole moment of a SWM4-NDP water molecule being close to the experimental gas phase value. The results thus show that the dependence of T(m) on the molecular quadrupole, observed for the three- and four-site water models, is generally lost if polarization is explicitly included. The study also shows that adding polarizability to a planar three-charge water model increases orientational disorder in hexagonal ice. In addition, analysis of the tetrahedral order in bulk ice reveals a correlation between the pre-existing degree of orientational disorder in ice simulated using different polarizable and nonpolarizable models and the melting temperature of the models. Our findings thus suggest some new considerations regarding the role of polarization forces in a crystalline solid that may guide future development of reliable polarizable water models for ice.     

A non-polarizable model of water that yields the dielectric constant and the density anomalies of the liquid: TIP4Q

A four-site rigid water model is presented, whose parameters are fitted to reproduce the experimental static dielectric constant at 298 K, the maximum density of liquid water and the equation of state at low pressures. The model has a positive charge on each of the three atomic nuclei and a negative charge located at the bisector of the HOH bending angle. This charge distribution allows increasing the molecular dipole moment relative to four-site models with only three charges and improves the liquid dielectric constant at different temperatures. Several other properties of the liquid and of ice Ih resulting from numerical simulations with the model are in good agreement with experimental values over a wide range of temperatures and pressures. Moreover, the model yields the minimum density of supercooled water at 190 K and the minimum thermal compressibility at 310 K, close to the experimental values. A discussion is presented on the structural changes of liquid water in the supercooled region where the derivative of density with respect to temperature is a maximum.     

Quantum path integral simulation of isotope effects in the melting temperature of ice Ih

The isotope effect in the melting temperature of ice Ih has been studied by free energy calculations within the path integral formulation of statistical mechanics. Free energy differences between isotopes are related to the dependence of their kinetic energy on the isotope mass. The water simulations were performed by using the q-TIP4P/F model, a point charge empirical potential that includes molecular flexibility and anharmonicity in the OH stretch of the water molecule. The reported melting temperature at ambient pressure of this model (T=251?K) increases by 6.5±0.5 and 8.2±0.5?K upon isotopic substitution of hydrogen by deuterium and tritium, respectively. These temperature shifts are larger than the experimental ones (3.8 and 4.5 K, respectively). In the classical limit, the melting temperature is nearly the same as that for tritiated ice. This unexpected behavior is rationalized by the coupling between intermolecular interactions and molecular flexibility. This coupling makes the kinetic energy of the OH stretching modes larger in the liquid than in the solid phase. However, the opposite behavior is found for intramolecular modes, which display larger kinetic energy in ice than in liquid water.     

Classical interaction model for the water molecule

The authors propose a new classical model for the water molecule. The geometry of the molecule is built on the rigid TIP5P model and has the experimental gas phase dipole moment of water created by four equal point charges. The model preserves its rigidity but the size of the charges increases or decreases following the electric field created by the rest of the molecules. The polarization is expressed by an electric field dependent nonlinear polarization function. The increasing dipole of the molecule slightly increases the size of the water molecule expressed by the oxygen-centered sigma parameter of the Lennard-Jones interaction. After refining the adjustable parameters, the authors performed Monte Carlo simulations to check the ability of the new model in the ice, liquid, and gas phases. They determined the density and internal energy of several ice polymorphs, liquid water, and gaseous water and calculated the heat capacity, the isothermal compressibility, the isobar heat expansion coefficients, and the dielectric constant of ambient water. They also determined the pair-correlation functions of ambient water and calculated the energy of the water dimer. The accuracy of theirs results was satisfactory.     

Isotope effects in ice Ih: a path-integral simulation

Ice Ih has been studied by path-integral molecular dynamics simulations, using the effective q-TIP4P/F potential model for flexible water. This has allowed us to analyze finite-temperature quantum effects in this solid phase from 25 to 300 K at ambient pressure. Among these effects we find a negative thermal expansion of ice at low temperatures, which does not appear in classical molecular dynamics simulations. The compressibility derived from volume fluctuations gives results in line with experimental data. We have analyzed isotope effects in ice Ih by considering normal, heavy, and tritiated water. In particular, we studied the effect of changing the isotopic mass of hydrogen on the kinetic energy and atomic delocalization in the crystal as well as on structural properties such as interatomic distances and molar volume. For D(2)O ice Ih at 100 K we obtained a decrease in molar volume and intramolecular O-H distance of 0.6% and 0.4%, respectively, as compared to H(2)O ice.     

An application of flexible constraints in Monte Carlo simulations of the isobaric--isothermal ensemble of liquid water and ice Ih with the polarizable and flexible mobile charge densities in harmonic oscillators model

The method of flexible constraints was implemented in a Monte Carlo code to perform numerical simulations of liquid water and ice Ih in the constant number of molecules, volume, and temperature and constant pressure, instead of volume ensembles, using the polarizable and flexible mobile charge densities in harmonic oscillators (MCDHO) model. The structural and energetic results for the liquid at T=298 K and rho=997 kg m(-3) were in good agreement with those obtained from molecular dynamics. The density obtained at P=1 atm with flexible constraints, rho=1008 kg m(-3), was slightly lower than with the classical sampling of the intramolecular vibrations, rho=1010 kg m(-3). The comparison of the structures and energies found for water hexamers and for ice Ih with six standard empirical models to those obtained with MCDHO, show this latter to perform better in describing water far from ambient conditions: the MCDHO minimum lattice energy, density, and lattice constants were in good agreement with experiment. The average angle HOH of the water molecule in ice was predicted to be slightly larger than in the liquid, yet 1.2% smaller than the experimental value.     

Quantum effects in ice Ih

Quantum and classical simulations are carried out on ice Ih over a range of temperatures utilizing the TIP4P water model. The rigid-body centroid molecular dynamics method employed allows for the investigation of equilibrium and dynamical properties of the quantum system. The impact of quantization on the local structure, as measured by the radial and spatial distribution functions, as well as the energy is presented. The effects of quantization on the lattice vibrations, associated with the molecular translations and librations, are also reported. Comparison of quantum and classical simulation results indicates that shifts in the average potential energy are equivalent to rising the temperature about 80 K and are therefore non-negligible. The energy shifts due to quantization and the quantum mechanical uncertainties observed in ice are smaller than the values previously reported for liquid water. Additionally, we carry out a comparative study of melting in our classical and quantum simulations and show that there are significant differences between classical and quantum ice.     

Temperature dependence of crystal growth of hexagonal ice (I(h))

The transformations between water and ice have many implications across numerous fields of study. A better understanding of this process would benefit many areas of science and technology such as medicine, biology, and atmospheric and material sciences. In the present work the temperature dependence of the rate of growth (melting) of the basal face of hexagonal ice I(h) and the effect of system size are investigated in molecular dynamics simulations. Using an effective pair potential model of water, systems are studied over temperatures ranging from T(M) - 40 to T(M) + 16 K, where T(M) is the melting temperature of the model. It is found that the growth rates reach a maximum value of 0.7 ? ns(-1) (7 cm s(-1)) at about 12 K below the melting temperature. A noticeable effect of the system size on the melting temperature and ice growth rates is observed; it is shown that the size effect arises in smaller systems due to the artificial ordering under periodic conditions. The decrease in melting entropy in the smallest system by 0.4 J (mol K)(-1) relative to the largest system results in an up-shift in the melting temperature by about 2 K. An almost 60% increase in the maximum growth rate is observed for the smallest system.     

Molecular dynamics simulations of the thermal conductivity of methane hydrate

Nonequilibrium molecular dynamics simulations with the nonpolarizable SPC/E (Berendsen et al., J. Phys. Chem. 1987, 91, 6269) and the polarizable COS/G2 (Yu and van Gunsteren, J. Chem. Phys. 2004, 121, 9549) force fields have been employed to calculate the thermal conductivity and other associated properties of methane hydrate over a temperature range from 30 to 260 K. The calculated results are compared to experimental data over this same range. The values of the thermal conductivity calculated with the COS/G2 model are closer to the experimental values than are those calculated with the nonpolarizable SPC/E model. The calculations match the temperature trend in the experimental data at temperatures below 50 K; however, they exhibit a slight decrease in thermal conductivity at higher temperatures in comparison to an opposite trend in the experimental data. The calculated thermal conductivity values are found to be relatively insensitive to the occupancy of the cages except at low (T相似文献   

A potential model for the study of ices and amorphous water: TIP4P/Ice

The ability of several water models to predict the properties of ices is discussed. The emphasis is put on the results for the densities and the coexistence curves between the different ice forms. It is concluded that none of the most commonly used rigid models is satisfactory. A new model specifically designed to cope with solid-phase properties is proposed. The parameters have been obtained by fitting the equation of state and selected points of the melting lines and of the coexistence lines involving different ice forms. The phase diagram is then calculated for the new potential. The predicted melting temperature of hexagonal ice (Ih) at 1 bar is 272.2 K. This excellent value does not imply a deterioration of the rest of the properties. In fact, the predictions for both the densities and the coexistence curves are better than for TIP4P, which previously yielded the best estimations of the ice properties.     

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.     

Molecular dynamics simulations of liquid methanol and methanol-water mixtures with polarizable models

A polarizable model for simulation of liquid methanol, compatible with the COS/G2 water model, has been developed using the Charge-on-Spring (COS) technique. The model consists of three point charges, with one polarizable center on the oxygen atom. The Lennard-Jones parameters on the oxygen atom together with the molecular polarizability were varied to reproduce the experimental heat of vaporization and density of liquid methanol at ambient conditions. We examined the energies of various methanol dimers in the gas phase and compared them with values obtained from ab initio calculations. The model was then used to study the thermodynamic, dynamic, structural, and dielectric properties of liquid methanol as well as of a methanol-water mixture. A microscopic picture of the structure of pure liquid methanol and of the methanol-water mixture is provided. Good agreement was found between the results from our model simulations and available experimental and ab initio calculation data. In particular, the experimental dielectric permittivity of 32 could be reproduced, which had been shown to be difficult when using nonpolarizable models.     

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.