A Theoretical Study of the Hydration of Methane, from the Aqueous Solution to the sI Hydrate-Liquid Water-Gas Coexistence

Monte Carlo and molecular dynamics simulations were done with three recent water models TIP4P/2005 (Transferable Intermolecular Potential with 4 Points/2005), TIP4P/Ice (Transferable Intermolecular Potential with 4 Points/ Ice) and TIP4Q (Transferable Intermolecular Potential with 4 charges) combined with two models for methane: an all-atom one OPLS-AA (Optimal Parametrization for the Liquid State) and a united-atom one (UA); a correction for the C–O interaction was applied to the latter and used in a third set of simulations. The models were validated by comparison to experimental values of the free energy of hydration at 280, 300, 330 and 370 K, all under a pressure of 1 bar, and to the experimental radial distribution functions at 277, 283 and 291 K, under a pressure of 145 bar. Regardless of the combination rules used for σC,O, good agreement was found, except when the correction to the UA model was applied. Thus, further simulations of the sI hydrate were performed with the united-atom model to compare the thermal expansivity to the experiment. A final set of simulations was done with the UA methane model and the three water models, to study the sI hydrate-liquid water-gas coexistence at 80, 230 and 400 bar. The melting temperatures were compared to the experimental values. The results show the need to perform simulations with various different models to attain a reliable and robust molecular image of the systems of interest.


Introduction
Gas hydrates are compounds formed by the inclusion of gas molecules in cavities of the crystal lattice of water, and they can exist at elevated pressures for temperatures somewhat above the melting point of hexagonal ice (ice Ih) . Depending on the properties of the guest gas molecules and the details of the hydrate formation procedure, different structures can be obtained, the most important being structure I (sI), structure II (sII) and structure H (sH) [1,2]. Typically, smaller gas molecules (such as methane, ethane and carbon dioxide) tend to form sI hydrates, while larger molecules preferentially form sII (propane, iso-butane) and sH (cyclohexane, cycloheptane) hydrates. These structures differ in the size of the cavities in the clathrate network of water molecules, as well as in the number of cavities of different types in the unit cell. The smallest cavity found in gas hydrates is the pentagonal dodecahedral cage (5 12 ) comprising twelve pentagons (average radius < r > = 3.95 Å). Larger cavities include solutes enhanced the structure of water [35]; however, from a more recent comparison of 13 C chemical shifts obtained from magic-angle-spinning nuclear magnetic resonance (MAS NMR) [36] for methane in the hydrate and in the aqueous phase, it was concluded that n H = 20, arguing that this value is indisputable, albeit with a dynamic aqueous methane hydration shell where water molecules might continuously enter and leave the hydration sphere. The same number is reported in a recent molecular dynamics study [33]. Nonetheless, the hydration number n H = 16 determined from neutron diffraction [34] still poses a problem of interpretation.
Theoretical studies of the methane-water systems complement the information gathered from experiments, by providing both interpretations at the molecular level and an inexpensive means to assess the feasibility and even the economic cost of using a certain method to impede hydrate formation. The reliability of the predictions obtained from numerical simulations depends on the accuracy of the molecular models that are employed. To be able to study the formation and the melting of hydrates, these models should ideally perform equally well over a range of thermodynamic conditions ample enough to comprise the three phases present in a pipeline: the gaseous methane, its aqueous solution and the crystalline solid. Unfortunately, this is currently out of the question: no water model exists to date that is capable of describing equally well the ices and the liquid. The best model for the ices, TIP4P/Ice [22], fails to reproduce the equation of state ρ(T) of the liquid, whereas TIP4P/2005 [37] is probably close to the best description of water that can be achieved with a non-polarizable model described with a single Lennard-Jones (LJ) site, and three charges though cannot reproduce the static dielectric constant (T) and produce a too-low melting temperature for ice Ih. This last feature is common to the more recent TIP4Q [38], which improved the agreement with experimental data of the liquid, especially the dielectric constant (T). Though the strongest interactions in hydrates are the same as those for ices, namely hydrogen bonding between water molecules, the size of the cages and, especially, their occupancy also depend on the gas-water interactions [39], methane in this case. The methane molecule can be modeled either considering all of the hydrogens, the all-atom approach (for instance the OPLS-AA [40]) or with an electrically-neutral single site, the united-atom (UA) approximation [30]; in both cases, the interaction with water has been modeled with a standard Lennard-Jones (LJ) potential. Whereas the non-zero charges of the AA approach were obtained from quantum calculations, the parameters of the LJ potential, and σ, were fitted to reproduce the methane-water interaction with a specific water model, the original TIP4P. The use of these methane models with different water models can be done either with combination rules or with a re-parametrization of the methane-water potential. This was the subject of a study by the group of Vega [41], who concluded that a 7% increase in the C,O parameter sufficed to reproduce the solubility of the gas and the properties of the methane hydrate, for the UA methane model combined with the TIP4P/2005 water model. However, in a more recent study of the three-phase coexistence, the same group used the TIP4P/Ice model with the original UA methane, but without the 7% correction [42]. Jensen et al. [43] and Michalis et al. [44] have also calculated this phase equilibrium line quite rigorously for methane hydrates directly from molecular simulation.
It becomes then relevant to compare the predictions of the different models on the behavior of the systems of interest, as this allows one to assess the robustness of the various conclusions that can be attained. Thus, the purpose of this study is to apply the different techniques of numerical simulations to compare the performance of the rigid models of water TIP4P/2005, TIP4P/Ice and TIP4Q in reproducing the experimental data of the hydration of methane. Therefore, in the present work, we present the results of Monte Carlo (MC) and molecular dynamics (MD) simulations of the diluted aqueous solution of methane, the sI hydrate and the methane gas-liquid water-sI hydrate coexistence, performed with the three previously-mentioned water models, combined with the OPLS-AA all-atom model for methane [40] and a more recent united-atom (UA) model [41]. The comparison to experimental data is made with the hydration free energies, the coordination properties, the sI hydrate thermal expansivity and the gas-liquid-hydrate coexistence conditions.

Free Energies of Hydration
A minimum requirement for an empirical model intended to correctly describe the interaction of methane with water is the reproduction of the hydration free energy ∆ hyd G at infinite dilution, ideally at various different temperatures [41]. The very low solubility of hydrophobic molecules, a ratio of about 1/4000 waters under ambient conditions [45], poses a problem for simulations with a much smaller number of water molecules, in this case, the ratio being 1/241. However, the very low energy of the methane molecule with water, and even lower with the other methane molecule, allows one to obtain quantitative agreement with experiments from the higher simulated concentration.
The hydration free energy ∆ hyd G was computed as described in Section 3 for the model combinations in Table 1, and the results are depicted in Figure 1, along with those reported in [41] for the experimental data and the values for combinations 2005-2 and 2005-3. It can be seen that when the Bennett's Acceptance Ratio [46] (BAR) was used, the combinations 2005-1, 2005-2, Q-1 and Q-2 yielded values in very close agreement with the experiment, whereas the TIP4P/Ice water model produced somewhat larger deviations. With this method, all of the estimates for ∆ hyd G with the 7% correction on the C−O interaction were underestimated, with TIP4P/Ice increasing the discrepancy at lower temperatures. As it turns out from using the BAR method, the 7% correction worsens the agreement with experimental data; thus, the conclusion is opposite that in [41]. However, the computation of ∆ hyd G with equal acceptance resulted in agreement with the data that were obtained from MC with the Widom insertion method [47] in [41]. These differences highlight a common problem of all empirical molecular models, viz. the dependence of the parametrization on the methods used to compute the target experimental data. It is worth mentioning that, to the best of our knowledge, none of the more recent simulations on the formation and the melting of hydrates [42,43,[48][49][50] employs any correction to the so-called Lorentz-Berthelot combination rules.  Comparison to the experimentally determined hydration free energy of methane (as reported in Reference [41]) of those obtained from MD simulations with Bennett's Acceptance Ratio [46] (g_bar) and evenly spaced (g_energy) thermodynamic integration. OPLS-AA: All-atom force-field from Reference [40]; UA: United-atom model from Reference [41]; UA (7%): Same model with a 7% correction for the C−Øinteraction. The symbols were given sizes slightly larger than the corresponding standard deviations.

Radial Distribution Functions and Coordination Numbers n H
In the present work, Monte Carlo simulations were performed on a system with one methane molecule in 343 water molecules, which amounts to an order of magnitude larger than the solubility of methane [45]. The sampling was done on the isothermal-isobaric (NpT) ensemble, using isotropic pressure, the analytical model potentials and the thermodynamic conditions that are described in Section 3. A spherical cutoff of 1 nm was used, and long-range interactions were handled with Ewald sums. One MC step comprised 5000 trials, divided into the following fractions: 0.003 for CH 4 moves, 0.994 for H 2 O moves and 0.003 for attempts to change the volume. The molecular displacements and rotations, as well as the volume changes, were adjusted to yield a 50% acceptance: trial ratio. The simulation of each system started from an arbitrary configuration, and an initial run of 3 × 10 4 MC steps was used for equilibration. Production runs comprised 7 × 10 4 MC steps, and their statistical significance was assessed with the blocking method [51], whence an average and a standard deviation were assigned, for instance, to the densities. The standard deviation was the same for all runs, ±2 kg· m −3 . The TIP4P/Ice model systematically produced lower densities, albeit only slightly (no more than 1%).
The methane-water rdfs g COw (r) and g CHw (r) are shown in Figures 2-4, for the TIP4P/2005, TIP4P/Ice and TIP4Q water models, respectively. All models predict a more ample cavity for methane, which is inferred from the ca. 0.25 Å shift to the right of the first peak in both rdfs. In general, the OPLS-AA produced less structure than the UA model. While the 7% correction to C,O flattened the rdfs when used with both TIP4P/2005 and TIP4P/Ice, it had a much smaller (and opposite) effect with TIP4Q. In fact, it can be noticed in Figure 4 that TIP4Q is less sensitive to the choice of CH 4 model.
Because the first minimum of the rdfs does not attain zero in any case, there is not a clear-cut first hydration shell; nonetheless, all of the g COw (r) first minima occur around r = 5.5 Å. Hence, instead of integrating g COw (r), a histogram was made with the number of water molecules at a maximum distance of r = 5.5 Å, sampled each MC step, to estimate the coordination number n H . All histograms turned out to be normal distributions. All of the simulations yield n H ∼ 20 ± 3, in agreement with the MAS NMR data of [34], but with distributions that range from n H = 11 to n H = 30; that is to say very dynamic.

The sI Hydrate
The behavior of the unit cell length of the sI hydrate as a function of temperature that resulted from the simulations described in Section 3 is depicted in Figure 5, and the numerical values are presented in Table 2. It can be seen that the TIP4Q gives the best agreement with the experimental values, whereas TIP4P/2005 has the correct trend, but slightly underestimated in some 0.03 Å. On the other hand, the TIP4P/Ice overestimated the experimental values in some 0.02 Å.  Table 1 and experimental data.

The Gas-Liquid-Hydrate Coexistence
This three-phase coexistence has already been studied with different molecular models [42][43][44]48,50]. Because TIP4P/Ice yields the best reported reproduction of the phase diagram of the ices [22], the combination Ice-2 in Table 1 was used in [42,43,50], but somewhat different results were obtained, which have been ascribed to the different area of the contact surfaces. This discrepancies corroborate the observation made in Section 2.1 about the dependence on the simulation methods used to compute the target values, of the parametrizations of empirical molecular models, thus supporting the main idea of the present work, that different models have to be used to attain reliable molecular images of the systems of interest. Thus, the two water models TIP4P/2005 and TIP4P/Ice are used in this work for a comparison with those previous studies, and the TIP4Q potential is added to check on its performance, all combined with the UA model methane.
The first set of simulations was made at 230 bar, and the evolution in time of the potential energy at various different temperatures is shown in Figure 6 for TIP4P/2005; in Figure 7, for TIP4P/Ice, and in Figure 8, for TIP4Q, all with System A. In Figure 9, we show the results for the three models, all with System B. The resulting three-phase coexistence temperatures, T 3 , are presented in Table 3 and in Figure 10; a good agreement was found with the data in [42,44] for the TIP4P/2005 and TIP4P/Ice models, and the value obtained for TIP4Q is close to that of TIP4P/2005. It can be seen that the three water models give the same coexistence temperature for System A and for System B at 230 bar, with the difference that System B crystallizes faster than System A with the same water potential. For example, using the TIP4Q water model, System B crystallizes in around 50,000 ps, while System A crystallizes in 150,000 ps. With the TIP4P/Ice water model, System B crystallizes in around 40,000 ps, while System A in 50,000 ps. Additionally, with the TIP4P/2005, System B crystallizes in around 40,000 ps, while System A in 150,000 ps. Opposite the findings in [52], the simulation of a larger system affected solely the rate at which crystallization occurred, but not the predicted coexistence temperature. It can be seen that both the TIP4Q and TIP4P/2005 models underestimate the experimental values of T 3 for all pressures. This is more clearly seen in Figure 11, where the logarithm of the pressure has been plotted as a function of temperature. The experimental data have been taken from [1]. The TIP4P/Ice water model is the best of the three models considered in this study, in reproducing the coexistence temperature for the range of pressures studied in this work. The direct coexistence method possesses an inherent degree of stochasticity [44,53], which is why our result of the TIP4P/Ice model at 400 bar is slightly different from the results of Conde    The method used in this work to determine T 3 has been criticized [54] because it considers melting as an isothermal process, whereas it is more closely adiabatic in the real system, with significant spatial and temperature gradients. While taking them into account does modify the rate and the mechanism of decomposition [24,52,54], no critical effect on T 3 has been reported [70].

Methods
In the simulations in this study, we used the rigid/non-polarizable TIP4P/Ice, TIP4P/2005 and TIP4Q water potentials. The TIP4P/Ice and TIP4P/2005 water models have an LJ interaction site located on the oxygen atom, positive charges located at the positions of the H atoms and a negative charge located at a distance d OM from the oxygen along the H-O-H bisector; whereas the TIP4Q water model has an LJ interaction site located on the oxygen atom, positive charges located at the positions of the H atoms, a positive charge located at the position of the O atom and a negative charge located at a distance d OM from the oxygen along the H-O-H bisector. The parameters of all of the molecular models used in this work are shown in Tables 4 and 5. When the OPLS-AA model was used for methane, the Lennard-Jones parameters for the C−O and H−O interactions resulted from the following geometric averages: which is the default combination for the OPLS-AA force-field, whereas with the UA methane, the arithmetic mean was used for σ i,j = 1 2 (σ i,i + σ j,j ). Furthermore, the 7% correction recommended in [41] was also applied to C,O , thus yielding nine model combinations that will henceforth be referred to as labeled in Table 1.
The GROMACS 4.5.1 package [55][56][57][58] was used for MD simulations to compute the hydration free energy of methane with the thermodynamic integration method [46,[59][60][61][62] built in it. A system with one methane molecule in 241 water molecules was simulated under constant pressure and temperature (NpT ensemble), at 1 bar, controlled with the Berendsen barostat and at temperatures of 280, 300, 330 and 370 K, controlled by using the stochastic dynamics (sd) integrator with a 2-fs time-step, as described in the GROMACS user manual [63] (for the thermodynamic integration and BAR [46] methods, see Section 3.12.2 in user manual).  4 3.73 147.5 0.0 A homemade Monte Carlo program was used with a system of one methane molecule in 343 water molecules, to simulate the dilute aqueous solution of methane under a pressure of 145 bar and at 277, 283 and 291 K, to obtain the coordination numbers n H and the radial distribution functions (rdfs) and to compare them to the experimentally-determined data under the same conditions [34,36]. This was done because in Monte Carlo, the pressure and temperature controls are exact, and the structural information can be obtained somewhat more readily.
The GROMACS 4.5.1 code was also used to simulate the fully-occupied sI hydrate in the NpT ensemble: the unit cell, with a side of length 1.203 nm, was built according to the X-ray crystallographic data [64]; the hydrogen atoms of the water molecules were distributed randomly, but following the Bernal-Fowler rules and changing the orientations until achieving a near-zero total dipole moment. The unit cell was then replicated 2 times in each orthogonal direction (2 × 2 × 2) to form a cubic cell of side 2.406 nm, constituted by 368 molecules of water and 64 of methane. The pressure was kept at 30 bar by means of the isotropic Parrinello-Rahman barostat [65,66] with time constant for coupling tau_p = 2, and the temperature was successively fixed at 175, 200, 225, 250 and 270 K with a Nosé-Hoover thermostat [67,68]. For the long-range Coulombic interaction, the Particle Mesh Ewald (PME) algorithm was used with a cut-off radius of 0.9 nm; an LJ interaction was implemented with a cut-off radius of 0.9 nm; and the Lorentz-Berthelot mixing rules were implemented in both cases. The simulations were implemented as NpT MD simulations using three-dimensional periodic boundary conditions. The system was run for 200 ps at each temperature, with a 3-fs time-step, to compute the thermal expansivity of the sI hydrate model. This time span was proven to yield statistically-meaningful averages for the unit cell length in [41].
The same MD code with the direct phase coexistence method [69] was used to calculate the three-phase gas-liquid-hydrate coexistence temperature. We used the same time-step, the same thermostat and the same barostat, but allowing each of the three orthogonal directions to fluctuate independently, as described in [42]. It is worth mentioning that this procedure has been acknowledged to lead to an appropriate computational prediction of the phase diagram [70]. We made the coexistence analysis on two systems; System A comprised the same initial sI hydrate in contact with two other cubic boxes of the same size, one with 368 waters in the liquid phase to one side and another with 64 methanes in the gas phase to the other side, that yielded a computational cell of size 2.406 nm × 2.406 nm × 7.218 nm, with the contact interfaces perpendicular to the Z-axis. Additionally, System B was comprised of System A replicated 2 × 1 × 1 times; this means that the hydrate phase of System B comprised 736 water molecules and 128 methane molecules; the liquid phase comprised 736 water molecules; and the gaseous phase comprised 128 methane molecules; which yielded a computational cell of size 4.812 nm × 2.406 nm × 7.218 nm ( Figure 12). To equilibrate the initial configurations with each model, a short 50-ps simulation was performed at 250 K under pressures of 80, 230 and 400 bar for System A and 230 bar for System B, which did not result in either melting or crystallization. Temperature scans were then performed at 80, 230 and at 400 bar for System A and 230 bar for System B, using the same controls and time-step for all simulations of each model. The systems were simulated at each temperature while analyzing the evolution of the potential energy as a function of time. The increase in potential energy indicated melting, whereas its decrease indicated crystallization. The equilibrium phase coexistence temperature was taken as an average of the lowest temperature at which the hydrate melted and the highest temperature at which the system froze, as in [42,48].

Conclusions
A series of MC and MD simulations were performed to calibrate the predictions that can be obtained from different models on the hydration of methane. Whereas TIP4P/Ice has a better performance for the ice phase diagram, TIP4P/2005 and TIP4Q provide a better description of liquid water; thus, it was no surprise that these two latter produced better agreement with the experimental data of the diluted aqueous solution of methane. The hydration free energy ∆ hyd G was computed in this work with a more robust algorithm than in previous studies, which resulted in agreement with experimental data for all six combinations of water-methane models with their original parameters, regardless of the combination rule used for σ C,O , opposite the conclusion in [41]. Finally, we have performed molecular dynamics simulations to estimate the three-phase (methane hydrate-water-methane) coexistence temperature T 3 at three different pressures (80, 230 and 400 bar) by using the direct coexistence method and the three water models (TIP4P/Ice, TIP4Q and TIP4P/2005) with the UA model methane. The results showed that the three-phase coexistence temperatures obtained with the TIP4P/Ice model were in good agreement with the experimental data. Results obtained by TIP4P/2005 and TIP4Q models were shifted to lower temperatures by about 20 and 25 K, respectively, with respect to the experimental data. A caveat is in order, as the method used to determine the coexistence temperatures relies solely on the behavior of the potential energy, disregarding other criteria to track the formation/melting of the hydrate. To be on the safe side, the final configurations of System A at P = 230 bar and two temperatures, T = 290 K and T = 300 K, are shown in Figure 12, where it can be seen that the former resulted in the complete formation of the hydrate, whereas the latter yielded a liquid-like geometry. We have observed that System A and System B give the same coexistence temperature and that System B crystallized faster than System A using the same potential and the same thermodynamics conditions; this could mean that using the direct coexistence method with bigger systems could result in better and faster results. The unexpected better performance of TIP4Q with regard to the thermal expansivity of the fully-occupied sI hydrate in the vicinity of the melting conditions suggests that the ability of TIP4P/Ice to produce a higher melting temperature of the hydrate, and perhaps also of ices, is related to the lower density of the model.