A method distinguishing between guest molecules that can form sI, sII, and sH hydrogen clathrates

A new method based on free energy calculations is proposed for discriminating between promoters that can form sI, sII, and sH hydrogen clathrates. The method is validated by comparing results of this computational approach to known experimental data reporting which type of clathrate structure is formed with the help of a particular promoter molecule. Good agreement is found. The free energy results confirm a well-known simple rule of thumb based on the van der Waals volume of the promoter molecules to distinguish between potential sH or sII hydrate formers. The method predicts an unusually strong temperature dependence of the stability of hydrogen clathrate hydrates for some hitherto untested promoter molecules, which can be exploited if easy release of hydrogen is required.


I. Introduction
Clathrate hydrates are molecular crystals consisting of a hydrogen-bonded water framework with enclathrated guest molecules 1 stabilizing the compounds. Clathrates have been proposed as potential systems for hydrogen storage, because they can hold gas at relatively high density in molecular form. However, extreme pressures of $220 MPa are required for the formation of stable pure hydrogen clathrates, 2,3 but the formation pressure can be signicantly reduced by putting a second guest molecule, a so-called promoter, in the large cavities of hydrogen clathrate hydrates. 4 There are three important structures of gas hydrates: cubic structure I (sI) with two small and six large cages, cubic structure II (sII) with sixteen small and eight large cages, and hexagonal structure H (sH) with three small, two medium and one large cavity in the unit cell. The latter structure is particular promising for hydrogen storage since it allows the highest density storage when both the small and medium cages are occupied by H 2 . The aim of this work is to present and validate a predictive screening method for discriminating between promoters that form hydrogen sI, sII, and sH clathrates. Since the most prospective candidates for hydrogen storage material are structures II and H, in this study we focus mainly on sII and sH hydrogen clathrate hydrates.
Many organic water-soluble and -insoluble compounds have been identied experimentally as efficient stabilizing promoter molecules for sI, sII, and sH hydrogen clathrate formation. Frankcombe et al. 5 presented a computational method to discriminate between guest molecules that form sII clathrate and molecules capable (together with a help gas) to form sH clathrate hydrate. The screening method is based on calculations of the difference between the energies of sII and sH clathrate formation. A clathrate formation energy is dened with reference to the total potential energies of the corresponding clathrate structure (sII or sH) and of the potential energy of hexagonal ice structure. The clathrate formation energy dened by Frankcombe et al. requires knowledge of the empty clathrate energy. The expressions of the energy of sII clathrate formation and the formation energy for sH clathrate formation consider only congurational, i.e., potential energy. Hence, the considerations are valid only at 0 K and 0 Pa. Moreover, the expressions do not take into account guest-guest interactions. We propose a new more general computational approach to distinguish between promoter molecules capable to form structure sI, sII hydrogen clathrate and those that can form structure H hydrogen clathrate. Our predictive method is based on Gibbs free energy calculations of sI, sII, and sH hydrogen clathrate and thus extends to nite temperatures and pressures.
In Section II we provide computational details. We then in Section III give the results and demonstrate that the van der Waals volume can be used to discriminate between promoters that form hydrogen sI, sII or sH hydrogen clathrate hydrates. We conclude with summary.

II. Computational method
The initial positions of oxygens in sI, sII and in sH are obtained from the experimental X-ray crystallographic clathrate structure as described in ref. [6][7][8][9]. The Monte Carlo algorithm developed by Buch 10 is employed to generate initial proton-disordered congurations according to the ice rules 11 and the conguration with dipole moment near zero is selected. The sI 2 Â 2 Â 2 supercell (24.06 Â 24.06 Â 24.06Å 3 initial dimensions) consists of 368 water molecules with 48 large and 16 small cages. The fully occupied sI hydrogen clathrate contains 432 molecules. The simulated cell of the sII clathrate hydrate is one unit cell (a ¼ 17.31Å) containing 136 water molecules in the host framework with 16 small cages and 8 large cages. The total number of molecules in the simulation box of the sII fully lled structure is 160. For the sH clathrate a 2 Â 2 Â 2 supercell (24.42 Â 24.42 Â 20.28Å 3 ) composed of 272 water molecules in the host lattice with 24 small, 16 medium, and 8 large cavities is used. The completely occupied sH clathrate consists of 320 molecules. Single occupancy of any cage of all types of clathrates is assumed. Double occupancy of a single cage is not considered. Complete occupancy of sI, sII, and sH clathrates means that the small and medium cages are occupied by a single H 2 molecule whereas the large cages are lled with a promoter molecule, i.e., there are no empty cages le. All molecules are treated as rigid.
Molecular interactions are pairwise additive and dened by the Lennard-Jones (LJ) potential for describing the van der Waals interactions. Coulombic long-range electrostatic interactions are treated using the Ewald summation method. 12 The intermolecular potential is given by where N is the number of LJ sites and M is the number of Coulomb charges. Long range corrections 12 for the truncation of the LJ potential are added to eqn (1). The cutoff of the Lennard-Jones potential and the cutoff for the real-space part of the Ewald sum equal half of the shortest simulation box side. Periodic boundary conditions are applied. The TIP4P/2005 (ref. 13) model is used for the water molecules in the host lattice. A linear rigid model with bond length of 0.7414Å is used for hydrogen molecule. 14 The geometries of all tested promoter molecules and partial charges on their atomic sites are taken from the ESI of ref. 5. The parameters of the Lennard-Jones potential for different promoter molecules are taken from Cornell et al., 15 For unlike site-site interactions, the Lorentz-Berthelot combining rule 3 ij ¼ ffiffiffiffiffiffiffiffi ffi 3 ii 3 jj p , s ij ¼ (s ii + s jj )/2 is employed. The Atomic and Bond Contributions method 16 was used to calculate van der Waals volume of molecules. Table 1 lists tested promoter molecules with corresponding abbreviations and known experimentally type of formed clathrate structure.
Consider the transformation of a fully occupied sII hydrogen clathrate hydrate with a certain promoter in the large cages and hydrogen in the small cages into a fully occupied sH hydrogen clathrate hydrate with the same promoter. The transformation reaction per water molecule is where S sII and S sH refer to the solid fully occupied sII and sH hydrogen clathrate, respectively, and H 2 and pr denote gaseous hydrogen and liquid promoter phase, respectively. Since the chemical potential is given by where f is a fugacity of a component, f 0 is a fugacity of the ideal gas reference state and m 0 ¼ 0 because in this study we are focused on the excess free energy calculations, the Gibbs free energy change per water molecule is where G sH and G sII denote the free energy per water molecule of solid fully occupied sH and sII hydrogen clathrate, respectively, and R is the gas constant, f pr and f H 2 are fugacities of liquid promoter phase and gas H 2 phase, respectively. The reaction proceeds to the right if DG form (sH, sII) < 0, to the le otherwise. Since bothS sII and S sH are solids the reaction proceeds completion at xed fugacities of the gaseous and the liquid components.
The fugacity of the promoter component f pr in the liquid phase is calculated using the expression where P* stands for the saturated vapour pressure calculated via the Antoine equation, 22 f* is the fugacity coefficient of the liquid phase component at saturation condition which is assumed to be equal to unity, x and V L are the mole fraction and molar volume of the promoter in liquid phase, and the activity coefficient g has been determined via the van Laar equation 22 with parameters calculated by tting innite dilution activity coefficient data 23 following the equation The fugacity of gaseous hydrogen phase has been obtained from the Virial equation with parameters taken from ref. 24: Calculations of the term (4) for all promoter compounds from Table 1 and found in ref. 23 demonstrate that the absolute value of the term is almost always less than 0.08 kcal mol À1 which is the error of the free energy difference G sH À G sII . Only for tetrahydrofuran (THF), the modulus of the term is greater than the error and equals to 0.12 kcal mol À1 . Therefore, it is a good approximation to neglect the term in eqn (4) and hence The expression takes into account all the host-guest and guest-guest interactions. Finite temperature and pressure are taken into consideration by the entropy included in the Gibbs free energy. The sign of the formation free energy in eqn (8) indicates the preference to the formation of either sII hydrogen clathrate or sH hydrogen clathrate. Bearing in mind that G sH and G sII are always negative, hydrogen structure H is energetically favoured to be formed when DG form is negative (G sH < G sII ). Similarly, when DG form (sH, sII) is positive (G sH > G sII ), then hydrogen sII hydrogen clathrate is preferably formed.
Similarly, the formation free energy discriminating between crystallization of sH and sI hydrogen clathrate is approximated by The formation free energy distinguishing between synthesis of sII and sI hydrogen clathrate is expressed as The self-referential (SR) technique 25 is employed for the free energy calculations. The method uses the principle that at xed volume and temperature the Helmholtz free energy is an extensive quantity and it doubles if the system size doubles. In the SR technique, the absolute Helmholtz free energy is calculated as the free energy difference between two systems with different number of molecules. The SR method consists of two stages: 'replication' and 'relaxation'. The replication stage calculates the free energy difference between the small singlesize system and a highly constrained double-size system produced by adding a replica to the original single-size system of interest. The double-size system has twice the Helmholtz free energy of the single-size system of interest. During the relaxation stage this self-similarity constraint applied at the replication stage is gradually relaxed until the molecules do not feel the constraint any more. The relaxation stage is governed by a constraint parameter a. This stage delivers the free energy difference between the fully constrained double-size system and fully relaxed double-size. The double-size system has twice the Helmholtz free energy of the single-size system of interest. Hence the Helmholtz free energy of the crystal of interest is the sum of the free energies differences obtained for the two stages. The Gibbs free energy is obtained via the Helmholtz free energy using the relation G ¼ F + PV.
Simulations in the isothermal isobaric ensemble (NPT) are performed to equilibrate a clathrate at given pressure and temperature. The Helmholtz free energy and congurational energy calculations of empty and fully lled sII and sH clathrate hydrates are conducted in the canonical (NVT) ensemble. Each Monte Carlo (MC) free energy run for a given value of constraint parameter a involves of 400 000-700 000 cycles for equilibration and 60 000-70 000 production cycles for obtaining ensemble averages. Translational and rotational moves are attempted at random with equal probability. Volume trial changes in the NPT simulations occur with frequency 1/N. More theoretical and computational details on the Helmholtz free energy calculations of molecular crystals are provided in ref. 25 and 26. Fig. 1 DG form (sH, sII) formation free energy for various tested promoter molecules at 1000 atm and 233 K, 253 K, and 273 K. The empty columns represent the formation free energy for enclathrated promoters known experimentally as sH clathrate formers. The solid filled columns refer to DG form (sH, sII) for entrapped guest molecule determined experimentally as sII (in some cases both sI and sII) clathrate promoter. The columns filled with wide downward diagonal lines are used for experimentally untested cyclononane (CNO), cyclodecane (CDE). Errors are to two standard deviations.

A. Distinguishing sII and sH clathrates formation
To validate the computational approach for differentiating between guest molecules that can form sH clathrates and molecules that can form sII clathrates, 37 experimentally known promoter molecules have been tested. Structure sII and sH are of special interest of this study because these clathrate types are prospective candidates for hydrogen storage material. The Gibbs free energies of sII and sH hydrogen clathrate hydrates were obtained from self-referential MC simulations at 1000 atm. All small cages of sII and sH clathrates as well as the medium cages of structure H are occupied by single H 2 . In all large cavities of sII and sH a promoter (experimentally known as either sII or sH former) is inserted to calculate the Gibbs free energies at 233 K, 253 K, and 273 K. The formation free energy DG form (sH, sII) is then calculated with eqn (8) for various guest molecules. DG form (sH, sII) is depicted in Fig. 1 in an ascending order at the three studied temperatures. Fig. 1 indicates that DG form (sH, sII) is always negative (i.e., G sH < G sII ) for promoters forming sH clathrate. However, it is worth noting that the formation free energy for sII hydrogen clathrate with cyclohexane (CHX) at 253 K and dimethylpropane (DMPH) and methylpropene (MPROPE) at 273 K is smaller than the error bar. The proposed computational approach for predicting the type of clathrate structure (either sII or sH) to be formed by a potential promoter molecule, demonstrates excellent agreement with experimental evidences.
The order of the formation free energies for different promoters seen in Fig. 1 cannot be altered with the change of pressure. The preference for formation of a particular type of clathrate (sII or sH) cannot be facilitated by pressure either. This is demonstrated at 233 K in Fig. 2

B. Dependence on volume
At all studied temperatures it has been observed that the van der Waals volume (V vwd ) of enclathrated promoter molecule is a very important parameter in screening potential molecules capable to form either sII or sH hydrogen clathrates. Fig. 3 demonstrates the formation free energy versus the van der Waals volume of known promoters forming sII and sH clathrates. Fig. 3 conrms computational study 5 and an experimental rule of thumb that at all studied temperatures there is a rough threshold of V vwd of the promoter which separates secondary guest molecules forming sII clathrates from those known as sH clathrate formers. This is found to be the case at all temperatures considered even though the relative stability of the hydrates of different promoters alters with temperature. Nonetheless, it is worth noting that as the van der Waals volumes of sH and sII formers overlap, there is no denite threshold which able to predict type of formed clathrate. The same conclusion has been reached in ref. 5.

C. Distinguishing sI and sII clathrates formation
To put this method further to the test we also tried to distinguish between formation of sI and sII clathrates. To determine whether sI type hydrogen clathrate can be formed by a specic promoter molecule, additional calculations of formation free energy DG form (sH, sI) and DG form (sII, sI) using eqn (9) and (10) are required. The clathrate with the lowest formation free energy would be the type predicted. Fig. 4 demonstrates DG form (sII, sI) for ten selected promoter molecules at three different temperatures. In comparison with DG form (sH, sII) and DG form (sH, sI), the formation free energies DG form (sII, sI) for these chosen ten promoter at studied temperatures have the lowest value and hence, determine whether either sII or sI type of clathrate would be crystallized. In Fig. 4 DG form (sII, sI) for cyclopentane (CP), propane (PROP), and sulfurhexauoride (SHF) demonstrates preference for sII type hydrogen clathrate formation at all temperatures, except at 233 K when it is positive, however, the formation free energy at that temperature is smaller than the error bar. DG form (sII, sI) is always positive for ethane (ETA), ethene (ETE), and uoromethane (FLU) that  conrms that sI hydrogen clathrates are synthesized with help of those promoters. Acetylene (ACET), cyclobutane (CBUT), methane (MET), and oxetane (OXE) compounds can form both sI and sII clathrates and DG form (sII, sI) can be both positive and negative to still agree with experiments. Fig. 4 shows very good agreement with experimental information on the clathrate type formed by a selected promoter. It is also encouraging to see that at the lowest temperature, DG form (sII, sI) of the clathrates which can form both sI and sII types mostly fall in between the values of sI and sII formers.
To the best our knowledge compounds cyclononane (CNO) and cyclodecane (CDE) have not been tested for synthesizing clathrate hydrates experimentally. Our calculations have shown that at all temperatures the formation free energy DG form (sH, sII) for these compounds is always negative (Fig. 1). We therefore propose CNO and CDE as potential candidates for sH hydrogen clathrates formers. In a previous computational study 27 we have shown that the stability of sH hydrogen clathrate hydrates with these two promoter molecules demonstrates an unusually strong temperature dependence, which can be exploited if easy release of hydrogen is required.

IV. Conclusions
We have presented and validated a computational method for predicting a type (either sI or sII or sH) of hydrogen clathrate which can be formed by a particular promoter. The method is based on the Gibbs free energy calculations and is valid at nite temperatures and pressures. It has been found that the van der Waals volume of a potential promoter is an important parameter which can determine a type (either sII or sH) of clathrate to be formed. It has been found that it is impossible to change a type of hydrogen clathrate structure by decreasing/increasing pressure imposed on the clathrate. At all studied temperatures there is a rough threshold of the van der Waals volume of entrapped promoter molecule which distinguishes promoter molecules capable to form sII clathrates from those forming sH clathrates. Our formation free energy calculations indicate that the promoters in Table 1 known as formers of pure sI or sII or sH clathrates formers can form corresponding hydrogen clathrates. However, this requires experimental verication.