Molecular dynamics study of crystalline water ices

The behavior of structures of H₂O crystalline ices Ih, Ic, XI, VII, VIII, VI is studied in molecular dynamics experiment using the potential offered by Poltev and Malenkov. The behavior of the system consisting of one of the two identical interpenetrating, but without any common hydrogen bonds, water frameworks comprising the ice VI structure is also simulated. As a result of simulations, the ice VII structure has collapsed, whereas other systems proved to be stable. The reasons of instability of the ice VII and previously studied ice IV structures in molecular dynamics experiments are discussed. Based on the simulation results of the above-mentioned ices and previous simulation of ices II, III, IX, IV, and XII, the general regularities of dynamic properties of water molecules in crystalline water ices are formulated. Unreliability of results obtained by molecular dynamics in the investigation of self-organizing processes in aqueous systems is shown.     

Calculation of heat capacities of light and heavy water by path-integral molecular dynamics

As an application of atomistic simulation methods to heat capacities, path-integral molecular dynamics has been used to calculate the constant-volume heat capacities of light and heavy water in the gas, liquid, and solid phases. While the classical simulation based on conventional molecular dynamics has estimated the heat capacities too high, the quantum simulation based on path-integral molecular dynamics has given reasonable results based on the simple point-charge/flexible potential model. The calculated heat capacities (divided by the Boltzmann constant) in the quantum simulation are 3.1 in the vapor H2O at 300 K, 6.9 in the liquid H2O at 300 K, and 4.1 in the ice Ih H2O at 250 K, respectively, which are comparable to the experimental data of 3.04, 8.9, and 4.1, respectively. The quantum simulation also reproduces the isotope effect. The heat capacity in the liquid D2O has been calculated to be 10% higher than that of H2O, while it is 13% higher in the experiment. The results demonstrate that the path-integral simulation is a promising approach to quantitatively evaluate the heat capacities for molecular systems, taking account of quantum-mechanical vibrations as well as strongly anharmonic motions.     

A classical polarizable model for simulations of water and ice

We develop a classical rigid polarizable model of water for molecular simulations of water and ice. The model uses the Rowlinson five-site geometry: oxygen bearing the Lennard-Jones interaction and linearly polarizable point dipole, two positively charged hydrogens, and two massless negative charges placed symmetrically off oxygen so that the experimental dipole moment is reproduced. The target properties are the densities of water and ice Ih, diffusivity, enthalpies of fusion and vaporization, and the ice Ih melting point. The surface tension at lower temperatures is by 7% underestimated whereas the dielectric constant by 6% overestimated. Diffusivity and viscosity worsen at higher temperatures, although the Stokes radius is overestimated only by 2-7%. The ice Ih melting temperature is 260 K and the temperature of maximum density is 269 K. Rescaling the charges by a factor of 1.01 and Lennard-Jones energy by 1.0201 improves the melting point and energy-related quantities but shifts the agreement of kinetic properties to higher temperatures. For the model we propose abbreviation POL4D.     

Dynamics of TIP5P and TIP4P/ice potentials

The dynamics of a thin film of ice Ih deposited on MgO (001) is studied through molecular dynamics simulations performed with two new potential models of ice. This system is chosen because it is possible to compare the results of the simulations to incoherent neutron quasielastic scattering experiments performed few years ago and to previous molecular dynamics simulations using the TIP4P potential model. The present simulations are performed to determine the evolution of the translational and orientational order parameters of the ice film upon temperature increase in the 250-280 K range. They are also used to calculate the translational and orientational diffusion coefficients of the water molecules in the supported film as a function of the temperature. When using the TIP5P potential, the present results show a better agreement with experimental data than those calculated with the TIP4P potential, especially regarding the temperature above which significant changes are obtained in the dynamics of the water film. Similar conclusions are obtained when using the TIP4P/ice potential, although this latter potential clearly underestimates the translational diffusion coefficients.     

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 solvation dynamics of probe molecules in methanol-doped ice and in liquid ethanol

We have studied the solvation statics and dynamics of coumarin 343 and a strong photoacid (pK* approximately 0.7) 2-naphthol-6, 8-disulfonate (2N68DS) in methanol-doped ice (1% molar concentration of methanol) and in cold liquid ethanol in the temperature range of 160-270 K. Both probe molecules show a relatively fast solvation dynamics in ice, ranging from a few tens of picoseconds at about 240 K to nanoseconds at about 160 K. At about 160 K in doped ice, we observe a sharp decrease of the dynamic Stokes shift of both coumarin 343 and 2N68DS. Its value is approximately only 200 cm-1 at approximately 160 K compared to about 1100 cm-1 at T >/= 200 K (at times longer than t > 10 ps). We find a good correlation between the inefficient and slow excited-state proton-transfer rate at low-temperature ice, T < 180 K, and the dramatic decrease of the solvation energy, as measured by the dynamic band shift, at these low temperatures. We find that the average solvation rate in ice is similar to its value in liquid ethanol at all given temperatures in the range of 200-250 K. The surprisingly fast solvation rate in ice is explained by the relatively large freedom of the water hydrogen rotation in ice Ih.     

Quantum effects in liquid water and ice: model dependence

This paper explores the influence of choice of potential model on the quantum effects observed in liquid water and ice. This study utilizes standard rigid models and a more formal context for the rigid-body centroid molecular dynamics methodology used to perform the quantum simulations is provided. Quantum and classical molecular dynamics simulations are carried out for liquid water and ice Ih at 298 and 220 K, respectively, with the simple point charge/extended and TIP4P-Ew water models. The results obtained for equilibrium and dynamical properties are compared with those recently reported on TIP4P [L. Hernandez de la Pena and P. G. Kusalik, J. Chem. Phys. 121, 5992 (2004); L. Hernandez de la Pena et al., J. Chem. Phys 123, 144506 (2005)]. For the liquid, an energy shift of about 8% and an average molecular uncertainty of about 11 degrees were found independently of the water model. The self-diffusion coefficient consistently increases by more than 50% when going from the classical to the quantum system and quantum dynamics are found to reproduce the experimental isotopic shifts with the models examined. The ice results compare remarkably well with those previously reported for the TIP4P water model; they confirm that quantum effects are considerable and that the quantum mechanical uncertainty and the energy shifts due to quantization are smaller in ice than in liquid water. The relevance of these findings in the context of the construction of water models is briefly discussed.     

Quasi-harmonic approximation of thermodynamic properties of ice Ih, II, and III

Several thermodynamic properties of ice Ih, II, and III are studied by a quasi-harmonic approximation and compared to results of quantum path integral and classical simulations. This approximation allows to obtain thermodynamic information at a fraction of the computational cost of standard simulation methods, and at the same time permits studying quantum effects related to zero-point vibrations of the atoms. Specifically, we have studied the crystal volume, bulk modulus, kinetic energy, enthalpy, and heat capacity of the three ice phases as a function of temperature and pressure. The flexible q-TIP4P/F model of water was employed for this study, although the results concerning the capability of the quasi-harmonic approximation are expected to be valid independently of the employed water model. The quasi-harmonic approximation reproduces with reasonable accuracy the results of quantum and classical simulations showing an improved agreement at low temperatures (T< 100 K). This agreement does not deteriorate as a function of pressure as long as it is not too close to the limit of mechanical stability of the ice phases.     

Two-dimensional infrared spectroscopy of isotope-diluted ice Ih

We present experimental 2D IR spectra of isotope diluted ice Ih (i.e., the OH stretch mode of HOD in D(2)O and the OD stretch mode of HOD in H(2)O) at T = 80 K. The main spectral features are the extremely broad 1-2 excited state transition, much broader than the corresponding 0-1 groundstate transition, as well as the presence of quantum beats. We do not observe any inhomogeneous broadening that might be expected due to proton disorder in ice Ih. Complementary, we perform simulations in the framework of the Lippincott-Schroeder model, which qualitatively reproduce the experimental observations. We conclude that the origin of the observed line shape features is the coupling of the OH-vibrational coordinate with crystal phonons and explain the beatings as a coherent oscillation of the O···O hydrogen bond degree of freedom.     

Hydration structure of water confined between mica surfaces

We report further molecular dynamics simulations on the structure of bound hydration layers under extreme confinement between mica surfaces. We find that the liquid phase of water is maintained down to 2 monolayer (ML) thick, whereas the structure of the K(+) ion hydration shell is close to the bulk structure even under D = 0.92 nm confinement. Unexpectedly, the density of confined water remains approximately the bulk value or less, whereas the diffusion of water molecules decreases dramatically. Further increase in confinement leads to a transition to a bilayer ice, whose density is much less than that of ice Ih due to the formation of a specific hydrogen-bonding network.     

Investigation of the hydrogen bonding in ice Ih by first-principles density function methods

It is a well recognized difficult task to simulate the vibrational dynamics of ices using the density functional theory (DFT), and there has thus been rather limited success in modelling the inelastic neutron scattering (INS) spectra for even the simplest structure of ice, ice Ih, particularly in the translational region below 400 cm(-1). The reason is partly due to the complex nature of hydrogen bonding (H-bond) among water-water molecules which require considerable improvement of the quantum mechanical simulation methods, and partly owing to the randomness of protons in ice structures which often requires simulation of large super-lattices. In this report, we present the first series of successful simulation results for ice Ih using DFT methods. On the basis of the recent advancement in the DFT programs, we have achieved for the first time theoretical outcomes that not only reproduce the rotational frequencies between 500 to 1200 cm(-1) for ice Ih, but also the two optic peaks at ～240 and 320 cm(-1) in the translational region of the INS spectra [J. C. Li, J. Chem. Phys 105, 6733 (1996)]. Besides, we have also investigated the impact of pairwise configurations of H(2)O molecules on the H-bond and found that different proton arrangements of pairwise H(2)O in the ice Ih crystal lattice could not alter the nature of H-bond as significantly as suggested in an early paper [J. C. Li and D. K. Ross, Nature (London) 365, 327 (1993)], i.e., reproducing the two experimental optic peaks do not need to invoke the two H-bonds as proposed in the previous model which led to considerable debates. The results of this work suggest that the observed optic peaks may be attributed to the coupling between the two bands of H-O stretching modes in H(2)O. The current computational work is expected to shed new light on the nature of the H-bonds in water, and in addition to offer a new approach towards probing the interaction between water and biomaterials for which H-bond is essential.     

Transient absorption of vibrationally excited ice Ih

The ultrafast dynamics of HDO:D2O ice Ih at 180 K is studied by midinfrared ultrafast pump-probe spectroscopy. The vibrational relaxation of HDO:D2O ice is observed to proceed via an intermediate state, which has a blueshifted absorption spectrum. Polarization resolved measurements reveal that the intermediate state is part of the intramolecular relaxation pathway of the HDO molecule. In addition, slow dynamics on a time scale of the order of 10-100 ps is observed, related to thermally induced collective reorganizations of the ice lattice. The transient absorption line shape is analyzed within a Lippincott-Schroeder model for the OH-stretch potential. This analysis identifies the main mechanism behind the strong spectral broadening of the v(OH)=1-->2 transition.     

Antifreezes Act as Catalysts for Methane Hydrate Formation from Ice

Contrary to the thermodynamic inhibiting effect of methanol on methane hydrate formation from aqueous phases, hydrate forms quickly at high yield by exposing frozen water–methanol mixtures with methanol concentrations ranging from 0.6–10 wt % to methane gas at pressures from 125 bars at 253 K. Formation rates are some two orders of magnitude greater than those obtained for samples without methanol and conversion of ice is essentially complete. Ammonia has a similar catalytic effect when used in concentrations of 0.3–2.7 wt %. The structure I methane hydrate formed in this manner was characterized by powder X‐ray diffraction and Raman spectroscopy. Steps in the possible mechanism of action of methanol were studied with molecular dynamics simulations of the Ih (0001) basal plane exposed to methanol and methane gas. Simulations show that methanol from a surface aqueous layer slowly migrates into the ice lattice. Methane gas is preferentially adsorbed into the aqueous methanol surface layer. Possible consequences of the catalytic methane hydrate formation on hydrate plug formation in gas pipelines, on large scale energy‐efficient gas hydrate formation, and in planetary science are discussed.     

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.     

Sticking of CO to crystalline and amorphous ice surfaces

We present results of classical trajectory calculations on the sticking of hyperthermal CO to the basal plane (0001) face of crystalline ice Ih and to the surface of amorphous ice Ia. The calculations were performed for normal incidence at a surface temperature Ts = 90 K for ice Ia, and at Ts = 90 and 150 K for ice Ih. For both surfaces, the sticking probability can be fitted to a simple exponentially decaying function of the incidence energy, Ei: Ps = 1.0e(-Ei(kJ/mol)/90(kJ/mol)) at Ts = 90 K. The energy transfer from the impinging molecule to the crystalline and the amorphous surface is found to be quite efficient, in agreement with the results of molecular beam experiments on the scattering of the similar molecule, N2, from crystalline and amorphous ice. However, the energy transfer is less efficient for amorphous than for crystalline ice. Our calculations predict that the sticking probability decreases with Ts for CO scattering from crystalline ice, as the energy transfer from the impinging molecule to the warmer surfaces becomes less efficient. At high Ei (up to 193 kJ/mol), no surface penetration occurs in the case of crystalline ice. However, for CO colliding with the amorphous surface, a penetrating trajectory was observed to occur into a large water pore. The molecular dynamics calculations predict that the average potential energy of CO adsorbed to ice Ih is -10.1 +/- 0.2 and -8.4 +/- 0.2 kJ/mol for CO adsorbed to ice Ia. These values are in agreement with previous experimental and theoretical data. The distribution of the potential energy of CO adsorbed to ice Ia was found to be wider (with a standard deviation sigma of 2.4 kJ/mol) than that of CO interacting with ice Ih (sigma = 2.0 kJ/mol). In collisions with ice Ia, the CO molecules scatter at larger angles and over a wider distribution of angles than in collisions with ice Ih.     

Homogeneous ice nucleation from supercooled water

Homogeneous ice nucleation from supercooled water was studied in the temperature range of 220-240 K through combining the forward flux sampling method (Allen et al., J. Chem. Phys., 2006, 124, 024102) with molecular dynamics simulations (FFS/MD), based on a recently developed coarse-grained water model (mW) (Molinero et al., J. Phys. Chem. B, 2009, 113, 4008). The calculated ice nucleation rates display a strong temperature dependence, ranging from 2.148 ± 0.635 × 10(25) m(-3) s(-1) at 220 K to 1.672 ± 0.970 × 10(-7) m(-3) s(-1) at 240 K. These rates can be fitted according to the classical nucleation theory, yielding an estimate of the effective ice-water interface energy γ(ls) of 31.01 ± 0.21 mJ m(-2) for the mW water model. Compared to experiments, our calculation underestimates the homogeneous ice nucleation rate by a few orders of magnitude. Possible reasons for the discrepancy are discussed. The nucleating ice embryo contains both cubic ice Ic and hexagonal ice Ih, with the fraction of each structure being roughly 50% when the critical size is reached. In particular, a novel defect structure containing nearly five-fold twin boundaries is identified in the ice clusters formed during nucleation. The way such defect structure is formed is found to be different from mechanisms proposed for the formation of the same defect in metallic nanoparticles and thin film. The quasi five-fold twin boundary structure found here is expected to occur in the crystallization of a wide range of materials with the diamond cubic structure, including ice.     

Thermodynamic model for water and high-pressure ices up to 2.2 GPa and down to the metastable domain

We propose a thermodynamic model of the properties of liquid water and ices I, III, V, and VI that can be used in the ranges of 0-2200 MPa and 180-360 K. This model is the first to be applicable to all H(2)O phases in these wide ranges, which exceed the stability domain of all phases. Developing empirical or semiempirical expressions for the specific volumes of liquid water or ices has been necessary. The model has been tested on available experimental data sets. The specific volume of liquid water is reproduced with an accuracy better than 1%. The error on the specific volume of ices remains within 2%. The model has also been used to describe the melting curves of high-pressure ice polymorphs and compared with new Simon equations fitting available data. Our calculations suggest a slight revision of the triple point positions in the H(2)O phase diagram. We have ensured the reliability of our model up to 1.5 GPa, and we have shown that it can be used with good confidence up to 2.2 GPa. In order to show the validity of this model in the low-temperature domains, the melting curve of ice Ih in the water-ammonia system has been modeled. This curve is reproduced with good accuracy down to 180 K, at a 1 bar pressure. It shows that this model can be used in further studies for modeling equilibriums involving liquid or solid phases of H(2)O under pressure and for investigating the effect of inhibitors in complex water-rich systems.     

Simulations of proton order and disorder in ice Ih

Computer simulations of ice Ih with different proton orientations are presented. Simulations of proton disordered ice are carried out using a Monte Carlo method which samples over proton degree of freedom, allowing for the calculation of the dielectric constant and for the examination of the degree of proton disorder. Simulations are also presented for two proton ordered structures of ice Ih, the ferroelectric Cmc2(1) structure or ice XI and the antiferroelectric Pna2(1) structure. These simulations indicate that a transition to a proton ordered phase occurs at low temperatures (below 80 K). The symmetry of the ordered phase is found to be dependent on the water potential. The stability of the two proton ordered structures is due to a balance of short-ranged interactions which tend to stabilize the Pna2(1) structure and longer-range interactions which stabilize the Cmc2(1) structure.