J. Phys. Chem. 1994, 98, 11948-11961
11948
Investigation of the Dynamics of Benzene in Silicalite Using Transition-State Theory Randall Q. Snurr, Alexis T. Bell, and Doros N. Theodorou* Department of Chemical Engineering, University of Califomia, and Center for Advanced Materials, Lawrence Berkeley Laboratory, Berkeley, Califomia 94720 Received: May 18, 1994@
The dynamics of benzene in silicalite at low loading was investigated using transition-state theory. Benzene was found to diffuse by infrequent hops between preferred adsorption sites. Potential energy minima and saddle points were located using an atomistic model, and diffusion paths connecting pairs of minima were constructed through each saddle point (transition state). The intrinsic reaction coordinate (IRC) approach was used to construct the diffusion paths in six dimensions. The IRC equations are presented for the motion of a rigid body (benzene) through a static potential field (silicalite). A rate constant for each transition between minima was calculated using a harmonic approximation to the potential energy function. From the rate constants, the self-diffusivity was computed with a dynamic Monte Carlo simulation. An activation energy of 36.7 kJ/mol was calculated. This is larger than the experimental value, and the predicted diffusivities are 1-2 orders of magnitude smaller than experiment. Likely reasons for this discrepancy are the harmonic approximation invoked in calculating the rate constants and our neglect of zeolite flexibility in the calculations. The predicted time scales for local motions within the channel intersections agree well with spectroscopic resulk. Many of these motions correspond to rotations of the benzene molecule about its c6 axis.
Introduction Medium-pore zeolites such as ZSM-5 are widely used in the petroleum and petrochemical industries as catalysts for xylene isomerization, catalytic dewaxing, and conversion of methanol to gasoline.' Z S M J is well suited for these applications because its pore size inhibits coke formation while still admitting larger molecules such as monocyclic aromatics. Additionally,the pore size of ZSM-5 enables desirable shape selectivity. For example, xylene isomerization catalyzed by ZSM-5 is biased toward production of the most valuable isomer, p-xylene, because the diffusivity of p-xylene in ZSM-5 is much greater than that of m-xylene of o-xylene, thus allowing the para product to diffuse selectively out of the zeolite crystal.2 There is much interest in understanding the transport mechanisms of aromatics in ZSM-5 and other zeolites to spur development of new shape-selective processes and to allow further selectivity improvements in existing processes. A large body of experimental work has provided information on the diffusion of aromatics in ZSM-5 and its all-silica form known as silicalite. The experiments range from gravimetric measurements of diffusion coefficients and activation energies3 to spectroscopic measurements of the time scales for rotation of an adsorbed benzene molecule about its c6 axis4 Interpretation of spectroscopic studies is, however, not straightforward, and the integration of the microdynamics information with the macroscopic transport picture is not always easy. Molecular simulation based on an atomistic description of zeoliteisorbate interactions is a powerful tool for gaining insight into the transport mechanisms of adsorbed molecules and for predicting the diffusivity and activation energy for diffusion. Nowak et al.5 have used an atomistic model of benzene in silicalite to address the energetics of diffusion in a simple way by calculating potential energy maps. A pairwise-additive Lennard-Jones potential between atomic sites was employed. The maps show the variation in potential energy as a benzene
* To whom correspondence should be addressed. @
Abstract published in Advance ACS Abstracts, October 15, 1994.
0022-3654/94/2098-11948$04.50/0
molecule moves along the center of a silicalite channel. At each position, the change in potential energy was computed as a function of the rotation about a single axis. The potential maps allow an estimate of the activation energy for diffusion to be made, as well as the identification of possible sorption sites. An improvement on this technique was made by Pickett and co-workers.6 They traced out minimum-energy paths for diffusion by shifting a benzene molecule along the channel axis in steps of 0.1 A; at each step, they minimized the energy with respect to movement along the remaining five translational and rotational degrees of freedom. The activation energy for diffusion can be estimated from the minimum-energy paths. Talu7 has used such paths to explore the reorientations that a benzene molecule undergoes as it passes through the straight and sinusoidal channels of the silicalite pore system. VignkMaeder and Jobic8 calculated minimum-energy paths for diffusion of benzene in silicalite, using a different potential model that included electrostatic and polarization terms in addition to dispersive and repulsive interactions. They also calculated the barriers for rotation about several axes for molecules located in the potential energy minima. Schroder and Sauer9calculated minimum-energy paths for diffusion of single benzene molecules through silicalite and then applied transition-state theory (TST) to obtain the diffusivity. Diffusion through the sinusoidal channels was neglected, and it was assumed that a single minimum in the intersection was separated from a minimum in the straight channel by a single barrier. Librational and translational frequencies of the benzene molecule were calculated at the minima and the transition state to obtain the activation entropies and enthalpies, and then transition-state theory was used to calculate the diffusivity. Schroder and Sauer point out that many minimum-energy paths exist through the straight and sinusoidal channels. Also, it is rather difficult to follow exactly one of the paths; depending on the algorithm used for stepping along the path, the molecule may experience sudden unphysical reorientations at points along the predicted path. That is, the molecule suddenly switches from one (smooth) path to another. Starting from a given configuration,
0 1994 American Chemical Society
Dynamics of Benzene in Silicalite the path generated by the method of Pickett et al. is not guaranteed to go to the lowest energy minimum or to go through a low barrier in traversing between minima. Molecular dynamics (MD) simulation is another technique often used to obtain information from an atomistic model about the dynamics of species adsorbed in zeolites.' The classical Newtonian equations of motion are integrated forward in time, and the positions of the molecules are tracked. From the meansquared displacement, the self-diffusivity is calculated using the Einstein equation.lo MD studies of diffusion in silicalite have been performed by a number of investigators for simple sorbates such as methane."-14 More recently, the diffusion of larger molecules, such as butane and hexane, has been studied using MD simulation^.^^,^^ The simulations also provide predictions of the heats of sorption, siting locations of molecules within the channel system, and details of the sorbate motion. The diffusivities and other properties predicted from the MD simulations to data have generally been in good agreement with experimental measurements. Unfortunately, even with today's fastest supercomputers, MD simulations cannot probe time scales longer than nanoseconds. From the experimental diffusivity of benzene in silicalite, one can estimate that the time required for a benzene molecule to diffuse several unit cells is orders of magnitude longer than time scales accessible to MD. From a variety of experimental evidence,I6-l9 it is thought that there are energetically favored sorption sites for benzene in silicalite. Recent Monte Carlo simulation^^^-^^ reveal cleanly-defined sorption sites in the straight channels, the sinusoidal channels, and the channel intersections. According to the simulations, virtually all adsorbed benzene molecules reside in these sites; that is, the probability of finding a molecule between the sites is extremely low. It would seem reasonable then to view the transport of benzene in silicalite as proceeding by infrequent hops of molecules from site to site. Such a model is consistent with deuterium NMR studies of Zibrowius et aLZ4 In this paper, transition-state theory (TST) is systematically applied to the low occupancy diffusion of benzene in silicalite. TST is particularly appropriate for the estimation of rates of infrequent dynamical events, such as the hopping of a molecule from one state to another.25 TST provides a framework for calculating the rate constants for the hopping events, and from the rate constants the diffusivity can be obtained. Because it can focus attention on the events that actually contribute to translational motion of the molecules (that is, the hopping motions), TST is a much more computationally efficient technique for studying infrequent event processes than MD. A molecular dynamics simulation of benzene in silicalite would "waste" all of its time tracking the local movements of the benzene molecules within the sorption sites, but it would be unable to sample the long-time motion between states. Like MD, the TST approach described in this paper is a predictive method that starts with a detailed atomistic model of the system and attempts to arrive at the macroscopically-observed selfdiffusion coefficient. In addition, a wealth of information about the mechanism and rates of sorbate motion inside the zeolite is gleaned over a very wide spectrum of time scales. Transition-state theory has been applied by several authors to the motion of small molecules between the a cages of zeolite A.26-29 The transition state was taken to be in the windows that connect the a cages. Kiirger et aLZ9calculated a diffusivity for methane in 5A zeolite that is in good agreement with the experimental value from pulsed field gradient Nh4R measurements. Properly accounting for the rotational degrees of freedom of methane was found to be important in obtaining the correct diffusivity value. More recently, June et aL30 used transition-state theory to investigate the low-temperature dif-
J. Phys. Chem., Vol. 98, No. 46, 1994 11949 TABLE 1: Potential Parameters for BenzendSilicalite Interaction@
co HO C H 0 Si
~
1.335 x lo6 1.599 x lo5
1.807 x lo3 5.108 x lo2 -0.15
+0.15 -1 +2
Partial charges are given in units of the electron charge le/.
fusion of the spherical sorbates xenon and sulfur hexafluoride in silicalite. For the more complicated silicalite structure, the locations of the transition states are not as intuitively obvious as in type A zeolites, but algorithms exist for locating saddle points on the zeolite/sorbate potential energy hypersurface. In the work of June et al., these saddle points were taken to be the transition states. Minimum-energy diffusion trajectories connecting sorption sites were found by constructing two steepest descent paths from each transition state. The paths terminate at potential energy minima, and sorption states were constructed around the minima. Rate constants for transitions between sorption states were determined from the evaluation of the configuration integrals in the states and on the dividing surfaces between states. Dynamical correction factors to the rate constants, accounting for recrossings of the dividing surfaces and correlated multistate jumps, were computed through shorttime MD simulations initiated on the dividing surfaces. The rate constants were then used in a continuous-time/discrete-space Monte Carlo calculation to determine the self-diffusivity of the sorbate.
Model The atomistic representation and the potential parameters for the zeolite/sorbate interactions are exactly the same as those used in our calculations of adsorption thermodynamics of benzene in silicalite.20-22 The orthorhombic form of silicalite having space group Pnma is considered here, with the atomic coordinates taken from the X-ray diffraction work of Olson et aL31 A number of crystal phase transitions are known to occur in silicalite systems. (See refs 18, 21, and 22 and references therein.) For low loadings of benzene we expect to find silicate in the Pnma form. Silicalite has two sets of intersecting channels, each about 5.5 8, in diameter. A set of straight channels run in the y direction and are intersected by a set of sinusoidal channels that run in the x direction. The silicalite structure is considered here to be rigid and to be comprised purely of silicon and oxygen with no aluminum. Benzene is modeled as a rigid, planar molecule with explicit carbon and hydrogen atoms. The benzene and silicalite are assumed to interact through a pairwise-additive potential between atoms of the adsorbate and atoms of the zeolite framework. The atom-atom interactions are modeled with a Lennard-Jones plus point-charge potential:
where i and j are atoms of the sorbate molecule and of the silicalite structure, respectively, and rv is the distance between them. Au and Bu are the Lennard-Jones constants, and qi and q, are the partial charges of the atoms. The values of Aq, Bu, qi, and % are given in Table 1. The basis for choosing these values has been presented earlier.20,21A cutoff radius of 13 A is applied to the Lennard-Jones interactions, and the long-range electrostatic interactions are calculated using the Ewald summation technique.I0 The benzene/silicalite model described here
11950 J. Phys. Chem., Vol. 98, No. 46, 1994 has been used in earlier simulationsof low-occupancy adsorption thermodynamics.20 The model yields Henry’s constants and isosteric heats of adsorption that are in excellent quantitative agreement with experimental data.
Methodology and Calculations Overview. This section describes the application of transition-state theory to the diffusion of a rigid body in an external potential field, specifically to the diffusion of a rigid, planar benzene molecule in the potential field of a static silicalite pore system. A rigid benzene molecule has six degrees of freedom. These can be taken as the center-of-mass coordinates of the molecule, ro, describing its position, and three Eulerian angles, Y, describing its ~ r i e n t a t i o n . We ~ ~ will define as a state or microstate the region in the 6-dimensional configuration space around a local minimum of the potential energy qro, Y). We will define as a macrostate a collection of neighboring microstates separated by energy barriers that are low relative to kT. The molecules move rapidly among the microstates within a particular macrostate but only infrequently cross the highenergy barriers separating the macrostates. We assume that when a molecule enters a new macrostate it quickly thermalizes there, and a further jump of the molecule into another macrostate is uncorrelated with the previous jump. The evolution of the system in time is viewed as a Poisson process consisting of successiveuncorrelated jumps between neighboring macrostates. An ensemble of such systems can thus be described by the following master equations:
p ~ ( tis) the probability of finding a system of the ensemble in macrostate A at time t, and kA-B is a first-order constant for a jump from A to B . The equilibrium occupation probabilities, piq,must obey the condition of microscopic reversibility eq
(3) Piq = kB-,4 PB The self-diffusivity of the adsorbate can be obtained from the solution of the master equation under prescribed initial conditions. In the rest of this section, the various steps leading up to the calculation of the diffusivity are presented in detail. First, the minima and saddle points of the potential energy hypersurface are found. Next, the intrinsic reaction corrdinate (IRC) approach of F ~ k u isi ~applied ~ to determine all minimum-energy paths connecting pairs of adjacent minima, each path going through a saddle point. This serves two purposes: first, it establishes the connectivity of the microstates; second, the IRC paths provide information on the most likely motions of the adsorbed molecules. Note that, for the case under consideration here, these are paths in six dimensions, tracking both translational and rotational movements of the benzene molecule. Next, rate constants are calculated for all possible transitions (IRC paths) in the system. Then, based on the values of the rate constants, the microstates connected by fast transitions are lumped together into macrostates. Occupation probabilities for the macrostates are calculated, as well as overall effective rate constants between the macrostates. Finally, the self-diffusivity is calculated by solving the master equation with a continuous-time/discretespace Monte Carlo simulation. Note the hierarchical nature of the approach. The detailed atomistic model is used to identify the sorption sites and the rate constants for moving between them, but the diffusivity is obtained from a more coarse-grained lattice Monte Carlo simulation. The rate constants serve as the connecting link between the two levels of the hierarchy, as they kA-B
Snurr et al. are obtained as outputs of the atomistic calculations and then used as inputs to the lattice Monte Carlo simulation. Identification of Minima and Saddle Points. Minima of the potential energy hypersurface were found using a quasiNewton minimization algorithm employing the BFGS Hessian matrix updating scheme. Due to the symmetry of the silicalite structure, only the asymmetric unit of the silicalite unit cell had to be searched for minima. Degenerate configurations were only counted once. (Rotating benzene by 60” about its axis gives rise to a degeneracy.) An exhaustive search for all minima was conducted as follows. Initital configurations were generated randomly. Those with energies below a cutoff of 0 kJ/mol were used as starting points for minimization; higher-energy initial configurations were discarded. Minimizations from 10 000 starting configurations yielded 27 unique minima. Minimizations from an additional 50 000 starting configurations yielded no additional minima. The search for fist-order saddle points of the potential energy used a Newton-like algorithm due to Baker.34 Recall that at a first-order saddle point the gradient of the potential energy is zero, and the Hessian matrix of second derivatives has (exactly) one negative eigenvalue. An exhaustive search for all transition states in the asymmetric unit was conducted in a manner similar to the search for minima. The algorithm requires the Hessian matrix; this was calculated by finite difference from the analytic first derivatives. The Baker algorithm can also be used to find minima. As a check, we performed the search for minima a second time using the Baker algorithm instead of the BFGS method. The same 27 minima were obtained. Diffusion Paths. The intrinsic reaction coordinate (IRC) is a classical path in configuration space that connects a saddle point with two local minima.33 To obtain the classical equations of motion for a rigid body, we start with the kinetic energy expres~ion:~~
+ +
+
9-= V2M(X2+ Y2 Z2) 1/2(zlw12 z2w;
+ Z3w?)
(4)
where M is the mass of the body; X , Y, and Z are the centerof-mass coordinates, with time derivatives X , Y, and Z; ZI,12, and 1 3 are the principal moments of inertia of the body, defined in a set of axes fixed with respect to the rigid body; and w1, w2, and w3 are the angular velocities about the body-fixed principal axes. The Eulerian angles specify how the body-fixed axes are oriented with respect to the laboratory (or zeolite) coordinate system. We will use the convention employed by Goldstein for the Eulerian angles 4, q, and The angular velocities are related to the time derivatives of the Eulerian angles as follows:
4sin e sin q + 8 cos q w2 = 4sin e cos - 6 sin q w1 =
w g = &OS
e+
(5)
Substituting (5) into (4),we obtain
+ + Z2>+ 1/z~l(4’ sin2 e sin’ q +
T= 1 / 2 ~ ( X 2Y2
8’ cos’ q + 248 sin e sin q cos q) + 1/2~2(42 sin2 e cos2 q 8’ sin2 q - 248 sin e cos q x sin q) + 1/2Z3(42 cos’ 8 + +2 + 244 cos 0) ( 6 )
+
The equations of motion are obtained from Lagrange’s equations: (7)
where y i s the Lagrangian, equal to T-
R and qj represents
J. Phys. Chem., Vol. 98, No. 46, 1994 11951
Dynamics of Benzene in Silicalite each of the elements of the generalized coordinate vector (X, Y, Z, q5, q ,0). For example, for qj = X,
Following Fukui’s approach, we obtain the IRC by considering motion from a given point with an infinitesimal velocity. Integrating forward in time for a very small time step, we have for qj = X
(9) The initial condition is XO = 0 at to = 0 because we want the path of infinitesimal velocity or a path with continuous dissipation of the kinetic energy. Being a path with infinitesimal velocity, the IRC can be regarded as a “center line” about which actual motions will occur. Applying this initial condition to eq 9 yields
or
a .ti ax
MdX=--tdt
(11)
Applying this approach to all six degrees of freedom produces the following:
a 7’dt ax a .t) - -t dt
with respect to the mass-weighted Cartesian c o o r d i r ~ a t e s ) . ~ ~ ~ ~ ~ The rest of the path would simply be a steepest-descent trajectory in the mass-weighted Cartesian coordinates. For the case of rigid-body dynamics, where we are working in the generalized coordinates q, we can obtain a similar prescription for the first step. Near the saddle point, the gradient of Y’can be expressed as
where H,, is the Hessian matrix of second derivatives with respect to the generalized coordinates and dq is the distance vector from the saddle point. Substituting eq 14 into eq 13 yields -dq 1 = G-IH,, d q (15) dz Thus, the first step dq is along the eigenvector corresponding to the negative eigenvalue of the matrix G-IH,,. Equation 13 provides the equations for the diffusion paths. Some details of the computation will now be discussed. The steps along the diffusion path were taken with constant step size of ds = 0.001. With distances measured in angstroms, angles measured in radians, and masses measured in g/mol, the following equation is used to calculate a length in the 6-dimensional configuration space:32
M dX = - -t MdY=
ay
a az
.ti
M d Z = - -tdt
+
+
(IIsin2 0 sin2 q I2 sin2 0 cos2 II, I3 cos2 0) dq5 (II sin 0 sin q cos I) - Z, sin 8 cos sin )I
(IIcos2 11,
+ Z2 sin2 qj d 0 + (I,sin 0 sin 11, cos
+ +
-
$J
a .ti dt ao
Z2 sin 8 sin q cos q)dq5 = - -t
(12)
Equations 12 provide a differential description for a line in sixdimensional space, with parameter t = -t2/2. We can use them to step along this line (the IRC) as follows. Set - t dt equal to some small scaling factor d t ; calculate the gradient, V,U’= (a Ylax, a Yiay, a Ylaz, a Yia& a Ylaq, a Yia0jT; calculate dq = (dx,dY, dZ,d4, d q , d0jTfrom eqs 12, and step along the IRC. Note that the equations for dq5, d v , and df3 are coupled. Equations 12 can be written more compactly as35 d q = G-IV, Y’dt
(13)
where G is the metric tensor of the configuration space32but also includes the appropriate mass and moment of inertia factors. Each IRC is initiated at a saddle point and descends in two directions. The path in each direction terminates at a minimum. The prescription just given, however, does not tell us how to take the first step. At the saddle point, the gradient is zero, and eqs 12 cannot be used to calculate dq. If one were working in mass-weighted Cartesian coordinates, the first step would be taken in the direction of the eigenvector corresponding to the negative eigenvalue of the Hessian matrix (derivatives taken
where again Gg is an element of the metric tensor of the space. Steps of constant size were assured by the choice of dz at each step. With a step size of ds = 0.001, the values of V, Y’at the beginning and end of each step were close enough that the sequence of rectilinear dqs constituted an excellent approximation to the continuous diffusion path. Visualizations of all resulting paths in three dimensions were smooth. It is verified that each step along the path resulted in a lowering of the potential energy and that each path ultimately terminated at one of the already-known minima. Diffusion paths were computed in the asymmetric unit of the silicalite structure, and the paths for the full unit cell were obtained by applying the appropriate symmetry operators. It is well known that whenever the Eulerian angle 0 approaches 0 or n,the rigid-body equations of motion di~erge.lO3~~ We got around this problem in a simple way. Each path was attempted in the Eulerian angles defined as mentioned above. If the value of 0 got too close to 0 or n, the path was aborted and rerun using a different convention for the Eulerian angles (the so-called xyz c ~ n v e n t i o n ) . ~ It ~is conceivable that a particular path could fail with both conventions, but fortunately no such paths were encountered. Calculation of Rate Constants. Vineyard37 has presented an expression for the TST rate constant in many degrees of freedom. In his formulation, the degrees of freedom used are the mass-weighted Cartesian coordinates of the atoms in the system, x. The rate constant for going from state i to state j is25
Q m is the partition function for the canonical ensemble
where the sum over r represents integrating over all coordinates and momenta of the system. 53% the Hamiltonian, 5T= T+ Y! The partition function in the numerator of eq 17 is evaluated
11952 J. Phys. Chem., Vol. 98, No. 46, 1994 at the dividing surface in configurational space between states i andj, and the partition function in the denominator is evaluated over state i. Note that the dividing surface between i and j contains the saddle point or transition state. After integrating the kinetic energy, k z can be written as the ratio of two configurational integrals:
To evaluate the configurational integrals in an approximate way, Vineyard expands the potential energy at the minima and at the transition states in a Taylor series to second order. In this harmonic approximation, the rate constant becomes N
m= 1
where N is the number of degrees of freedom in the system (3 times the number of atoms), Yh is the energy at the saddle point, L I ‘ I , is the energy at the minimum i, v k are the normal frequencies for vibration of the entire system at the starting point of the transition, and v i are the normal frequencies of the system constrained in the saddle point configuration. The frequencies v k are obtained as the square roots of the N eigenvalues (Am) of the Hessian matrix H,, evaluated at the minimum i, where the derivatives of H, are with respect to the mass-weighted Cartesian coordinates x. When evaluated at the saddle point, the Hessian H,, has N - 1 positive eigenvalues and one negative eigenvalue. The square roots of the positive eigenvalues are the frequencies vi. Note that the integral in the denominator of eq 19 is proportional to the equilibrium probability of occupying microstate i, pyq. We now seek to develop similar expressions for our model of benzene in silicalite. Vineyard analyzes a system in which all atoms move according to prescribed potentials, and the degrees of freedom are then taken as the mass-weighted Cartesian coordinates of the atoms, x. We wish to analyze the diffusion of a rigid benzene molecule within a rigid silicalite framework. Since the zeolite is taken as rigid, the zeolite will only enter the analysis as an external potential field in which the benzene molecule moves. Treating the benzene molecule as a rigid body allows us to use only six generalized degrees of freedom (N = 6) instead of the 36 atomic Cartesian coordinates. In Vineyard’s derivation, the kinetic energy part of the Hamiltonian can be integrated analytically because the kinetic energy can be written in the form
c N
T=
Snurr et al.
(22)
Going through a derivation analogous to Vineyard’s, we obtain again eq 20, where now the frequencies are calculated from the eigenvalues of H,,; that is, the derivatives are taken with respect to the mass-weighted generalized coordinates r. Working with r is inconvenient, but the eigenvalues of H,, are the same as those of G-lH,, (see Appendix), which provides a way for estimating the configurational integrals in the harmonic approximation and thus obtaining the rate constants. As in Vineyard’s work, we also calculate the microstate occupation probabilities within the harmonic approximation. Definition of Macrostates. After the rate constants are calculated as discussed above the transitions of benzene in silicalite between all identified microstates, one observes that the rate constants span many orders of magnitude. Not all of the time scales associated with such disparate rate constants are relevant for diffusion. Many of the individual transitions occur very rapidly with low energy barriers and are simply rotations of the molecule or small local displacements. We are more interested in the less frequent motions of the molecules through the high-energy bottlenecks in configuration space; it is these motions that determine the long-range diffusive mobility. To focus on these rate-limiting transitions, we wish to collect together groups of microstates that are separated by only lowenergy barriers (and thus mutually accessible by “fast” motions). We will call each of these collections a macrostate. We wish to calculate the occupation probabilities of the macrostates and the effective rate constants for transitions between the macrostates. These quantities will be calculated from the information already obtained for the microstates and the transitions between microstates. First we must determine a criterion for deciding which transitions we will call fast and which we will call slow. Looking at the large number of rate constants for the transitions between microstates, it is difficult to find the natural separation of time scales that should be used as the basis for identifying fast transitions. Rather than looking at all of the rate constants, we can find “natural” relaxation times of the system as described next;38 if a clear separation of time scales is seen in these relaxation times, we can use that as the criterion for distinguishing fast from slow motions. To accomplish this, we will construct the matrix k of rate constants. Each off-diagonal element k, will simply be the rate constant k,, associated with a transition from microstate i to microstatej. If more than one path exists from i to j , then we use the sum of the ki-J values for all such paths. To satisfy conservation of probability, the diagonal elements of the rate matrix are defined as
p,2
i=l
where the pi’s are the momenta conjugate to the mass-weighted Cartesian coordinates, x. We can get an expression for kTsT similar to Vineyard’s by starting with the kinetic energy expression of eq 4 rather than eq 6. Our “mass”-weighted generalized coordinates, r, are
r1 = f
ix
r2 = &MY
r3=&z
Given that detailed balance holds, it can be shown that k has one zero eigenvalue and that the other eigenvalues are real and negative.38 The eigenvector corresponding to the zero eigenvalue gives the equilibrium probability of occupying each microstate, p;’. Since we already know plq,this simply serves as a check on the calculation. The nonzero eigenvalues, A,, give the relaxation times of the system: 5, = While a separation of time scales could not be easily seen from k, the values of A show a very clean break of 3 orders of magnitude. At 200, 300, 400, and 500 K, we find that there are three eigenvalues corresponding to slow motions and then a break of
J. Phys. Chem., Vol. 98, No. 46, 1994 11953
Dynamics of Benzene in Silicalite
TABLE 2: Eigenvalues, Ai, of the Matrix k of Rate Constants at 300 K” -1.20 -4.84 -6.15 -7.27 -2.15 -4.33 -1.56 -1.64 -1.65
x 10-4 x 10’
103 x 103 x x x x x
lo6 107 10’ 10’ lo8
-6.47 -5.35 -3.11 -3.78 -1.13 -1.43 -1.80 -6.30 -7.60
x x x x x x x x x
10’ 109 10” 10’’ 10” 10” 10” 10” 10”
-9.10 -1.09 -1.25 -1.49 -1.64 -1.90 -1.93 -2.27 -2.63
of possible events that could occur next is then 1000
x 10” x 10” x 10” x 10”
x 10” x lo1*
x 10” x 10”
x 10’’
Values are given in units of s-l. The first eigenvalue is essentially zero as expected. The next three eigenvalues correspond to slow relaxations. These are separated from the faster motions by a gap of several orders of magnitude in the eigenvalue spectrum.
3-4 orders of magnitude followed by a near continuum of faster motions. The spectrum of eigenvalues at 300 K is given in Table 2. The reader may note the similarity of this approach to that of Wei and Prater in analyzing systems of linear chemical reactions; they diagonalize the matrix of rate constants to determine the relaxation times of the network.39 Microstates are grouped together into macrostates if the rate constants between all pairs of microstates are higher than the eigenvalue on the high side of the observed break in the spectrum of eigenvalues. Note that the grouping could be different at different temperatures. The equilibrium occupany probability of the macrostate is set equal to the sum of the probabilities of the microstates that constitute it.
The microstates within a macrostate are viewed as being in equilibrium with one another, so that when a molecule moves into the macrostate, the microstates instantaneously become populated with probability pgq/piq.40 To reflect this rapid equilibration (which occurs on a time scale much faster than that for jumping to another macrostate), the rate constants of transitions out of the macrostate are weighted by p‘”piq in coarse-graining the network of microstate^.^^ In the dynamic Monte Carlo simulations of diffusion detailed below, the macrostate is represented by a single point in three-dimensional space (within the zeolite) that is taken as the average of the positions of the constituent microstates, with weighting factor p;q/piq. Dynamic Monte Carlo Calculations. To obtain the diffusivity from the rate constants, the master equation, eq 2, is solved numerically with a stochastic s i m ~ l a t i o n . ~Molecules ~ , ~ ~ , ~ ~in the simulation move between adjacent sites with frequencies determined by the rate constants ~ A - B . Since each possible event A B has a first-order rate constant associated with it, the time between events of this type is an exponentially distributed random variable
-
-
The chain of events constitutes a Poisson process. The average time (tAB) between A B moves is calculated to be 1lkA-B from
The simulations were run with 1000 noninteracting (“ghost”) molecules to obtain good statistics. At the beginning of the simulation, the molecules were distributed among the macrostates according to the equilibrium probability distribution (piq}. At a given time during the simulation, each of those 1000 molecules has a number, ni, of transitions available to it depending on what state the molecule is in. The total number
nT=
C nj
(27)
i= 1
Given n~ possible moves, with the time that must elapse for event i to occur being represented by the random variable ti, the time interval for the first event to occur is tmin = min(tl,t2,...,t,J. The time interval for the first event is also an exponentially distributed random variable, and the overall rate constant is equal to the sum of the individual rate constants
The probability that the ith event will occur first is (kA-B)i/@. The algorithm for the continuous-time/discrete-space Monte Carlo proceeds as follows. For each move, the overall transition rate e for the system in its current configuration is calculated from eq 28, and the time interval z to the first event is selected according to the Poisson statistics, eq 25:
z=
ln(1 - 5)
-Q
(29)
where 5 is a random number between 0 and 1. The time t is added to the simulation clock, and the actual event to occur is chosen from the probability distribution (kA-B)j/@. The mean-squared displacements are accumulated as the system evolves in time, and the Einstein equation is used to calculate the diffusivity in each of the x, y , and z directions:
Dxx= lim p- 2t The overall average self-diffusivity at low loading is then given by
Results A total of 27 unique minima and 100 unique transition states were identified for benzene in the asymmetric unit of silicalite. The energies of the minima ranged from -62.1 to -1.4 kJ/ mol, and the transition states ranged from -61.1 to +86.4 kJ/ mol. Figure 1 shows part of the silicalite unit cell. A straight channel is seen running up through the center of the figure, with sinusoidal channels running left to right. The large green spheres represent the minima, and the medium-sized red spheres represent the transition states. The small dots represent points along the diffusion paths, with dark blue signifying paths within a macrostate and white signifying paths between macrostates. Note that the minima, the transition states, and all points along the diffusion paths are actually points in a 6-dimensional configuration space that includes the molecule’s rotational degrees of freedom. Figure 1 just shows a 3-dimensional projection, indicating the positions of the center of mass. Also note that all paths are shown; some of them have very high energies and are seldom traversed. Figure 2 is similar to Figure 1 but focuses on the region near one of the channel intersections. Figures 1 and 2 clearly show that the macrostates correspond to adsorption sites in the straight channels (S), the sinusoidal channels (Z), and the channel intersections (I). One might expect such macrostates based on the geometry of the silicalite pore system. It should be emphasized, however, that we have not imposed this particular grouping of the microstates into macrostates. Rather, this grouping emerged naturally from the
11954 J. Phys. Chem., Vol. 98, No. 46, 1994
Snurr et al.
Figure 1. Diffusion paths followed by the center of mass for benzene in silicalite. The silicalite straight channel can be seen in the center of the figure, with the sinusoidal channels running left to right. The large green spheres represent the minima, the medium-sized red spheres the transition states, and the small dots the diffusion paths. Dark blue paths are within macrostates, and white paths connect minima in different macrostates.
systematic procedure outlined above for defining the macrostates based on the rate constants for transitions between microstates and the natural separation of time scales. Having adsorption sites I, S, and Z is consistent with previous experimental and modeling work in this area. (See references in Snurr et al.22) Evaluating the configurational integrals with the harmonic approximation yields occupation probabilities for the I, S, and Z sites of 0.9989, 0.0008, and 0.0003, respectively, at 300 K. These values can be compared with values of 0.9963, 0.0023, and 0.0014 at 293 K obtained with the same atomistic model from full Monte Carlo integrations.22 As discussed below, the discrepancy provides some measure of the error introduced by the harmonic approximation. The occupation probabilities obtained by Li and tal^^^ (0.873, 0.021, and 0.060 for I, S, and Z sites, respectively) are considerably different; a different potential representation and, conceivably, a different definition of the sites were used in that work. Some of the microstates in
the S and 2 sites have energies comparable to microstates in the I sites, but the intersections are favored entropically, leading to the predicted occupation probabilities. As discussed previously,20-22our prediction that benzene is in the channel intersections at low loading agrees with the experimental evidence in the literature. The molecules spend most of their time in the channel intersections and occasionally move to another intersection by way of either a straight channel or a sinusoidal channel. The rate constants for hops between macrostates are strongly temperature dependent, as they correspond to transitions with high activation barriers. Many of the transitions between microstates withm the same macrostate have very low activation energies, and the associated rate constants are relatively insensitive to temperature. These transitions are just small local displacements and rotations of the benzene molecule. This division between slow transitions with high activation energies
.
Dynamics of Benzene in Silicalite
J. Phys. Chem., VoE. 98, No. 46, 1994 11955
Figure 2. Diffusion paths near the channel intersection region. Macrostates (in dark blue) corresponding to straight channels, sinusoidal channels, and the channel intersection are clearly seen.
and fast transitions that are insensitive to temperature agrees with the picture obtained from deuterium N M R studies= and quasi-elastic neutron scattering (QENS),19 where it was seen that fast rotational motions persist down to temperatures around 100 K. Gjven the macrostate description of the system, the dynamic Monte Carlo simulations yield plots of mean-squared displacement versus time like the one in Figure 3, which is for benzene in silicdite at 300 K. The simulations were run for 10 million moves with 1000 noninteracting molecules. Each Monte Carlo on took about an hour on a workstation. Note that the time step during the simulation depends on the temperature (because the rate constants depend on temperature) and is inversely proportional to the number of molecules used. At 300 K with 1000 molecules, the average time step is 2.8 x loy7 s, while at 500 K with 1000 molecules it is 7.0 x average times between moves of a single molecule are thus 2.8 x and 7.0 x lo-’ s at 300 and 500 K, respectively.
The TST methodology automatically chooses the appropriate average time step for each temperature. The mean-squared displacement in Figure 3 corresponds to a distance of about 20 unit cells. The linearity of the plot indicates that the system is diffusive and that the self-diffusivity can be correctly obtained from eqs 30 and 31. As was seen for other sorbates in silicalite,15v30 benzene diffuses most readily in the y direction, the direction of the straight channels. There is also significant movement in the x direction, corresponding to the sinusoidal channels, but much less in the z direc is no direct channel system than runs in the z direction, but motion in z is possible by a series of moves between the straight and sinusoidal channels. The self-diffusivities calculated from eqs 30 and 31 are reported in Table 3. The nondiagonal components of D fluctuate about zero. Again, the anisotropy of the diffusivity tensor is evident. K&geQ2recently derived a simple correlation between the principal values of the diffusivity
Snurr et al.
11956 J. Pkys. Ckem., Vol. 98, No. 46, 1994
3
t 5
150000
t 10-l~
Eioi4 -
A
A 0
1 ''-I5
"
0
0.5
1.0
1.5
2.0
2.5
10-16
3.0
time (seconds)
tensor in ZSM-5 zeolites:
D,
a2
D,
+-Db2,
where a, b, and c are the unit cell lengths (20.07, 19.92, and 13.42 A, respectively). The correlation was derived by assuming that the time constant for the loss of "memory" of the diffusants is much smaller than the time of migration from one channel intersection to an adjacent one. The values of D,: predicted from eq 32 using our calculated values of D,, and Dyy are given in the last column of Table 3. The correlation yields values (LIZ?) that are in excellent agreement with our calculations (D Z J . The orientationally-averaged diffusivity calculated from eq 31 is compared with low-loading experimental data43-47 in Figure 4 in the form of an Arrhenius plot. The experimental results shown are for silicalite samples with very low aluminum content. Diffusivities of benzene in ZSM-5 (higher aluminum content) are reported to be a factor of 2-3 smaller than in ~ i l i c a l i t e .The ~ ~ diffusivities predicted by our TST calculations are 1-2 orders of magnitude smaller than the experimental values. Fitting the TST results to the expression
D = Do exp(-E,/RT)
+TST I " 0.002
"
'
"
0.003
I
' 0.004
"
' 0.005
1/ T (K')
Figure 3. Mean-squared displacement versus time for benzene in silicate at low loading and 300 K. Benzene diffuses more rapidly in the y direction (straight channels) than in x or z.
-c2_ -
0
Wu(1983) Shah(1988) Eic (1989) Van-Den-Begin (1989) Shen(1991)
(33)
yields a preexponential factor DOof 2.54 x cm2/s and an activation energy E, of 36.7 kJ/mol. Doing a least-squares fit of the data of Eic, Van-Den-Begin, and Shen gives a DOvalue of 5.33 x cm2/s and an activation energy of 25.2 kJ/mol, in line with the commonly cited value of Ea.3 If Shah's data are included in the least-squares fit, the values become 4.02 x cm2/s and 32.5 kJ/mol for DO and Ea, respectively. Compared to either analysis of the experimental diffusivities, our value of DO is too small and our activation energy is too large. The predicted activation energies in the x, y. and z directions are 37.8, 36.5, and 37.5 kJ/mol, respectively. The continuous-time/discrete-space Monte Carlo scheme does not require that one invoke a macrostate picture; it can be used with all of the microstates included. Such a simulation would
Figure 4. Comparison of TST predictions and experimental results for the self-diffusivity of benzene in silicalite at low loading.
both provide information on the short-time dynamics of the system and serve as a check on the macrostate picture. Because there are so many fast local transitions that occur very frequently, the average time step of such a simulation becomes very short, and the simulation spends most of its time tracking these fast local motions. The simulation must be run for very many time steps to explore long-range diffusive motion. This is reminiscent of the problem with using MD to study this system. The Monte Carlo scheme is still faster than MD, however, even if all microstates are included, because the Monte Carlo algorithm does not require any time-consuming potential energy evaluations. We performed a dynamic Monte Carlo simulation for benzene in silicalite at 500 K with all microstates included. The plot of mean-squared displacement versus time is shown in Figure 5. Only 100 (noninteracting) molecules were used, resulting in the increased noise compared to Figure 3. This simulation of 100 molecules for 6 ,us took about a week of computer time on a workstation, compared to an hour for the simulation with the macrostate picture, even though the macrostate simulation at 300 K used 1000 molecules and covered a time of 6000 p s . Simulations not invoking macrostates are not feasible at lower temperatures, as the separation of time scales between the fast motion and the slow motions becomes too large. The mean-squared displacements of Figure 5 appear to be large enough (compared to the relevant distance scales on the silicalite pore network) to allow a meaningful diffusivity to be extracted. The result is 3.4 x cm2/s, in good agreement with the value extracted from the macrostate picture and reported in Table 3. In addition to the macroscopic diffusion coefficient, the nature and mechanism of benzene motion within the zeolite are also of interest. To elucidate the nature of the short-time motions, computer graphics were used to visualize some of the calculated diffusion paths for motions within a channel intersection. (Recall that, at low loading, the molecules are thought to reside primarily in the intersection regions.) The translational and rotational trajectories of benzene along the diffusion paths were followed on the computer screen. As interpretation of NMR and QENS spectra suggest^,*^,^^ a number of paths resulted in rotations of the benzene molecule about its c6 axis. Rotations
J. Phys. Chem., Vol. 98, No. 46, 1994 11957 125 -
. ..
100
. ..
Total
--- Y
-
-.-.-
a -
5 c c
E -Q
75-
P
.-'0ln
-0
p!
a
50-
Ea
I
2 25 -
....,..._,......,...
c
I
Y
0
1
2
3
4
5
8
7
lm
lime (micro9econd$1
Figure 5. Mean-squared displacement for benzene in silicalite at low loading and 500 K. This simulation used only 100 molecules, resulting in increased noise compared to Figure 3. The macrostate picture was not invoked for this simulation; transitions between all microstates were included explicitly.
TABLE 3: Components of the Self-Diffusivity of Benzene in Silicalite at Low Loading.* 200 6.5 x 300 1.1 x 400 500
4.0 x
3.7 x
2.3 x 8.0 x 1.9 x lo-" 2.0 x
1.6 x
9.1 x
2.4 x
2.9 x
9.3 x 10-l'
6.9 x 7.0 x
8.4 x
9.0 x 2.7 x lo-" 7.1 x 7.2 x lo-"
t
"Values are reported in cm2/s. The values in the last column, 0:". were calculated from the correlation of Kiirger (eq 32) using the
U
values 0, and Dyy.
of 15", 30", and 60" were observed; most of these paths also involved some other slight motions of the molecule such as a retilting of the plane of the aromatic ring withm the intersection. One such path is illustrated in Figure 6; the initial and final configurations of the path are shown. Both of these configurations are local potential energy minima. The view is looking down the straight channels. One of the benzene hydrogen atoms is colored purple to show the rotation of the molecule about the c6 axis. The time scale for the c6 rotational motions is an s. T h s is consistent with the correlation the order of time for c 6 rotational motion measured by Jobic et aLi9 of 3.2 x s at 300 K at low loading. Several flips of the aromatic ring itself were also observed. Flips of 90" and 180" were seen within the channel intersection. Two of the 180" flips also included a 60" c6 rotation. These motions occur on a time scale of 10-10-10-9 s. Such motions are slower than the QENS experiments of Jobic et al. but would be considered fast on NMR time scales. Some of the paths corresponding to the infrequent transitions between macrostates were also visualized. One of the primary paths through the straight channels is shown in Figure 7. The view of the figure looks down the sinusoidal channels, and the straight channels run top to bottom in this view. In Figure 7a, the benzene molecule starts in a potential energy minimum in a channel intersection ( T I = -62.1 kT/mol). The molecule proceeds to the transition state shown in Figure 7b (T)= -24.0 kJ/mol) and then to the minimum shown in Figure 7c (Q'=
t I t 4 'II A -
m
U
1 Figure 6. Initial (a, top) and final
L
of the diffusion paths. The benzen intersection and undergoes * motion of the aromatic ring The view is down the straight cham atoms is colored purple to show the rotation of the molecule about the c 6 axis.
-39.9 kJ/mol). The molecule can continue into the straight channel by going through another transition state (not shown;
11958 J. Phys. Chem., Vol. 98, No. 46, 1994
Snurr et al.
I
x
-+,
Figure 7. Path leading from a channel intersection to a straight channel. The molecule starts at a potential energy minimum in the channel intersection (a, top left) and goes over a high-energy transition state (b, top right) to enter the straight channel (c, bottom left). The molecule can continue into the straight channel by going over another transition state (not shown) and then to another minimum (d, bottom right). See text for more details.
9"= -37.4 kJ/mol) and then to the minimum shown in Figure 7d
(V= -55.8 kJ/mol). Note that going over the energetic
barrier in Figure 7b requires a significant rotation of the benzene ring. In the channel intersection, the benzene c 6 axis points
Dynamics of Benzene in Silicalite generally in the z direction. (See also Figure 6.) As the 6 axis molecule moves through the straight channel, the c realigns toward the x direction. For the paths through the sinusoidal channels that we visualized, the c6 axis points generally in the y direction. Rotational motions appear to play a significant role in the long-range diffusive movement of benzene in silicalite.
Discussion There are several possible explanations for the large discrepancy between the experimental values of the diffusivity and the TST predictions. One possibility is a problem with the potential parameters. These same potential parameters, however, give good agreement with experiment for the thermodynamic properties of benzene in silicalite over the full range of sorbate loading.21 At infinite dilution, the predicted Henry's constant and isosteric heat match experimental values quite well, thus providing confidence in the potentials. One should therefore examine other possible sources of error first. A second possibility is errors introduced by using the harmonic approximation in evaluating the configurational integrals for obtaining the rate constants and microstate occupation probabilities. Errors in the evaluation of the configurationalintegrals would affect the preexponential factor of the rate constants and of the diffusivity. Above, however, we noted that our prediction of the preexponential factor, DO,is only off from the experimental value by a factor of 2. Below we will suggest computational techniques that are more accurate than the harmonic approximation. A very likely source of large errors in our modeling is the assumption that the zeolite is rigid. Zeolite framework flexibility has been shown to have only a small effect on the predicted diffusivity of small molecules in ~ i l i c a l i t e , but ' ~ ~for ~~ systems where the adsorbates fit very tightly in the pores, breathing motions or random vibrations of the zeolite might have an important influence on sorbate diffusion. Deem et al.48 have used a potential model for zeolite frameworks to simulate the dynamics of several zeolites, including silicalite. They find that the oxygen-oxygen distances across the pores generally show periodicity with time. The amplitude of the pore breathing is a function of temperature. For silicalite the distances across a pore varied up to 0.5 8, at 300 K. Such a variation at one of the transition states for diffusion could clearly have a large effect on the activation energy. The work of Deem et al. considered only the empty zeolite framework; the presence of a large adsorbate like benzene may change the breathing motions of the silicalite structure. Flexibility of the zeolite could be included within a TST calculation by using an approach introduced by Gusev and Suter for studying the dynamics of small penetrants in dense polymers.49 In this approach, atoms of the sorbent are assumed to move isotropically and independently of one another. Gusev and Suter explain how this allows for an efficient way of computing the potential of mean force experienced by the penetrant. The TST calculation is then carried out much as described above but using the potential of mean force instead of the potential energy. A critical parameter in this approach is the amplitude of the fluctuations experienced by the zeolite atoms. These values could possibly be extracted from the work of Deem et aL4*or through the approach used by Gusev and S ~ t e r When . ~ ~ this approach was applied to diffusion through polymers, diffusivities were predicted within a factor of 2 of the experimental values. This was a tremendous improvement over considering the polymer network to be rigid in the TST calculation, which leads to trapping of penetrants in some
J. Phys. Chem., Vol. 98, No. 46, 1994 11959 regions and diffusivities that are underestimated by 3-5 orders of magnitude. It is possible that the diffusion of benzene in silicalite depends on concerted openings of the bottlenecks in the pores rather than random fluctuations of the zeolite atoms. Applying the method of Gusev and Suter would capture effects of the latter but not of the former. Concerted breathing motions could, however, be modeled by incorporating degrees of freedom of the zeolite atoms explicitly in the TST calculations. Such an approach is currently being developed for diffusion in polym e r ~ . The ~ ~ searches , ~ ~ for minima and transition states would allow for relaxation of the zeolite atoms in addition to translation and rotation of the benzene molecule. This would, of course, require a potential model for intrazeolite interactions in addition to zeolitehorbate interactions. Tracing the diffusion paths and calculating the TST rate constants would also incorporate zeolite degrees of freedom. Above we mentioned the differences in the equilibrium macrostate occupation probabilities, piq, as calculated by the harmonic approximation at each microstate or by Monte Carlo integration over the macrostate. The significant differences in the values of pzq from the two methods are indicative of a fairly large discrepancy in the values of the configurational integrals from the two methods. Errors in calculating the configurational integrals will affect not only piq but also the rate constants. For systems of low dimensionality (e.g., a spherical penetrant in a rigid zeolite, having three degrees of freedom), the configurationalintegrals of eq 19 can be computed directly by Monte Carlo integration. June et al.30did this for xenon in silicalite. The three-dimensional configuration space was tessellated into voxels, and all voxels were assigned to one of the states. The configurational integral over a state was then evaluated by Monte Carlo integration over the voxels of that state. The configurational integral over the dividing surface between two states (numerator of eq 19) was evaluated by a Monte Carlo integration over the common faces of voxels on the border between two states. This technique is not easily extended to systems with more than three degrees of freedom, but other techniques could be used.36 For example, in a free energy perturbation s~heme,~O,~l a series of Monte Carlo or MD simulations could be carried out along the diffusion path, each constrained to a particular hyperplane normal to the diffusion path. This provides a free energy profile along the diffusion path and can be used to obtain the TST rate constants. Calculations with the harmonic approximation are much faster and should be carried out first to identify important paths where the more rigorous calculations should concentrate. Transition-statetheory assumes that all dynamical trajectories crossing the transition state in the direction from i to j undergo a successful crossing into state j . In reality, the system may fail to thermalize in the new state and may cross back into the original state. Another possibility is that the molecule may undergo a multistate transition and move directly to a second product state. The transition-state theory rate constant, k,,,TST, can, however, be corrected to account for these occurrences by computing a dynamical correction factor,
x-j.30,50
(34)
fi-,
is estimated from short MD runs that are initiated at the dividing surface between i and j . The fraction of the MD trajectories that ultimately thermalize in the state toward which they were moving at the beginning of the MD run gives the dynamical correction factor. The MD simulations required to estimate fi-, are not nearly as long as those that would be
11960 J. Phys. Chem., Vol. 98, No. 46, 1994 required to calculate the diffusivity from MD alone. The time scale for thermalization once the system is at the dividing surface is much shorter than the time scale between transitions. While concerted breathing motions of the zeolite or correlated multistate jumps may enhance diffusion of benzene through silicalite, random vibrations of zeolite atoms may decrease the dynamical correction factor and thus lower the diffusivity. The relative importance of these various effects should be examined in future work.
Conclusions We have presented a series of calculations that investigate the dynamics of benzene in silicalite at low loading using transition-state theory. Benzene diffuses through silicalite by infrequent hops between preferred adsorption sites, thus making transition-state theory applicable. The TST analysis can investigate the long time scales relevant to diffusion in this system. Predictions include the diffusion coefficient and the dynamical trajectories followed by the sorbate molecules. The calculations proceed through a series of stages. The first stage is the identification of all potential energy minima and saddle points. Diffusion paths connecting pairs of minima are constructed through each saddle point (transition state). The intrinsic reaction coordinate (IRC) approach was used to construct the diffusion paths in six dimensions using the equations of motion for a rigid body (benzene) through a static potential field (silicalite). To the best of our knowledge, this is the first time that such an approach has been applied to the diffusion of a molecule with internal structure through a microporous medium. A rate constant for each transition between minima (microstates) spanned by a diffusion path was calculated using a harmonic approximation to the potential energy function. The rate constants at a given temperature were found to span many orders of magnitude. Microstates mutually accessible by fast motions were lumped together to form macrostates, and effective rate constants for the slower motions between macrostates were then calculated. The macrostates were found to correspond to the diffusant residing in the straight channels, the sinusoidal channels, and the channel intersections. Transitions between microstates within a macrostate are orders of magnitude faster than transitions between macrostates, allowing the microstates within a given macrostate to be viewed as being in equilibrium with one another. The fastest motions were found to be relatively insensitive to temperature, while the rate constants for transitions between macrostates are strongly temperature dependent. These observations agree with experimental studies of this system. The predicted time scales for local motions within the channel intersections agree well with spectroscopic results. Many of these motions correspond to rotations of the benzene molecule about its axis. Given the rate constants, the self-diffusivity was computed with a dynamic Monte Carlo simulation. Benzene was seen to diffuse most readily along the straight channels, but diffusion in the sinusoidal channels is not negligible. Self-diffusivities at low occupancy were computed at 200, 300,400, and 500 K, and an activation energy of 36.7 kJ/mol was obtained. The predicted activation energy is larger than the experimental value, and the predicted diffusivities are 1-2 orders of magnitude smaller than the experimental values. Possible reasons for this discrepancy were discussed. One likely reason is our assumption that the zeolite is rigid; breathing modes of the zeolite pores might significantly enhance diffusion of the tightly-fitting benzene molecule. Additional calculations to investigate this point have been suggested, including TST calculations incor-
Snurr et al. porating a flexible representation of the zeolite and more rigorous calculations that do not necessitate invoking a harmonic approximation in evaluating the configurational integrals. If future work can predict the diffusivity for benzene in silicalite in better agreement with experiment, an important extension would be to apply transition-state theory to diffusion of the isomers of xylene through silicalite. This would allow a truly molecular-level look at the mechanisms involved in shapeselective transport and catalysis.
Acknowledgment. Support for this work was provided by the Advanced Industrial Concepts Division of the U S . Department of Energy under Contract DE-AC03-76SF00098 and by W. R. Grace and Co.-Connecticut. We thank Dr. Tim Robinson for his kind assistance with the computer graphics work and NIH Grant S10 lUtO5651-01 for support of the graphics facilities. Helpful conversations with M. L. Greenfield and J. S . Shaffer are gratefully acknowledged. Appendix To evaluate the configurational integrals for the TST rate constant in the harmonic approximation, one needs the eigenvalues of the Hessian matrix H,, where the derivatives are taken with respect to the mass-weighted generalized coordinates r (eq 22). Unfortunately, working with r is not convenient. Rather than calculating the eigenvalues of H,,, one might intuitively expect that one could use the eigenvalues of the matrix G-lH,,, since the eigenvectors of G-lH,, were used to locate the first step of the diffusion path starting at the saddle point. If these two approaches are equivalent, then the same eigenvalues E. must satisfy the following two equations simultaneously:
H,, d r = il d r
(35)
G-'H,, dq = il dq To show that these two equations are indeed satisfied simultaneously by the same eigenvalues and eigenvectors, we can transform eq 35: (37) or dq = il dq
($)-'Hrr($)
where adaq represents the Jacobian of transformation from q to r coordinates. The postulated simultaneous validity of eqs 36 and 38 with the same eigenvalues and eigenvectors would imply (39) The Jacobian can be obtained from the coefficients of the following equations (see eq 5): d d sin 8 sin ly
+ d e cos I$ = dQl
d 4 sin 8 cos ly - d e sin I$ = dS2, d 4 cos 8
+ dly = dS2,
Thus we obtain eq 41 (Chart 1).
(40)
Dynamics of Benzene in Silicalite
J. Phys. Chem., Vol. 98, No. 46, 1994 11961
CHART 1
ax
3 ae1
=I
0
0
o
o
ZrMO o A s i n e sin+ o
h c o s q
Using these expressions, it can indeed by shown that eq 39 is true and that H, and G - ' q q have the same eigenvalues. Thus, we can use the eigenvalues of G-*H,, to calculate the frequencies needed for evaluating kHA. Another approach with the same final result is suggested by Greenfield and T h e o d ~ r o u . Rather ~~ than having six degrees of freedom, the benzene molecule is described by the massweighted Cartesian coordinates, x, of the atoms of the molecule ( N = 36). Stiff harmonic potentials govern the bond lengths and bond angles. The Vineyard analysis now applies directly, and eq 20 gives kHA. The stiff degrees of freedom are next taken to the limit of infinite stiffness. This causes all but six of the eigenvalues of the Hessian H,, (i-e., all eigenvalues associated with the stiff degrees of freedom) to approach positive infinity in the same way at both the saddle point and the minimum; they thus cancel out the expression for kHA. The remaining frequencies are the square roots of the eigenvalues of the matrix G - l q q , giving the same result as derived above. Greenfield and Theodorou make a harmonic approximation in the mass-weighted Cartesian coordinates of the atoms, while above we made the harmonic approximation in the mass (or moment of inertia)-weighted generalized coordinates r. We only expect to get the same final result for the case of a rigid body with no intemal degrees of freedom. If the molecule has flexible torsional degrees of freedom, then regarding the stiff degrees of freedom as truly rigid is not the same as treating them as having stiff potentials in the limit of infinite stiffness."
References and Notes (1) van Bekkum, H., Flanigen, E. M., Jansen, J. C., Eds. Introduction to Zeolite Science and Practice; Elsevier: Amsterdam, 1991. (2) Wei, J. J. Catal. 1982, 76, 433-439. (3) K&ger, J.; Ruthven, D. M. Diffusion in Zeolites and Other Microporous Solids; Wiley-Interscience: New York, 1992. (4) Caro, J.; Jobic, H.; Biilow, M.; Kirger, J.; Zibrowius, B. In Advances in Catalysis;Eley, D. D., Pines, H., Weisz, P. B., Eds.; Academic Press: New York, 1993; Vol. 39, pp 351-414. ( 5 ) Nowak, A. K.; Cheetham, A. K.; Pickett, S. D.; Ramdas, S . Mol. Simul. 1987, 1, 67-77. (6) Pickett, S. D.; Nowak, A. K.; Thomas, J. M.; Cheetham, A. K. Zeolites 1989, 9, 123-128. (7) Talu, 0. Mol. Simul. 1991, 8, 119-132. (8) Vignt-Maeder, F.; Jobic, H. Chem. Phys. Letr. 1990, 169, 31-35.
(9) Schraer, K.-P.; Sauer, J. Z. Phys. Chem. (Leipzig)1990,271,289296. (10) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (11) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1990, 94, 8232-8240. (12) Demontis, P.; Fois, E. S.; Suffritti, G. B.; Quartieri, S. J . Phys. Chem. 1990, 94, 4329-4334. (13) Goodbody, S . J.; Watanabe, K.; MacGowan, D.; Walton, J. P. R. B.; Quirke, N. J . Chem. Soc., Faraday Trans. 1991, 87, 1951-1958. (14) Maginn, E. J.; Bell, A. T.; Theordorou, D. N. J. Phys. Chem. 1993, 97, 4173-4181. (15) June, R. L.; Bell, A. T.; Theodorou, D. N. 1.Phys. Chem. 1992, 96, 1051-1060. (16) Thamm, H. J . Phys. Chem. 1987, 91, 8-11. (17) Forste, C.; Gemanus, A,; Kiger, J.; Pfeifer, H.; Caro, J.; Pilz, W.; Zikanova, A. J . Chem. Soc., Faraday Trans. 1 1987,83, 2301-2309. (18) van Koningsveld, H.; Tuinstra, F.; van Bekkum, H.; Jansen, J. C. Acta Crystallogr. 1989, B45, 423-431. (19) Jobic, H.; Bee, M.; Dianoux, A. J. J . Chem. SOC.,Faraday Trans. I 1989, 85, 2525-2534. (20) Snurr, R. Q.; Bell, A. T.; Theodorou, D. N. In Proceedings from the Ninth International Zeolite Conference, Montreal 1992; von Ballmoos, R., Higgins, J. B.; Treacy, M. M. J., Eds.; Buttenvorth-Heinemann: Boston, 1993; Vol. II, pp 71-78. (21) Snurr, R. Q.; Bell, A. T.; Theodorou, D. N. J . Phys. Chem. 1993, 97, 13742-13752. (22) Snurr, R. Q.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1994, 98, 5111-5119. (23) Li, J.; Talu, 0. J . Chem. SOC.,Faraday Trans. 1993, 89, 16831687. (24) Zibrowius, B.; Caro, J.; Pfeifer, H. J . Chem. Soc., Faraday Trans. I 1988, 84, 2347-2356. (25) Hill, T. L. An Introduction to Statistical Thermodynamics; Dover: New York, 1986. (26) Sargent, R. W. H.; Whitford, C. J. In Molecular Sieve Zeolites II; Gould, R. I., Ed.; Advances in Chemistry Series 102; American Chemical Society: Washington, DC, 1971. (27) Ruthven, D. M.; Derrah, R. I. J . Chem. SOC.,Faraday Trans. 1 1972, 68, 2332-2343. (28) Betemps, M.; Jutard, A. J. Phys. D:Appl. Phys. 1980, 13, 423. (29) Kirger, J.; Pfeifer, H.; Haberlandt, R. J. Chem. SOC.,Faraday Trans. 11980, 76, 1569-1575. (30) June, R. L.; Bell, A. T.; Theodorou, D. N. J . Phys. Chem. 1991, 95, 8866-8878. (31) Olson, D. H.; Kokotailo, G. T.; Lawton, S. L.; Meier, W. M. J. Phys. Chem. 1981, 85, 2238-2243. (32) Goldstein, H. Classical Mechanics, 2nd ed.; Addison-Wesley: Reading, MA, 1980. (33) Fukui, K. Acc. Chem. Res. 1981, 14, 363-368. (34) Baker, J. J. Comput. Chem. 1986, 7, 385-395. (35) Greenfield, M. L.; Theodorou, D. N. Manuscript in preparation. (36) Sevick, E. M.; Bell, A. T.; Theodorou, D. N. J . Chem. Phys. 1993, 98, 3196-3212. (37) Vineyard, G. J. Phys. Chem. Solids 1957, 3, 121-127. (38) Oppenheim, I.; Shuler, K. E.; Weiss, G. H. Stochastic Processes in Chemical Physics: The Master Equation; MIT Press: Cambridge, MA, 1977. (39) Wei, J.; Prater, C. D. In Advances in Catalysis; Eley, D. D., Selwood, P. W.; Weisz, P. B., Eds.; Academic Press: New York, 1962; Vol. 13, pp 203-392. (40) Shaffer, J. S.; Chakraborty, A. K. Macromolecules 1993,26,11201136. (41) Tsikoyiannis, J. G.; Wei, J. Chem. Eng. Sci. 1991, 46, 233-255. (42) K&ger, 5. J . Phys. Chem. 1991, 95, 5558-5560. (43) Wu, P.; Debebe, A.; Ma, Y. H. Zeolites 1983, 3, 118-122. (44) Shah, D. B.; Hayhurst, D. T.; Evanina, G.; Guo, C. J. AIChE J . 1988, 34, 1713-1717. (45) Eic, M.; Ruthven, D. M. In Zeolites: Facts, Figures, Future; Jacobs, P. A., van Santen, R. A,, Eds.; Elsevier: Amsterdam, 1989; pp 897-905. (46) Van-Den-Begin, N.; Rees, L. V. C.; Caro, J.; Biilow, M. Zeolites 1989, 9, 287-292. (47) Shen, D.; Rees, L. V. C. Zeolites 1991, 11, 666-671. (48) Deem, M. W.; Newsam, J. M.; Creighton, J. A. J . Am. Chem. SOC. 1992, 114, 7198-7207. (49) Gusev, A. A,; Suter, U. W. J . Chem. Phys. 1993,99,2228-2234. (50) Theodorou, D. N. In Difision in Polymers; Neogi, P., Ed.; Marcel Dekker: New York, in press. (51) Elber, R. J. Chem. Phys. 1990, 93, 4312-4321.
