DFT calculation of the potential energy landscape topology and Raman spectra of type I CH₄ and CO₂ hydrates

Abstract

CO₂ and CH₄ clathrate hydrates of type I were studied by means of DFT and QTAIM, in order to better understand their properties at the molecular level. Sub-cells of type I hydrates were modeled as independent rigid cages, both empty and containing guest molecules. Interaction potentials of guest molecules inside each cage, and moving from a cell to the adjacent one, were calculated using the DFT approximation B3LYP/6-311+g(d,p), considering the cases with and without long range Coulombic corrections. The selected theory level was validated by comparison of the simulated Raman spectra with the experimental ones, for the case of type I lattice at full occupation of CO₂ and CH₄, respectively. For this comparison, Fermi resonances of CO₂ were taken into account by transforming experimental bands to the corresponding theoretical non-mixed states. On the one hand, our results confirm the validity of the theory level selected for the model. We have shown the high anisotropy of the guest–cell interaction potential of the molecules analyzed, which has implications in the formulation and use of equations of state, and in the study of transport properties as well. On the other hand, our results suggest that the concentration of guest species inside type I hydrates could be computed from the comparison of experimental and predicted Raman spectra, although there are non-trivial experimental limitations to get over for that purpose.

---

1 Introduction

Clathrate hydrates (or gas hydrates) are non-stochiometric inclusion solids where water molecules form a crystalline regular network through hydrogen bonding, leaving cage structures that may enclathrate small guest molecules such as CO₂ or CH₄, arranged into periodic lattices, at low temperatures or moderately high pressures. Gas hydrates have attracted large research efforts from different scientific and technological communities due to the number of processes and applications where they are involved,¹ as for instance their potential to capture and retain molecules involved in the atmospheric greenhouse effect,^2,3 with special emphasis on potential CO₂ sequestration.^4,5 In addition, the well stated abundance of CH₄ hydrates inside certain strata of the oceanic seabed, or in permafrost substrates, has been regarded as a very promising source of natural gas,^6,7 whose exploitation is not well developed yet.

Although many of the geological and structural characteristics of clathrate hydrates are well known, other relevant microscopic aspects such as the mechanism of cage formation, or the dynamics of the lattice, are not fully understood yet. On the other hand, hydrate phase diagrams can be experimentally obtained but, again, there is a lack of a complete fundamental explanation from a microscopic point of view and, subsequently, obtaining accurate estimations through modelling remains a scientific challenge. Traditionally, the calculation of phase diagrams is done using equations of state (EoS). From this point of view, the van der Waals–Platteeuw^8,9 (vdWP) method has provided a extremely useful tool to compute the thermodynamic phase equilibrium conditions for hydrates, combined with either classical cubic EoS or more sophisticated approaches such as, for instance, the molecular based SAFT EoS.^10–14 From a different perspective, molecular simulation (MS) techniques such as Monte Carlo (MC) and molecular dynamics (MD) have also been widely used to study hydrates, as reviewed recently by Barnes and Sum.¹⁵ Using these approaches, the details of the molecular structure and intermolecular interaction potential are key aspects determining the predictive ability of the models. vdWP approaches often consider the guest molecule–hydrate interaction to be described using spherically symmetric Kihara potentials¹⁶ or even more simplified square well (SW) potentials.¹³ Nevertheless, the limitations of the Kihara potential used for this purpose have been already pointed out,¹⁷ as well as the approximation of spherically symmetric potentials, denying the angular dependence.¹⁸ When performing molecular simulations, as usual only effective pair potentials are considered, which are sometimes limited to describe the complexity of the interactions involved. These considerations provide good results as a first approximation, but further refinements might be implemented if some microscopic considerations about the structure and interactions need to be incorporated into the theories.

Related with this, there are comparatively less studies in the literature concerning first principles and ab initio clathrate calculations, and the information about the detailed atomic structure, interaction energies and different structure dependent vibrational modes is very relevant in the discussion of the refinement of thermodynamic models and adequate understanding of the dynamics of the hydrate structures. As representative examples of this, Hiratsuka et al.¹⁹ studied the vibrational spectra of CH₄ molecules in type I hydrates using ab initio MD, emphasizing the differences observed for the two different types of cages appearing in this structure. Using a similar approach, Wang et al.²⁰ used a combined strategy performing MD and electronic density functional theory (DFT) calculations to study hydrogen type II hydrates, discussing in detail both atomic structure and computed Raman spectra displacements. Alavi et al.²¹ also studied hydrogen type II hydrates, discussing the guest migration between adjacent cages, considering the different face geometries. Recently, Martos-Villa et al.²² analyzed the properties of CH₄ and CO₂ hydrates using a combination of ab initio and MD calculations. The authors used a solid periodic network approximation in their DFT calculations, and discussed the geometric parameters of the hydrates, including radial distribution functions (RDFs) and the determined shifts in the main vibrational modes. Velaga and Anderson²³ also used quantum mechanics (QM) ab initio Hartree Fock calculations to determine the CO₂–H₂O intermolecular potential within a hydrate structure. With this parametrization, they computed the thermodynamic coexistence conditions for this system and the cage occupancy. These examples show how the ab initio computational methods are currently being used to determine molecular scale parameters such as interaction force fields, structural variables, and vibrational and rotational spectra, which are extremely useful to understand the basic aspects about hydrate formation, stability and thermodynamics. This type of calculations also provide parameterizations not solely based on regression of experimental data, which are later incorporated into thermodynamic models, EoS or MS, providing a description of hydrates from different alternative and mutually enriching perspectives.

In this context, the objective of this work is to apply QM electronic DFT calculations to type I hydrates of CH₄ and CO₂, with two different objectives. First, the electronic density landscapes have been computed in order to determine the guest–host interaction potential energy inside each one of the two different types of cells and the probability of the host molecule to move from one cell to another. The objective is to evaluate the differences depending on the guest molecule and the anisotropy resulting from the cell internal geometry. This serves to refine the existing guest–cage interacting models, and also to evaluate the possibility of a future diffusion theory for the guest molecules across the crystalline network depending on the space directions and different connectivity between adjacent cages. The second objective is the theoretical determination of the Raman spectra of CH₄ and CO₂ hydrates. It is well known that hydrate Raman spectra can be measured with great accuracy^24,25 and represent an extremely useful tool to determine hydrate structures and study their formation and dissociation processes, as well as their molecular dynamics. In addition, the study of the main vibrational frequency shift of hydrates in comparison with that of pure guest molecules is a valuable tool to identify hydrates, especially when thinking of the case of mixed hydrates, as reported recently by Qin and Kuhs.²⁶

2 Methods

The type I hydrate structure can be subdivided into cages of two kinds: a regular dodecahedron (D) and a truncated hexagonal trapezohedron (T) (Fig. 1). Another usual notation is 5¹²6² for T and 5¹² for D, alluding to the number of faces of each type. One cage of each type was modeled with the structure having the H-bond ordering lower energy.²⁷ Due to the fact that the water molecules must be arranged in this type of hydrate structure with proton disorder, the number of possible conformations of H-bond orientations in a T cage satisfying the Bernal–Fowler²⁸ rules is as large as 3 043 836, making it impossible to consider all of them in energy calculations, as described by Yoo et al.²⁹ However, the difference between the lower energy conformation and the immediately higher is only about 2.5 kJ mol⁻¹.

Fig. 1 Geometrical structure of the T cage (A) and D cage (B) of type I hydrates. Red spheres (larger) correspond to oxygen atoms and white to hydrogen. Blue spheres (smaller) represent the ring critical points at the center of the faces, and blue thin lines are the lower density paths connecting them with the center of the cage. Above: top view of corresponding cells, below: side view.

Structures of cages T and D were minimized without the guest molecule by means of electronic DFT.³⁰ In particular, the B3LYP/6-311+g(d,p) approximation, as implemented in Gaussian 09³¹ computational chemistry software, was employed. As usual, B3LYP stands for Becke³² three parameter Lee–Yang–Parr³³ hybrid functional, and 6-311+g(d,p) is the chosen Pople-type³¹ basis set. In addition, with the aim to consider specifically the guest transition between cages, a dual cage TT was also modeled at the same level of theory.

Topological features of electronic density inside the T and D cages were investigated using the quantum theory of atoms in molecules³⁴ (QTAIM) approach, applying the recently released open source software Multiwfn³⁵ to wavefunctions previously calculated using Gaussian 09. In particular, bond paths and minimum density paths connecting the center of the cage with the ring critical point of each face were calculated, after the optimized wave-function described above was obtained. Minimum electronic density paths were confirmed to be orthogonal to each face in their surroundings, passing through the geometric center of it, so for the calculation of the most favorable energy barrier to enter (or exit from) the cage, we have followed these paths.

The interaction potential of the guest molecule passing through every cage face was calculated using DFT at the same level of approximation used for structure minimization. First, all the system was minimized with the guest molecule at the geometrical center, and then, the energy was recalculated at a series of guest positions following the exiting path from the host structure determined in the previous step, with the aim to obtain the corresponding energy profile. In addition, and for the sake of comparison, this calculation was repeated using the CAM-B3LYP³⁶ (Coulomb attenuating method – B3LYP) functional for a better description of the long range interactions, which were expected to be relevant in this type of systems.

We have studied the potential profile for the guest molecule considering the cases of CH₄ and CO₂ molecules. CO₂ being a linear molecule that belongs to the point group D_∞h, the potential of interaction will be different depending on the direction taken by the molecule with respect to the hexagonal or pentagonal face of the polyhedron. Although carbon dioxide is electrically neutral, the partial charge distribution between the atoms, leading to a large quadrupole moment, determines to a great extent the molecule interaction, especially when condensed phases are considered, as shown in a recent work.³⁷ To simplify the study we have considered two complementary geometries. In the first case, the molecule crosses the face with its longitudinal axis being parallel to the plane defined by the oxygen atoms in the hexagonal face. In the second case, the molecule axis was considered to be perpendicular to this plane. This last configuration minimizes the interaction energy, and thus constitutes the preferential way of crossing.

The next objective of this work is the calculation of the Raman spectra for both CH₄ and CO₂ type I hydrates. In order to make a reliable computational spectroscopic study using ab initio methods, the first step is to choose an adequate level of theory. We decided to scan through the most common basis sets which include various combinations of diffuse and polarization functions, and to consider several of the most common calculation methods: Hartree–Fock (based on Eckhart theorem, which does not include electron correlation) and post-Hartree–Fock (Moller–Plesset method 2, MP2, based on the perturbation theorem), QCISD (configuration interaction with truncation at the quadratic term) and a method based on the DFT using the B3LYP hybrid functional. With the aim of comparing the characteristic molecular variables with those retrieved from the NIST database, single point and geometry optimization calculations were performed. A comparison of the average percent deviations (APD) associated with geometrical properties, bond strength and main vibrations provided within those methods was done in order to estimate the corresponding Raman spectra. Taking into consideration the computational cost and the detail of theory needed to do further calculations, the method with the three-parameter B3LYP hybrid functional was selected again as the optimal procedure, since it has been shown to predict successfully a wide range of properties of both systems: guest molecules and the empty hydrate cell. For instance, Pele et al.³⁸ have emphasized the superiority of the B3LYP method for the calculation of Raman spectra when compared with other available alternatives. Also, a comparative analysis of harmonic vibrational calculations made using different popular functionals was done by Cappelli and Biczysko,³⁹ revealing that B3LYP results are the most reliable among them, although it being one of the less computationally demanding functionals. The basis set used here was 6-31+g(d,p) since results showed that this degree of detail provides a good estimation of bond distances and angles, and also the main features of the vibrational and rotational molecular spectra. Raman spectroscopy is successfully applied to the analysis of a wide range of systems. Additionally, this kind of spectroscopy provides useful information about molecular vibrations which can be used for qualitative and quantitative applications through spectral interpretation. The visualization of the molecular structures was obtained using VMD,⁴⁰ and the presented figures were rendered with the help of the included Tachyon⁴¹ program. Analysis of molecular structures and spectroscopic data was performed using Gabedit.⁴²

3 Results

For the investigation of the potential affecting the guest molecule inside the hydrate, the electronic density corresponding to individual cages of type T and D was calculated, and a QTAIM analysis was performed afterwards. In Fig. 1, paths connecting the center of the cage with the ring critical points, and the corresponding critical points, are depicted. These are usually called ring paths (RP). A section of the electronic density of the T cage, considering the projection in Fig. 1A (bottom), is illustrated in Fig. 2, including isocontour curves, gradient field lines, critical points and ring paths also. From a probabilistic point of view, the lower the density at a given distance, the higher the probability of finding a guest molecule at that distance, so the ring paths are the ways which would follow the guest moving in or out of the cell. They are the lines central to all the water basins in a face, which implies that they have a gradient flux equal to zero.

Fig. 2 Topological analysis of the electronic density of cage T in the YZ plane (X = 0), which passes through the center of both hexagonal faces and the center of two of the pentagonal faces also. The projection is the same as that of Fig. 1A (below). Tick labels correspond to distances in Å. Ring paths (same as in Fig. 1) are shown in black, connecting the cell center (green) with the ring critical points (orange) of the hexagonal faces (top and bottom) and two of the pentagonal faces (sides). The blue stroke represents a density value of 0.001, which is generally regarded as the molecular surface. Contours (black lines) correspond to values 0.002, 0.004, 0.008, 0.02, 0.04, 0.08, and so on, starting from the surface. Gradient lines (grey) are also shown, to visualize the attraction basins.

The type I hydrate, as previously mentioned, is composed of only two cage types: T and D. There are also two types of faces only, pentagonal and hexagonal. T cells share pentagonal faces with neighbour T cells and with D cells also, but D cells share pentagonal faces with T cells only. Actually, the I lattice can be regarded as an arrangement of T cells, with D holes in the interstices. We will, therefore, center the following discussion on T cells.

QTAIM analysis provides valuable information about the structure. Fig. 1 and 2 illustrate the general nonlinearity of pentagonal RPs, in both T and D cells, in contrast to hexagonal RPs which are in fact linear. A closer inspection shows that the pentagonal RP is always orthogonal to the ring plane in the vicinity of the ring critical point, and the curvature of the RP is simply a consequence of the combination of a radial line starting from the center and the corresponding orthogonal line through the center of the ring, so it might be explained using almost exclusively geometrical considerations.

Using this fact, we have calculated the interaction potential along the straight portion of the path around the ring of pentagonal RP and hexagonal RP in a T cell for a guest CH₄ molecule (Fig. 3), using the B3LYP functional as detailed in Section 2. Calculations were repeated using CAM-B3LYP for a better description of long-range interactions, but the energy difference obtained when including the CAM method was found to be not significant for the present study (Fig. 4). For example, dissociation energy of C–C bonds is typically referred as 4.4 eV, so quantitative values in Fig. 3 are unrealistically high. This is due to the rigidity of the cell in the calculation, while in reality the H-bond lattice is expected to have a considerable flexibility. On the one hand, structure flexibility depends strongly on the number of layers considered, making the realistic full calculation taking flexibility into account inviable by means of ab initio or DFT techniques. On the other hand, results are essentially correct for topological analysis, face comparison purposes and for describing the potential perceived for the guest molecule inside the cages. Care must be taken, however, of not considering these energy values as experimentally realistic due to the assumptions described.

Fig. 3 Interaction potential, in eV, of CH₄ approaching the center of pentagonal (hollow squares) and hexagonal (filled circles) faces from inside the T cage. Increasing negative values, in Å, correspond to distances inside the cage and 0 corresponds to the face.

Fig. 4 Interaction potential in eV of CH₄ approaching the center of the pentagonal face from inside the T cage, calculated using the B3LYP functional (hollow squares) and CAM-B3LYP (solid circles). Increasing negative values on the x axis, in Å, correspond to distances inside the cage and 0 corresponds to the face.

Taking into account the particular structure of type I, where the hexagonal faces are aligned forming straight channels, these results suggest some kind of preferential gas transport along the corresponding aligned hexagonal RPs. It is expected to see three independent and orthogonal priority ways where the movement of the molecules may take place. Guest trajectories would be in contact with pentagonal faces in each cell, but these exhibit a much lower probability of passing through them, and their contribution could be neglected in a first approximation. Moreover, the guest molecular displacement should be independent for each channel, so it might be studied using a unidimensional model of transport in each direction. From these considerations, the type I hydrate is expected to show anisotropic transport for molecules of the size of methane and smaller, in contrast with II and H types.

II and H type hydrate structures can be described considering a greater variety of cells, but in all cases, only three types of faces are present, hexagonal, pentagonal, and square, with all of them having the dimensions constrained by the H-bonds of the water molecules.³ The probability of transition through a square face will be even lower than through pentagonal one due to its reduced size. Therefore, the previous considerations can be generalized in a straightforward way to all hydrate structures.

All these considerations are only qualitative due to the intrinsic difficulty of the proper calculation at molecular level of all the mechanisms involved in flexible transitions. One possibility, proposed by Peters et al., is that guest transport is supported by a water-hopping type mechanism.⁴³ This entails breaking completely the hydrogen bonds of a water molecule in one of the faces of a cell, while the elongation of the other bonds in the ring remains unaffected during the transition. Their MD calculations show that this process will give a higher probability of transition through hexagonal faces, implying transport along the hexagonal channels, in accordance with our results. Peters arrives at a diffusion coefficient value of D ≈ 7 × 10⁻¹⁵ m² s⁻¹, which seems to be compatible with experimental observations, although it might arguably be an overestimation due to the assumptions implied in a water-hopping type mechanism. Another study in this line⁴⁴ considers enhancement of diffusivity through the presence of the so-called “help-gas”, an additional guest species that can form defects in the water networks, enhancing the main guest mobility, by means of MD and nuclear magnetic resonance. On the other hand, in a hexagonal face, the flexibility of the bonds shared in equal proportion across all of them was taken into account by Román-Pérez et al., giving a value of about 1 eV for the methane transition barrier.⁴⁵ These two approximations lay in the extremes of a series of reasonable processes: ranging from one molecule supporting the complete H-bond network distortion, to the equitable sharing of the bond elongations. If we consider also the possibility of further layers of neighbor waters contributing to the processes, then the diversity of detailed molecular mechanisms involved will rapidly grow unmanageable. But, from a macroscopic perspective of transport phenomena, all of them collapse into an average probability for each type of transition. From this point of view, the rigid approximation adopted in the present work is not only computationally affordable, but reliable, reproducible, and easily comparable with different theoretical models.

Once having analyzed the CH₄ hydrate, we have repeated the calculation for CO₂ through an hexagonal RP, considering orientations parallel and perpendicular to the hexagonal face. Results are shown in Fig. 5 in comparison with the corresponding ones of methane. It can be observed that in the case of CO₂, there is a huge difference between orientations: the less favorable one can be guessed as forbidden, whilst the more favorable height is as low as one-third of that of methane. The orientation, thus, plays a key role in the movement of CO₂ molecules through the type I lattice. When its axis is aligned with the hexagonal channels, diffusion will be very feasible, but, at the same time, inside each T cell it could be easily trapped by the dipolar interactions of its bonds. As a consequence, the type of transport could be very dependent on the occupation of D cells (even being only marginal) and additional research is necessary in this direction.

Fig. 5 Comparison of the interaction potential (in eV) of CH₄ and CO₂ approaching the center of the hexagonal face along the HRP. Increasing negative values correspond to distances, in Å, inside the cage, and 0 corresponds to the face. Triangles: CH₄, open circles: CO₂ with the axis oriented perpendicularly to the plane of the ring, open squares: CO₂ parallel to the ring.

Many thermodynamic models and EoS commonly used in the simulation of interfacial and other thermodynamic properties are based on a basic model known as particle in a box. The infinite square well potential is an approach that provides results that to a great extent can be adjusted to experimental data. Nevertheless, it is a very rough potential, so given these results, we can observe that the classical square well may be improved by smoothing the walls. If the bottom of the calculated potential can be considered flat, as in the square well, there is a finite slope toward the walls, and also to be considered is that there is a finite maximum in significant regions along the surface of the cage.

As mentioned in the introduction, Velaga and Anderson²³ used QM ab initio Hartree Fock calculations to determine the CO₂–H₂O intermolecular potential over a grid of points. These data were used in order to fit the energy results calculated over the grid to a intermolecular potential. Martos-Villa et al.²² analyzed the properties of CH₄ and CO₂ hydrates using a combination of ab initio and MD calculations. The present work follows a different strategy, as the calculations of the potential between a methane and a carbon dioxide molecule and the hydrate were done using the QTAIM results (minimum energy paths). Our approach can be regarded as complementary to the one used in Martos-Villa paper. The authors used a solid state lattice approach (plane wavefunctions) in their calculations, and we will use a molecular approach (localized wavefunctions). Some geometric parameters of the hydrates obtained using both approaches are compared in Tables 1 and 2.

Table 1 Main average geometrical parameters of cage D, comparing experimental results taken from the literature,⁴⁶ab initio calculations in this work (DFT6-31+g(d,p) theory level) and Martos-Villa et al.²² work using the CVFFH force field and GGA/PBE (optimized using SIESTA). Distances are tabulated in Å and angles in degrees^a

Features	Experimental^46	This work	HF (HF/6-311G**)^22	CVFFH^22	GGA/PBE^22
a d(C–Hm): bond distance between carbon and methane hydrogens, d(C⋯O): distance between methane carbon and oxygen atom from water, d(Hw⋯O): distance between water hydrogens and water oxygens, d(Hm⋯O): distance between methane hydrogens and water oxygens, a(H–C–H): angle in the methane molecule.
d(C–Hm)	1.070	1.093	0.955	1.089	1.072
d(C⋯O)	3.860	3.878	4.011	4.006	3.855
d(Hw⋯O)	1.772	1.783	1.952	1.709	1.710
d(Hm⋯O)	2.922	2.938	2.867	2.985	2.789
a(H–C–H)	109.5	109.5	109.5	109.5	109.5

---

Table 2 Comparative table of the main geometrical features for cage T^a

Features	Experimental^46	This work	HF (HF/6-311G**)^22	CVFFH^22	GGA/PBE^22
a d(C–Hm): bond distance between carbon and methane hydrogens, d(C⋯O): distance between methane carbon and oxygen atom from water, d(Hw⋯O): distance between water hydrogens and water oxygens, d(Hm⋯O): distance between methane hydrogens and water oxygens, a(H–C–H): angle in the methane molecule.
d(C–Hm)	1.098	1.090	1.084	1.090	1.085
d(C⋯O)	4.505	4.293	4.100	4.485	4.545
d(Hw⋯O)	1.772	1.796	2.050	1.765	1.770
d(Hm⋯O)	3.300	3.265	3.095	3.675	3.675
a(H–C–H)	109.5	109.5	109.7	109.5	109.5

---

The next stage in this study is the estimation of the corresponding Raman spectra for both CH₄ and CO₂ hydrates. Having all stationary points characterized by means of harmonic analysis performed immediately after the final optimization step of the iterative equation-solving process, normal modes that correspond to positive frequencies were plotted, fixing a temperature value of 298.15 K. By doing that, we were able to compare the Raman spectra of the isolated cages T and D, and then with both of them filled with guest CH₄ and CO₂ molecules. In each case only single occupancy was considered. The change in spectra is obvious depending on whether the structure is full or empty. There are also some significant variations depending on the guest molecule which would allow one to identify experimentally the type of guest that fills the structure in the case of an unknown gas hydrate structure. In the calculated Raman spectra of an empty D cell, a very noisy spectrum is obtained in the range of frequency shifts between 500 cm⁻¹ and 2000 cm⁻¹, which is probably due to the free movement of external hydrogens. The situation is not the same in the case of T cell, because no detectable signals appear up to 2600 cm⁻¹. When the hydrate is filled with CH₄ in the case of D cell, the 6-31G+(d,p) basis set with the B3LYP method provides a very clean spectrum in the region below 2750 cm⁻¹. The only visible and noteworthy peak in the intensity signal is a small double peak of the CH₄ molecule around 1500 cm⁻¹. The other characteristic signal appears at a wavenumber slightly above 3000 cm⁻¹. In this case, and for identification purposes, it is preferable to use the peak located around 2900–3000 cm⁻¹, due to its higher spectral intensity. In a recent paper, Komatsu et al.⁴⁷ have overlapped and compared the bands of gas CH₄ with those of the hydrate, and also those of CH₄ in aqueous solution. For the case of CH₄ in solution, we have used the polarizable continuum model (PCM) of Tomasi et al. and an explicit model of solvation introducing a grid of water molecules by hand to compute the Raman spectra. PCM is a commonly used method to model solvation effects, due to the much lower computational cost when compared with explicit solvent simulations. Modeling H₂O as a continuum made the computation feasible in a much shorter time. The molecular free energy in solution is computed using this method as the sum over the dispersion–repulsion components, electrostatic components and the cavitation energy. The cavity is defined through interlocking van deer Waals-spheres centered at atomic positions. The effect of polarization of the solvent continuum is computed by numerical integration rather than an approximation to the analytical form used in the Onsager model. According to Tulk et al.²⁵ the vibrational frequencies of the symmetric stretching C−H vibration of methane must be different on cages T and D on type I hydrates. With the objective to provide additional experimental information to this affirmation, Qin and Kuhs²⁶ determined the spectra of deuterated and normal hydrate filled with CH₄ as guest using synchrotron X-ray powder diffraction and Raman spectroscopy. The most important result was that, although there were no significant differences on relative cross sections between deuterated and hydrogenated hydrate, the Raman scattering cross sections of CH₄ in cages T and D were close but not the same. Also, Kortus et al.⁴⁸ reported this split in two bands, giving an experimental difference among them of 11 cm⁻¹, while their theoretical estimations were of only 3 cm⁻¹. Our calculations (Table 4) show that this different molecular environment causes slightly different vibrational frequencies on cages T and D, and the obtained distance among them matches closely the declared experimental shifts.

In order to understand the vibrational spectra of CH₄ in gas hydrates, it is necessary to analyze the molecular interactions between guest and host molecules and the resulting vibrational modes. These Raman signals have been studied by several different approaches such as classical MD,^49–51ab initio MD⁵² or the semi-empirical “loose cage–tight cage” model.⁵³ The last treats the anharmonic term of potential energy as a small perturbation of the harmonic potential of the free molecule using a matrix of interactions to get the vibrational frequency shifts.⁵⁴ One of the first references on this topic is that of Hiratsuka et al.¹⁹ where they used Car–Parrinello molecular dynamics (CPMD)⁵⁵ running a simulation of the order of 30 ps in the microcanonical ensemble, analyzing five different autocorrelation functions in order to make the Fourier transform to get the vibrational spectra. CH₄ has got symmetric and asymmetric stretching modes and also rocking and bending vibrational modes. For identification and structure analysis, only the stretching modes can be used, because vibrational frequencies of bending and rocking modes overlap in both cages. Some calculations such as those by Hiratsuka et al.¹⁹ show that the area of the asymmetric stretching mode is larger than that of the symmetric mode but overlap in the spectral region with the O–H stretching vibration of water. This fact actually makes it impossible to measure it experimentally nowadays. The results show that the stretching vibrational frequencies for CH₄ in T cages are lower than those in D cages. According to Hiratsuka et al.,¹⁹ the C–H length (stretching) and the C–H–C angle (rocking) of the vibration spectra provide the clearest intramolecular modes distinguishing the occupancy of CH₄ in T and D cages of I hydrates. It is now well known that the symmetric C–H stretching vibrational frequencies of CH₄ in the T cage are lower than those in D cages.^19,22,52 This is also true for this work. The microscopic interpretation for these results arises from the fact that the cages are not perfectly symmetric, but there is an elongation because water molecules in hexagonal faces are closer to the centre than the ones of the elongated direction. This asymmetry makes the CH₄ molecule not to be located in the center and creates an anisotropic interaction of the guest molecule with water molecules in the network. The methane C–H bond equilibrium distance in D cells is slightly smaller than that in T cells, in which it is longer due to elongation of T cells. This work has obtained values of 1.070 Å for the C–H bond in the small and 1.098 Å for the bigger cell. These distances indicate different degrees of interaction with the water molecules of the cell and the anisotropy of the elongated cell interaction as discussed above. All these results are compared in Table 3. Although all methods have their advantages and limitations, such as the treatment of the fictitious electron mass in the case of Car–Parrinello dynamics or the van der Waals correction in the DFT method, we can observe some general results: the trend of vibrational shifts is similar regardless of the method employed, even with the “loose cage–tight cage” semi-empirical method.

Table 3 Comparison of stretching frequencies (cm⁻¹) and C–H bond lengths (Å) as calculated by different methods

	Cage type	This work	Martos-Villa^22 HF 6-31G**	Martos-Villa SIESTA GGA/PBE	Hiratsuka^19	Tse^52	Experimental
Symmetric stretch	D	2917	2920	2920	2871	3014	2915.2914
	T	2908	2901	—	2863	2967	2905.2901
Asymmetric stretch	D	3056	3021	3066	2976	3167	—
	T	3021	3008	—	2960	3152	—
C–H bond length	D	1.070	—	1.072	1.097	1.107	1.098
	T	1.098	1.083	1.085	1.098	1.112	1.148

---

Table 4 Comparison between vibrational symmetric stretching frequency values of CH₄ as a guest molecule in the hydrate, CH₄ in aqueous solution (PCM model) and gas CH₄, all in cm⁻¹, APD in percentage

	Cell D	Cell T	Methane in water	Methane (gas)
Original band	3030	3021	3021.78; 3029.24	3030.68
Scaled band	2917	2908	2911.78; 2919.24	2920
Literature	2915	2905	2910; 2917	2918.48
APD	0.069	0.103	0.061; 0.061	0.052

---

From a quantitative point of view, ab initio vibrational frequencies are quite often larger than the ones observed experimentally. Errors appear mainly because of the incomplete incorporation of electron correlation, neglect of anharmonicity effects and use of finite basis sets. A good agreement between experimental and theoretical frequencies can be obtained using a convenient scaling factor, which is usually calculated by comparing frequency calculations of a large amount of molecules in gas phase with experimental results. In this case, we found that calculated frequencies overestimate experimental values. We have considered that, since the CH₄ hydrate is a crystalline inclusion compound, it would not be adequate to apply a scaling factor obtained from gas phase. A better solution is to calculate a hydrate-specific one. In order to do so, we calculated the ratios between various computationally obtained frequencies of vibration of CH₄ hydrates and the corresponding experimental values. The scale factor would then take a value of 0.963 ± 0.001, considering the limitations in the accuracy of its calculation, as pointed out by Teixeira et al.⁵⁶ Taking this scaling factor into account, Table 4 shows as well the Raman shifts for the rescaled spectrum. The obtained APD errors show that the agreement between experimental bands and ab initio calculations is accurate using the appropriate scaling factor, and the calculation method used to estimate the Raman spectra presents the same relative distances between the characteristic bands in filled T and D cells, aqueous, and gas CH₄, respectively. This result emphasizes the reliability of the method used and the level of calculation selected in this case.

Concerning the Raman spectra of type I CO₂ hydrate, Nakano et al.⁵⁷ reported two experimental characteristic bands placed at 1280 and 1386 cm⁻¹. This characteristic doublet, denoted as Fermi dyad,⁵⁸ was also reported by Qin and Kuhs.²⁶ These two bands correspond to the lower (ν_c−) and upper (ν_c+) levels of Fermi resonance corresponding to vibrational modes with very similar energies, the symmetric stretching vibration (ν₁) and the first overtone of the bending modes (2ν₂). This quantum resonance effect, which has been described in detail for the particular case of CO₂ by McCluskey and Stoker,⁵⁹ produces a shift of both original frequencies from their original values and an increase in the intensity of the weaker vibrational mode band. The Raman spectra computed in this work, as shown in Fig. 7 and 8, do not show the (ν_c+) nor the (ν_c−) bands, and the reason for this is that their description needs the use of anharmonic vibration models in the calculation of the Raman spectrum, but anharmonic analysis is currently available for only IR frequencies in Gaussian but not for Raman frequencies. The absence of these Fermi dyad bands deserves a more detailed analysis, as it could lead to the conclusion that the level of theory selected is not describing properly the vibrational states of the molecule. Actually, the estimated spectra show clearly all the Raman shifts corresponding to the expected vibrational modes, which are the symmetric (ν₁) and asymmetric (ν₃) stretching, and two bending (ν_2a and ν_2b) modes. Table 6 shows the obtained shifts (not rescaled in this case) for each of these modes, for isolated CO₂, and also for the molecule inside T and D cells. Both bending modes, which are Raman inactive in the isolated molecule but do appear in the hydrate case due to the coupled interaction with the host molecules, are very close to each other, so they merge into a single band from an experimental perspective, and it is the first overtone of these combined bending frequencies along with the symmetric stretching that produces the Fermi resonance, increasing their original low intensity, an effect which cannot be captured with conventional DFT, as noted before.

In order to compare the CO₂ spectral bands obtained experimentally with the DFT counterparts, we have taken an alternative approach. Instead of implementing an ab initio or DFT anharmonic approximation, the corresponding pure harmonic excitations from the experimental bands were computed. For that aim, the assumption that the real system can be described using a perturbative treatment must be made, as the result of the interaction of the harmonic unperturbed states. This is a standard analytical procedure with negligible computational cost and rather accurate results, even when compared with other computationally demanding alternatives such as iterative numerical procedures (see e.g. the references by Andersson et al.⁶⁰ and Pinto et al.⁶¹). With the approach followed, an interaction Hamiltonian (Ĥ′) accounting for the difference between harmonic and anharmonic models is then assumed, therefore assigning the unperturbed Hamiltonian (Ĥ⁰) to the harmonic oscillator and the perturbed (Ĥ = Ĥ⁰ + Ĥ′) to the anharmonic one. In the following, spectroscopic units (cm⁻¹) will be considered as energy units, and states will be labeled with the corresponding energy. With the fundamental state energy E₀₀₀ being equal to 0, spectroscopic signals can be directly assigned to excitation energies, denoting as E₊ the excited state of transition ν_c+, E₋ the excited state of ν_c−, E₁₀₀ identified then as ν₁ and E₀₂₀ as 2ν₂. Therefore, the anharmonic states E₊ and E₋ will be the result of the mixing interaction of the unperturbed states E₁₀₀ and E₀₂₀. Under the described scheme, the following secular determinant of the characteristic equation is obtained:

where E_F stands for the Fermi dyad energy and E_int stands for the mixing interaction energy. Some calculations over the determinant show that the solutions of E_F, E₋ and E₊ must satisfy the following two conditions:

E₋ + E₊ = E₁₀₀ + E₀₂₀

(E₋ − E₊)² = (E₁₀₀ − E₀₂₀)² + 4E_int², (2)

or, equivalently, the spectral bands have to hold:

ν_c− + ν_c+ = ν₁ + 2ν₂

(ν_c− − ν_c+)² = (ν₁ − 2ν₂)² + 4E_int². (3)

These relations allow one to calculate the Fermi dyad if harmonic levels are known, and the other way round, it is also possible to calculate the pure harmonic levels corresponding to the Fermi states, along with the mixing (interaction) energy E_int. We are interested in the last option, in order to obtain experimental harmonic levels to compare with. To do so, we have first considered experimental values available in the literature for CO₂ in gas state: from the IR band⁶²ν₂ = 667 cm⁻¹ and Raman levels⁶³ν_c− = 1284 cm⁻¹ and ν_c+ = 1381 cm⁻¹, and this way we have estimated ν₁ = 1338 cm⁻¹ and E_int = −51.998 cm⁻¹. Considering another different source of experimental values for Raman levels,²⁶ν_c− = 1280 cm⁻¹ and ν_c+ = 1386 cm⁻¹, and the same IR band used before, we have obtained alternative values of ν₁ = 1332 cm⁻¹ and E_int = −52.998 cm⁻¹. The ν₁ values are in remarkably good agreement with the standard experimental reference⁶² of ν₁ = 1333 cm⁻¹. Both estimated values of E_int are also in very good agreement with E_int = −51.232 cm⁻¹ reported by McCluskey and Stoker.⁵⁹

This last value, the energy interaction, was supposed to be constant due to the fact that the experimental Fermi bands are placed 103 cm⁻¹ apart for CO₂ in gas phase and 104 cm⁻¹ in type I hydrate.⁶³ Hence, we used the calculated interaction energy from CO₂ (gas), together with experimental hydrate values of ν_c− and ν_c+, to estimate ν₁ and ν₂ inside the hydrate. These results are shown in Table 5, in comparison with the corresponding DFT calculations of the present work previously discussed. The agreement between both sets of values is satisfactory.

Table 5 Experimental Raman signals (expressed in cm⁻¹) for the Fermi levels (ν_c−, ν_c+) and the corresponding harmonic bands (ν₁, ν₂) of CO₂ calculated in this work by means of DFT. The superscript (^a) denotes harmonic modes calculated from experimental Fermi bands using eqn (2). The bending band actually corresponds to two orthogonal modes, which can be resolved by DFT, and whose values are shown under ν₂ in those cases

	ν c−	ν c+	ν 1	ν 2
CO2 gas
DFT			1347	679; 679
Chazallon^63	1284	1381
Qin and Kuhs^26	1280	1386
Person and Zerbi^62			1333	667
CO2 hydrate (sI)
DFT, D cell			1374	644; 646
DFT, T cell			1376	659; 660
Chazallon^63	1277	1381	1329^a	664^a
Qin and Kuhs^26	1275	1379	1327^a	664^a

---

Table 6 Raman shifts of CO₂ vibrational modes as an isolated molecule and also as guest in T and D cells of type I hydrates

Type	Shift (cm^−1)	Intensity
Isolated CO2
ν 1	1347.09	20.2106
ν 3	2417.27	0.0000
ν 2a	678.76	0.0000
ν 2b	678.76	0.0000
T cell
ν 1	1376.14	17.6293
ν 3	2420.66	0.0329
ν 2a	659.26	0.0981
ν 2b	659.76	0.0301
D cell
ν 1	1374.17	19.0925
ν 3	2424.56	0.0258
ν 2a	643.77	0.0981
ν 2b	646.27	0.0301

---

Taking into account that no scaling factor was applied to DFT spectra, we can observe that there is a good agreement in absolute values between our computational results and the experimental estimations under the harmonic model. Nevertheless, there is a discrepancy in the tendency when comparing the obtained ν₁ DFT values in the gas and hydrate, as the values are higher in the hydrate, contrary to the trend of the experimental Fermi dyad frequencies, which are always slightly smaller in the hydrate case. This discrepancy may be attributed to the fact that this approach uses a number of simplifications, such as the fact of discarding anharmonic effects or considering a limited number of molecules in the hydrate. In contrast, estimated ν₂ values follow the correct trend, showing a small variation in gas state when compared with the hydrate (again, considering the mean value of the two ν₂ bending modes) for both experimental references.

An additional difficulty arises from the analysis of the different type of cages. Experimental results cannot be obtained independently for each type of cage, but this is not the case of computational approaches. DFT calculations performed are able to resolve the difference in the detected ν₁ symmetric stretching band corresponding to the T and D cages, as shown in Fig. 6.

Fig. 6 Comparison of the symmetric stretch C–O. Dashed black line: D cage. Continuous red line: T cage.

Fig. 7 Calculated Raman spectra for T cell (without scaling), where the first image corresponds with the empty hydrate, the second is the T cell filled with methane and the third one has carbon dioxide as the guest molecule. The enlarged part of the figure corresponds with an interesting region to identify the guest molecule (see Results for more information). Wavenumber is in cm⁻¹, and intensity is expressed in arbitrary units.

Fig. 8 Calculated Raman spectra for D cell (without scaling) where the first image corresponds with the empty hydrate, the second one is the D cell filled with methane and the third one has carbon dioxide as the guest molecule. Wavenumber is in cm⁻¹, and intensity is expressed in arbitrary units.

This shift between T and D cell frequencies has been a matter of debate, and for instance Sum et al.²⁴ affirmed that no split peaks could be resolved in their experimental Raman spectra. If the guest is filling both cavities, the peak should be split, so they concluded that in their experiments the CO₂ was only present in the large cages. In their work, Cao et al.⁶⁴ also focused on the same band detected here, the C–O symmetric stretching mode, and reported a difference of 4 cm⁻¹ between computed CO₂ bands in the small and large cage of the I hydrate using ab initio DFT quantum calculations with the B97-D functional. Our calculations, as reported in Table 6, yield a similar separation value, 1.97 cm⁻¹, and both bands have been represented in Fig. 6. Nevertheless, and from a practical perspective, Qin and Kuhs²⁶ report that these two peaks cannot be resolved experimentally due to the typical resolution of the spectrometers. In addition to spectrometer resolution, in our opinion other factors may also play a role here, such as the fact that experimental peaks have a certain width, which may be comparable to the device resolution, producing overlapping due to the proximity of the bands. As a result, the shift between T and D cell bands seems to be not useful as a practical tool to discuss the relative occupancy of both types of cells, although from a theoretical point of view it can be estimated with great accuracy. This is coincident with the conclusions obtained from the experimental studies of Murphy and Roberts,⁶⁵ Kumar et al.⁶⁶ and Beeskow-Strauch and Schicks.⁶⁷

A final consideration can be mentioned concerning the obtained spectral dependence of the CO₂ molecular orientation inside the cell. Udachin et al.⁶⁸ used X-ray diffraction data to determine the preferential orientations of the molecule inside the I hydrate cells. For the case of T cells, the authors describe, at low temperatures (173 K), a set of preferential orientations ranging from the molecule lying in the plane parallel to both hexagonal faces, to a tilt angle of the O–C–O axis of 14.4 degrees with this plane. The authors declare also that these preferential orientations vanish at higher temperatures due to the increased dynamics of the H₂O molecules. In our case, the calculation of the Raman spectra was performed using the optimized structure obtained with the first stage QM calculations, and imposing a temperature of 298.15 K. As the orientation of the molecule might affect the Raman displacements, a further verification was then necessary to quantify this effect under the conditions tested. With this objective, a new Raman calculation was performed, fixing the orientation tilt of the molecule to 14.4 degrees reported by Udachin and coworkers. The results show that the displacement of the spectral bands obtained is lower than 1 cm⁻¹, so the effect of this molecular orientation on the results presented here may be considered negligible.

The cited characteristic bands placed around 1500 cm⁻¹ can be outlined as distinctive features of both hydrates studied. For the CH₄ case, the mark is going to appear above this value while for CO₂ it is going to be placed below. In the case of coexistence (mixed CH₄ and CO₂ structures) the expected result is a combination band centered around 1500 cm⁻¹. Depending on the relative concentration, the band should be displaced above or below this value (CH₄ or CO₂ predominance respectively), which could serve to determine the relative abundance of both types of guest molecules in the hydrate structure.

4 Conclusions

In the present work some advances in the elucidation of the properties at the molecular level of CH₄ and CO₂ type I hydrates were made. Using DFT and QTAIM computational tools, topological features of the electronic density in both type D and T cells were analyzed. In addition, the energy profile of the guest molecule transition from one cell to the adjacent one was calculated, considering a rigid lattice. The main conclusion is that for the guest molecules under consideration the potential inside the clathrate cavities is highly asymmetric, having important implications. Namely, the square well and spherical approximations usually considered in EoS calculations are insufficient in the solid phase, with our results suggesting the convenience of smoothing the hard walls of a square well (or a Lennard-Jones type repulsion) and also considering an anisotropic interaction potential depending on the application. This conclusion can easily be generalized to other neutral molecules of similar size, and to type II and type H hydrate lattices due to the type of cells in which they can be subdivided. For transport phenomena, our results suggest that one preferential direction has to be considered: the central axis of hexagonal faces of T cells. In type I hydrates, this fact translates into three mutually orthogonal directions of unidimensional channels for the guest molecules.

With Raman spectroscopy being an invaluable and ubiquitous experimental tool for the study of the molecular properties, we have considered to test the computational model employed in the DFT study with available experimental Raman data. This strategy provides us with two convenient tools for the future: the necessary bridge between the theoretical-computational and the empirical points of view for the study of molecular structures and processes and, implicitly, the experimental validation of the chosen computational models. DFT at the B3LYP/6-311+g(d,p) level of theory was confirmed to be a convenient compromise between good experimental prediction and moderate computational resources for the systems under consideration. Spectra of T and D type cells, both empty and containing one guest molecule (CO₂ or CH₄), were calculated, showing good agreement with experimental values of clathrates at full occupation. Fermi resonance of CO₂ containing systems was tested using standard harmonic DFT calculations, by comparison of computational results with theoretical unmixed states obtained from the experimental Fermi bands. Finally, our results suggest the possibility of theoretically calculating the concentration of each species inside the hydrate, by comparing the calculated and the experimental spectra, but experimental limitations might make this hardly viable nowadays; therefore, more research in this direction will be necessary to overcome these difficulties.

Acknowledgements

The authors acknowledge project Ref. FIS2012-33621 (cofinanced with EU FEDER funds), from the Ministerio de Economía (Spain). E. Péré and M. Rérat from IPREM lab. (UPPA, France), A. Desmedt from ISM lab. (University Bordeaux 1, France) and B. Chazallon from PhLAM Lab. (University Lille 1, France) are thanked for their valuable scientific comments and suggestions.

References

1. A. K. Sum, C. A. Koh and E. D. Sloan, Ind. Eng. Chem. Res., 2009, 48, 7457–7465.
2. E. D. Sloan, Nature, 2003, 426, 353–359.
3. E. D. Sloan and C. Koh, Clathrate Hydrates of Natural Gases, CRC Press, New York, 3rd edn, 2008.
4. K. Shin, Y. Park, M. Cha, K.-P. Park, D.-G. Huh, J. Lee, S.-J. Kim and H. Lee, Energy Fuels, 2008, 22, 3160–3163.
5. G. Ersland, J. Husebo, A. Graue and B. Kvamme, Energy Procedia, 2009, 1, 3477–3484.
6. M. D. Max, A. H. Johnson and W. P. Dillon, Economic Geology of NaturalGas Hydrate, Springer, 2006.
7. Natural Gas Hydrate in Oceanic and Permafrost Environments, ed. M. D. Max, Springer, 2011.
8. J. C. Platteeuw and J. H. van der Waals, Mol. Phys., 1958, 1, 91–96.
9. J. H. van der Waals and J. C. Platteeuw, Adv. Chem. Phys., 1959, 2, 1–57.
10. X. S. Li, H. J. Wu and P. Englezos, Ind. Eng. Chem. Res., 2006, 45, 2131–2137.
11. A. Martín and C. J. Peters, J. Phys. Chem. B, 2009, 113, 7548–7557.
12. P. Paricaud, J. Phys. Chem. B, 2011, 115, 288–299.
13. S. Dufal, A. Galindo, G. Jackson and A. J. Haslam, Mol. Phys., 2012, 110, 1223–1240.
14. H. Jiang and H. Adidharma, Chem. Eng. Sci., 2012, 82, 14–21.
15. B. C. Barnes and A. K. Sum, Curr. Opin. Chem. Eng., 2013, 2, 184–190.
16. W. R. Parrish and J. M. Prausnitz, Ind. Eng. Chem. Process Des. Dev., 1972, 11, 26–35.
17. R. Sun and Z. H. Duan, Geochim. Cosmochim. Acta, 2005, 69, 4411–4424.
18. J. Klauda and S. Sandler, J. Phys. Chem. B, 2002, 106, 5722–5732.
19. M. Hiratsuka, R. Ohmura, A. K. Sum and K. Yasuoka, J. Chem. Phys., 2012, 136, 044508.
20. J. W. Wang, H. L. Lu, J. A. Ripmeester and U. Becker, J. Phys. Chem. C, 2010, 114, 21042–21050.
21. S. Alavi and J. A. Ripmeester, Angew. Chem., Int. Ed., 2007, 46, 6102–6105.
22. R. Martos-Villa, M. Francisco-Márquez, M. P. Mata and C. I. Sáinz-Díaz, J. Mol. Graphics Modell., 2013, 44, 253–265.
23. S. C. Velaga and B. J. Anderson, J. Phys. Chem. B, 2014, 118, 577–589.
24. A. K. Sum, R. C. Burruss and E. D. Sloan, J. Phys. Chem. B, 1997, 101, 7371–7377.
25. C. A. Tulk, J. A. Ripmeester and D. D. Klug, Gas Hydrates: Challenges for the Future, New York, 2000, pp. 859–872.
26. J. Qin and W. F. Kuhs, AIChE J., 2013, 59, 2155–2167.
27. S. Yoo, M. V. Kirov and S. S. Xantheas, J. Am. Chem. Soc., 2009, 131, 7564–7566.
28. J. D. Bernal and R. H. Fowler, J. Chem. Phys., 1933, 1, 515.
29. S. Yoo, M. V. Kirov and S. S. Xantheas, J. Am. Chem. Soc., 2009, 131, 7564–7566.
30. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140, A1133.
31. Gaussian 09 Revision D.01, Gaussian Inc., Wallingford CT, 2009.
32. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
33. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785.
34. R. Bader, Chem. Rev., 1991, 91, 893.
35. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
36. T. Yanai, D. Tew and N. Handy, Chem. Phys. Lett., 2004, 393, 51–57.
37. G. Pérez-Sánchez, D. González-Salgado, M. M. Piñeiro and C. Vega, J. Chem. Phys., 2013, 138, 084506.
38. L. Pele, J. Sebek, E. O. Potma and R. B. Gerber, Chem. Phys. Lett., 2011, 515, 7–12.
39. C. Cappelli and M. Biczysko, Time-independent approach to vibrational spectroscopy, in Computational Strategies for Spectroscopy, ed. V. Barone, Willey, 2012, pp. 309–360.
40. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
41. J. Stone, MSc thesis, Computer Science Department, University of Missouri-Rolla, 1998.
42. A.-R. Allouche, J. Comput. Chem., 2011, 32, 174–182.
43. B. Peters, N. E. R. Zimmermann, G. T. Beckham, J. W. Tester and B. L. Trout, J. Am. Chem. Soc., 2008, 130, 17342–17350.
44. A. Torres-Trueba, M. C. Kroon, C. J. Peters, I. L. Moudrakovski, C. I. Ratcliffe, S. Alavi and J. A. Ripmeester, J. Chem. Phys., 2014, 140, 214703.
45. G. Román-Pérez, M. Moaied, J. M. Soler and F. Yndurain, Phys. Rev. Lett., 2010, 105, 145901.
46. A. Klapproth, PhD thesis, Universität zu Göttingen, 2002.
47. H. Komatsu, M. Ota, R. L. Smith and H. Inomata, J. Taiwan Inst. Chem. Eng., 2013, 44, 517–537.
48. J. Kortus, G. Irmer, J. Monecke and M. R. Pederson, Modell. Simul. Mater. Sci. Eng., 2000, 8, 403–411.
49. H. Itoh and K. Kawamura, Ann. N. Y. Acad. Sci., 2000, 912, 693–701.
50. J. A. Greathouse, R. T. Cygan and B. A. Simmons, J. Phys. Chem. B, 2006, 110, 6428–6431.
51. F. Castillo-Borja, R. Vazquez-Roman and U. I. Bravo-Sanchez, Mol. Simul., 2010, 36, 229–239.
52. J. S. Tse, J. Supramol. Chem., 2002, 2, 429–433.
53. G. C. Pimentel and S. W. Charles, Pure Appl. Chem., 1963, 7, 111.
54. S. Subramanian and E. D. Sloan, J. Phys. Chem. B, 2002, 106, 4348–4355.
55. R. Car and M. Parrinello, Phys. Rev. Lett., 1985, 55, 2471–2474.
56. F. Teixeira, A. Melo and M. N. D. S. Cordeiro, J. Chem. Phys., 2010, 133, 114109.
57. S. Nakano, M. Moritoki and K. Ohgaki, J. Chem. Eng. Data, 1998, 43, 807–810.
58. E. Fermi, Z. Phys., 1931, 71, 250–259.
59. C. W. McCluskey and D. S. Stoker, 2006, arxiv:physics/0601182.
60. M. P. Andersson, P. Uvdal and A. D. MacKerell, J. Phys. Chem. B, 2002, 106, 5200–5211.
61. A. S. S. Pinto, R. B. de Barros, M. N. D. S. Cordeiro, J. A. N. F. Gomes, A. R. Garcá and L. M. Ilharco, Surf. Sci., 2004, 566–568, 965–970.
62. Vibrational Intensities in Infrared and Raman Spectroscopy, ed. W. B. Person and G. Zerbi, Elsevier, Amsterdam, 1982.
63. B. Chazallon, CO₂ Capture by Gas Hydrate Crystallization: Investigation of Equilibrium and Compositional Properties of CO₂-N₂ Hydrates by Micro-Raman Spectroscopy, in Physics and Chemistry of Ice, ed. Y. Furukawa, et al., Hokkaido Univ. Press, 2011, pp. 173–181.
64. X. Cao, Y. Su, Y. Liu, J. Zhao and C. Liu, J. Phys. Chem. A, 2014, 118, 215–222.
65. P. J. Murphy and S. Roberts, Geochim. Cosmochim. Acta, 1995, 59, 4809–4824.
66. R. Kumar, S. Lang, P. Englezos and J. Ripmeester, J. Phys. Chem. A, 2009, 113, 6308–6313.
67. B. Beeskow-Strauch and J. M. Schicks, Energies, 2012, 5, 420–437.
68. K. A. Udachin, C. I. Ratcliffe and J. A. Ripmeester, J. Phys. Chem. B, 2001, 105, 4200–4204.

---