Three-phase equilibria of hydrates from computer simulation. I. Finite-size effects in the methane hydrate.

Clathrate hydrates are vital in energy research and environmental applications. Understanding their stability is crucial for harnessing their potential. In this work, we employ direct coexistence simulations to study finite-size effects in the determination of the three-phase equilibrium temperature (T3) for methane hydrates. Two popular water models, TIP4P/Ice and TIP4P/2005, are employed, exploring various system sizes by varying the number of molecules in the hydrate, liquid, and gas phases. The results reveal that finite-size effects play a crucial role in determining T3. The study includes nine configurations with varying system sizes, demonstrating that smaller systems, particularly those leading to stoichiometric conditions and bubble formation, may yield inaccurate T3 values. The emergence of methane bubbles within the liquid phase, observed in smaller configurations, significantly influences the behavior of the system and can lead to erroneous temperature estimations. Our findings reveal finite-size effects on the calculation of T3 by direct coexistence simulations and clarify the system size convergence for both models, shedding light on discrepancies found in the literature. The results contribute to a deeper understanding of the phase equilibrium of gas hydrates and offer valuable information for future research in this field.

Three-phase equilibria of hydrates from computer simulation. II. Finite-size effects in the carbon dioxide hydrate.

In this work, the effects of finite size on the determination of the three-phase coexistence temperature (T3) of the carbon dioxide (CO2) hydrate have been studied by molecular dynamic simulations and using the direct coexistence technique. According to this technique, the three phases involved (hydrate-aqueous solution-liquid CO2) are placed together in the same simulation box. By varying the number of molecules of each phase, it is possible to analyze the effect of simulation size and stoichiometry on the T3 determination. In this work, we have determined the T3 value at 8 different pressures (from 100 to 6000 bar) and using 6 different simulation boxes with different numbers of molecules and sizes. In two of these configurations, the ratio of the number of water and CO2 molecules in the aqueous solution and the liquid CO2 phase is the same as in the hydrate (stoichiometric configuration). In both stoichiometric configurations, the formation of a liquid drop of CO2 in the aqueous phase is observed. This drop, which has a cylindrical geometry, increases the amount of CO2 available in the aqueous solution and can in some cases lead to the crystallization of the hydrate at temperatures above T3, overestimating the T3 value obtained from direct coexistence simulations. The simulation results obtained for the CO2 hydrate confirm the sensitivity of T3 depending on the size and composition of the system, explaining the discrepancies observed in the original work by Míguez et al. [J. Chem Phys. 142, 124505 (2015)]. Non-stoichiometric configurations with larger unit cells show a convergence of T3 values, suggesting that finite-size effects for these system sizes, regardless of drop formation, can be safely neglected. The results obtained in this work highlight that the choice of a correct initial configuration is essential to accurately estimate the three-phase coexistence temperature of hydrates by direct coexistence simulations.

Three-phase equilibria of hydrates from computer simulation. III. Effect of dispersive interactions in the methane and carbon dioxide hydrates.

In this work, the effect of the range of dispersive interactions in determining the three-phase coexistence line of the CO2 and CH4 hydrates has been studied. In particular, the temperature (T3) at which solid hydrate, water, and liquid CO2/gas CH4 coexist has been determined through molecular dynamics simulations using different cutoff values (from 0.9 to 1.6 nm) for dispersive interactions. The T3 of both hydrates has been determined using the direct coexistence simulation technique. Following this method, the three phases in equilibrium are put together in the same simulation box, the pressure is fixed, and simulations are performed at different temperatures T. If the hydrate melts, then T > T3. Conversely, if the hydrate grows, then T < T3. The effect of the cutoff distance on the dissociation temperature has been analyzed at three different pressures for CO2 hydrate: 100, 400, and 1000 bar. Then, we have changed the guest and studied the effect of the cutoff distance on the dissociation temperature of the CH4 hydrate at 400 bar. Moreover, the effect of long-range corrections for dispersive interactions has been analyzed by running simulations with homo- and inhomogeneous corrections and a cutoff value of 0.9 nm. The results obtained in this work highlight that the cutoff distance for the dispersive interactions affects the stability conditions of these hydrates. This effect is enhanced when the pressure is decreased, displacing the T3 about 2-4 K depending on the system and the pressure.

Scaled charges for ions: An improvement but not the final word for modeling electrolytes in water.

In this work, we discuss the use of scaled charges when developing force fields for NaCl in water. We shall develop force fields for Na+ and Cl- using the following values for the scaled charge (in electron units): ±0.75, ±0.80, ±0.85, and ±0.92 along with the TIP4P/2005 model of water (for which previous force fields were proposed for q = ±0.85 and q = ±1). The properties considered in this work are densities, structural properties, transport properties, surface tension, freezing point depression, and maximum in density. All the developed models were able to describe quite well the experimental values of the densities. Structural properties were well described by models with charges equal to or larger than ±0.85, surface tension by the charge ±0.92, maximum in density by the charge ±0.85, and transport properties by the charge ±0.75. The use of a scaled charge of ±0.75 is able to reproduce with high accuracy the viscosities and diffusion coefficients of NaCl solutions for the first time. We have also considered the case of KCl in water, and the results obtained were fully consistent with those of NaCl. There is no value of the scaled charge able to reproduce all the properties considered in this work. Although certainly scaled charges are not the final word in the development of force fields for electrolytes in water, its use may have some practical advantages. Certain values of the scaled charge could be the best option when the interest is to describe certain experimental properties.

The Madrid-2019 force field for electrolytes in water using TIP4P/2005 and scaled charges: Extension to the ions F-, Br-, I-, Rb+, and Cs.

In this work, an extension of the Madrid-2019 force field is presented. We have added the cations Rb+ and Cs+ and the anions F-, Br-, and I-. These ions were the remaining alkaline and halogen ions, not previously considered in the Madrid-2019 force field. The force field, denoted as Madrid-2019-Extended, does not include polarizability and uses the TIP4P/2005 model of water and scaled charges for the ions. A charge of ±0.85e is assigned to monovalent ions. The force field developed provides an accurate description of aqueous solution densities over a wide range of concentrations up to the solubility limit of each salt studied. Good predictions of viscosity and diffusion coefficients are obtained for concentrations below 2 m. Structural properties obtained with this force field are also in reasonable agreement with the experiment. The number of contact ion pairs has been controlled to be low so as to avoid precipitation of the system at concentrations close to the experimental solubility limit. A comprehensive comparison of the performance for aqueous solutions of alkaline halides of force fields of electrolytes using scaled and integer charges is now possible. This comparison will help in the future to learn about the benefits and limitations of the use of scaled charges to describe electrolyte solutions.

Phase diagram of the NaCl-water system from computer simulations.

NaCl aqueous solutions are ubiquitous. They can crystallize into ice, NaCl, or NaCl · 2H2O depending on the temperature-concentration conditions. These crystallization transitions have important implications in geology, cryopreservation, or atmospheric science. Computer simulations can help understand the crystallization of these solids, which requires a detailed knowledge of the equilibrium phase diagram. We use molecular simulations in which we put at contact the solution with the solid of interest to determine points of the solid-solution coexistence lines. We follow two different approaches, one in which we narrow down the melting temperature for a given concentration and the other in which we equilibrate the concentration for a given temperature, obtaining consistent results. The phase diagram thus calculated for the selected model (TIP4P/2005 for water molecules and Joung-Cheatham for the ions) correctly predicts coexistence between the solution and ice. We were only able to determine NaCl · 2H2O-solution coexistence points at higher temperatures and concentrations than in the experiment, so we could not establish a direct comparison in this case. On the other hand, the model underestimates the concentration of the solution in equilibrium with the NaCl solid. Our results, alongside other literature evidence, seem to indicate that ion-ion interactions are too strong in the model. Our work is a good starting point for the improvement of the potential model and for the study of the nucleation kinetics of the solid phases involved in the phase diagram.

Spontaneous NaCl-doped ices Ih, Ic, III, V and VI. Understanding the mechanism of ion inclusion and its dependence on the crystalline structure of ice.

Direct coexistence simulations on a microsecond time scale have been performed for different types of ice (Ih, Ic, III, V, and VI) in contact with a NaCl aqueous solution at different pressures. In line with the previous results obtained for ice Ih [Conde et al., Phys. Chem. Chem. Phys., 2017, 19, 9566-9574], our results reveal the spontaneous growth of a new ice doped phase and the formation of a brine rejection phase in all ices studied. However, both the preferential incorporation of ions into the ice lattice and the inclusion mechanisms depend on the crystalline structure of each ice. This work shows the inclusion of Cl- and Na+ ions in ice from salt using molecular dynamics simulation, in agreement with the experimental evidence found in the literature. The model used for water is TIP4P/2005. For NaCl we employ a set of potential parameters that uses unit charges for the ions.

Freezing Temperatures, Ice Nanotubes Structures, and Proton Ordering of TIP4P/ICE Water inside Single Wall Carbon Nanotubes.

A very recent experimental paper importantly and unexpectedly showed that water in carbon nanotubes is already in the solid ordered phase at the temperature where bulk water boils. The water models used so far in literature for molecular dynamics simulations in carbon nanotubes show freezing temperatures lower than the experiments. We present here results from molecular dynamics simulations of water inside single walled carbon nanotubes using an extremely realistic model for both liquid and icy water, the TIP4P/ICE. The water behavior inside nanotubes of different diameters has been studied upon cooling along the isobars at ambient pressure starting from temperatures where water is in a liquid state. We studied the liquid/solid transition, and we observed freezing temperatures higher than in bulk water and that depend on the diameter of the nanotube. The maximum freezing temperature found is 390 K, which is in remarkable agreement with the recent experimental measurements. We have also analyzed the ice structure called "ice nanotube" that water forms inside the single walled carbon nanotubes when it freezes. The ice forms observed are in agreement with previous results obtained with different water models. A novel finding, a partial proton ordering, is evidenced in our ice nanotubes at finite temperature.

Spontaneous NaCl-doped ice at seawater conditions: focus on the mechanisms of ion inclusion.

Molecular dynamics simulations on the microsecond time scale have been performed on an aqueous solution of TIP4P/2005 water and NaCl by using the direct coexistence technique to study the ice growth and the ice/liquid interface water. At ambient pressure, for temperatures above the eutectic point of the salt and at seawater concentrations the brine rejection phenomenon and the spontaneous growth of an ice slab doped by the salt are obtained, as found in natural terrestrial and planetary environments. Experiments indicate that Cl- goes via substitution to ice sites. In line with this evidence we find a new result: the Cl- ion included in the lattice always substitutes not one but two water molecules, leaving the surrounding ice structure not distorted. The Na+ ion shows a lower probability of being included in the ice and it occupies an interstitial site, causing a local distortion of the lattice. No signs of significant ion diffusion are observed in the lattice.

Structural properties and fragile to strong transition in confined water.

We derive by computer simulation the radial distribution functions of water confined in a silica pore modeled to reproduce MCM-41. We perform the calculations in a range of temperatures from ambient to deep supercooling for the subset of water molecules that reside in the inner shell (free water) by applying the excluded volume corrections. By comparing with bulk water we find that the first shell of the oxygen-oxygen and hydrogen-hydrogen radial distribution functions is less sharp and the first minimum fills in while the oxygen-hydrogen structure does not significantly change, indicating that the free water keeps the hydrogen bond short range order. The two body excess entropy of supercooled water is calculated from the radial distribution functions. We connect the behavior of this function to the relaxation time of the same system already studied in previous simulations. We show that the two body entropy changes its behavior in coincidence with the crossover of the relaxation time from the mode coupling fragile to the strong Arrhenius regime. As for bulk water also in confinement, the two body entropy has a strict connection with the dynamical relaxation.

High precision determination of the melting points of water TIP4P/2005 and water TIP4P/Ice models by the direct coexistence technique.

An exhaustive study by molecular dynamics has been performed to analyze the factors that enhance the precision of the technique of direct coexistence for a system of ice and liquid water. The factors analyzed are the stochastic nature of the method, the finite size effects, and the influence of the initial ice configuration used. The results obtained show that the precision of estimates obtained through the technique of direct coexistence is markedly affected by the effects of finite size, requiring systems with a large number of molecules to reduce the error bar of the melting point. This increase in size causes an increase in the simulation time, but the estimate of the melting point with a great accuracy is important, for example, in studies on the ice surface. We also verified that the choice of the initial ice Ih configuration with different proton arrangements does not significantly affect the estimate of the melting point. Importantly this study leads us to estimate the melting point at ambient pressure of two of the most popular models of water, TIP4P/2005 and TIP4P/Ice, with the greatest precision to date.

Revisiting the thermodynamic modelling of type I gas-hydroquinone clathrates.

Under specific pressure and temperature conditions, certain gaseous species can be engaged in a host lattice of hydroquinone molecules, forming a supramolecular entity called a gas hydroquinone clathrate. This study is devoted to the thermodynamic modelling of type I hydroquinone clathrates. The gases considered in this work are argon, krypton, xenon, methane, nitrogen, oxygen and hydrogen sulphide. The basic van der Waals and Platteeuw model, which is, for example, not able to predict well the phase equilibrium properties of such clathrates at high temperature, is modified and extended by considering first the solubility of the guest in solid HQ and then the mutual interactions between the gaseous molecules inside the clathrate structure (i.e. guest-guest interactions). Other improvements of the basic theory, such as the choice of the reference state, are proposed, and a unique set of thermodynamic parameters valid for all the studied guests are finally calculated. Very good agreement is obtained between the model predictions and the experimental data available in the literature. Our results clearly demonstrate that the highest level of theory is necessary to describe well both the triphasic equilibrium line (where the HQ clathrate, the native hydroquinone HQα and the gas coexist), the occupancy of the guest in the clathrate, and the intercalation enthalpy.

Molecular dynamics simulation of CO2 hydrates: Prediction of three phase coexistence line.

The three phase equilibrium line (hydrate-liquid water-liquid carbon dioxide) has been estimated for the water + carbon dioxide binary mixture using molecular dynamics simulation and the direct coexistence technique. Both molecules have been represented using rigid nonpolarizable models. TIP4P/2005 and TIP4P/Ice were used for the case of water, while carbon dioxide was considered as a three center linear molecule with the parameterizations of MSM, EPM2, TraPPE, and ZD. The influence of the initial guest occupancy fraction on the hydrate stability has been analyzed first in order to determine the optimal starting configuration for the simulations, paying attention to the influence of the two different cells existing in the sI hydrate structure. The three phase coexistence temperature was then determined for a pressure range from 2 to 500 MPa. The qualitative shape of the equilibrium curve estimated is correct, including the high pressure temperature maximum that determines the hydrate re-entrant behaviour. However, in order to obtain quantitative agreement with experimental results, a positive deviation from the classical Lorentz-Berthelot combining rules must be considered.

Determining the phase diagram of water from direct coexistence simulations: the phase diagram of the TIP4P/2005 model revisited.

Direct coexistence simulations between the fluid and solid phases are performed for several ices. For ices Ih and VII it has already been shown that the methodology is successful and the melting point is in agreement with that obtained from free energy calculations. In this work the methodology is applied to ices II, III, V, and VI. The lengths of the direct coexistence runs for the high pressure polymorphs are not too long and last less than 20 ns for all ices except for ice II where longer runs (of about 150 ns) are needed. For ices II, V, and VI the results obtained are completely consistent with those obtained from free energy calculations. However, for ice III it is found that the melting point from direct coexistence simulations is higher than that obtained from free energy calculations, the difference being greater than the statistical error. Since ice III presents partial proton orientational disorder, the departure is attributed to differences in the partial proton order in the water model with respect to that found in the experiment. The phase diagram of the TIP4P/2005 model is recalculated using the melting points obtained from direct coexistence simulations. The new phase diagram is similar to the previous one except for the coexistence lines where ice III is involved. The range of stability of ice III on the p-T plot of the phase diagram increases significantly. It is seen that the model qualitatively describes the phase diagram of water. In this work it is shown that the complete phase diagram of water including ices Ih, II, III, V, VI, VII, and the fluid phase can be obtained from direct coexistence simulations without the need of free energy calculations.

Determining the three-phase coexistence line in methane hydrates using computer simulations.

Molecular dynamics simulations have been performed to estimate the three-phase (solid hydrate-liquid water-gaseous methane) coexistence line for the water-methane binary mixture. The temperature at which the three phases are in equilibrium was determined for three different pressures, namely, 40, 100, and 400 bar by using direct coexistence simulations. In the simulations water was described by using either TIP4P, TIP4P/2005, or TIP4P/Ice models and methane was described as simple Lennard-Jones interaction site. Lorentz-Berthelot combining rules were used to obtain the parameters of the cross interactions. For the TIP4P/2005 model positive deviations from the energetic Lorentz-Berthelot rule were also considered to indirectly account for the polarization of methane when introduced in liquid water. To locate the three-phase coexistence point, two different global compositions were used, which yielded (to within statistical uncertainty) the same predictions for the three-phase coexistence temperatures, although with a somewhat different time evolution. The three-phase coexistence temperatures obtained at different pressures when using the TIP4P/Ice model of water were in agreement with the experimental results. The main reason for this is that the TIP4P/Ice model reproduces the melting point of ice I(h).

Heat capacity of water: A signature of nuclear quantum effects.

In this note we present results for the heat capacity at constant pressure for the TIP4PQ/2005 model, as obtained from path-integral simulations. The model does a rather good job of describing both the heat capacity of ice I(h) and of liquid water. Classical simulations using the TIP4P/2005, TIP3P, TIP4P, TIP4P-Ew, simple point charge/extended, and TIP5P models are unable to reproduce the heat capacity of water. Given that classical simulations do not satisfy the third law of thermodynamics, one would expect such a failure at low temperatures. However, it seems that for water, nuclear quantum effects influence the heat capacities all the way up to room temperature. The failure of classical simulations to reproduce C(p) points to the necessity of incorporating nuclear quantum effects to describe this property accurately.

The phase diagram of water at negative pressures: virtual ices.

The phase diagram of water at negative pressures as obtained from computer simulations for two models of water, TIP4P/2005 and TIP5P is presented. Several solid structures with lower densities than ice Ih, so-called virtual ices, were considered as possible candidates to occupy the negative pressure region of the phase diagram of water. In particular the empty hydrate structures sI, sII, and sH and another, recently proposed, low-density ice structure. The relative stabilities of these structures at 0 K was determined using empirical water potentials and density functional theory calculations. By performing free energy calculations and Gibbs-Duhem integration the phase diagram of TIP4P/2005 was determined at negative pressures. The empty hydrates sII and sH appear to be the stable solid phases of water at negative pressures. The phase boundary between ice Ih and sII clathrate occurs at moderate negative pressures, while at large negative pressures sH becomes the most stable phase. This behavior is in reasonable agreement with what is observed in density functional theory calculations.

The phase diagram of water at high pressures as obtained by computer simulations of the TIP4P/2005 model: the appearance of a plastic crystal phase.

In this work the high pressure region of the phase diagram of water has been studied by computer simulation by using the TIP4P/2005 model of water. Free energy calculations were performed for ices VII and VIII and for the fluid phase to determine the melting curve of these ices. In addition, molecular dynamics simulations were performed at high temperatures (440 K) observing the spontaneous freezing of the liquid into a solid phase at pressures of about 80,000 bar. The analysis of the structure obtained lead to the conclusion that a plastic crystal phase was formed. In the plastic crystal phase the oxygen atoms were arranged forming a body center cubic structure, as in ice VII, but the water molecules were able to rotate almost freely. Free energy calculations were performed for this new phase, and it was found that for TIP4P/2005 this plastic crystal phase is thermodynamically stable with respect to ices VII and VIII for temperatures higher than about 400 K, although the precise value depends on the pressure. By using Gibbs-Duhem simulations, all coexistence lines were determined, and the phase diagram of the TIP4P/2005 model was obtained, including ices VIII and VII and the new plastic crystal phase. The TIP4P/2005 model is able to describe qualitatively the phase diagram of water. It would be of interest to study if such a plastic crystal phase does indeed exist for real water. The nearly spherical shape of water makes possible the formation of a plastic crystal phase at high temperatures. The formation of a plastic crystal phase at high temperatures (with a bcc arrangements of oxygen atoms) is fast from a kinetic point of view occurring in about 2 ns. This is in contrast to the nucleation of ice Ih which requires simulations of the order of hundreds of ns.

What ice can teach us about water interactions: a critical comparison of the performance of different water models.

The performance of several popular water models (TIP3P, TIP4P, TIP5P and TIP4P/2005) is analyzed. For that purpose the predictions for ten different properties of water are investigated, namely: 1. vapour-liquid equilibria (VLE) and critical temperature; 2. surface tension; 3. densities of the different solid structures of water (ices); 4. phase diagram; 5. melting-point properties; 6. maximum in the density of water at room pressure and thermal coefficients alpha and KT; 7. structure of liquid water and ice; 8. equation of state at high pressures; 9. self-diffusion coefficient; 10. dielectric constant. For each property, the performance of each model is analyzed in detail with a critical discussion of the possible reason of the success or failure of the model. A final judgement on the quality of these models is provided. TIP4P/2005 provides the best description of almost all properties of the list, the only exception being the dielectric constant. In second position, TIP5P and TIP4P yield a similar performance overall, and the last place with the poorest description of the water properties is provided by TIP3P. The ideas leading to the proposal and design of the TIP4P/2005 are also discussed in detail. TIP4P/2005 is probably close to the best description of water that can be achieved with a non-polarizable model described by a single Lennard-Jones (LJ) site and three charges.