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

CO2 and CH4 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 CO2 and CH4, respectively. For this comparison, Fermi resonances of CO2 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.


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 2 or CH 4 , 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, 1 as for instance their potential to capture and retain molecules involved in the atmospheric greenhouse effect, 2,3 with special emphasis on potential CO 2 sequestration. 4,5In addition, the well stated abundance of CH 4 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).1][12][13][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. 15 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 moleculehydrate interaction to be described using spherically symmetric This journal is © the Owner Societies 2015 Kihara potentials 16 or even more simplified square well (SW) potentials. 13Nevertheless, the limitations of the Kihara potential used for this purpose have been already pointed out, 17 as well as the approximation of spherically symmetric potentials, denying the angular dependence. 18When 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. 19 studied the vibrational spectra of CH 4 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. 20 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. 21also studied hydrogen type II hydrates, discussing the guest migration between adjacent cages, considering the different face geometries.Recently, Martos-Villa et al. 22 analyzed the properties of CH 4 and CO 2 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 23 also used quantum mechanics (QM) ab initio Hartree Fock calculations to determine the CO 2 -H 2 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 4 and CO 2 , with two different objectives.First, the electronic density landscapes have been computed in order to determine the guesthost 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 4 and CO 2 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. 26

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 12 6 2 for T and 5 12 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. 27Due 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 28 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. 29 However, the difference between the lower energy conformation and the immediately higher is only about 2.5 kJ mol À1 .
Structures of cages T and D were minimized without the guest molecule by means of electronic DFT. 30In particular, the B3LYP/6-311+g(d,p) approximation, as implemented in Gaussian 09 31 computational chemistry software, was employed.As usual, B3LYP stands for Becke 32 three parameter Lee-Yang-Parr 33 hybrid functional, and 6-311+g(d,p) is the chosen Pople-type 31 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 34 (QTAIM) approach, applying the recently released open source software Multiwfn 35 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 36 (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 4 and CO 2 molecules.CO 2 being a linear molecule that belongs to the point group D Nh , 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. 37To 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 4 and CO 2 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. 38 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, 39 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, 40 and the presented figures were rendered with the help of the included Tachyon 41 program.Analysis of molecular structures and spectroscopic data was performed using Gabedit. 42ines, 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.
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 4 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.
Taking into account the particular structure of type I, where the hexagonal faces are aligned forming straight channels, 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.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. 3The 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. 43This 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 E 7 Â 10 À15 m 2 s À1 , which seems to be compatible with experimental observations, although it might arguably be an overestimation due to the assumptions implied in a waterhopping type mechanism.Another study in this line 44 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 Roma ´n-Pe ´rez et al., giving a value of about 1 eV for the methane transition barrier. 45These 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 4 hydrate, we have repeated the calculation for CO 2 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 2 , 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 2 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.
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 23 used QM ab initio Hartree Fock calculations to determine the CO 2 -H 2 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. 22 analyzed the properties of CH 4 and CO 2 hydrates using a combination of ab initio and MD calculations.The present work follows a different strategy, as the calculations of the This journal is © the Owner Societies 2015 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.
The next stage in this study is the estimation of the corresponding Raman spectra for both CH 4 and CO 2 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 4 and CO 2 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 À1 and 2000 cm À1 , 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 À1 .When the hydrate is filled with CH 4 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 À1 .The only visible and noteworthy peak in the intensity signal is a small double peak of the CH 4 molecule around 1500 cm À1 .The other characteristic signal appears at a wavenumber slightly above 3000 cm À1 .In this case, and for identification purposes, it is preferable to use the peak located around 2900-3000 cm À1 , due to its higher spectral intensity.In a recent paper, Komatsu et al. 47 have overlapped and compared the bands of gas CH 4 with those of the hydrate, and also those of CH 4 in aqueous solution.For the case of CH 4 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 2 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. 25 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 26 determined the spectra of deuterated and normal hydrate filled with CH 4 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 4 in cages T and D were close but not the same.Also, Kortus et al. 48reported this split in two bands, giving an experimental difference among them of 11 cm À1 , while their theoretical estimations were of only 3 cm À1 .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 4 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][50][51] ab initio MD 52 or the semi-empirical ''loose cage-tight cage'' model. 53The 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. 54One of the first references on this topic is that of   Hiratsuka et al. 19 where they used Car-Parrinello molecular dynamics (CPMD) 55 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 4 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. 19 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 4 in T cages are lower than those in D cages.According to Hiratsuka et al., 19 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 4 in T and D cages of I hydrates.It is now well known that the symmetric C-H stretching vibrational frequencies of CH 4 in the T cage are lower than those in D cages. 19,22,52This 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 4 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.
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 4 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 4 hydrates and the corresponding experimental values.The scale factor would then take a value of 0.963 AE 0.001, considering the limitations in the accuracy of its calculation, as pointed out by Teixeira et al. 56 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 4 , 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 2 hydrate, Nakano et al. 57 reported two experimental characteristic bands placed at 1280 and 1386 cm À1 .This characteristic doublet, denoted as Fermi dyad, 58 was also reported by Qin and Kuhs. 26These two bands correspond to the lower (n cÀ ) and upper (n c+ ) levels of Fermi resonance corresponding to vibrational modes with very similar energies, the symmetric stretching vibration (n 1 ) and the first overtone of the bending modes (2n 2 ).This quantum resonance effect, which has been described in detail for the particular case of CO 2 by McCluskey and Stoker, 59 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 (n c+ ) nor the (n 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 (n 1 ) and asymmetric (n 3 ) stretching, and two bending (n 2a and n 2b ) modes.Table 6 shows the obtained shifts (not rescaled in this case) for each of these modes, for isolated CO 2 , 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 2 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. 60 and Pinto et al. 61 ).With the approach followed, an interaction Hamiltonian (H ˆ0) accounting for the difference between harmonic and anharmonic models is then assumed, therefore assigning the unperturbed Hamiltonian (H ˆ0) to the harmonic oscillator and the perturbed (H ˆ= H ˆ0 + H ˆ0) to the anharmonic one.In the following, spectroscopic units (cm À1 ) will be considered as energy units, and states will be labeled with the corresponding energy.With the fundamental state energy E 000 being equal to 0, spectroscopic signals can be directly assigned to excitation energies, denoting as E + the excited state of transition n c+ , E À the excited state of n cÀ , E 100 identified then as n 1 and E 020 as 2n 2 .Therefore, the anharmonic states E + and E À will be the result of the mixing interaction of the unperturbed states E 100 and E 020 .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: or, equivalently, the spectral bands have to hold: 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 .his last value, the energy interaction, was supposed to be constant due to the fact that the experimental Fermi bands are placed 103 cm À1 apart for CO 2 in gas phase and 104 cm À1 in type I hydrate. 63Hence, we used the calculated interaction energy from CO 2 (gas), together with experimental hydrate values of n cÀ and n c+ , to estimate n 1 and n 2 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.
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 n 1 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 n 2 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 n 2 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 n 1 symmetric stretching band corresponding to the T and D cages, as shown in Fig. 6.
This shift between T and D cell frequencies has been a matter of debate, and for instance Sum et al. 24 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 2 was only present in the large cages.In their work, Cao et al. 64 also focused on the same band detected here, the C-O symmetric stretching mode, and reported a difference of 4 cm À1 between computed CO 2 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 À1 , and both bands have been represented in Fig. 6.Nevertheless, and from a practical perspective, Qin and Kuhs 26 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, 65 Kumar et al. 66 and Beeskow-Strauch and Schicks. 67 final consideration can be mentioned concerning the obtained spectral dependence of the CO 2 molecular orientation inside the cell.Udachin et al. 68 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 2 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 À1 , so the effect of this molecular orientation on the results presented here may be considered negligible.
The cited characteristic bands placed around 1500 cm À1 can be outlined as distinctive features of both hydrates studied.For the CH 4 case, the mark is going to appear above this value while for CO 2 it is going to be placed below.In the case of coexistence (mixed CH 4 and CO 2 structures) the expected result is a combination band centered around 1500 cm À1 .Depending on the relative concentration, the band should be displaced above or below this value (CH 4 or CO 2 predominance respectively), which could serve to determine the relative abundance of both types of guest molecules in the hydrate structure.

Conclusions
In the present work some advances in the elucidation of the properties at the molecular level of CH 4 and CO 2 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 2 or CH 4 ), were calculated, showing good agreement with experimental values of clathrates at full occupation.Fermi resonance of CO 2 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.

Fig. 1
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.

Fig. 3
Fig. 3 Interaction potential, in eV, of CH 4 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
Fig.4Interaction potential in eV of CH 4 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.

Fig. 5
Fig. 5 Comparison of the interaction potential (in eV) of CH 4 and CO 2 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 4 , open circles: CO 2 with the axis oriented perpendicularly to the plane of the ring, open squares: CO 2 parallel to the ring.
C-H m ): bond distance between carbon and methane hydrogens, d(CÁ Á ÁO): distance between methane carbon and oxygen atom from water, d(H w Á Á ÁO): distance between water hydrogens and water oxygens, d(H m Á Á ÁO): distance between methane hydrogens and water oxygens, a(H-C-H): angle in the methane molecule.

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

Fig. 7
Fig.7Calculated 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 À1 , and intensity is expressed in arbitrary units.

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

Table 3
Comparison of stretching frequencies (cm À1 ) and C-H bond lengths (Å) as calculated by different methods

Table 4
Comparison between vibrational symmetric stretching frequency values of CH 4 as a guest molecule in the hydrate, CH 4 in aqueous solution (PCM model) and gas CH 4 , all in cm À1 , APD in percentage Cell D Cell T Methane in water Methane (gas) This journal is © the Owner Societies 2015 62 do so, we have first considered experimental values available in the literature for CO 2 in gas state: from the IR band 62 n 2 = 667 cm À1 and Raman levels63n cÀ = 1284 cm À1 and n c+ = 1381 cm À1 , and this way we have estimated n 1 = 1338 cm À1 and E int = À51.998cmÀ1.Considering another different source of experimental values for Raman levels,26n cÀ = 1280 cm À1 and n c+ = 1386 cm À1 , and the same IR band used before, we have obtained alternative values of n 1 = 1332 cm À1 and E int = À52.998cmÀ1.The n 1 values are in remarkably good agreement with the standard experimental reference62of n 1 = 1333 cm À1 .Both estimated values of E int are also in very good agreement with E int = À51.232cm À1 reported by McCluskey and Stoker.

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

Table 5
Experimental Raman signals (expressed in cm À1 ) for the Fermi levels (n cÀ , n c+ ) and the corresponding harmonic bands (n 1 , n 2 ) of CO 2 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 n 2 in those cases aThis journal is © the Owner Societies 2015