Computer simulation of two new solid phases of water: Ice XIII and ice XIV

NpT Monte Carlo simulations have been performed for two recently discovered solid phases of water which have been denoted as ice XIII and ice XIV C. G. Salzmann et al. [Science311, 1758 (2006)]. Several potential models of water were considered, namely, the traditional SPC/E, TIP4P, and TIP5P and the more recent TIP5P-E, TIP4P-Ew, TIP4P/Ice, and TIP4P/2005 models. Significant differences in density and oxygen-oxygen radial distribution functions are found between the predictions of the SPC/E, TIP5P, and the models of the TIP4P family. The models TIP4P/Ice and TIP4P/2005 provide the best estimates of the density.     

Predicting the melting temperature of ice-Ih with only electronic structure information as input

The melting temperature of ice-Ih was calculated with only electronic structure information as input by creating a problem-specific force field. The force field, Water model by AFM for Ice and Liquid (WAIL), was developed with the adaptive force matching (AFM) method by fitting to post-Hartree-Fock quality forces obtained in quantum mechanics∕molecular mechanics calculations. WAIL predicts the ice-Ih melting temperature to be 270 K. The model also predicts the densities of ice and water, the temperature of maximum density of water, the heat of vaporizations, and the radial distribution functions for both ice and water in good agreement with experimental measurements. The non-dissociative WAIL model is very similar to a flexible version of the popular TIP4P potential and has comparable computational cost. By customizing to problem-specific configurations with the AFM approach, the resulting model is remarkably more accurate than any variants of TIP4P for simulating ice-Ih and water in the temperature range from 253 K and 293 K under ambient pressure.     

Simulating vapor-liquid nucleation of water: A combined histogram-reweighting and aggregation-volume-bias Monte Carlo investigation for fixed-charge and polarizable models

The method of histogram-reweighting was integrated with a recently developed approach using aggregation-volume-bias Monte Carlo and self-adaptive umbrella sampling to develop the AVUS-HR algorithm that allows for exceedingly efficient calculations of nucleation properties over a wide range of thermodynamic conditions. Simulations were carried out for water using both fixed-charge and polarizable force fields belonging to the TIP4P family (namely, TIP4P, TIP4P-FQ, TIP4P-pol2, and TIP4P-pol3) to investigate the nucleation of water over a temperature range from 200 to 300 K and the concentration of water clusters in the atmosphere at elevations up to 15 km. It was found that the nucleation free energy barriers and atmospheric concentrations are extremely sensitive to the force field, albeit all of the models investigated in this study support the following general conclusions: (i) ice nucleation is not present under the thermodynamic conditions and cluster-size range investigated here, i.e., the critical nuclei possess liquidlike structures, and (ii) the atmospheric concentrations of water clusters under homogeneous conditions are very low with the mole fraction of hexamers being about 10(-10), a number probably too low to influence the solar radiation balance. Compared to the experimental data, the TIP4P-pol3 model yields the most accurate nucleation results, consistent with its excellent performance for the second virial coefficient and the minimum cluster energies.     

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.     

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.     

Melting and superheating of sI methane hydrate: molecular dynamics study

Melting and decay of the superheated sI methane structure are studied using molecular dynamics simulation. The melting curve is calculated by the direct coexistence simulations in a wide range of pressures up to 5000 bar for the SPC/E, TIP4P/2005 and TIP4P/Ice water models and the united-atom model for methane. We locate the kinetic stability boundary of the superheated metastable sI structure that is found to be surprisingly high comparing with the predictions based on the classical nucleation theory.     

Physical properties of small water clusters in low and moderate electric fields

Likely candidates for the lowest minima of water clusters (H(2)O)(N) for N ≤ 20 interacting with a uniform electric field strength in the range E ≤ 0.6 V/? have been identified using basin-hopping global optimization. Two water-water model potentials were considered, namely TIP4P and the polarizable Dang-Chang potential. The two models produce some consistent results but also exhibit significant differences. The cluster internal energy and dipole moment indicate two varieties of topological transition in the structure of the global minimum as the field strength is increased. The first takes place at low field strengths (0.1 V/? 10) usually forming helical structures.     

Kinetic aspects of the thermostatted growth of ice from supercooled water in simulations

In experiments, the growth rate of ice from supercooled water is seen to increase with the degree of supercooling, that is, the lower the temperature, the faster the crystallization takes place. In molecular dynamics simulations of the freezing process, however, the temperature is usually kept constant by means of a thermostat that artificially removes the heat released during the crystallization by scaling the velocities of the particles. This direct removal of energy from the system replaces a more realistic heat-conduction mechanism and is believed to be responsible for the curious observation that the thermostatted ice growth proceeds fastest near the melting point and more slowly at lower temperatures, which is exactly opposite to the experimental findings [M. A. Carignano, P. B. Shepson, and I. Szleifer, Mol. Phys. 103, 2957 (2005)]. This trend is explained by the diffusion and the reorientation of molecules in the liquid becoming the rate-determining steps for the crystal growth, both of which are slower at low temperatures. Yet, for a different set of simulations, a kinetic behavior analogous to the experimental finding has been reported [H. Nada and Y. Furukawa, J. Crystal Growth 283, 242 (2005)]. To clarify this apparent contradiction, we perform relatively long simulations of the TIP4P/Ice model in an extended range of temperatures. The temperature dependence of the thermostatted ice growth is seen to be more complex than was previously reported: The crystallization process is very slow close to the melting point at 270 K, where the thermodynamic driving force for the phase transition is weak. On lowering the temperature, the growth rate initially increases, but displays a maximum near 260 K. At even lower temperatures, the freezing process slows down again due to the reduced diffusivity in the liquid. The velocity of the thermostatted melting process, in contrast, shows a monotonic increase upon raising the temperature beyond the normal melting point. In this case, the effects of the increasing thermodynamic driving force and the faster diffusion at higher temperatures reinforce each other. In the context of this study, we also report data for the diffusion coefficient as a function of temperature for the water models TIP4P/Ice and TIP4P/2005.     

Hydration thermodynamic properties of amino acid analogues: a systematic comparison of biomolecular force fields and water models

We present an extensive study on hydration thermodynamic properties of analogues of 13 amino acid side chains at 298 K and 1 atm. The hydration free energies DeltaG, entropies DeltaS, enthalpies DeltaH, and heat capacities Deltac(P)() were determined for 10 combinations of force fields and water models. The statistical sampling was extended such that precisions of 0.3, 0.8, 0.8 kJ/mol and 25 J/(mol K) were reached for DeltaG, TDeltaS, DeltaH, and Deltac(P)(), respectively. The three force fields used in this study are AMBER99, GROMOS 53A6, and OPLS-AA; the five water models are SPC, SPC/E, TIP3P, TIP4P, and TIP4P-Ew. We found that the choice of water model strongly influences the accuracy of the calculated hydration entropies, enthalpies, and heat capacities, while differences in accuracy between the force fields are small. On the basis of an analysis of the hydrophobic analogues of the amino acid side chains, we discuss what properties of the water models are responsible for the observed discrepancies between computed and experimental values. The SPC/E water model performs best with all three biomolecular force fields.     

Effects of external electric fields on double proton transfer kinetics in the formic acid dimer

Molecules can be exposed to strong local electric fields of the order of 10(8)-10(10) V m(-1) in the biological milieu. The effects of such fields on the rate constant (k) of a model reaction, the double-proton transfer reaction in the formic acid dimer (FAD), are investigated. The barrier heights and shapes are calculated in the absence and presence of several static homogenous external fields ranging from 5.14 × 10(8) to 5.14 × 10(9) V m(-1) using density functional theory (DFT/B3LYP) and second order M?ller-Plesset perturbation theory (MP2) in conjunction with the 6-311++G(d,p) Pople basis set. Conventional transition state theory (CTST) followed by Wigner tunneling correction is then applied to estimate the rate constants at 25 °C. It is found that electric fields parallel to the long axis of the dimer (the line joining the two carbon atoms) lower the uncorrected barrier height, and hence increase the raw k. These fields also flatten the potential energy surface near the transition state region and, hence, decrease the multiplicative tunneling correction factor. The net result of these two opposing effects is that fields increase k(corrected) by a factor of ca. 3-4 (DFT-MP2, respectively) compared to the field-free k. Field strengths of ～3 × 10(9) V m(-1) are found to be sufficient to double the tunneling-corrected double proton transfer rate constant at 25 °C. Field strengths of similar orders of magnitudes are encountered in the scanning tunneling microscope (STM), in the microenvironment of a DNA base-pair, in an enzyme active site, and in intense laser radiation fields. It is shown that the net (tunneling corrected) effect of the field on k can be closely fitted to an exponential relationship of the form k = aexp(bE), where a and b are constants and E the electric field strength.     

Polarizable model of water with field-dependent polarization

The polarizable charge-on-spring model of water with three Gaussian charges developed by the present authors [A. Baranyai and P. T. Kiss, J. Chem. Phys. 133, 144109 (2010)] was studied. We introduced an analytic function for the polarizability in terms of the local electric field. Following theoretical suggestions, the polarizability decreases from its experimental gas-phase value, in our approach, toward a high-field threshold. Using this modified polarizability, we reparameterized the model by calculating its dielectric constant and obtained good estimates of density and internal energy for ambient water, hexagonal ice, and water cluster properties. Mimicked by the new model, we studied liquid water under the impact of homogeneous static electric field in the rage of 0-2.5 V/?. Both the density and the average dipole moment increase with the strength of the electric field. However, the internal energy shows a minimum at ~0.35 V/?. At this field strength, the model starts ordering into a crystal structure. At higher fields the liquid forms a crystalline structure which is a special version of cubic ice.     

Computer simulation study of metastable ice VII and amorphous phases obtained by its melting

Molecular dynamics simulations of metastable ice VII and cubic ice Ic are carried out in order to examine (1) the ability of commonly used water interaction potentials to reproduce the properties of ices, and (2) the possibility of generating low-density amorphous (LDA) structures by heating ice VII, which is known to transform to LDA at approximately 135 K at normal pressure [S. Klotz, J. M. Besson, G. Hamel, R. J. Nelmes, J. S. Loveday, and W. G. Marshall, Nature (London) 398, 681 (1999)]. We test four simple empirical interaction potentials of water: TIP4P [W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983)], SPC/E [H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. B 91, 6269 (1987)], TIP5P [M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112, 8910 (2000)], and ST2 [F. H. Stillinger and A. Rahman, J. Chem. Phys. 60, 1545 (1974)]. We have found that TIP5P ice VII melts at 210 K, TIP4P at 90 K, and SPC/E at 70 K. Only TIP5P water after transition has a structure similar to that of LDA. TIP4P and SPC/E have almost identical structures, dissimilar to any known water or amorphous phases, but upon heating both slowly evolve towards LDA-like structure. ST2 ice VII is remarkably stable up to 430 K. TIP4P and SPC/E predict correctly the cubic ice collapse into a high-density amorphous ice (HDA) at approximately 1 GPa whereas TIP5P remains stable up to approximately 5 GPa. The densities of the simulated ice phases differ significantly, depending on the potential used, and are generally higher than experimental values. The importance of proper treatment of long-range electrostatic interactions is also discussed.     

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.     

The melting point of ice Ih for common water models calculated from direct coexistence of the solid-liquid interface

In this work we present an implementation for the calculation of the melting point of ice I(h) from direct coexistence of the solid-liquid interface. We use molecular dynamics simulations of boxes containing liquid water and ice in contact. The implementation is based on the analysis of the evolution of the total energy along NpT simulations at different temperatures. We report the calculation of the melting point of ice I(h) at 1 bar for seven water models: SPC/E, TIP4P, TIP4P-Ew, TIP4P/ice, TIP4P/2005, TIP5P, and TIP5P-E. The results for the melting temperature from the direct coexistence simulations of this work are in agreement (within the statistical uncertainty) with those obtained previously by us from free energy calculations. By taking into account the results of this work and those of our free energy calculations, recommended values of the melting point of ice I(h) at 1 bar for the above mentioned water models are provided.     

The series of molecular conductors and superconductors ET4[AFe(C2O4)3]·PhX (ET = bis(ethylenedithio)tetrathiafulvalene; (C2O4)2- = oxalate; A+ = H3O+, K+; X = F, Cl, Br, and I): influence of the halobenzene guest molecules on the crystal structure and superconducting properties

An extensive series of radical salts formed by the organic donor bis(ethylenedithio)tetrathiafulvalene (ET), the paramagnetic tris(oxalato)ferrate(III) anion [Fe(C(2)O(4))(3)](3-), and halobenzene guest molecules has been synthesized and characterized. The change of the halogen atom in this series has allowed the study of the effect of the size and charge polarization on the crystal structures and physical properties while keeping the geometry of the guest molecule. The general formula of the salts is ET(4)[A(I)Fe(C(2)O(4))(3)]·G with A/G = H(3)O(+)/PhF (1); H(3)O(+)/PhCl (2); H(3)O(+)/PhBr (3), and K(+)/PhI (4), (crystal data at room temperature: (1) monoclinic, space group C2/c with a = 10.3123(2) ?, b = 20.0205(3) ?, c = 35.2732(4) ?, β = 92.511(2)°, V = 7275.4(2) ?(3), Z = 4; (2) monoclinic, space group C2/c with a = 10.2899(4) ?, b = 20.026(10) ?, c = 35.411(10) ?, β = 92.974°, V = 7287(4) ?(3), Z = 4; (3) monoclinic, space group C2/c with a = 10.2875(3) ?, b = 20.0546(15) ?, c = 35.513(2) ?, β = 93.238(5)°, V = 7315.0(7) ?(3), Z = 4; (4) monoclinic, space group C2/c with a = 10.2260(2) ?, b = 19.9234(2) ?, c = 35.9064(6) ?, β = 93.3664(6)°, V = 7302.83(18) ?(3), Z = 4). The crystal structures at 120 K evidence that compounds 1-3 undergo a structural transition to a lower symmetry phase when the temperature is lowered (crystal data at 120 K: (1) triclinic, space group P1 with a = 10.2595(3) ?, b = 11.1403(3) ?, c = 34.9516(9) ?, α = 89.149(2)°, β = 86.762(2)°, γ = 62.578(3)°, V = 3539.96(19) ?(3), Z = 2; (2) triclinic, space group P1 with a = 10.25276(14) ?, b = 11.15081(13) ?, c = 35.1363(5) ?, α = 89.0829(10)°, β = 86.5203(11)°, γ = 62.6678(13)°, V = 3561.65(8) ?(3), Z = 2; (3) triclinic, space group P1 with a = 10.25554(17) ?, b = 11.16966(18) ?, c = 35.1997(5) ?, α = 62.7251(16)°, β = 86.3083(12)°, γ = 62.7251(16)°, V = 3575.99(10) ?(3), Z = 2; (4) monoclinic, space group C2/c with a = 10.1637(3) ?, b = 19.7251(6) ?, c = 35.6405(11) ?, β = 93.895(3)°, V = 7128.7(4) ?(3), Z = 4). A detailed crystallographic study shows a change in the symmetry of the crystal for compound 3 at about 200 K. This structural transition arises from the partial ordering of some ethylene groups in the ET molecules and involves a slight movement of the halobenzene guest molecules (which occupy hexagonal cavities in the anionic layers) toward one of the adjacent organic layers, giving rise to two nonequivalent organic layers at 120 K (compared to only one at room temperature). The structural transition at about 200 K is also observed in the electrical properties of 1-3 and in the magnetic properties of 1. The direct current (dc) conductivity shows metallic behavior in salts 1-3 with superconducting transitions at about 4.0 and 1.0 K in salts 3 and 1, respectively. Salt 4 shows a semiconductor behavior in the temperature range 300-50 K with an activation energy of 64 meV. The magnetic measurements confirm the presence of high spin S = 5/2 [Fe(C(2)O(4))(3)](3-) isolated monomers together with a Pauli paramagnetism, typical of metals, in compounds 1-3. The magnetic properties can be very well reproduced in the whole temperature range with a simple model of isolated S = 5/2 ions with a zero field splitting plus a temperature independent paramagnetism (Nα) with the following parameters: g = 1.965, |D| = 0.31 cm(-1), and Nα = 1.5 × 10(-3) emu mol(-1) for 1, g = 2.024, |D| = 0.65 cm(-1), and Nα = 1.4 × 10(-3) emu mol(-1) for 2, and g = 2.001, |D| = 0.52 cm(-1), and Nα = 1.5 × 10(-3) emu mol(-1) for 3.     

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.