8232
J . Phys. Chem. 1990, 94, 8232-8240
Molecular Dynamics Study of Methane and Xenon in Silicalite R. Larry June, Alexis T. Bell, and Doros N. Theodorou* Department of Chemical Engineering, University of California, Berkeley, Berkeley, California 94720 (Received: December 18, 1989; In Final Form: May 11, 1990) The dynamic and equilibrium properties of two different sorbate molecules, xenon and methane, in the pentasil zeolite silicalite were studied with molecular dynamics simulations. A rigid polyatomic model was employed for methane, while xenon was represented as a Lennard-Jones sphere. The zeolite lattice was assumed rigid. The effect of intracrystalline occupancy and temperature on the structure and dynamics of the sorbed phase was investigated. Predicted self-diffusivities are in good agreement with N M R measurements. Except at the lowest temperature studied, 200 K, the sorbate self-diffusivities were found to decrease monotonically with occupancy. At 200 K, the self-diffusivities exhibit a weak maximum before decreasing at high loadings. This maximum is attributed to the trapping of sorbate molecules in small potential wells formed within the pore network of the zeolite at low occupancy. The single particle density distribution function, p'(r), shows that at low loadings sorbate molecules are localized preferentially in the sinusoidal channels and avoid the energetically less favorable channel intersections and straight channels. At higher loadings, partial order is observed in the intracrystalline fluid. Introduction Diffusion of guest molecules and atoms inside the pores of zeolite crystals can have a significant effect on the performance of zeolites as sorbents and catalysts. For example, for very rapid chemical reactions, the rate of reaction within the zeolite may be severely limited by the rate of intrazeolite sorbate diffusion. The intracrystalline diffusivity depends on the size and shape of the sorbate molecule relative to that of the zeolite pore, the energetic interactions between the zeolite and the sorbate molecule, and the connectivity and dimensionality of the zeolite pore structure. Since the experimental determination of sorbate diffusivities is time-consuming, it would be highly desirable to be able to predict the transport rates of sorbate molecules a priori, from fundamental knowledge of the structure and energetics of a zeolite-sorbate system. A promising means for arriving at this goal is molecular dynamics (MD) simulation. Such simulations can also be used to explore characteristics of sorbate motion within the zeolite pores and to determine thermodynamic properties. Several reports have appeared in the literature concerning the use of M D simulations to characterize the motion of atoms and molecules in idealized porous media, as well as zeolites. Diffusion of a hard-sphere fluid contained within pores of width less than I O molecular diameters was first studied by Subramanian and Davis.' Diffusivities were found to decrease by a factor of 2 when reflection by the walls was changed from specular to diffuse. A liquid of Lennard-Jones particles was simulated in structureless slit pores of varying width by Magda et al., to investigate sorbate structure and self-diffusion.* The density profile was found to be oscillatory near the walls of the pore. In spite of the density variations within the pore, the calculated diffusivities were found to be relatively insensitive to position within the pore. MacElroy and Suh3 have compared the behavior of a monatomic Lennard-Jones fluid in a cylindrical pore with and without a structured surface, to investigate the effect of the wall architecture on sorbate diffusivity. Profound differences were found in the axial diffusivity of the sorbate depending on the roughness of the wall. Schoen et al. have reported similar results, but for a slit pore ~ y s t e m . ~ M D simulations of sorbate motion within zeolites have been discussed very recently by several groups. Trouw and Iton, Cohen de Lara et al., and Yashonath et al. have studied the motion of ,~ and sodium-A.' More commethane in f a ~ j a s i t esilicalite,6 ( I ) Subramanian, G.; Davis, H. T. Mol. Phys. 1979, 38, 1061-1066. (2) Magda, J. J.: Tirrell, M.; Davis, H. T. J . Chem. Phys. 1985, 83,
1888-1901. (3) MacElroy. J . M. D.; Suh, S.-H. Mol. Simulat. 1989, 2 , 313-351. (4) Schoen, M.; Cushman, J. H.; Diestler, D. J.; Rhykerd, C. L. J. Chem. Phys. 1988,88, 1394-1406. (5) Yashonath. S.: Demontis. P.; Klein, M. L. Chem. Phys. Left. 1988, 153, 551-556. (6) Trouw, F. R.; Iton, L. E. In Zeolites for the Nineties, Recent Research Reports; Jansen. J. C., Moswu, L., Post, M. F. M., Eds. Presented at the 1989 International Zeolites Conference, Amsterdam, The Netherlands, June 1989; pp 309-3 I O .
0022-3654/90/2094-8232$02.50/0
plicated systems have been addressed by LehertC et who investigated the dynamics of water in ferrierite, and by Demontis et al.,9 who studied the diffusion of benzene in sodium-\(. These simulations of sorbed molecules have produced estimates of the self-diffusivities that are in good agreement with experiment and have provided a large amount of detailed information regarding the microscopic pathways by which diffusion occurs within the zeolite pore network. In addition to the dynamical information gathered from M D simulations, Monte Carlo (MC) techniques have proven useful for elucidating thermodynamic properties, such as heats of sorption and sorption isotherms, in zeolitic systems. The grand canonical ensemble M C simulations of Soto and Myerslo and Woods and Rowlinson" have accurately predicted isotherms for small molecules in various faujasite structures. A recent work by Leggetter and Tildesley'* suggests that MC may be more efficient than MD in gathering thermodynamic properties. In this paper we report on the dynamic, structural, and thermodynamic properties of xenon and methane in silicalite over a range of occupancies and temperatures, as obtained by M D calculations. Experimental data, available for methane and xenon, provide an excellent test of the predictive power of M D simulations and of the assumptions invoked in the molecular representation. In the case of xenon, the simplicity of the model system enables the study of long-time dynamics (simulations exceeding 1000 ps). Calculated diffusivities for both systems compare well with those determined experimentally by N M R pulsed-field gradient techniques. The predicted structure and thermodynamic properties of the intracrystalline fluid at low occupancy are consistent with our earlier work.I3 The present study demonstrates that, with increasing occupancy, changes occur in the spatial organization and in the mechanism of motion of sorbate molecules. Model Representation: Structure and Interaction Potential The crystal lattice of silicalite is isostructural with ZSM-5 and contains virtually no aluminum atoms. The lattice parameters for the silicalite unit cell have been determined from X-ray diffraction studies and are reported as 20.07, 19.92, and 13.42 A in the Pnma space group (orth~rhombic).'~A change in sym(7) Cohen de Lara, E.; Kahn, R.: Goulay, A. M. J . Chem. Phys. 1989, 90, 7482-7491. (8) Lehertb, L.; Varcauteren, D. P.; Derouane, E. G.; Lie. G. C.; Clementi, E.; Andre, J.-M. In Zeolites: Facts, Figures, Future; Jacobs, P. A,, van Santen, R. A., Eds.; Elsevier: Amsterdam, 1989; pp 773-783. (9) Demontis, P.; Yashonath, S.; Klein, M. L. J. Phys. Chem. 1989, 93, 5016-501 9. (10) Soto, J. L.; Myers, A. L. Mol. Phys. 1981, 42, 971-983. ( I 1) Woods, G. B.; Rowlinson, J. S. J . Chem. Soc., Faraday Trans. 2 1989, 85. 765-781. ( I 2) Leggetter, S.; Tildesley, D. J. Mol. Phys. 1989, 68, 5 19-546. (13) June, R. L.; Bell, A. T.; Theodorou, D. N. J . Phys. Chem. 1990, 94, 1508-1516. (14) Olson, D. H.; Kokotailo, G.T.; Lawton, S. L.; Meier, W. M. J . Chem. Phys. 1981, 85. 2238-2243.
0 1990 American Chemical Society
The Journal of Physical Chemistry, Vol. 94, No. 21, 1990 8233
Methane and Xenon in Silicalite
TABLE I: van der Waals Radii ( r o ) ,Effective Number of Electrons (ne),and Atomic Polarizabilities (a)Used in Modeling All Interactions in the Silicalite-Methane and Silicalite-Xenon Systems atom
fl, A
ne
0 C H
1.575 1.8 1.3 2.32
7.0 5.0 0.9 32.4
Xe
ff,
A'
0.85 0.92 0.42 4.01
forces, B . . in eq 1 is calculated by the Slater-Kirkwood relationshipIr 3e2a01~2aiaj B..= (cgs units) (2) V (a/ne)i'/2 (a/ne)j'/2
+
while the short-range repulsion coefficient, A,, is obtained as
/.
d
/
>i
1-1 a
= 2007
A
6'
Figure 1. An outline of the pore network of the pentasil zeolite silicalite. The straight and sinusoidal channel systems are directed along the Yand
X directions, respectively. metry of the unit cell from orthorhombic to monoclinic can be induced by the replacement of some of the silicon in the lattice with aluminum, by the sorption of large sorbate molecules, or by For the simulations varying the temperature of the performed here, only the phase with orthorhombic symmetry was considered. The silicalite structure contains two sets of interconnected pores that are slightly elliptical in cross section and about 5.5 8, in diameter. One set of pores, directed along the [OIO] axis, follows a linear path through the crystal. A second set of pores follows a zigzag path, the average direction of which is along the [loo] axis. Motion in the [OOI] direction is also possible through a tortuous sequence of successive displacements along different segments of the interconnected pores. Figure 1 presents a schematic outline of the pore network, which emphasizes the connectivity and topology of the silicalite structure, as well as the coordinate system utilized in this work. For the simulations performed here, the silicalite crystal is assumed to have a perfect lattice with all framework atoms rigidly fixed in space and to contain no defects or aluminum atoms. The simulation box consists of a group of unit cells. Periodic boundary conditions were used along all three coordinate directions at positions corresponding to the unit cell boundaries. Depending upon the desired loading of sorbate molecules in a given simulation, two to four unit cells were used along each of the principal coordinate axes to form the simulation box. The simulation box size was chosen so that the total number of sorbate molecules remained fixed at a number that was computationally tractable and at the same time allowed relevant density fluctuations to occur. Structurally, methane was modeled as a five-center LennardJones molecule. Its bond angles and bond lengths were constrained to 109.5' and 1.1 A, respectively. Considered as a rigid body, our methane model has six degrees of freedom (three translational and three rotational). Xenon was represented as a single Lennard-Jones center with only three degrees of freedom. All interactions are treated atomistically by using a LennardJones potential to describe the interaction energy between any pair of atoms:
(3) Aij = '/zBij(r: + r:)6 where rij is the distance between interacting centers i and j , ai the atomic polarizability of atomic species i, e the charge of an electron, a, the Bohr radius, and n, an effective number of electrons in an atom. The expression for the repulsion coefficient, eq 3 , is determined by setting the force, -VU, between two atoms equal to zero at a separation equal to the sum of the van der Waals radii for the pair (r: and r:). Coulombic interactions are neglected, as the zeolite is assumed to be completely siliceous and the sorbate molecules under investigation are nonpolar. In zeolites with high aluminum contents, forces arising from electrostatic interactions involving permanent or induced charges in the sorbate molecules are significant and cannot be neglected. Examples of simulations in such systems are found in refs 5 and 7-9. To reduce the number of potential parameters, it was assumed that the framework silicon atoms do not interact significantly with the sorbate molecule^.'^ This assumption is justified because the sorbate molecule is surrounded by a latticework of framework oxygen atoms; silicon atoms, located at the centers of the SiOo tetrahedra, are recessed from the outer edges of the pores. By making the indicated assumption, all uncertainty in the lattice interaction potentials is absorbed into the oxygen parameters. All of the potential parameters were taken from the literature. The polarizability of the oxygen atoms was taken from Kiselev.19 The effective number of electrons and an optimal van der Waals radius of oxygen were taken from our previous work on the thermodynamics of the low-occupancy sorption of alkanes in silicalite.13 Values of the parameters ne, cy, and fl for carbon and hydrogen were taken from the work of Theodorou and SuterZ0 on glassy polymer systems. For xenon, the polarizability, a,was taken from Kiselev.19 The van der Waals radius, f', was determined from the Lennard-Jones u reported by Hirschfelder et using the well-known expression
p=
2116g
(4)
An effective number of electrons for xenon was obtained by
equating the Slater-Kirkwood expression for Bij to the Bi.derived from the values of t and cr listed by Hirschfelder et aLzl botential parameters for the sorbate atoms and the lattice oxygen atoms are presented in Table I. It should be noted that the parametrization for the methane-silicalite system is identical with that used in our previous work.I3 The total potential energy of the model system is determined by summing all interactions between pairs of sorbate molecules, as well as all sorbate interactions with the oxygen atoms of the zeolite lattice:
u = ljlattice + v a i r
(5)
The characteristic coefficient of long-range London dispersion
( 1 8) Bondi, A. Physical Properties of Molecular Crystals, Liquids, and Glasses; Wiley: New York, 1968. (19) Kiselev, A. V.; Lopatkin, A. A.; Shulga, A. A. Zeolites 1985, 5,
(15) Klinowski, J. T.;Carpenter, T. A.; Gladden, L. F. Zeolites 1987, 7 , 73-78. (16) Fyfe, C. A.; Kennedy, G. J.; Kokotailo, G. T.; Lyerla, J. R.; Fleming, W. W. 1. Chem. SOC.,Chem. Commun. 1985, 740-742. (17) Hay, D. G.; Jaeger, H. J . Chem. Soc., Chem. Commun. 1984, 1433.
26 1-267. (20) Theodorou, D. N.; Suter, U. W. Macromolecules 1985, 28, 14571478; 1986, 19, 139-154, 319-387. (21) Hirschfelder, J. 0.; Curtiss, C. F.; Bird, R. B. Molecular Theory of Gases and Liquids; Wiley: New York, 1964.
8234
The Journal of Physical Chemistry, Vol. 94, No. 21, 1990 N
.”
=
aJ(rIi)
Oatlice
i=l/E:
(7)
In eqs 5-7, the index I moves over all interacting oxygen atoms in the lattice and the indexes i a n d j move over all interacting pairs of sorbate atoms. To reduce computational effort, a cutoff radius of 13 A was employed in the evaluation of both the sorbatesorbate and the sorbate-zeolite potential summations. The error in the mean lattice energy, L/latticc,due to potential truncation, was estimated to be -0.2 kJ/mol by assuming that the oxygen atoms at a radius greater than the cutoff distance form a structureless background of uniform density equal to the mean oxygen density of the zeolite lattice. The assumption of a rigid zeolite lattice allows the sorbate atom-zeolite potentials to be evaluated once and pretabulated on a fine three-dimensional grid of 0.2-A spacing over the asymmetric unit of the unit ce11.I3 The potential energy between a sorbate atom and the zeolite lattice was determined rapidly and accurately by interpolation from the pretabulated information using a three-dimensional cubic hermite polynomial interpolation algorithm.22 The interpolation algorithm is different from the three-dimensional spline23used in our previous work.I3 The advantages of the new algorithm, which was designed specifically for M D simulations, include the use of analytical derivatives in fitting the hermite polynomial basis functions, continuity of derivatives across boundaries of the asymmetric unit of the unit cell. and improved vectorization of the interpolation code. Calculations The simulations for both methane and xenon were carried out in the standard M D NVE ensemble. The M D method creates a set of classical trajectories by integrating the equations of motion for a group of molecules within a fixed volume of zeolite. The equations of motion were formulated in Cartesian coordinates for all sorbate atoms:
+
mi: = FNewton Fc
The constraint forces, F,, arise from the fact that each methane molecule is constrained to a prescribed geometry. The formulation of eqs 8 and 9 allows the equations of motion to be integrated in Cartesian coordinates. For xenon, constraint forces are not employed. The constraint force, F,, for methane was accounted for by implementing the constraint algorithm of Edberg, Evans, and morris^.*^ The resulting equations of motion were integrated with a fifth-order GearZ5predictor-orrector algorithm. Standard Newtonian equations of motion for xenon were integrated with the Verlet algorithmz6 A Verlet neighbor list was employed in the evaluation of forces in the xenon simulations. Following the fast rotational motion of the hydrogen atoms in the methane molecule required a I-fs simulation time step, while for xenon a 40-fs time step was adequate. The time steps were adjusted in both systems so that energy drift was limited to 0.01 5%. Positions and velocities were recorded to magnetic tape every 50 fs for methane and every 200 fs for xenon. Analysis of the simulation data was performed by a postprocessing code at a later time. The total duration of the equilibrium portion of a run over which dynamical information was recorded was 1000 ps for xenon and 50 ps for methane. CPU requirements for the highest occupancy (22) Schultz. M . Spline Analysis; Prentice-Hall: Englewood Cliffs, NJ, 1973. (23) Sathyamurthy, N.; Raff, R. M . 1D/2D/3D (QCPE 322). Department of Chemistry, Oklahoma State University, Stillwater, OK. (24) Edberg, R.: Evans, D. J.; Morriss, G . P. J . Chem. Phys. 1986, 84, 6933-6939. (25) Gear, C. W . Numerical Initial Value Problems in Ordinary Differential Equotions; Prentice-Hall: Englewood Cliffs, NJ. 1971. (26) Verlet, L. Phys. Rec. 1967. 159. 98-103.
June et al. simulations vaned from about 100 min on a Sun 4/330 workstation for the xenon system to about 210 min on a Cray X / M P for the methane system. Sorbate loadings were set at 0, 4, 8, 12, and 16 molecules/unit cell by varying the number of unit cells contained within the simulation box while maintaining the number of sorbate molecules at either 96 or 128. Infinite dilution was simulated by turning off sorbatesorbate interactions entirely, so that sorbate molecules could move independently of each other. Modeling the infinitedilution state point in this manner has advantages in that the sampling of phase space is computationally more efficient than that achieved by using a series of single molecule runs.’ Temperatures of 200, 300, and 400 K were studied, corresponding to the range over which experimental data on diffusion are available. A good initial configuration is very important for reducing the time spent in equilibrating the system by MD. Here, the initial configuration was obtained by random placement of the molecules within the pores of the zeolite followed by 100 Metropolis Monte Carlo displacement moves per molecule to prevent excessively large accelerations in the initial steps of the integration of the equations of motion in the MD simulation. Center-of-mass (COM) velocities were then randomly assigned to the molecules from a Maxwell-Boltzmann distribution. Further equilibration of the system, combined with velocity scaling to achieve the desired run temperature, was carried out for 3000-6000 integration time steps with the M D algorithm before any data were recorded. Equilibration in the methane-silicalite system was evidenced by the equality of two different measures of temperature within the system: N d
“~lr
CmFtevJr2 l=l J=I
Tatomc =
kb(6Nmole)
(10)
In the definition of eqs 10 and 11, the index i runs over all molecules, index j runs over all atoms within a molecule, and kb is the Boltzmann constant. Each Lennard-Jones site within the methane molecule has a velocity vji and a mass mytc;the total mass of a molecule is mmole.Methane molecules have a welldefined COM, which corresponds to the position of the carbon atom; vi represents the velocity of the center of mass. The total number of sites within a molecule is nsite,and the total number Equivalence of the of molecules in the simulation box is NmOle. two temperatures, for methane, shows that kinetic energy has been properly distributed among the translational and rotational degrees of freedom. In contrast to bulk fluid simulations, the total momentum of the sorbate phase is not conserved, due to the potential field exerted by the rigid zeolite lattice. However, at the beginning of each simulation, care was taken to ensure that the initial velocity distribution did not give rise to a net drift of molecules. The total momentum of all sorbate molecules fluctuated about zero over the course of a simulation run. Results The dynamical properties of the two sorbate-zeolite systems are remarkably similar; the only differences lie in the time scales associated with the different molecular masses and interaction forces of methane and xenon. The size and shape of the two molecules are almost identical, resulting in comparable trends in the concentration and temperature dependence of the diffusivities. Sorbate self-diffusivitiescan be computed from the trajectories created by the MD simulations from either the Einstein or the Green-Kubo relationships. For simulations of limited duration, however, the Einstein relationship2’
Methane and Xenon in Silicalite
The Journal of Physical Chemistry, Vol. 94, No. 21, 1990 8235
Diffusivity Methone/Silicolite 2.5
7.
I
I'otentinl Contour - XZ Projection
I N M R Doto of Coro, et 01. * T = 300 K N M R Doto of Jobic. et 01. * T = 200 K Simulotion = 400 K = 300 K = 200 K -
0-OT O-OT
I
2
' l &-AT o
1.0
v
n
0.5 -
,
0.0 X Potential Contour - YZ Projection
1.0 I
I A
0.8
I
In
T=293 K (NMR Dato of Heink, et 01.) 0-OT = 200 K A-AT = 300 K 0-UT = 400 K
\ Figure 4. Contours of the zeolite-xenon potential surface in two planes
'A
v)
a
0.2
-
0.0
0
4
8
12
16
20
molecules/unit cell
Figure 3. Computed self-diffusivities for xenon in silicalite compared with very recent N M R measurements of Heink et aLZ9 Predicted selfdiffusivities show a weaker dependence on concentration than do the experimental data. Overall agreement is good, especially in the highoccupancy region where diffusivities plateau.
provides a more reliable (i.e., less error prone) picture of long-time diffusive behavior than the Green-Kubo relationship. For each simulation run, self-diffusivities were calculated by monitoring mean-square molecular displacements along each coordinate axis as a function of elapsed time from a series of statistically uncorrelated time origins. The ensemble-averaged mean-square displacement in the Z coordinate is given by
where Nmolcis the number of molecules and No is the ~ ~ m bofe r independent time origins used in the run. Similar expressions can be written for the X and Y coordinates. Methane diffusivitieswere based on center-of-mass displacements. The time origins for methane were spaced by 0.5-ps intervals; 2 . 0 - p intervals were used for xenon. Figures 2 and 3 show the Dredicted self-diffusivities for methane and lenon, respectively, as'a function of sorbate occupancy and temperature. Also plotted in the figures are the experimental data obtained by pulsed-field gradient NMR s p e c t r o s ~ o p y . ~ *For -~~ (27) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1987. (28) Caro, J.; Biilow, M.; Schirmer, W.; Kiirger, J.; Heink, W.;Pfeifer, H.; Zdanov. S . J . Chem. Soc., Faraday Trans. 1985,8/, 2541-2550.
passing through the unit cell. In (a) the XZ plane that bisects the sinusoidal channels is contoured, and in (b) the YZ plane that bisects the straight channels is contoured. The potential is lowest in the central region of each of the channels segments and highest at the channel intersections. The minimum contour level is at -22 kJ/mol, and the contour lines are drawn in 2 kJ/mol increments.
equivalent conditions, the predicted values of the self-diffusion coefficient are in good agreement with the experimental values. The values of the self-diffusivity of methane predicted by Trouw and Iton are slightly smaller than those calculated here, but still in reasonable agreement with experiment.6 It should be stressed that no adjustable parameters are involved in our simulations. Also, it should be noted that the observed deviations between theory and experiment cannot be viewed as a failure of the simulation model; uncertainties in the experimental data are as high as 50%. Figures 2 and 3 show that at 400 K the diffusivities for methane and xenon decrease monotonically with increasing sorbate occupancy of the zeolite. This occupancy dependence is weaker at 300 K, and at 200 K a weak maximum is observed in the predicted self-diffusivity. At 200 K and an occupancy of about 4 molecules/unit cell, a similar maximum has be& okerved very recently in the diffusion coefficients derived from frequency response measurements by v ~ ~ - D and ~ co-workers ~ - B ~in the ~ meth~ ~ ane-silicalite s y s t e m . 3 ~ The maxima in the diffusion coefficients predicted at 200 K can be explained as follows. Differences in sorbate-lattice interaction energy of approximately 4-6 kJ/mol are calculated (29) Pfeifer, H.; Karger, J. Personal communication, Karl Marx University, Leipzig, GDR. (30)Jobic, H.; BCe, M.; Caro, J.; KBrger, J. In Zeolites f o r the Nineties, Recent Research Reports; Jansen, J. C., Moscou, L., Post, M. F. M., Eds. Presented at the International Zeolite Conference, Amsterdam, The Netherlands. June 1989: DD 305-306. (31)'Van-Den-Be'&, N.; Rees, L. V. C.; Caro, J.; Btilow, M. Zeolites 1989, 9, 312-317.
8236 The Journal of Physical Chemistry, Vol. 94, No. 21, 1990 between the channel segments and the channel intersections, the former being energetically more favorable for the sorbate molecules. Potential energy contours for a single xenon atom in silicalite at infinite dilution in the XZ and YZ planes are presented in Figure 4. a and b, respectively. It can be seen that channel intersections present potential barriers to diffusion and that the barriers are higher in the sinusoidal channels. At low temperature and low occupancy, a significant portion of the sorbate molecules will be “trapped” in individual channel segments, and consequently, only those atoms or molecules possessing translational energies in excess of the potential barriers will be able to move freely between channel segments. With increasing temperature, the fraction of the sorbate molecules having energies in excess of the barrier increases, and diffusion through the zeolite pore space becomes more facile. As the concentration of sorbate molecules increases, intermolecular collisions act to push a sorbate molecule from a channel interior past the intersection and into an adjacent channel segment. Collisions between sorbate molecules occupying the same channel segment have two effects: the first is an exchange of energy during a collision, and the second is a reduction of mean free path for motion along the axis of the channel. The first of these effects leads to an increase in the rate of barrier crossings between channel segments, whereas the second effect leads to a decrease. The degree to which the first or the second effect dominates depends on temperature and sorbate occupancy. Reference to Figures 2 and 3 indicates that, at 200 K and occupancies below 4-8 sorbate molecules/unit cell, the effects of energy exchange are stronger than the effects of scattering collisions between sorbate pairs, and so the diffusivity increases with increasing sorbate occupancy. At higher loadings and temperatures, the effects of scattering collisions dominate and act to reduce the sorbate self-diffusivityat high occupancies. It should be noted, however, that our simulation does not allow for energy exchange during a sorbate molecule-wall collision. This enhances the probability of molecules being “trapped“ within the channel segments at low temperature and low occupancy. In reality, some energy exchange between the “trapped” molecule and the zeolite lattice will occur, and hence, the extent of trapping might, in fact, be less than that predicted here. The apparent activation energy for diffusion, determined from linear regression on a plot of In D,vs TI,is concentration dependent. In the infinte-dilution limit, it is about 5.6 kJ/mol for methane and 6.6 kJ/mol for xenon. The apparent activation energy decreases and eventually plateaus at higher occupancies, reaching values of 1.8 and 2.8 kJ/mol for methane and xenon, respectively. at I6 molecules/unit cell. Experimental values for the apparent activation energy of methane are in the range 4-5 kJ/mol.jO The reduction in the apparent energy barrier with increasing occupancy can be attributed to the role of sorbatesorbate collisions. Additional information regarding the dynamics of sorbed methane or xenon can be obtained by separately examining motion along the X, Y, and Z directions of the lattice. Figure 5 presents the mean-square displacements of a methane molecule along each of the coordinate axes at infinite dilution and a temperature of 400 K. The linearity of the mean-square displacements at long times indicates that the diffusional process within the zeolite is well-described by Fickian phenomenology. It can be seen in Figure 5 that the X component of diffusivity is only slightly larger than the Y component, while the component along the Z axis is much smaller. As seen in Figure 1, there is no direct path formed by the channel system in the Z direction and, as a result, diffusivities are much smaller in this direction. It follows that diffusive motion projected along the X and Y coordinate axes should be faster than along the Z axis, particularly at low intracrystalline concentrations, where sorbate-sorbate collisions that would scatter diffusing sorbate molecules are absent. All motion in the Z direction is necessarily accompanied by many collisions between the diffusing molecules and pore walls. At higher concentrations, however, sorbate-sorbate collisions become more common; this lowers the diffusivity along the X and Y coordinate axes and significantly reduces the anisotropic nature of diffusion.
June et al. Mean Squore Displocement
- Methone/Silicolite
80 /
0 molecuies/unit cell; 400 K
/
/ /
,
4
0
...
20
/
. .
i 0
4
/ /
/
0 2
0
4
6
Time (ps)
Figure 5. Computed mean-square displacement for methane in silicalite at 400 K in the infinite-dilution limit. Displacements along the X and Y axes are similar; displacements along the 2 coordinate axis are much smaller. This is because there is no channel system directed along the 2 coordinate.
Mechanistic aspects of the sorbate motion are most easily explored by examining the velocity autocorrelation function (VACF) and its Fourier transform. The VACF is calculated as
1
(cZ(t) u Z ( 0 ) ) =
No N d
E C N m o ~ f l o(,,=I i = l
of(t+t,) uf(to)
(15)
where the averaging takes place over time origins, to, spaced at a I-ps interval and over all molecules, Nmolc.For methane, the center-of-mass velocity was utilized for of. Similar expressions are employed in the X and Ydirections. To examine the frequency components of the VACF, a Fourier transform is performed. The Fourier transform of a real and even function, such as the VACF, reduces to a cosine transform2’ over real frequencies, v
and similarly for the X and Y directions. The spectra Ctu(v), ~ J u ) and , C&(v) were computed by using the Fourier transform algorithm from the subroutine FILONC provided in ref 27. The dominant frequencies, or time scales, of sorbate motion can be identified from the Fourier transform of the VACF. Figure 6, a, b, and c, presents the Fourier transform of the VACF for methane in silicalite at loadings of 0, 8, and 16 molecules/unit cell, respectively. Two peaks are visible in the infinite-dilution spectrum. The approximate length scales of motion corresponding to each of these two peaks can be determined assuming that the root-mean-square molecular velocity component in any coordinate direction is given by = ( k b T / m m o l , ) 1 / 2i ;= x, y , z
(17)
The low-frequency peak, located at about 0.4 ps-], is associated with motion in the Xand, to a lesser extent, the Ydirection. This peak corresponds to long-range motion occurring over length scales on the order of 8 A, very near the IO-A length of a straight or sinusoidal channel segment. As seen in Figure 4a,b, the interior of the channel segments is energetically more favorable than the channel intersections by at least 4 kJ/mol. At low occupancy, sorbate molecules tend to “shuttle” between channel intersections in either channel segment type; however, as is apparent from Figure 6a, shuttling motion is most dominant in the X direction. The second peak of the infinite-dilution spectrum, at a frequency of about 1.25 ps-’, corresponds to motion over a length scale of about 2.5 A; it originates in a “rattling” motion of sorbate molecules between the walls of the zeolite pore. The low-frequency peak is occupancy dependent; as the sorbate concentration within the zeolite rises, the nearest-neighbor distance between sorbate
The Journal of Physical Chemistry, Vol. 94, No. 21 1990 8237
Methane and Xenon in Silicalite
~
Rototionol Diffusivity - Methone/Silicalite
Fourier Transform Of CVv(7) - Methone/Silicolite 2.5
0 molecules/unit cell; 2 0 0 K
O-OT
-x
= 400 K
-
--Y
Z
-
v
n
I
I /
I
I
1
1.5
D,= sec-1 for liq CH4 Q 98 K 0.0
0.5
10
1.5
2.0
2.5
Frequency ( p s - ’ )
b
12
16
Figure 7. Rotational diffusivity of methane in silicalite at 400 K. The rate of rotational decorrelation increases as occupancy increases due to the influence of sorbate pair collisions.
8 molecules/unit cell; 2 0 0 K
0.4
3 0.3 v
> 0 ’
0.2
0.1 0.0
0.0
I
8
4
molecules/unit cell
Fourier Tronsform Of C V V ( ~ )- Methone/Silicolite
0.5
1.0 4 0
0.5
1 .o
1.5
Fourier Transform Of CVv(7)
-
2.0
2.5
)
Frequency (ps-
the spectra at high occupancy reflects the fact that sorbatesorbate pair collisions dominate the dynamics of the intracrystalline fluid. On the contrary, at low occupancy each channel system has a distinct spectrum, indicating that sorbate-zeolite collisions dominate the dynamics of the sorbate molecules. The rotational mobility of a methane molecule can be determined from the decorrelation in orientation of a unit vector, u(t), that is fixed in the coordinate frame of the molecule and rotates about the center of mass of the methane molecule with the hydrogen atoms. The decorrelation of u(t) is then defined by (u(O).u(t)). To relate (u(O).u(t)) to a rotational diffusion coefficient, a Fickian phenomenology is postulated for the diffusion of the orientational degrees of freedom, Q, of a spherical top D,V&
Methone/Silicalite
0.5
16 molecules/unit cell; 2 0 0 K
=
ac/at
(18)
where C is the van Hove self-correlation function for the orientational degrees of freedom and D, is the rotational diffusion coefficient.32 The solution to eq 18 can be represented as a series and expansion involving the product of rotation matrices, dmtm, an exponential time decay, F,(t):
F,(t) = bm,me-K/+I)&‘ (20) Sears has shown that at long times, where the phenomenology implied by eq 18 is valid, F , ( t ) and F 2 ( t ) can be equated with Legendre polynomials that express the decorrelaton of the unit vector u(t) defined above 0.0 L 0.0
0.5
1 .o
1.5
2.0
1 2.5
Frequency ( p s - ’ )
Figure 6. Fourier iransform of a series of velocity autocorrelation functions. In the low-occupancyspectrum (a), the low-frequency peak, at approximately 0.4 ps-I, is attributable to a ”shuttling”motion of the molecules within a channel segment between channel intersections. This behavior is less pronounced in (b) at 8 molecules/unit cell and essentially absent in (c) at a loading of 16 molecules/unit cell, where sorbate-sorbate collisions dominate the sorbate dynamics. All spectra show similar behavior at high frequencies.
molecules decreases, and increasingly frequent collisions reduce the likelihood of uninhibited “shuttling” motion; hence, the intensity of the low-frequency peak drops. Since there is no channel system that can allow direct motion in the Z direction, the ~ J Y ) spectrum has no low-frequency component; motion in the Z direction is composed mainly of “rattling” within the pores. The suppression of the low-frequency peak is clearly visible in Figure 6b,c. The great similarity between the individual components of
Pl(d = (u(O).u(t))
(21)
P2(0 = Y2P((u(0)*u(t))2)- 11
(22)
where P l ( t ) and P 2 ( t ) are the first and second Legendre polynomials. As implied by eqs 20, 21, and 22, fitting an exponential expression to the values of P l ( t ) and P 2 ( t )obtained at long times (0.2-1.0 ps) from MD simulations allows an estimate of the rotational diffusivity to be made, If the van Hove self-correlation function, G, obeys the simple phenomenology represented by eq 18, the ratio of the time constants for the decay of the first and second Legendre polynomials should be equal to 3. Rotational diffusivities extracted from the long-time behavior of P l ( t ) in our simulation are shown in Figure 7 at 400 K and a series of loadings. Error bars, in this and all subsequent plots, denote 95% confidence limits. The increase in rotational diffusivity with occupancy results from the increasing rate of sorbate-sorbate collisions, which act to decorrelate the molecular orientation. The (32) Sears, V. F. Can. J . Phys. 1965, 45, 231-254.
8238 The Journal of Physical Chemistry, Vol. 94, No. 21, 1990
June et al.
Energy of Sorption - Methone/Silicolite
-14
-15
-18
by I
0
4
o-oT
.-AT o-oT
8
12
= 200 K = 300 K = 400 K
16
20
molecules/unit cell
Figure 8. Occupancy and temperature dependence of the integral sorption energy of methane in silicalite.
rotational diffusivities predicted here compare well with the values of Gordon,33determined by infrared and Raman spectroscopy of liquid methane at 98 K, and with the value of Thorel et al., determined by quasi-elastic neutron scattering of methane adsorbed on graphite at 55 K.34 However, our values of 0,are much larger than the value of 3 X 1OIo s-l obtained by Jobic et al.35in neutron scattering experiments of methane sorbed in ZSM-5 at 200 K. Furthermore, for the system studied here, the ratio of the time constants in eq 20 averages 1.47;this is in excellent agreement with the experimental ratio of 1.46, measured in liquid methane at 98 K.32*33Therefore, the long-time decay of the orientational van Hove self-correlation function is seen not to be truly diffusive in nature. In addition to the dynamic properties discussed, the MD simulations are also useful in elucidating structural and thermodynamic properties of the intracrystalline sorbate phase. Two equilibrium quantities were calculated: the mean potential energy of the sorbate molecules, ( U ) ,and the single-particle density distribution function, pl(r). The mean potential energy (Or) is an integral energy of sorption, due to sorbate-zeolite and sorbate-sorbate interactions. In a thermodynamic formulation, ( U ) corresponds to the total internal energy of the solid phase (one mole of sorbate plus the zeolite containing it) minus the internal energy of the pure zeolite, minus the internal energy that the one mole of sorbate would possess if it were an ideal gas at the temperature of the system. A related differential property is the defined as the change in partial isosteric heat of sorption,I3 Q,,, molar enthalpy of the sorbate upon desorption. The relation between ( U ) and Q,,can be written as
The integrations in eq 23 are performed at constant temperature along the equilibrium isotherm. The quantity p represents the intracrystalline occupancy, or density, along an equilibrium isotherm. Q,,is the isosteric heat at temperature T and occupancy p’. @ is the molar enthalpy of the sorbate in a gas phase at equilibrium with the solid phase at temperature Tand occupancy p’, and HIG is the ideal gas molar enthalpy of the sorbate at temperature T . At low occupancies, in the Henry’s law region, eq 23 reduces toI3 lim ( U ) = - lim Q,,+ R T P’o
(24)
P’o
Values of the mean potential energy are shown in Figures 8 and 9 as a function of temperature and occupancy for methane and (33) Gordon, R. G. J . Chem. Phys. 1965, 43, 1307-1312. (34) Thorel, P. H.; Coulomb, J . P.; Bienfait. M. Surf. Sci. 1982, 114. L43-L47. (35) Jobic, H.; BEe, M.; Kearley, G . J. Zeolites 1989, 9, 312-317.
Figure 9. Occupancy and temperature dependence of the integral sorption energy of xenon in silicalite.
xenon, respectively. Favorable sorbate-sorbate interactions are evidenced by a slight increase in -( U ) as a function of sorbate loading. The value of 17.5 kJ/mol for the isosteric heat computed from molecular dynamics in the infinite-dilution limit via eq 24 for methane compares well against the value of 18 kJ/mol determined in our previous work with configurational integral techniquesI3 and is in good agreement with the experimentally determined value of 20 k J / m ~ l . ~ ~ The distribution of the sorbate molecules within the pores of the zeolite can be represented by the single particle density distribution function, ~ ‘ ( r ) ~ ’
where the integrals are taken over all external degrees of freedom (positional and orientational) of the sorbate molecules. The single particle density distribution function was computed from the simulated trajectories by spatially binning the sorbate molecules into a three-dimensional histogramI3 and taking advantage of all translational and crystallographic symmetries available within the simulation cell:
In eq 26, po is the overall density of sorbate molecules within the volume of the simulation cell, No is the total number of system snapshots accumulated in the distribution, Nuallis the number (=8) is the number of of unit cells in the simulation box, Nasym asymmetric units contained within the zeolite unit cell, and N (X,Y,Z,to)is the number of molecules contained within a volume AXAYAZ, centered at ( X , Y , Z ) ,at time to. The volume of each of the bins in the histograms, AXAYAZ, was held constant at 0.0079 A3.For methane No was 1000 snapshots, and for xenon No was 5000. The ratio pl(r)/po is a measure of the intensity of sorption at a position r; it indicates how much more concentrated the sorbate molecules are in a given region of the pore system than in a homogeneous fluid containing as many molecules per unit volume as the zeolite crystal. For the simulations performed here, this ratio ranged between 5 and about 70. Three-dimensional contours of constant pl(r)/po are shown in Figures 10-14. These surfaces were generated by using the subroutine 1SOSRF3* and are represented as a net comprised of bold lines. To provide a visual reference, the zeolite pore structure is represented as a net of thinner lines by contouring the xenon(36) Chiang, A. S.; Dixon, A. G.; Ma, Y. H.Chem. Eng. Sci. 1984, 39, 145 1-1459, I46 1-1468, (37) Hansen, J. P.; McDonald, 1. R. Theory of Simple Liquids; Academic Press: New York, 1986. (38) Wright, T. J. J . Comput. Graphics 1979, 13, 182-189; subroutine ISOSRF of the SCD graphics system of the National Center for Atmospheric Research, Boulder, CO..
Methane and Xenon in Silicalite SilicaliteiXenon 0 molecules/U.C., 400 K
Figure 10. Normalized three-dimensional sorbate single particle density distribution function, pl(r)/po, for xenon in silicalite at a temperature of 400 K and an occupancy of 0 molecules/unit cell. The distribution function is contoured at a level of pl(r)/po equal to 20% of its maximum value for the simulation. The sorbate molecules distribute themselves over broad regions of the channel segments but stay away from channel intersections.
The Journal of Physical Chemistry, Vol. 94, No. 21, I990 8239 S ~ I ~ ~ ~ i l i l c i X cI n h omno l c c u l c ~ i C l , 400
K
Figure 12. Contour plot of the normalized sorbate single particle density distribution function at an occupancy of 16 molecules/unit cell and a temperature of 400 K. The contours are again drawn at a level of pl(r)/p,, equal to 20% of its maximum value. Significant localization of the sorbate molecules is evidenced by the fact that the broad regions contoured in Figure IO have broken up into discrete domains of roughly spherical shape. Individual contour regions correspond to regions of highest probability for finding sorbate molecules.
SilicaIi~ciXcnonX molecule\/b.C . 400 K SilicaliteiXenon 0 molecules/U.C., 400 K
Figure 11. Contoured single particle density distribution function, p ' (r)/pO,for xenon in silicalite at a temperature of 400 K and an occupancy of 8 molecules/unit cell. The distribution function is contoured at a level of pl(r)/po equal to 20% of its maximum value for the simulation. The
sorbate structure is very similar to that of the infinite-dilutionstate point.
Figure 13. Single particle density distribution data under the same conditions as Figure 10, contoured at a level of pl(r)/po y u a l to 50%of its maximum value. The core of the sorbate distribution is revealed. Molecules are still seen to be delocalized in the interior of the channels.
zeolite potential at a fixed value. Only half of the unit cell is presented to enhance the clarity of the results. A mirror plane of symmetry, running through the X Z plane and located in the sinusoidal channels, is utilized to bisect the unit cell. As the results for both the methane and xenon systems are very similar and significantly more data were available from the xenon simulations, only results for xenon simulations are presented here. Figures 10, 1 I , and 12 compare the single particle density distribution function for xenon in silicalite a t 400 K at loadings of 0, 8, and 16 molecules/unit cell. The single particle density distribution function is contoured at 20% of the maximum value of the ratio p'(r)/p,,; the corresponding contoured volumes contain about 50% of the sampled centers of mass. At infinite dilution, the sorbate molecules sample a large fraction of the sites within the relatively energetically homogeneous pore structure of silicalite.
Figure 1 1 is drawn at an intermediate occupancy of 8 molecules/unit cell; the distribution appears similar to that of Figure IO. At 16 molecules/unit cell, however, the development of some structure in the sorbate phase is evident; sorbate molecules confine themselves to smaller regions of the channel system (Figure 12). This confinement stems from a compromise between sorbatesorbate interactions and a competition for energy minima within the zeolite. At all concentrations, the sorbate molecules avoid the energetically less favorable channel intersections-a result that is consistent with our infinite-dilution simulation^.'^ Figures 13 and 14 display the core of the single particle density distribution, at a level of p l ( r ) / p , , equal to 50% of its maximum value; the contour volumes contain about 25% of the sampled centers of mass. Figure 13 corresponds to 0 molecules/unit cell at 400 K, and Figure 14 corresponds to 16 molecules/unit cell at 400 K. Clearly,
8240 The Journal of Physical Chemistry, Vol. 94, No. 21, 1990
Figure 14. Single particle density distribution data under the same conditions as in Figure 12, contoured at a level of pl(r)/po equal to 50% of its maximum value. The core of the distribution reveals well-resolved adsorption sites in the sinusoidal channels. Molecules in the straight channel are less densely packed; therefore, they do not appear when contoured at this level.
the sorbate molecules prefer to reside in the sinusoidal channels. At the lower temperature of 200 K and infinite dilution, the contours disappear from the interior of the straight channels at a sorption intensity level, pl(r)/po, equal to 50% of its maximum value. In general, p ' ( r ) / p o becomes less diffuse at lower temperatures but still retains the qualitative features discussed here for the higher temperatures.
Conclusions The results presented here indicate that M D simulations can successfully describe the dynamics of small molecules sorbed in
June et al. zeolites. The assumption of a rigid lattice appears to be reasonable when the sorbate molecules are small with respect to the pores of the zeolite and interact weakly with the atoms of the lattice. Various properties, such as the sorbate self-diffusivities, isosteric heats, and the sorbate local density distributions, were computed over a wide range of occupancies and temperatures. Where experimental evidence is available, e.g., for the self-diffusivities, the comparison with experiment is very favorable. For some quantities for which no experimental evidence is available, such as the sorbate local density distribution, results generated here are consistent with results from our previous configurational integral work in the infinite-dilution limit. Different factors were found to influence the properties of the sorbate molecules at different levels of occupancy and temperature. At low occupancy, energy barriers presented by the channel intersections act to contain the sorbate molecules within the channel segments and induce a low-frequency "shuttling" motion. As occupancy rises, collisions between sorbate molecules decrease the mean free path within the zeolite and sorbate "rattling" is the only dominant mode of motion. Sorbate-sorbate collisions act to effectively lower energy barriers to sorbate motion, as is evidenced by the decrease in activation energy with increasing occupancy. increasing the occupancy of sorbate molecules also makes the density distribution less diffuse. As temperature is lowered, sorbate molecules are found to preferentially populate the sinusoidal channels over the straight channels, because the former are energetically more favorable. Furthermore, small nonpolar molecules such as methane and xenon do not frequent the channel intersections at any of the temperatures investigated here.
Acknowledgment. Support for this work was provided by the National Science Foundation through Grant NSF-CBT-87 18027 and by W. R. Grace and Co. Conn. D.N.T. further acknowledges the National Science Foundation for a Presidential Young Investigator Award, No. DMR-8857659. Finally, all of the authors thank the San Diego Supercomputer Center for a generous allocation of computer time on the Cray X-MP/48 and Professors H. Pfeifer and J. Karger for providing the experimental values of the self-diffusion coefficient of xenon in silicalite. Registry No. Xe, 7440-63-3; CH,, 74-82-8.
