Ab initio energetics and kinetics study of H_2 and CH_4 in the SI Clathrate Hydrate
Qi Li, Brian Kolb, Guillermo Roman-Perez, Jose M. Soler, Felix Yndurain, Lingzhu Kong, D. C. Langreth, T. Thonhauser
AAb initio energetics and kinetics study of H and CH in the SI Clathrate Hydrate Qi Li, Brian Kolb, Guillermo Rom´an-P´erez, Jos´e M. Soler, FelixYndurain, Lingzhu Kong, D. C. Langreth, and T. Thonhauser ∗ Department of Physics, Wake Forest University, Winston-Salem, NC 27109, USA Departamento de F´ısica de la Materia Condensada,Universidad Aut´onoma de Madrid, Cantoblanco, 28049 Madrid, Spain Department of Physics and Astronomy, Rutgers University, Piscataway, New Jersey 08854, USA (Dated: November 1, 2018)We present ab initio results at the density functional theory level for the energetics and kineticsof H and CH in the SI clathrate hydrate. Our results complement a recent article by some of theauthors [G. Rom´an-P´erez et al. , Phys. Rev. Lett. , 145901 (2010)] in that we show additionalresults of the energy landscape of H and CH in the various cages of the host material, as well asfurther results for energy barriers for all possible diffusion paths of H and CH through the waterframework. We also report structural data of the low-pressure phase SI and the higher-pressurephases SII and SH. PACS numbers: 66.30.je, 71.15.Mb, 84.60.Ve, 91.50.Hc
Clathrate hydrates are crystalline, ice-like structuresformed out of water molecules. The water frameworkcreates cavities in which gas molecules—typically O ,H , CO , CH , Ar, Kr, Xe—can be trapped, which sta-bilize the framework. The existence of clathrates wasfirst documented in 1810 by Sir Humphry Davy, andclathrates became the subject of intensive studies in the1930s, when oil companies became aware that clathratescan block pipelines. Nowadays, clathrate hydrates areof particular interest for two reasons: (i) they are formednaturally at the bottom of the ocean, where they areoften filled with CH . These deposits mean a tremen-dous stock pile of energy, while—at the same time—representing a possible global warming catastrophe if re-leased uncontrolled into the environment through melt-ing; (ii) clathrate hydrates can be used to store H inits cavities and can be a viable hydrogen-storage ma-terial (albeit with moderate hydrogen-storage density). For both cases, an understanding of the interaction be-tween the guest molecule and the host framework is cru-cial for their formation and melting processes, which arestill understood poorly. In this brief report, we presentresults that elucidate this crucial guest-molecule/host-framework interaction and complement a recent paper bysome of the authors. We show additional results of theenergy landscape of H and CH in the various cages ofthe host material, and we show further results for energybarriers for all possible diffusion paths of H and CH through the water framework. We also report structuraldata of the phases SI, SII, and SH.At low pressure, the methane filled clathrate forms thestructure SI, consisting of two types of cages. The smallercage is built of water molecules on the vertices of 12 pen-tagons with a diameter of approximately 7.86 ˚A, and werefer to this as 5 cage, or alternatively as D cage. Thelarger cage is built of 12 pentagons and two hexagonswith a diameter of approximately 8.62 ˚A, and we call it5 or T cage. The unitcell has cubic symmetry andconsists of two 5 and six 5 cages, with a total of 46 water molecules. At 250 MPa, the structure SI trans-forms into a new cubic phase SII, consisting of sixteen 5 and eight 5 cages, containing 136 water molecules inits unitcell. When the pressure is increased to 600 MPa,the structure undergoes another phase transition to thehexagonal phase SH. This phase has a smaller unit-cell of three 5 , two 4 , and one 5 cages, withonly 34 water molecules. Very nice graphical representa-tions of the different cages and structures can be foundin Refs. [2, 4, 6, and 8]. While other clathrate-hydratestructures exist, structure SI, SII, and SH are the mostcommon ones. Guest molecules such as H and CH in the cavities ofthe clathrate hydrates interact with the water frameworkthrough van der Waals forces. But even the water frame-work itself, i.e. the interaction of water molecules throughhydrogen bonds, has a van der Waals component. Tocapture these effects, we perform here density functionaltheory (DFT) calculations utilizing the truly non-localvdW-DF functional, which includes van der Waals inter-actions seamlessly into DFT.
We implemented vdW-DF using a very efficient FFT formulation into the lat-est release of PWscf , which is a part of the
Quantum-Espresso package. For our calculations we used ul-trasoft pseudopotentials with a kinetic energy cutoff forwave functions and charge densities of 35 Ry and 280 Ry,respectively. A self-consistency convergence criterion ofat least 1 × − Ry was used. All structures were fullyoptimized with respect to volume and atom positions,and the force convergence threshold was at least 10 − Ry/a.u. for SI and SH. We have also performed structuralcalculations on SII, but—due to the large unit cell with136 water molecules, i.e. 408 atoms—we used a slightlyless tight force convergence criterium of 5 × − Ry/a.u.For SI and SH we used a 2 × × while for SII we performed Γ-point calculationsonly.The empty cages are experimentally not stable, butthey have been shown to be a good starting point for a r X i v : . [ c ond - m a t . m t r l - s c i ] O c t TABLE I. Calculated and experimental lattice constants a and c for the SI, SII, and SH clathrate hydrates. In addition,calculated and experimental average nearest-neighbor and next-nearest-neighbor distances are given, as well as bond angles.Standard deviations are provided in square brackets. Experimental values for the lattice constants are taken from Ref. [8]for methane-filled cages. Experimental values for the averaged quantities are calculated from the structures given in thesupplemental materials of Ref. [8]. The experimental distances d nnO − H ( d nnnO − H ) seem to be underestimated (overestimated), mostlikely due to the difficulty of accurately determining H positions in X-ray experiments. For SI, neutron scattering experimentssuggest d nnO − H = 0 .
97 ˚A and d nnO − O = 2 .
755 ˚A. Also note that there is some variation in the experimental results for the latticeconstants in Refs. [1, 4, and 8]. a [˚A] c [˚A] d nnO − H [˚A] d nnnO − H [˚A] d nnO − O [˚A] (cid:54) H − O − H [ ◦ ] (cid:54) O − O − O [ ◦ ]SI calc. 11.97 — 0.994 [0.001] 1.790 [0.014] 2.781 [0.013] 107.1 ◦ [1.0] 108.6 ◦ [4.0]exp. 11.88 — 0.861 [0.031] 1.911 [0.022] 2.761 [0.017] 109.3 ◦ [3.0] 108.7 ◦ [3.7]SII calc. 17.35 — 0.994 [0.001] 1.792 [0.016] 2.784 [0.016] 107.1 ◦ [0.6] 109.2 ◦ [4.3]exp. 17.19 — 0.812 [0.016] 1.959 [0.025] 2.768 [0.013] 109.5 ◦ [2.0] 109.3 ◦ [4.0]SH calc. 12.32 10.01 0.994 [0.001] 1.793 [0.015] 2.782 [0.011] 107.2 ◦ [0.9] 108.4 ◦ [8.5]exp. 12.33 9.92 0.781 [0.040] 1.955 [0.022] 2.775 [0.005] 108.9 ◦ [5.1] 108.4 ◦ [8.3] calculations like ours. We have calculated the optimizedlattice parameters for the SI, SII, and SH structures andthe results are collected in Table I. We have also cal-culated the structures when filled with methane (onemethane molecule per cage) and filled with hydrogen(up to four H per cage), but the lattice parameters ex-pand less than 0.1% upon filling, such that we have usedthe parameters for the empty cages henceforth. Over-all, our optimized lattice constants agree well with pre-vious calculations and experiment. In Table I we fur-ther analyze the structure of the host materials by cal-culating the average nearest-neighbor and next-nearest-neighbor distances and important bond angles. In gen-eral, the calculated average distances and angles varyonly insignificantly amongst SI, SII, and SH, whereasthey show a slightly larger spread for some experimentalvalues. As a side-note, for a single water molecule we cal-culate d O − H = 0 .
973 ˚A and (cid:54) H − O − H = 104 . ◦ , in goodagreement with the experimental numbers of 0.958 ˚A and104.5 ◦ . Note that vdW-DF is known to give slightly toolarge binding distances.
Small deviations are visiblein the distances d nnO − H and d nnnO − H , which in sum mostlycancel to give very good agreement with the experimen-tal O–O distances. Reference [8] also gives the O–O dis-tances for all structures explicitly as between 2.725 ˚A and2.791 ˚A, in remarkable agreement with our calculations.Also, our calculated angles (cid:54) O − O − O agree very well withexperiment. However, the good agreement between oxy-gen distances and angles—which describe the structureas a whole—is closely related to the agreement for thelattice constants.We next focus on the binding energies of guestmolecules in the SI structure. In particular, we studythe binding energies of CH and H in the D and T cages as a function of the number of molecules; resultsare depicted in Fig. 1 for calculations where moleculesare added to only one cage in the unitcell, while all othercages are kept empty. Here, we define the binding energyas the energy difference between the “water-framework +guest-molecules” system minus the energy of the single Number of Molecules -0.4-0.200.2 B i nd i ng E ne r g y / M o l e c u l e [ e V ] CH in TCH in DH in TH in D FIG. 1. Binding energies per molecule for different numbers ofCH and H molecules in the D and T cages of SI. While the D and T cages can store only one CH molecule, the D and T cage can store up to two and four H molecules, respectively. constituents. In case of n -fold occupied cages, we sub-tract n times the energy of the single molecule. Methaneis a large molecule compared to the cage sizes and it canbe seen that in both D and T only one methane moleculecan be stored. Upon adding another methane molecule,the binding energy increases drastically. The situationis different for the much smaller H molecules. In thesmaller D cage we can store up to two H molecules, butincreasing the number to three or four results in a positivebinding energy, i.e. work is required to place more thantwo molecules into this cage. Note that the binding en-ergy that we find for double H occupancy is rather small,i.e. − D cages (see e.g. Ref. [20]),there are also reports that propose double occupancy orthat find inconclusive evidence. On the other hand,the larger T cage can store four H molecules. If a fifthmolecule is added, it escapes through one of the hexago-nal faces into the neighboring, empty T cage. Our calcu- FIG. 2. Energy landscape in meV for a rotating methanemolecule in a D cage. The x and y axes correspond to ro-tations about two mutually perpendicular axis. At the (0 , D cage. The differencebetween the minimum and maximum energy is 22 meV. lated H storage capacity of four molecules in the T cagesis in agreement with experiment. The binding energy forone H molecule compares well with quantum-chemistrycalculations on isolated cavities, which give –0.123 eV. Overall, our binding energies are slightly smaller than theones in Ref. [6]. Note that quantum motions have beenneglected in our approach, which may play an importantrole in the binding process and when determining thecage occupancy. A more precise treatment requires thecomputation of the corresponding thermodynamic parti-tion function, as for example shown in Ref. [26]. Nev-ertheless, we consider our calculations for the bindingenergy an important first step that already reveals im-portant information.It is also interesting to study where and how the H andCH molecules bind in the cages. If only one molecule ispresent in the cages, it binds in the center of the cage.Rotations and small displacements of H in that situa-tion are on an energy scale of approximately 1 meV andapproach the accuracy of our calculations. At room tem-perature, such perturbations are thus easily thermally ac-tivated. Since the methane molecule is larger, it cannotmove/rotate as easily. We have studied the rotation of asingle methane molecule centered in the D and T cages asa function of rotation about two mutually perpendicularaxes. The energy landscape for this rotation is depictedin Fig. 2 for the D cage. The D cage with its twelvepentagons and the methane molecule have a related sym-metry, which allows us to choose the (0 ,
0) point of theplot such that all methane hydrogens point exactly to anoxygen of the host lattice. At this point, hydrogen bondsare created and the total energy is the lowest. Upon ro-tation of the methane, the hydrogen bonds break andthe energy increases. The difference between the lowestand highest point of this energy landscape is 22 meV,suggesting thermal activation of rotations at room tem-perature, and quantifying an experimental assumption. Progress along Path [%] R e l a t i v e E ne r g y [ e V ] H T → T p H T → T h H T → D p CH T → T p CH T → T h CH T → D p FIG. 3. Barriers to diffusion of H and CH through the waterframework along different paths, as a function the relativeprogress along the path. The plots are labeled as A x → B ,where A is the type of the starting cage, B is the type of theending cage, and x refers to either pentagon (p) or hexagon(h), indicating the opening being used for traversing. Thesymbols are the results from 12-image NEB calculations andthe lines are fitted cubic splines, serving as a guide to the eye. We have also studied the rotation of a methane moleculein a T cage and the results are very similar to the resultspresented in Fig. 2, with the difference that the maximalenergy barrier is slightly smaller, i.e. 18 meV, not sur-prising as the T cage is slightly larger and the methanemolecule can rotate more easily. Calculations for bothrotation-energy landscapes have independently also beenperformed using Siesta and give essentially identicalresults.Finally, we present results for barriers to diffusionthrough the water framework in the SI structure. In thecase of H and hydrogen storage this is of much interest,as practical storage solutions require fast kinetics, i.e. lowbarriers. In the case of CH , the barriers can help us un-derstand the natural formation of the filled clathrates.We have calculated the barriers to diffusion with nudgedelastic band (NEB) calculations, using 12 images alongthe path from the center of one cage to the center of thenext cage; the results for all possible paths are plottedin Fig. 3. Note that the relaxation of the host latticeis crucial to obtain accurate barrier energies, and NEBcalculations allow for such relaxations perpendicular tothe path automatically. The plots are labeled as A x → B ,where A is the type of the starting cage, B is the typeof the ending cage, and x refers to either pentagon (p)or hexagon (h), indicating the opening being used fortraversing. Note that for the path T p → D there is onlyone choice of opening, i.e. a pentagon, and the path isnot symmetric as the distance from the center of T toits edge is longer than the corresponding distance in the D cage. Furthermore, this path’s end energy is differ-ent from its starting energy, since the guest moleculesare binding with different binding energies in the cages T and D , as already evident from Fig. 1. The lowestbarriers for H and CH diffusion agree well with previ-ous calculations. But, our H diffusion barrier is anoverestimation with respect to a recent NMR experiment,which gives 0.03 eV and warrants further investigation. The barriers are in general smaller for diffusion throughhexagons, simply because these openings are larger.For hydrogen-storage applications, the low barrier of ∼ T cages (going through a hexagon) isimportant. Through these T -cage channels, which threadthrough the material in all three dimensions, the hydro-gen can quickly be absorbed or released. However, toachieve the material’s full storage potential, some hy-drogen molecules will also have to get into the D cages,with a much higher barrier of ∼ ∼ and CH in hydrate clathrates. We havealso shown first results for the difficult-to-model, high-pressure phase SII with a large unit cell, finding goodagreement with experiment. While we have used vdW-DF for our study, it is conceivable that its successor,vdW-DF2, may further improve upon our results. Weencourage additional studies of the hydrate clathrates us-ing vdW-DF2, also including other types of cages, andmore detailed studies of the SII phase, which is one of themore promising phases amongst the hydrate clathratesfor hydrogen-storage applications.We would like to dedicate this report to the memory of Prof. David Langreth , who passed away just days beforeit was submitted—he is the “father” of vdW-DF and hisresearch inspired many. All calculations were performedon the WFU DEAC cluster. ∗ E-mail: [email protected] E. D. Sloan and C. A. Koh,
Clathrate Hydrates (CRC,Boca Ration, 2008). W. L. Mao, C. A. Koh, and E. D. Sloan, Physics Today , 42 (2007). B. Buffett and D. Archer, Earth Planet. Sci. Lett. ,185 (2004). V. V. Struzhkin, B. Militzer, W. L. Mao, H.-k. Mao, andR. J. Hemley, Chem. Rev. , 4133 (2007). S. Gao, W. House, and W. G. Chapman, J. Phys. Chem.B , 19090 (2005). G. Rom´an-P´erez, M. Moaied, J. M. Soler, and F. Yn-durain, Phys. Rev. Lett. , 145901 (2010). Here, we define the diameter as twice the average distancebetween oxygen atoms and the cage center. M. T. Kirchner, R. Boese, W. E. Billups, and L. R. Nor-man, J. Am. Chem. Soc. , 9407 (2004). B. Kolb and T. Thonhauser, Phys. Rev. B , 045116(2011). M. Dion, H. Rydberg, E. Schr¨oder, D. C. Langreth, andB. I. Lundqvist, Phys. Rev. Lett. , 246401 (2004); ,109902(E) (2005). T. Thonhauser, V. R. Cooper, S. Li, A. Puzder, P.Hyldgaard, and D. C. Langreth, Phys. Rev. B , 125112(2007). D. C. Langreth, B. I. Lundqvist, S. D. Chakarova-K¨ack,V. R. Cooper, M. Dion, P. Hyldgaard, A. Kelkkanen, J.Kleis, L. Kong, S. Li, P. G. Moses, E. Murray, A. Puzder,H. Rydberg, E. Schr¨oder, and T. Thonhauser, J. Phys.:Condens. Mattter , 084203 (2009). G. Rom´an-P´erez and J. M. Soler, Phys. Rev. Lett. ,096102 (2009). P. Giannozzi et al. , J. Phys.: Condens. Matter , 395502(2009). H. J. Monkhorst and J. D. Pack, Phys. Rev. B , 5188(1976); J. D. Pack and H. J. Monkhorst, Phys. Rev. B , 1748 (1977). C. Gutt, B. Asmussen, W. Press, M. R. Johnson, Y. P.Handa, and J. S. Tse, J. Chem. Phys. , 4713 (2000). G. Graner, E. Hirota, T. Iijima, K. Kuchitsu, D. A. Ram-say, J. Vogt, and N. Vogt,
2. Inorganic Molecules. Part4 , K. Kuchitsu (ed.), SpringerMaterials – The Landolt-B¨ornstein Database. DOI: . T. Thonhauser, A. Puzder, and D. C. Langreth, J. Chem.Phys. , 164106 (2006). K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist, and D. C.Langreth, Phys. Rev. B , 081101(R) (2010). K. A. Lokshin, Y. Zhao, D. He, W. L. Mao, H.-K. Mao, R.J. Hemley, M. V. Lobanov, and M. Greenblatt, Phys. Rev.Lett. , 125503 (2004). S. Patchkovskii and J. S. Tse, Proc. Nat. Acad. Sci. ,14645 (2003). T. M. Inerbaev, V. R. Belosludov, R. V. Belosludov, M.Sluiter, and Y. Kawazoe, Comp. Matt. Sci. , 229 (2006). F. Sebastianelli, M. Z. Xu, Y. S. Elmatad, J. W.Moskowitz, and Z. Bacic, J. Phys. Chem. C , 2497(2007). L. J. Florusse, C. J. Peters, J. Schoonman, K. C. Hester, C.A. Koh, S. F. Dec, K. N. Marsh, and E. D. Sloan, Science , 469 (2004). J. Wang, H. Lu, J. A. Ripmeester, and U. Becker, J. Phys.Chem. C , 21042 (2010). S. Patchkovskii and S. N. Yurchenko, Phys. Chem. Chem.Phys. , 4152 (2004). J. M. Soler, E. Artacho, J. D. Gale, A. Garcia, J. Jun-quera, P. Ordej´on, and D. S´anchez-Portal, J. Phys. Con-dens. Matter , 2745 (2002). P. Ordej´on, E. Artacho, and J. M. Soler, Phys. Rev. B ,R10441 (1996). S. Alavi and J. A. Ripmeester, Angew. Chem. , 6214(2007). T. Okuchi, I. L. Moudrakovski, and J. A. Ripmeester,Appl. Phys. Lett.91