Density functional theory studies of hydrogen bonding vibrations in sI gas hydrates

To analyze the vibrational modes of water and methane in structure I gas hydrates, we constructed a 178-atom supercell with two small cages of type 512 and six large cages of type 51262. We applied the density functional theory method to simulate the vibrational spectrum and normal modes of methane hydrates. In accord with our previous studies, we confirmed that two groups of hydrogen bond (H-bond) peaks (at around 291 and 210 cm−1) in the translational bands come from two kinds of intermolecular H-bond vibrational modes. This is the first investigation of H-bond vibrations in methane hydrates. The partial modes of CH4 were extracted. We found that the CH4 phonons in the translational region are below 180 cm−1 so that the influence of methane on the H-bond is insignificant. We proposed a new method to decompose gas hydrates via direct application of terahertz radiation to the H-bonds. Herein, we confirmed that CH4 molecules do not absorb this energy.


Introduction
With the depletion of oil and coal, gas hydrates will play a vital role in energy use and environmental protection. Sloan [1] estimated that the current energy reserves of methane hydrates are twice the total fuel fossil reserves, Lunine and Stevenson [2] wrote that methane hydrates make up the major satellite nebulae. And Takeya et al [3] revealed that the potential of methane hydrate used in actual methane storage. In addition, solutions to the blockage of methane hydrates in natural gas pipelines [4] and the accelerating effects of climate change [5] are highly valued. Traditional methods for extraction of natural methane hydrates include thermal stimulation, pressure reduction, and catalysts. However, commercial exploitation has not yet been realized due to economic and safety factors [6,7]. Therefore, it is imperative to develop an efficient and sustainable method to dissociate methane hydrates.
Gas hydrates are crystalline solids that are often composed of cage-like structures of host molecules surrounding various nonpolar guest molecules. Gas hydrates are divided into various structures according to the cage numbers and cage types. Among these, structure I (sI) is the most common structure found in nature, first presented by Muller and von Stackelberg in 1951 [8], in which methane constitutes the overwhelming majority of guest molecules [9]. Hoshikawa et al [10] examined hydrogen of deuterated methane hydrate by means of the neutron powder diffraction data. An sI-type crystal comprises two small cages of type 5 12 (pentagonal dodecahedrons) and six large cages of type 5 12 6 2 (a type of tetrakaidecahedral cavity), and only small molecules (0.4-0.55 nm) [11] such as methane and carbon dioxide can be encaged.
Water molecules are entrapped into cages by hydrogen bonding, whereas the host and guest molecules are connected by van der Waals forces.
Ice XVII, one kind of sI-type lattice, is derived from the phase diagram of the first H 2 -H 2 O compounds prepared in the laboratory [12,13]. Ice XVII has the same structural definition as sI: a composition of pentagons linked with other pentagons via shared vertices [14]. For sI, two additional hexagons are required to maintain the hydrogen-bond strain. However, ice XVII can only form pentagons and therefore adopts a pipe structure rather than a cage structure like sI.
Since Anthonsen [15] first presented the Raman scattering spectroscopy of halogen hydrates, many researchers have focused on the effects of encapsulated methane in the band of C-H stretching via Raman spectroscopy [15][16][17]. Sum et al [16] first assigned the peaks at 2905 and 2915 cm −1 to the large cage and small cage, respectively. The two peaks' intensity ratio of 3:1 was interpreted as the number ratio of these two cages. This phenomenon in the vibrational spectrum was also confirmed by Greathouse et al [17]. In another study, Kortus et al [18] used melanophlogite, which has the same framework structure as sI, to discuss the influence of structure on the Raman spectrum of methane.
In contrast, Schicks et al [19] explored the influence of host molecules by Raman spectroscopy rather than guest molecules and found similarities in the spectra between sI and ice Ih in O-H stretching and H-bond translational vibrations. Hanson and Berg [20] indicated that metastable methane hydrates must have the ability to exist in different cavities. Tulk et al [21] compared the host lattice of methane hydrates with ice Ih in the region of coupled O-H vibrations.
In theoretical simulations, John [22] simulated the sI structure with an ab initio method to explore the relationship between the stretching vibrational frequencies and the CH bond length. Ramya et al [23] investigated the effects of methane in large and small cages, and the properties of empty cages.
Although many experimental and theoretical studies of methane hydrates have been carried out, the molecular vibrations in the far-infrared (IR) region have yet to be investigated due to the complexities of mixed host and guest phonons. From the view of physics, to explore an efficient method of the exploitation of gas hydrates for energy purposes, one should investigate the interactive mechanisms of H-bonds inside them.
In this study, we simulated the vibrational modes of H-bonds in sI methane hydrates. We present the calculated IR and Raman spectra and phonon density of states (PDOS). Compared with ice XVII, an sI-type clathrate hydrate without guest molecules [24], the encapsulation of methane was found to have little effect on H-bonds. This finding confirmed our prediction that a new method for the dissociation of gas hydrates could be explored [25].

Computational details
We performed calculations on sI methane hydrates using the CASTEP [26] code, a first-principles density functional theory module in the Materials Studio platform. The first difficult task was to construct a hydrogen disordered sI crystal structure, a supercell with 178 atoms. We generated several zero-total-dipole-moment sI clathrate hydrate structure files by GenIce [27] script with various random seeds. We then selected the structure that possessed the most even distribution for this study. The structure is illustrated in a supplementary file (https://stacks.iop.org/NJP/22/093066/mmedia) of supporting information.
We chose the revised Perdew-Burke-Ernzerhof functional [28], a generalized gradient approximation exchange-correlation functional, for the quantum mechanics calculation. The self-consistent field and energy tolerance was set at 1.0 × 10 −9 eV/atom. The energy cutoff was set at 750 eV, and we calculated the ω(q) relationship at the gamma point. The reduced Brillouin zone of a large supercell is very small, resulting in a tiny ω(q) dispersion. The norm-conserving pseudopotential was used for calculation of phonons. The hydrostatic pressure was set at 1 MPa. Figure 1 depicts the simulated spectra of Raman, IR, PDOS (sI), and PDOS (CH 4 ), and the correlated data are listed in table 1. Theoretically, the optical vibrational modes are all Raman and IR active due to disordered hydrogen. Due to the broad scale of the intensity in various regions, we adjusted their proportions and displayed the results in four parts for comparison. In the translation region, because the polarizability does not change much, the strong group of H-bond vibrations disappears in the IR and Raman spectrum. All phonons herein can be seen from the PDOS curve. The two groups are not particularly regular due to mixing with CH 4 vibrational modes. Li et al [29] first noted in 1989 that many  ice phases have two main H-bond peaks. Li and Ross [30] suggested that ice has two different H-bond strengths; however, this model has not been widely accepted [31][32][33]. In 2017, we found two intrinsic H-bond vibrational modes in ice Ic [34]. For one H 2 O molecule in an ice Ic lattice, the strong mode includes vibration of four linked H-bonds along the bisector of the HOH angle, which is called a four-bond mode. In the weak mode, only two H-bonds vibrate while the other two remain almost stationary, which is called a two-bond mode. Later, we found that these two kinds of H-bond vibrational modes that comprise the two main peaks in the translation region are a general rule among ice phases [34][35][36][37][38][39][40], and we proved that this originated from the local tetrahedral structure of ice [25]. Regarding the hydrogen-ordered ice Ic as the ideal model, the vibrational frequency ratio of these two modes is √ 2. According to the simulated normal modes of H-bonds in clathrate ice in this work, the strongest frequency of the four-bond mode is at 291 cm −1 , and we selected a typical two-bond mode at 210 cm −1 for comparison. Figure 2 illustrates the vibration processes of these two modes. (For the dynamic processes, please see the supplementary material.) Celli et al [41] reported inelastic neutron scattering (INS) experiments of various clathrate hydrates. For their sI structure, the guest molecule was xenon, so the PDOS curve came mainly from phonons of ice. Two sharp peaks could be seen at 294 and 208 cm −1 . The revised  Perdew-Burke-Ernzerhof functional slightly underestimates the hydrogen bonding, so our simulated H-bond length is greater than that from the experiment by Klapproth [46]. Table 2 shows that the calculated average length of the H-bonds is 1.871 Å, which is 5.6% greater than the 1.772 Å length from the experimental measurement. The translational region includes a dormant Raman peak at 179 cm −1 , and only one obvious IR peak at 178 cm −1 . These results are consistent with experiments in which the strong H-bond peak was difficult to detect with a photon-scattering method [30]. Schicks et al [19] presented a Raman scattering peak at 206 cm −1 in this region. Tse et al [42] used the Perdew-Burke-Ernzerhof generalized gradient approximation functional to propose a peak at 217 cm −1 . To verify the effects of the guest molecules, methane, the partial PDOS of CH 4 is presented in figure 1. In the translation region, the main contributions of methane are below 180 cm −1 . Because no H-bonds exist between clathrate ice and methane, most of the phonons in this region come from the H-bonds of ice.

Results and discussion
Ice XVII, a kind of sI-type clathrate ice structure, was investigated in our previous study, and its vibrational spectrum is presented in figure 3 for comparison [47]. The two triangle peaks of ice XVII are obvious, similar to ice Ih [30]. However, they are less apparent in the sI structure due to mixing with guest molecules.
We compare the INS spectra of sI and liquid water in the far-IR region in figure 4 [41,48]. The spectrum of sI has two main H-bond peaks at approximately 26-37 meV [31]. However, no such H-bond vibration modes are found in liquid water because it lacks a rigid tetrahedral structure. For liquid water, the strength of molecular rotation modes is lower and leaves an only weak absorption valley in this area. We discovered this phenomenon by comparing liquid water and ice Ih [25]. Although the H-bond absorption band in sI has a slight redshift relative to ice Ih, one can see that the obvious H-bond vibrational frequencies of clathrate ice in sI gas hydrates also fall in the valley of liquid water. The application of terahertz radiation in this area to decompose gas hydrates in water surroundings may be an efficient method. This study further confirms our previous suggestion for energy resource exploitation [25].
As shown in figure 1, the spectrum of the region of libration shows three groups due to the three types of molecular rotation modes: rocking (highest peak, 598 cm −1 ), wagging (highest peak, 745 cm −1 ), and twisting (highest peak, 975 cm −1 ). The encapsulated CH 4 molecules do not contribute any phonons in this band. The simulated PDOS herein agrees well with ice XVII, as shown in figure 3. In contrast, the simulated Raman peaks only manifest one group, the main peak is found at 963 cm −1 . And the simulated IR spectrum only shows the wagging band, and the highest absorption peak is found at 854 cm −1 .
The third part in figure 1 shows the intramolecular bending modes of CH 4 and H 2 O. The two weak peaks at 1303 and 1523 cm −1 are CH 4 bending and rocking modes, which agrees well with previous reports [22, 23, 43-45, 49, 50]. There are 25 modes at around 1303 cm −1 . Figure 5(a) shows one of the vibrational modes at 1303 cm −1 , where the four hydrogens in the gold-colored molecule vibrate in the same direction,   [41] (the guest molecule is xenon) and liquid water [48]. H-bond absorption ranges from 20 to 40 meV. which can be named CH 4 bending. In contrast, the 16 modes around 1523 cm −1 are CH 4 rocking modes. Due to limited change in the dipole moment and polarizability for a CH 4 molecule, the intensity of mode at 1523 cm −1 is very weak for both Raman and IR spectra. However, Chazallon et al [44,45] observed an overtone at 3053 cm −1 and is corresponding to CH 4 rocking. Another overtone at 2570 cm −1 was assigned to bending of methane [44,50]. For the CH 4 bending modes, the changes in vibrational strength within cages of various sizes are very limited. The wavenumbers from 1619 to 1716 cm −1 are the bending modes of H 2 O. For any IR or Raman spectrum, the peaks are too small to be detected experimentally. The strengths of H 2 O bending are also remained stable in various ice phases, and figure 3 shows that the band width is nearly the same as that of ice XVII.
The PDOS curve for the CH 4 stretching band ranges from 2950 to 3125 cm −1 . The stronger peaks correspond to asymmetric stretching, and the weaker peaks correspond to symmetric stretching. Both bands show obvious splitting due to the two cage sizes. Compared with the high intensities of OH stretching, these peaks are nearly invisible in the Raman and IR spectra.  literature [16,17], the splitting of the CH 4 stretching mode is 20 cm −1 for the two sizes of cages. This observation is in accord with simulations from Atilhan et al [51] and Hiratsuka et al [49]. The splitting of the CH 4 stretching band is related to the difference of the CH bond length in the large and small cages.
The vibration band from 3144 to 3416 cm −1 represents OH stretching, as shown in the fourth part of figure 1. The IR spectrum shows two very active peaks at 3211 and 3282 cm −1 , which represent more symmetric stretching. However, the strongest Raman peak is at 3415 cm −1 , which represents asymmetric stretching modes. As shown in table 2, the OH bond length of sI is shorter than that of ice XVII so that the intramolecular OH stretching band presents a 40 cm −1 blueshift, as shown in figure 3, so the wavenumbers of the sI H-bonds have a redshift from ice XVII.

Conclusions
Based on the first-principles density functional theory method, we calculated the normal modes of sI gas hydrates, and the simulated Raman scattering and IR absorption spectra are presented for comparison. Because INS experiments can detect phonons throughout the reduced Brillouin zone, previous INS results can be compared with our PDOS curve integrated by ω(q) dispersions. With good agreement with the literature, we focused on the H-bonds in the far-IR molecular translation band. Although the two groups of H-bonds of clathrate ice do not present two distinct triangle shapes such as in ice XVII, we can still confirm that the phonons in this region come from two kinds of H-bond vibration modes. Because the contributions of CH 4 in this area are below 180 cm −1 , we confirm that the effects of methane on these two main peaks are negligible.
The interactions between methane and water cages can be ignored in all regions except the far-IR region. We observed similarities between sI and ice XVII. For the phonons of methane, the splitting of CH stretching can be observed due to the influence of differences in cage size.
We have proposed a new method to decompose gas hydrates via direct application of terahertz radiation to the H-bonds. After a comparison with the partial PDOS of CH 4 in this band, we confirmed that the CH 4 molecules do not absorb this energy. Further experimental measurements are needed.