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 H2 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. H2 molecules enter the small cavities, whereas THF molecules enter the large cavities. H2 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 H2 storage density by tuning the THF composition in the liquid phase from which hydrates form. The H2 + THF hydrates were formed at 120 bar and 270 K.
However, in recent years, there developed a controversy regarding the occupancy of H2 in the small cavities of pure H2 hydrates and in the small cavities of the double hydrates of THF and H2. Different theoretical and experimental studies yielded different occupancy numbers for the small cavity. Mao et al. [3] obtained the mole ratio H2 : H2O 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 H2 Raman vibron intensities in the sII clathrates with the same amount of H2 in the gas phase to obtain the ratio H2 : H2O as . Both of the H2 Raman vibron intensities were obtained under identical conditions. Both of these results compare well with the theoretical value of 0.47 for 4H2(L) + 2H2(S) case. Lee et al. [6] synthesized double hydrates of THF and H2 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 H2 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 H2 molecules doubly occupy small cavities [6].
Florusse et al. [5] synthesized H2 + 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 D2 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 H2 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 H2 vapor from the hydrates.
Alavi et al. [11, 12] performed molecular dynamics simulations of pure H2 hydrates and H2 + THF double hydrates. They were able to confirm the single occupancy of H2 in pure H2 hydrate small cavities, but not in the small cavities of THF-H2 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 H2 hydrates. They found that H2 molecules doubly occupy the small cages.
In view of the above controversy surrounding the stable H2 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 H2 molecules will be placed in the small cavity and four H2 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 (512) and one large cavity (51264) 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 H2 occupancy of the small cavity while having the two different guest molecules (4H2 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.520 and 109.50. 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 H2 molecules in the large cavity are assumed to form a tetrahedral cluster as in the previous study [13] and the H2–H2 distances vary between 2.90 Å and 3.13 Å. For the two H2 molecules in the small cavity, the distance between H2 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 + H2 binary hydrates, one THF molecule replaces the 4H2 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: 1H2(S) + 4H2(L), 2H2(S) + 4H2(L), 1H2(S) + THF(L), and 2H2(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 (2H2(S) + 4H2(L)) is more favorable than the single occupancy (1H2(S) + 4H2(L)) for pure hydrogen hydrates. However, the opposite result occurs for the THF-H2 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 H2 and double H2 + THF hydrates. In Figure 3, the binding energies of the pure H2 hydrate are plotted for different small cage occupancies using various basis sets. It is seen that the configuration 1H2(S) + 4H2(L) is slightly more favorable than 2H2(S) + 4H2(L). The binding energies of THF + H2 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 H2 hydrates and THF + H2 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 + H2 sII hydrate and pure H2 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 H2 hydrates and binary THF + H2 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 H2 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 + H2 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 H2 hydrates and THF + H2 binary hydrates with the highest level of basis set, which may be in line with recent experimental results.