Phase stability of the ice XVII-based CO2 chiral hydrate from molecular dynamics simulations

We computed the phase diagram of CO 2 hydrates at high pressure (HP), from 0.3 to 20 kbar, by means of molecular dynamics simulations. The two CO 2 hydrates known to occur in this pressure range are the cubic structure I (sI) clathrate and the HP hydrate, whose water frame-work is the recently discovered ice XVII. We investigated the stability of both hydrates upon heating (melting) as well as the phase changes upon compression. The CO 2 -filled ice XVII is found to be more stable than the sI clathrate and than the mixture of ice VI and dry ice at pressure values ranging from 6 to 18 kbar and in a wide temperature range, although a phenomenological correction suggests that the stability should more realistically range from 6.5 to 13.5 kbar. Our simulation results support the current hypothesis that the HP hydrate is stable at temperatures above the melting curve of ice VI.


I. INTRODUCTION
Clathrate hydrates, also known as gas hydrates, are a class of host-guest materials where guest molecules reside within, and help stabilize, crystalline cages formed by the hydrogen-bond network of water molecules. 1 Naturally occurring clathrate hydrates, mainly biogenic methane hydrate, are present in large amounts on the sea floor (estimates vary between 100 and 63.000 Gt of carbon 2 ), where methane hydrate is stable below 300 m, 3 making it possibly the largest source of hydrocarbon on Earth. While clathrate formation is of high practical relevance as it is responsible for severe pipeline clogging, 4,5 much of the interest in clathrates revolves around the possibility of using them as energy source, 6 as well as a possible means of CO 2 sequestration and long-time storage, 7 whereas there are still important open questions regarding the role of clathrates in the climate change feedback. 3,8 Hydrogen hydrates have also been used as alternative routes to obtain new ice polymorphs, as in the case of ices XVI and XVII, 9,10 showing how these structures keep prompting new exciting findings, more than 200 years after their discovery. 11,12 Since the first determination by von Stackelberg and Müller of the crystal structure of a clathrate hydrate using X-Ray diffraction, 13 several additional structures have been identified, the most common being denoted as the cubic structure I (sI), 14 the cubic structure II (sII), 15 and the hexagonal structure (sH). 16 In most cases, clathrates that are in the sI or sII form, at low pressures, transform upon compression to the sH and, eventually, to the methane hydrate-III form (MH-III), 17 possibly passing through the tetragonal structure (ST). 18 In 2010, however, Hirai and co-workers found that at pressure higher than 6 kbar, the CO 2 clathrate hydrate, which starts at low pressure in the sI form, changes into a previously unknown form of hydrate (named by Hirai high-pressure phase, or HP 19 ), only to decompose into dry ice and ice VI above 10 kbar. A series of experiments 10,20 showed that the water network associated with these new types of hydrates is in fact ice XVII (also known as Sχ) and is characterized by open, helicoidal channels, which give the structure a chiral nature.
The question of the stability field boundaries of this new type of hydrate is still open, and computer simulations can provide a viable alternative [21][22][23][24][25][26][27][28][29][30] to probe those regions of the phase diagram, which are difficult to investigate experimentally. In this work, we employ

II. METHODS
The water and CO 2 molecules are modeled, respectively, using the TIP4P/ice potential developed by Vega and co-workers to reproduce the coexistence lines of different ice forms 31 and the TraPPE potential developed in 2001 by Potoff and Siepmann to describe mixtures containing alkanes, CO 2 , and nitrogen. 32 Both models treat the molecules as rigid, planar (water), or linear (CO 2 ) objects and make use of interaction centers located not only at the atomic positions. The TIP4P/ice model uses 4 interaction centers, namely, one oxygen atom OW, two hydrogen atoms HW, and one virtual site MW located along the bisector of the HW-OW-HW angle.
To model the CO 2 molecule, we make use of five interaction centers, namely, three massless (virtual) sites bearing charge and Lennard-Jones interactions (one carbon atom C and two oxygen atoms O) and two real sites X bearing no interaction but the mass needed to reproduce the inertia tensor of the original three-site model. The use of virtual sites is necessary to simulate the linear molecule using Lagrange multiplier-base methods, which would be otherwise not possible given the linear dependency of the C-O and C-C vectors.
The parameters of the models are summarized in Table I. All simulations were performed using the GROMACS molecular dynamics simulation package, version 2018.2, 33,34 simulating the isothermal-isobaric ensemble by imposing the temperature and pressure equilibrium values by means of the Nosé-Hoover thermostat 35,36 (1 ps relaxation time) and Parrinello-Rahman barostat 37 (2 ps relaxation time), respectively. The pressure coupling scheme adopted is the anisotropic one, where the diagonal elements of the pressure tensor are relaxed independently from each other, if one solid phase is present in the system. When simulating the water/CO 2 liquid mixture, we adopt the slab configuration, and in this case, we keep the box edges parallel to the liquid-liquid  The equations of motion were integrated using the leapfrog algorithm with a time step of 1 fs, and the molecular structure of water and CO 2 was kept rigid using the SETTLE 38 and LINCS 39 algorithms, respectively. Both the electrostatic and the dispersion interactions were calculated using the smooth variant of the Particle Mesh Ewald (sPME) method 40 using a real space mesh with a spacing of 0.15 nm and a relative accuracy of the potential energy, calculated at the cutoff of 1 nm, of 10 −5 and 10 −3 for the Coulomb and dispersion interactions, respectively. The initial structures of sI hydrate, ice VI, and ice XVII were generated using the Genice program, 41 taking care of rearranging the position of the hydrogen atoms in the crystalline lattice in order to satisfy (as much as possible) the ice rules, as both ice VI and ice XVII are hydrogen-disordered. The single phase simulations, namely, sI clathrate, HP hydrate, ice VI, and CO 2 ice I, were composed of 432, 493, 640, and 864 molecules, respectively. For the explicit coexistence simulations, instead, we employed larger systems, composed of 2300, 2238, and 2787 molecules for the sI/liquid, HP/liquid, and ice VI/liquid cases, respectively. In Fig. 1, we report snapshots from a trajectory of the CO 2 sI/liquid and HP/liquid interfaces, as well as of the sI and HP hydrates. The water/CO 2 ratio used in the simulation of these systems is discussed later on in the text.
It is worth noticing that, since the structure of ices and liquid water is mainly determined by hydrogen bonds, the effect of the long range part of the dispersion interaction is expected to play a minor role for ice, similarly to what happens in liquid water. 42,43 The effect of the long-range part of the dispersion interaction between the CO 2 molecules (which have only a quadrupole moment) and water, on the contrary, is expected to be relevant because CO 2 hydrates are known to be stabilized by the CO 2 /water dispersion interactions. The combination of TIP4P/Ice and TraPPE potentials (with cutoff Lennard-Jones interaction) has in fact been used previously to describe the CO 2 sI clathrate hydrates by Vega and co-workers 28 who, in order to reproduce the experimental melting line of the sI hydrate, modified the mixing rules that determine the Lennard-Jones interaction between water and CO 2 . As will be shown in the discussion of the results, the inclusion of the long-range part of the dispersion interaction moves the melting line of the CO 2 sI clathrate closer to the experimental values without resorting to modifying the mixing rules.
To map the phase boundaries between the different phases taken into consideration (sI, HP, ice VI, CO 2 ice I, liquid CO 2 /water), (i) we determined the melting temperature of hydrates sI and HP, and that of ice VI, in contact with the CO 2 -saturated liquid, by means of the three-phase explicit coexistence method; 44 (ii) starting from the set of melting points that could be determined more precisely, we followed the corresponding solid/liquid coexistence lines by means of the Gibbs-Duhem integration method; 45 (iii) from the crossing of pairs of melting curves (sI/liquid, HP/liquid, ice VI/liquid), we identified the corresponding quadruple points; (iv) starting from the quadruple points, we performed additional Gibbs-Duhem integrations moving away from the liquid phase, that is, along the HP/sI/CO 2 ice I and HP/ice VI/CO 2 ice I coexistence lines, thus completing the tracing of the phase boundaries. of Chemical Physics

III. RESULTS
To determine suitable starting points for the Gibbs-Duhem integration, we performed a series of simulations at different pressures and temperatures of the explicit interface between the CO 2saturated liquid phase and either of the hydrates or ice VI. The explicit three-phase interface is realized by juxtaposing, in slab configuration, a fragment of the solid phase of interest (sI, HP, or ice VI) to a pre-equilibrated, partially mixed solution of water and liquid CO 2 , at, or close to, the thermodynamic point of interest. This procedure is needed because the kinetics of CO 2 mixing could be slower than the onset of melting or crystallization at the boundary of the solid phase. It is particularly important to equilibrate the saturated solution when the pressure is modified, as the miscibility gap depends more pronouncedly on the latter rather than on the temperature, in the range of parameters considered here.
Points on the melting curve can be estimated with good accuracy by identifying the temperature at which the time evolution of the potential energy switches from an increasing trend (melting) to a decreasing one (crystallization). With simulation runs of up to 50 ns (the closer to the melting point, the longer the simulation needs to be in order to notice any trend), we were able to determine the melting temperatures with an accuracy of up to 1 K by interpolating the slope of the energy trend in the initial 20 ns. In Fig. 2, we report the time evolution of the potential energy of the hydrate HP at 10 4 atm, blockaveraged over consecutive windows of 150 frames, for four different temperatures (280, 290, 300, and 305 K), plus that of control run at 295 K. The lower panel of Fig. 2 shows the slope of the energy trend obtained from a linear fit of the energy as a function of temperature, which is used to improve the accuracy of the coexistence temperature estimate. The slope of the energy drift is becoming zero at the coexistence line; therefore, it is possible to estimate the latter by an interpolation method of choice. We use simply the equation of the straight line passing through the points at 290 and 300 K. Error propagation allows us to estimate a coexistence temperature of 295 ± 1 K. The result from a successive control run at 295 K shows that this is a good estimate as the energy trend at this temperature has a negligible slope. Furthermore, to check that the system is indeed en route to either the melting or the crystallization, we prolonged the runs until the full transition was seen (305, 300, and 280 K at around 15, 70, and 60 ns, respectively). We did not observe a complete crystallization at 290 K, but expect it at about 100 ns, whereas the control The coexistence temperature can be estimated as the one where the interpolated slope crosses zero. The square represents the slope at the control temperature of 295 K.

ARTICLE
scitation.org/journal/jcp system at 295 K did not show any appreciable change in the fraction of crystal and liquid phases. This method requires, evidently, very long simulation runs, and in some cases, either because of thermal fluctuations or because of a very slow kinetics of crystal growth, there could be no clear-cut difference between an increasing and a decreasing trend. For this reason, we decided to use the melting points that were determined more clearly (10 3 and 10 4 atm for sI and HP, respectively) as starting points for a Gibbs-Duhem integration. 45 The coexistence line between hydrates sI or HP and the CO 2 -saturated solution is traced by integrating the Clapeyron equation in the form using the classical (fourth order) Runge-Kutta method, where δh and δv represent the molar enthalpy and the volume difference between the solid phase (either sI or HP) and the liquid phase. For each value of pressure and temperature, we performed short simulation runs of the homogeneous solid phase with full anisotropic pressure coupling and of the liquid phase in a slab configuration with a CO 2 -containing, water-rich phase and a water-containing, CO 2 -rich phase, with pressure coupling applied only in the direction normal to the interface. Each simulation consisted of 200 ps of relaxation and 200 ps of equilibrium sampling, which was long enough to determine enthalpy and volume within 3% and 1%, respectively. In addition, along the coexistence line, we sampled the number of CO 2 molecules per water in the water-rich, CO 2 -containing liquid, which is reported in Fig. 3 as a function of pressure, and compared with the phenomenological fit of Diamond and Akinfiev. 46 The number of water molecules per CO 2 in the CO 2 -rich phase was always below 1/1000. The global composition in the liquid phase has been set equal to the composition of the corresponding hydrate, in order to compute directly the quantity δv/δh from each pair of simulations using a weighted difference of the average volume and enthalpy of the two simulation boxes. At the chosen molar fraction, the liquid phase was always within the miscibility gap; therefore, it always maintained a separation between a water-rich, CO 2 -containing phase and a CO 2rich, water-containing one. Note also that as a starting configuration for every new thermodynamic point along the coexistence line, we always used the last configuration generated at the previous one. In this way, the 200 ps relaxation phase was enough to ensure that the local composition in the two phases (water-rich and CO 2 -rich) could reach its equilibrium value.
As hydrates are known to be stable only in a very narrow compositional range, 47,51,52 here we used the working hypothesis that the hydrates are always filled with CO 2 (5.75 and 3.52 water molecules for every CO 2 molecule for the sI and HP hydrates, respectively). In addition, the solubility of CO 2 in ice VI and that of water in CO 2 ice I were not considered. To trace the melting of ice VI, we used a slightly different procedure as the relevant coexistence is that between ice VI, the saturated liquid, and CO 2 ice I. Therefore, at each state point, we simulated these three phases separately and computed the ratio δv/δh using the appropriate stoichiometric coefficients. In Fig. 4, we report the result of the Gibbs-Duhem integrations, as well as some control points computed by means of the direct coexistence, 44 along with the respective error bars. The part of a melting curve drawn with a continuous line marks the region where the corresponding solid phase is the most stable, while the dashed part marks the region where it is not. The quadruple points are found at the crossing between the melting curves (sI and HP, and HP and ice VI). As these points are always located below the melting line of CO 2 ice I (dashed red line 53 ), any excess CO 2 will be found in the CO 2 ice I state. Therefore, at the two quadruple points, one has coexistence of sI, HP, CO 2 ice I, and saturated liquid, and HP, ice VI, CO 2 ice I, and saturated liquid, respectively. Starting from the quadruple points and removing one phase (in our case, the saturated liquid) gives the system one degree of freedom, which can be used to trace the coexistence between the remaining three solid phases (sI/HP/CO 2 ice I and HP/ice VI/CO 2 ice I). In this case, we applied the Gibbs-Duhem scheme by using the temperature as the independent variable and calculating from independent simulations of the three solid phases the ratio δh/δv, again by using the appropriate stoichiometric coefficients. In this way, we were able to locate the boundaries of maximum stability of the sI and HP clathrates, which are reported by the shaded areas in Fig. 4. Note that because the system has two components, it is still possible to find, depending on the composition, two coexisting solid phases within one region of stability. For example, if the total ratio of water to CO 2 molecules is larger than 3.52 at 10kbar, then the HP hydrate will coexist with ice VI (below its melting curve) or, in the case of a smaller ratio, with CO 2 ice I.
In Fig. 4, lower panel, we report a detail of the stability field of the hydrates, overlaid with the experimental data from Refs. 19 and 47-50, as well as the experimental (dashed-dotted thin lines) and TIP4P/ice (thin dashed lines) phase boundaries of water and of ice polymorphs. 31 As already mentioned, the introduction of long range dispersion forces brings the computed melting temperature closer to the experimental one, as compared to the values computed with the cutoff Lennard-Jones potential (inset of Fig. 4 The other melting curve which is experimentally known is that of ice VI in equilibrium with the CO 2 -saturated liquid. 47 In this case, the simulation results show a depression of the melting temperature with respect to the pure water case, which is consistent with experimental observations. 47 Between the stability fields of ice VI and sI clathrate lies that of the HP hydrate, which is delimited by its melting curve at a high temperature and by the coexistence curves with the sI clathrate at low pressure and with that of ice VI/CO 2 ice I at high pressures, respectively. The stability field of the HP hydrate is the experimentally least well-known. The high-and low-pressure boundaries are known to be weakly temperature-dependent, 19,47 but experimental evidence of its high-temperature stability limit is missing, with no HP hydrate found at temperatures higher than 280 K, although it is believed, according to thermodynamic arguments and because of the topology of the phase diagram, that its melting curve should be located at temperatures slightly higher than the melting curve of ice VI. 47 The present simulation results confirm that the solid-solid coexistence curves show only a weak dependence on the temperature and support the hypothesis that the HP hydrate is stable in a temperature range that extends above the melting curve of ice VI. The large temperature gap between the melting point of ice VI and that of the HP hydrate and the corresponding width in pressure of the stability field of the HP hydrate are, however, in quantitative disagreement with the corresponding measurements (stability range between 6 and 10 kbar) or expectations (HP hydrate melting line close to that of ice VI). 47 However, one should consider that the pressure shift (500 bars) with respect to experimental values of the sI hydrate melting curve (which could presumably be the same for the HP hydrate) is opposite in sign to that of the melting curve of ice VI. The latter could be brought close to the experimental values by shifting the temperature by 15 K and the pressure by about −1.3 kbar. By applying these shifts to the simulation results, the gap between the ice VI and HP hydrate melting curves reduces noticeably (the maximum gap at the quadruple point decreases from 32 to 15 K) and the stability range in pressure changes from 6-18 kbar to 6.5-13.5 kbar, bringing the last figures more in line with the experimental findings.

IV. CONCLUSIONS
By simulating ternary equilibria between the CO 2 sI clathrate, the CO 2 -filled ice XVII (HP hydrate), ice VI, CO 2 ice I, and the CO 2water solution at saturation, we traced the boundaries of the stability fields for the two CO 2 hydrates at high pressure. The melting curve of the sI clathrate, simulated by taking into account the longrange contribution to the dispersion forces, is shifted by 8 K and 500 bars with respect to the well-determined experimental values. Ice VI in equilibrium with the CO 2 -saturated liquid shows, in accordance with experimental findings, a depression of its melting temperature with respect to the one of the neat system. From the crossing between these two melting curves and that of the HP hydrate, we determined the stability limits of the HP hydrate at high and low pressure values, which extend from 6 to 18 kbar and up to 308 K. This pressure range can be corrected in order to take into account the fact that the melting curves of the sI clathrate and of ice VI are displaced by different amounts with respect to the experimental data. The correction brings the pressure stability range in the range 6.5-13.5 kbar, up to the maximum temperature of 310 K. Independent from the application of this phenomenological correction, our simulations confirm the hypothesis that the HP hydrate must be stable above the melting curve of ice VI. In addition, we showed that the inclusion of the longrange part of the dispersion interaction improves the estimate of the