Dynamics of n-Butaneâ Methane Mixtures in Silicalite, Using

J. Phys. Chem. B 2000, 104, 5541-5552

5541

Dynamics of n-Butane-Methane Mixtures in Silicalite, Using Quasielastic Neutron Scattering and Molecular Dynamics Simulations Leonidas N. Gergidis,† Doros N. Theodorou,† and Herve´ Jobic*,‡ Department of Chemical Engineering, UniVersity of Patras, and Institute of Chemical Engineering and High-Temperature Chemical Processes, GR 26500 Patras, Greece, and Institut de Recherches sur la Catalyse, CNRS, 2 AVenue Albert Einstein, 69626 Villeurbanne, France ReceiVed: January 4, 2000; In Final Form: April 3, 2000

The transport of n-butane-methane mixtures in the zeolite silicalite was studied using molecular dynamics simulations and quasielastic neutron scattering experiments over a range of loadings and compositions at 200 K. Self-diffusivities are seen to decrease monotonically with loading of either species. Self-diffusivity values calculated from the dynamics simulations are in excellent agreement with the experimental measurements from quasielastic neutron scattering if one takes into account the errors associated with both techniques. We also studied the detailed dynamical behavior of the system. The dependence on the wave vector of the halfwidth at half-maximum of the incoherent dynamic structure factor, describing self-correlations in the motion of sorbate molecules, is indicative of a jump diffusion process. By monitoring and analyzing the molecular motion in the simulation, we confirmed that diffusion takes place through successive jumps between the interiors of adjacent channel segments. Precise and quantitative calculations mapping the MD trajectories onto a coarse-grained jump model reveal mechanistic aspects of the motion. Distributions of jump lengths and rate constants are accumulated for the various jump types executed by each sorbate species. Jump lengths are widely distributed between 0 and 15 Å, the mean jump length being a decreasing function of the loading. The mean time between jumps is actually smaller at higher occupancies, because there short jumps prevail, occurring back and forth in a highly correlated fashion.

1. Introduction Zeolites are widely used as catalysts in the petroleum and petrochemical industries in a number of processes, such as xylene isomerization, catalytic dewaxing, and conversion of methanol to gasoline. They are also used in industrial separations and have assumed an important role in water cleaning and softening as selective adsorbents.1,2 In recent years, molecular simulation techniques, and especially molecular dynamics, have proven to be a powerful tool for predicting zeolite structureproperty relationships and elucidating the detailed molecular processes shaping these relationships.3 Furthermore, molecular dynamics is tailored for calculating transport properties, such as self-diffusivities. Past simulation studies for a variety of zeolite-sorbate systems are in very good agreement with microscopic techniques such as pulsed field-gradient NMR (PFG-NMR) and quasielastic neutron scattering (QENS). The disagreement often observed with more macroscopic experimental techniques used for the calculation of transport properties, for example, membrane permeation and uptake rate experiments, usually originates from the fact that these macroscopic measurements are rate-limited by processes other than intracrystalline diffusion. Although the good agreement between the three microscopic techniques, PFG-NMR, QENS, and MD, is relatively well established for pure sorbates,4,5 relatively few studies exist in the literature presenting direct comparisons between those techniques for mixed sorbates. A direct com* To whom correspondence should be addressed. E-mail: jobic@ catalyse.univ-lyon1.fr. † University of Patras. ‡ Institut de Recherches sur la Catalyse.

parison between MD simulations and experimental self-diffusivities from PFG-NMR for a binary mixture of methanetetrafluoromethane was made by Snurr and Ka¨rger.6 The reported self-diffusivities from the simulations were in good agreement with the experimental measurements. These authors observed that the diffusivities for both components, at constant total loading, decreased as the fraction of the larger and less mobile CF4 increased. Also very recently, Jost et al.7,8 have reported MD and pulsed field gradient NMR studies for the diffusion of methane-xenon mixtures over a wide range of loadings and compositions. This study showed that xenon diffusivity was unaffected by the presence of methane molecules at various concentrations. By varying the concentration of xenon, a monotonic decrease of the diffusivity of both species was observed. Here, we present a first attempt to compare results obtained from molecular dynamics simulations and quasielastic neutron scattering experiments on a binary mixture of sorbates in a zeolite and to elucidate the nature of intracrystalline motion probed by the two techniques. The present investigation focuses on mixtures of n-butane and methane in silicalite. This system was chosen because of its practical applications and because its experimental and computational study is both worthwhile and feasible. Previous MD studies of the same system at 300 K have revealed interesting behavior.9 Self-diffusivities were seen to decrease monotonically with the loading of either species. Raising the loading of n-butane from 2 to 9 molecules per unit cell caused the diffusivity of methane to drop by a factor of 60. The spatial distribution of molecules of the two coadsorbed species was investigated, showing that, at high occupancies, n-butane molecules are forced to populate pref-

10.1021/jp0000073 CCC: $19.00 © 2000 American Chemical Society Published on Web 05/23/2000

5542 J. Phys. Chem. B, Vol. 104, No. 23, 2000

Gergidis et al.

erentially the gauche conformation. An anomalous diffusion regime was identified for both species at higher loadings. Anomalous effects were more pronounced for methane than for n-butane in all three directions, but most strongly in the z direction, along which no direct channel pathway exists. Crossover to normal “Fickian” diffusion was found to occur at times on the order of nanoseconds. Visualization of trajectories from the dynamic simulations revealed a jumplike character of intracrystalline motion. Exploring the nature of these jumps is one of the major objectives of the present experimental and simulation study. We will investigate how the motion of methane is affected by coadsorbed n-butane at 200 K, over a variety of intracrystalline occupancies.

We performed molecular dynamics simulations in the canonical statistical ensemble (NVT) using the Nose´-Hoover12 thermostat. As in our previous work,9 the n-butane molecules were represented in terms of four united-atom methyl and methylene groups (interaction sites). Methane was modeled as a single interaction site. For the n-butane molecules, we assigned the correct deuterated masses of 18 g/mol and 16 g/mol to methyl and methylene, respectively. Details of the simulation algorithm can be found elsewhere.9 Molecular dynamics simulations conducted in this investigation had a duration ranging from 12 to 20 ns. Production runs were carried out on SGI Origin 200 and DIGITAL 433au workstations.

2. Neutron Scattering Experiments, Molecular Dynamics Simulations The experiments were performed at the Institut LaueLangevin, Grenoble. We used the time-of-flight (TOF) spectrometer IN5 with an incident wavelength of λ ) 9 Å, giving an elastic resolution of ca. 18 µeV (full width at half-maximum) which was fairly constant over the whole range of scattering angles (10° < 2 θ < 125°). Because of the large incoherent cross section of hydrogen, the scattering from methane and the methyls and methylenes of n-butane is largely incoherent. For probing the motion of one of the components, it is necessery to turn off the signal from the other. This can be established by deuterating one of the two molecules, because deuterium has a much smaller incoherent cross section. In our case, we used perdeuterated n-butane (99.7% isotopic enrichment), and background was recorded with the zeolite loaded with n-butane-d10 at various concentrations. The coherent scattering from the zeolite could be avoided as much as possible by taking selected groupings of the detectors between the Bragg peaks of the zeolite. The Bragg peaks of the zeolite were determined on the spectrometer by comparing the intensities of the elastic peaks at each angular position with those obtained from a standard vanadium plate (vanadium is an almost totally incoherent scatterer and is used to calibrate the detector’s efficiencies and to measure the elastic resolution). Experiments took place at a temperature of 200 K in order to reach high loadings at reasonable pressures (P < 1 atm). The ZSM-5 sample is the same as in previous QENS studies on alkane diffusion (C1C8).10,11 It has a Si/Al ratio of 36 and a good crystallinity with a low density of silanol groups. The zeolite was activated by heating to 770 K under flowing oxygen. After cooling, the sample was pumped to 10-4 Pa while heating again up to 770 K. The zeolite (7 g) was transferred inside a glovebox into a slab-shaped aluminum container (diameter 50 mm). The container was connected to a gas inlet system, which allowed sorption in situ, even at low temperature. The adsorbed amounts of methane and n-butane-d10 were determined volumetrically during the neutron experiment. When n-butane was adsorbed first at 200 K (with a loading of 4 molecules per unit cell), the equilibrium pressure over the sample was very small (less than 1 mbar), indicating that n-butane is strongly adsorbed and that no molecules are present in the gas phase. Increasing concentrations of methane were then adsorbed onto the sample at the same temperature, and we assume that n-butane does not desorb upon methane adsorption. On the other hand, when methane was adsorbed first, the equilibrium pressure rose progressively upon n-butane adsorption (e.g., from 20 mbar to 28 mbar after adsorption of 2.5 n-butane molecules per unit cell). This shows that a small fraction of methane molecules is driven out of the zeolite, so that the real concentration of methane is slightly less than the quoted value.

3. Theory To connect the “theoretical” information, extracted from the MD simulation trajectories, with the experimental measurements, it is necessary to calculate the self-part of the van Hove correlation function using13

Gs(r,t) )

1 N

∑i ∫ 〈δ(r - r′ + Ri(0))δ(r′ - Ri(t))〉 d3r′

(1)

The angular brackets denote a statistical or thermal average, and the sum is over all particles of a given type (e.g., methane molecules). The self-part of the van Hove correlation function represents the conditional probability density for finding a particle at time t at a distance r from the position it occupied at time t ) 0. Since Gs(r,t) denotes a probability density, its integral over volume equals 1:

∫Gs(r,t) d3r ) 1

(2)

As a next step, we have calculated the intermediate incoherent scattering function Is(Q,t). This is defined as the spatial Fourier transform of the van Hove autocorrelation function:

Is(Q,t) )

∫Gs(r,t)exp(iQ‚r)d3r

(3)

In the QENS experiment, pQ is the neutron momentum transfer. It is defined by Q ) k - k0 where k and k0 are the scattered and incident wave vectors, respectively. Similarly, pω ) E - E0 is the neutron energy transfer, with ω denoting the angular velocity. In the hydrodynamic limit (Q w 0, t w ∞) where Fick’s law describes diffusive motion, the orientationally averaged Is(Q,t) becomes equal to exp(-DsQ2 t) with Ds the self-diffusivity. In a next step, we calculated the incoherent (self) scattering function (or incoherent dynamic structure factor) Sinc(Q,ω), defined as the time Fourier transform of the intermediate scattering function

Sinc(Q,ω) )

1 2π

∫ Is(Q,t)exp(-iωt) dt

(4)

Sinc(Q,ω) is proportional to the measured scattering intensity in an incoherent neutron scattering experiment. The intensity scattered by the sample is related to self-correlations in space and time of the hydrogen atoms under the effect of the different kinds of molecular motions. The total scattering function is usually set equal to a convolution product of the individual scattering laws, corresponding to the different motions of translation, rotation, and vibration:

Sinc(Q,ω) ) STinc(Q,ω)XSRinc(Q,ω)XSVinc(Q,ω)

(5)

n-Butane-Methane Mixtures in Silicalite

J. Phys. Chem. B, Vol. 104, No. 23, 2000 5543

For the energy range studied here ((1 meV), the vibrational term in eq 5 consists only of a frequency-independent DebyeWaller factor, which has the effect of decreasing the total intensity of the spectra for increasing momentum transfer values. The rotational term of methane is well modeled by an isotropic rotational diffusion.14 The broadening due to rotation is small at small momentum transfers (its relative intensity is 5% at 0.35 Å-1), but this contribution cannot be neglected at higher Q values. At very small momentum transfer (Q ≈ 0.1 Å-1), translational motion obeys Fick’s law, and the broadening of the shape of the orientationally averaged Sinc(Q,ω) (half-width at half-maximum or HWHM) is proportional to DsQ2, according to the following expression for the dynamic structure factor:

DsQ2 1 Sinc(Q,ω) ) π ω2 + (D Q2)2

(6)

s

When the momentum transfer increases, the profile of the translational scattering function is still Lorentzian, but the HWHM varies in a more complicated way. Here, only the translational motion will be considered. 4. Jump Models Whereas Fickian diffusion is obtained at small Q values, jump diffusion implies that at large Q (short distances), a deviation will be seen from the straight-line dependence of the HWHM on Q2 expected from eq 6. The interpetation of the QENS spectra at larger Q values requires a model that contains as parameters the characteristic times and lengths of the elementary steps. One such model is the model proposed by Chudley and Elliott (CE).13,15 It is based on the hypothesis that a particle (here, a sorbate molecule) vibrates on a given site during time τ, and, after this time, it jumps to another site situated at a distance d, with the additional requirement that the time taken for the jump is much shorter than the residence time. Under the above hypotheses, the incoherent scattering function is still a Lorentzian function:

Sinc(Q,ω) )

∆ω(Q) 1 2 π ω + (∆ω(Q)2)

(7)

If the jumps have a fixed length d, and if they can occur in any direction, the HWHM is then

∆ω(Q) )

[

]

sin(Qd) 1 1τ (Qd)

(8)

The limit of eq 8 at small Q values is

∆ω(Q) )

Q2d2 6τ

) DsQ2

(9)

Another model proposed by Jobic16 is based on an isotropic jump length distribution of the form

F(r) )

(

)

(r - d0)2 r exp d0r0(2π)1/2 2r02

(10)

where d0 is the distance between two sites and r0 is a measure of the delocalization of the molecule on its site. The mean square jump length corresponding to this distribution is

〈r2〉 )

∫0∞ r2 F(r) dr ) d02 + 3r02

(11)

and the following analytical expression is obtained for the HWHM of the Lorentzian:

∆ω(Q) )

[

1 1 1τ Qd0r0(2π)1/2

(

∫0∞ sin(Qr) exp -

[

)] ( )]

(r - d0)2 2r02

dr )

sin(Qd0) Q2r02 1 1exp τ Qd0 2

(12)

It must be noted that all jump models share the same linear variation at low Q values (long distances), where Fick’s law holds, with variations on the large Q region. In fact, since all the models are based on the same postulates (Markovian random walk with a mean jump rate), one can always obtain the width through the following equation:

∆ω(Q) )

[

1 1τ

∫0∞

sin(Qr) F(r) dr Qr

]

(13)

It is only the jump distribution that varies. The Singwi-Sjo¨lander model13,16,17 can also be defined from eq 13, and its HWHM is given by

∆ω(Q) )

Q2〈r2〉 1 6τ 1 + Q2〈r2〉/6

(14)

MD simulation is a powerful tool for extracting details of the molecular motion at different time and length scales. Direct comparison with the QENS experiments and predictions using theoretical jump diffusion models is one of the objectives of this work. 5. Results 5.1. Quasielastic Neutron Scattering Experiments. Energy spectra obtained for a fixed n-butane concentration and varying methane loadings are shown in Figure 1 for Q ) 0.39 Å-1 (the signal from the zeolite and from n-butane has been subtracted). At this Q value, the contribution from the rotation is small (6%), so that the quasielastic broadenings are essentially due to translational motion. Some of the spectra obtained for a fixed methane loading and increasing perdeuterated n-butane concentrations are shown in Figure 2. The experimental spectra for methane obtained at the different Q values were first fitted individually with a Lorentzian function assigned to the translation, convoluted with an isotropic rotation and with the instrumental resolution. The resolution function of the spectrometer was found to be better simulated with a Gaussian than with a triangle, as is traditionally done on this instrument. It appears from Figures 1 and 2 that the agreement between experimental and fitted spectra is good. The selfdiffusivities for methane derived from the dependence of the HWHM of the energy spectra on Q in the low-Q region, where Fickian diffusion is followed (see section 5.4) are shown in Figures 3 and 4. Similar values of the self-diffusivities were obtained from analyzing the HWHM in the whole Q range with a jump diffusion model, as discussed in section 5.4. The diffusivities of n-butane could also be estimated from the experiments, although the accuracy is less than for methane because the signal is weaker. At low concentrations of n-butane, a value DsC4D10 ) 4 × 10-6 cm2/s is obtained from the QENS. The diffusivity measured by QENS for n-butane was found to increase with the loading, because it is in fact the transport diffusivity that is measured for a deuterated compound.18


Gergidis et al.

Figure 1. QENS energy spectra at various methane loadings. The methane occupancy, in molecules per unit cell, is shown in the figure. The n-butane loading is held constant at 4 molecules per unit cell. All spectra are at a value of the wave vector transfer Q ) 0.39 Å-1. Fits to the spectra are presented with straight lines and experimental points with crosses. See text for details.

Figure 2. QENS energy spectra at various perdeuterated n-butane loadings. The methane loading is held constant at 4 molecules per unit cell. The value of wave vector transfer is Q ) 0.39 Å-1. Fitted spectra are presented with straight lines and experimental points with crosses.

5.2. Diffusion Coefficients from MD. The orientationally averaged value of the self-diffusivity, Ds, was extracted from the long-time, linear (“Fickian”) part of the mean squared displacement curves obtained from the MD trajectory for each species:

d 1 Ds(k) ) lim 〈(r(k)(t) - r(k)(0))2〉 6 tw∞ dt

(15)

where the index k corresponds to the different species. Direct comparisons of the self-diffusivity coefficients of methane from simulations and from QENS measurements as functions of the methane and n-butane loading are shown in Figures 3 and 4, respectively. Semilogarithmic coordinates were used to make the form of the occupancy dependence more apparent. Selfdiffusivities of methane molecules decrease with increasing loading of either component. The dependence of the methane self-diffusivity on n-butane occupancy is much stronger than that of the methane diffusivity on methane loading: As the n-butane loading rises from 0 to 7.5 molecules per unit cell, DsCH4 decreases by a factor of 20 (both for QENS and MD); whereas DsCH4 falls by a factor of 9 for QENS and a factor of 4 for MD as the methane loading increases from 2 to 10 molecules per unit cell. The MD results are consistent and in very good qualitative agreement with the behavior of the system studied at 300 K.9 The agreement between the experimental measurements and the computed values is extremely good, given the uncertainty in the potential parameters of the MD and the experimental errors inherent in the QENS. The estimated errors

Figure 3. Orientationally averaged self-diffusivities of methane in mixtures of methane and perdeuterated n-butane sorbed in silicalite at 200 K, shown as functions of the methane loading. The n-butane loading is held constant at 4 molecules per unit cell.

for the MD simulations are in the range 20 to 50%; for the QENS measurements, they range from 50 to 70%. The diffusivity values obtained from the MD simulations are somewhat higher than the diffusivity values obtained from the experiments. This is expected, because in the simulation we assume a perfectly rigid crystal lattice, without defects and impurities that would reduce mobility of the molecules. Mutual molecular interference plays a key role in the diffusion process. The



Figure 4. Orientationally averaged self-diffusivity of methane in mixtures of methane and perdeuterated n-butane sorbed in silicalite at 200 K, shown as functions of the n-butane loading. The methane loading is held constant at 4 molecules per unit cell.

TABLE 1: Self Diffusivities of n-butane as Obtained from the MD Simulationsa CH4 loading

Ds,C4D10, 10-6 cm2/s

C4D10 loading

Ds,C4D10, 10-6 cm2/s

2 4 6 10

4.1(4.0) 4.3 3.7 1.9

4 5 7.5

4.3(4.0) 2.7 0.8

a The n-butane loading is held constant at 4 molecules per unit cell while varying the methane concentration. Correspondingly, the methane loading was kept constant at 4 molecules per unit cell while varying the n-butane loading. The QENS values are in parentheses.

increase of the loading in both methane and n-butane imposes extra energy barriers on the motion of the confined molecules. The discrepancy at the high methane loadings in Figure 3 can be attributed to the larger error associated with the neutron scattering experiments for the smallest diffusivities. Self-diffusivities of n-butane, as computed by MD at various loadings, are presented in Table 1. They are seen to decrease with the loading of both n-butane and methane molecules. At low loadings, where Ds for the deuterated n-butane molecules can be extracted from the QENS measurements, excellent agreement is seen between the experimental and simulation values of DsC4D10. 5.3. Van Hove Correlation-Intermediate Scattering Functions. The self-part of the van Hove correlation function Gs(r,t), defined in eq 1 correlates the motion of a particle in space and time. Here, we will be interested in the translational motion of molecules belonging to each one of the sorbate species that is associated with their self-diffusive translational process. Typical calculations for the components Gs(ri,t) (ri ) x, y, z), which represent the probability density for a molecule to experience a displacement ri along coordinate axis i in time t relative to its position at time t ) 0, are shown in Figures 5 and 6 for methane and n-butane, respectively. The loading is 4 methane and 4 n-butane molecules per unit cell for both figures. In the case of n-butane, molecular centers of mass were used to accumulate Gs. From Figures 5 and 6, it is clear that the Gs(y,t) curve decays rapidly for both species (faster for methane) in comparison to those for the x and, especially, the z direction. Also, along all three directions, one observes fine structure (periodic maxima and minima). This structure results from the

Figure 5. Self-part of the van Hove correlation function for methane in silicalite at 200 K at a loading of 4 methane and 4 perdeuterated n-butane molecules per unit cell in x, y, and z directions.

fast local motion of the penetrants within the channel interiors (sorption sites), which are periodically distributed in space. This fast motion is intercepted by infrequent jumps between the sites. As shown by Gaub et al.,19 Gs(ri, t) along each axis can be approximated very well as a convolution of the Gaussian expression corresponding to self-diffusion in the hydrodynamic limit and a spatial autocorrelation function of the equilibrium occupancy probability along the i direction. For the silicalite crystal, which consists of a regular network of channels, the correlation in space and time actually depends on both starting and final positions and not only on their difference. The curves shown in Figures 5 and 6 represent equilibrium averages over all starting positions. In the following calculations, we will focus on the orientationally averaged Gs(r, t) for each species, which is obtained from Gs(r,t) and related to the component Gs(ri, t) as discussed by Gaub et al.19 This “isotropic” Gs(r,t) corresponds to the QENS measurements. Gs(r,t) does not display the fine structure exhibited by Gs(ri,t), as those are washed out during the orientational averaging.19 To check the consistency of our MD computations, we calculated the self-diffusivity from the intermediate scattering function Is(Q, t) in the limit of small Q and long times t, where


Gergidis et al.

Figure 7. Plot of the logarithm of the intermediate scattering function Is(Q,t) against time for methane in silicalite at 200 K at a loading of 2 methane and 4 n-butane molecules per unit cell. The shape of the function is shown for Q ) δQ, 2δQ, and 3δQ, with δQ ) 1.82 × 10-2 Å-1. Self-diffusivity calculation is shown for Q ) δQ and compared with the value of DMSD obtained from the slope of the s mean-square displacement with respect to time.

Figure 6. Self-part of the van Hove correlation function for n-butane in silicalite at 200 K at a loading of 4 methane and 4 perdeuterated n-butane molecules per unit cell in x, y, and z directions.

the relation Is(Q, t) ) exp(-DsQ2 t) holds. Figure 7 shows the linear dependence of ln Is(Q, t) on time for methane. The selfdiffusivity value Ds extracted in this way for methane is in excellent agreement with the value obtained for the slope of the mean square displacement (MSD) as a function of time in the hydrodynamic limit (see inset in Figure 7). The n-butane self-diffusivity from Is(Q, t) is also in good agreement with the MSD value but is a little more difficult to obtain, because n-butane’s diffusive process is slower; hence, the statistics in the Q w 0 limit are not adequate. 5.4. Scattering Function Line shape Analysis. For molecules performing long-range self-diffusive motion, a Lorentzian shape is expected for the scattering function (eq 6). The spectra Sinc(Q, ω), calculated from the MD simulation, are shown in Figures 8 and 9 for methane and n-butane, respectively, at 200 K. These figures show the characteristic broadening of the shape of Sinc(Q, ω) with increasing wave vector Q; the angular velocity is given in units of energy pω. Neither the convolution with instrumental resolution nor the decay of intensity due to the Debye-Waller factor was taken into account in these calculations. The broadening of the curves is more expressed for methane than for n-butane molecules, indicating that the diffusivity of methane is larger than the diffusivity of n-butane. The height

Figure 8. Incoherent dynamic structure factor (scattering function) S(Q,ω) for methane in silicalite at 200 K at a loading of 4 methane and 4 n-butane molecules per unit cell, as computed from MD.

of methane’s and n-butane’s Sinc(Q, ω) is a decreasing function of Q (Figures 8, 9). Figure 10 displays the HWHM as a function of Q2 for a given loading for both MD simulation and experiment. For small values of Q2 (long distances), the HWHM is proportional to Q2, as expected for “Fickian” diffusion. The slope of the “Fickian” line of the MD points is slightly higher than the slope of the line between the experimental points because the diffusivity predicted by the simulation is slightly higher than the self-diffusivity value from the QENS experiments. The agreement between experimental points and simulation predictions is good, as already discussed in conjunction with the diffusivity results (Figures 3 and 4). The large Q2 range (short distance region) is sensitive to details of the molecular motion. For this loading of 2 methane and 4 n-butane molecules per unit cell, the experimental HWHM goes through a maximum and then reaches a plateau value for large Q vectors. This behavior is captured well by the MD simulation, although the second hump found at about 0.6 Å-2 is not clearly observed in the experiment. In Figure 10, no experimental points are available between 0.24 and 0.56 Å-2 because Bragg peaks of the zeolite fall in this range. Along with the QENS and MD


Figure 9. Incoherent dynamic structure factor (scattering function) S(Q,ω) for n-butane in silicalite at 200 K at a loading of 4 methane and 4 n-butane molecules per unit cell, as computed from MD.

Figure 10. Half-width at half-maximum of the incoherent scattering function plotted against Q2 for methane at a loading of 2 methane and 4 n-butane molecules per unit cell in silicalite at 200 K. Both simulation predictions and points from QENS measurements are shown. Also plotted are fits to the experimental data with the Chudley-Elliott, SingwiSjo¨lander, and Jobic models.

results in Figure 10 are shown fits to the QENS points with three jump models: the Chudley-Elliott (eq 8), the SingwiSjo¨lander (eq 14), and the Jobic model introduced in section 4 (eq 12). The Singwi-Sjo¨lander model does not capture the maximum in HWHM seen in both experiments and simulation. The Chudley-Elliott model does a better job in describing the experimental points, but the small hump observed by QENS and MD at Q2 of about 0.15 Å-2 is better reproduced by the Jobic model (with parameters d0 ) 9.57 Å and r0 ) 2.5 Å). It is noteworthy that the MD predictions, without any adjustable parameters, provide as good a description of the experimental line shape as any of the fits to theoretical jump models. The picture of the HWHM as a function of Q changes upon increasing the loading. At a loading of 4 methane and 4 n-butane molecules per unit cell, the HWHM displays two humps (Figure 11). This behavior seems to be exhibited by the experimental points as well, although to a lesser degree, and even less by the fits to the theoretical models of Chudley-Elliot and of Jobic. This complicated variation found from the molecular dynamics study cannot be checked with our QENS experiment. One would have to repeat the measurements with a higher Q resolution. 5.5. Trajectory Analysis: Jump Dynamics. Previous molecular dynamics investigations9 have shown that molecular


Figure 11. As in Figure 10 but at a loading of 4 methane and 4 n-butane molecules per unit cell. Both simulation predictions and points from QENS measurements are shown. Also plotted are fits to the experimental data with the Chudley-Elliott and Jobic models.

displacements result from a combination of short oscillatory motions and of infrequent, abrupt jump-like processes that are completed within short time periods. It was apparent from our MD trajectories that molecules are temporarily trapped within energetically favorable sorption “sites”, between which jumps are executed. “Jump” lengths between neighboring sites change when the concentration and composition of the mixture change. It was evident, for example, that when the loading is increased, the mean jump length is reduced significantly. Calculations of the slope of the HWHM as a function of Q2 (section 5.4) and the fine structure in Gs(r,t) (section 5.3) also show that jump diffusion is justified for the description of the motion. The restriction of the sorbed molecules in the confining potential landscape of the zeolitic crystal gives rise to motions such as shuttling and rattling.20 These motions are primarily oscillatory in character. At low occupancy, these oscillations have different amplitudes and frequencies, set by the spacings between the pore walls and between channel intersections, which represent the natural energy barriers to the molecular mobility imposed by the crystal structure. At high occupancies, they are strongly affected by sorbate-sorbate interactions within the channels. In the present investigation, we have developed an approach for the definition of jumps along the MD trajectory of penetrant molecules, which eliminates fast local motions such as rattling and shuttling. A “jump” is considered as leading to a change of the intracrystalline environment (channel segment) where the molecule resides. Such events must be completed within “short” time scales, ensuring the abrupt character inherent in the nature of a jump. For identifying the change in environment, it is necessary to label each channel segment contained in our simulation box. In our investigations we have used a simulation box containing 16 unit cells (2×2×4 along the x, y, and z directions, respectively). The identification of channel segments was achieved by using the points representing the local minima of the potential energy hypersurface for xenon in silicalite, as identified in the transition-state theory investigations of June et al.21 and used in the hierarchical modeling work of Maginn et al.22 There are two minimum energy points in each straight channel segment, three in each sinusoidal channel segment, and one in each channel intersection. Coordinates of those minimum energy points that are located in the asymmetric unit are given in Table 2 and a picture of the unit cell of silicalite displaying a set of


Gergidis et al.

TABLE 2: Coordinates of the Minimum Energy Points Within the Asymmetric Unit of a Silicalite Unit Cell environment

x (Å)

y (Å)

z (Å)

straight (S) sinusoidal (Z) sinusoidal (Z) sinusoidal (Z) intersection (I)

0.36 4.00 5.85 1.43 0.46

0.93 4.98 4.98 4.98 4.98

6.11 12.13 10.88 11.56 6.06

TABLE 3: Numbers of Transitions and Transition Rate Constants for Methane 2CH4-4C4D10

10CH4-4C4D10

transition iwj

jumps

kiwj, ns-1 a

kiwjb

jumps

kiwj, ns-1 a

kiwjb

SwS SwZ ZwZ ZwS

6386 4308 2270 4316

41.64 28.09 16.86 32.05

33.45 27.80 13.50 25.79

35690 16181 14044 16182

52.59 23.84 18.45 21.254

29.70 15.50 14.5 13.81

a From enumeration of jumps. b From distribution of residence times between jumps.

TABLE 4: Numbers of Transitions and Transition Rate Constants for n-Butane 2CH4-4C4D10

transition iwj

jumps

SwS SwZ ZwZ ZwS

2140 1196 397 1175

a

kiwj,

ns-1 a

10CH4-4C4D10 kiwj

jumps

kiwj, ns-1 a

kiwjb

6.9 3.2

3115 2612 846 2608

10.76 9.02 2.95 9.10

7.9 7.0

b

8.38 4.68 1.26 3.66

2.2

7.2

b

From enumeration of jumps. From distribution of residence times between jumps.

Figure 12. Points representing energy minima in the three silicalite environments. Sinusoidal channel points are represented in black, with straight points in gray and intersection points with very light gray. The silicalite crystal unit cell is represented with light gray for reference.

minimum energy points is shown in Figure 12. The 4×2 + 4×3 + 4 ) 24 points shown in Figure 12 are all obtained from the points listed in Table 2 by application of the symmetry operations of the zeolite. To assign a molecule to a particular segment, we calculate the distances between the position of the molecule and the positions of all minimum energy points in the lattice, and we select that minimum point for which this distance becomes minimal. The molecule is assigned to the channel segment or intersection containing this point. With this designation, the actual trajectory of a molecule can be reduced to a coarse-grained trajectory consisting only of jumps between channel segments. Our algorithm follows the actual dynamical trajectory for each molecule in time and tracks the environment (particular straight, sinusoidal channel segment, or intersection) in which the molecule finds itself at each time. As long as the molecule remains in the same channel segment or moves into an intersection region, we do not have any jump event. As soon as the molecule moves into a new channel segment, we mark a jump. The position assigned to the molecule between any two jumps is its average position over the interval that intervened between the jumps. The jump length is taken as the difference between the average position in the “old” environment and the average position in the “new” environment. With this approach, there are four types of jumps: jumps between adjacent straight channel interiors (S w S), jumps between adjacent sinusoidal channel segments (Z w Z), and jumps between straight and sinusoidal channel segments (S w Z, Z w S), which “turn a corner” in the zeolite. It was decided to absorb intersections into the channel states and not consider motions originating or terminating at intersections as separate jumps. This was done to keep the jump model simple, in view of the fact that intersections correspond to high-energy reaction intermediates21 and are therefore visited for only short periods of time. With these definitions, we are able to fit the actual trajectory of each sorbate molecule with a succession of steps. Such

trajectory “fits” are shown in Figures 13 and 14 for methane and n-butane, respectively, at two different occupancies. The jumplike character of the motion is clear from these pictures. It is also clear that our stepwise fits to the trajectories follow the actual trajectories very closely, providing justification for the coarse-graining algorithm adopted here. The pictures in Figures 13 and 14 are quite similar to the pictures from our previous investigation9 at 300 K although now, at the lower temperature of 200 K, methane’s motion is clearly more jumplike. A first observation from Figures 13 and 14 is that, at low occupancy, the jumplike character is more pronounced for the n-butane than for the methane molecules. A physical explanation for this phenomenon is that, at low occupancy, the methane molecules can move through the different silicalite landscapes more easily due to their lower interaction energies with the zeolite. Methane molecules are not strongly localized in the channel interiors, as happens for the n-butane molecules, due to a different balance between energetic and entropic effects.9 At high occupancies, however, other methane molecules and, especially, the more slowly moving n-butane molecules impose additional obstacles to methane mobility. A second observation from Figures 13 and 14 is a drastic decrease in the jump length as occupancy rises for both species. The molecules are trapped in the channel segments for long periods of time and move much less whenever a jump occurs. Of course, this affects the diffusivity values, as expected (see section 5.2). For each transition from initial state i to destination state j, the calculations should obey microscopic reversibility:

pi kiwj ) pj kjwi

(16)

where pi is the equilibrium probability of occupancy of state i and kiwj is the rate constant for a i w j jump. In all the calculations, the above condition was fulfilled. This is characteristically seen in Tables 3 and 4, which give the total numbers of transitions of each type for each species, recorded during the 9-ns-long simulation runs. The distributions of the jump lengths for methane at two extreme loadings are shown in Figure 15. It should be stressed that these jumps are performed within a time step of 3 ps. This explains why the peak of the distribution in the low occupancy case is centered at about 6 Å, whereas the experimental value is 9.57 Å, because the time scale of the QENS measurement is


Figure 13. Distances traversed by a methane molecule (thin and dashed lines) and trajectory fits with successions of jumps (thick, straight lines) shown as functions of time for two different sorbate loadings. The upper curves are for a loading of 4 n-butane and 2 methane molecules per unit cell; the lower ones are for 4 n-butane and 10 methane molecules per unit cell. Origins for the measurement of distance are arbitrary.

Figure 14. Distances traversed by an n-butane molecule (thin and dashed lines) and trajectory fits with successions of jumps (thick, straight lines) shown as functions of time for two different sorbate loadings. The upper curves are for a loading of 4 n-butane and 2 methane molecules per unit cell; the lower ones are for 4 n-butane and 10 methane molecules per unit cell. Origins for the measurement of distance are arbitrary.

much longer (at least 40 ps). Increasing the methane loading to 10 methane molecules per unit cell shifts this peak to approximately 2.5 Å. The tail of the distribution is also shifted to shorter distances by increasing the methane loading. We are able to resolve the total jump length distribution into contributions from the individual types of transitions between the channel segments. In Figure 16, we show the jump length distributions for the particular transitions of methane. The first peak (around 1.5 Å) originates in our having absorbed intersection states into the (mainly straight) channel segment states. Crossing motions from a straight channel interior to an intersection or the reverse can occur over a distance of 1.5 Å (this kind of transition displays the larger population of jumps at high occupancy, for both methane and n-butane molecules). The Z w Z transition distributions extend to somewhat larger lengths because the geometric distance, in three dimensions, between two neighboring sinusoidal channel segments is larger than that between two straight or between two adjacent straight and sinusoidal channel


Figure 15. Distributions of jump lengths for methane molecules at different loadings. Straight and broken lines correspond to loadings of 2 methane and 4 n-butane and 10 methane and 4 n-butane molecules per unit cell, respectively.

segments. The characteristics of S w Z and Z w S transitions are identical, as expected, and display shorter tails. Increasing the loading results in a decrease of the jump length for all transitions, but most dramatically for the transitions between straight channels. At the higher loading, methane molecules are pushed to the edges of straight channels.9 Jump length distributions for the other three types of transitions are almost indistinguishable. Figure 17 shows the jump length distributions for the centers of mass of n-butane molecules for two different methane loadings. The picture is very different from that seen in the case of methane. There are at least two peaks in the n-butane distributions. Figure 18 clearly shows the origin of the peaks. The first peak comes, again, from short jumps between intersection regions (absorbed in the straight channels) and adjacent straight channels. The peak around 5 to 6 Å is from n-butane molecules undergoing S w Z or Z w S transitions. There is another “turning the corner” peak due to S w Z and Z w S jumps around 8 Å. The S w S distributions display one extended peak at approximately 9 to 10 Å, corresponding to jumps between the interiors of neighboring channel segments. The Z w Z distribution again extends to longer distances in comparison to the distributions for other transitions for the reason discussed above. It displays a sharp peak around 10.5 Å. When the methane loading is increased to 10 molecules per unit cell (the n-butane being held constant at 4 molecules per unit cell, Figure 18b), practically all peaks beyond 7 Å vanish, indicating that the likelihood of a direct jump between the interiors of adjacent straight or sinusoidal channels is very small. At this extreme occupancy, there is a dramatic increase of the peak at 1.5 Å, which is indicative of fast crossings in the neighborhood of the intersection regions. The picture extracted from the jump length distributions is consistent with the picture derived from direct monitoring of the real- and coarse-grained trajectories. Tables 3 and 4 present the rate constants of individual jumps, calculated for each species by the following relation:

kiwj ) number of kiwj jumps/(number of molecules × probability of state i × simulation time) (17) Interestingly, increasing the loading of methane from 2 to 10 molecules per unit cell leads to an increase of the rate


Gergidis et al.

Figure 17. Distributions of jump lengths for n-butane molecules at different loadings. Straight and broken lines correspond to loadings of 2 methane and 4 n-butane and 10 methane and 4 n-butane molecules per unit cell, respectively.

Figure 16. Distributions of jump lengths for methane molecules at different loadings, differentiated according to the jump type: (a) 2 methane and 4 n-butane and (b) 10 methane and 4 n-butane molecules per unit cell. S w S and Z w Z transitions are shown as thick and thin straight lines, respectively. S w Z and Z w S transitions are shown as thick and thin broken lines, respectively.

constants for S w S and Z w Z transitions, primarily because of an abundance of fast, short-jump processes involving the intersection regions (peak around 1.5 Å in Figures 15 and 16). For n-butane, increasing the occupancy leads to an increase in the rate constants of all transitions, again due to fast local crossings at the intersections. It is interesting that the abundance of “straight-through” S w S transitions is suppressed relative to that of “turning the corner” S w Z transitions at high occupancy, as is also clear from the jump length distributions of Figure 18. An alternative way to estimate the transition rate constants kiwj from the MD trajectories is to accumulate the distribution of waiting times, tR,iwj that a molecule spends in state i before a jump to state j occurs. For a Poisson process consisting of independent jump events, the probability density of each tR,iwj should be exponential,23 of the form e-kiwj tR,iwj. Plotting the logarithm of the probability density versus tR,iwj would thus yield the rate constant kiwj as the negative of the slope. Actual distributions of the waiting times tR,iwj for each transition show a rapidly dropping, convex initial portion that is a signature of fast, correlated recrossing events, wherein a molecule leaving environment i to enter environment j does not thermalize in j but quickly returns back to i. The jumps over ca. 1.5-Å distances

discussed in connection with Figures 15-18 belong largely to this category. At long tR,iwj, the distributions of waiting times display the exponential decay characteristic of a Poisson process. Estimates of the rate constants kiwj extracted from the longtime slopes of waiting time distributions are given in Tables 3 and 4 along with the corresponding values from enumerating the transitions, discussed above (eq 17). For very infrequent transitions, the waiting time distributions accumulated over the duration of the simulation are subject to large statistical error; no estimates are given in Tables 3 and 4 for these cases. As seen in Tables 3 and 4, kiwj values extracted from the waiting time distributions are lower than the corresponding values from enumerating transitions, as the latter are not corrected for dynamical recrossing events. From the transition rate constants and equilibrium state probabilities, one can calculate a mean waiting time for a molecule between jumps as

〈tR〉 )

1

∑i ∑j pi kiwj

(18)

In our case, i, j ∈ {S, Z}. In the above equation, if kiwj is estimated by enumerating jumps (eq 17), 〈tR〉 is simply the total simulation time divided by the total number of jumps of any kind per molecule for the molecular species considered. Mean waiting times obtained through eq 18 from such enumerations of jumps are given in Tables 5 and 6. In the same table are given the mean jump lengths 〈d〉 for each species (over all kinds of jumps), as calculated from the overall jump length distributions of Figures 15 and 17. One can attempt a crude estimate of the self-diffusivity Ds through the equation

Ds )

〈d〉2 6〈tR〉

(19)

This equation would be strictly valid for a random walk of constant step size 〈d〉 and time intervals 〈tR〉 between successive jumps. Estimates of Ds based on eq 19 are also given in Tables 5 and 6, along with the true self-diffusivities extracted from the long-time slope of the mean square displacement curves.


J. Phys. Chem. B, Vol. 104, No. 23, 2000 5551 Summary and Conclusions

Figure 18. Distributions of jump lengths for n-butane molecules at different loadings, differentiated according to jump type: (a) 2 methane and 4 n-butane and (b) 10 methane and 4 n-butane molecules per unit cell. S w S and Z w Z transitions are shown as thick and thin straight lines, respectively. S w Z and Z w S transitions are shown as thick and thin broken lines, respectively.

TABLE 5: Mean Jump Lengths, Mean Residence Times between Jumps, and Self-Diffusivity Values for Methane loading

〈d〉 Å

〈tR〉, ps

, Dcrude s 10-5cm2/s

2CH4-4C4D10 10CH4-4C4D10

6.08 3.54

16.6 17.5

3.7 1.2

Dmsd s , 10-5cm2/s 2.0 0.6

TABLE 6: Mean Jump Lengths, Mean Residence Times between Jumps, and Self-Diffusivity Values for n-Butane loading

〈d〉 Å

〈tR〉, ps

, Dcrude s 10-6 cm2/s

Dmsd s , 10-6 cm2/s

2CH4-4C4D10 10CH4-4C4D10

5.91 2.91

108.08 57.6

5.4 2.5

4.0 1.9

One sees that eq 19 overestimates the self-diffusivity by roughly a factor of 2. The reasons for this are that (a) there is a whole distribution of jump lengths and jump times, so an estimate based on single, average jump length and jump time values, eq 19, is too crude; and (b) there is an abundance of fast, correlated jumps over short distances, as a result of which estimates of kiwj from single enumeration of jumps are too high; hence, (eq 18) estimates of 〈tR〉 are too short, leading to an overestimation of Ds through eq 19.

Results from quasielastic neutron scattering experiments and molecular dynamics simulations of binary mixtures of n-butane and methane in silicalite at various loadings and compositions have been reported. Calculations of the self-part of the van Hove correlation function for both species from the MD along the three principal directions of the crystal clearly indicate different molecular mobilities in the different crystal directions and exhibit a “fine structure” due to spatial correlations resulting from the specific siting preferences of the molecules in the pores. The fine structure disappears when the orientationally averaged Gs(r,t) is considered.19 Methane molecules are more mobile than their n-butane counterparts at low loadings. At higher loadings, especially for the high loadings of mixture in n-butane, the selfdiffusivities of the two species approach each other, as the more sluggish n-butane molecules inevitably block the translational progress of methanes as well. Calculation of the intermediate scattering function Is(Q,t) revealed an exponential dependence on time t at low Q, yielding the same diffusivity values as those obtained from the mean square displacement. Going one step further, the incoherent dynamic structure factors Sinc(Q,ω), calculated by Fourier transformation of Is(Q,t), show the expected Lorentzian shape, whereas the broadening of the curves with time (proportional to self-diffusivity at low Q) is more pronounced for methane than for n-butane molecules. A direct comparison of the broadening of the scattering function, as this can be characterized through the half-width at half-maximum (HWHM) from QENS and from MD at an occupancy of 2 methane and 4 n-butane molecules per unit cell, yielded good agreement for all Q examined. Thus, MD simulations are seen to track not only the long-range diffusive motion but also the local motions of the sorbate molecules, in very good qualitative and quantitative agreement with the QENS experiments. For the low Q values (in the linear, “Fick’s law” regime), the agreement is excellent. In the region of shorter distances, (large Q values) an undulating behavior leading to a plateau value is seen in both the MD simulations and in the experiments. This behavior indicates the existence of a jump diffusion process. Theoretical models of jump diffusion, such as the Chudley-Elliott and the more sophisticated Jobic model, were shown capable of providing a good representation of the experimental broadening. We used a coarse-graining approach to reduce the actual trajectories of the sorbate molecules into sequences of jumps between different intracrystalline environments [straight (S) and sinusoidal (Z) channel segments]. With this approach, we were able to calculate jump distances and residence times (jump rates) for S w S, S w Z, Z w Z, and Z w S jumps. Coarse-grained trajectories were seen to track the actual dynamical trajectories very well. Our jump analysis has shown a dramatic decrease in jump lengths, as well as interesting shifts in the jump rate constants with increasing occupancy. The mean waiting time between jumps actually becomes somewhat shorter with increasing occupancy, reflecting mainly an increase in fast, shortlengthscale recrossing events in the neighborhood of intersection regions. On the other hand, long S w S and Z w Z jumps across an intersection become less abundant in relation to “turning the corner” S w Z jumps for n-butane at high loadings. The calculation of the orientationally averaged self-diffusion coefficient from simulations ranging in duration from 12 ns (for the lower loadings) to 20 ns (for the higher loadings) has shown that the self-diffusivity decreases monotonically with concentration of either component. Direct comparison of the MD

5552 J. Phys. Chem. B, Vol. 104, No. 23, 2000 predictions with the QENS measurements for the methane selfdiffusivity showed excellent agreement, considering the errors inherent in each technique. Also, a comparison of the selfdiffusivity of n-butane with the corresponding estimate that can be extracted from QENS at low occupancy gave equally satisfactory agreement. It is clear that simulation techniques, such as molecular dynamics, can not only predict the trends and the functional forms describing the dependence of the selfdiffusivity on loading, but are also in very good quantitative agreement with QENS measurements at all length scales; at the same time, they provide a concrete interpretation of the observed dynamics. The combination of dynamic simulations and quasielastic neutron scattering experiments thus emerges as a powerful tool for elucidating the mechanism and rates of dynamics in confinement. Acknowledgment. We thank A. F. Terzis, M. Doxastakis, V. Mavrantzas and Prof. M. Beé for helpful discusions. We would like also to thank the Institute of Chemical Engineering and High Temperature Chemical Processes in Patras for generously providing part of the computer resources used in this work. L.N.G. gratefully acknowledges the General Secretariat of Research and Technology of Greece for financial support in the form of a PENED′95 program, No. 218, and the Institut Laue-Langevin, Grenoble, for the hospitality shown to him during his visit there. References and Notes (1) Barrer, R. M. Zeolites and Clay Minerals; Academic Press: London, 1978. (2) Ka¨rger, J.; Ruthven, D. M. Diffusion in Zeolites and Other Microporous Solids; Wiley-Interscience: New York, 1992.

Gergidis et al. (3) Bell, A. T.; Maginn, E. J.; Theodorou, D. N. In Handbook of Heterogeneous Catalysis; Ertl, G., Kno¨zinger, H., Weitkamp, J., Eds.; VCH: Weinheim, 1997; Vol. 3, p 1165. (4) Jobic, H.; Ernst, H.; Heink, W.; Ka¨rger, J.; Tuel, A.; Beé M. Micropor. Mesopor. Mater. 1998, 26, 67. (5) Jobic, H.; Beé, M.; Kearley, G. J. Phys. Chem. 1994, 98, 4660. (6) Snurr, R. Q.; Ka¨rger, J. J. Phys. Chem. B. 1997, 101, 6469. (7) Jost, S.; Fritzsche, S.; Haberlandt, R. Chem. Phys. Lett. 1997, 219, 385. (8) Jost, S.; Bar, N. K.; Fritzsche, S.; Haberlandt, R.; Ka¨rger, J. J. Phys. Chem. B. 1998, 102, 6375. (9) Gergidis, L. N.; Theodorou, D. N. J. Phys. Chem. B 1999, 103, 3380. (10) Stepanov, A. G.; Shubin, A. A.; Luzgin, M. V.; Jobic, H.; Tuel, A. J. Phys. Chem. B 1998, 102, 10860. (11) Millot, B.; Me´thivier, A.; Jobic, H.; Moueddeb, H.; Beé, M. J. Phys. Chem. B 1999, 103, 1096. (12) Frenkel, D.; Smit, B. Understanding Molecular Simulation: from Algorithms to Applications; Academic Press: San Diego, 1996. (13) Beé, M. Quasielastic Neutron Scattering; Adam Hilger: Bristol, 1988. (14) Jobic, H.; Beé, M.; Kearley, G. J. Zeolites 1989, 9, 312. (15) Chudley, C. T.; Elliott, R. J. Proc. Phys. Soc., London 1961, 77, 353. (16) Jobic, H. In Recent AdVances in Gas Separation by Microporous Membranes; Kanellopoulos N., Ed.; Elsevier: New York, in press. (17) Singwi, K. S.; Sjo¨lander, A. Phys. ReV. 1960, 119, 863. (18) Jobic, H.; Ka¨rger, J.; Beé, M Phys. ReV. Lett. 1999, 82, 4260. (19) Gaub, M.; Fritzshe, S.; Haberlandt, R.; Theodorou, D. N. J. Phys. Chem. B 1999, 103, 4721. (20) June, L. R.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1992, 96, 1051. (21) June, L. R.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1991, 95, 8866. (22) Maginn, E. J.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1996, 100, 7155. (23) Crame´r, H. The Elements of Probability Theory and Some of Its Applications; Krieger: New York, 1955.

Dynamics of n-Butaneâ Methane Mixtures in Silicalite, Using

Recommend Documents

Dynamics of n-Butaneâ Methane Mixtures in Silicalite, Using