The thermal properties of binary structure sI clathrate hydrate from molecular dynamics simulation

ABSTRACT In this work, we present temperature dependence of lattice parameter and normalised lattice parameter in the atmospheric pressure and 120 bar and also pressure dependence of unit cell volume and normalised unit cell volume at 150 and 250 K for variety guests with different size, polarity and guest–host hydrogen bonding capability such as trimethylene oxide (TMO), ethylene oxide (EO), formaldehyde (FA), cyclobutane (CB), cyclopropane (CP) and ethane (Et) in the large cages with CH4 in small cages of sI clathrate hydrates by molecular dynamics simulations. The obtained values of lattice parameters for the guest species are compatible with the experimental values. These clathrate hydrates are simulated with TIP4P/ice four-site water potential. Herein, isobaric thermal expansivity and isothermal compressibility are calculated at a temperature range of 50–250 K and a wide pressure range. These structural properties have been compared for guests which they are isoelectronic and have similar masses but with different size and polarity. We use molecular dynamics simulations to relate microscopic guest properties, like guest–host hydrogen bonding to macroscopic sI clathrate hydrate properties. The temperature dependence of thermodynamic properties such as constant-volume and constant-pressure heat capacity is presented in the atmospheric pressure for these guest species.


Introduction
Clathrate hydrates are non-stoichiometric crystalline inclusion compounds which the host in the crystalline structure contains cavities, holes or channels that guest atoms or molecules with appropriate size are trapped. The structure of clathrate hydrate is formed by the size of the guest molecules which according to this, the clathrate hydrates have three known common structures; structure I (sI), structure II (sII) and structure H (sH) and five cavities, pentagonal dodecahedron (5 12 ), tetrakaidecahedron (5 12 6 2 ), hexakaidecahedron (5 12 6 4 ), irregular dodecahedron (4 3 5 6 6 3 ) and icosahedron (5 12 6 8 ) [1].
The knowledge of the structural and thermodynamic properties of clathrate hydrates such as thermal expansion, isothermal compressibility coefficient and heat capacity are important for increasing our understanding of intermolecular interactions between guest-host and guest-guest molecules and knowledge of these properties is needed to understand the details of the spectroscopic studies. The research results have been shown most of the properties of clathrate hydrate such as Raman and infrared spectra [2,3] are similar to those of ice I h due to the similar hydrogen bonded but clathrate hydrate and hexagonal ice are different in some of the properties such as thermal conductivity and thermal expansivity [4]. The thermal conductivity of hydrates is lower than that of ice Ih. The guest molecules in cages of clathrate hydrate have interaction with the host water lattice that this guesthost interaction is responsible for vibrations of guest molecules with a large anharmonic. Anyway many studies were linked that demonstrate the thermal conductivity of gas hydrates depend on the complicated nature of the hydrate lattice and guest type [5][6][7].
Tse et al. have compared the thermal expansivity of different guests in clathrate hydrate with ice I h and have reported the thermal expansivity is much larger than that of ice I h especially below 200 K due to the presence of guest molecules in clathrate hydrate and guest-host interactions between them [8][9][10][11][12]. The guest molecules size and shape, the rotational motions of guests, water cage occupancy percentage and guest-host hydrogen bonding capability affect the lattice parameter and thermal expansivity. Many experimental research [13][14][15] and theoretical studies have been carried out in order to investigate effect of type and size guest and number of the guest molecules in the cages on the lattice parameter and thermal expansivity clathrate hydrate. The lattice parameters and thermal expansivity have been shown for four Structure I (C 2 H 6 , CO 2 , 47% C 2 H 6 + 53% CO 2 , and 85% CH 4 + 15% CO 2 ) and seven Structure II (C 3 H 8 , 60% CH 4 + 40% C 3 H 8 , 30% C 2 H 6 + 70% C 3 H 8 , 18% CO 2 + 82% C 3 H 8 , 87.6% CH 4 + 12.4% i-C 4 H 10 , 95% CH 4 + 5% C 5 H 10 O, and a natural gas mixture) by experimental measurement [13]. Hester et al reported that there is a general correlation between in the thermal expansivities for a given hydrate structure. The thermal expansion for clathrate hydrates of Ar, Kr and C 3 H 8 of cubic sII, Xe and CH 4 hydrates of cubic sI and also for empty lattices of sI and sII calculated within framework of lattice dynamics approach in quasiharmonic approximation [16]. The MD simulation performed for measurement the lattice parameter and thermal expansivity of the binary THF and H 2 sII clathrate hydrate [17], the CH 4 and CO 2 sI clathrate hydrate with using TIP4P/2005 and TIP4P/ice water models [18,19] and large molecule guest substance (LMGS) like NH in sH clathrate hydrate [20] to investigation effect the size of guest molecules and guest molecule occupancies on the thermal expansion. These studies show clathrate hydrate formed by guest molecules with large van der Waals radius lead to expanding of the host lattice and also the thermal expansion coefficients for these molecules are less than hydrate including small guest molecule. The larger the size of the guest molecules, the vibrational anharminicity and thermal expansivity decrease. However, the effect of hydrogen bonding of the guest molecules and cage water molecules are not fully investigated on the lattice parameter and thermal expansivity. The isothermal compressibility is one of the guest-dependent structural properties. The powder neutron diffraction be used for determining structural properties of N 2 , O 2 , and air clathrate hydrate such as isothermal compressibility that have specified significant differences in the compressibilities of N 2 and O 2 clathrate hydrates [21]. The unit cell parameters and isothermal compressibilities for O 2 -clathrate and N 2 -clathrate that N 2 and O 2 molecules have almost identical size appear to be determined not only by the size and shape of the guest molecules but also by the strength of molecular interaction. The unit cell volume and isothermal compressibility for The O 2 -clathrate are larger than N 2 -clathrate that due to the O 2 -clathrate has stronger guest-host interaction. The temperature dependence of polyhedral cage volumes of trimethylene oxide (TMO) sI and sII clathrate hydrates are calculated from atomic positions determined from neutron powder-diffraction data. This study reveals the guest-host interaction in certain cases, and differences between the dodecahedral cages (small cages) of sI and sII containing TMO molecule [22]. Molecular dynamics simulations have been performed on sI methane hydrate for considering the effect of the amount of methane contained in the hydrate. These MD simulations determine that the presence of methane increases slightly the value of the unit cell and decreases slightly the compressibility of the structure [23]. However, studies of thermal expansivity, and isothermal compressibility of the guest species that are different in the ability to form guest-host hydrogen bonding with water lattice by the MD simulation method are fewer. Alavi et al. determined the effect of the guest species (cyclic ethers) like THF, 1,3-dioxolane, THP, and p-dioxane in sII clathrate on the structure, thermal expansivity, and isothermal compressibility [24] but this effect has not been studied for three-member and four-member heterocyclic ethers and small molecules with carbonyl functional groups in sI clathrate hydrate.
The heat capacity of clathrate hydrate is the most important property. The heat capacity of sI ethylene oxide clathrate hydrate [25], the pure and KOH-doped argon and tetrahydrofuran clathrate hydrates [26,27] and structure I hydrate of trimethylene oxide have been measured by experimental methods [28]. However, measuring the heat capacity of clathrate hydrate by calorimetric method is also rather difficult. We recently studied the probability of guest-host hydrogen bonding formation trimethylene oxide (TMO), ethylene oxide (EO) and formaldehyde (FA) as polar guests compared with non-polar guest molecules with similar structure cyclobutane (CB), cyclopropane (CP) and ethane (Et), respectively for TMO, EO, and FA using MD simulations with the SPC/E and TIP4P/ice force field [29]. The guest-water hydrogen bonding probability increases with temperature in all of guests and at each temperature, the TMO molecules have the largest probability of hydrogen bonding with cage waters. In this paper, we show that the guest-host interaction plays an important role in the determination of the macroscopic properties of gas hydrates at high temperature.
Herein, in order to study the effect of the size and nature of guest molecules on the structural properties of clathrate hydrate, we plot lattice parameter versus temperature and pressure for guest species such as TMO, EO, FA, CB, CP and Et in large cages and CH 4 help gas in small cages by molecular dynamics simulation on a wide range of pressure and temperature. Using these plots, we calculate thermal expansivity and isothermal compressibility of the clathrate phase for the guest molecules. The heat capacities at constant volume (C V ) and constant pressure (C P ) are estimated respectively from MD simulations in NVT and NPT ensembles at 50-250 K.

Computational details
In the molecular dynamics simulations, the initial coordinates of the water oxygen atoms in the sI clathrate hydrate cubic unit cell (a = 12.09 Å) are taken from X-ray crystallography [30] and the initial positions of the water hydrogen atoms about the oxygen atoms were assigned to be consistent with the ice rules. In the simulations, 3 × 3 × 3 replica of the sI clathrate hydrate unit cell with 46 water molecules is used. The initial guest molecule structures were determined by quantum chemical structural optimisation at DFT level using (B3LYP parameterisation) with the 6-311++G(d,p) basis set using the Gaussian 09 suite of programs [31]. The intermolecular van der Waals potentials between atoms i and j on different molecules are modelled as a sum of Lennard-Jones (LJ) and electrostatic point charge interactions: where σ ij and ε ij are the distance and energy parameters of the ij pair of atoms separated by a distance of r ij , and q i and q j are the electrostatic point charges on the atoms. For potentials between unlike atoms, standard combination rules, 1 ij = (1 ii 1 jj ) 1/2 and s ij = (s ii + s jj )/2, are used. Partial electrostatic charges on the atoms of the structurally optimised guest molecules were calculated using the 'charges from electrostatic potential grid' (CHELPG) method [32]. The values for the parameters σ ii , ε ii and q i for selected atom types used in the simulations are given in Table S1 of the Supporting Information (SI). The large cage guest molecules (TMO, EO, FA, CB, CP) were modelled with the general AMBER force field (GAFF) [33] and the Murad-Gubbins potential [34] and the OPLS force field [35] were chosen for methane and ethane molecules, respectively. The water molecules in the clathrate were modelled using the TIP4P/ice four-charge model [36].
The molecular dynamics simulations on periodic simulation cells were performed using the DL_POLY software program version 2.18 [37] with pressure and temperature regulated using the modified Nosé−Hoover thermostat-barostat algorithm [38][39][40] and thermostat and barostat relaxation times of 0.1 and 1.0 ps, respectively. A cutoff of 15 Å is applied to the short-range Lennard-Jones interactions and the electrostatic interactions are evaluated through the Ewald sum method. In all MD simulations a Verlet leap-frog algorithm with a time step of 1 fs was used for the integration of the equations of motion. At the desired temperature, the NPT or NVT simulations with 270,000 steps is used. The first 70,000 steps are used for equilibration and not included in the averaging. During the equilibration stage, the calculated clathrate hydrate structures and energies were seen to converge. To get good statistics in these calculations, at every temperature, at least five NPT or NVT simulations are performed using different starting configurations for the guest molecules and the average results of all simulations have been used. In this way it was ensured that an ensemble average is taken, which provides sufficient accuracy.

Thermal expansivity
In our work, we considered the cell vector length dependence on temperature for TMO, EO, FA, CB, CP and Et as guest molecules with using TIP4P/ice water model from temperatures 50-250 K at ambient pressure and 120 bar that have shown in Figure 1. Experimental lattice parameters for TMO, EO and Et clathrates are also plotted for comparison. The fourmember TMO and CB molecules have a larger molecular size than the three-member ring EO and CP guest molecules and these molecules, in turn, have a larger size than FA and Et molecules. Since the size of guest molecules decreases as ring-CH 2group are replaced by oxygen atom, TMO, EO and FA molecules have smaller size than CB, CP and Et molecules, respectively. In all of case except TMO, we observed that with increase size molecule and temperature, unit cell length increases. At high temperature for the TMO clathrate, the effect of the guest polarity dominates effect of the guest molecule size on the lattice parameter of the clathrate.
As shown in Figure 1, the lattice parameters of the TMO, EO and Et hydrates as a function of temperature at atmospheric pressure have been compared with the experimental values. The lattice parameter obtained for Et is in better agreement with the experimental results. The TIP4P/ice model exhibits larger value of the lattice parameters for TMO and EO hydrate than experimental values.
The slopes of the lattice parameters as a function of temperature appear similar. The slope of these figures is related to thermal expansivity (α) where this property was obtained according to the following definition, where n is the number of moles in a system and a is lattice parameters, Å: The temperature dependence of α is assumed to be: By separation of variables and integration we obtain Equation (4): Expansion coefficients of Equation (3)  (T-T 0 ) for variety guest species in sI clathrate hydrate at 1 bar and 120 bar have been shown in Fig. S1 of SI.
As shown in Figure S1, the correlation between data and fitted line of the form of Equation (4) is very good for all of guest molecules. The coefficients a 1 , a 2 , a 3 for all of guest species in clathrate hydrate at atmosphere pressure and 120 bar are given in Tables S2 and S3 of the SI, respectively. The thermal expansion of different guest molecules for both pressure have been calculated using the expansion coefficients given in Tables S2 and S3 and Equation (3) and the results are shown in Figure 2.
We have showed in our previous work that the TMO, EO and FA guest molecules form guest-host hydrogen bonds (see Figure S2 of SI) in contrast with non-polar guests with analogous structures, CB, CP and Et, respectively. The guest-host hydrogen bonding probability increases with temperature and so the guest-host hydrogen bonding between the TMO molecule and water is stronger than the other guest molecules. Figure 2(a) indicates that in the low temperature region (100-150 K) the effect of guest molecules size on the value of the thermal expansivity of the mixed hydrates overcomes the effect guest-host hydrogen bonding. Therefore, at high temperatures, the larger value of thermal expansion for TMO, EO and FA compared to CB, CP and Et can be attributed to guest-host hydrogen bonding with water lattice. Figure 2(a) also shows dependence of the sI hydrates thermal expansion on temperature at atmospheric pressure using the experimental data of Hester et al. [13] For comparison, the calculated isobaric thermal expansion for methane sI hydrate using MD simulation with TIP4P/ice water model at atmospheric pressure by Costandy et al. [19] are presented in Figure 2. We observe that the obtained values of thermal expansion in this work are closer to the experimental result.

Isothermal compressibility
The pressure dependence of the unit cell volume for different sI clathrates in the range of pressure 0.01-5 kbar at 150 and 250 K are shown in Figure 3. As seen from Figure 3, with an increase in the size of guest molecules, the unit cell volume of hydrate increases because the large molecule stretches lattice hydrate and has larger unit cell volume. The guest molecules at low temperature (150 K) have less effect on the unit cell volume than high temperature (250 K) because the guest-host interaction at low temperature is weak. The unit cell volume of different sI clathrate hydrate was used to determine the isothermal compressibilities (β) which it is defined as: where V represents unit cell volume, Å 3 . We suppose, β is as a function of pressure that have shown in Equation (6): By separation of variables and integration, Equation (7) is obtained: where V 0 is lattice volume at a reference pressure (P 0 ). Expansion coefficients of Equation (6) [13]. The MD data of thermal expansion for methane sI hydrate using TIP4P/ice force field for water at 1 bar by Costandy et al. [19] (green stars).
The ln (V/V 0 ) against (P-P 0 ) for variety guest species in sI clathrate hydrate has been shown at 150 and 250 K in Figure S3 of SI.
As shown in Fig. S3, the correlation fits of the data on the base Equation (7) is very good for all of the guest molecules. The expansion coefficients for TMO, EO, FA, CB, CP and Et clathrate hydrate at 150 and 250 K can be seen in Table S4 and S5 of SI, respectively. The isothermal compressibility values for different guest molecules have been calculated by expansion coefficient given in Tables S4 and S5 and Equation (6). These values for both of temperatures are shown at Figure 4. Figure 4 demonstrates that β is different for variety type of guest molecules at lower pressures and is almost similar at higher pressure. Figure 4(a) shows that the size and geometry of the guest molecules in the large cages at low pressures have a greater effect on the compressibility than at high pressures. The TMO guest molecule is tethered to the water by strong hydrogen bonds which have short average life times at high temperatures with the hydrogen bonds forming and breaking at short time intervals. Therefore, with increasing temperature the vibrational anharmonicity due to the formation and breaking hydrogen bonds increases and the compressibility for this molecule is larger than CB clathrate at 250 K.
In Figure 4(b), the values of the isothermal compressibility are compared with experimental and other MD simulation results for methane and ethane sI hydrate. Helgerud et al. [41] have reported experimental data for the bulk modulus of CH 4 sI hydrate in the temperature range of −20-15°C and pressure range of 30.5-97.7 MPa. We have extracted the isothermal compressibility values using the equation where B T is the bulk modulus. Manakov et al. [15] obtained the experimental isothermal compressibility of ethane sI hydrate between 0 and 2 MPa at room temperature. The isothermal compressibility of fully occupied methane sI hydrate was calculated using the MD simulations at 250 K with TIP4P/2005 water model by Docherty et al. [23] and at 271.15 with TIP4P/2005 water model by Ning et al. [18] and also 250 K with TIP4P/ice by Costandy et al. [19]. In Figure  4(b), we find that our isothermal compressibility results for ethane hydrate are similar to those of Docherty et al. on slope.

Heat capacity at constant volume and pressure
The heat capacity by using molecular dynamics simulation, Parallel-tempering isothermal-isobaric Monte Carlo simulations for Xe hydrate [42] [45], methane and ethane hydrates [46,47] have been calculated. Using the experimental method, the heat capacity per unit volume have been measured for structure II clathrate hydrates of 1,3-dioxolane (DO) and cyclobutanone (CB), and compared the results with similar previous measurements for tetrahydrofuran (THF).
In this work, molecular dynamics (MD) simulation was performed to evaluate the specific heat capacities for various guests in the sI clathrate hydrate. By MD simulations in NVT and NPT ensembles, the heat capacities at constant-volume (C V ) and at constant-pressure (C P ) for binary clathrate hydrate are estimated between 50-250 K and 1 bar respectively. In order to calculation C p and C V , we have obtained the fourth degree polynomial fit of E NPT and E NVT as a function of temperature (see S4 and S5 Figs. in SI), respectively and we have used from 8 and 9 equations and differentiate the fitted function to obtain heat capacity. The E NPT and E NVT are enthalpy of system and total internal energy of system, respectively.
The constant-pressure heat capacity of EO, TMO, FA, CB, CP and Et guest molecules in sI clathrate hydrate as a function of temperature at 1 bar can be fitted to Equation (10) as follows: The value of fit parameters a, b, c and d are given in Table 1.
The calculated C P of guest molecules in sI clathrate hydrate is plotted in J kg −1 K −1 in Figure 5. For different guest molecules in clathrate hydrate with increasing temperature, the constant-pressure heat capacities increase slightly. The calculated C P of EO and Et hydrate are plotted in Figure 5 compared to publish experimental results of C P for EO·6.86H 2 O sI hydrates [25], EO·6.89H 2 O sI hydrates [45] and Et·7.67H 2 O sI hydrates [47]. Although there is different between heat capacities calculated and experimental data especially at low temperatures but the trend of C P for both plots (experimental and this work) increases with increasing temperature. The cages in the gas hydrates in the lab and in nature may not be fully occupied by guest molecules and the hydrate samples may have impurity that these empty cages and impurity may also affect the experimental value of the heat capacity. The cage occupancy depends upon the temperature, pressure, guest concentration and the guest and cage types. The results of heat capacities of methane and CO 2 hydrates obtained  [41]. Experimental data for the isothermal compressibility of ethane sI hydrate at room temperature are from Manakov et al. [15] (orange circles). The MD simulation data for methane hydrate isothermal compressibility as a function of pressure at 250 K using the TIP4P/2005 water model by Docherty et al. [23] (violet stars), at 271.15 using the TIP4P/2005 water model by Ning et al. [18] (cyan multiplication sign) and 250 K using the TIP4P/ice water model by Costandy et al. [19] (yellow triangle right). from Ning et al. [18] are also higher than those of the experimental reports at low temperatures. We have obtained the following fourth-order polynomial fit of volume-constant heat capacity of EO, TMO, FA, CB, CP and Et guest molecules in sI clathrate hydrate as a function of temperature at 1 bar which the value of fit parameters are given in Table 2.
The calculated C V of guest molecules in sI clathrate hydrate is plotted in J kg −1 K −1 in Figure 6. As expected, the constantpressure heat capacity of variety guest molecules is larger than corresponding to the constant-volume heat capacity.
The results for the constant-volume heat capacity of methane sI hydrate are presented using MD simulation by English [48] and the results are also in good agreement with experiment. The obtained isochoric heat capacity using the fluctuation formula and reaction field method for 3×3×3 systems of methane hydrate at 265 K and 1 bar is equal to 2.24 kJ kg −1 K −1 . Figures 5 and 6 that corresponding to C P and C V for a variety guest molecules show that there is a relation between guest species size and specific heat capacity. According to Figures 5  and 6, at lower temperatures the largest values of C P and C V belongs to the FA and ethane clathrate hydrates where the guests have the smallest size, and the lowest value of C P and C V belongs to the guest molecules that have large size such as TMO and CB.
In our previous work, we indicated the guest-host interaction increases with temperature for TMO, EO and FA and the guest-host hydrogen bonds are weaker than host-host hydrogen bonds in the lattice. The enthalpy of hydrogen bond formation (kJ mol −1 ) has been calculated for the TMO, EO and FA guests with water. In TMO, EO and FA hydrates, the guest-water hydrogen bond formation is endothermic which implies that the hydrogen bonds the TMO-water, EOwater and FA-water hydrogen are weaker than bonds between water molecules of clathrate lattice [18]. The smaller value of the hydrogen bond formation enthalpy for TMO molecule indicates that the hydrogen bond between the TMO molecule and water is stronger than the other guest molecules. Therefore, the smaller value of heat capacity of TMO molecule compared with EO and FA can be interpreted by the smaller enthalpy of hydrogen bond formation for TMO than the other guest molecules. As seen in Figures 5 and 6, TMO heat capacity increases extremely with temperature compared with EO and FA which can be related to the probability of formation of hydrogen bonding for this guest is greater and also the numbers of hydrogen bonds that break and form are higher at high temperature.

Conclusions
In this work, we performed molecular dynamics simulation to study volumetric and thermal properties of sI hydrates including variety guest molecules. We presented the temperaturedependent of lattice parameters and thermal expansivities and also pressure-dependent volume of unit cell and isothermal compressibilities for guest molecules that are isoelectronic but have different size and polarity. For this purpose, the oxygen atom replaced with the CH 2 group such as TMO and CB, EO and CP, FA and Et which the cell vector length and unit cell volume increases. The temperature-dependent of thermal expansion of guest species in sI clathrate hydrate in the atmosphere pressure showed that guest molecules with large size and strong guest-host interactions (like TMO) have higher values of thermal expansion at high temperatures because these interactions are important at high temperature. The pressuredependent isothermal compressibilities at high temperature showed that guest-host interaction is one of the important factors that can be affected on compressing of clathrate hydrate. The stronger this interaction (like TMO molecule by hydrogen bonding with host lattice) the easier the clathrate hydrate compressed. The isothermal compressibilities and the isobaric thermal expansion were compared to reported experimental measurements and other MD simulation results for sI hydrates. The presented results exhibited the similar trends with the experimental data and other MD simulation results. The MD simulation have been performed for variety guest in cage of sI hydrates by TIP4P/ice water model to determine the effect of guest molecules size and nature on the specific constantpressure and constant-volume heat capacities at different temperature and pressure. The results at low temperatures indicated that C P and C V of the TMO, EO, FA, CB, CP and Et hydrate decreases with increasing size of guest molecules. The calculated C P (J kg −1 K −1 ) for the EO and Et molecules with experimental measurements is much less than the C P calculated in this paper which may be due to the difficulty of compressing clathrate hydrate and the presence of gases within their cavities.