13044
J. Phys. Chem. C 2007, 111, 13044-13052
Molecular Dynamics Simulations of Type-sII Hydrogen Clathrate Hydrate Close to Equilibrium Conditions Terry J. Frankcombe* and Geert-Jan Kroes† Leiden Institute of Chemistry, Gorlaeus Laboratories, Leiden UniVersity, P.O. Box 9502, 2300 RA Leiden, The Netherlands ReceiVed: February 5, 2007; In Final Form: July 5, 2007
Molecular dynamics simulations of hydrogen-containing type-sII clathrate hydrate slabs have been performed. A range of water-water interaction potentials were tested. Three potentials that yield water ice melting points close to 273 K gave very different clathrate stabilities. The TIP5P potential was found to give the closest match to experimentally observed clathrate stabilities. Using this potential, we investigated molecular hydrogen mobility within the clathrate. It was found that hydrogen molecules within the large clathrate cages were highly mobile, migrating rapidly from cage to cage and giving instantaneous cage occupations of up to six hydrogen molecules per cage. Self-diffusion coefficients for hydrogen migration were calculated and found to obey an Arrhenius relation. Hydrogen in the small cages was not mobile on the time scale of the simulations. Simulations including the known promoter molecule tetrahydrofuran demonstrated enhanced stability of the clathrate, consistent with experimental data.
1. Introduction Using hydrogen as an energy carrier is an attractive alternative to the current global dependence on hydrocarbon-based fuels, particularly for vehicular transport. Storing pure hydrogen on board passenger cars is unattractive from an engineering point of view, primarily because of the low density of hydrogen gas and the effort required to compress it to a reasonable volume or liquefy it for cryogenic storage at 20 K. Although catalyzed complex metal hydride systems remain one of the most popular potential high-density storage materials, hydrogen clathrate hydrates have recently emerged as an alternative. Crystalline inclusion compounds in which hydrogen molecules occupy the spaces inside well-ordered cages of water molecules, type-sII hydrogen-containing clathrates1 were first observed under high pressure.2 Although the hydrogen density within sII hydrogen clathrates is not yet firmly established (with single small-cage occupation demonstrated by neutron diffraction studies of D2-containing clathrate3 and double small-cage occupation suggested by Raman spectroscopy of H2-containing clathrate2), pure hydrogen-water clathrates yield material hydrogen densities of 3.9-5.0 wt % and 31-43 kg H2 m-3, comparable to the densities achieved in pressurized gas systems. At temperatures not far below 0 °C, pressures required to form the clathrate from a water and hydrogen mix are on the order of 1000 bar.4 The requirement of these high pressures is a significant impediment to using pure hydrogen clathrate hydrates in hydrogen storage systems. It has been shown recently that introducing small amounts of the “promoter” molecule tetrahydrofuran (THF, C4H8O) significantly reduces the pressures required to form a clathrate.5 The THF-H2-H2O clathrate forms at the accessible pressure of 50 bar at 280 K. This reduction of required pressure comes * Present address: Department of Chemistry, Physical Chemistry, Go¨teborg University, SE-412 96 Go¨teborg, Sweden. E-mail:
[email protected];
[email protected]. † E-mail:
[email protected].
at the cost of reduced hydrogen storage capacity, with potential hydrogen storage volume taken up by the promoter molecule. It has been proposed that the hydrogen storage capacity of THF-H2-H2O clathrate can be tuned by varying the THF concentration, with hydrogen densities over 4 wt % being measured using Raman spectroscopy.6 However, efforts to confirm these results with direct gas content measurements, NMR spectroscopy, and neutron diffractometry7,8 found a hydrogen content of at most 1 wt % in similar THF-H2-H2O clathrate samples, with all large clathrate cages being filled with THF. Nonetheless, the possibility remains that other promoter molecules may interact with water in such a way as to stabilize the clathrate structure without fully occupying the large cages of the sII clathrate, or allow an sH clathrate to form with a lower promoter-molecule-occupied cage fraction. Thus, investigating other promoter molecules could lead to binary clathrate systems with a more advantageous balance between clathrate formation conditions and hydrogen storage capacity. Given the large range of potential promoter compounds, insight gained through computational modeling has the potential to significantly speed the development of promoted hydrogencontaining clathrate hydrates. Of principle interest are the characteristics of a good promoter molecule, in terms of both its interactions with water and hydrogen and the mechanics of the processes that lead to promoted clathrate formation and stabilization. Hydrogen clathrates have been modeled recently by Alavi et al. using molecular dynamics (MD) methods. These authors have published studies of the effect of multiple cage occupation9 and promotion with THF10 in type-sII clathrates, and type-sH clathrates promoted with 2-methoxy-2-methyl-propane (methyltert-butylether).11 These investigations were focused primarily on the effect of multiple H2 occupation on the calculated total potential energy, unit cell volume, and thermal expansion coefficients of the clathrate. Model potentials were used to describe the interactions between the molecular components of the clathrate. Generic model potentials were used for the
10.1021/jp071006e CCC: $37.00 © 2007 American Chemical Society Published on Web 08/16/2007
Simulations of Type-sII Hydrogen Clathrate Hydrate interactions with the promoter molecules. The SPC/E model12 was used for the water-water interaction, which has known deficiencies in describing liquid water and water ice.13-15 There was no report of validation of these potentials for describing pure or promoted hydrogen clathrate hydrates or investigation of the sensitivity of the results to the nature of the potentials. Thus, one of the main focuses of this work was to test model interaction potentials in order to identify potentials that do a good job at reproducing experimentally observed properties of hydrogen-containing clathrates. Direct microscopic modeling of the stabilizing effect of THF in hydrogen-containing clathrate hydrates have not been published previously. (Though Alavi et al.10 did examine the energetic and volumetric effects of varying the THF fraction, the stability of the clathrate was not directly assessed in that work.) Thus, we have performed preliminary modeling demonstrating the stabilizing effect. Additionally, we have explored hydrogen diffusion through the clathrate structure. This paper is structured as follows. The following section describes some of the computational methods we have used. Section 3 describes calculations designed to test and validate the intermolecular potential functions used in this work. The diffusion of hydrogen through the clathrate is studied in Section 4, followed by some observations on the mechanism of the collapse of the clathrate structure in Section 5. Modeling including THF as a promoter is described in Section 6, before a concluding section. 2. Computational Methods All MD calculations were performed using the DLpoly2 package.16 3D periodic NPT simulations were performed using the anisotropic, stress-relaxing Berendsen thermostat and barostat,17 with time constants of 0.5 ps for both. A time step of 0.001 ps was used throughout. Except where noted, the system being modeled was an sII clathrate slab immersed in pressurized hydrogen gas, with the intent of mimicking the conditions in the experiments of Lokshin and Zhao.4 The slab was constructed with a 2 × 2 surface cell and was two unit cells thick (thus, the simulation cell contained a 2 × 2 × 2 replication of the conventional Fd3hm unit cell of the ideal sII clathrate). The water molecules were arranged with proton disorder that minimized the slab net dipole moment while obeying the so-called ice rules.18 The small (512) clathrate cavities were initially populated with single randomly oriented H2 molecules, whereas the large (51264) cavities were each populated with randomly oriented clusters of four H2 molecules, in accord with occupations measured directly by neutron diffractometry.3 The surfaces of the slab were constructed by simply cleaving at the surfaces of the Fd3hm unit cell, with water molecules left intact. Incomplete cages on this {100} surface were left unmodified. The 8 complete unit cells comprised 1088 water molecules. The lattice parameter in the z (surface normal) direction was approximately 68 Å (twice the slab thickness, giving initial c/a and c/b ratios of 2). The space between slab images was filled with H2 molecules sampled from a simulation of hydrogen gas at the relevant temperature and pressure. The force field used was a sum of Coulombic and other nonbonded terms. The Coulomb interaction was applied between a set of partial charges located on a number of sites of the rigid water and hydrogen molecules. A consistent set of charges was used throughout the simulation, with the contribution of the Coulomb terms to the total energy and forces calculated by Ewald summation.19 The remaining interaction was modeled by a sum of radial functions, described below.
J. Phys. Chem. C, Vol. 111, No. 35, 2007 13045 Seven different model potentials (and corresponding H2O geometries) were used to describe the water-water interaction, namely, the SPC/E,12 TIP3P,20 TIP4P,20 TIP4P/ice,21 TIP4P/ 2005,22 and TIP5P23 models and the six-site model of Nada and van der Eerden24 (hereafter referred to as “6 site”). These popular models specify Lennard-Jones interactions between the oxygen atoms of the water molecules for the non-Coulombic van der Waals interaction, as well as specifying different sets of point charges located at various points around the rigid water molecule. The potentials for interactions involving H2 were similar, with charges to reproduce the quadrupole moment of H2 being placed on the atom sites and at the center of the molecule.25 These charges were not altered between the simulations. Thus, while the H2 -H2 interaction remained constant, the electrostatic part of the water-H2 potential varied as the charges on the water sites were varied. The isotropic potential of Buck et al.26 was applied between the H2 molecules. Lennard-Jones-like 8-6 potentials were applied between the center point of the H2 molecules and the water oxygen and hydrogen atom sites and a site midway between the two water H atoms.25 Although morerecent water-H227,28 and H2-H229-34 potentials have been published, these are not in the radial site-site form preferred for MD simulations, usually being expressed as spherical harmonic expansions. The H2-H2 potential adapted30 from that of Silvera and Goldman35 that is popular for solid hydrogen simulations could have been used in place of the Buck potential. No corrections for energetic and entropic quantum effects (such as the quantisation of molecular rotations) have been applied. The work of Wang et al.36 indicates that at the temperatures considered in this work quantum effects in gasphase hydrogen are sufficiently small to be ignored. Although quantum effects in water ice are significant,37,38 no procedure to correct MD simulations for these effects have been established. The truncated parabola tunneling correction of Fermann and Auerbach39,40 was estimated on approximate hydrogen migration paths in the context of discussing diffusion rates. This is described in full in Section 4. 3. Assessment of Water-Water Potentials Simulations of the clathrate slabs were performed in order to assess the suitability of each of the considered water-water potential models. These tests were based on reproducing the observed phase behavior of the hydrogen and water system under excess hydrogen conditions. A pressure of 1500 atm (152 MPa) was selected. The work of Lokshin and Zhao4 shows that under hydrogen gas clathrates should be stable at 265 K and decompose at 285 K. Reproducing this behavior for the modeled clathrate slab was chosen as the criteria against which the model potentials were assessed. A method reminiscent of the “direct coexistence” method of determining melting temperatures41-43 was used, in which MD simulations were performed at a number of temperature and pressure conditions, starting with a solid clathrate in contact with hydrogen gas. The results of the simulations were analyzed to determine whether or not the clathrate decomposed. At the temperatures of these simulations, the instantaneous positions of the 1088 water molecules did not clearly reveal the clathrate structure upon casual inspection. Thus, a number of measures were developed to characterize the clathrate structure and assess whether the clathrate remained stable or decomposed. No measures directly based on calculating radial distribution functions were found to be successful at distinguishing clathrate structures from nonclathrate structures. In addition,
13046 J. Phys. Chem. C, Vol. 111, No. 35, 2007
Frankcombe and Kroes
Figure 2. Average coordination number of water around isolated, internal H2 molecules. 265 K, 1500 atm. Clathrate (20), liquid water (16), and threshold (18.5) values are indicated.
Figure 1. z distribution of hydrogen through the clathrate slab as a function of time. Yellow indicates high H2 density, moving through red, purple, to black as the H2 density decreases. (a) TIP4P potential, 265 K and (b) 6-site potential, 285 K, both at 1500 atm pressure.
no Fourier-based methods have proven to be successful in studies of hexagonal water ices.44,45 One rather qualitative but successful clathrate structure measure was the z distribution of molecular hydrogen through the simulation as a function of time. To calculate this measure, the simulation cell was divided into a number of thin slices parallel to the (001) face of the simulation cell (thus parallel to the nominal slab surface). For each time step, hydrogen molecules within each slice were counted to build up the distribution. The distribution as a function of time can be plotted as a color intensity plot. Examples are shown in Figure 1 for the TIP4P potential at 265 K and for the 6-site potential at 285 K. In this plot, z coordinates corresponding to high hydrogen density are shown in yellow, with red, purple, blue, and black indicating successively lower hydrogen density. Note that only the central 56 Å of the approximately 68 Å z coordinate of the simulation cell is shown. In the first part of the representation in Figure 1a, corresponding to times less than about 50 ps, the hydrogen clearly remains localized to particular levels within the slab, suggesting that a well-ordered clathrate structure exists. At times beyond 50 ps, the well-ordered hydrogen structure breaks up, appearing to some extent to exclude the hydrogen from the central, watercontaining slab region. One can conclude from this that the clathrate slab is not stable beyond the first 50 ps at 265 K with the TIP4P water-water interaction. Alternatively, the portion of the hydrogen density evolution shown in Figure 1b (covering 1 ns from 4 ns of simulation time) indicates that with the 6-site potential the clathrate remained stable for at least 5 ns at 285 K. Examination of these hydrogen density profiles indicated that the clathrate slabs described using the SPC/E, TIP3P, and TIP4P potentials were extremely unstable at 265 K. Slabs simulated using the SPC/E and TIP4P potentials lost their clathrate character after around 50 ps. The TIP3P potential resulted in even less-stable slabs, losing hydrogen z structure in less than 10 ps. The tested modifications of the TIP4P potential, TIP4P/ 2005 and TIP4P/ice, yielded slabs that decomposed in 250 and 850 ps, respectively. Only the slab described by the 6-site potential showed no sign of decomposing in 5 ns, with the TIP5P
potential giving a slab that decomposed after 1.4 ns. However, 5-ns simulations at 285 K with the 6-site potential showed an essentially identical evolution of the hydrogen z density, indicating that the modeled clathrate was stable at this temperature. 285 K is above the experimentally observed decomposition temperature of the clathrate, which lies below 280 K at this pressure. At 285 K, the TIP5P slab decomposed in 270 ps. Several other features of the hydrogen distribution are evident in these plots. The first of these is that the outlying peaks of the initial hydrogen distribution (the peak closest to the slab surface on each side, located at around z ) (15 Å at zero time in Figure 1a) disperse in the first few picoseconds of the simulation. This can be interpreted as the water cages making up the clathrate cavities just below the surface being disrupted as the surface adjusts to the crude cleavage applied when constructing the slab, allowing the contained hydrogen to escape. Second, from a few picoseconds into the simulation there appears to be an accumulation of hydrogen adhered to the surface of the clathrate, around z ) (20 Å. This surface layer of higher hydrogen concentration was on the order of 5 Å thick. A third feature is a region at the surface of the slab with low hydrogen density. This region lies around z ) (15 Å and forms once the initial outlying hydrogen peaks have dispersed very early in the simulations. The hydrogen has been excluded from a water layer between the well-ordered clathrate of the slab and the hydrogen interface layer. These features were general, being observed in every simulation. A second measure of the clathrate structure used was the coordination number of water around hydrogen molecules. The 512 cages of sII clathrate comprise 20 water molecules at a radius of around 4 Å. Recent MD and Monte Carlo simulations46 indicate that the coordination of water molecules around H2 suspended in bulk water is 16. Thus, assuming single occupation of H2 molecules in the small cages of the sII clathrate, the number of water molecules within 5.0 Å of isolated H2 molecules provides a discriminant between the clathrate structure and hydrogen-containing water. It is not clear what the value of this parameter would be for hydrogen somehow contained within hexagonal ice. Figure 2 shows the calculated average coordination number as a function of time for the simulations at 265 K using the seven different water-water potentials. For each time step, the coordination number around a hydrogen molecule was included in the calculation of the average coordination number if the center of mass of the hydrogen molecule was at least 4 Å from
Simulations of Type-sII Hydrogen Clathrate Hydrate
Figure 3. Change in the total energy of the simulation cell under constant temperature conditions. 265 K, unless indicated otherwise, and 1500 atm.
all other hydrogen molecules and within 14 Å of the original middle of the clathrate slab. Because the coordination number data exhibited a high degree of noise (typically (0.4 in regions where it had fallen significantly below the ideal clathrate value of 20), Gaussian smoothing with a width parameter of 6 ps has been applied to the data shown in Figure 2. Note that simulations were not continued if it was obvious that the clathrate slab had collapsed. For all potentials except the TIP5P and 6-site potentials, the calculated average coordination number started close to the ideal clathrate value of 20 and dropped toward the water value of 16. With the TIP5P potential, the coordination number stabilized around 18. The evolution of the calculated coordination number gave a more well-defined measure of the destruction of the clathrate slab than the z distribution of hydrogen discussed above. Considering a threshold value of 18.5 for the coordination number indicating collapse of the clathrate, the two measures gave remarkably similar lifetimes for the clathrate structures. Only using the TIP5P and 6-site potentials yielded clathrate slabs that remained stable for more than 1 ns. For the TIP5P potential, the slab remained intact for almost 1.5 ns at 265 K and decomposed after just 170 ps at 285 K. The coordination number confirmed that the clathrate slab modeled using the 6-site potential was very stable. At 285 K, the behavior of the average coordination number was essentially identical to that shown in Figure 2 for 265 K, oscillating around 19.85 (around 0.05 lower than the corresponding value at 265 K) with no sign of decrease after 5 ns of simulation time. A third measure of the clathrate structure, the mean squared deviation of the oxygen-oxygen distances within the known hydrogen-bonding pattern of the sII clathrate, proved useful. This measure shall be discussed further in Section 5. From the point of view of the stability of the clathrate slabs modeled with the different water-water potentials, this mean squared deviation measure correlated well with the hydrogen z distribution and average coordination number measures discussed above. The fourth measure of the clathrate structure considered was the total energy (potential plus kinetic) of the unit cell. This is the measure commonly used in the direct coexistence method of determining melting points from MD simulations and relies on the fact that a melting structure converts kinetic energy into potential energy. The simulation thermostat then effectively adds kinetic energy to the system in the course of maintaining the temperature, leading to an increase in the total energy for a melting system.43 Figure 3 shows the change in total energy per unit cell for each of the simulations at 265 K and 1500
J. Phys. Chem. C, Vol. 111, No. 35, 2007 13047 atm. With this measure, the slightly higher stability of the TIP5P clathrate indicated by the other measures was not obvious, with the total energy increasing after initiation at a rate similar to that exhibited by the TIP4P/ice potential. Figure 3 shows the total energy of the 6-site potential clathrate in simulations at 265 and 285 K. Even at the higher temperature where the clathrate should decompose, the total energy of the 6-site potential simulation showed a slight tendency to decrease, indicating a stable solid phase.43 The preceding results clearly show that the SPC/E, TIP3P, TIP4P, TIP4P/2005, and TIP4P/ice potentials yield clathrates that are significantly less stable than the real clathrate, decomposing quickly under conditions at which the real clathrate is stable. The TIP5P potential also appears to yield unstable clathrates at 265 K, but these are more stable than those for the other potentials mentioned above. Furthermore, 5-ns simulations of the clathrate slab at 245 K using the TIP5P potential showed the clathrate to be stable, indicating that the decomposition temperature of the clathrate modeled with TIP5P lies in the range 245-265 K. Of course, one must bear in mind that it is known that in general thin slabs do melt at slightly lower temperatures than the corresponding bulk material,47,48 though the effect is difficult to quantify using the MD approach used in this work. The 6-site model appears to yield too stable clathrates, with the modeled clathrate slab showing no sign of decomposing after 5 ns at 285 K, above the clathrate stability region identified experimentally.4 In principle, the MD calculations performed in this work can only show potentials to yield insufficiently stable clathrate structures. Starting from a clathrate structure, when a clathrate has decomposed under conditions at which it should be stable one can definitely conclude that the potential is no good for modeling the clathrate. However, with MD simulations (as opposed to Monte Carlo simulations) one cannot conclude that clathrates that do not decompose in any given time frame (accessible to MD simulations) are stable. They may decompose if given sufficient time. This is particularly the case for H2water clathrate systems, which exhibit significant freezingmelting hysteresis.49 The fact that the “melting line” of the phase diagram presented by Lokshin and Zhao4 is based on clathrate formation data strengthens the conclusion that the clathrates described by potentials other than TIP5P and the 6-site potential are far too unstable: in reality, at 265 K the clathrate forms from a H2-water mix despite the hysteresis, so they should definitely not decompose if initialized as a clathrate in the simulations. The stability ordering implied by Figures 2 and 3 (TIP3P is the least stable, followed by SPC/E, TIP4P, and TIP4P/2005) correlates well with the ice-melting temperatures of these models, from 145 K for TIP3P to 250 K for TIP4P/ 2005.22,24,43,50,51 This correlation does not hold for the remaining models, with TIP4P/ice, TIP5P and the 6-site model all giving ice-melting temperatures within a few degrees of the 273 K melting point of real water,22,24,43,50,51 yet demonstrating widely different clathrate stabilities. Presumably, this is due to differences in the medium-range behavior of the potentials. Although each water molecule participates in four hydrogen bonds in both the ice Ih and clathrate structures and thus feels a similar immediate environment, the voids in the clathrate structure mean that water molecules in clathrate interact with a lower number of water molecules at medium distances than those in water ice, arranged with a different distribution of orientations. Differences in the potentials in these regions are thus reflected in the different behaviors for water ice and clathrate hydrate.
13048 J. Phys. Chem. C, Vol. 111, No. 35, 2007
Frankcombe and Kroes
Although one cannot rule out the possibility that the 6-site potential yields a clathrate that would decompose above the decomposition temperature of the real clathrate given sufficient time, the present indications are that it is far too stable. Thus, we conclude that the TIP5P potential, though not perfect and yielding a decomposition temperature that is too low for the clathrate at 1500 atm, is the most suitable for clathrate simulations, of those tested. The fact that this potential yields clathrates that decompose slightly too fast may be considered an advantage when studying the effect of promoter molecules on clathrates close to equilibrium conditions and the mechanisms of clathrate decomposition. 4. Hydrogen Diffusion MD simulations were performed in order to study the diffusion of molecular hydrogen within the clathrate structure. The methodology was similar to that described above, except simulations were performed using a single unit cell with 3D periodic boundary conditions but without gas-filled gaps in the simulation cell, thus simulating the bulk material. The TIP5P water-water potential was used. Simulations were performed with an external pressure of 1500 atm at temperatures of 200, 245, 255, and 265 K. Two types of simulations were performed, with the clathrate cavities either initially fully occupied with hydrogen according to the neutron diffraction occupations3 (four hydrogens in large cages and one in small cages) or with full occupation except a single large cage completely empty of hydrogen. For each temperature and initial occupation setup, a number of simulations were performed, for up to 4 ns of simulation time. Every 1 ps of simulation time, the hydrogen occupation of the clathrate cavities was analyzed. From this analysis, hydrogen jump eventsswhen a hydrogen molecule moved from one cavity to anotherswere identified. Hydrogen diffusion was observed to be sufficiently slow that the 1-ps time resolution of this analysis was sufficient. A number of unexpected results arose from this molecular diffusion analysis. In all of the simulations performed, not a single change in the occupation of the singly occupied small cages was observed. At the temperatures being modeled, with the current set of potentials, the potential barrier presented to a hydrogen molecule traversing through one of the pentagonal faces of the 512 water cage was insurmountable in the time scale of the simulations. Although the H2-H2 interaction energy could contribute to the effective barrier in the case of attempting to doubly occupy the small cavities,9 no such energy penalty need be paid to move a hydrogen molecule from a small cage into a large cage of low occupation, implying that it was indeed the H2-(H2O)5 energy involved in moving hydrogen through the pentagonal faces that prevented altering the small-cage occupation. Among the large cages, the hydrogen was much-more mobile than expected, even in the “fully occupied” case. Cage occupation was highly variable. Occupations of six hydrogen molecules in a large cage were observed frequently, with lifetimes measured in hundreds of picoseconds. The Einstein relation was used to determine the self-diffusion coefficient of hydrogen moving within the clathrate. Treating the clathrate as a network of sites to be occupied, the Einstein relation can be written as
D ) gl2Γ
(1)
TABLE 1: Average Lattice Parameter, Hydrogen Cage Jump Frequency, and Resultant Hydrogen Diffusion Coefficient for Simulations of Fully Occupied sII Clathrates temp (K)
〈a〉 (Å)
Γ (ns-1 cage-1)
D (× 10-11 m2s-1)
200 245 255 265
16.78 16.88 16.90 16.92
0.05 ( 0.08 1.1 ( 0.3 1.8 ( 1.4 3.6 ( 0.7
0.7 ( 1.1 15 ( 4 24 ( 18 49 ( 10
in which l is the site to site migration distance, Γ is the hopping frequency, and g is the geometry factor.52 Although the geometry factor for bulk diffusion is usually 1/6, the tetrahedral coordination of large cages in the sII clathrate (sharing their hexagonal faces) leads to a geometry factor of 1/4 in this case. Treating the migration distance as the distance between the centers of the large cages (which are located at the 8b positions of the Fd3hm space group) gives, on average
l)
x3a 4
(2)
where a is the lattice constant. The lattice constant was taken as the time-average over the simulation. The hopping frequency was determined by accumulating the inverses of the times between hopping events (cage occupation changes). After dividing by eight (the number of large cages in the unit cell), the mean of these frequencies were calculated for each temperature, yielding the hopping frequencies shown in Table 1. The uncertainties quoted are plus/minus one standard deviation of the observed inverse periods. In the simulations at 245 K or above, no induction period was evident. The self-diffusion coefficients could then be simply calculated via eq 1. Clearly, the hydrogen migration at 200 K was so slow that it made the statistics attained in the 1-ns simulations performed insufficient for an accurate determination of the hopping frequency. A larger number of simulations were performed for 245 and 265 K, leading to the more-accurate determination of the hopping frequencies in these cases than those at 255 K. Indeed, the results from the simulations at 200 and 255 K have been included for completeness only. Considering only the more precise 245 and 265 K mean hopping frequency data, the diffusion coefficients can be fit to an Arrhenius form
D ) D0 exp
(-Q kT )
(3)
where k is the Boltzmann constant and T is the temperature, with an activation energy Q ) 32.3 kJ/mol and preexponential factor D0 ) 1.15 × 10-3 m2 s-1. Compared to the quoted uncertainty, this fit passes very close to the calculated 255 K value, giving D ) 27.8 × 10-11 m 2 s-1. Considering the one standard deviation uncertainty in the calculated diffusion coefficients introduces significant uncertainty into these values, to Q ) 32 ( 12 kJ/mol and D0 ) 10-3(2 m2 s-1. In the simulations initially with an empty large cage, there appeared to be a tendency toward slower hydrogen migration. The calculated diffusion coefficients are shown in Table 2. No obvious tendency toward initially filling the empty cage was observed. Indeed, no spatial differences in the rate of cage to cage jumping could be discerned. The average lattice parameter of the simulations with lower hydrogen content was slightly smaller than that in the corresponding fully occupied simulation, though the average potential energy per unit cell was slightly higher. This is consistent with previous simulations.9
Simulations of Type-sII Hydrogen Clathrate Hydrate TABLE 2: Average Lattice Parameter, Hydrogen Cage Jump Frequency, and Resultant Hydrogen Diffusion Coefficient for Simulations of sII Clathrates Initially with an Empty Large Cage temp (K)
〈a〉 (Å)
Γ (ns-1cage-1)
D (× 10-11 m2s-1)
245 265
16.86 16.90
0.9 ( 0.5 1.7 ( 1.0
12 ( 7 21 ( 13
An explicit calculation of the energy barriers along the hydrogen migration paths was performed in a static approximation. The potential energy of a clathrate hydrate with a single hydrogen guest molecule was calculated, with the hydrogen guest molecule progressively moved along a straight line from the center of a large cage to the center of an adjacent cage. The clathrate was held fixed at the ideal structure with lattice parameter 16.9 Å. That is, no relaxation was allowed as the hydrogen molecule passed through the cage face. Two paths were considered, one to an adjacent large cage (through a hexagonal face) and one to an adjacent small cage (through a pentagonal face). Calculations were performed with the hydrogen molecule oriented with its molecular axis pointing along the migration path and with the hydrogen molecule oriented with its molecular axis in an arbitrary direction perpendicular to the migration path (lying in a plane roughly parallel to the plane of the face through which it was passing). The calculated energies are shown in Figure 4. These calculated energies are consistent with the observed diffusion behavior. The barriers to migration through a pentagonal face (to or from a small cage) were found to be more than 7 times higher than those for migration between large cages through a hexagonal face. The calculated barriers for migration through hexagonal faces (27.6 kJ/mol and 23.8 kJ/mol) were remarkably similar to the 32 ( 12 kJ/mol activation energy calculated for the temperature dependence of the diffusion coefficient. If one takes the calculated barrier for migration through the pentagonal faces as the Arrhenius activation energy and assume the same D0 value as for migration through the hexagonal faces, then eq 3 yields a vanishingly small diffusion coefficient on the order of 10-42 m2 s-1 at 265 K. Clearly, any reduction in the barrier height caused by relaxation and thermal distortion of the clathrate water structure as the hydrogen passed through the pentagonal face would greatly increase this value, in accord with the exponential
Figure 4. Potential energy change as a single hydrogen molecule moves through a fixed clathrate structure, in terms of a normalized migration coordinate that takes the value 0 at the center of the originating large cage and 1 at the center of the adjacent large or small cage. The axis of the hydrogen molecule is aligned with (solid) and perpendicular to (dashed) the migration path.
J. Phys. Chem. C, Vol. 111, No. 35, 2007 13049 dependence of the diffusion coefficient on the negative of the barrier height. The calculated energies along static migration paths could be used to calculate an estimate of the magnitude of hydrogen tunneling. Using the truncated parabola model of Fermann and Auerbach,39,40 the increase in the rate of migration due to tunneling was estimated as
[ (
ktunnelling 1 ∆E - θ0 + ) exp kclassical 2 kT
)
∫-∞θ exp 0
( )
]
∆Eθ sech2θ dθ kTθ0 (4)
where ∆E is the calculated barrier height and
θ0 )
π∆E p|ω|
(5)
with p being the Planck constant divided by 2π and ω being the imaginary vibrational frequency corresponding to motion over the barrier. Fitting a quadratic to the maximum region of the lowest calculated static migration barrier gave a harmonic frequency of 2700i cm-1. Equation 4 then yields tunneling factors decreasing from 1.52 at 200 K to 1.26 at 265 K. These factors act as crude estimates of the effect of tunneling on the calculated diffusion coefficients listed in Table 1. Thus, including this estimate of tunneling effects, the diffusion coefficient for hydrogen migration at 265 K would be (62 ( 10) × 10-11 m2 s-1. For migration through the pentagonal faces of the water cases, the calculated harmonic frequency of 5000i cm-1 gave much-larger estimates of tunneling effects, but the higher calculated migration barrier led to a stronger temperature effect. The calculated tunneling factors were 10.2, 3.2, and 2.6 at 200, 245, and 265 K, respectively. However, the extremely small diffusion coefficient due to the high barrier (see above) means the tunneling effect would need to be many orders of magnitude larger to explain significant diffusion to and from the small clathrate cages. The high hydrogen mobility observed in the 3D periodic calculations described in this section was also observed in the slab calculations of the previous section. In all cases, the hydrogen content within the confines of the clathrate slab was seen to decrease with time, with the hydrogen migrating into the gaseous region. This is consistent with the nonideal occupations at higher temperatures measured by Lokshin et al.3 Counting occupations within spheres of radius 3.9 Å centered on the ideal cage positions, Figure 5 shows the total number of molecules contained in large cages. In Figure 5, these are divided into molecules contained in the inner five, three, and one layers of large cages within the slab (of the seven layers of complete large cages in the ideal clathrate slab). Contrary to expectations, the 265 K slab modeled with the 6-site potential exhibited a slower migration of the hydrogen out of the slab than a slab modeled with the TIP5P potential at the lower temperature of 245 K. The apparent rigidity of the slab described by the 6-site potential slowed the hydrogen migration through the oscillating hexagonal faces of the clathrate cages. Although the occupations in the 265 K, 6-site potential slab stabilized at values consistent with three molecules per cage between 4 and 5 ns of simulation time, occupations within the 245 K, TIP5P potential slab reached this level after around 1.5 ns and continued to fall, stabilizing at lower values after around 3 ns of simulation time. 5. Mechanism of Degradation In simulations of unstable clathrate slabs, the slabs decomposed from the surface inward, cage by cage. The surface-first
13050 J. Phys. Chem. C, Vol. 111, No. 35, 2007
Figure 5. Total number of H2 molecules contained in large clathrate cages, split into slab regions. Data from the occupation of the central layer of eight cages is shown in red, that from the inner three layers is shown in green, and that from the inner five layers is shown in blue. TIP5P potential at 245 K solid lines, 6-site potential at 265 K dashed lines. 1500 atm pressure. The expected values for average occupation 4 and 3 are indicated as thin dashed lines.
Frankcombe and Kroes In Section 3 above, it was suggested that the water molecules from cleaved cages form a liquid-like layer on the surface of the clathrate slab. Even for stable clathrate slabs, this is confirmed by growing D values when the outer 2-3 Å of the slab is included, indicating mobile water molecules. Both hightemperature slabs at 285 K stable by virtue of being described by the 6-site potential and TIP5P slabs at 245 K were found to exhibit this mobile water layer. Despite the lower temperature, the latter was found to have a much-greater degree of mobility, in terms of both the speed of the increase in D and the rate at which the mobility dropped off moving into the slab. Examination of the structure during the decomposition process indicated that the clathrate structure was not lost uniformly in the final, first-order decomposition of the central regions of the slab. Individual water cages, encapsulating some hydrogen molecules, remained intact after other water cages around them had lost their structure. These cages were not stable and decomposed in due course. They were not sufficiently stable to migrate from the region of their original location within the structure before collapsing. In contrast, during the initial surface disordering no preference for delocalising the water molecules belonging to particular clathrate cages could be discerned. In particular, no initial surface pitting of the solid clathrate surface could be identified under the mobile water layer, even on length scales consistent with individual near-surface cages collapsing. 6. Effect of THF
Figure 6. Average squared hydrogen-bond displacement for TIP5P clathrate slab (1500 atm, 265 K), considering only molecules within the central 8, 16, 24, and 28 Å of the slab
nature of the decomposition can be demonstrated easily by examining the changes in the oxygen-oxygen distances between initially hydrogen-bonded water molecules in the clathrate. Taking the mean square displacement of these distances from their ideal values gives a direct geometric measure of the clathrate structure, here called D. This geometric measure can be restricted to water molecules initially lying within certain regions of the clathrate, such as a central slice of a given thickness. Examples of the clathrate slab described with the TIP5P potential at 265 K are shown in Figure 6. Thermal motion means that even for stable clathrates D fluctuates around a small nonzero value, D0. The fact that the clathrate slab decomposes from the surface inward is clearly indicated by the fact that the onset of growth in D from D0 occurs at successively later times for successively thinner slices (which incorporate only regions successively further from the original surface). This trend only applies for the outer regions of the slab, however. The innermost regions of the slab disorder at approximately the same time. This is consistent with previous work on surface-induced melting which shows that slabs tend to initially undergo a continuous transition to quasi-liquid near the surface, before the entire slab melts in a first-order transition.48
Preliminary calculations have been performed to directly model the effect of THF on the stability of the hydrogencontaining clathrates. The methodology used was similar to that reported in Section 3 above. Using the TIP5P water-water interaction, simulations were performed at temperatures of 265-305 K under a pressure of 500 atm (51 MPa) for clathrate slabs containing only hydrogen guests and with the hydrogen contained in the large clathrate cages replaced by THF molecules, with one THF molecule per cage. For the THF-H2O and THF-H2 interactions, the model potentials used previously by Alavi et al.10 were used. At 265 K, 500 atm is outside the region of stability of the hydrogen-only clathrate.4 However, under these conditions the stable state for the system is hydrogen gas and water ice. The solid-phase clathrate to water ice transition occurs much-more slowly than the transition from the clathrate to liquid water and hydrogen. Thus, in the simulations of the clathrate at 265 K and 500 atm the clathrate slab was not observed to decompose in the timescales accessible to the MD simulations. The results are summarized in Figure 7, which shows the coordination number of water around isolated hydrogen within the originally clathrate region. Like in Figure 2, Gaussian smoothing has been applied to the data, in this case with a larger width parameter of 16 Å to counteract the greater rapid variations in the average coordination number observed in the THF-containing simulations. In both the hydrogen-only and hydrogen with THF clathrate cases, the qualitative behavior was similar to the higher pressure simulations of section 3, indicated in Figure 2. At lower temperatures, the average coordination number remained stable at values of 19.5 or greater. As the temperature increased, thermal motion distorted the cage structure more, reducing the average coordination to 19 or less. At still higher temperatures, the clathrate structure collapsed, indicated by the rapid drop of the average coordination number toward the liquid water value of 16. The presence of THF molecules in the liquid state resulted in a lower average coordination number of less than 14 because the THF molecules
Simulations of Type-sII Hydrogen Clathrate Hydrate
Figure 7. Average coordination number of water around isolated, internal H2 molecules for simulations of H2 (dashed) and H2 +THF (solid) clathrate slabs. 500 atm, 265-305 K. Clathrate (20) and liquid water (16) values are indicated.
reduced the water density and displaced water molecules from the solvation shells around the hydrogen molecules. Including THF in the simulations resulted in an increase in the temperature required to cause the collapse of the clathrate structure within the 4 ns of the simulations by around 20 K. Although the pure-hydrogen clathrate decomposed in around 1 ns at 285 K, similar simulations of THF-containing clathrate maintained an average coordination number around 19.5, indicating that the clathrate structure was maintained. At 295 K, the THF-containing clathrate exhibited lower average water coordination around hydrogen of around 19, consistent with greater thermal effects, but it did not decompose in 4 ns of simulation time. Only in the simulation at 305 K did the THFcontaining clathrate decompose on a time scale similar to that of the hydrogen-only clathrate at 285 K (around 1 ns). 7. Summary and Conclusions In this work, we have performed molecular dynamics simulations of hydrogen clathrate hydrates, mostly as slabs under hydrogen gas. A range of water-water interaction potentials were tested. Four different measures of clathrate structure were developed (based on bulk-structural, geometric, and energetic considerations), which showed a high degree of correlation. Potentials that give ice melting temperatures that are too low also yield clathrates that decompose at temperatures that are too low. However, three potentials that give reasonable and very similar ice melting temperatures (the TIP4P/ice, TIP5P, and 6-site potentials) yield clathrates of very-different stabilities. Although no water-water interaction was found to reproduce the behavior expected at the tested temperature and pressure conditions, the TIP5P potential seemed to give the closest match to the stability behavior expected given the published hydrogen clathrate phase diagram. Thus, we recommend that the TIP5P water-water potential should be used for future microscopic modeling of hydrogen clathrates. Of course, water-water interactions alone cannot describe a clathrate structure because the water structure is not stable at anything but the lowest temperatures without the presence of a guest gas filling the clathrate cages. Calculations modeling hydrogen-containing clathrate structures show that the waterH2 potential strongly influences properties of the clathrate, whereas the H2-H2 interactions are of secondary importance.53 Because different water-water potentials usually imply different partial charges on the water sites, they also alter the water-H2
J. Phys. Chem. C, Vol. 111, No. 35, 2007 13051 potential through the electrostatic component. Therefore, in investigating the dependence of clathrate stability on the waterwater pair potential model, we also investigate the dependence of the stability on the water-H2 potential. Previous modeling of hydrogen-containing clathrates has utilized the SPC/E potential, apparently successfully.9-11 The simulations performed in the present work indicate that at the temperatures used in the previous work the SPC/E potential yields an unstable clathrate structure, with slabs of the material decomposing quickly. The fact that the clathrates modeled by Alavi et al. were not observed to decompose can readily be explained by the lack of a melt-nucleating material surface in those simulations. Conclusions drawn from the previous modeling should be assessed carefully in light of the fact that the water-water potential used does not give a satisfactory description of the stability of the clathrates. In both slab and bulk simulations, two different levels of hydrogen mobility were observed, depending on whether the hydrogen molecules in question were contained in the small or large clathrate cages. Although no migration of the small-cage hydrogen was observed in multiple nanosecond simulations, the hydrogen in large cages was remarkably mobile. These two levels of hydrogen mobility were supported by the calculation of energy barriers along migration paths, with migration through hexagonal faces between large cages passing a barrier less than 30 kJ/mol and the migration through pentagonal faces required for changes to the small-cage occupation having calculated barriers over 200 kJ/mol. However, hydrogen must migrate to and from the small cages on longer timescales because neutron diffraction experiments of D2-containing clathrates have observed changes to the small-cage occupation as the temperature of a single sample was changed.3 An estimate of hydrogen tunneling through the migration barrier associated with hydrogen motion through the hexagonal clathrate cage face indicated an increase in the calculated migration rate due to tunneling of 25% at the highest temperature considered. A large degree of variability of the occupation of the large cages was observed, both spatially and temporally. In the slab calculations, this intercage hydrogen mobility lead to migration of the hydrogen out of the clathrate phase into the gas phase. The mobility of the hydrogen within the clathrate would presumably be greatly affected by the presence of promoter molecules, such as THF, occupying the large cages. Even in concentrations below that leading to full occupation of the large clathrate cages by promoters, the presence of promoter molecules would likely block the hydrogen migration paths. Though not directly investigated in this work, one may speculate about the role of reduced hydrogen mobility in the stabilizing effect of promoters like THF, not allowing hydrogen to migrate out of the clathrate structure by blocking the migration paths through the large cages. The hydrogen contained within the large cages would thus maintain “mechanical” resistance against clathrate cage collapse, in addition to that provided directly by the promoter molecules. In preliminary calculations, we have directly demonstrated the stabilizing effect of THF. In our simulations at 500 atm, the clathrate decomposition temperature was increased by an amount on the order of 20 K. This work lays some groundwork for investigations seeking insight into the role of promoter molecules in hydrogencontaining clathrates, aiming to identify promoter molecules with optimal properties. These investigations, potentially leading to materials optimized for hydrogen storage applications, are ongoing.
13052 J. Phys. Chem. C, Vol. 111, No. 35, 2007 Acknowledgment. Some of the calculations reported in this paper were performed by Sebastiaan Nooij and Johan Blanker. This work was supported by an STW/NWO grant and with computer time from NCF/NWO. References and Notes (1) Sloan, E. D. Clathrate Hydrates of Natural Gases, 2nd ed.; CRC Press: Boca Raton, FL, 1998. (2) Mao, W. L.; Mao, H. K.; Goncharov, A. F.; Struzhkin, V. V.; Guo, Q.; Hu, J.; Shu, J.; Hemley, R. J.; Somayazulu, M.; Zhao, Y. Science 2002, 297, 2247. (3) Lokshin, K. A.; Zhao, Y.; He, D.; Mao, W. L.; Mao, H. K.; Hemley, R. J.; Lobanov, M. V.; Greenblatt, M. Phys. ReV. Lett. 2004, 93, 125503. (4) Lokshin, K. A.; Zhao, Y. Appl. Phys. Lett. 2006, 88, 131909. (5) Florusse, L. J.; Peters, C. J.; Schoonman, J.; Hester, K. C.; Koh, C. A.; Dec, S. F.; Marsh, K. N. ; Sloan, E. D. Science 2004, 306, 469. (6) Lee, H.; Lee, J.; Kim, D. Y.; Park, J.; Seo, Y. T.; Zeng, H.; Moudrakovski, I. L.; Ratcliffe, C. I.; Ripmeester, J. A. Nature 2005, 434, 743. (7) Hester, K. C.; Strobel, T. A.; Sloan, E. D.; Koh, C. A.; Huq, A.; Schultz, A. J. J. Phys. Chem. B 2006, 110, 14024. (8) Strobel, T. A.; Taylor, C. J.; Hester, K. C.; Dec, S. F.; Koh, C. A.; Miller, K. T.; Sloan, E. D. J. Phys. Chem. B 2006, 110, 17121. (9) Alavi, S.; Ripmeester, J. A.; Klug, D. D. J. Chem. Phys. 2005, 123, 024507. (10) Alavi, S.; Ripmeester, J. A.; Klug, D. D. J. Chem. Phys. 2006, 124, 014704. (11) Alavi, S.; Ripmeester, J. A.; Klug, D. D. J. Chem. Phys. 2006, 124, 204707. (12) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (13) Ba´ez, L. A.; Clancy, P. J. Chem. Phys. 1994, 101, 9837. (14) Motakabbir, K. A.; Berkowitz, M. J. Phys. Chem. 1990, 94, 8359. (15) Chialvo, A. A.; Cummings, P. T. J. Phys. Chem. 1996, 100, 1309. (16) Smith, W.; Forester, T. R. J. Mol. Graphics 1996, 14, 136. (17) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (18) Bernal, J. D.; Fowler, R. H. J. Chem. Phys. 1933, 1, 515. (19) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (20) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (21) Abascal, J. L. F.; Sanz, E.; Ferna´ndez, R. G.; Vega, C. J. Chem. Phys. 2005, 122, 234511. (22) Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2005, 123, 234505. (23) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910. (24) Nada, H.; van der Eerden, J. P. J. M. J. Chem. Phys. 2003, 118, 7401.
Frankcombe and Kroes (25) Hixson, H. G.; Wojcik, M. J.; Devlin, M. S.; Devlin, J. P.; Buch, V. J. Chem. Phys. 1992, 97, 753. (26) Buck, U.; Hulsken, F.; Kohlhase, A.; Otten, D. J. Chem. Phys. 1983, 78, 4439. (27) Hodges, M. P.; Wheatley, R. J.; Schenter, G. K.; Harvey, A. H. J. Chem. Phys. 2004, 120, 710. (28) Faure, A.; Valiron, P.; Wernli, M.; Wiesenfeld, L.; Rist, C.; Noga, J.; Tennyson, J. J. Chem. Phys. 2005, 122, 221102. (29) Wind, P.; Røwggen, I. Chem. Phys. 1992, 167, 263. (30) Cui, T.; Cheng, E.; Alder, B. J.; Whaley, K. B. Phys. ReV. B 1997, 55, 12253. (31) Diep, P.; Johnson, J. K. J. Chem. Phys. 2000, 112, 4465. (32) Diep, P.; Johnson, J. K. J. Chem. Phys. 2000, 113, 3480. (33) Boothroyd, A. I.; Martin, P. G.; Keogh, W. J.; Peterson, M. J. J. Chem. Phys. 2002, 116, 666. (34) Carmichael, M.; Chenoweth, K.; Dykstra, C. E. J. Phys. Chem. A 2004, 108, 3143. (35) Silvera, I. F.; Goldman, V. V. J. Chem. Phys. 1978, 69, 4209. (36) Wang, Q.; Johnson, J. K.; Broughton, J. Q. J. Chem. Phys. 1997, 107, 5108. (37) Herna´ndez de la Pen˜a, L.; Gulam Razul, M. S.; Kusalik, P. G. J. Chem. Phys. 2005, 123, 144506. (38) Herna´ndez de la Pen˜a, L.; Kusalik, P. G. J. Chem. Phys. 2006, 125, 054512. (39) Fermann, J. T.; Auerbach, S. J. Chem. Phys. 2000, 112, 6787. (40) Peters, B.; Bell, A. T.; Chakraborty, A. J. Chem. Phys. 2004, 121, 4461. (41) Ladd, A. J. C.; Woodcock, L. V. Chem. Phys. Lett. 1977, 51, 155. (42) Bryk, T.; Haymet, A. D. J. Mol. Sim. 2004, 30, 131. (43) Ferna´ndez, R. G.; Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2006, 124, 144506. (44) Kroes, G. J. Surf. Sci. 1992, 275, 365. (45) Kroes, G. J. and Clary, D. C. J. Phys. Chem. 1992, 96, 7079. (46) Sabo, D.; Rempe, S. B.; Greathouse, J. A.; Martin, M. G. Mol. Sim. 2006, 32, 269. (47) Weissmann, M.; Ramı´rez, R.; Kiwi, M. Phys. ReV. B 1992, 46, 2577. (48) Sakai, H. Surf. Sci. 1996, 348, 387. (49) Mao, W. L.; Mao, H. K. Proc. Natl. Acad. Sci. 2004, 101, 708. (50) Vega, C.; Abascal, J. L. F.; Nezbeda, I. J. Chem. Phys. 2006, 125, 034503. (51) Vega, C.; Sanz, E.; Abascal, J. L. F. J. Chem. Phys. 2005, 122, 114507. (52) Kizilyalli, M.; Corish, J.; Metselaar, R. Pure Appl. Chem. 1999, 71, 1307. (53) Patchkovskii, S.; Tse, J. S. Proc. Nat. Acad. Sci. 2003, 100, 14645.
