#### Abstract

Structure II clathrate hydrates of pure hydrogen and binary hydrates of are studied using *ab initio* calculations to determine the stable occupancies of small cavities. *Ab initio* calculations are carried out for a double cavity consisting of one dodecahedron (small cavity) and one hexakaidecahedron (large cavity). These two cavities are attached to each other as in sII hydrates to form a double cavity. One or two molecules are placed in the small cavity and one THF (or 4 molecules) molecule is placed in the large cavity. We have determined the binding energies of the double cavities at the MP2 level using various basis sets (3-21G, 3-21G(2p), 3-21G(2p), 6-31G, 6-31G(2p), and 6-31G(2p)). Different basis sets yield different stable occupancies of the small cavity. The results from the highest basis set (6-31G(2p) with zero point energy corrections) indicate that the single occupancy is slightly more favorable than the double occupancy in both the cases of pure hydrates and THF + double hydrates.

#### 1. Introduction

Gas hydrates are crystalline molecular complexes of water and low molecular weight gases such as methane, ethane, propane, iso-butane, nitrogen, hydrogen, carbon dioxide, hydrogen sulfide, and [1]. The water molecules form cavities through hydrogen bonding [2]. The gas molecules are enclathrated in the interstitial cavities. Guest molecules and host water molecules interact through van der Waals dispersion forces. The clathrate hydrates are existing in 3 different structures: Structures I, II, and H [2]. Structures I and II are more common than structure H. One general observation is that guests with diameters of 0.4–0.55 nm form structure I and guests with diameters of 0.6–0.7 nm or less than 0.4 nm form structure II, whereas structure H is formed from mixtures of small guest molecules and larger guest molecules (0.8–0.9 nm). The occupancy of the cavities, in general, is one. But recently H_{2} clathrate hydrates have been synthesized with multiple occupancies in the cavities [3]. Noble gases are also known to form hydrates with multiple occupancies [4].

Tanaka et al. [3] synthesized pure hydrogen hydrates at 2000 bar and 250 K. Its gravimetric storage density is 5.0 wt% with two and four hydrogen molecules in the small and large cavities of sII hydrates. It can be preserved at 140 K under ambient pressure. To reduce the hydrate formation pressure, Florusse et al. [5] used tetrahydrofuran (THF) as a promoter guest molecule. They made hydrogen hydrates at much lower pressures (around 50 bar and 280 K). The hydrates formed are of sII type. H_{2} molecules enter the small cavities, whereas THF molecules enter the large cavities. H_{2} molecules also enter the large cavities. But the storage capacity is less compared to that of pure hydrogen hydrates: 2 wt% versus 5.0 wt%. Subsequently, Lee et al. [6] reported 4 wt% of the H_{2} storage density by tuning the THF composition in the liquid phase from which hydrates form. The H_{2} + THF hydrates were formed at 120 bar and 270 K.

However, in recent years, there developed a controversy regarding the occupancy of H_{2} in the small cavities of pure H_{2} hydrates and in the small cavities of the double hydrates of THF and H_{2}. Different theoretical and experimental studies yielded different occupancy numbers for the small cavity. Mao et al. [3] obtained the mole ratio H_{2} : H_{2}O to be 0.45 from volumetric measurements, which is consistent with double occupancy in small cavity and quadruple occupancy in large cavity. They also compared the integrated H_{2} Raman vibron intensities in the sII clathrates with the same amount of H_{2} in the gas phase to obtain the ratio H_{2} : H_{2}O as . Both of the H_{2} Raman vibron intensities were obtained under identical conditions. Both of these results compare well with the theoretical value of 0.47 for 4H_{2}(L) + 2H_{2}(S) case. Lee et al. [6] synthesized double hydrates of THF and H_{2} at 277.3 K and 120 bar. They gradually changed the initial THF concentration from 5.56 moL% THF to 0.15 moL% THF to obtain H_{2} composition of 4 wt% in hydrate phase at 270 K and 100 bar. From direct gas release measurements and Raman measurements, it was concluded that H_{2} molecules doubly occupy small cavities [6].

Florusse et al. [5] synthesized H_{2} + THF double hydrates at 50 bar and 279.6 K. Based on NMR measurements, they suggested that some small cavities could be doubly occupied whereas other small cavities are singly occupied. Lokshin et al. [7] determined the D_{2} clathrate structure using neutron diffraction. From their neutron diffraction study, they determined the small cage occupancy to be one, and the large cage occupancy to vary between 2 and 4. Hester et al. [8] and Strobel et al. [9] reported the single occupancy of H_{2} molecule in the small cavity using high-resolution neutron diffraction and NMR. Anderson et al. [10] obtained the same result using volumetric measurements of released H_{2} vapor from the hydrates.

Alavi et al. [11, 12] performed molecular dynamics simulations of pure H_{2} hydrates and H_{2} + THF double hydrates. They were able to confirm the single occupancy of H_{2} in pure H_{2} hydrate small cavities, but not in the small cavities of THF-H_{2} double hydrates. They observed that the double hydrates’ small cavities can be singly or doubly occupied. From *ab initio *calculations of the small cavity and the large cavity, Patchkovskii et al. [13] determined the stable occupancy in the small cavity to be 2 and in the large cavity to be 4. Sluiter et al. [14] did a density functional theory study of pure H_{2} hydrates. They found that H_{2} molecules doubly occupy the small cages.

In view of the above controversy surrounding the stable H_{2} occupancy in the small cavity, we will perform *ab initio* calculations firstly with a double cavity. We will investigate the stability of one double cavity consisting of one dodecahedron (small cavity) and one hexakaidecahedron (large cavity) as in sII hydrates at a higher level of electron basis set or correlation than the previous works [13, 14]. One or two H_{2} molecules will be placed in the small cavity and four H_{2} molecules or one THF molecule is placed in the large cavity. We will determine the binding energies of the singly occupied and the doubly occupied small cavity of the double cavity using various basis sets at the MP2 level.

#### 2. Computational Details

To determine the stable small cavity occupancy, we have combined one small cavity (5^{12}) and one large cavity (5^{12}6^{4}) to form a double cavity. The cavities share a pentagon as their common face because only 12 pentagons are available sides in the small cavity and the large cavity also has 12 pentagons with 4 hexagons. The double cavity is chosen in our ab initio calculations because we want to understand the stability of the H_{2} occupancy of the small cavity while having the two different guest molecules (4H_{2} or THF) in the large cavity. We speculate that the stability of the occupancy will be different from the previous single cavity simulation case [13].

PCGAMESS 7.0 [15] is used to carry out the *ab initio* calculations. The coordinates of O atoms of the small cavity and the large cavity are obtained from an available website [16]. The H atoms are placed in such a way that the H-bonding is maximized. The O–H bond length is 1.0 Å as used in the previous work [13] and its length then is optimized in the structure optimization of the combined double cavity. The H–O–H bond angle varies between 104.52^{0} and 109.5^{0}. We have placed 1 or 2 hydrogen molecules in the small cavity and 4 hydrogen molecules in the large cavity. H–H bond length is initially 0.74 Å and relaxed for the structure optimization of the combined double cavity. The 4 H_{2} molecules in the large cavity are assumed to form a tetrahedral cluster as in the previous study [13] and the H_{2}–H_{2} distances vary between 2.90 Å and 3.13 Å. For the two H_{2} molecules in the small cavity, the distance between H_{2} molecules is initially 2.58 Å [13] and is optimized for the structural optimization. These initial distance parameters of the tetrahedral cluster and the bihydrogen cluster were obtained by Patchkovskii et al. [13] from DFT (Density Functional Theory) calculations.

When considering THF + H_{2} binary hydrates, one THF molecule replaces the 4H_{2} cluster in the large cavity. The THF coordinates are obtained from NIST’s Chemistry Webbook [17]. DFT calculations are used to optimize the complex geometries using B3LYP functional and 3-21G basis set 1. The geometric optimizations are carried out while all atoms are allowed to move with cavity oxygen atoms’ coordinates fixed. We have then calculated the binding energies for the following four optimized configurations: 1H_{2}(S) + 4H_{2}(L), 2H_{2}(S) + 4H_{2}(L), 1H_{2}(S) + THF(L), and 2H_{2}(S) + THF(L). Here L indicates large cavity and S, small cavity. The binding energy is obtained using the following equation:

where E is the total energy and the subscripts of complex, EDC, guest (L), and guest (S) represent the complex of guest molecules and cavities, empty double cavity, guest molecules in a large cavity, and guest molecules in a small cavity, respectively.

The small cavity and large cavity are shown in Figure 1 and are generated by *MOLEKEL* software [18, 19]. The combined double cavity formed from the small cavity and the large cavity is shown in Figure 2, which is also generated by *MOLEKEL* software [18, 19]. As mentioned before, the two cavities share a pentagon as a common face. We have performed the calculations in PCGAMESS 7.0 [15] using six different basis sets: 3-21G, 3-21G(2p), 3-21++G(2p), 6-31G, 6-31G(2p), and 6-31++G(2p). The electron correlation is treated at the MP2 level. The calculations are done using a 32-bit computer node with a clock speed of 2.0 GHz and 2 GB RAM. It takes 6–8 days for each complex’s total energy calculation. Empty double cavity takes around 5 days for computing its total energy. The total energy calculations for guest molecules take few (0.5–6.0) hours on this machine.

**(a) Small cavity.**

**(b) Large cavity.**

**(a)**

**(b)**

#### 3. MP2 Results

The binding energy values without including the structure optimization are given in Table 1. The units are J/moL based on the complex double cavity. Using the initial coordinates described in the previous section, we carried out the binding energy calculations. The double occupancy in the small cavity (2H_{2}(S) + 4H_{2}(L)) is more favorable than the single occupancy (1H_{2}(S) + 4H_{2}(L)) for pure hydrogen hydrates. However, the opposite result occurs for the THF-H_{2} binary hydrates; that is, the single occupancy is slightly more stable than the double occupancy.

Performing MP2-level calculations based on the structure optimization with complex double cavities presents the results of binding energies in Table 2. These binding energies are plotted in Figures 3 and 4 for pure H_{2} and double H_{2} + THF hydrates. In Figure 3, the binding energies of the pure H_{2} hydrate are plotted for different small cage occupancies using various basis sets. It is seen that the configuration 1H_{2}(S) + 4H_{2}(L) is slightly more favorable than 2H_{2}(S) + 4H_{2}(L). The binding energies of THF + H_{2} double hydrates are shown in Figure 4. There is no definite trend in binding energy values with the level of the basis sets used and the single occupancy is as stable as the double occupancy. From the result of the highest level of basis set used (6-31++G(2p)), the configuration of single occupancy for both cases is more stable than the configuration of double occupancy, which may be consistent with the recent experimental works [7–10].

#### 4. Zero Point Energy Calculations

Molecules are not stationary at 0 K. Molecules vibrate even at 0 K. Because of this, molecules possess a vibrational energy. This energy is called the zero point vibrational energy (ZPVE) or zero point energy (ZPE). The calculations in the previous section treated the nuclei as a frozen core. By adding the ZPE to the frozen nuclei energy, we will obtain the total energy of the molecule at 0 K. ZPE calculations involve geometry optimization and frequency calculations [20]. We have obtained optimized geometries of the complexes from DFT calculations using B3LYP exchange-correlational functional and 3-21G basis set as we mentioned before. The geometric optimization is very expensive in terms of calculation time and it takes more than 3 weeks using the same PC (2.0 GHz clock speed and 2 GB RAM) used in the previous section of MP2 calculations. Thus, we carry out the geometric optimization at the high performance computing cluster available at the Texas A&M University, Kingsville, and it takes around 3 days.

The ZPEs determined by the optimized geometries are then added to the frozen nuclei energies 2. All vibrational motions of all molecules are considered in the zero-point energy calculations. From these energies, binding energies of various configurations are calculated and compared. It is valid to add the ZPE obtained from the lower level calculations to the frozen nuclei energies obtained from the higher-level calculations (MP2) [20]. These energies are called :

In finding out the stable occupancies of gas hydrates, binding energies of various configurations are compared. The results are given in Table 3. The binding energies with ZPE corrections are plotted in Figures 5 and 6. With ZPE corrections, the stable occupancy of small cavity is one for both pure H_{2} hydrates and THF + H_{2} double hydrates. This is a concurrent observation with the experimental analyses [7–10] and MD simulations [11, 12]. However, the difference of binding energies between single occupancy and double occupancy is not much significant.

Two cases in pure hydrogen hydrates in Table 3 showed positive values of binding energies. To eliminate the positive binding energies, we tried one of these two calculations with counterpoise corrections with 6-31++G(2p) basis set but the positive values are more pronounced (e.g., 55839.9 J/moL versus 8866.79 J/moL) even if the trend is still valid (single occupancy is more stable than double occupancy in the small cavity). On the other hand, we tried with free hydrogens as guest molecules instead of bi-hydrogen and tetra-hydrogen clusters. The binding energies move toward more positive (34435.08 J/moL versus 8866.79 J/moL) but the same trend is obtained that the single occupancy is more favorable than the double occupancy in the small cavity.

#### 5. Conclusions

We have tried to determine the stable small cavity occupancies of binary THF + H_{2} sII hydrate and pure H_{2} hydrates using *ab initio* calculations. The systems were studied using a double cavity. The double cavity was made up of one small cavity and one large cavity. Binding energies of pure H_{2} hydrates and binary THF + H_{2} hydrates were compared to determine the stable occupancies of the small cavity. Without ZPE corrections, the single occupancy was more preferable than the double occupancy for the small cavities of pure H_{2} hydrates. However, this was opposite without any structural optimization of the double cavity. The single occupancy was as stable as the double occupancy for THF + H_{2} binary hydrates with respect to all of basis sets regardless the structural optimization. With ZPE corrections based on the structural optimization, the single occupancy was more favorable than double occupancy. In these two cases, the single occupancy was slightly preferable in both pure H_{2} hydrates and THF + H_{2} binary hydrates with the highest level of basis set, which may be in line with recent experimental results.