THE JOURNAL OF CHEMICAL PHYSICS 131, 214506 共2009兲

Augmented stability of hydrogen clathrate hydrates by weakly polar
molecules
Takato Nakayama, Kenichiro Koga, and Hideki Tanakaa兲
Department of Chemistry, Faculty of Science, Okayama University, 3-1-1 Tsushima-naka,
Kitaku, Okayama 700-8530, Japan

共Received 7 October 2009; accepted 13 November 2009; published online 4 December 2009兲
Thermodynamic stability of hydrogen clathrate hydrates has been examined in a wide range ofpressure based solely on the intermolecular interactions involved. We show that the stability is
indeed augmented by a second guest species 共here acetone兲 called a promoter, a consequence of
which is notable reduction in the dissociation pressure of the hydrates encaging hydrogen alone.
This evaluation is made by extension of the van der Waals–Platteeuw theory combined withsemi-grand-canonical Monte Carlo 共GCMC兲 simulations where the number of hydrogen molecules
is allowed to vary while those of host water and promoter acetone molecules are fixed. The GCMC
simulations then provide various types of cage occupancies of hydrogen from single to quadruple,
from which the chemical potential of water in the clathrate hydrate is obtained as a function of thecage occupancy by acetone and the pressure. These occupancies are used to calculate the chemical
potential of water in the clathrate hydrate. The stability is estimated by comparison of the chemical
potential of water in the clathrate hydrate with that in hexagonal ice. We show the extent to which
the dissociation pressure is reduced with increasing the occupancy of the larger cages by acetone.
© 2009 American Institute of Physics. 关doi:10.1063/1.3271341兴
I. INTRODUCTION

Clathrate hydrates are nonstoichiometric compounds
made from water and guest molecules.1 Clathrate hydrates
are stable only when guest molecules occupy most of thecages. It has long been believed that each cage accommodates at most one guest molecule.1 The thermodynamic stability of the clathrate hydrates has been accounted for by van
der Waals and Platteeuw 共vdWP兲 under the assumption of
single occupancy.2 Dissociation pressures of many clathrate
hydrates have been successfully predicted by the vdWP
theory with accurate intermolecular interactions.1 An assumption in the original vdWP theory that host water lattice
is not affected by inclusion of guest molecules is reasonable
for small guest molecules but it is not so for guest molecules
large enough to alter the host vibrational motions
significantly.3–5 It was shown that a simple modification of
the vdWP theory enables to treat multiple occupancy in asingle cage.6
Since the discovery of hydrogen clathrate hydrate, it attracted considerable attention due to its potential use as hydrogen storage.7,8 A hydrogen clathrate hydrate is identified
with structure II, which is made by both 12- and 16-hedra.
Considering the size of a hydrogen molecule, multiple occupancies are expected at a high pressure where a hydrogenclathrate hydrate is formed.9 In fact, each larger cage accommodates maximally four hydrogen molecules while a smaller
one encages at most two molecules according to the originalreport. It is, however, pointed out later that the highest occupancy in the smaller cage seems to be limited to be one.
The cage occupancy of hydrogen clathrate hydrate hasa兲Electronic mail: htanakaa@cc.okayama-u.ac.jp.

0021-9606/2009/131共21兲/214506/10/$25.00

been examined by grand canonical/NPT Monte Carlo 共GC/
NPT MC兲 simulations for wide ranges of temperature andpressure.10 The simulations were carried out with a fixed
number of water molecules and a fixed chemical potential of
the guest species so that hydrogen molecules can be createdor removed in the clathrate. This hybrid type of the GCMC
simulation turned out to be very useful to calculate the cageoccupancy correctly, in which it is allowed to adjust itself
under a preset pressure so as to take account of compression
or expansion by a hydrostatic pressure. It was found that thesmaller cage in structure II is practically incapable of accommodating more than a single guest molecule even at pressures as high as 500 MPa, which agrees with the recent experimental investigations.8,11,12 The larger cage was found toencapsulate maximally four hydrogen molecules but its occupancy is dependent significantly on the chemical potentialof hydrogen, which is related to the pressure of the bulkhydrogen outside.
One of the key issues toward the practical use of a hydrogen clathrate is to explore possibilities for reduction of its
dissociation pressure. It is stable only at pressures above 200
MPa and its cage occupancy by hydrogen is almost saturated
around at 500 MPa. Such high pressures make it difficult tohandle this substance by simple containers. To overcome thisinconvenience, some molecules called promoters are added
to increase its stability. They are mostly large molecules accommodated in the larger cages so as to stabilize the clathrate while hydrogen molecules occupy both the smaller andthe remaining larger cages. For example, the dissociationpressure can be fairly reduced by accommodating tetrahydrofuran 共THF兲.13–17 This is advantageous to its use as hydrogenstorage because a moderate pressure such as 15 MPa suffices

131, 214506-1© 2009 American Institute of Physics

Downloaded 19 Jan 2010 to 150.46.100.222. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

214506-2

J. Chem. Phys. 131, 214506 共2009兲

Nakayama, Koga, and Tanaka

TABLE I. Geometries and parameters of intermolecular interactions for water, acetone, and hydrogen. Theinteractions are described by charge 共q兲, LJ size 共␴兲, and energy 共␧兲 parameters. The coordinates are those of
the sites with the origin of the center of gravity.
PotentialSite␴ / 0.1 nm

␧ / kJ mol−1

q/e

x / 0.1 nm

y / 0.1 nm

z / 0.1 nm

TIP5P

O
H
ChargeC
O
CH3H
Center

3.12

0.669

3.75
2.96
3.91

0.43960.87930.66993.038

0.285200.241⫺0.2410.300⫺0.4240.062
0.4932⫺0.986400⫾0.57200⫾1.2864000⫾0.757000.0893
1.3112⫺0.696⫾0.370700.06510.5859
0.404
0
0
0
0
0

AcetoneHydrogenfor its stabilization.13–17 However, it accompanies a serious
disadvantage. That is, occupation by a second guest necessarily reduces the occupancy of hydrogen. Thus, there must
be a kind of tradeoff between reduction in the dissociationpressure and its efficiency. There have been several attempts
to seek promoter species and to determine optimal occupancy by the promoter that achieves a low dissociation pressure and high occupancy by hydrogen.18–20
It is highly desirable to establish a relation between reduction of the dissociation pressure and the efficiency as hydrogen storage in a wide range of pressure with varying composition of a promoter species. To this end, we should find asimple way to evaluate the thermodynamic stability of clathrate hydrates and its cage occupancy from intermolecularinteractions currently available under a fixed occupancy of apromoter guest. Instead of THF, we choose acetone as a promoter. Acetone seems to be a little harder to be encaged instructure II clathrate hydrate.21 However, the resultant stability does not depend on the promoter species as far as it isencaged. We first extend the vdWP theory to a fixed numberof promoter molecules, without taking account of an equilibrium of a promoter with that outside explicitly. This condition is justified by the fact that acetone is much more preferentially encaged than hydrogen. It is also helpful toperform GCMC simulations at high pressures where multiple
occupancy is expected.
We propose a method to estimate the stability ofhydrogen-acetone binary clathrate hydrates by the extended
vdWP theory combined with GC/NPT MC simulations. This
seems to be the only way to implement numerical evaluationof the free energy for multiple cage occupancy since that for
triple or quadruple is actually impossible to calculate. Various cage occupancies are obtained from GC/NPT MC simulations in the presence of acetone. Once a certain amount ofpromoter species is encapsulated, its stability is dependent on
its occupancy but is irrelevant to the choice of it in neglecting the guest-guest interaction. Thus, the thermodynamic stability is evaluated with the chemical potential values of theempty hydrate and ice in equilibrium with the hydrate.
The present paper is organized as follows. Models andsimulation methods used in the present study are described inSec. II. The cage occupancies for hydrogen clathrate at various temperatures and pressures are shown and a relation be-

tween the vdWP theory and the GC/NPT hybrid ensemble
MC simulation is discussed in Sec. III. The conclusion is
given in Sec. IV.
II. SIMULATION METHOD

The choice of intermolecular interactions and the methods of calculating the chemical potential of hydrogen fluid,
generating the crystalline structure of hydrate II, and implementing GC/NPT MC simulations, are briefly outlined here;
the details were described elsewhere.10 The methods for calculating various free energies are also described.
A. Intermolecular interactions and structureof clathrate hydrate II

The thermodynamic stability and the cage occupancy ofclathrate hydrate are evaluated with intermolecular potentials
which are reliable but as simple as possible. In the present
study, all interactions are assumed to be pairwise additive.
The water-water intermolecular interaction is described by
the TIP5P potential.22 This potential has often been used in
the study of water. It is believed to be the most reliable
within the framework of pair potential in reproducing phasebehavior of pure water.23 The TIP5P model consists of fiveinteraction sites; a positive charge, q, on each proton, two
negative charge 共−q兲 sites extended from the oxygen toward
the remaining tetrahedral directions, also a Lennard-Jones
共LJ兲 potential between oxygen nuclei. Choice of the waterwater interaction, however, does not make a significant difference in cage occupancy under moderate compression,
which does not lead to collapse of the original clathrate hydrate structure or to significant hydrogen bond cleavage. An
acetone molecule is treated as a rigid rotor with four interaction sites, two methyl groups, a carbonyl carbon, and an
oxygen atom.24 The individual sites interact with each othervia LJ and Coulombic potentials. A hydrogen molecule ismodeled by a rigid rotor with two sites, 0.074 14 nm apart,
having a positive charge and one site, at the midpoint, having
an LJ potential center and a negative charge.25 For all theunlike site-site interactions, the Lorentz–Berthelot rule is applied. The parameters on the interaction potentials and themolecular geometries are tabulated in Table I. The interaction potentials for all pairs of molecules are truncatedsmoothly at 0.8655 nm.26 The correction term arising from

Downloaded 19 Jan 2010 to 150.46.100.222. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

214506-3

J. Chem. Phys. 131, 214506 共2009兲

Hydrogen clathrate hydrates

the truncation of all the LJ interactions is taken into accountin calculating the free energy of hydrogen fluid and the occupancy in clathrate hydrate.
A method to generate clathrate hydrate II structure wasdescribed previously.3 The unit cell of the hydrate II is cubicand experimental lattice parameters are used in the initialstage of GC/NPT MC simulation; a = b = c = 1.731 nm with␣ = ␤ = ␥ = ␲ / 2.1 We generate 32 different structures, each ofwhich is of proton-disordered form and has nearly zero net
dipole moment. The basic cell in the present MC simulations
is composed of 8 unit cells of clathrate hydrate II, in which
1088 water molecules are placed on the lattice sites.

B. GC/NPT Monte Carlo simulations

Accommodation of a guest molecule can be regarded as
an adsorption in each cavity. Thus, the number of guest molecules at a given pressure 共at a given chemical potential of
guest species in its fluid phase兲 can be evaluated in the same
fashion as usual adsorption process by GCMC simulation.
This simulation is carried out with the fixed parameters of
the temperature, the volume of the clathrate hydrate, V, andthe chemical potential of the guest species, ␮h 共or ␮a兲. The
chemical potential of the guest molecule is calculated fromthe equation of state as shown below. We assume that thehydrostatic pressure, P, on the hydrate is equal to that of thehydrogen fluid. The pressure covered in this study ranges
from 1 to 200 MPa, which may lead to alternation of the celldimension. Then, there is a more appropriate way to calculate the cage occupancy at a given temperature and pressure
as a function of the pressure of hydrogen in equilibrium withthe clathrate hydrate in question. That is, a standard MC
technique with a constant pressure is applied with the fixed
numbers of water and guest molecules. This meets the Gibbs
phase rule. A system is treated as an isothermal-isobaric
共NPT兲 ensemble with respect to water and guest at a given
instant and the volume change is attempted according to
NPT MC simulation scheme in addition to the grand canonical procedures. Thus, our simulation is a hybrid one, called
GC/NPT MC simulation. The volume is dependent on thehydrostatic pressure. A guest hydrogen molecule is inserted
or deleted according to the standard technique of GCMC
simulation at a given temperature T and pressure P. A trial
insertion is made at an arbitrary position with random distributions of polar and azimuthal angles. Such an insertion accompanied by the interaction energy w with all other molecules including Nh hydrogen molecules in a volume, V, is
accepted for a polar angle ␪ with a probability

再冎Vz sin ␪
exp共− ␤w兲 ,min 1,
共Nh + 1兲

共1兲

where z stands for the fugacity of guest molecule, and ␤ is
1 / kBT with the Boltzmann constant kB; z is calculated withthe thermal de Broglie wavelength ⌳ asid
兲 + ln ⌳3 ,ln z = ␤共␮h − ␮rot

共2兲id
are the chemical potential of hydrogen andwhere ␮h and ␮rot
its ideal rotor part, respectively. A deletion of an arbitrarily

chosen molecule out of Nh guests in the system is attempted.
The trial is accepted with a probability

再min 1,

冎Nhexp共␤w⬘兲 ,Vz sin ␪

共3兲

where w⬘ is the interaction energy of the chosen guest withall other molecules. The intermolecular interaction is truncated smoothly as in GCMC simulations, but its correction to
the pressure is taken into consideration with an assumptionthat molecular distribution is uniform beyond the truncation.
In the present GC/NPT MC simulations, a single MC
step consists of either a trial creation or annihilation 共withthe equal probability兲 followed by five trial moves of molecules, and a subsequent volume change. Each MC simulation at given temperature and pressure includes 4 ⫻ 107 steps.

C. Chemical potential for hydrogen fluid

The equation of state for the model hydrogen fluid
whose constituent molecules interact with the prescribed potential is unknown. We thus calculate the pressure-density
isotherms from NVT MC simulations and then obtain thechemical potential as a function of pressure for each temperature. The number of molecules used in the simulation istaken to be in a range between 128 and 512, depending on
the density, so that the volume falls into a certain range. The
density ranges from 0.2 to 64 mol L−1, which is enough forthe present purpose. The chemical potential of the imperfectgas 共fluid兲 for a linear rigid rotor is calculated as follows. At
a given temperature, we perform NVT MC simulations ofpure hydrogen with various volumes 共number densities, ␳兲.
The long-range interaction is again taken into account in thecalculation of the pressure. Thus, we obtain an isotherm of
the pressure.
The free energy per molecule, ah, at T and ␳ for hydrogen with the mass, m, and the moment of inertia, I, is given
by sum of the ideal and the nonideal parts as冋 冉 冊

ah = kBT ln ␳
+

冕

冉 冊册

␤h2 3/2
␤h2− 1 + ln2␲m
4 ␲ 2I␳p − ␳ k BT
d␳ ,
␳2
0

共4兲

where h is the Planck constant. Since we deal with a onecomponent bulk fluid system, the chemical potential, ␮h, issimply given by␮h = ah + P␳−1 ,共5兲

where ␳ at a given pressure P is calculated from the interpolation of the isotherm.

D. Calculation of free energy of solid and liquidstates and free energy of cage occupation

The free energy of ice is calculated by its individualcomponents, the interaction energy at 0 K, Uq共V兲, the harmonic free energy, Fh共T , V兲, the anharmonic free energy,
Fa共T , V兲, and the residual entropy term, Sc, as

Downloaded 19 Jan 2010 to 150.46.100.222. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

214506-4

J. Chem. Phys. 131, 214506 共2009兲

Nakayama, Koga, and Tanaka

A共T,V兲 = Uq共V兲 + Fh共T,V兲 + Fa共T,V兲 − TSc .

共6兲

It should be noted that the interaction energy at 0 K is afunction of volume alone, which is obtained from conventional NPT MC simulations. The harmonic vibrational energy is calculated from a set of the normal frequencies, ␯ as
Fh共T,V兲 = kBT 兺 ln共␤h␯i兲.

共7兲i

3 ⫻ 106 MC steps 共1 MC step= N trial move兲 for equilibration and 3 ⫻ 106 MC steps for averaging. This method isapplied to liquid acetone.
The free energy of cage occupation is calculated by asimple integration on coordinates and orientations of an acetone inside a larger cage as
f a = − kBT ln ␴

The anharmonic free energy is calculated as
Fa共T,V兲
=−
T

冕

T

Ua
2 dT,
0 T

共8兲

where Ua is the anharmonic energy given by
共9兲

Ua = U − Uq − Uh ,with Uh = 3共Nw − 1 / 2兲kBT for Nw water molecules. The residual entropy term is given by NwkB ln共3 / 2兲 for both ice andempty clathrate and does not actually contribute to thechemical potential difference. The Gibbs free energy affectedby pressure is evaluated by
共10兲

G = A + PV.

All those properties depend on the volume, and thus the pressure. However, their pressure dependences are, here, expected to be small in the pressure range examined. It is assumed that the compressibilities for ice and clathrate aresimilar to each other and the pressure dependence of thechemical potential difference between ice and clathrate arising from the compressibility difference can be neglected.
The free energy difference of a fluid with a potential
energy function ⌽N from the corresponding ideal gas state
Aid at the same density is calculated by thermodynamic integration,冕冕冓 ⬘ 冔
1⌬A = A − Aid =

dA共␭兲
d␭
d␭
0
1

=

0

⬇兺
k

⳵ ⌽N共␭兲⳵␭冓⳵ ⌽N⬘ 共␭k兲
⳵␭

冔

⌬␭,

共11兲

N,T,V,␭k␤h2 3/2−12␲m

冋 兿 冉 冊 册冊
␲
␴

− ln

8␲ I j
␤h2
2

j

v

⍀

exp关− ␤␾共r,⍀兲兴drd⍀ ,共13兲

where ␾共r , ⍀兲 is the interaction potential of an acetone withthe surrounding molecules. This method is applicable only toa single occupancy of a nonlinear rigid molecule at thepresent stage. In calculating the above quantity, it is foundthat the interaction between water and acetone in the present
model is not so strong as to form a hydrogen bond breaking
a hydrogen bond between water molecules; the lowest valueof ␾共r , ⍀兲 is found to be around −47 kJ mol−1, the magnitude of which is less than the half of the magnitude of thecorresponding energy for water.
E. Grand canonical MC simulation for acetone
clathrate hydrate

In order to perform GCMC simulation for a large molecule such as acetone, a different kind of GCMC simulation
method has been proposed and applied to propane.5 Thismethod is briefly outlined here. The GCMC simulation iscarried out by choosing one of the cages that can accommodate a guest molecule. If it is empty, a creation trial is attempted by designating position and orientation where aweighting function, ␹, is multiplied to facilitate convergenceand then accept-reject scheme is applied; a creation is accepted with a probability,
共14兲

where vc is the volume of each cage but does not appear in apractical GCMC simulation. The weighting function is chosen to bed␭
N,T,V,␭

冉 冉 冊
1/2

⫻

2␲m 3/2
2␲I j 1/2兿2␤h2
j=1 ␤hmin兵1,zvc sin ␪ exp共− ␤w兲␹−1其,where the perturbation potential is given by ⌽N⬘ 共␭兲 = ␭k⌽N
with k set to 8. In numerical calculation, the increment in ␭ ischosen to be 0.02. The ideal free energy is given by
Aid = NkBT ln ␳

再冉 冊 冉 冊
冕冕
冎
l

−1

冉

␹共r,⍀兲 = vc共␤␬/2␲兲3/2 exp − ␤␬

冊

r2
.
2

共15兲

The parameter ␬ is introduced to sample the central region ofa cage more preferentially and is set to 2␲kBT, where r is in
0.1 nm. If a chosen cage is occupied, a deletion trial is accepted with a probability,
min兵1,exp共␤w兲␹/zvc sin ␪其.

共16兲

1/2

,

共12兲

where ␴ is the symmetry number and I j stands for jth diagonalized component of the inertia tensor. In the actual calculation, NPT MC simulations for N = 300 are performed withthe steps of 3 ⫻ 106 共1 MC step= N trial move +1 volume
change兲 to determine the equilibrium volume, which arefollowed by NVT MC simulations at each ␭ value with

III. THEORY, SIMULATION RESULTS,
AND DISCUSSION

We extend the vdWP theory to binary systems where thenumber of promoter molecules is either allowed to alter or
fixed to a constant. Combined with GC/NPT MC simulations, the phase equilibrium for the binary clathrate hydrateis examined on the basis of the theory developed.

Downloaded 19 Jan 2010 to 150.46.100.222. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

214506-5

J. Chem. Phys. 131, 214506 共2009兲

Hydrogen clathrate hydrates

TABLE II. Free energy components of hexagonal ice and empty clathrate II
and chemical potential in kJ mol−1.
Uq

Fa

Fh

A

PV

␮

Ice
⫺57.4958 ⫺1.7918 8.9552 ⫺51.2609 0.1843 ⫺51.0766
Clathrate ⫺56.5531 ⫺1.7302 8.7946 ⫺50.4173 0.2158 ⫺50.2014
Difference⫺0.8752

FIG. 1. Temperature dependences of the chemical potential of acetone 共open
circle兲 and the free energy of cage occupation 共filled circle兲 at pressure 0.1
MPa.

stable than ice Ih. The obtained difference, ⌬␮ = ␮i − ␮0c= −0.88 kJ mol−1, is similar to that obtained previously
for a different water model.5 Its temperature dependence isroughly estimated from the entropy difference ⌬s = 共⌬h− ⌬␮兲 / T, where ⌬h is the enthalpy difference obtained by
the conventional NPT MC simulation. A similar estimation ismade for pressure dependence using the volume difference⌬v.

A. Free energy of cage occupation and chemical
potential of ice and empty hydrate

The free energy of acetone in its pure liquid is calculated
at several temperatures and pressures, which is converted tothe chemical potential. The free energy of the larger-cage
occupation by acetone is calculated by the numerical integration at several temperatures. Its temperature dependence isplotted in Fig. 1 along with the chemical potential at 0.1
MPa. It indicates that the chemical potential is considerably
higher than that of the free energy. This indicates that a great
number of the larger cages is occupied by acetone molecules,
which leads to stability of the acetone clathrate hydrate. 共Its
spontaneous formation is not confirmed from only this estimation but it must be done by comparison of the chemical
potentials of water in two phases.兲 Pressure dependence ofthe chemical potential of acetone at 273 K is plotted inFig. 2. It is evident from the slope of the curve that thevolume of pure acetone is almost constant over the pressure
range from 0 to 100 MPa.
The free energies of ice Ih and empty clathrate hydrate
structure II are calculated at 273 K and 10 MPa. The individual components are listed in Table II. The chemical potential of ice Ih is lower than that of the empty clathrate,
which ensures that structure II empty clathrate hydrate is less

B. Phase equilibrium and the van der Waals
and Platteeuw theory for binary clathrates

In the present study, occupancies of hydrogen and thepromoter acetone are examined. That is, we consider theequilibria among three phases 共p = 3兲 in the binary clathrate
hydrate with the components of water, hydrogen, and acetone 共the number of components, c = 3兲. The clathrate hydrate is in equilibrium with either aqueous solution or ice
coexisting with hydrogen fluid. It is known that the larger
cage accommodates maximally four hydrogen molecules
共m = 4兲 or at most one acetone molecule and the smaller one
contains zero or one hydrogen molecule. If both guest species are variable being subject to constant chemical potential
values according to the original vdWP theory, the grand partition function with the free energy of the empty clathrate Aw0
is given byNl⌶ = exp共− ␤Aw0 兲 兺

na=0

冦 冢
m

Nl

⫻ 兿兺
j=1 n j冉冊Nlexp关␤na共␮a − f a兲兴naj−1

Nl − na − 兺 nk
k=1

nj

冧 冉冊
Ns

⫻exp关␤n j共j␮h − f j兲兴冣Nsexp关␤ns共␮h − f s兲兴, 共17兲
兺
n
s
n =0
s

FIG. 2. Pressure dependence of the chemical potential of acetone at temperature 273 K.

where Nl and Ns stand for the numbers of the larger andsmaller cages in the clathrate hydrates and the free energies
of cage occupancy by acetone and by j hydrogen molecules
in a larger cage, and that by a hydrogen molecule in asmaller cage are denoted by f a, f j, and f s, respectively. Finally a simple form is obtained as

Downloaded 19 Jan 2010 to 150.46.100.222. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

214506-6

J. Chem. Phys. 131, 214506 共2009兲

Nakayama, Koga, and Tanaka

再molecule, and the smaller-cage occupancy xs by a singlehydrogen molecule are given, respectively, as⌶ = exp共− ␤Aw0 兲 1 + exp关␤共␮a − f a兲兴m+ 兺 exp关␤共j␮h − f j兲兴j=1

冎

x j = exp关␤共j␮h − f j兲兴

Nl

兵1 + exp关␤共␮h − f s兲兴其 .
Ns

共18兲

One of the independent variables, V, is changed to the pressure P applied on the clathrate hydrate as
⌼=

1
v0

冕

⬁

⌶ exp共− ␤ PV兲dV,

共19兲⳵ ln ⌼
⳵ 共␤␮h兲

再
再

= Nl

m

j exp关␤共j␮h − f j兲兴
兺
j=1

冎
m

⫻ 1 + exp关␤共␮a − f a兲兴 + 兺 exp关␤共j␮h − f j兲兴
j=1+ Ns exp关␤共␮h − f s兲兴兵1 + exp关␤共␮h − f s兲兴其−1

再m⫻ 1 + exp关␤共␮a − f a兲兴 + 兺 exp关␤共j␮h − f j兲兴⳵ ln ⌼
⳵ Nw

冎

−1

共20兲

j=1

冎

−1

共21兲

再

+ 兺 exp关␤共j␮h − f j兲兴
j=1

j=1

冎

冎

−1

,
共23兲

xs = exp关␤共␮h − f s兲兴兵1 + exp关␤共␮h − f s兲兴其−1 .
The occupancies, x, are inserted in the above equation to findthe equilibrium chemical potential of water given by
m

j=1

冊

We first examine an equilibrium among the clathrate, theaqueous solution, and the hydrogen fluid, although the equilibrium will not be treated explicitly in the present paper.
Since there are three phases 共p = 3兲 and three components
共c = 3; water, acetone, and hydrogen兲, the number of independent thermodynamic variables is given as f = c + 2 − p = 2 according to the phase rule. Here they are taken to be T and P.
The aqueous solution is assumed, hereafter, to be ideal forthe sake of simplicity 共its extension to nonideal solution isstraightforward introducing the activity兲. The chemical potentials of the aqueous solution and of the clathrate hydrate
are then described, respectively, by

␮w共T, P兲 = ␮w0 共T, P兲 + kBT ln共1 − y a兲,共25兲

␮a共T, P兲 = ␮0a共T, P兲 + kBT ln y a ,

共26兲

.

= ␤␮0c − ␣l ln 1 + exp关␤共␮a − f a兲兴
m

再m⫻ 1 + exp关␤共␮a − f a兲兴 + 兺 exp关␤共j␮h − f j兲兴

,

共24兲

Since the clathrate hydrate is treated as an NPT ensemble
with respect to water, the chemical potential of water is calculated as␤␮c = −

xa = exp关␤共␮a − f a兲兴冉⳵ ln ⌼
⳵ 共␤␮a兲= Nl exp关␤共␮a − f a兲兴冎−1␮c = ␮0c + ␣lkBT ln 1 − xa − 兺 x j + ␣skBT ln共1 − xs兲.

and具na典 =

j=1

0

where v0 is the unit volume.
The mean numbers of hydrogen and acetone molecules
are expressed as
具nh典 =

再m⫻ 1 + exp关␤共␮a − f a兲兴 + 兺 exp关␤共j␮h − f j兲兴− ␣s ln兵1 + exp关␤共␮h − f s兲兴其.
共22兲

The larger-cage occupancy 共fraction兲 x j by j hydrogenmolecules, the corresponding occupancy xa by an acetone

under the condition ␮w ⬍ ␮i, where y a is the mole fraction ofacetone in aqueous solution. The occupancies at equilibrium
are obtained by solving the equation ␮c共T , P兲 = ␮w共T , P兲
combined with Eq. 共23兲 assuming the chemical potential ofhydrogen is equal to the pure fluid as ␮h共T , P兲 ⬇ ␮0h共T , P兲.
Next, equilibrium among ice, clathrate, and hydrogen
fluid phases is considered. The number of components andthe number of phases are the same as in the previous equilibrium and so is the number of independent variables f = 2.
Those are again taken to be T and P. The chemical potentialof water in aqueous solution is substituted for that of ice. The
occupancies of hydrogen and acetone 共xh and xa兲 are obtained from

␮i共T, P兲 = ␮c共T, P兲,

共27兲

with ␮h共T , P兲 ⬇ ␮0h共T , P兲 under ␮i ⬍ ␮w. The relation of thechemical potential between ice and water in aqueous solution
may impose a restriction on the range of xa.

Downloaded 19 Jan 2010 to 150.46.100.222. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

214506-7

J. Chem. Phys. 131, 214506 共2009兲

Hydrogen clathrate hydrates

具nh典 =

⳵ ln ⌼
⳵ 共␤␮h兲= 共Nl − na兲

再

再

m

j exp关␤共j␮h − f j兲兴
兺
j=1

m

⫻ 1 + 兺 exp关␤共j␮h − f j兲兴
j=1

冎

冎

−1

+ Ns exp关␤共␮h − f s兲兴兵1 + exp关␤共␮h − f s兲兴其−1 , 共30兲␤␮c = −

= ␤␮0c + ␣l ln共1 − xa兲

再FIG. 3. Chemical potential difference between ice and clathrate at 273 K
and various pressures. Filled circle: calculated from the free energy
关left hand side of Eq. 共28兲兴. Open circle: calculated from the occupancy
共right hand side兲.

␮i共T, P兲 − ␮0c 共T, P兲 = ␣lkBT ln共1 − xa兲.

共28兲Both sides of the above equation are plotted in Fig. 3 at
273.15 K and several pressures. It is shown that the righthand side is always lower. Since the concentration of acetone
has an upper limit, this result simply indicates an acetone
clathrate is more stable than ice. However, the most stablephase can be an aqueous solution where the chemical potential of water is lower than that in clathrate. Therefore, thisleads to a conclusion that the combination of the intermolecular interactions suffice to reproduce acetone clathrate formation.

D. Phase equilibrium and the van der Waals
and Platteeuw theory for restricted binary clathrates

If the number of acetone molecules in the larger cages isfixed to na, the relevant grand partition function, the occupancy of hydrogen, and the chemical potential of water in theclathrate hydrate are given by⌶ = exp共− ␤Aw0 兲

再冉冊

⫻ 1 + 兺 exp关␤共j␮h − f j兲兴
j=1⫻兵1 + exp关␤共␮h − f s兲兴其Ns ,

冎

j=1

冎

− ␣s ln兵1 + exp关␤共␮h − f s兲兴其冉= ␤␮0c + ␣l ln

冊m

1 − xa − 兺 x j + ␣s ln共1 − xs兲,
j=1

共31兲

with
xj =

共1 − xa兲exp关␤共j␮h − f j兲兴.
兺k=0exp关␤共k␮h − f k兲兴

共32兲

Even when the number of acetone molecules is fixed to a
constant, the chemical potential of water can be expressed
apparently by the same equation. In case of an equilibrium
between hydrate and ice, there are three phases 共p = 3; clathrate, ice, and hydrogen fluid兲 with the components of two
species 共c = 2; water and hydrogen兲. From the phase rule, f= c + 2 − p = 1, and T can be an independent thermodynamic
variable while the occupancy xa 共=na / Nl兲 is kept fixed and isregarded as a parameter. The equilibrium between ice andclathrate for a given xa is obtained by solving冉␮i共T, P兲 = ␮0c 共T, P兲 + ␣lkBT lnm

1 − xa − 兺 x j
j=1

冊+ ␣skBT ln共1 − xs兲.

共33兲This equation is used for evaluation of stability of clathrate
by GC/NPT MC simulations with a fixed number of acetone
molecules. It is evident from Eq. 共32兲 that the smaller xa
gives rise to partial decomposition of the clathrate hydrate
into ice.
In passing, it can be shown that a set of xa and y a should
satisfy in order for a clathrate to be in equilibrium with aqueous solution,冉Nlexp共− ␤na f a兲
na

m

m

− ␣l ln 1 + 兺 exp关␤共j␮h − f j兲兴

C. GCMC simulations of acetone clathrate hydrates

GCMC simulations of acetone clathrate hydrates serve
to examine whether these potentials of water and acetone areappropriate to model the clathrate hydrates. This is actually
made by calculating the mean occupancy at several pressurestentatively using the chemical potential of pure liquid acetone at a given temperature and pressure. 共This condition isan artificial one since the high concentration of acetone inhibits formation of the clathrate.兲 The occupancy, xa, is related to the chemical potential difference between ice and theempty clathrate hydrate as⳵ ln ⌼
⳵ Nwm冊

␮0c 共T, P兲 + ␣lkBT ln 1 − xa − 兺 x j + ␣skBT ln共1 − xs兲

Nl−na

j=1

= ␮w0 共T, P兲 + kBT ln共1 − y a兲,
共29兲

共34兲

subject to the condition

Downloaded 19 Jan 2010 to 150.46.100.222. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

214506-8

J. Chem. Phys. 131, 214506 共2009兲

Nakayama, Koga, and Tanaka

FIG. 4. Pressure dependence of various types of hydrogen occupancies at
250 K in the larger cage without acetone. Dotted: single. Dashed: double.
Solid: triple. Dot-dash: quadruple.

冉

m

冊

FIG. 6. Pressure dependence of occupancies of hydrogen at 250 K in thesmaller cages with acetone occupancy of 0 共dotted兲, 0.5 共dashed兲, and 1
共solid兲.

Occupation by acetone prevents hydrogen moleculesfrom entering the larger cages, which thus undermines its
efficiency. Utilization of a binary clathrate hydrate as a hydrogen storage is solely dependent on the balance between

the extent of depression of the dissociation pressure and thechange in the occupancy of hydrogen. One of the goals in thepresent study is to evaluate the stability of hydrogen clathrate
hydrate as a function of the fraction of acetone in the larger
cages at various pressures, thereby making it possible to predict, via intermolecular interactions alone, whether it isformed at given conditions. To this end, calculation of thefree energy for multiple cage occupancy is required. It is,however, impossible or impractical to calculate f 3 and f 4, forthis involves a considerable numerical integration task. To
circumvent this difficulty, we take advantage of the fact that
evaluation of the chemical potential of water in the clathrate
hydrate requires only the various types of cage occupancies.
The latter quantities can be obtained from GC/NPT MC
simulations with the fixed number of acetone molecules. The
larger-cage occupancies of various types for xa = 0 and 0.5
are given in Figs. 4 and 5, and the smaller-cage occupancies
for xa = 0, 0.5, and 1 are shown in Fig. 6. The occupancy inthe larger cage decreases monotonically with the occupancyof acetone. A detailed comparison shows that the relative
occupancies of hydrogen to the available larger cages otherthan those encapsulating acetone remain almost constant irrelevant to the occupancy of acetone. However, a small en-

FIG. 5. Pressure dependence of various types of hydrogen occupancies at
250 K in the larger cage with half the larger cages occupied by acetone.
Dotted: single. Dashed: double. Solid: triple. Dot-dash: quadruple.

FIG. 7. Volume of the hydrogen clathrate at 1 MPa 共circle兲, 100 MPa
共square兲, 150 MPa 共diamond兲, and 200 MPa 共triangle兲 and at 250 K. Filled
symbols indicate clathrate is more stable than ice.

f a共T, P兲 + kBT ln xa − kBT ln 1 − xa − 兺 x j
j=1= ␮0a共T, P兲 + kBT ln y a ,共35兲

the left hand side of which is derived from ␮a= −kBT共⳵ ln ⌼ / ⳵na兲. There are possibly two intersections in
the chemical potentials, which limit the stable concentrationy a since the chemical potential of water decreases with increasing acetone concentration. This also poses the upperand lower limits on xa; the stable phase in the region between
the two limits is the clathrate hydrate alone. Thus, an acetone
concentration in an initial aqueous solution to prepare aclathrate hydrate cannot be too high for a clathrate hydrate tobe the only stable phase.
E. Occupancy of hydrogen with fixed numberof acetone molecules

Downloaded 19 Jan 2010 to 150.46.100.222. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

214506-9

J. Chem. Phys. 131, 214506 共2009兲

Hydrogen clathrate hydrates

TABLE III. Hydrogen occupation ratio to water 共H2 / H2O兲 in hydrogen clathrate hydrate with acetone at 250 K
and at various pressures and acetone occupancies. No occupation ratio is given to an unstable state.
P / MPa
xa

1

10

50

100

150

200

0
1/2
3/4
7/8
15/16
31/32
63/64
1

¯
¯
¯
¯
¯
¯
¯
0.008

¯
¯
¯
¯
¯
¯
¯
0.025

¯
¯
¯
¯
¯
¯
¯
0.075

¯
¯
¯
¯
¯
¯
¯
0.100

¯
¯
¯
¯
0.115
0.113
0.111
0.109

0.231
0.174
0.144
0.129
0.121
0.117
0.115
0.113

hancement in occupancy of hydrogen by promoter acetone isobserved in the smaller cage. This seems to stem from thefact that incorporation of acetone expands the volume 共Fig.
7兲 and also makes the amplitudes of the thermal vibrationalmotions in the host lattice smaller, which prevents the cages
from being deformed. This leads to easier incorporation ofhydrogen molecules although the increase in the occupancy
is small.
It is possible to evaluate the stability of the clathrate
hydrate at a given pressure with a fixed occupancy of acetone
from the chemical potential difference between ice andempty clathrate hydrate, Eq. 共33兲. This evaluation is madewith the various types of hydrogen occupancies. Its stability
relative to ice at 250 K is examined. Table III tabulates themean hydrogen ratio to water for the stable clathrate at agiven pressure. Unless there is a special correlation betweenhydrogen and acetone occupations, it is obvious from Eq.
共32兲 that the larger xa decreases x j at an equilibrium. A stablebinary clathrate hydrate must have a larger xa under a low
pressure of hydrogen. Thus, a clathrate containing a small
amount of acetone under high pressure of hydrogen becomes
unstable when it is decompressed. It is decomposed into ice
and clathrate with a larger xa.
IV. CONCLUDING REMARKS

Augmented stability of hydrogen clathrate hydrates by apromoter species has been investigated by an extension of
the vdWP theory. Phase equilibria of hydrogen clathrates
with either ice Ih or an aqueous solution are discussed. Acetone is chosen in the present study as a molecule to stabilizehydrogen clathrates but our theoretical treatment can be applied to any kind of binary clathrate hydrates with promoter
molecules including THF. First, it is shown that acetone isindeed encaged to form clathrate hydrate by GCMC simulation and the free energy of cage occupation. New GC/NPTMC simulations are carried out with a fixed number of acetone molecules and a fixed chemical potential of guest hydrogen molecules. Various types of the cage occupancies inhydrogen clathrate hydrate have been obtained by these
GCMC simulations in a wide range of hydrogen pressure
since the triple and quadruple occupancies are otherwise impossible to evaluate. These occupancies enable to calculate
the chemical potential of water in clathrate hydrate. Thus, we
can discuss augmented thermodynamic stability of hydrogenclathrates with the existence of a promoter guest species. It isrigorously shown that the augmented stability is due to theoccupation of cages by the promoter molecules and that apparent weak and indirect attraction between acetone and hydrogen in the smaller cage is observed, which seems to becaused by stabilization of the cage structure. The present
method is advantageous in establishing the efficiencystability relation with the simple intermolecular interactions.
We show a condition that a clathrate hydrate with a fixed
number of promoter molecules becomes stable. The partial
occupancy of the larger cage is indeed effective to reduce thedissociation pressure. Since the occupancy depends solely on
the affinity of promoter molecules to a clathrate, it is worthwhile discussing the stability of such a clathrate made ofwater and the promoter species in terms of the property in its
aqueous solution.

ACKNOWLEDGMENTS

This work was supported by grant-in-aid and the Next
Generation Super Computing Project, Nanoscience Program,
MEXT, Japan and also by Okayama Foundation for Scienceand Technology.
1

E. D. Sloan and C. A. Koh, Clathrate Hydrates of Natural Gases, 3rd ed.
共CRC, Boca Raton, 2008兲.
2
J. H. van der Waals and J. C. Platteew, Adv. Chem. Phys. 2, 1 共1959兲.
3
H. Tanaka and K. Kiyohara, J. Chem. Phys. 98, 4098 共1993兲.
4
H. Tanaka and K. Kiyohara, J. Chem. Phys. 98, 8110 共1993兲.
5
H. Tanaka, J. Chem. Phys. 101, 10833 共1994兲.
6
H. Tanaka, T. Nakatsuka, and K. Koga, J. Chem. Phys. 121, 5488 共2004兲.
7
W. L. Mao, H.-K. Mao, A. F. Goncharov, V. V. Struzhkin, Q. Guo, J. Shu,
R. J. Hemley, M. Somayazulu, and Y. Zhao, Science 297, 2247 共2002兲.
8
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. 93, 125503 共2004兲.
9V. V. Struzhkin, B. Militzer, W. L. Mao, H. Mao, and R. J. Hemley,
Chem. Rev. 共Washington, D.C.兲 107, 4133 共2007兲.
10
K. Katsumasa, K. Koga, and H. Tanaka, J. Chem. Phys. 127, 044509
共2007兲.
11
K. C. Hester, T. A. Strobel, E. D. Sloan, C. A. Koh, A. H. Hug, and A. J.
Schultz, J. Phys. Chem. B 110, 14024 共2006兲.
12
T. A. Strobel, C. J. Taylor, K. C. Hester, S. F. Dec, C. A. Koh, K. T.
Miller, and E. D. Sloan, J. Phys. Chem. B 110, 17121 共2006兲.
13
L. J. Florussse, C. J. Peters, J. Schoonman, K. C. Hester, C. A. Koh, S. F.
Dec, K. N. Marsh, and E. D. Sloan, Science 206, 2247 共2002兲.
14
T. A. Strobel, C. A. Koh, and E. D. Sloan, J. Phys. Chem. B 112, 1885
共2008兲.

Downloaded 19 Jan 2010 to 150.46.100.222. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

214506-10
15

J. Chem. Phys. 131, 214506 共2009兲

Nakayama, Koga, and Tanaka

H. Lee, J. W. Lee, D. Y. Kim, J. Park, Y. T. Seo, H. Zeng, I. L. Moudrakovski, C. I. Ratcliffe, and J. A. Ripmeester, Nature 共London兲 434, 743
共2005兲.
16
D. Y. Kim and H. Lee, J. Am. Chem. Soc. 127, 9996 共2005兲.
17
R. Anderson, A. Chapoy, and B. Tohidi, Langmuir 23, 3440 共2007兲.
18
S. Alavi, J. A. Ripmeester, and D. D. Klug, J. Chem. Phys. 297, 024507
共2005兲.
19
N. I. Papadimitriou, I. N. Tsimpanogiannis, A. T. Papaionnou, and A. K.
Stubos, J. Phys. Chem. C 112, 10294 共2008兲.
20
N. I. Papadimitriou, I. N. Tsimpanogiannis, C. J. Peters, A. Th. Papaio-annou, and A. K. Stubos, J. Phys. Chem. B 112, 14206 共2008兲.
O. Yamamuro, T. Matsuo, H. Suga, W. I. F. David, R. M. Ibberson, andA. J. Leadbetter, Physica B 213–214, 405 共1995兲.
22
M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112, 8910 共2000兲.
23
Y. Koyama, H. Tanaka, G.T. Gao, and X. C. Zeng, J. Chem. Phys. 121,
7926 共2004兲.
24
A. Bródka and T. W. Zerda, J. Chem. Phys. 104, 6319 共1996兲.
25
I. F. Silvera and V. Goldman, J. Chem. Phys. 69, 4209 共1978兲.
26
I. Ohmine, H. Tanaka, and P. G. Wolynes, J. Chem. Phys. 89, 5852
共1988兲.
21

Downloaded 19 Jan 2010 to 150.46.100.222. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp