J. Phys. Chem. B 2006, 110, 501-506
501
Molecular Dynamics-based Approach to Study the Anisotropic Self-Diffusion of Molecules in Porous Materials with Multiple Cage Types: Application to H2 in Losod Annemieke W. C. van den Berg,† Edwin Flikkema,‡ Sander Lems,§ Stefan T. Bromley,*,| and Jacobus C. Jansen† Ceramic Membrane Centre “The Pore”, Delft UniVersity of Technology, Julianalaan 136, 2628 BL Delft, The Netherlands, Cambridge UniVersity Centre for Computational Chemistry, Department of Chemistry, Lensfield Road, Cambridge CB2 1EW, U. K., Delft Institute for Sustainable Energy, Delft UniVersity of Technology, Department of Chemistry, Julianalaan 136, 2628 BL Delft, The Netherlands, and Department Quimica Fisica and Centre Especial de Recerca en Quimica Teorica, UniVersitat de Barcelona and Parc Cientific, c/Marti i Franques 1, E-08028 Barcelona, Spain ReceiVed: September 6, 2005; In Final Form: NoVember 3, 2005
The anisotropic self-diffusion of molecular hydrogen in the multiple cage clathrasil losod (LOS) is modeled by means of molecular dynamics (MD) simulations of up to 1 µs for the temperature range 900-1200 K while treating the framework as fully flexible. The LOS diffusion tensor is calculated employing an analytical method based on hopping rates. The diffusion in the c-direction of the unit cell is found to be approximately two times more rapid than in the a- and the b-directions, a characteristic of importance for the application of LOS as a membrane. The overall diffusion is based on five different hop types for which the individual hopping rates and diffusion barriers are calculated separately. We show explicitly that the shape and volume of the cages have a significant effect on the hopping rates and further that even small deformations of the circular Si6O6 apertures have a large influence on the energetic barrier for hydrogen diffusion. Compared to the single cage clathrasils dodecasil 3C (MTN) and sodalite (SOD), LOS has a lower diffusion rate. However, from a technical point of view this rate (at 573 K) is still fast enough for LOS to be interesting as a sizeselective membrane or as a hydrogen-adsorption medium.
Introduction When assessing the technical feasibility of a material as a hydrogen-storage medium or as a molecular sieve for the sizeselective separation of H2 from impurities, it is important to know the diffusion rate of H2 through the material as it is this rate which effectively determines the hydrogen uptake and release rates. Many computational studies on gas diffusion and adsorption in zeolites show that the current quality of computational methods is sufficient to obtain valuable insights into the performance of zeolites as adsorbents and membranes.1-5 In previous computational works,5-10 we have assessed the feasibility of sodalite (SOD) and dodecasil 3C (MTN) as potential hydrogen-storage materials and found promising results. These two structures belong to the group of zeolitic sixring clathrates of which 17 members are currently known to exist;11 see Table 1. The six-ring clathrates are porous crystals in which the pores are constituted by rings with six T-atoms (T ) Si, Al, P, Ge, etc.) and six oxygen atoms, which are separated by cages (see for example Figure 1). The accessible diameter of these six-rings is only slightly smaller than the LennardJones diameter of a H2 molecule (∼3 Å) and is relatively flexible.12 These factors, combined with the compliant interaction of the hydrogen molecule with the six-ring, allow for H2 to cross the six-rings with a modest activation barrier. This activated diffusion is possible at relatively elevated temperatures * To whom correspondence should be addressed. E-mail: s.bromley@ qf.ub.es. † Ceramic Membrane Centre “The Pore”. ‡ Cambridge University Centre for Computational Chemistry. § Delft Institute for Sustainable Energy. | Universitat de Barcelona and Parc Cientific.
while at ambient conditions the hydrogen is locked inside the cavities between the six-rings. This phenomenon, known as encapsulation,13 hypothetically enables hydrogen storage at ambient conditions with an easily controllable uptake and release mechanism, either by temperature14,15 or by external force16 to deform larger rings. Furthermore, the energetic cost of such a storage mechanism is low as H2 is stored as a molecule and no dissociation is required. Although they possess these favorable characteristics, of all the six-ring clathrates, the only one that has received much attention in the literature is SOD. Experimentally measured hydrogen-storage capacities in SOD can be found in refs 14 and 15, and calculated capacities are reported in ref 9. Regarding diffusion of gases in SOD, a few papers4,6,7,17-19 report computational studies while no experimental work is currently known. The reasons for this strong focus on SOD are probably its high symmetry (there exists only one type of six-ring and cage type) and its very accessible structure, making the description of diffusion in its interior straightforward and enabling the sampling of a statistically viable number of cage hops in a feasible computational time. Recently, we also reported8 on hydrogen diffusion in MTN, which has the comparable features of a single cage type, single six-ring type, and relatively high symmetry. Table 1, however, shows that, unlike MTN and SOD, most zeolitic clathrates possess pore systems consisting of more than one type of six-ring. It is thus important to study diffusion in multiple ring/cage-type systems in order to understand how such differences affect the diffusion characteristics in these materials. From the list given in Table 1, we have selected the clathrate losod (LOS) as a suitable test case, because it has a reasonable number of distinct six-ring types combined with relatively small
10.1021/jp055033l CCC: $33.50 © 2006 American Chemical Society Published on Web 12/09/2005
502 J. Phys. Chem. B, Vol. 110, No. 1, 2006
van den Berg et al.
TABLE 1: Overview of All 17 Clathrate Structures with Six-Rings as Largest Accessible Openings name
structure code
packinge
number of six-ring types
average number of six-rings per accessible cage
afghanite AlPO4-16 dodecasil 1H franzinite giuseppettite liottite losod linde type N marinellite melanophlogite MCM-61a ZSM-39b nonasil RUB-10c sigma-2 sodalite IM-10d
AFG AST DOH FRA GIU LIO LOS LTN MAR MEP MSO MTN NON RUT SGT SOD UOZ
(4665)3(46617)1 (46612)1(46)1 (51268)1(435663)2(512)3 (46611)1(4665)1(4668)3 (46623)1(4668)2(4665)5 (46611)1(46617)1(4665)4 (4665)1(46611)1 (4668)3(4665)4(4126886)1(476881)6(4662)2 (46617)1(4668)2(4665)3 (512)1(51262)3 (68)1(46620)1(4662)2(4264)3 (512)2(51264)1 (5464)2(4158)2(58612)1 (445462)1(485868)1 (51268)1(4356)2 (4668)1 (46)2(4264)1(410620)1
7 1 2 12 8 8 4 7 8 1 7 1 2 4 2 1 3
8.0 12.0 4.7 8.0 8.0 8.0 8.0 6.5 8.0 2.0 6.3 4.0 6.7 5.0 8.0 8.0 12.0
a MCM ) mobile composition of matter. b ZSM ) zeolite socony mobil. c RUB ) Ruhr University Bochum. d IM ) Institut Franc¸ ais du Pe´troleMulhouse. e Each cage-type has its own description between the brackets, and the ratio between the different cage-types is given by the superscripts outside the brackets. Between brackets the ring-types (normal text) and their occurrence (superscript) in a cage-type are given.20
Figure 1. Schematic representation of the LOS structure in which the shapes of the smaller cage type A, (4665), the larger cage type B, (46611), and the connectivities between the cages are clearly visible. The vertexes are Si-atoms while O-atoms are not shown but lie on the middle of an edge between two vertexes.
cages. In previous work,8 we have shown that a smaller cage corresponds to a larger hopping rate, thus, taking the structure with the smallest cages reduces the computational time needed for a viable number of cage hops, making this simulation possible in a reasonable time. A schematic representation of the LOS structure is given in Figure 1 showing that LOS consists of two types of cages. The smaller cage consists of six four-rings and five six-rings, (4665), which we name type A. The larger cage consists of six fourrings and 11 six-rings, (46611), which we name type B. The ratio between the cages is 1:1, so LOS can be described20 as (4665)1(46611)1. The symmetry of the framework allows classification of the six-ring portals into five distinct types FG () AA, AB, BA, BB1, and BB2). Cage A contains two six-ring apertures of type AA (each shared with another A-cage) and three of type AB (each shared with a B-cage). Cage B contains three six-ring exits of type BA (each shared with an A-cage), two of type BB1 (each shared with a neighboring B-cage in the c-direction of the unit cellsthe vertical direction in Figure 1),
and six of type BB2 (each shared with a neighboring B-cage in the ab-plane of the unit cellsthe horizontal direction in Figure 1). Since the B-cages in the ab-plane are shifted compared with each other, the vector connecting their centers also has a c-component; LOS contains thus straight and sinusoidal pores. In LOS, there exist 16 distinct vectors for hopping from one cage center to a neighboring cage center, which we call the hop vectors: AA1, AA2, AB1, AB2, AB3, BA1, BA2, BA3, BB11, BB12, BB21, BB22, BB23, BB24, BB25, and BB26. In this study, the anisotropic self-diffusion of molecular hydrogen in the all-silica LOS framework is modeled by means of force-field-based molecular dynamics (MD) simulations for the temperature range 900-1200 K. Since our previous study12 showed that framework flexibility is a very important factor in determining the energetic barrier for diffusion through six-rings and, thus, for the diffusion rate, the LOS framework was allowed to be fully flexible in all calculations. Simulations of up to 1 µs were performed to obtain enough statistical data to estimate the hopping rates for all five nearest-neighbor hop types (i.e., through ring portals AA, AB, BA, BB1, and BB2); based on these rates, a diffusion tensor could be derived. The obtained results are compared to the previous results6,8 for SOD and MTN, and the observed differences are rationalized in terms of six-ring shape, cage shape, and cage volume. The results are discussed with respect to their importance in understanding diffusion in complex multicage/ring porous systems and the application of such materials as selective membranes and hydrogen-storage media. Computational Methodology The force field used to model the LOS clathrasil framework utilizes a Buckingham potential form for Si-O and O-O interactions, a harmonic O-Si-O three-body term, and a spring constant to define a negative shell around a positive oxygen core.21 The interactions between the hydrogen molecule and the clathrasil framework atoms employ specifically parametrized Lennard-Jones potentials.22 The combination of these potentials is described and validated in previous studies,9,12 and it was demonstrated by means of DFT calculations that the H2 molecule behaves as an effectively neutral and centrosymmetric particle inside the clathrasil even in the largely repulsive environment of the six-ring center.12 The cutoff for all nonCoulombic interatomic potentials used was set to 10 Å.
Anisotropic Self-Diffusion of H2 in LOS
J. Phys. Chem. B, Vol. 110, No. 1, 2006 503
LOS is modeled as a 2 × 2 × 2 unit cell system with periodic boundary conditions, Si192O384: 16 (4665) cages and 16 (46611) cages, employing no symmetry constraints, and with one hydrogen molecule inside. The LOS geometry is taken from the International Zeolite Association (IZA) database.11 Since this geometry is based on experimental measurements on the silica-alumina form of LOS filled with sodium cations and water, the pure silicate framework structure was first energy minimized at a constant pressure of 1 bar with the BFGS optimization algorithm,23 using the general utility lattice program (GULP),24 in order to obtain a geometry that gives a better representation of the empty all-silica structure. The optimized supercell has the following dimensions: a, b, c ) 24.652, 24.652, 20.390 Å and R, β, γ ) 90, 90, 120°. All MD simulations were performed using the computer code DLPOLY.25 The starting position for the hydrogen molecule was always located near the center of an arbitrary A- or B-cage, and an equilibration time of 200 ps was employed. Simulations were performed at 900, 1000, 1100, and 1200 K with the NVT Evans ensemble26 using a time step of 0.001 ps. These temperatures are sufficiently high to assume negligible quantum effects on the H2 molecule.7 Simulations at lower temperatures could not be performed due to prohibitively demanding computational simulation times needed to obtain sufficient hop statistics. For each temperature, the simulation time was always chosen long enough to include at least 200 hops in order to ensure a reasonably low statistical error in the hopping rate. The MD-output data was analyzed with an in-house hop counting program.7,8 All hops that happened within a few picoseconds of each other were checked manually for recrossings and subsequent hops. Theory Hydrogen diffusion through the LOS framework is only possible via the six-rings (see Figure 1). Due to the nature of the framework, the diffusion of hydrogen in this structure will be anisotropic and, thus, the diffusion coefficient will be a tensor27 rather than a scalar as shown in eq 1.
B J ) -D∇φ
(
Dxx Dyx Dzx with D ) Dxy Dyy Dzy Dxz Dyz Dzz
)
(1)
〈(R(t) - R(0))(β(t) - β(0))〉 2t
(2)
In this equation, the numerator gives the mean-square displacement in the Rβ-direction and t is time. However, generating MD trajectories sufficiently long to be able to apply eq 2 directly would be computationally too demanding. Instead, the diffusion tensor can be approximated with eq 3, if the hopping events are uncorrelated.
DRβ )
1
R β bFGi ∑ NFGibFGi
2t FGi
hopping possibility no.
hop vector FGi
x bFGi
y bFGi
z bFGi
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
AA1 AA2 AB1 AB2 AB3 BA1 BA2 BA3 BB11 BB12 BB21 BB22 BB23 BB24 BB25 BB26
0 0 7.2 -3.6 -3.6 7.2 -3.6 -3.6 0 0 7.2 -3.6 -3.6 7.2 -3.6 -3.6
0 0 0 6.2 -6.2 0 6.2 -6.2 0 0 0 6.2 -6.2 0 6.2 -6.2
5.1 -5.1 0 0 0 0 0 0 10.2 -10.2 5.1 5.1 5.1 -5.1 -5.1 -5.1
Here, t is the total simulated time starting at zero when the first hop occurs, NFGi is the number of hops with vector-type FGi observed during the simulation (hop at t ) 0 not included), R β and bFGi and bFGi are the components R and β of the corresponding vector. For LOS, the summation has to be performed over the 16 different hop vectors already described in the Introduction. The coordinates system of the framework is chosen such that the z-direction coincides with the c-direction of the unit cell and the xy-plane coincides with the ab-plane. The resulting hop vectors are given in Table 2 with the origin of each vector placed in the center of the cage from where it starts. Since self-diffusion has no driving force, NFGi can be assumed to be equal for all hop vectors of a certain hop type, that is, vectors for which the FG is equal (NAA1 ) NAA2, NBB11 ) NBB12 * NBB21, etc.). Combining eq 3 with Table 2 shows that all nondiagonal elements of the diffusion tensor are zero and that Dxx ) Dyy * Dzz. The physical explanation for this simple diffusion tensor can be found in the high symmetry of the LOS structure. The Dzz represents the diffusion in the c-direction of the LOS unit cell (Dc), and the Dxx and Dyy both represent the diffusion in the ab-plane (Da ) Db, in short Dab). An alternative form of eq 3 expressing more specific information about the diffusion behavior of H2 in LOS is given by eq 4.
DRβ )
In this equation, B J is the flux, D is the diffusion tensor, ∇φ is the concentration gradient, and DRβ (with R and β both being x, y, or z) is the diffusion coefficient for the Rβ-component. The diffusion tensor is related to the tensor of the mean-square displacements as shown in eq 2.
DRβ )
TABLE 2: Vector Components of the 16 Distinct Hop Possibilities in LOS
(3)
1
R β pFkFGibFGi bFGi ∑ 2 FGi
(4)
Here, pF is the probability of hydrogen being present in a cage type F (with F being A or B for LOS) and kFGi is the hopping rate (hops/s) of the FGi hops. The pF values for LOS can be extracted directly from the MD data (assuming that the dataset is large enough to be representative for the complete diffusion behavior, vide infra) by determining the fraction of time that the H2 molecule is present in cage type A and in cage type B
pA )
tA t
pB )
tB t
and (5)
With tA being the total time that the H2 molecule is present in cage type A and tB being the total time that H2 is present in
504 J. Phys. Chem. B, Vol. 110, No. 1, 2006
van den Berg et al.
TABLE 3: Main Results of the MD-simulations: Total Simulation Times, Number of Parallel Simulations, Total Number of Observed Hops, and Total Number of Subsequent Hops
TABLE 4: Specification of the Hops Observed during the MD Simulationsa
number of parallel total number total number of total simulation simulations of hops Nb subsequent hopsc T (K) time t (ns)a 900 1000 1100 1200
1000 447 298 148
25 15 10 10
201 208 202 211
2 5 2 0
Excluding equilibration time. b Including the first hops and the subsequent hops, rapid recrossings were not observed during the simulations. c All subsequent hops were of the type A f A f A, except for one at 1000 K which was B f B f A.
cage type B. The hopping rates can be calculated from the average residence time prior to a hop τ(F) and the relative frequencies of hopping through the different ring-types as determined from the MD data
1
NFGi
τ(F)
NFGi ∑ Gi
-∆EFGi/kBT
(7)
As an illustration of the technological relevance of the results, the calculated diffusion constants are interpreted in terms of a hydrogen uptake time for a spherical crystal. The hydrogen uptake as a function of time can be estimated by means of eq 8 where Dtrans represents the transport diffusion coefficient. At low loadings, Dtrans equals the self-diffusion coefficient D.27 Furthermore, d is the diameter of the spherical crystal, 〈CH2〉 is the average hydrogen concentration in the crystal, and CH∞2 is the hydrogen-saturation concentration.
〈CH2〉 CH∞2
)1-
6
∞
∑
1
π2 n)1 n2
(
exp -n
2
NAB
NBA
NBB1
NBB2
900 1000 1100 1200
72 83 34 55
12 13 13 16
4 7 11 12
9 7 11 9
75 73 119 109
Excluding all first hops and all subsequent hops.
TABLE 5: Average H2 Residence Times Per Cage Type and the Probabilities of Finding the H2 Molecule in Cage Type A or Ba T (K)
τ (A) (ps)
τ (B) (ps)
pA
pB
900 1000 1100 1200
2328.3 1131.2 580.9 288.7
4406.9 2511.3 1384.9 860.7
0.34 0.33 0.12 0.15
0.66 0.67 0.88 0.85
a Excluding (the time prior to) the first hop of each simulation and all subsequent hops, see Table 3.
(6)
Here, ∑Gi NFGi is the total number of hops departing from cage type F. Since the residence time in cage F does not depend on the following hop type, it can simply be calculated with τ(F) ) tA/∑Gi NFGi. Previously, we already assumed that NFGi is equal for all hop vectors of one hop type FG, therefore, the same holds for kFGi. To rationalize the temperature dependence of the hopping rates, an Arrhenius-type representation in terms of an effective energetic barrier (∆EFGi) and a preexponential hopping rate (AFGi) was used for each FGi hop
kFGi ) AFGie
NAA
a
a
kFGi )
T (K)
)
4π2Dtranst d2
(8)
Since the diffusion coefficient in eq 8 is a scalar, it is only possible to determine limiting values for the hydrogen uptake time in LOS by using the distinct elements of the diffusion tensor, Dab and Dc. Results and Discussion The data obtained from the MD simulations are summarized in Tables 3-5. A limited number of correlated events were observed in the MD trajectories in the form of subsequent hops occurring within 0.5 ps, but their statistical significance is found to be too low to warrant inclusion of these correlated events in the model (see Table 3).
Figure 2. Arrhenius plots of the hydrogen hopping rates for each sixring type.
The hopping rates are calculated with eq 6 under the assumption that the observed hops of type FG (Table 4) are equally distributed over the i hop vector types (NAAi ) NAA/2, NABi ) NAB/3, NBAi ) NBA/3, NBB1i ) NBB1/2, NBB2i ) NBB2/6) and are given in Table 6 together with their Poisson errors. For the calculation of these statistical errors, we assumed no correlation in the data points, an assumption previously shown to be correct for the H2 self-diffusion in the clathrate structures SOD and MTN.6,8 The error in the hopping rates then equals 100%/xNFG. Since the AA and BB2 hop types occur the most frequently, their errors are relatively low, while the errors in the less frequent events AB, BA, and BB1 are rather high. However, to obtain smaller errors for these events we would have to run simulations in the order of a millisecond which, with the current computer speed, would take years. Furthermore, the overall diffusion tensor mainly depends on the frequent hopping events and their (smaller) errors. The volumes of the A- and B-cage in LOS are 183 and 488 Å3, respectively. In previous work,8 using transition-state theory, we have shown that a smaller cage volume leads to a higher hopping rate out of that cage. This trend is also observed in Table 6 for the cages in LOS. Figure 2 shows the Arrhenius plots of the hopping rates and a linear fit through the points (see also eqs 9a-e). Combining these fits with eq 7 gives the values for the energetic barrier for each hop type (i.e., each six-ring type) and the preexponential rate factors (Table 7). As can be expected from the errors in the hopping rates and the exploding effect of the exponential relation, the errors in the pre-factors are very large. Fortunately, the values for the energetic barriers are more reliable and allow
Anisotropic Self-Diffusion of H2 in LOS
J. Phys. Chem. B, Vol. 110, No. 1, 2006 505
TABLE 6: Hopping Rates for Each FGi Hop and Their Statistical Errors for the Temperature Range 900-1200 K T (K)
kAAi/1 × 107 (hops/s)
err. (%)
kABi/1 × 107 (hops/s)
err. (%)
kBAi/1 × 107 (hops/s)
err. (%)
kBB1i/1 × 107 (hops/s)
err. (%)
kBB2i/1 × 107 (hops/s)
err. (%)
900 1000 1100 1200
18.41 38.21 62.27 134.16
11 11 17 13
2.05 3.99 15.87 26.02
28 26 28 25
0.34 1.07 1.88 3.58
41 33 29 29
1.16 1.60 2.82 4.02
29 35 29 33
3.22 5.57 10.16 16.24
11 11 9 9
TABLE 7: Pre-exponential Rate Factors and the Energetic Barriers for Hopping through the Five Six-Ring Exit Types Together with Their Mean Least-Square Errors A (hops/s)
error (hops/s)
T (K)
Dab (m2/s)
Dc (m2/s)
DSOD (m2/s)
DMTN (m2/s)
3.87 × 1011 8.75 × 1011 3.56 × 1010 1.80 × 109 2.11 × 1010
2.28 × 1011 10.9 × 1011 1.98 × 1010 0.86 × 109 0.56 × 1010
900 1000 1100 1200
2.0 × 10-11 3.6 × 10-11 8.1 × 10-11 1.3 × 10-10
4.1 × 10-11 7.3 × 10-11 1.1 × 10-10 2.0 × 10-10
6.8 × 10-10 8.8 × 10-10 1.2 × 10-9 1.8 × 10-9
1.5 × 10-10 2.6 × 10-10 2.9 × 10-10 4.9 × 10-10
hop type ∆E (kJ/mol) error (kJ/mol) AA AB BA BB1 BB2
57.5 80.7 68.6 38.3 48.8
5.1 10.7 4.8 4.1 2.3
TABLE 9: Distinct Elements of the Diffusion Tensor of Hydrogen in LOS and the Diffusion Coefficients in SOD and MTN for the Temperature Range 900-1200 K
TABLE 8: Distances (Å) of the Opposing O-atoms in the Six-Rings of SOD, MTN, and LOS
TABLE 10: Ratios between the Uptake Times of LOS (both directions), MTN, and SOD for the Temperature Range 900-1200 K
SOD
MTN
AA
AB and BA
BB1
BB2
T (K)
LOS (ab)
LOS (c)
SOD
MTN
5.05 5.05 5.06
5.15 5.15 5.20
5.07 5.07 5.07
5.00 5.00 4.75
5.08 5.08 5.08
4.98 4.98 5.30
900 1000 1100 1200
34 24 15 14
17 12 11 9
1 1 1 1
5 3 4 4
for a comparison with the barriers found for SOD6 (32.9 ( 1.5 kJ/mol) and MTN8 (33.8 ( 2.8 kJ/mol). The linear fits in Figure 2 are given by eqs 9a-e.
ln(kAA) ) -6920.1/T + 26.681
with R2 ) 0.98 (9a)
ln(kAB) ) -9704.5/T + 27.498
with R2 ) 0.96 (9b)
ln(kBA) ) -8255.7/T + 24.297
with R2 ) 0.99 (9c)
ln(kBB1) ) -4601.8/T + 21.313
with R2 ) 0.97
(9d)
ln(kBB2) ) -5870.3/T + 23.772
with R2 ) 0.99
(9e)
The energetic barriers of the six-rings in LOS are all higher than those of MTN and SOD. In part, this can be explained by the higher adsorption energy of H2 in the smaller A-cage caused by the larger number of framework atoms that can attract the H2 simultaneously. This effect increases the energy needed to escape from the A-cage. For the same reason, the barrier for the hop from cage A to cage B is higher than that for the reverse hop. Additionally, the elongated shape of the B-cage also has a small increasing effect on the adsorption energy although weaker than that observed for the A-cage. A second important effect is illustrated by Table 8 showing the distances between opposing O-atoms in the six-rings of the optimized empty clathrasil structures. These distances indicate how circular the six-ring is. The SOD, MTN, AA, and BB1 rings have three similar values and are thus completely circular, while the AB and BB2 rings have one shorter and two longer distances indicating a flattened ring (see also Figure 1). Since a circular form is the most favorable shape for passage (having the most symmetrically open form for a centrosymmeric molecule), the deformed six-rings will have a higher energetic barrier for H2 passage. This is consistent with the barrier values in Table 7. The diffusion tensor is calculated with eq 3 for the temperature range 900-1200 K (see Table 9). The diffusion coefficients of SOD and MTN are also given for comparison.
Table 9 shows that the diffusion in the c-direction, that is, via the AA, BB1, and BB2 hops, is approximately two times faster than that in the ab-plane, that is, via the AB, BA, and BB2 hops. This is mainly a result of the lower energetic barrier of AA and BB1 compared to AB and BA. Furthermore, the overall diffusion in LOS is slower than that in SOD and MTN, due mainly to the higher average energetic barrier for diffusion in LOS. Additionally, there may be a small influence from the earlier mentioned difference in adsorption energy and hopping rate caused by different shapes and volumes of the LOS cages (A ) 183 Å3 and flattened, B ) 488 Å3 and elongated) compared to the SOD (337 Å3 and spherical) and MTN cage (461 Å3 and spherical). The resulting relative uptakes (eq 8) in the three clathrasils are compared in Table 10 showing that uptake in LOS is slower than that in SOD and MTN and that the difference increases with decreasing temperature. The uptake times in SOD and MTN at a technically reasonable temperature of 573 K are 1.2 and 5.0 s.6,8 Extrapolating the values in Table 9 shows that the uptake of LOS at 573 K is about 40-220 times slower than that for SOD, that is, an uptake time in the order of minutes. Still indicating, however, that from a technical point of view the LOS structure is (at 573 K) also an interesting candidate for applications requiring a fast H2 uptake such as hydrogen adsorption (for storage purposes) and separation/purification via a size-selective membrane. The difference in diffusion rate and uptake time for the different directions in LOS is a potentially very important fact for the application of LOS as a membrane to be used, for example, for the purification of hydrogen. Our calculations predict that a LOS membrane grown with the c-direction perpendicular to the support material gives the quickest permeance and is thus the most desirable form. Such utilization of zeolitic orientation, demonstrating that the principle has a strong influence on technological performance, has already been shown experimentally for an MFI membrane for organic vapor separation.28 Conclusions The anisotropic self-diffusion tensor of molecular hydrogen in LOS is calculated for the temperature range 900-1200 K
506 J. Phys. Chem. B, Vol. 110, No. 1, 2006 by means of molecular dynamics simulations allowing for full framework flexibility. The energetic barriers for the five different hop types of H2 in LOS are all higher than the barriers found for the six-rings of SOD and MTN due to the slight deformation in two of the four LOS six-ring types and the higher adsorption energy in the small (4665) LOS cages. Due to these higher barriers, LOS has a considerably lower H2 diffusion rate than MTN and SOD. The corresponding H2 uptake time in LOS at 573 K is in the order of several minutes, meaning that the LOS uptake rate is still fast enough at this temperature to make the framework technically interesting as a hydrogen-adsorption medium or as a size-selective membrane. The diffusion in the c-direction is also found to be approximately two times as fast as that in the ab-plane indicating that a LOS membrane with its c-orientation perpendicular to the support material will show the best technological performance in terms of permeance speed for H2 separation. Acknowledgment. The authors thank Dr. Martijn A. Zwijnenburg for creating Figure 1 and Dr. Olaf Delgado-Friedrichs for providing his 3D-Tiler visualization software. E.F. wishes to acknowledge that this research was supported by a Marie Curie Intra-European Fellowship (MEIF-CT-2005-011551) within the 6th European Community Framework Programme. References and Notes (1) Kar, S.; Chakravarty, C. Chem. Phys. Lett. 2002, 352, 294. (2) Jia, W.; Murad, S. J. Chem. Phys. 2005, 122, 234708. (3) Akten, E. D.; Siriwardane, R.; Sholl, D. S. Energy Fuels 2003, 17, 977. (4) De Luca, G.; Pullumbi, P.; Barbieri, G.; Fama, A. D.; Bernardo, P.; Drioli, E. Sep. Purif. Technol. 2004, 36, 215. (5) van den Berg, A. W. C.; Bromley, S. T.; Wojdel, J. C.; Jansen, J. C. Microporous Mesoporous Mater. 2005, 87, 235. (6) van den Berg, A. W. C.; Bromley, S. T.; Flikkema, E.; Wojdel, J.; Maschmeyer, T.; Jansen, J. C. J. Chem. Phys. 2004, 120, 10285.
van den Berg et al. (7) van den Berg, A. W. C.; Bromley, S. T.; Flikkema, E.; Jansen, J. C. J. Chem. Phys. 2004, 121, 10209. (8) van den Berg, A. W. C.; Flikkema, E.; Bromley, S. T.; Jansen, J. C. J. Chem. Phys. 2005, 122, 204710. (9) van den Berg, A. W. C.; Bromley, S. T.; Jansen, J. C. Microporous Mesoporous Mater. 2005, 78, 63. (10) van den Berg, A. W. C.; Bromley, S. T.; Wojdel, J. C.; Jansen, J. C. Phys. ReV. B 2005, 72, 155428 (11) Baerlocher, C.; Meier, W. M.; Olson, D. H. Atlas of Zeolite Framework Types; Elsevier: Amsterdam, The Netherlands, 2001 (updates available at http://www.iza-structure.org/). (12) van den Berg, A. W. C.; Bromley, S. T.; Ramsahye, N.; Maschmeyer, T. J. Phys Chem. B 2004, 108, 5088. (13) Breck, D. W. Zeolites molecular sieVes, structure, chemistry and use; John Wiley & Sons: New York, 1974; pp 623-628. (14) Weitkamp, J.; Fritz, M.; Ernst, S. Int. J. Hydrogen Energy 1995, 20, 967. (15) Weitkamp, J.; Schmid, D.; Fritz, M.; Cubero, F.; Ernst, S. Wasserst. Energietraeger, Kolloq. 1994 Sondersforschungsbereichs 270; University Stuttgart: Stuttgart, Germany, 1994; pp 287-300. (16) Maschmeyer, T.; Jansen, J. C. EP 2002, 1188477. (17) Kopelevich, D. I.; Chang, H.-C. J. Chem. Phys. 2001, 115, 9519. (18) Murphy, M. J.; Voth, G. A.; Bug, A. L. R. J. Phys. Chem. B 1997, 101, 491. (19) Nada, R.; Nicholas, J. B.; McCarthy, M. I.; Hess, A. C. Int. J. Quantum Chem. 1996, 60, 809. (20) Zwijnenburg, M. A.; Bromley, S. T.; Foster, M. D.; Bell, R. G.; Delgado-Friedrichs, O.; Jansen, J. C.; Maschmeyer, T. Chem. Mater. 2004, 16, 3809. (21) Sanders, M. J.; Leslie, M.; Catlow, C. R. A. J. Chem. Soc., Chem. Commun. 1984, 19, 1271. (22) Watanabe, K.; Austin, N.; Stapleton, M. R. Mol. Simul. 1995, 15, 197. (23) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes; Cambridge University Press: Cambridge, U.K., 1992. (24) Gale, J. D. J. Chem. Soc., Faraday Trans. 1997, 93, 629. (25) Forester, T. R.; Smith, W. The DL_POLY_2 User Manual; CCLRC Daresbury Laboratory: Daresbury, U.K., 1997. (26) Evans, D. J.; Morriss, G. P. Comput. Phys. Rep. 1984, 1, 297. (27) Karger J.; Ruthven, D. M. Diffusion in Zeolites and Other Microporous Solids; J. Wiley & Sons: New York, 1992. (28) Lai, Z.; Bonilla, G.; Diaz, I.; Nery, J. G.; Sujaoti, K.; Amat, M. A.; Kokkoli, E.; Terasaki, O.; Thompson, R. W.; Tsapatsis, M.; Vlachos, D. G. Science 2003, 300, 456.
