Delving into guest-free and He-filled sI and sII clathrate hydrates: a first-principles computational study

Abstract

The dynamics of the formation of a specific clathrate hydrate as well as its thermodynamic transitions depend on the interactions between the trapped molecules and the host water lattice. The molecular-level understanding of the different underlying processes benefits not only the description of the properties of the system, but also allows the development of multiple technological applications such as gas storage, gas separation, energy transport, etc. In this work we investigate the stability of periodic crystalline structures, such as He@sI and He@sII clathrate hydrates by first-principles computations. We consider such host water networks interacting with a guest He atom using selected density functional theory approaches, in order to explore the effects on the encapsulation of a light atom in the sI/sII crystals, by deriving all energy components (guest–water, water–water, guest–guest). Structural properties and energies were first computed by structural relaxations of the He-filled and empty sI/sII unit cells, yielding lattice and compressibility parameters comparable to experimental and theoretical values available for those hydrates. According to the results obtained, the He enclathration in the sI/sII unit cells is a stabilizing process, and both He@sI and He@sII clathrates, considering single cage occupancy, are predicted to be stable whatever the XDM or D4 dispersion correction applied. Our results further reveal that despite the weak underlying interactions the He encapsulation has a rather notable effect on both lattice parameters and energetics, with the He@sII being the most energetically favorable in accord with recent experimental observations.

---

1 Introduction

Clathrate hydrates are crystalline compounds formed by a host water skeleton containing polyhedral cavities of different size and shape, which can be occupied, fully or partially, by diverse guest molecules (noble gases, H₂, N₂, CH₄, CO₂, H₂S,…) depending mostly on the thermodynamic conditions (low temperature and high pressure).^1–3 These ice-like lattices are naturally occurring in the ocean depths, as well as in permafrost regions of Earth.^4,5 Moreover, they are also likely to exist in other planetary bodies, such as Pluto's ocean,⁶ Mars’ ice caps⁷ or Titan's surface.⁸ The stability of these complexes arises from the hydrogen bonds established between water molecules and the van der Waals (vdW) dispersion forces between host and guest species.^9,10

The three most common clathrate hydrates structures adopted naturally are the cubic sI and sII lattices and the hexagonal sH. What distinguishes them is the size, shape and number of water cages present in each of the structures, as well as the periodic crystal structure. Thus, the cubic sI host lattice consists of 46 water molecules per unit cell, forming two pentagonal dodecahedron (small cavity formed by 20 water molecules, 5¹²) and six ellipsoidal-shaped tetrakaidecahedron cages (large cavity formed by 24 water molecules, 5¹²6²). The cubic sII unit cell comprises 136 water molecules located in sixteen small 5¹² cages and eight hexakaidecahedron large cages (composed by 28 water molecules, 5¹²6⁴). Finally, the hexagonal sH unit cell is made up of 34 water molecules distributed in three 5¹², two irregular dodecahedra 4³5⁶6³ (with 20 water molecules) small cages and one large icosahedron 5¹²6⁸ (with 36 water molecules) cage.^5,11,12 Less common clathrate hydrates structures, such as the C₀, C₁, C₂, etc. have been also reported.^13–18 Generally, there is only one guest atom/molecule within each small cage, whereas two or more species can be encapsulated in the larger cages as reported in different multiple-occupancy studies.^19–25

Due to the flexible character of the H-bonds, water presents a rich phase diagram with a great diversity of crystalline phases, both clathrate-like and ice-like. The formation of one structure or another strongly depends, among other things, on the conditions of pressure (P) and temperature (T). Therefore, transitions between different phases can occur or even several structures can coexist.³ Is it important to emphasize the need to understand the stability of these systems, taking into account nuclear quantum effects,^26,27 as well as the mechanisms involved in their formation/dissociation,²⁸ since new structures can be found by both filling with guest molecules²³ or by emptying^16,29 the corresponding hydrate cavities or even the inter-cage molecule migration of one cage to the neighboring ones,^30–32 under certain P–T conditions.³³ This may be useful and relevant, not only from a fundamental point of view, given that new regimes of the phase diagram of water can be investigated, such as low-density ice structures,^34,35 but also from a practical point of view, since promising applications of gas hydrates can be further explored, such as gas storage materials,^36–38 potential energy sources,^39,40 gas separation,^41–43 energy transport^44,45 and desalination techniques.^46–48

So far, a great number of clathrate hydrates filled with different guest species have already been studied.^49–55 In this context, noble gas hydrates are a valuable target due to the simple electronic structure of these inert atoms, together with their unique stability and properties.^10,56 These features make them perfect to create prototype models which may serve to study, for example, the energy variation with the guest size or with the cage occupancy, as well as to predict the stability of other similar systems.^57,58 In this regard, the clathrate hydrates of Ar, Kr, and Xe were discovered long time ago,^59–63 whereas the clathrate hydrates of the lighter noble gases (He and Ne) have been synthesized more recently.^23,29 Focusing our attention on the He atom, it has been observed that it is more likely to form ice-like structures.^64–67 However, as mentioned before, Falenty et al. were able to stabilize the eagerly awaited He@sII clathrate hydrate.²³ Helium also forms binary hydrates with molecules such as hydrogen or tetrahydrofuran in sII clathrate structure.^68,69 Recently, valuable information on benchmark data,^10,70 as well as on the stability of each of the clathrate cavities of the sI, sII structures⁷⁰ have been reported, indicating that the reactions associated with the encapsulation of the He atom inside the 5¹², 5¹²6² and 5¹²6⁴ cages are thermodynamically favored in the range of experimental conditions.²³ Therefore, given the shortage of information available on these systems, and following our previous studies on He–water complexes of different size, in this work we aim to go one step further by studying from a theoretical point of view the stability of periodic crystalline structures of the He@sI and He@sII clathrate hydrates from first-principles approaches.

Since its inception, computational chemistry has made remarkable progress, becoming a useful tool in the study of diverse chemical species, from simple atoms and molecules to complex systems. It yields valuable information about understanding and predicting trends in structures and properties of compounds, as well as about reaction mechanisms or dynamics.⁷¹ The continuous development of quantum theories combined with progressing technology steadily increases the variety of chemical problems that can be treated computationally in terms of reasonable accuracy and affordable computational time.^72,73 Moreover, the recent advances in data science techniques, such as in machine-learning,^74–79 offer new opportunities to efficiently evaluate the properties of molecular systems by successfully accelerating computational algorithms and amplifying the knowledge available from computational chemistry techniques.⁸⁰ Most of the theoretical investigations on noble gas clathrate hydrates have been reported from molecular dynamics simulations, and ab initio quantum chemical approaches.^10,13,56,58 In this context, density functional theory (DFT) is nowadays a suitable and an appropriate tool for studying molecular interactions in condensed phase systems. Modern DFT-D approaches including dispersion corrections improve the description of the nonstandard, non-covalent guest–host interactions, and such schemes are important when one is dealing with gas hydrates. Thus, in this work we performed DFT-D calculations on He-filled and empty sI and sII clathrate hydrates in order to obtain valuable structural and energetic insights, by evaluating their binding through the computation of the two main, host–host and guest–host, interactions. Comparisons with the individual cages,⁷⁰ as well as with experimental and previous theoretical data available in the literature for similar systems are also considered and discussed.

2 Computational details

In Fig. 1 the He@sI and He@sII crystal unit cells together with their corresponding building block cages are shown. Both sI and sII empty crystal structure was initially taken from the 3D crystalline framework given in ref. 81. The oxygen atoms positions have been extracted from the X-ray diffraction experiments, while the hydrogen atoms orientations have been determined from TIP4P water model optimizations, satisfying both the ice rules and give a net zero dipole moment for the unit cell configuration.⁸¹ The cubic sI unit cell has 46 water molecules and Pm3̄n space group symmetry, while the cubic sII hydrate unit cell has 136 water molecules and Fd3̄m space group symmetry. The He atoms are positioned at the center of each cages in the cell, with a H₂O : He ratio, called hydration number, of 46 : 8 and 136 : 24 in the He@sI and He@sII structures, respectively, using the DENEB software package.⁸²

Fig. 1 Individual clathrate-like sI and sII cages (left panel) and bulk He@sI and He@sII clathrate hydrate crystal structures (right panel), with the black box indicating their unit cells. Red color corresponds to oxygen atoms, gray color to hydrogen atoms and magenta color to He atoms.

Electronic structure DFT calculations were performed for both the empty and He-filled sI and sII periodic structures, using the PW86PBE functional as implemented in the LIBXC library⁸³ of the Quantum Espresso (QE) code.^84–86 Dispersion corrections were also considered using the exchange dipole moment (XDM) scheme,^87,88 as implemented in QE, and the D4 term^89,90 as a postprocessing through the DFTD4⁹¹ program. The selection of the PW86PBE functional and the specific dispersion correction schemes is made on the basis of previous investigations on the assessment of various GGA, meta-GGA and hybrid DFT functionals for finite-size clathrate-like cages,^3,10,70,92 as well as on similar gas hydrate systems,^3,12,93 in comparison with benchmark data from DF-MP2, DLPNO-MP2 and DLPNO-CCSD(T) approaches, taking into consideration the computational resources available and the cost of the calculations. Further, we should also note that non-local dispersion effects are also explored through the vdW-DF and vdW-DF2 functionals,^94,95 as implemented in the QE code, although we found an expansion in the lattice parameter for both He-filled sI and sII structures, observed previously in similar clathrate systems,^2,51 and thus they were not considered any further.

The plane-wave/pseudopotential approach within the projector-augmented-waves (PAW) method⁹⁶ was employed using the standard PAW pseudopotentials supplied within QE, with PBE-based pseudopotentials used for the PW86PBE^97,98 computations. Several single-point runs were carried out to check the convergence for the energy cutoff for the plane-wave expansion of the wavefunctions, Ecut_wfc, and for the charge density, Ecut_rho, (see Fig. S1–S4 of the ESI†), as well as the number of k-points in the DFT calculations using the PW86PBE without and with XDM and D4 dispersion corrections for each He-filled and empty sI and sII structures. In this way the Ecut_wfc was selected at 90/80 Ry (1224/1088 eV) and for the Ecut_rho at 480/360 Ry (6530/4898 eV) for the empty and He-filled sI/sII structures, respectively, while a Monkhorst–Pack 2 × 2 × 2/1 × 1 × 1 k-point grid⁹⁹ in the reciprocal space was used per unit cell in the sI/sII systems. As regards the DFT calculations for the isolated H₂O molecule and He atom, they were performed at the Γ-point, in a cubic simulation cell of volume 20 × 20 × 20 ų. Regarding the structure optimizations, the Broyden–Fletcher–Goldfarb–Shanno (BFGS) quasi-newton algorithm was chosen to perform the corresponding geometry relaxation until components of energies and forces on all atoms are below 10⁻⁴ Ry and 10⁻³ Ry bohr⁻¹, respectively.

The cohesive energies per water molecule for the fully occupied He@sI/He@sII structures and the empty sI/sII hydrates were computed as:

where E^He@sI/sII_opt are the total energies of the fully occupied He@sI/sII clathrate hydrates' unit cells obtained by structural relaxation at each lattice constant a, E^H₂O/E^He are the total energies of the isolated ground-state H₂O molecules and He atoms, respectively, and E^sI/sII_opt(a) are the optimized total energies of the empty sI/sII hydrate unit cells at the different values of lattice constant a. N corresponds to 46/136 water molecules in the sI/sII systems, while M stands for 8/24 He atoms in the fully filled sI/sII systems.

Thus, the binding energy per He atom encaged in the empty sI/sII systems is given by,

Finally, the saturation energy corresponding to the total energy difference when the sI/sII hydrates are fully occupied, with one helium per cage, is given by

ΔE_sat = E^He@sI/sII_opt(a₀) − E^sI/sII_opt(a₀) − E^M·He_opt(a₀) (4)

at the corresponding full geometry-optimized configuration for a₀.

3 Results and discussion

3.1 Structural empty and He-filled sI/sII lattice properties

We started by carrying out full geometry optimizations on the fully He-filled and empty sI/sII periodic structures, at a variety of lattice constant, a, values. Total energies E^He@sI/sII_opt and E^sI/sII_opt as a function of a are shown in in Fig. S5 of the ESI,† while the resulting cohesive energies per water molecule as a function of a using the PW86PBE functional with the D4 and XDM dispersion correction schemes (see Tables S1 and S2 for numerical results, ESI†) are plotted in Fig. 2.

Fig. 2 Cohesive energies, in kcal per mol per water molecule, of the fully occupied He@sI clathrate and the empty sI hydrate (left panel), and of the fully occupied He@sII clathrate and the empty sII hydrate (right panel), as a function of the lattice constant a in Å. Solid lines corresponds to the MEOS fits, while the values of equilibrium lattice constants are drawn with dashed (this work) and dotted (DMC value from ref. 51 for sI, and experimental ones from ref. 23 and 29 for He@sII and sII) lines.

The curves present similar behavior in both panels, although there are certain differences that should be mentioned. First of all, the lattice parameters associated with the sI crystal are smaller than those corresponding to the sII crystal, as a consequence of the smaller size of the sI lattice. In both systems, the PW86PBE-XDM functional predicts more energetically favorable lattices than the same functional with the D4 dispersion correction. Nevertheless, the energy difference between the He-filled (He@sI/He@sII) and empty (sI/sII) systems, when considering both D4 and XDM schemes, is higher for the sI crystals than for the sII ones. For example, at a = 11.8 Å (sI/He@sI systems) and 16.8 Å (sII/He@sII systems) this difference is around 0.88 and 0.83 kcal mol⁻¹ for the He@sI and sI curves, respectively, and 0.67 and 0.66 kcal mol⁻¹ for the He@sII and sII, respectively. Therefore, we can extract two conclusions from this analysis: on the one hand, the energy difference between the D4 and XDM energies is higher for the He-filled and empty sI systems, and on the other hand, the curves corresponding to the empty structures for both sI and sII systems are energetically closer than those of the corresponding He-filled lattices. Another interesting feature that should be mentioned is that the inclusion of the He atom in both sI and sII lattices generates full He-filled structures with lower energy than the sI and sII ones. In particular, an energy lowering by around 0.09 kcal per mol per H₂O molecule is obtained for both He@sI and He@sII from the PW86PBE-XDM calculations, while smaller values of around 0.03–0.06 kcal per mol per H₂O molecule is estimated from the PW86PBE-D4 ones. Thus, such energy differences (per water molecule) between the He-filled and empty systems curves are approximately equal when one considers the XDM dispersion scheme, but higher for the sII–He@sII curves when the D4 correction is included.

The next step in our analysis consists in obtaining the equilibrium lattice constants, a₀, and the corresponding ΔE^He@sI/sII_coh(a₀) and ΔE^sI/sII_coh(a₀) energies. The calculated total E^He@sI/sII(a) and E^sI/sII(a) energies as a function of the lattice parameter are fitted to the Murnaghan's (see solid and dashed lines in Fig. 2) and Vinet's equations of state (MEOS, VEOS).^100,101 Similar trends are found for both MEOS and VEOS fits, with the resulting equilibrium lattice parameters, as well as the bulk moduli (bulk modulus, B₀, and bulk modulus pressure derivative, , at zero pressure) listed in Table 1.

Table 1 Murnaghan's/Vinet's equation of state fitting results obtained from the PW86PBE-XDM and -D4 calculations for both fully He-filled and empty sI/sII clathrate hydrates

	PW86PBE-XDM				PW86PBE-D4
	He@sI	sI	He@sII	sII	He@sI	sI	He@sII	sII
a 0 (Å)	11.741/11.745	11.753/11.757	16.903/16.914	16.912/16.923	11.840/11.844	11.838/11.841	16.950/16.961	16.943/16.953
B 0 (GPa)	12.14/12.79	11.90/12.50	15.85/16.26	15.59/15.96	10.63/11.34	10.69/11.41	15.08/15.57	15.29/15.87
	6.05/5.71	5.89/5.58	8.38/7.45	8.35/7.42	6.09/5.79	6.13/5.83	8.86/7.97	9.33/8.34

---

In the present DFT calculations, the a₀ values are computed through potential energy optimizations for the empty and fully filled He@sI/He@sII crystal cells, and thus temperature or pressure, as well as vibrational zero-point effects are not taken into account. However, since the experimental data available in the literature reported at specific P–T conditions, we should only consider comparisons with those results at zero pressure and temperature. Thus, for the empty sI clathrate we compare the a₀ obtained from the present MEOS fits with the recent theoretical value of a₀ = 11.77 Å reported from MEOS fits to diffusion Monte Carlo (DMC) calculations, while for the empty sII hydrate an equilibrium lattice parameter a₀ = 17.1251 Å has been obtained by fitting to a polynomial T-dependence expression the data from the neutron diffraction experiments on the deuterated hydrate.²⁹ If we compare these results with the present values, we get relative errors of 0.144% and 0.578% at PW86PBE-XDM and PW86PBE-D4 levels, respectively, for the empty sI system; whereas the error increases to 1.244% and 1.063%, respectively, for the sII crystal. We should note that the deuterated hydrates have somewhat weaker D-bonds that could result to slightly longer cell parameters than those of the H-bonding ones, as one can observe such isotopic H/D effect between the calculated and experimentally measured a₀ value for the sII hydrates. As we can also see, the equilibrium lattice parameter computed including XDM dispersion correction is closer to the reference theoretical value for the sI hydrate, while the calculated value with the D4 dispersion has the smallest relative error for the sII structure.

Regarding the case of the He-filled sI and sII crystals, until recently studies were limited to He atom encapsulation in ice-frameworks. Nevertheless, in 2018 the He@sII structure has been synthesized in the laboratory,²³ and for such partially filled clathrates an a₀ value of 17.0763 Å has been reported for samples at 100 MPa, and a₀ = 17.1163 Å at ambient pressure for a temperature of 80 K. We should highlight here that the changes reported²³ in the lattice parameters of the empty sII and partially filled He@sII count 0.01 Å, at the same P–T conditions (both structures with deuterated water framework). Unfortunately, there is no more information on lattice parameters at low temperatures for such He clathrates, although an extrapolation by considering the analytical expression of the thermal expansivity is available for the Ne@sII in ref. 29. The polynomial, a(T) = A₀ + A₁T + A₂T² + A₃T³ + A₄T⁴, expression²⁹ extrapolates the measured lattice parameter from 145 K to 0 K, and estimates of a₀ = 17.1025 and 17.1251 Å have been reported for the Ne@sII and sII, respectively. The a₀ value for the Ne@sII is shorter by less than 0.1% from the experimental value at 80 K, thus we may suggest that a similar change might be expected for the He@sII system. Indeed, such experimental value is very close to the 16.903/16.950 Å obtained from the PW86PBE-XDM/D4 calculations for the He@sII, and it is by 0.01 Å smaller than those of the empty sII at the PW86PBE-XDM level of calculation and 0.01 Å larger at the PW86PBE-D4 (see Table 1), with the XDM predicting the lattice changes in the same direction with the experimental data. However, we should be rather cautious as multiple or partial occupancy observed in the experimental study, as well as deuterated effects could also affect these results.

In the absence of further He data, we also compare our results with those reported in the literature for similar systems. For example, values of a₀ = 17.103 Å for Ne@sII,²⁹ and a₀ = 11.835 Å for Xe@sI¹⁰² for rare-gas clathrates. Continuing with small molecules clathrates, a₀ = 17.083 and 11.811 Å for the N₂@sII²⁹ and N₂@sI,¹⁰³ respectively, while for systems encapsulating CO, such as CO@sI/sII a₀ = 11.819/17.063 Å have been reported² from neutron powder diffraction experiments. We can also compare our outcomes with those available for larger encapsulated molecules, such as a₀ = 11.835 Å for ethylene oxide sI hydrate, and a₀ = 17.130 Å for tetrahydrofuran sII hydrate.^104,105 In most cases, the difference with respect to the computed He@sI and He@sII lattice parameters is less than 0.08 and 0.17 Å, respectively. One can see that all these lattice parameters are larger than those of He-filled systems, as expected from species with stronger interactions.

In turn, we have studied the effects of pressure on both empty sI/sII and fully filled He@sI/sII structures. The corresponding MEOS of each structure are displayed in the a–P diagrams of Fig. 3. The bulk moduli, both B₀ and were initially estimated from the corresponding MEOS fits, assuming linear dependence of the bulk modulus with pressure, , with B₀ being constant and valid for the range . Using an initial guess of , the B₀ values listed in Table 1 are obtained. The PW86PBE-XDM predicts higher values of B₀ and for the He-filled structures compared with empty hydrates, thus indicating that the crystal without guest molecules is more compressible, whereas the PW86PBE-D4 predicts lower values for the same parameters in the He-filled structures, concluding that the empty systems are more resistant to compression. Therefore, the PW86PBE functional with XDM and D4 dispersion corrections produces opposite results. To clarify the results related to the structural parameters, we have also fitted the values of the total E^He@sI/sII_opt(a) and E^sI/sII_opt(a) energies considering the VEOS, which is applicable for all classes of solids in compression.¹⁰¹ With this idea, we have found that PW86PBE-XDM predicts that empty systems (sI and sII) are more compressible than fully occupied clathrates (He@sI/sII) having a smaller bulk modulus, while, PW86PBE-D4 predicts that filled systems have higher compressibility than empty hydrates. Although, the trend observed using the MEOS is confirmed by using the VEOS (see corresponding values in the Table 1), it should be noted that such overall claims on the bulk modulus are weakened by accuracy issues in the theoretical bulk modulus, as testified by the fact that by changing the EOS form their value difference is larger that the differences obtained between filled and empty hydrates systems employing XDM or D4 dispersion corrections. Next, in Fig. 3, one can see that the higher the pressure, the shorter the lattice parameter or what is the same, the He-filled and empty sI and sII hydrates suffer greater compression as pressure increases. The curves corresponding to He-filled and empty clathrates are almost overlapped when XDM or D4 dispersion corrections are applied.

Fig. 3 Pressure effects on the lattice constant a of the He-filled and empty sI (left panel) and sII clathrate hydrates (right panel). Solid lines correspond to the MEOS fits (see Table 1). The zero-pressure lattice parameter values are indicated for the sI/He@sI, while the values at ambient pressure (10⁻⁴ GPa) are reported for the sII/He@sII.

The bulk modulus obtained for the He@sI/sII clathrates at PW86PBE-XDM/D4 level (Table 1) are smaller than those obtained for other systems.^3,12 For example, for the CO₂@sI hydrate values of 13.64/12.53 GPa for the bulk modulus have been recently reported, while in the He@sI the values obtained are 12.140/10.632 GPa for the same functionals. In the case of CO₂@sII, these values are 17.149/16.053 GPa, whereas for the He@sII hydrate we obtained 15.851/15.084 GPa. According to this, it is possible to see that for both structures (sI and sII), the bulk modulus values for the clathrate hydrates with helium encapsulated are smaller than in the cases with the larger CO₂ molecule, that indicates the higher compressibility of the He@sI/sII clathrates compared to the CO₂@sI/sII ones. We should note that all comparisons of the structure parameters and compressibility of He clathrates with experimental data on other clathrates is informative on the strength of the cage–host interaction, it is not meant as a validation of the calculations.

3.2 Stability of the He@sI and He@sII clathrates

The He clathrate hydrate has been recently observed only in the sII structure,²³ and prediction of stability issues involves the calculation of guest–host and host–host interactions in the such crystal lattices. All energies were computed for a single occupancy of all cages with hydration number of 5.75 and 5.67 for the He@sI and He@sII clathrates, respectively, using the PW86-XDM/D4 DFT-D functionals.

Total and cohesive energy values are reported above, and in order to better understand these results, we have decomposed the total cohesive energy into contributions arising from the He atom binding to the empty sI or sII hydrate (see eqn (3)), ΔE^He@sI/sII, the guest–host interaction part, and the cohesive energy of the empty sI and sII hydrate (see eqn (2)), ΔE^sI/sII_coh, corresponding to the host–host interaction component. In Fig. 4 these binding energy values of He to the sI/sII hydrates, as well as the cohesive energies of the empty sI/sII crystals at their corresponding equilibrium lattice constant, a₀, obtained from the present DFT-D computations, are displayed. The binding energies (per He atom) characterizing the guest–host interactions, and the cohesive energies per water molecule are consistent with energies previously reported for the individual He-filled and empty 5¹² (present in sI and sII), 5¹²6² (present in sI) and 5¹²6⁴ (present in sII) cages from DLPNO-CCSD(T)/AVTZ reference calculations,⁷⁰ and are plotted in Fig. 4 (see vertical dashed and dashed-dotted lines), together with the reference DMC value for the bulk sI hydrate, ΔE^sI_coh, (see vertical solid line) available in the literature.⁵¹

Fig. 4 Binding (per He atom) and cohesive (per water molecule) energies, in kcal mol⁻¹, of the He-filled and empty sI and sII clathrate hydrates, respectively, obtained from the present PW86PBE-XDM/D4 computations. Vertical dashed/dashed-dotted lines correspond to interaction and cohesive (per water molecule) energies, in kcal mol⁻¹, of the individual He-filled/empty sI and sII clathrate-like cages, respectively, computed at DLPNO-CCSD(T)/AVTZ level,⁷⁰ while the reference DMC value from ref. 51 for the empty sI hydrate is displayed by a solid color (cyan) line.

One can see that the PW86PBE-XDM/D4 values underestimate the bulk sI formation energies compared to the references DMC data,⁵¹ while the DLPNO-CCSD(T) values for finite-size cages should be considered as lower energy limits in all bulk hydrates studied here. The PW86PBE-XDM calculations overestimate binding energies of the He@sII clathrates, as well as formation energies of empty sI and sII hydrates, and underestimate binding of the He@sI structure compared to the PW86PBE-D4 results. Another aspect which stands out is the fact that for the He@sII crystal the binding is stronger than for the He@sI one by approximately 17 kcal mol⁻¹, whatever the dispersion correction is. On the contrary, the empty sI hydrate is predicted to have lower cohesive energy than the sII one by approximately 3 kcal mol⁻¹. Although the difference in the cohesive energies of the water framework is smaller on a per molecule basis, compared to those in the He binding energies, apparently the H₂O : He ratio means that small errors/differences in water–water interactions are much affected compared to the larger values in binding per He atom. Nevertheless, one can see that the PW86PBE functional with both XDM and D4 dispersion schemes exhibit the same qualitative behavior.

In turn we also present the results of energetics involved in the total energy release after saturation of all 8 and 24 cages in the bulk sI and sII hydrates. The saturation energy, ΔE_sat, is obtained as the total energy difference between the empty sI/sII and free He atoms to full occupied He@sI/sII clathrate (see eqn (4)). We found that single He filling of all 24 sII cages is energetically more favorable than filling all 8 cavities of the sI structure. The values of the corresponding saturation energies are −3.99/−1.45 and −11.95/−7.22 kcal mol⁻¹ for He@sI and He@sII, respectively, calculated at PW86PBE-XDM/D4 (see Fig. 5). Again we can see that results obtained from both XDM and D4 dispersion energies are qualitatively similar, predicting that the full single He occupation stabilizes the sII and sI structures, He@sII is clearly more stable than He@sI, although the XDM saturation energies being much lower than those obtained from the D4 scheme.

Fig. 5 Saturation energies, ΔE_sat, in kcal mol⁻¹, for fully single cage occupation He@sI and He@sII clathrate hydrates, as computed from the present PW86PBE-XDM/D4 calculations.

Such conclusion could be supported by considering the size and number of cages forming the bulk sI and sII unit cells. The fact that the sII unit cell has a larger number of small 5¹² than the sI one, and further the size of the 5¹²6² and 5¹²6⁴ cages forming the sI and sII clathrates, respectively, could affect the He binding in these cages, as have been also found from the DLPNO-CCSD(T) calculations of these individual cages⁷⁰ (see also Fig. 4). Similar results have been also reported for the CO and CO₂ molecules and their stabilizing role in the sII and sI structures,^2,3,12 considering both single and multiple cages occupancy. Overall, our calculations show that the full single He occupation of the sII hydrate is energetically more favorable than the sI one, although entropy and enthalpy effects need to be considered too.

4 Summary and conclusions

By performing first-principles DFT-D calculations, we study the effects on the encapsulation of a He atom in sI/sII clathrate hydrates considering the entire periodic 3D crystal lattices. Such light rare-gas@hydrate are challenging targets, as they are suitable for checking the accuracy of computational quantum chemistry approaches on nonstandard non-covalent interactions, and at the same time, the output could contribute to the current emerging research on stabilization of new low-density ice polymorphs, such as ice-XVI or sIII. Thus, on the basis of previous benchmark studies, we have employed the PW86PBE functional with the XDM and D4 schemes for describing the water–water (host–host) and He–water (guest–host) interactions, and evaluating the role of dispersion forces in the complete periodic lattices. Cohesive and binding energies, through structural geometry relaxation DFT-D calculations of the fully He-filled and empty sI/sII unit cells were determined, and their structural properties and stability were discussed. Our results show that the stability could be properly determined by taking into account dispersion corrected schemes such as D4 or XDM ones, with the PW86PBE-D4/XDM functional yielding lattice constants and compressibility parameters comparable to recent experimental and previous theoretical values available for the bulk sII/He@sII and sI structures, respectively. We have highlighted that for an overall reasonable description both the H-bonded water framework and the vdW bound He interactions should be accurately represented, emphasizing the importance of H₂O : He ratio in the error analysis of water–water/He–water interactions.

The present computations reveal that both processes for full He occupation of the sI and sII hydrates cages are energetically favored, with the encapsulation of a single He atom per cage playing a stabilizing role in the formation of He@sI and He@sII clathrates, whatever XDM or D4 dispersion scheme used. As shown, guest–host interactions are far from being negligible contributing to the stabilization of such He-filled structures. Total saturation energies of −3.99/−1.45 and −11.95/−7.22 kcal mol⁻¹ were obtained from the present PW86PBE-XDM/D4 calculations for the He@sI and He@sII crystal structures, respectively.

A complete description of the formation of He clathrate structures requires taking into account thermodynamic factors, such as temperature and pressure, as well as nuclear quantum effects.^26,27 The incorporation of enthalpy and entropy effects could determine whether the calculated stabilization energies are related to observed phase diagram transformations.^{13,33,67,70,106,107} Understanding the underlying factors involved in the stability of such hydrates allows the extension of this work in this direction. Thus the next step consists to consider multiple cage occupancy similar to those established by diffraction experiments²³ with up to 4 He atoms in the large sII cages and one guest in the small ones, for a more direct comparison with measurements on structural properties. Admittedly, the large sII unit cell of the He clathrate structure as the number of guests increases could become computationally expensive, although tractable for few of them. Moreover, investigations involving the filling with heavier rare-gas, such as Ne, atoms, as well as the stability of other clathrate structures (e.g. He@sH, He@C₀,…) or ice-polymorphs could further contribute to delve deeper into the subject, and will be addressed in future studies.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We would like to thank Raúl Rodríguez-Segundo for assistance with DENEB-Atelgraphics software. The authors thank to Centro de Cálculo del IFF, SGAI (CSIC) and CESGA for allocation of computer time. This work has been supported by MINECO grant no. FIS2017-83157-P, MICINN grant no. PID2020-114654GB-I00, Comunidad de Madrid grant no. IND2018/TIC-9467, and COST Action CA18212(MD-GAS).

Notes and references

1. S. Patchkovskii and J. S. Tse, Proc. Natl. Acad. Sci. U. S. A., 2003, 100, 14645–14650.
2. C. Pétuya, L. Martin-Gondre, P. Aurel, F. Damay and A. Desmedt, J. Chem. Phys., 2019, 150, 184705.
3. A. Cabrera-Ramírez, R. Yanes-Rodríguez and R. Prosmiti, J. Chem. Phys., 2021, 154, 044301.
4. F. J.-A. L. Cruz and J. P.-B. Mota, Phys. Chem. Chem. Phys., 2021, 23, 16033.
5. J. P. Heindel, K. M. Herman, E. Aprá and S. S. Xantheas, J. Phys. Chem. Lett., 2021, 12, 7574–7582.
6. S. Kamata, F. Nimmo, Y. Sekine, K. Kuramoto, N. Noguchi, J. Kimura and A. Tani, Nat. Geosci., 2019, 12, 407–410.
7. D. Ambuehl and M. E. Madden, Icarus, 2014, 234, 45–52.
8. T. H. Vu, M. Choukroun, C. Sotin, Y. Muñoz Iglesias and H. E. Maynard-Casely, Geophys. Res. Lett., 2020, 47, e2019GL086265.
9. J. Lee and J. W. Kenney III, in Solidification, ed. A. E. Ares, IntechOpen, Rijeka, 2018, ch. 7.
10. R. Yanes-Rodríguez, D. J. Arismendi-Arrieta and R. Prosmiti, J. Chem. Inf. Model., 2020, 60, 3043–3056.
11. E. Kim and Y. Seo, Chem. Eng. Sci., 2019, 359, 775–778.
12. A. Cabrera-Ramírez, D. J. Arismendi-Arrieta, A. Valdés and R. Prosmiti, ChemPhysChem, 2021, 22, 359–369.
13. P. Teeratchanan and A. Hermann, J. Chem. Phys., 2015, 143, 154507.
14. D. Londono, J. L. Finney and W. F. Kuhs, J. Chem. Phys., 1992, 97, 547–552.
15. J. S. Loveday and R. J. Nelmes, Phys. Chem. Chem. Phys., 2008, 10, 913–1068.
16. L. del Rosso, M. Celli and L. Ulivi, Nat. Commun., 2016, 7, 13394.
17. T. A. Strobel, M. Somayazulu, S. V. Sinogeikin, P. Dera and R. J. Hemley, J. Am. Chem. Soc., 2016, 138, 13786–13789.
18. D. Amos, M.-E. Donelly, P. Teeratchanan, C. Bull, A. Falenty, W. Kuhs, A. Hermann and J. Loveday, J. Phys. Chem. Lett., 2017, 8, 4295–4299.
19. H. Itoh and J. S. Tse, J. Chem. Phys., 2001, 115, 9414.
20. W. L. Mao, H.-k Mao, A. F. Goncharov, V. V. Struzhkin, Q. Guo, J. Hu, J. Shu, R. J. Hemley, M. Somayazulu and Y. Zhao, Science, 2002, 297, 2247–2249.
21. H. Lu, J. Wang, C. Liu, C. I. Ratcliffe, U. Becker, R. Kumar and J. Ripmeester, J. Am. Chem. Soc., 2012, 134, 9160–9162.
22. X. Cao, Y. Su, Y. Liu, J. Zhao and C. Liu, J. Phys. Chem. A, 2014, 118, 215–222.
23. W. F. Kuhs, T. C. Hansen and A. Falenty, J. Phys. Chem. Lett., 2018, 9, 3194–3198.
24. R. Ma, H. Zhong, J. Liu, J. Zhong, Y. Yan, J. Zhang and J. Xu, Processes, 2019, 7, 699.
25. C. Métais, C. Petuya, S. Espert, J. Ollivier, L. Martin-Gondre and A. Desmedt, J. Phys. Chem. C, 2021, 125, 6433–6441.
26. A. Witt, F. Sebastianelli, M. E. Tuckerman and Z. Bacić, J. Phys. Chem. C, 2010, 114, 20775–20782.
27. A. Valdés, D. J. Arismendi-Arrieta and R. Prosmiti, J. Phys. Chem. C, 2015, 119, 3945–3956.
28. J. R. Zhong, X. Y. Zeng, F. H. Zhou, Q. D. Ran, C. Y. Sun, R. Q. Zhong, L. Y. Yang, G. J. Chen and C. A. Koh, Sci. Rep., 2016, 6, 38855.
29. A. Falenty, T. C. Hansen and W. F. Kuhs, Nature, 2014, 516, 231–233.
30. H. Cao, N. J. English and J. M.-D. MacElroy, J. Chem. Phys., 2013, 138, 094507.
31. C. Burnham and N. J. English, J. Phys. Chem. C, 2016, 120, 16561–16567.
32. J. R. Cendagorta, A. Powers, T. J.-H. Hele, O. Marsalek, Z. Bacić and M. E. Tuckerman, Phys. Chem. Chem. Phys., 2016, 18, 32169–32177.
33. K. V. Gets, R. K. Zhdanov, Y. Y. Bozhko and V. R. Belosludov, J. Phys. Chem. C, 2021, 125, 15659–15663.
34. Y. Liu and L. Ojamäe, Phys. Chem. Chem. Phys., 2018, 20, 8333–8340.
35. Y. Huang, C. Zhu, L. Wang, X. Cao, Y. Su, X. Jiang, S. Meng, J. Zhao and X. C. Zeng, Sci. Adv., 2016, 2, e1501010.
36. V. R. Belosludov, Y. Y. Bozhko, O. S. Gets, K. V. Subbotin and Y. Kawazoe, J. Phys.: Conf. Ser., 2018, 1128, 012031.
37. C. Cheng, F. Wang, Y. Tian, X. Wu, J. Zheng, J. Zhang, L. Li, P. Yang and J. Zhao, Renewable Sustainable Energy Rev., 2020, 117, 109492.
38. A. K. Both, Y. Gao, X. C. Zeng and C. L. Cheung, Nanoscale, 2021, 13, 7447.
39. Z. R. Chong, S. H.-B. Yang, P. Babu, P. Linga and X. S. Li, Appl. Energy, 2016, 162, 1633–1652.
40. Y. H. Ahn, S. Moon, D. Y. Koh, S. Hong, H. Lee, J. W. Lee and Y. Park, Energy Storage Mater., 2020, 24, 655–661.
41. V. R. Avula, P. Gupta, R. L. Gardas and J. S. Sangwai, Asia-Pac. J. Chem. Eng., 2017, 12, 709–722.
42. P. Linga, R. Kumar and P. Englezos, J. Hazard. Mater., 2007, 149, 625–629.
43. P. Babu, P. Linga, R. Kumar and P. Englezos, Energy, 2015, 85, 261–279.
44. J. S. Gudmundsson, Method for production of gas hydrates for transportation and storage, US Pat., 5536893, 1996.
45. N.-J. Kim, J. H. Lee, Y. S. Cho and W. Chun, Energy, 2010, 35, 2717–2722.
46. P. Babu, A. Nambiar, T. He, I. A. Karimi, J. D. Lee, P. Englezos and P. Linga, ACS Sustainable Chem. Eng., 2018, 6, 8093–8107.
47. J. Lee and Y. Seo, Chem. Eng. J., 2019, 375, 121974.
48. N. Thakre and A. K. Jana, Renewable Sustainable Energy Rev., 2021, 135, 110150.
49. S. Y. Willow and S. S. Xantheas, Chem. Phys. Lett., 2012, 525–526, 13–18.
50. J. Zhu, S. Du, X. Yu, J. Zhang, H. Xu, S. C. Vogel, T. C. Germann, J. S. Francisco, F. Izumi, K. Momma, Y. Kawamura, C. Jin and Y. Zhao, Nat. Commun., 2014, 5, 4128.
51. S. J. Cox, M. D. Towler, D. Alfè and A. Michaelides, J. Chem. Phys., 2014, 140, 174703.
52. S. Alavi and J. A. Ripmeester, Mol. Simul., 2017, 43, 808–820.
53. U. Ranieri, M. M. Koza, W. F. Kuhs, S. Klotz, A. Falenty, P. Gillet and L. E. Bove, Nat. Commun., 2017, 8, 1076.
54. J. Ghosh, R. R.-J. Methikkalam, R. G. Bhuin, G. Ragupathy, N. Choudhary, R. Kumar and T. Pradeep, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 1526–1531.
55. J. Zheng, Z. R. Chong, M. F. Qureshi and P. Linga, Energy Fuels, 2020, 34, 10529–10546.
56. S. Mondal and P. K. Chattaraj, Phys. Chem. Chem. Phys., 2014, 16, 17943–17954.
57. S. P. Kaur and C. N. Ramachandran, Mol. Phys., 2018, 116, 54–63.
58. P. H.-B. Brant Carvalho, A. Mace, O. Andersson, C. A. Tulk, J. Molaison, A. P. Lyubartsev, I. M. Nangoi, A. A. Leitão and U. Häussermann, Phys. Rev. B, 2021, 103, 064205.
59. P. Villard, C. R. Hebd. Seances Acad. Sci., 1896, 123, 377–379.
60. P. Villard, C. R. Hebd. Seances Acad. Sci., 1896, 123, 377–379.
61. R. de Forcrand, C. R. Hebd. Seances Acad. Sci., 1923, 176, 355–358.
62. R. M. Barrer and A. V.-J. Edge, Proc. R. Soc. London, Ser. A, 1967, 300, 1–24.
63. Y. A. Dyadin, E. Y. Aladko, A. Y. Manakov, F. V. Zhurko, T. V. Mikina, V. Y. Komarov and E. V. Grachev, J. Struct. Chem., 1999, 40, 790–795.
64. D. Londono, W. F. Kuhs and J. L. Finney, Nature, 1998, 332, 141–142.
65. C. Lobban and J. L. Finney, J. Chem. Phys., 2002, 117, 3928.
66. A. V. Ildyakov, A. Y. Manakov, E. Y. Aladko, V. I. Kosyakov and V. A. Shestakov, J. Phys. Chem. B, 2013, 117, 7756–7762.
67. R. V. Belosludov, Y. Y. Bozhko, O. S. Subbotin, V. R. Belosludov, H. Mizuseki, Y. Kawazoe and V. M. Fomin, J. Phys. Chem. C, 2014, 118, 2587–2593.
68. N. I. Papadimitriou, I. N. Tsimpanogiannis, A. K. Stubos, A. Martin, L. J. Rovetto and C. J. Peters, J. Phys. Chem. Lett., 2010, 1, 1014–1017.
69. U. Ranieri, M. M. Koza, W. F. Kuhs, R. Gaal, S. Klotz, A. Falenty, D. Wallacher, J. Ollivier, P. Gillet and L. E. Bove, J. Phys. Chem. C, 2019, 123, 1888–1903.
70. R. Yanes-Rodríguez and R. Prosmiti, Phys. Chem. Chem. Phys., 2022, 24, 1475–1485.
71. R. A. Friesner and B. J. Berne, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 6648–6653.
72. Y. Ozaki, K. B. Beć, Y. Morisawa, S. Yamamoto, I. Tanabe, C. W. Huck and T. S. Hofer, Chem. Soc. Rev., 2021, 50, 10917–10954.
73. I. Tuvi-Arad, Isr. J. Chem., 2021, 61, 1–10.
74. M. Gastegger, J. Behler and P. Marquetand, Chem. Sci., 2017, 8, 6924–6935.
75. G. Carleo, Y. Nomura and M. Imada, Nat. Commun., 2018, 9, 5322.
76. F. Noé, A. Tkatchenko, K.-R. Müller and C. Clementi, Annu. Rev. Phys. Chem., 2020, 71, 361–390.
77. J. Westermayr and P. Marquetand, Chem. Rev., 2021, 121, 9873–9926.
78. Q. Ai, V. Bhat, S. M. Ryno, K. Jarolimek, P. Sornberger, A. Smith, M. M. Haley, J. E. Anthony and C. Risko, J. Chem. Phys., 2021, 154, 174705.
79. M. Ceriotti, C. Clementi and O. A. von Lilienfeld, Chem. Rev., 2021, 121, 9719–9721.
80. J. A. Keith, V. Vassilev-Galindo, B. Cheng, S. Chmiela, M. Gastegger, K. R. Müller and A. Tkatchenko, Chem. Rev., 2021, 121, 9816–9872.
81. F. Takeuchi, M. Hiratsuka, R. Ohmura, S. Alavi, A. K. Sum and K. Yasuoka, J. Chem. Phys., 2013, 138, 124504.
82. DENEB 1.30 beta: the Nanotechnology Software by Atelgraphics, 2021, https://www.atelgraphics.com.
83. M. A.-L. Marques, M. J.-T. Oliveira and T. Burnus, Comput. Phys. Commun., 2012, 183, 2272–2281.
84. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
85. P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Kücükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. Otero-de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu and S. Baroni, J. Phys.: Condens. Matter, 2017, 29, 465901.
86. P. Giannozzi, O. Baseggio, P. Bonfá, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru and S. Baroni, J. Chem. Phys., 2020, 152, 154105.
87. A. D. Becke and E. R. Johnson, J. Chem. Phys., 2007, 127, 154108.
88. A. Otero de la Roza and E. R. Johnson, J. Chem. Phys., 2012, 136, 174109.
89. E. Caldeweyher, C. Bannwarth and S. Grimme, J. Chem. Phys., 2017, 147, 034112.
90. E. Caldeweyher, S. Ehlert, A. Hansen, H. Neugebauer, S. Spicher, C. Bannwarth and S. Grimme, J. Chem. Phys., 2019, 150, 154122.
91. Universität-Bonn, D4- A Generally Applicable Atomic-Charge Dependent London Dispersion Correction, https://www.chemie.uni-bonn.de/pctc/mulliken-center/software/dftd4.
92. A. Cabrera-Ramírez, D. J. Arismendi-Arrieta, A. Valdés and R. Prosmiti, Chem. Phys. Chem., 2020, 21, 2618–2628.
93. F. Izquierdo-Ruiz, A. Otero-de-la Roza, J. Contreras-Garciía, O. Prieto-Ballesteros and J. M. Recio, Materials, 2016, 9, 777.
94. M. Dion, H. Rydberg, E. Schröder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2004, 92, 246401.
95. K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist and D. C. Langreth, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 081101.
96. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
97. J. P. Perdew and W. Yue, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 8800.
98. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
99. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Solid State, 1976, 13, 5188.
100. F. D. Murnaghan, Proc. Natl. Acad. Sci. U. S. A., 1944, 30, 244–247.
101. P. Vinet, J. Ferrante, J. Smith and J. Rose, J. Phys. C: Solid State Phys., 1986, 19, L467.
102. T. C. Hansen, A. Falenty and W. F. Kuhs, J. Chem. Phys., 2016, 144, 054301.
103. V. Vins, A. Jäger, S. Hielscher, R. Span, J. Hrubý and C. Breitkopf, EPJ Web Conf., 2017, 143, 02141.
104. J. S. Tse, W. R. McKinnon and M. Marchi, J. Phys. Chem., 1987, 91, 4188–4193.
105. J. S. Tse, J. Phys., Colloq., 1987, 48, 543–549.
106. A. Vítek, D. J. Arismendi-Arrieta, R. Rodríguez-Cantano, R. Prosmiti, P. Villarreal, R. Kalus and G. Delgado-Barrio, Phys. Chem. Chem. Phys., 2015, 17, 8792–8801.
107. A. Vítek, D. J. Arismendi-Arrieta, M. Sarmanová, R. Kalus and R. Prosmiti, J. Phys. Chem. A, 2020, 124, 4036–4047.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2cp00701k

---