12748
J. Phys. Chem. B 2004, 108, 12748-12756
Transport Diffusivity of N2 and CO2 in Silicalite: Coherent Quasielastic Neutron Scattering Measurements and Molecular Dynamics Simulations George K. Papadopoulos,‡ Herve´ Jobic,† and Doros N. Theodorou*,‡ Institut de Recherches sur la Catalyse, CNRS, 2 AVenue Albert Einstein, 69626 Villeurbanne, France, and School of Chemical Engineering, National Technical UniVersity of Athens, GR 15780 Zografou Campus, Athens, Greece ReceiVed: February 18, 2004; In Final Form: June 14, 2004
A direct comparison of experimental measurements (coherent quasielastic neutron scattering) and simulation (molecular dynamics) is carried out for the first time to investigate the concentration dependence of transport diffusivity of nitrogen and carbon dioxide in silicalite-1. Theoretical elaboration of the simulation results on the basis of the quasichemical mean field approximation and a model of surface diffusion due to D. A. Reed and G. Ehrlich (Surf. Sci. 1981, 102, 588) leads to the conclusion that the corrected (Darken) diffusivity D0 can exhibit significant dependence on the loading θ of the sorbed phase, varying with the strength of interactions between the sorbate molecules. As a consequence of sorbate-sorbate interactions being considerably more attractive in CO2/silicalite-1 than in N2/silicalite-1, D0(θ) is a decreasing function in CO2/silicalite-1 but displays a weak maximum in N2/silicalite-1.
1. Introduction The literature contains numerous reports on experimental and theoretical studies of confined fluids (sorbates) in either amorphous or crystalline porous materials (refs 1-4 and references therein). The interpretation of equilibrium and transport data in such materials is often a complicated task because of the underlying complexity of their pore network. This is especially true of zeolitic materials, where the regular geometry, orientation, and connectivity of nanopores and the possible existence of defects within crystals (intracrystalline morphology), along with intergrowths and boundaries between crystals and meso- and macropores (intercrystalline morphology), give rise to complicated variations in the potential field experienced by sorbed molecules.1 Consequently, classical theoretical approaches to the transport under high confinement are difficult to apply because of the discontinuities introduced in the fluid by the strong, spatially varying potential field of the porous medium. Computer simulation studies of microporous solids at a molecular level, given a model for all the interactions present, have contributed significantly to a better understanding of many aspects of the equilibrium and kinetics of sorbed phases. Measurements of isothermal sorption of a variety of sorbates over a range from low to high occupancies have fully validated simulation predictions from Monte Carlo techniques.3-6 In contrast, the interpretation of the transport mechanism in zeolites and other nanoporous materials based on experimental measurements has proved to be a much more arduous task than the prediction of sorption. This difficulty, which is more pronounced in zeolitic media possessing structure over a variety of length scales (e.g., membranes and beds), appears primarily as a * To whom correspondence should be addressed. E-mail: doros@ central.ntua.gr. Fax: +30 21 07723112. Phone: +30 21 07723157. ‡ National Technical University of Athens. † Institut de Recherches sur la Catalyse.
consequence of the intergrowths and defects in real zeolite crystals, as well as of the interplay between zeolitic, Knudsen, and molecular diffusion in the nano-, meso-, and macropores of these materials. A related complication is the sometimes poor agreement between experimental diffusivity values obtained by different techniques that probe sorbate motion over different length scales (e.g., micro- vs macroscopic diffusion measurements).1 Moreover, the time dependence of phenomena involved in the measurement of diffusivity by a particular method (e.g., transient or steady state) can also cause noticeable variations in the measured values.1,6 Previous simulation studies have been concerned mainly with sorption and diffusion of several gases either in microporous amorphous substrates by modeling a part of their pore structure (e.g., single pores of various geometries or canonical networks thereof)9,12 or in zeolites by atomistic representation of their crystal unit cell.5,6,10 In the above studies, molecular interactions within the sorbed phase have been represented by several effective (usually pair) potentials (e.g., Lennard-Jones or Buckingham), with the long-range electrostatic interactions due to the presence of multipoles being modeled by Coulomb potentials and calculated via the Ewald summation method.7 Interactions of sorbate molecules with the solid phase have been similarly modeled by pairwise additive analytical expressions for the potential energy of interaction between individual atoms (sites) on the solid and individual sites on the sorbate molecule;11 an alternative strategy involves the suitable calibration of the model surface with respect to experimental data such as the energy of adsorption of a given fluid, collected under the conditions of the simulation experiment.12 Atomistic computer simulation of transport belongs to the category of “microscopic” techniques.1 Consequently, its results can be compared directly to those from pulsed field gradient nuclear magnetic resonance (PFG-NMR)1,13 and quasielastic neutron scattering (QENS) measurements.8 Comparisons between atomistic molecular dynamics (MD) simulation and PFG-
10.1021/jp049265g CCC: $27.50 © 2004 American Chemical Society Published on Web 08/04/2004
Transport Diffusivity in Silicalite
J. Phys. Chem. B, Vol. 108, No. 34, 2004 12749
NMR or QENS experiments are usually conducted at the level of the self-diffusivity, Ds, leading to excellent results.4 In practical separations and catalytic applications, however, it is the transport diffusivity, Dt, that is of greater importance. A recent simulation study6 has given excellent agreement between predicted and experimental sorption properties of carbon dioxide and nitrogen through silicalite-1. Attempts to estimate transport diffusivity in the same systems via the Darken relation, however, failed in reproducing absolute experimental permeability values in silicalite-1 membranes under nonequilibrium conditions. The latter result has stimulated an investigation of sorbate transport in relation to the macroscopic structure of zeolitic membranes, as described in ref 14. At the molecular level, however, the dependence of the transport diffusivity Dt and the so-called corrected diffusivity D0 on intracrystalline occupancy needs to be resolved. The definition of Dt involves the response to an intracrystalline concentration (chemical potential) gradient; therefore, its direct determination calls for nonequilibrium experiments. Thanks to linear response theory, however, Dt and D0 can be determined directly under equilibrium conditions, where chemical potential gradients are absent. Experimentally, this objective can be accomplished by coherent QENS, which probes the collective motion of sorbate molecules at equilibrium.39 With simulation, the same objective can be accomplished using equilibrium molecular dynamics (EMD). The objectives of the present work are (a) to obtain the transport (Dt) diffusivities of N2 and CO2 sorbed in silicalite-1 as functions of intracrystalline occupancy by QENS measurements, (b) to predict the transport and the corrected (D0) diffusivities from EMD computer experiments, and (c) by taking advantage of the opportunity of a direct comparison between computed and measured Dt coefficients to understand better the molecular factors that govern the occupancy dependence of D0, a “pure” kinetic parameter that is largely free of the influence of the isotherm.17 In previous studies,5,14,15 D0 has often been assumed to be independent of sorbate concentration and approximately equal to the self-diffusivity, Ds, in the limit of zero occupancy; this assumption will be reassessed in light of the results obtained here. To our knowledge, this is the first time that a direct comparison of Dt values obtained directly from simulation and from a microscopic experimental method is being presented. 2. Transport Theory Activated Diffusion. The diffusion mechanism in the confined spaces of zeolite crystals depends basically on the concentration and the ratio of pore width over sorbate molecular size; when these sizes become comparable, the diffusing molecules remain subject to the strong potential field of the crystal at all times. Since this field is usually spatially inhomogeneous, transport occurs chiefly via random jumps of thermally activated molecules between “sites”, or regions around potential minima, separated by energy barriers Ea; this behavior is frequently met inside the lattice of nanoporous zeolitic crystals and referred to as an intracrystalline diffusion process being subject to a temperature dependence of the Arrhenius type, i.e.,
( )
D ) D∞ exp -
Ea RT
(1)
where R is the gas constant and T is the absolute temperature.
Macroscopic. A general expression for the isothermal steadystate single-component sorbate flow in the absence of external force fields can be written in Fickian form as follows:
J ) -Dt(F) ∇F
(2)
In eq 2, the vector J is the macroscopic flux, expressed as the mean number of molecules flowing per unit cross-sectional area of the nanoporous material due to a concentration (number density F) gradient, and Dt is the transport (Fickian) diffusion coefficient, which in general can be a function of concentration. Under the irreversible thermodynamics formalism, wherein the gradient of chemical potential over thermal energy acts as the driving force for diffusion of the sorbate species 1, eq 2 can be equivalently written in terms of the phenomenological coefficients L according to Onsager notation (assuming that the zeolitic phase 2 is stationary so that L12 ) 0 and hence L11 ) L) as follows:4
J ) -β L(F) ∇µ
(3)
where β ) 1/(kBT), kB ) R/NL is the Boltzmann constant, and NL is Avogadro’s number. The chemical potential gradient can be reexpressed as
∇µ )
1 ∇f βf
(4)
with f being the fugacity of the sorbed fluid at the considered position. After manipulation and combination with eq 4, eq 2 reads
J ) -Dt(F) βF
∂ ln F ∇µ ∂ ln f
(5)
Comparison of eqs 3 and 5 leads to the following relation
Dt(F) )
L(F) ∂ ln f ∂ ln f ) D0(F) F ∂ ln F ∂ ln F
(6a, 6b)
Equation 6b, with D0(F) replaced by Ds, is known as the Darken approximation. Moreover, by eqs 6a and 6b the transport coefficient L and the corrected diffusivity D0 are related to the more commonly used transport diffusivity Dt. However, the experimental determination of Dt in terms of eqs 6a and 6b has not been practiced up to now, since neither L nor D0 is a directly measurable quantity in the majority of experimental techniques. In the literature, relations of structure analogous to that of eq 6b whose derivation is based more or less on kinetic arguments have been proposed;16 nevertheless, these models also suffer from inadequacy in providing a direct prediction of Dt. Statistical Mechanics. The elaboration of molecular dynamics results toward the estimation of transport coefficients relies on linear response theory (LRT).18 This provides the bridge between equilibrium time correlation functions and nonequilibrium response to weak perturbations. Onsager, through his regression hypothesis, namely, that the relaxation of weak nonequilibrium disturbances follows the same laws as the regression of the spontaneous fluctuations at equilibrium, made the first attempt toward such a bridge before LRT formalized it.33 LRT proves that the time integrals of autocorrelation functions are related to transport coefficients, via relations known as Green-Kubo, of the following type:
Ds )
1 d0
∫0∞ dt 〈vi(t)‚vi(0)〉
(7a)
12750 J. Phys. Chem. B, Vol. 108, No. 34, 2004
Papadopoulos et al.
Equation 7a relates the self-diffusivity Ds (averaged over a d0 dimensional space) and the time integral of the autocorrelation function of the translational (center of mass) velocity v of the sorbate molecule i. In MD simulations, the above equation is usually used in the averaged form of eq 7b, which allows accumulating statistics over all the N molecules of the system:
Ds )
1
∫0
∞
d0
dt
〈
N
1
∑vi(t)‚vi(0) N i)1
〉
(7b)
An alternative, mathematically equivalent relation to eq 7b is the Einstein equation, where the average is taken over time for the mean squared displacement of the center of mass position vectors r of all the molecules in the system, i.e.,
Ds ) lim tf∞
〈∑
d 1 1
〉
N
dt 2d0 N i)1
[ri(0) - ri(t)]2
(8)
Thus, the self-diffusivity Ds expresses the displacement of a tagged molecule among the remaining N - 1 untagged molecules. It is the property measured using microscopic experimental techniques such as PFG-NMR, NMR lifetime, and incoherent QENS.1 According to LRT, the nonequilibrium macroscopic flux J described by eq 3, which develops in response to a chemical potential gradient under steady-state conditions, can be written in terms of the time integral of the autocorrelation function of the instantaneous flux j(t) under equilibrium conditions (absence of chemical potential gradient) as follows:18
J ) 〈j〉neq ) -
V d0
∫0∞ dt 〈j(t)‚j(0)〉β∇µ
(9)
where V is the volume occupied by the system and the instantaneous flux j(t) is defined as
j(t) )
1
N
∑vi(t) V i)1
where vi ) r3 i is the center-of-mass velocity of molecule i. Comparing eqs 3, 6, and 9, the microscopic expression of D0 is identified with
D0 )
N
1
N
∞ dt ∑∑〈vi(t)‚vj(0)〉 ∫ 0 dN
) Ds +
N
1
N d0
N
∞ dt ∑∑〈vi(t)‚vj(0)〉 ∫ 0 dN
(10b)
i)1 j)1 j*i
0
)
(10a)
i)1 j)1
0
∫0∞ dt 〈u(t)‚u(0)〉
(10c)
where u(t) is the center-of-mass velocity of the swarm of N molecules, which can be identified with the macroscopic streaming velocity of the system, i.e.,
u(t) )
1
N
∑ vi(t)
N i)1
Clearly, u(t) ) R4 (t) with R being the center-of-mass of the swarm of N molecules.
By analogy to the set of eqs 8 and 7, eq 10c may equivalently be written in terms of the mean squared displacement of the center of mass position of the swarm, i.e.,
D0 ) lim tf∞
d N 〈[R(0) - R(t)]2〉 dt 2d0
(11)
An alternative way of deriving the theoretical expression of D0 in the Fourier space in combination with Onsager’s hypothesis has been recently employed by Briels and co-workers.27 From eq 10b it is clear that as sorbate loading tends to infinite dilution, cross-correlation functions between velocities of different molecules approach zero and the autocorrelation terms (self-diffusivity) predominate. Also, the sorption isotherm approaches the linear (Henry’s law) regime. In this limit,
lim D0(F) ) lim Dt(F) ) lim Ds(F) Ff0
Ff0
Ff0
The statistical mechanics formulation of transport through eqs 6, 7, and 10 clearly indicates the physical character of each diffusion coefficient used in the study of sorbate flow. Thus, Ds is a pure kinetic property of individual molecules, D0 is a pure kinetic property19 that has a collective nature, since in addition to Ds it comprises the sum of the off-diagonal crosscorrelation terms in the array of the velocity dot products (eq 10b). On the other hand, Dt, which is also a collective property, embodies the isotherm shape (eq 6b), namely, direct dependence on the sorption thermodynamics. All three properties are in general concentration-dependent, exhibiting different concentration dependences. An alternative to EMD for obtaining transport diffusivity is the utilization of nonequilibrium molecular dynamics techniques (NEMD), either beyond the linear response in the form of transient time correlation functions,20 which can be considered as the nonlinear generalization of the Green-Kubo formulas, or within the linear response regime by calculating the sorbate flux caused by a (macroscopic) chemical potential gradient under transient5,21 or steady-state flow conditions.5,22,23 The work of Maginn et al.5 was the first application of NEMD in order to predict transport diffusivity in zeolites. They used two methods for the study of methane transport inside a silicalite-1 crystal; one was based on a gradient relaxation molecular dynamics technique that mimics macroscopic experiments, while the other applied an external force generated by a fictitious “color” field, which in that method plays the role of the driving force for diffusion. A disadvantage of the first method is the difficulty in establishing an initial macroscopic concentration profile and in ensuring that the system is maintained in the linear response regime. Both these difficulties are overcome by the color field method. In the latter, however, very weak fields may be required to keep the system within the linear response regime, resulting in high statistical noise; reliable values of the macroscopic flux can only be accumulated through very long runs. An extensive comparison of EMD, NEMD, and hybrid methods for direct determination of transport diffusivity such as the dual control Volume grand canonical molecular dynamics under constant chemical potential gradient22 is found in the work of Maginn and co-workers.23 Sholl and co-workers24 have used standard equilibrium MD techniques for the prediction of nonequilibrium transport coefficients of various sorbates in silicalite-1 and carbon nanotubes;26 also, Krishna and co-workers24 have performed studies toward the prediction of transport diffusivities for pure sorbates and mixtures in zeolites on the basis of the MaxwellStefan theory. Sanborn and Snurr25 have studied binary mixtures
Transport Diffusivity in Silicalite
J. Phys. Chem. B, Vol. 108, No. 34, 2004 12751
in faujasite using MD in order to obtain structural and transport properties of the system. Briels and co-workers27 have performed a study of transport diffusivity of Ar sorbed in AlPO4-5 based on equilibrium MD runs by employing the Green-Kubo approach in Fourier space; a wave-vector-dependent diffusivity was obtained, allowing them to assess the correct linear-response result. 3. Experimental Details 3.1. QENS. The QENS experiments were performed at the Institut Laue-Langevin, Grenoble, using the spectrometer IN6. The energy resolution depends on the incident neutron energy; in this work, the incident energy was taken as 3.12 meV (5.12 Å). The elastic energy resolution was measured with a vanadium plate. The line shape was of a Gaussian form with a half-width at half-maximum (hwhm) increasing with increasing scattering angle (10° < 2θ < 115°). The hwhm of the resolution varied typically from 40 µeV for small scattering angles to about 50 µeV for the largest scattering angle. The range of momentum transfer was 0.22 Å-1 < pQ < 2.07 Å-1, where Q ) k0 - k1 is the scattering vector, 2pπ is the Planck constant, and k0 and k1 are the incident and final wave vectors, respectively. The silicalite-1 sample was activated by heating to 773 K under oxygen flow. From the nitrogen adsorption isotherm measured at 77 K, a maximum concentration of 31 molecules per unit cell was obtained, which is in good agreement with previous experimental and computed occupancies.6 The zeolite was contained in a slab-shaped aluminum container (diameter 50 mm) connected to a gas inlet system. The N2 and CO2 loadings were determined from the equilibrium pressures and from the simulated adsorption isotherms calculated at the same temperatures. Since CO2 is relatively strongly adsorbed in silicalite-1, the loading dependence of the diffusivity could be measured at 300 K, whereas for N2 a temperature of 200 K was used to reach high loadings at a reasonable pressure (e1 atm). Both CO2 and N2 molecules exhibit coherent scattering so that the collective motions of the system are predominantly probed. For coherent scattering, the scattering function depends on the total van Hove correlation function:18
G(r,t) ) Gs(r,t) + Gd(r,t)
(12)
The self-part Gs(r,t) in eq 12 reflects the self-diffusivity term, and the distinct part Gd(r,t) correlates the motion of different atoms as a function of time. Therefore, G(r,t) gives the correlation between the densities at two different positions and different times:
1 G(r,t) ) 〈F(0,0) F(r,t)〉 F
(13)
The total correlation function is directly connected with the density fluctuations in the adsorbed phase. The transport diffusivity, Dt, can be obtained from QENS experiments because eq 13 gives the evolution of local concentration gradients around equilibrium. Provided that the length scale probed is large enough, the collective diffusion coefficient derived from experiment, D(Q), becomes the macroscopic transport diffusivity. In a recent MD simulation study of Briels and co-workers,27 the transport diffusivity was also extracted from the fluctuations in an equilibrium simulation; in the 1D system of that work, D(Q) was found to reach the macroscopic transport diffusion coefficient for Q ≈ 0.05 Å-1.
Figure 1. Intensities derived for N2 in silicalite-1 at 200 K at different loadings: (. 0.4 N2/uc (uc ) unit cell), (0) 0.85 N2/uc, (O) 1.7 N2/uc, (4) 2.6 N2/uc, (]) 3.75 N2/uc, (3) 5.0 N2/uc, (+) 6.45 N2/uc.
At low Q values, the coherent scattering function for isotropic motion has the form
Scoh(Q,ω) )
Dt(θ) Q2 S(Q) π ω2 + (D (θ) Q2)2
(14)
t
where the coherent structure factor is defined by the integral
S(Q) )
∫-∞∞ Scoh(Q,ω) dω
(15)
In eq 14, the concentration dependence of Dt has been expressed on the basis of loading (fractional occupancy) defined as θ ) F/Fm, with Fm being the maximum capacity of zeolite for the particular sorbate molecule at the temperature of interest. The first difference with incoherent scattering is that the intensity varies as S(Q), apart from a trivial Debye-Waller factor. Therefore, one expects a nonmonotonic variation of the total intensity as a function of Q. In the small Q domain, the line shape of the coherent scattering function is still Lorentzian, with an fwhm of
∆ω(Q)coh ) 2Dt(θ) Q2 At zero momentum transfer (macroscopic density fluctuations), the limiting value of the structure factor S(Q) is proportional to the isothermal compressibility, itself related to the inverse of the thermodynamic factor (eq 6b):
∂ ln p 1 ) ∂ ln F S(0) Thus, the width of the coherent signal at very low Q can also be written as
∆ω(Q)coh )
2D0(θ) Q2 S(Q)
(16)
On a time-of-flight instrument, the conversion of the scattered intensities to S(Q) is not straightforward because the measurements are made at constant scattering angles and not constant Q. Furthermore, a series of corrections have to be made, and the integration over energy transfers in eq 15 is restricted experimentally. The determination of the intensities scattered in the quasielastic domain (Figure 1) yields, however, a quantity that is related to the structure factor, after normalization with respect to the number of adsorbed molecules. This is shown in Figure 2. The number of data points at small Q values is not sufficient to determine S(0). However, one can check in Figure 2 that the inverse of the structure factor for small momentum
12752 J. Phys. Chem. B, Vol. 108, No. 34, 2004
Papadopoulos et al.
Figure 2. Structure factors derived for N2 in silicalite-1 at 200 K at different loadings (same symbols as in Figure 1). The lines indicate extrapolation to small Q values.
transfers becomes larger when the loading increases; i.e., the thermodynamic factor increases. 3.2. Simulation. Models. A perfect silicalite-1 crystal is modeled in its orthorhombic form (Pnma space group), with the oxygen and silicon atoms bearing partial charges of -1e and +2e, respectively, and placed in a rigid framework according to crystallographic data28 (model LJCB.JBTLC29 in the notation of ref 6). The nitrogen molecule was represented as a quadrupole with three partial charges on the molecular axis, which was considered rigid, and two Lennard-Jones (LJ) sites on the atoms (model 2LJ3CB.MSKM30 in the notation of ref 6). The linear triatomic molecule of carbon dioxide was modeled similarly in terms of an LJ site on each atom and three partial charges distributed along the molecular axis (model 3LJ3CB.EPM231 in the notation of ref 6). Short-range sorbate-sorbate and sorbate-zeolite interactions were calculated by a 12-6 LJ potential, and the longrange electrostatic interactions were summed by the Ewald technique. All the parameters and model abbreviations used in the present work are identical to those explained in detail in our previous detailed study of the same systems6 and will not be repeated here. Molecular Dynamics. Grand canonical Monte Carlo (GCMC) simulations were conducted to sample state points at the desired T and F. Translation, rotation, creation, and destruction trials were employed according to the Metropolis scheme, as implemented by the Adams algorithm.32 Configurations from the GCMC (up to 5 million microstates with 54-216 molecules in the simulation box) were used as initial configurations for the MD runs.6 Equilibrium MD simulations were carried out in the microcanonical (N, V, E) ensemble. The LEN algorithm due to Fincham6,32 was used for the time evolution of the linear molecules. The simulation box, after testing for system size effects, consisted of 8 unit cells in the case of GCMC simulations and 27 unit cells in the case of MD simulations. As in the work of ref 6, the larger box size in the case of MD was chosen to ensure collection of good statistics and to avoid artificial correlations in the motion as a consequence of periodic boundary conditions. The runs have been carried out on a Beowulf Linux cluster consisting of 12 Athlon-MP CPUs communicating through a gigabit switch in order to be capable of performing parallel runs. An initial run lasting about 100 ps was sufficiently long to establish the desired temperature. Long simulation times up to 20 ns were necessary in order to ensure good statistics for calculating D0 especially at high occupation densities; as seen
Figure 3. Isotherms of N2 at 200 K and CO2 at 300 K predicted by GCMC. An experimental CO2 isotherm at 300 K is also given. Comparisons of N2 results against experiment have been included in ref 6. Lines display the results of nonlinear regression of the simulated θ(ln f) dependence using the QC model for the isotherm. Optimal values of the dimensionless sorbate-sorbate interaction energy βw obtained through the regression are indicated (see also Table 1).
from a comparison of eqs 8 and 11, D0, being a collective property, suffers from inherently much poorer statistics than Ds. 4. Results and Discussion 4.1. QENS MD. Figure 3 shows the predicted isotherms from GCMC simulation. These predictions have been shown to be in excellent agreement with experiment6 and have been used as a means to determine the conditions under which containers should be filled in order to achieve the desired loadings in the QENS experiments. They were also used to go from D0 to Dt values. From Figure 4, the Arrhenius plots of transport diffusivities from QENS and simulation for both sorbates at the same occupancy are seen to be in excellent agreement. Predicted activation energy values are very close to those experimentally determined. The results for the transport diffusivity measured directly from QENS measurements and our MD predictions for N2 and CO2 are presented in Figure 5. The corrected diffusivities in the QENS graphs have been produced via Darken factors obtained from the simulated isotherms. The increasing trend for the MD transport diffusivity with occupancy for both gases is in excellent agreement with the QENS experiments, indicating a slightly concentration-dependent D0 over the sorbate loading range used in the experimental and simulation work. However, particular attention should be paid to the corrected diffusivity, which, as discussed in the theoretical part of this paper, is largely free of thermodynamic influences (eq 6b), reflecting this way only the rate at which sorbate molecules escape from the potential minima due to the atomic structure of the sorbent. Figure 5 shows that the measured corrected diffusivity of CO2 decreases as sorbate density increases, being in good agreement with the predictions from MD. In the case of N2, the QENS measurements of D0 exhibit a shallow maximum at a loading of roughly two molecules per unit cell; on the other hand, MD results show a rather occupancy-independent behavior of D0 values. Small differences between predicted and measured values of Dt and D0 can be attributed to the force fields used. The reader
Transport Diffusivity in Silicalite
J. Phys. Chem. B, Vol. 108, No. 34, 2004 12753
Figure 4. Temperature variation (Arrhenius plots; see eq 1) of the transport diffusivity of N2 (top) and CO2 (bottom) as obtained from QENS measurements (left) and MD simulations (right) under the same conditions of loading. The activation energies extracted from the plots are given for comparison.
Figure 5. Measured (QENS) and predicted (MD) transport (Dt) and corrected (D0) diffusivities in silicalite-1 shown as functions of loading.
is reminded that no parameter adjustment was implemented at any point during calculations. The dependence of D0 on the sorbate concentration is studied in detail in the next sections, in light of an existing theoretical model. 4.2. Quasichemical Approach to Adsorption. An effort to understand the dependence of D0 on θ theoretically on the basis of sorbate-sorbate interactions has been made by Reed and Ehrlich (R-E).35 It is instructive to compare our simulation results against the predictions of this theory. Since the R-E theory invokes the quasichemical (QC) approximation, we begin
by mapping our sorption thermodynamics predictions, as obtained from GCMC simulations, onto the quasichemical model. This mapping will provide us with the parameters invoked by the R-E theory. The quasichemical mean field approximation due to Guggenheim for the study of the thermodynamics of nonideal solutions (equivalent to Bethe-Peierls approximation for the Ising model of ferromagnetism)34 can be adapted to analyze a “solution” of NA adsorbate molecules distributed among a total of NS adsorption sites arranged on a simple lattice of uniform
12754 J. Phys. Chem. B, Vol. 108, No. 34, 2004
Papadopoulos et al.
coordination number z. Each lattice site is capable of accommodating at most one molecule. Each sorbate molecule interacts with its nearest neighbors. The total number of nearest neighbor pairs is denoted by NAA; the contribution from each pair of nearest-neighbor sorbed molecules to the overall potential energy is wAA. In the following, we will use the quantity w ) zwAA/2 for convenience. With these definitions, the potential energy UAA of sorbate-sorbate interactions reads
UAA ) NAA
2w z
(17)
For a given number of sorbate molecules, NA, on the Ns sites, there are
( )
Ns! Ns NA ) NA!(Ns - NA)!
arrangements on the lattice. The value of NAA, and hence the energy UAA, may vary among these configurations. We now introduce a mean field approximation wherein all configurations of NA molecules are assigned the same energy, 〈UAA〉 ) 2w〈NAA〉/z. Furthermore, we consider the system in the grand canonical ensemble at temperature T at equilibrium with a bulk fluid phase where the chemical potential of the sorbate is µ. The grand partition function Ξ(µ,V,T) of the adsorbed molecules in the lattice can be written as follows:
( )
NS
Ξ(µ,V,T) )
N
exp(-β〈UAA〉) Ns ∑ A N )0 A
(
exp(NAβµ)[qs(T)]NA (18)
2 - 2θ ζ+1
)
(19)
where θ ) NA/Ns ) F/Fm is the fractional occupancy, as above, and ζ is defined by the expression
ζ ) {1 - 4θ(1 - θ)[1 - exp(-2βw/z)]}1/2
(20)
Therefore, for the total potential energy U(θ) per adsorbate molecule the following relation is obtained:
(
)
∂ ln q (T) 2(1 - θ) β U(θ) ) T + βw 1 NA ∂T ζ+1 s
(21)
Equation 21 at infinite sorbate dilution reduces to s
∂ ln q (T) β U(0) ) T NA ∂T
(22)
From eqs 21 and 22,
U(θ) - U(0) 2 - 2θ ) βw 1 β NA ζ+1
(
CO2, T ) 300 K N2, T ) 200 K
βw
Fm (molecules/uc)
Fm (GCMC) (molecules/uc)
-0.1 0.7
24.03 19.66
23.10 20.17
a Fm and Fm (GCMC) denote the maximum loading capacity, as estimated from fitting and as calculated from GCMC simulation, respectively. uc ) unit cell.
Through the grand partition function (eq 18), the isotherm of the system can be derived, resulting in the relation
βΛ3qs(T) 2 - 2θ z θ f ) bf ) q 1 - θ ζ + 1 - 2θ
(
)
(23)
)
(24)
where q is the partition function over the internal degrees of freedom in the ideal gas, Λ is the thermal wavelength, and the parameter b is related to the Henry’s constant KH for sorption through the expression
b ) lim ff0
θ F 1 1 ) lim ) K f Fm ff0 f Fm H
Obviously, for no sorbate-sorbate interactions (w ) 0) the QC isotherm (eq 24) degenerates into the simpler Langmuir isotherm:
bf )
where qs(T) is the partition function over the internal degrees of freedom of an adsorbed molecule on the lattice, referenced to the lowest-energy internal state of the molecule in the gas phase at infinite separation. The mathematical method for estimating the mean pair energy 〈UAA〉 from eq 18 is due to Bethe,36 resulting eventually in the relation
〈UAA〉 ) wNA 1 -
TABLE 1: Results from the Fit of Eq 24 to the Simulated Isothermsa
θ 1-θ
To map our atomistic sorption thermodynamics simulation onto the QC model (Figure 3), (i) the simulated isotherms were compared to experimental results for both sorbates and are found to be in excellent agreement (for N2, isotherm predictions from GCMC were validated against experiment at 300 and 77 K in ref 6). (ii) The theoretical ln f versus θ relationship from the QC model (eq 24) was fitted to the isotherms predicted by GCMC, using βw and Fm as adjustable parameters. Results from this fitting procedure are shown in Figure 3. The estimated pair interaction energies βw and the maximum sorbate capacity values Fm resulting from this fitting, assuming a lattice of coordination number z ) 8, are presented in Table 1. From the values obtained for the two sorbates, one is clearly led to the conclusion that intermolecular interactions in the sorbed phase of N2 are weaker (less attractive) than those of CO2. Clearly, use of a constant coordination number z is an oversimplification for our systems because sorbate molecules do not arrange on a regular lattice. At saturation, each of the three environments (straight channel, sinusoidal channel, and channel intersections) of silicalite-1 contains roughly two molecules; on the basis of the spatial arrangement of these environments, one would estimate five neighbors for channel interiors and nine for channel intersections. In each case, one or two neighbors would be closer than the rest. The value z ) 8 was found to give the best fit of the QC model to the isotherms among the values z ) 2, 4, 6, 8, 10 tried. Sorbate-sorbate energetics can also be probed through the quantity [U(θ) - U(0)]/NA (total energy per sorbate molecule, referenced to the zero occupancy limit). This quantity was estimated directly as a difference of ensemble averages accumulated in the GCMC simulations carried out at occupancy θ and at very low occupancy. A theoretical expression based on the QC theory is given by eq 23. In Figure 6 is shown the quantity [U(θ) - U(0)]/NA as a function of θ, as obtained from simulation and as estimated from eq 23 for z ) 8 and various
Transport Diffusivity in Silicalite
J. Phys. Chem. B, Vol. 108, No. 34, 2004 12755
Figure 6. Energy of the sorbed phase per sorbed molecule referenced to the zero-occupancy limit as obtained from MD simulations. Lines are graphs of the energy of QC model calculated for various values of pair interaction energy w and lattice coordination number z ) 8.
values of βw. The values of βw determined from the best fit to the isotherms do not capture the simulated values [U(θ) - U(0)]/ NA. This emphasizes that the QC model is but a very crude representation of the actual zeolite systems. Nevertheless, it is clear from Figure 6 that sorbate-sorbate interactions for CO2 at 300 K are significantly more attractive than for N2 at 200 K as obtained from fitting the simulated isotherms. 4.3. Occupancy Dependence of Corrected Diffusivity. The work37 of June et al., Snurr et al., and Saravanan et al. has employed simulation studies on lattices in order to investigate the concentration dependence of Ds. As already mentioned, an attempt to relate sorbate-sorbate energetics in the adsorbed phase to the occupancy dependence of D0 has been made by Reed and Ehrlich.35 Their model has originally been developed for surface diffusion on a regular lattice of sites. It can be adapted to our zeolite sorbate systems, however, through the mapping procedure discussed in section 4.2. In the R-E work, the overall effective jump rateG(θ) of adsorbed molecules on a simple lattice is written as a function of G (i), the individual jump rate of a molecule surrounded by i adsorbate molecules at nearest-neighbor sites, and of the probability P(i) that such a molecule is surrounded by i neighbors. The latter quantity is estimated by invoking the QC approximation (section 4.2), giving
P(i) )
()
z [ζ - 1 + 2θ]i[2 - 2θ]z-i[ζ + 1]-z i
(25)
In expressing G(i), it is assumed that any attempt to jump into an already occupied site is unsuccessful. The overall effective jump rate emerges as
G(θ) ) G(0)
[2ζ-+2θ1 ] [1 + ζ -2 -1 +2θ2θ exp(2βw/z)] -z
z-1
(26)
where G (0), the jump rate on an empty lattice, is the limiting value of eq 26 as θ approaches zero. From the calculation of the flux vector of the adsorbed molecules executing jumps of constant length λ on the lattice, with the chemical potential gradient acting as a driving force, and from a comparison with eq 5, the following relation for the transport diffusivity is derived:31
Dt(θ) ) G(θ) λ2
∂ ln f ∂ ln θ
(27)
Figure 7. Normalized corrected diffusivities from MD simulation and corresponding QC theoretical curves obtained from the R-E model35 (see eq 29) for various pair interaction energy values, including those extracted from the fitting procedure of Figure 3.
In the limit of very low loadings, by combining eqs 6b, 26, and 27, one obtains the relation
Dt(0) ) D0(0) ) G(0)λ2
(28)
Combination of eqs 6b and 26-28 leads to the following expression for the occupancy dependence of the corrected diffusivity on a lattice of coordination number z.
D0(θ) ) D0(0)
[2ζ-+2θ1 ] [1 + ζ -2 -1 +2θ2θ exp(2βw/z)] -z
z-1
(29)
For w ) 0 (no sorbate-sorbate interaction apart from the requirement of no multiple occupancy of sites), eq 29 gives
D0(θ) ) D0(0) (1 - θ)
(30)
Our simulation estimates for the diffusivities D0(θ) for CO2 and N2 in silicalite-1, normalized by the corresponding D0(0) values, are presented in Figure 7. In the same figure are given the D0(θ)/ D0(0) curves obtained from the R-E model, eq 29, for the values of βw and Fm extracted from fitting the sorption thermodynamics (Table 1). Plots of eq 29 for a higher and lower value of βw have also been included to display the overall behavior predicted by the R-E model. The R-E model provides a satisfactory representation of the occupancy dependence of the corrected diffusivity. This supports the idea that the physical origin of the different D0(θ) dependences seen in N2/silicalite-1 at 200 K and CO2/silicalite-1 at 300 K lies in the higher (more attractive) sorbate-sorbate interactions in the latter system. 5. Conclusions Our predicted sorption thermodynamics in the systems N2/ silicalite-1 at 200 K and CO2/silicalite-1 at 300 K, based on GCMC simulation, is in excellent agreement with previous experimental results. Mapping the predicted thermodynamics onto the quasichemical model, which assumes adsorption on a lattice of discrete sites and invokes a mean field approximation to estimate the number of nearest-neighbor pairs, led to the conclusion that sorbate-sorbate interactions are considerably
12756 J. Phys. Chem. B, Vol. 108, No. 34, 2004 stronger (more attractive) in CO2/silicalite-1 than in N2/silicalite1. The same conclusion is supported by direct calculation of the potential energy, referenced to its zero-occupancy limit, from the GCMC simulation of the two systems. Coherent quasielastic neutron scattering measurements of the transport diffusivity Dt as a function of occupancy θ were predicted successfully by our molecular dynamics simulations. These measurements pointed out a significant difference in the occupancy dependence of the corrected diffusivity, D0(θ), between the two systems investigated. While in CO2/silicalite-1 at 300 K, D0(θ) is clearly a decreasing function, and in N2/ silicalite-1 at 200 K it eventually stays constant, exhibiting a weak maximum. Our MD simulations captured this difference satisfactorily. To understand the physical origin of the different D0(θ) dependency, we have invoked the Reed and Ehrlich35 model of surface diffusion, with parameters obtained from fitting the simulated isotherms. The model is but a crude representation of actual transport in the zeolite pores. Nevertheless, with the mapping procedure adapted in our work, it captures the different simulated D0(θ) dependencies satisfactorily. Differences in sorbate-sorbate energetics emerge as the physical reason for the different occupancy dependence of the corrected diffusivity. Strongly attractive sorbate-sorbate interactions in CO2/silicalite-1 at 300 K (βw negative in the R-E model) lead to clearly decreasing D0(θ), while weaker attraction in N2/silicalite-1 at 200 K (βw positive in the R-E model) leads to an almost constant D0(θ). Recent neutron scattering experiments in the system D2/NaX resulted in a measured D0 that is an increasing function of θ,8 whereas in the systems CF4/silicalite-1 and benzene/NaY the measured D0(θ) exhibits a clear decreasing behavior.38 These findings strengthen the aforementioned conclusions about the effect of sorbate-sorbate interactions on the corrected diffusivity. While the R-E model is successful in predicting a simple, qualitative explanation of the D0(θ) dependence, there is much incentive for the development of more refined theoretical models of D0(θ), based on a more faithful representation of zeolite/ sorbate systems such as CO2/silicalite-1 and N2/silicalite-1. Such models could involve multiple coordination numbers and interaction energy parameters depending on the local environment (e.g., straight channels, sinusoidal, interconnections) of the zeolite. They should also allow for the ability of sorbate molecules to bypass each other in certain regions (e.g., intersections) of the zeolite. Acknowledgment. The neutron experiments were performed at the Institut Laue-Langevin, Grenoble, France, using the IN6 spectrometer. We are grateful to Prof. R. Krishna for pointing our attention to the Reed and Ehrlich work. We thank Dr. H. Schober for his help during the measurements and Dr. K. Makrodimitris for some calculations. Financial support provided by the TROCAT project (Grant No. G5RD-CT-2001-00520) of the European Commission is gratefully acknowledged. References and Notes (1) Ka¨rger, J.; Ruthven, D. M. Diffusion in Zeolites and Other Porous Materials; Wiley-Interscience: New York, 1992. Rouquerol F.; Sing, K. S. W. Adsorption by Powders and Porous Solids; Academic Press: London, 1999. (2) Bearman, R. J.; Kirkwood, J. G. J. Chem. Phys. 1958, 28, 136. Bearman, R. J. J. Chem. Phys. 1959, 31, 751. Mason, E. A.; Viehland,
Papadopoulos et al. L. A. J. Chem. Phys. 1978, 68, 3562. Mason, E. A.; Lonsdale, S. J. Membr. Sci. 1990, 51, 1. (3) Auerbach, S. M. Int. ReV. Phys. Chem. 2000, 19, 155. Fuchs, A. H.; Cheetman, A. K. J. Phys. Chem. B 2001, 105, 7375. (4) Theodorou, D. N.; Snurr, R. Q.; Bell, A. T. In Molecular Dynamics and Diffusion in Microporous Materials; Alberti, G., Bein T., Eds.; Comprehensive Supramolecular Chemistry, Vol. 7; Pergamon: New York, 1999. (5) Maginn, E. J.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1993, 97, 4173. (6) Makrodimitris, K.; Papadopoulos, G. K.; Theodorou, D. N. J. Phys. Chem. B 2001, 105, 777 and references therein. (7) Frenkel, D.; Smit, B. Understanding Molecular Simulation, 2nd ed.; Academic Press: San Diego, CA, 2002. (8) Jobic, H.; Ka¨rger, J.; Bee, M. Phys. ReV. Lett. 1999, 82, 4260. (9) Magda, J. J.; Tirrell, M.; Davis, H. T. J. Chem. Phys. 1985, 83, 1888. MacElroy, J. M. D.; Suh, S.-H. Mol. Phys. 1987, 60, 475. (10) Cracknell, R. F.; Gubbins, K. E. Langmuir 1993, 9, 824. (11) Steele, W. A. Surf. Sci. 1973, 36, 317. Steele, W. A. The Interaction of Gases with Solid Surfaces; Pergamon Press: Oxford, 1974. (12) Papadopoulos, G. K.; Suh, S.-H.; Nicholson, D. Mol. Simul. 1999, 22, 237. (13) Ka¨rger J.; Pfeifer, H. Zeolites 1987, 7, 90. (14) Makrodimitris, K.; Papadopoulos, G. K.; Philippopoulos, C.; Theodorou, D. N. J. Chem Phys. 2002, 117, 5876. (15) Papadopoulos, G. K. J. Chem. Phys. 2001, 114, 8139. (16) Ruthven, D. M.; Derrah, R. I. J. Chem. Soc., Faraday Trans. 1 1972, 68, 2332. (17) Papadopoulos, G. K.; Petropoulos, J. H. J. Membr. Sci. 1995, 101, 27. (18) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids; Academic Press: London, 1986. Reichl, L. E. A Modern Course in Statistical Physics; Edward Arnold: Kent, U.K., 1991. (19) Papadopoulos, G. K.; Petropoulos, J. H. J. Chem. Soc., Faraday Trans. 1996, 92 (17), 3217. (20) Visscher, V. M. Phys. ReV. A 1974, 10, 246. Evans, D. J.; Morriss, G. P. Phys. ReV. A 1988, 38, 4142. Petravic, J.; Evans, D. J. Phys. ReV. Lett. 1997, 78, 1199. (21) Arya, G.; Maginn, E. J.; Chang, H.-C. J. Chem. Phys. 2000, 113, 2079. (22) Heffelfinger, G. S.; van Swol, F. J. Chem. Phys. 1994, 100, 7548. MacElroy, J. M. D. J. Chem. Phys. 1994, 101, 5274. (23) Arya, G.; Chang, H.-C.; Maginn, E. J. J. Chem. Phys. 2001, 115, 8112. (24) Skoulidas, A. I.; Sholl, D. S. J. Phys. Chem B. 2002, 106, 5058. Skoulidas, A. I.; Sholl, D. S.; Krishna, R. Langmuir 2003, 19, 7977. (25) Sanborn, M. J.; Snurr, R. Q. Sep. Purif. Technol. 2000, 20, 1. (26) Skoulidas, A. I.; Ackerman, D. M.; Johnson, J. K.; Sholl, D. S. Phys. ReV. Lett. 2002, 89, 185901. (27) Hoogenboom, J. P.; Tepper, H. L.; van der Vegt, N. F. A.; Briels, W. J. J. Chem. Phys. 2000, 113, 6875. Tepper, H. L.; Briels, W. J. J. Chem. Phys. 2002, 116, 9464. (28) Olson, D. H.; Kokotailo, G. T.; Lawton, S. L. J. Phys. Chem. 1981, 85, 2238. (29) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1990, 94, 1508. (30) Murthy, C. S.; Singer, K.; Klein, M. L.; McDonald, I. R. Mol. Phys. 1980, 41, 1387. (31) Harris, J. G.; Yung, K. H. J. Phys. Chem. 1995, 99, 12021. (32) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (33) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (34) Fowler, R.; Guggenheim, E. A. Statistical Thermodynamics; University Press: Cambridge, 1949; p 421. (35) Reed, D. A.; Ehrlich, G. Surf. Sci. 1981, 102, 588. (36) Bethe, A. Proc. R. Soc., Ser. A 1935, 150, 552. (37) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1991, 95, 8866. Snurr, R. Q.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1994, 98, 5111. Saravanan, C.; Jousse, F.; Auerbach, S. M. J. Chem. Phys. 1998, 108, 2162. (38) Jobic, H.; Skoulidas, A. I.; Sholl, D. S. J. Phys. Chem. B 2004, 108, 10613. Jobic, H. et al., in preparation. (39) Jobic, H.; Makzodimitris, K.; Papadopoulos, G. K.; Schober, H.; Theodorou, D. N. In Proceedings of the 14th International Zeolite Conference; van Steen, E. et al., Ed.; Cape Town, 2004; p 2056.
