J. Phys. Chem. 1996, 100, 4571-4581
4571
Cage-to-Cage Jumps of Xenon in Zeolite NaY: Transmission Coefficients from Molecular Dynamics Simulations Thomas Mosell,† Gerhard Schrimpf,†,§ Christoph Hahn,†,⊥ and Ju1 rgen Brickmann*,†,‡ Institut fu¨ r Physikalische Chemie I and Darmsta¨ dter Zentrum fu¨ r Wissenschaftliches Rechnen, Technische Hochschule Darmstadt, Petersenstrasse 20, D-64287 Darmstadt, Germany ReceiVed: September 8, 1995X
The diffusion of xenon in zeolite NaY at low temperatures is studied by molecular dynamics simulations of individual trajectories started at the cage-to-cage boundary. Transmission coefficients are obtained from the normalized reactive-flux correlation function. Recrossings of the cage-to-cage boundary by xenon atoms at very low temperatures are attributed to instantaneous shifts of the maximum of the potential energy along the diffusion path. At higher temperatures there are increasingly many recrossings due to the slow energy transfer from guest atoms to the zeolite lattice: The xenon atom can pass several energy barriers during one activation period. The individual jumps are investigated by introducing an absorbing barrier surrounding the two supercages on both sides of the central cage-to-cage boundary.
I. Introduction Zeolites are microporous solids. They are used in the chemical industry as molecular sieves and as heterogeneous catalysts.1 The diffusion of guest molecules in zeolites has attracted a lot of attention in the literature.1,2 In addition to experimental techniques for the study of diffusion processes, the method of molecular dynamics (MD) simulations3 is increasingly applied (see, e.g., ref 4). Reports on the simulation of alkanes5-8 and rare gases9,10 in silicalite, xenon in zeolite NaY,11-14 argon in zeolite NaCaA,13,15 methane in zeolite ZK4,16,17 and benzene in purely siliceous faujasite18 are among the more recent publications. The time scale that is accessible by conventional MD simulations limits the applicability of this method to guest molecules whose diffusion is relatively fast. Bull et al.18 estimated that a simulation time of about 200 ns is necessary for a proper determination of the diffusion coefficient of benzene in zeolite NaY at room temperature. If a guest molecule remains for a considerable time at the same adsorption site, it is possible to describe the diffusion process as a sequence of jump events, i.e., the diffusion process can be described by a hopping model. Moreover, if subsequent jumps can be considered as independent, transition-state theory19 can suitably be applied to calculate the probability for a single jump.20,21 In conventional transition-state theory, it is assumed that every trajectory remains on the side of the products once it has passed through the transition state.22,23 Not all molecular systems comply with this condition, and the results obtained from transition-state theory are often corrected by the introduction of a transmission coefficient which takes into account the possibility of recrossing of trajectories to the educt site. The transmission coefficient can be determined from MD simulations by averaging over a large number of trajectories.24 Conventionally, the trajectories are started in the transition-state †
Institut fu¨r Physikalische Chemie. Darmsta¨dter Zentrum fu¨r Wissenschaftliches Rechnen. Present address: IBM Deutschland Informationssysteme GmbH, Wissenschaftliches Zentrum Heidelberg, Vangerostrasse 18, D-69115 Heidelberg, Germany. ⊥ Present address: AKZO Nobel Central Research Obernburg, D-63784 Obernburg, Germany. * To whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, February 15, 1996. ‡ §
0022-3654/96/20100-4571$12.00/0
region. This approach has been successfully applied to chemical reactions. Proton-25,26 and electron27-transfer reactions, ion-pair interconversion,28 SN129 and SN230,31 reactions, and conformational isomerization processes32-34 have been investigated within this model approach. Furthermore, the approach has been applied to diffusion processes with high activation barriers, e.g., the diffusion of crystalline defects in solids35 and of small molecules in proteins.36 Working on the diffusion of xenon and SF6 in silicalite, June et al.37 extended the range of applications to the diffusion of guest molecules in zeolites. They calculated transition probabilities for jumps between adsorption sites by dynamically corrected transition-state theory and subsequently calculated macroscopic diffusion coefficients by a Monte Carlo simulation on a lattice. In this work, the results of simulations that determine the diffusion of xenon atoms in zeolite NaY are presented, and transmission coefficients are reported. The trajectories of single xenon atoms in zeolite NaY were started in the transition state, and transmission coefficients were determined for different temperatures. Additionally, the energy relaxation after the barrier passing and the frequency of jump events that consist of several barrier passages have been studied. The possibilities to divide the diffusion process into separate jump events are discussed in detail. There are several symmetrically identical adsorption sites for xenon inside a supercage. The barrier heights for jumps from supercage to supercage and for jumps between adsorption sites in the same supercage are different. Our recent conventional MD simulations14 can be analyzed with respect to mean residence times of the atoms at an adsorption site and in a supercage. For a loading of 1 xenon per supercage, the mean residence times at 300 K are 12 and 53 ps for the adsorption site and the supercage, respectively. Thus, the cage-to-cage jumps are the rate-determining steps in the diffusion process. In this work only jumps between supercages were investigated. The activation barrier for the diffusion of xenon in zeolite NaY is comparably small, and as a result, the range of temperatures where the diffusion can sensibly be described by a hopping model is restricted to low temperatures. Xenon in zeolite NaY was nevertheless chosen as the model system, because the diffusion of xenon in this zeolite has already been investigated by conventional MD simulations,11-14 and an © 1996 American Chemical Society
4572 J. Phys. Chem., Vol. 100, No. 11, 1996
Mosell et al. successfully to zeolite A.44 Schrimpf et al.41 modified the NaNa and Na-O interactions for zeolite NaY in order to avoid a migration of Na ions from site SII into the sodalite cage. In this work, the latter version of the force field was used. For a complete list of the parameter values, the reader is referred to ref 41. The potential energy of the system can be divided into a part Vzeo that depends only on the atoms of the zeolite lattice and a part Vxenon that contains all interactions connected with a guest molecules:
V ) Vzeo + Vxenon
(1)
The potential energy Vzeo of the zeolite lattice is calculated as the sum of two-body interactions, which are taken into account for the following combinations of atomic species:
Vzeo ) ∑VT-O + ∑VO-Na + ∑VO-(T)-O + ∑VNa-Na ij
Figure 1. Sketch of a supercage of zeolite NaY. The oxygen linkages between the tetrahedrally coordinated atoms are marked by lines. Cation positions are marked by black spheres. SI and SII sites are located at the center of the hexagonal prism and within the supercage (in front of the sodalite cage), respectively.
analysis of the transition-state behavior was presented by Yashonath and Santikary.12 Moreover, it will be possible to compare simulation results with experimental diffusion coefficients obtained for xenon in zeolite NaX by pulsed field gradient NMR.38,39 The paper is organized as follows. Section II contains the force-field representation of the model system. The determination of transmission coefficients and the definition of the transition state are described in section III. Some details on the calculations are reported in section IV. In section V, the results of the MD simulations of the cage-to-cage jumps are presented and discussed. Conclusions are stated in section VI.
ij
ij
(2)
ij
The T-O bonds, the interactions between Na and O, and the Urey-Bradley term for the interaction O-(T)-O of two oxygens separated by a T atom are modeled by anharmonic bonded potentials between atom i of type x and atom j of type y:
Vx-y(rij) ) V0xy + (kxy/2)(rij - r0xy)2 + (Axy/6)(rij - r0xy)3 (3) where rij is the distance between the two atoms. A cutoff at a distance rcut
rcut ) r0xy - 2kxy/Axy
(4)
is applied to these interactions. A repulsive term for the nonbonded Na-Na interactions is used in addition to the bonded potentials. The Na-Na interaction is suitably modeled by a Coulomb potential: 2
II. Model The framework of zeolite NaY consists of cuboctahedral sodalite cages that are linked by six-membered rings of oxygen atoms (see Figure 1). This structure contains large cavities called supercages. The supercages are tetrahedrally connected by windows framed by rings of 12 oxygen atoms and 12 tetrahedrally coordinated Si and A1 atoms. The cubic unit cell of zeolite NaY with space group Fd3hm and cell parameter a ) 24.85 Å40 contains eight supercages whose centers form a diamond lattice. The aluminum content of the zeolite lattice is variable. At a Si/Al ratio of 3:1, one unit cell consists of 144 Si, 48 Al, 384 O, and 48 Na atoms. The simulations are performed with a flexible zeolite lattice. The model used in these simulations has already been applied and described in detail before.41 No differentiation is made between Si and Al atoms, which are described by a joint atom type T with the mass of the silicon atoms. The simulations were performed at a Si/Al ratio of 3:1. Sites SI and SII were completely occupied by sodium ions. Other sites40 are not considered. Zeolite lattices are often kept rigid in simulations.4,5,8 In general, however, the trajectory of the xenon atom is influenced by the energy exchange with the lattice.41,42 Therefore, lattice flexibility is taken into account in our simulations. The atoms of the zeolite lattice may act as a heat bath.41 Demontis et al.43 developed a valence force field for natrolite and applied a slightly altered version of this force field
1 qNa VNa-Na ) 4π0 rij
(5)
The parametrization of this intermolecular interaction resulted in qNa ) 0.65e.41 The intermolecular interactions Vxenon between the xenon atom and the atoms of the zeolite lattice are modeled not only by Lennard-Jones potentials, but also by the induction energy Vpol that results from the strong polarizability of the xenon atoms:
Vxenon ) ∑VXe-NaLJ + ∑VXe-OLJ + Vpol ij
(6)
ij
All interactions between xenon and the T atoms of the zeolite lattice are neglected, because the tetrahedral atoms are wellshielded by the oxygen atoms. The parameters for the LennardJones potentials
VXe-y(rij) ) 4Xe-y[(σXe-y/rij)12 - (σXe-y/rij)6]
(7)
were taken from Kiselev and Du45 who had modified interaction parameters obtained from the Kirkwood-Mu¨ller equation in order to reproduce experimental Henry constants. Shifted-force potentials with a cutoff of 10 Å were used for the nonbonded interactions to reduce the CPU time for the calculation of the intermolecular interactions. The use of shifted-force potentials instead of the Ewald summation technique46 may seem critical for a periodic structures like the zeolite
Cage-to-Cage Jumps of Xenon in Zeolite NaY
J. Phys. Chem., Vol. 100, No. 11, 1996 4573
lattice. However, Dufner et al.47 showed that shifted-force potentials can be adequately used in simulations of zeolites. Point charges of qi ) +1e on the Na atoms and of qi ) -1/8e on the O atoms48 are taken for the calculation of the induction energy Vpol of the xenon atom:
where ri(t) is the position of the system on the reaction coordinate and rq is the point that defines the transition state. The number of systems in state A is given by
NA(t) ) ∑θ(rq - ri(t))
(14)
i
3 1 Vpol ) - ∑∑qiTβimuβ 2 i β
(8)
At the macroscopic level, an exponential decay of the correlation function eq 12 is expected:51,52
C(t) ) exp(-t/τ)
where µβ is the component of the induced dipole moment along coordinate β and Tβi is
1 1 ∇ Tβ ) 4π0 β|rXe - ri| i
(9)
(15)
In a two-state system, the relaxation time τ is linked to the rate constant k for transitions from A to B by
τ ) (1 - xA)/k
(16)
The dipole moment is calculated from
muβ ) -RXe∑Tβiqi
(10)
i
where RXe is the isotropic polarizability of xenon, for which the value RXe ) 4.01 Å3 of Bezus et al.48 was used. In this work, an iterative calculation of the induction energy, as suggested by A Ä ngya´n et al.49 and applied in ref 41, is not necessary, since only one xenon atom is contained in the simulation box. III. Simulation Technique In transition-state theory19 the rate constant kTST is determined solely on the basis of statistical information. Wigner22 showed that the basis of transition-state theory is equivalent to the assumption that no recrossings are possible after a trajectory has passed the transition state.23 Deviations from transitionstate theory can be taken into account by the introduction of a transmission coefficient κ:
k ) κkTST
(11)
where k is jump rate. III.1. Transmission Coefficients. Transmission coefficients can be determined from MD simulations of independent trajectories that are started in the transition state to avoid extensive sampling of nonreactive events.24 At least around 100 independent trajectories must be calculated to obtain transmission coefficients with acceptable accuracy.27,29,33,36 The calculation of the transmission coefficient can be based on an analysis of the history of each trajectory30 according to the stable-states picture of chemical reactions of Northrup and Hynes.50 However, this approach involves the a priori definition of a boundary between a transition-state region and stable states. A time-correlation-function formalism that was developed by Chandler51 is usually used instead.28,29 Chandler’s approach51 is based on the time-autocorrelation function C(t) for the fluctuations δNA(t) ) NA(t) - 〈NA〉 of the number of molecules NA in state A:
C(t) )
〈δNA(0)δNA(t)〉 〈(δNA(0))2〉
(12)
The brackets denote statistical averaging. The assignment of a molecule to state A or B of a two-state system can be formally described by the Heavyside function θ:
θ(rq ) -ri(t)) )
{
q
1 r g ri(t) 0 rq < ri(t)
(13)
where xA is the equilibrium value of the mole fraction of molecules in state A. Chandler derived an expression for the reactive flux k(t) in the form51
k(t) ) 〈r˘ (0) rˆδ(r(0) - rq) θ(r(t) - rq)〉/xA
(17)
where rˆ is the unit length of the reaction coordinate in reciprocal space, and δ denotes the Dirac delta function. The reactive flux k(t) initially decays, because a part of the trajectories recrosses the transition state. After some time, a plateau is obtained. When the plateau is reached, only the part of the flux that has resulted in the formation of products contributes to the reactive flux. The remaining reactive flux corresponds to the rate of transitions between the two states and, therefore, to the macroscopic rate constant k. The reactive flux k(t) can show a slow decay instead of a true plateau. This continuous decrease is a result of the occurrence of some independent rare events on the time scale of the simulations. According to Berne et al.53 the rate constant k may, nevertheless, be obtained, if the reactive flux converges against a linear asymptotic limit (LAL). Following these authors, the rate constant k is obtained in this case from the intersection between the vertical axis at t ) 0 and the extrapolation of the linear decay function LAL. If there is no separation between the time scale of recrossings during the same rare event and the time scale of independent rare events, the reactive flux k(t) will not show a linear behavior asymptotically. In this case, the two states model is no longer applicable.28,51,54 This limits the applicability of the approach to systems with activation barriers of at least a few kBT.28 The transmission coefficient can be determined from the extrapolation of the asymptotics of the normalized reactive-flux correlation function κ(t) for t ) 0 in the same manner as the transition rate.28,54 Since no recrossings can have occurred at t ) 0, k(0) and kTST must be equal and κ(t) can be obtained from the ratio k(t)/kTST.28,54 The rate constant kTST depends on quantities at t ) 0 only, and corresponds to51,52
kTST ) 〈r˘ (0) rˆδ(r(0) - rq) θ(r˘ (0))〉/xA
(18)
Finally, κ(t) is given by28,54
κ(t) )
〈r˘ (0) δ(r(0) - rq) θ(r(t) - rq)〉 〈r˘ (0) δ(r(0) - rq) θ(r˘ (0))〉
(19)
The velocity component r˘ (0) along the reaction coordinate in both numerator and denominator of eq 14 indicates that fluxweighted averaging has to be carried out. It is statistically more efficient to take this into account by using a flux-weighted
4574 J. Phys. Chem., Vol. 100, No. 11, 1996
Mosell et al.
distribution of r˘ (0), which leads to an equally weighted contribution of each trajectory to the results, than by sampling with a standard Boltzmann distribution, which would lead to a domination of a few trajectories.29,32 The coordinates of the starting configurations are determined from a simulation where the system is constrained in its transition state. In this simulation, which has been termed “constrained reaction coordinate dynamics” (CRCD),55 r˘ (0) constantly remains zero due to the application of a constraint for the reaction coordinate. Carter et al.55 pointed out that the constraining of the velocity along the reaction coordinate can change the probabilities for the different configurations in the transition state. However, this does not occur when only distance constraints are used.28 The values for r˘ (0) in the starting configuration are taken from a flux-weighted distribution. As a result, κ(t) is obtained from a sum over all trajectories:31
κ(t) )
1 N+
N+
1
(20)
-
where κ(t) has been divided into the contributions that arise from N+ trajectories with a positive and N- trajectories with a negative velocity component in the transition state. Note that the flux-weighted averaging results in a mean kinetic energy of 2kBT.59 III.2. Characterization of the Transition State. Our reaction coordinate for the cage-to-cage jumps in zeolite NaY points along the vector [111] in the unit cell, i.e., each point on the reaction coordinate corresponds to a plane with its normal vector n along [111]. A xenon atom at a position p is located on the plane
n‚(p - a) ) d
σ(p,a) ) n‚(p - a) - d ) 0
(22)
The stochastic thermalization method of Kast et al.57 was used to obtain a canonical distribution. After an equilibration period of 100 ps, the system was integrated for 1.1 ns, during which starting conditions for the individual trajectories were stored every 10 ps. Subsequently, the individual trajectories were followed for 35 ps (70 ps at 210 K) in both directions of time. The 220 half-trajectories were performed at constant energy.
N-
θ(r(t) - r ) - ∑θ(rq - r(t)) ∑ N i)1 i)1 q
simulations were performed at constant volume with periodic boundary conditions. Simulations with a box that contained one unit cell were performed for the temperatures 77, 100, 140, 160, 180, and 210 K. First, a CRCD simulation was performed where the system was constrained to the transition state. In this simulation, the xenon atom is confined to a plane due to the constraint σ:
(21)
that lies at a distance d from point a. The variable d and a constant point a are used to define the reaction coordinate. In this work, point a coincides with the center of the window that is located at the coordinates 1/2 1/2 1/2 of the crystal structure. As a result, the coordinate point d ) 0 Å coincides with the window plane. To consider the motion of the lattice atoms, the instantaneous coordinates of point a are obtained from the center of gravity of the lattice atoms in the simulation box. It should be noted here that the definition of the window plane by d ) 0 Å is not identical to that in ref 41 where the plane was calculated from the positions of oxygen atoms in the 12-T ring. The normal vector n is defined in the coordinate system of the simulation box. The definition within the coordinate system of the box is possible, because there is no rotation of the zeolite lattice with respect to the coordinate system of the box during the simulations. The window plane d ) 0 Å is used to define the transition state, since the boundary between two supercages is the natural choice for the transition-state plane for cage-to-cage jumps. Moreover, the maximum of the minimum-energy path along the reaction coordinate for xenon in zeolite NaY is located in this plane.14 At low temperatures, the potential of mean force along the reaction coordinate should not deviate too much from the corresponding potential-energy curve. As a result, the transition state is close to the plane that corresponds to the maximum of the potential of mean force along the reaction coordinate. IV. Details of the Simulations In the MD simulations, the equations of motion are integrated with the Verlet algorithm56 with a time step of ∆t ) 1 fs. The
V. Results and Discussion V.1. Transmission Coefficients. The normalized reactiveflux correlation function κ(t) for cage-to-cage jumps of xenon in zeolite NaY is depicted in Figure 2. The results for the temperatures from 77 to 180 K are displayed in Figure 2a. Even at the lowest temperature, the asymptotic linear limit is reached already after a period of 10 ps. The length of this period increases when temperature is raised and becomes significantly longer than the 10 ps period as was found by June et al.37 for their flux correlation functions. At 210 K, the period of 35 ps is not long enough to reach an asymptotic linear limit. Therefore, the simulations were extended to 70 ps at this temperature. The results are displayed in Figure 2b. A LAL is obtained after ca. 25-30 ps. A slow decay in κ(t) is observed for all temperatures except at 77 K. Some independent rare events occur during the simulations, and accordingly, one part of the xenon atoms will leave the two supercages that are on either side of the central transition state through one of the other windows. As a result, xenon atoms can pass through the (infinitely extended) transition-state plane while they are in other supercages. To avoid any influence of these passages on the transmission coefficient, trajectories are continually assigned to the same side of the transition-state plane once the xenon atom has left the twosupercage region. The temperature dependence of the transmission coefficient is displayed in Figure 3. This coefficient slowly decreases from κ ) 0.91 at 77 K to 0.77 at 210 K. At all temperatures, the transmission coefficient for cage-to-cage jumps of xenon in zeolite NaY is relatively high. Only a few trajectories return through the transition state. Therefore, transition-state theory turns out to provide the correct order-of-magnitude of the results. An analysis of MD simulations of xenon in zeolite NaY by reactive-flux correlation functions has already been presented by Yashonath and Santikary.12 They performed conventional MD simulations with a fixed lattice. Although a LAL is clearly to be seen in their correlation functions, they state that the mean residence time for a xenon atom in a supercage is too short for a separation of independent jump events from consecutive passages though the transition state during the same barrier passage. They determined a mean residence time from all events, whereas the condition that the mean residence time must be significantly larger than the time until a LAL is reached applies only to a mean residence time for events where a stable state was reached in the new minimum. Therefore, there are no indications that hopping models are invalid for xenon in zeolite NaY at low temperatures.
Cage-to-Cage Jumps of Xenon in Zeolite NaY
J. Phys. Chem., Vol. 100, No. 11, 1996 4575
Figure 2. Time dependence of the normalized reactive-flux correlation function κ(t) for cage-to-cage jumps of xenon in zeolite NaY (a, top) for temperatures in the range 77-180 K and (b, bottom) at 210 K. The determination of the transmission coefficient κ by extrapolation of the LAL to t ) 0 ps is indicated for the results at 210 K.
V.2. Minima in the Supercages. The history of some typical trajectories, which were all taken from the simulations at 140 K, has been depicted in Figure 4. Trajectories that remain on the “side of the products” are depicted in the upper and central part of Figure 4. The vibrations at d ≈ 3 Å and d ≈ 8 Å are caused by xenon atoms near the xenon adsorption sites in front of 4-T rings of the hexagonal prism. This location for the adsorption site is consistent with the results obtained in other simulations.11,14 At the minimum near d ≈ 8 Å, the vibrations extend to d ≈ 6 Å, because the xenon atoms move along a sequence of 4-T rings. The high mobility of xenon along this sequence has already been described earlier.14 These vibrations are invisible for the adsorption sites at d ≈ 3 Å, because the
sequence of 4-T rings is nearly parallel to the reaction coordinate in this case. Trajectory types A, B, and C correspond to simple jumps between adsorption sites in neighboring supercages. Trajectories that include an additional jump between two adsorption sites in the same supercage (as in trajectory D) occur frequently at higher temperatures. For that reason, a classification of trajectories as belonging to type A, B, or C becomes more and more ambiguous when the temperature is raised. An investigation into a preferential occupation of specific adsorption sites after cage-to-cage jumps should, therefore, be based on a different criterion. The ratio n1/n2 of xenon atoms near 3 Å to those near 8 Å after a period of 35 ps is listed in
4576 J. Phys. Chem., Vol. 100, No. 11, 1996
Mosell et al. TABLE 2: Number of Recrossings for Types F and G + Ha no. of trajectories
a
Figure 3. Temperature dependence of the transmission coefficient κ for cage-to-cage jumps of xenon in zeolite NaY.
Figure 4. Position d on the reaction coordinate as a function of time t for some typical trajectories A-H, which are used to characterize different trajectory types in the text. All trajectories in this figure were taken from the simulations at 140 K.
TABLE 1: Ratio n1/n2 of Xenon Atoms at Adsorption Sites near d ) 3 Å to Those at Adsorption Sites near d ) 8 Å (after 35 ps) T [K]
n1/n2
T [K]
n1/n2
77 100 140
4.2 2.5 2.2
160 180 210
1.5 1.6 1.8
Table 1. Trajectories that left the supercage, e.g., trajectory E, are ignored in this analysis. At low temperatures, the adsorption sites next to the transition-state plane are strongly favored. The ratio n1/n2 decreases with increasing temperature not only because of a larger number of trajectories of type D but also due to an increasing number of xenon atoms that reach one of the remote adsorption sites in the supercage directly. Thus, the tendency toward a uniform population of the adsorption sites becomes apparent on the time scale of the simulations at higher temperatures. This equilibration process within the supercage confirms the approach to
T [K]
type F
type G + H
77 100 140 160 180 210
7 8 3 0 1 2
3 6 15 15 24 28
The trajectory types are defined in the text.
restrict the investigation of jumps of xenon in zeolite NaY to cage-to-cage jumps. V.3. Recrossing Trajectories. Trajectories of xenon atoms that return through the transition state into their original supercage are displayed in the lower part of Figure 4. Trajectories of type F recross immediately after the initial passage of the transition state. Trajectories G and H get into the proximity of the adsorption sites before a second barrier transition takes place. The rate of occurrence of these trajectory types is listed in Table 2 where trajectories that penetrate more than 1 Å into the new supercage are classified as type G + H. No trajectory of type F got further than 0.26 Å on the reaction coordinate. Neither did any of these trajectories remain in the new supercage for more than 3 ps. Type F is mainly observed at very low temperatures. The type predominantes (among the recrossing trajectories) at 77 K but is rarely observed above 140 K. The predominance of type F at low temperatures is responsible for the fact that, during the first 5 ps, the reactiveflux correlation function κ(t) (Figure 2) decreases faster at low temperatures than at high temperatures, although the LAL value of κ(t) increases when the temperature is decreased. The recrossings of type G + H occur after collisions between the xenon atom and the lattice in the region of the adsorption sites. Some xenon atoms are reflected once and recross the transition state directly afterward (trajectory G), while others perform a few vibrations at one of the adsorption sites, before the recrossing occurs (trajectory H). The sharp boundary between type F and type G + H results from the increasing acceleration of the xenon atom when it proceeds along the reaction coordinate after passing through the transition state. The influence of the velocity component d˙ (0) and the force component Fd(0) along the reaction coordinate on the trajectory is shown in Figure 5 for the trajectories at 77 K. Recrossings of type F and type G + H are indicated by circles and rhombi, respectively. A positive value of the force denotes a force component Fd(0) that is parallel to the velocity component d˙ (0). Type F is characterized by a force component in the direction opposite from that of the velocity. Furthermore, the velocity component d˙ (0) is relatively small. The same behavior is observed for type-F trajectories at higher temperatures. The window plane coincides with the potential-energy maximum of the minimum-energy path along the reaction coordinate. However, the location of the potential-energy maximum on the reaction coordinate can be different for a xenon atom that does not follow the minimum-energy path. Moreover, the potential energy along the reaction coordinate can be influenced by vibrations of the lattice atoms. As a result, xenon atoms may have to surmount an additional energy barrier. The xenon atom will not pass the overall barrier, if its kinetic energy is insufficient to compensate for the additional barrier.30 The decrease in type-F trajectories at higher temperatures can be attributed to the increase of the kinetic energy at these temperatures.
Cage-to-Cage Jumps of Xenon in Zeolite NaY
J. Phys. Chem., Vol. 100, No. 11, 1996 4577
Figure 5. Occurrence of recrossings of type F (circles) and type G + H (rhombi) for different initial conditions (velocity component d˙ (0) and force component Fd(0) along the reaction coordinate) in the simulations at 77 K. Trajectories without a recrossing are marked by crosses.
TABLE 3: Mean Velocity Component 〈|d4 (0)|〉 along the Reaction Coordinate 〈|d˙ (0)|P) [m s-1]a T [K]
recross traj
leaving traj
all traj
77 100 140 160 180 210
30.8 36.6 78.3 74.1 79.7 71.9
58.4 67.6 82.3 85.7 82.7 87.7
49.4 58.2 71.8 75.2 76.1 80.1
a First column: recrossing trajectories. Second column: trajectories where xenon leaves the two supercages next to the central transitionstate plane. Third column: all trajectories.
The trajectories of type G + H do not show any apparent correlation with the starting conditions. Recrossings are observed even for large positive values of the force component Fd(0). In Table 3, the velocity components 〈|d˙ (0)|〉 of recrossing trajectories and trajectories that leave the supercage through one of the other windows are compared with the mean velocity component of all trajectories. While the mean velocity component 〈|d˙ (0)|〉 of the trajectories that leave the two-supercage region is enlarged, the differences between 〈|d˙ (0)|〉 of the recrossing trajectories and 〈|d˙ (0)|〉 of all trajectories do not exceed 10% in the temperature range from 140 to 210 K (where type G + H predominates among the recrossing trajectories). Moreover, the deviations are not systematic in this temperature range, and accordingly, the occurrence of a recrossing cannot be anticipated from the starting conditions. V.4. Energy Relaxation. The length of the period until a LAL is reached may depend on the rate of energy transfer from the xenon atom to the lattice. The time dependence of the average energy EXe
EXe ) Ekin + VLJ + Vpol
(23)
for a xenon atom after its passage through the transition state is depicted in Figure 6. This energy includes the kinetic energy Ekin of the xenon atom, the contributions from the induction energy Vpol, and the Lennard-Jones interactions VLJ between xenon and the atoms of the zeolite lattice. The decrease in EXe is a slow process and energy relaxation remains incomplete on the time scale of 35 ps. The curves of EXe differ only in the mean energy of the starting configurations. Thus, the lattice takes up the energy at a rate that is approximately independent
from temperature. The time scale of the energy relaxation is in accordance with results from simulations of kinetic-energy relaxation in nonequilibrium systems of xenon in zeolite NaY.41 A higher mean energy of the starting configurations, therefore, implies a longer period until the xenon atom has lost so much energy that a new barrier passage is energetically impossible. This leads to a greater probability for additional barrier passages and explains the decrease of the transmission coefficient at higher temperatures. The slow decrease of EXe is responsible for the length of the period until a LAL in κ(t) is reached. Therefore, the occurence of trajectories of type G + H is a consequence of the slow energy transfer to the lattice. The mean energy EXe in Figure 6 decreases monotonously because of the averaging over different trajectories. A typical individual trajectory is displayed in Figure 7. The xenon atom can repeatedly lose and gain energy in the collisions with lattice atoms. As a consequence, the mean energy EXe is not suitable quantity to anticipate the length of the period until deactivation for an individual trajectory and a more detailed investigation into correlations between the mean energy EXe and the length of the period until the LAL is reached is not appropriate. The different contributions to EXe are depicted in Figure 8 for the trajectories at 140 K. The Lennard-Jones VLJ interactions show a gradual, slow decrease during the complete time interval. In contrast, the induction energy Vpol decreases rapidly at the beginning, because the polarization of the xenon atom increases strongly near the adsorption sites due to the proximity of a sodium ion. Thus, the dynamic behavior of the system xenon in zeolite NaY shortly after a barrier crossing is strongly influenced by the high polarizability of the xenon atom. At first, a large part of the potential energy is converted to kinetic energy Ekin of the xenon atom. The kinetic energy increases strongly at the beginning. After 35 ps, the kinetic energy has decreased to approximately the value in the starting configurations. However, this is no equipartition of energy between all degrees of freedom, because the velocities of the xenon atoms were chosen according to a flux-weighted Boltzmann distribution, and the mean kinetic energy at t ) 0 ps is 2 kBt. Therefore, the relaxation of the kinetic energy Ekin is still incomplete at the end of the 35 ps simulations. V.5. Jump Events to Distant Minima. The slow energetic deactivation of the xenon atoms leads to the occurrence of trajectories with several barrier passages during a single activation period. Such a sequence of jumps, which we will call a “jump event”, can end in a more distant supercage. These jump events can be considered by calculating transition probabilities for every possible combination of minima.37,54 Alternatively, an absorbing boundary can be introduced to examine the occurrence of these jump events. For these investigations, we defined absorbing boundaries that are located in all windows that link the two supercages on either side of the central transition-state plane with other supercages. Straub and Berne58 applied an absorbing boundary in order to absorb trajectories that return through the transition state. Their approach has been criticized, because the conditions in the transition state during the first and second passage may be strongly correlated due to a slow relaxation of the local surroundings.31 This does not apply to an absorbing boundary that is located in another transition state, because the local surroundings (e.g., vibrations of the lattice atoms, position of xenon in the window plane) are uncorrelated. Therefore, the usage of an absorbing boundary to define a two-supercage region around the central transition state is possible, although the data in Table 3 revealed that the probability for a passage through
4578 J. Phys. Chem., Vol. 100, No. 11, 1996
Mosell et al.
Figure 6. Time dependence of the mean energy EXe of a xenon atom after the passage through the transition state for different temperatures.
Figure 7. Time dependence of the energy of a xenon atom (kinetic energy, induction energy and contributions from the Lennard-Jones interactions) for an individual trajectory taken from the simulations at 140 K.
the absorbing boundary depends on the velocity component d˙ (0) of the starting configuration. In this work, the occurrence of jump events to distant minima is investigated by considering a function χ(t) defined as the part of the trajectories that has been “absorbed” at the boundaries of the two-supercage region. If the diffusion process can be divided into separate jump events, this absorbing-boundary function χ(t) will contain a “plateau” that corresponds to the LAL of the normalized reactive-flux correlation function κ(t). We suggest to define an absorbing boundary factor χ that is obtained by extrapolating the LAL to t ) 0, creating an analogue to the transmission coefficient κ. The absorbing boundary factor χ states which part of the trajectories is absorbed during a barrier passage that has to be assigned to the original jump event. In Figure 9, the absorbing-boundary function χ(t) is displayed as a function of the time since the start of the simulation in the transition state. At the beginning, χ(t) remains at zero, because a xenon atom must travel the distance between two windows before it can contribute to χ(t). Subsequently, χ(t) increases strongly and reaches a LAL on the same time scale as κ(t). At
Figure 8. Time dependence of the contributions to the mean energy EXe of a xenon atom after the passage through the transition state at 140 K. These contributions are (from top) the kinetic energy Ekin, the Lennard-Jones interactions ELJ, and the induction energy Epol.
210 K, the occurrence of a LAL is verified by the simulation of 70 ps. Although it would in principle be possible to extend the examination of κ(t) for LALs to even higher temperatures, a further increase of temperature does not seem to warrant the necessary computer time. The absorbing-boundary factors χ that were determined from χ(t) are listed in Table 4. The values increase strongly when temperature is raised. However, the absorbing-boundary factor is significantly smaller than χ ) 1 at all temperatures, and accordingly, only a limited part of the trajectories includes a
Cage-to-Cage Jumps of Xenon in Zeolite NaY
J. Phys. Chem., Vol. 100, No. 11, 1996 4579
Figure 9. Absorbing-boundary function χ(t), representing the part of trajectories that has reached the absorbing boundary of the two-supercage region, (a, top) for temperatures in the range 77-180 K and (b, bottom) at 210 K. The determination of the absorbing-boundary factor χ by extrapolation of the LAL to t ) 0 ps is indicated for the results at 210 K.
TABLE 4: Absorbing-Boundary Factor χ for Xenon in Zeolite NaY T [K]
χ
T [K]
χ
77 100 140
0.04 0.10 0.15
160 180 210
0.19 0.22 0.34
jump event with several barrier passages. Therefore, the diffusion of xenon in zeolite NaY may be divided into separate activation periods. The applicability of a hopping model to the diffusion process of xenon in zeolite NaY is confirmed.
V.6. Investigation into a Centralized Mode for Cage-toCage Diffusion. Yashonath and Santikary12 suggested that there are two different pathways for the diffusion of xenon in zeolite NaY. In their simulations the xenon atoms constantly remained near the walls at low temperatures, whereas an increasing part of the xenon atoms passed from the window into the central part of the supercage at higher temperatures. Yashonath and Santikary concluded that not only a surface-mediated but also a centralized mode contributes to cage-to-cage diffusion. The numbers nCEN of half-trajectories that included a passage of the xenon atom through the central part of the supercage are
4580 J. Phys. Chem., Vol. 100, No. 11, 1996
Mosell et al.
TABLE 5: Passages through the Central Part of the Sueprcagea,b T [K]
nCEN
χCEN(35 ps)
χ(35 ps)
77 100 140 160 180 210
15 26 47 62 79 91
0.13 0.12 0.23 0.29 0.29 0.38
0.05 0.11 0.19 0.24 0.27 0.38
a The central part of the supercage is defined by a sphere with a radius of 3 Å around the center of the supercage. b nCEN ) number of the half-trajectories with a passage through the central part of the supercage; χCEN(t) ) absorbing-boundary function for half-trajectories with a passage through the central part of the supercage; χ(t) ) absorbing-boundary function for all half-trajectories. The absorbingboundary function is defined in the text.
listed in Table 5 for our simulations. Following Yashonath and Santikary, we defined the central part of the supercage by a sphere with the radius 3 Å. Only a few of the trajectories pass through the central part of the supercage at low temperatures and nCEN increases strongly when temperature is raised. Therefore, the tendency of the trajectories to pass through the central part of the supercage at higher temperatures is also found in our simulations. A xenon atom in an energetically deactivated state must be located near the adsorption site in proximity to the wall. Therefore, a separate diffusion mechanism through the central part of the supercage mainly concerns successive cage-to-cage jumps during the same activation period. As a result, the probability for further barrier passages in the same jump event should depend on the passage of the xenon atom through the central part of the supercage. The values of the absorbingboundary function χ(t) at the end of the simulations are listed in Table 5. The differences between χCEN(35 ps) for the halftrajectories that include a passage of the central part of the sueprcage and χ(35 ps) for all half-trajectories are small. An increase of the probability densities in the central part of the supercage at higher temperatures is to be expected due to the high potential energy in this region. In the period until deactivation, the energy of the xenon atom is sufficiently high to reach the central part of the supercage. Therefore, the temporary occupation of a high-energy position in the center of the supercage has to be seen as a side-effect of the slow energy transfer rather than as an indication of two separate diffusion pathways. The diffusion of xenon in zeolite NaY can be described by a single diffusion mechanism. Thus, the high value of 〈|d˙ (0)|〉, the mean velocity component along the reaction coordinate in the transition state, for trajectories that leave the two-supercage region (see Table 3) does not result from different diffusion mechanisms. The high value of 〈|d˙ (0)|〉 must be a consequence of deterministic movements. The faster the xenon atoms are moving toward a window, the less they can be deflected from this direction. Therefore, the high value of 〈|d˙ (0)|〉 for trajectories that leave the two-supercage region is an indication for correlation effects between different cage-to-cage jumps within the same jump event. VI. Conclusions The cage-to-cage diffusion of xenon in zeolite NaY was investigated by MD simulations of trajectories started in the window plane that forms the boundary between two supercages. Transmission coefficients were obtained from the normalized reactive-flux correlation function κ(t), which reaches its linear asymptotic limit (LAL) after a period of ca. 10-30 ps. In the
temperature range 77-210 K, the transmission coefficient decreases with increasing temperatures, but it does not drop below 0.75. Transition-state theory could, therefore, be used to obtain the order-of-magnitude of the diffusion coefficient. The typical low-temperature recrossing occurs when the kinetic energy is insufficient to surmount an additional potentialenergy barrier due to an instantaneous shift of the barrier maximum along the reaction coordinate. At higher temperatures, the typical recrossing takes place after collisions of the xenon atom with the walls in the new supercage. The latter recrossings are determined by the geometry of the supercage and do not depend on the starting conditions. The recrossings at high temperatures occur due to the slowness of the energy transfer from the xenon atom to the zeolite lattice. The energetic deactivation of the xenon atom takes place by a repeated transfer of energy in collisions with the lattice atoms. The higher the temperature, the longer the period until a trajectory becomes energetically deactivated for cage-to-cage diffusion. Therefore, the slow energy transfer can account for the time scale in which κ(t) approaches its LAL. The decrease in the transmission coefficient with increasing temperature can be explained from the increasing time period for deactivation. An absorbing-boundary function χ(t), representing the percentage of trajectories that has reached the boundary of the twosupercage region around the central transition-state plane, was used to investigate the tendency of xenon atoms to pass through several transition-states during one activation period. This function reaches its LAL on the same time scale as the reactiveflux correlation function. The fact that the LAL value is well below χ(t) ) 1 indicates that the majority of the trajectories becomes deactivated within the boundaries of the two-supercage region. Therefore, the diffusion process can be divided into separate jump events that consist of a limited number of barrier passages during a single activation period. We suggest to consider the occurrence of several barrier passages during one jump event in terms of an absorbing-boundary factor χ that is determined in the same way as the transmission coefficient is determined from the normalized reactive-flux correlation function. Since the diffusion process can be divided into separate jump events, diffusion can be described by a hopping model. This allows to use transition-state theory with dynamical corrections for the calculation of diffusion coefficients. However, a hopping-model approach to the calculation of diffusion coefficients for xenon in zeolite NaY will have to take into account not only the transmission coefficient κ but also the absorbingboundary factor χ. Further work on this topic is in progress. Acknowledgment. We thank the Fonds der Chemischen Ind., Frankfurt, for financial support of this work, particularly for a grant to T.M. References and Notes (1) Chen, N. Y.; Degnan, T. K.; Smith, C. M. Molecular Transport and Reaction in Zeolites; VCH Publishers: New York, 1994. (2) Ka¨rger, J.; Ruthven, D. M. Diffusion in Zeolites and Other Microporous Solids; John Wiley: New York, 1992. (3) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 1987. (4) Demontis, P.; Suffritti, G. B. In: Modelling of Structure and ReactiVity in Zeolites; Catlow, C. R. A., Ed.; Academic Press: London, 1992. (5) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1992, 96, 1051. (6) Demontis, P.; Suffritti, G. B.; Fois, E. S.; Quartieri, S. J. Phys. Chem. 1992, 96, 1482.
Cage-to-Cage Jumps of Xenon in Zeolite NaY (7) Kawano, M.; Vessal, B.; Catlow, C. R. A. J. Chem. Soc., Chem. Commun. 1992, 879. (8) Nicholas, J. B.; Trouw, F. R.; Mertz, J. E.; Iton, L. E.; Hopfinger, A. J. J. Phys. Chem. 1993, 97, 4149. (9) El Alrami, S.; Vigne´-Maeder, F.; Bigot, B. J. Phys. Chem. 1992, 96, 9417. (10) El Alrami, S.; Kolb, M. J. Chem. Phys. 1993, 98, 1509. (11) Santikary, P.; Yashonath, S.; Ananthakrishna, G. J. Phys. Chem. 1992, 96, 10469. (12) Yashonath, S.; Santikary, P. J. Phys. Chem. 1993, 97, 3849. (13) Yashonath, S.; Santikary, P. Mol. Phys. 1993, 78, 1. (14) Schrimpf, G.; Brickmann, J. J. Comput.-Aided Mater. Design 1995, 2, 49. (15) Yashonath, S.; Santikary, P. J. Phys. Chem. 1993, 97, 13778. (16) Fritzsche, S.; Haberlandt, R.; Ka¨rger, J.; Pfeifer, H.; Heinzinger, K. Chem. Phys. 1993, 174, 229. (17) Demontis, P.; Suffritti, G. B. Chem. Phys. Lett. 1994, 223, 355. (18) Bull, L. M.; Henson, N. J.; Cheetham, A. K.; Newsam, J. M.; Heyes, S. J. J. Phys. Chem. 1993, 97, 11776. (19) Eyring, H. J. Chem. Phys. 1935, 3, 107. (20) Ruthven, D. M.; Derrah, I. H. J. Chem. Soc., Faraday Trans. 1 1972, 68, 2332. (21) Ka¨rger, J.; Pfeifer, H.; Haberlandt, R. J. Chem. Soc., Faraday Trans. 1 1980, 76, 1569. (22) Wigner, E. Trans. Faraday Soc. 1938, 34, 29. (23) Miller, W. H. J. Chem. Phys. 1976, 65, 2216. (24) Keck, J. Discuss. Faraday Soc. 1962, 33, 173. (25) Allen, M. P.; Schofield, P. Mol. Phys. 1980, 39, 207. (26) Laria, D.; Ciccotti, G.; Ferrario, M.; Kapral, R. J. Chem. Phys. 1992, 97, 378. (27) Zichi, D. A.; Ciccotti, G.; Hynes, J. T.; Ferrario, M. J. Phys. Chem. 1989, 93, 6261. (28) Ciccotti, G.; Ferrario, M.; Hynes, J. T.; Kapral, R. J. Chem. Phys. 1990, 93, 7137. (29) Keirstaed, W. P.; Wilson, K. R.; Hynes, J. T. J. Chem. Phys. 1991, 95, 5256. (30) Bergsma, J. P.; Gertner, B. J.; Wilson, K. R.; Hynes, J. T. J. Chem. Phys. 1987, 86, 1356. (31) Gertner, B. J.; Wilson, K. R.; Hynes, J. T. J. Chem. Phys. 1989, 90, 3537. (32) Rosenberg, R. O.; Berne, B. J.; Chandler, D. Chem. Phys. Lett. 1980, 75, 162. (33) Wilson, M. A.; Chandler, D. Chem. Phys. 1990, 149, 11. (34) Depaepe, J. M.; Ryckaert, J. P.; Bellemans, A. Mol. Phys. 1993, 78, 1575.
J. Phys. Chem., Vol. 100, No. 11, 1996 4581 (35) Becker, K. D.; Hoheisel, C. J. Chem. Phys. 1982, 77, 5108. (36) Verkhivker, G.; Elber, R.; Gibson, Q. H. J. Am. Chem. Soc. 1992, 114, 7866. (37) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1991, 95, 8866. (38) Heink, W.; Ka¨rger, J.; Pfeifer, H.; Stallmach, F. J. Am. Chem. Soc. 1990, 112, 2175. (39) Ka¨rger, J.; Pfeifer, H.; Stallmach, F.; Feoktistova, N. N.; Zhdanov, S. P. Zeolites 1993, 13, 50. (40) Fitch, A. N.; Jobic, H.; Renouprez, A. J. Phys. Chem. 1986, 90, 1311. (41) Schrimpf, G.; Schlenkrich, M.; Brickmann, J.; Bopp, P. J. Phys. Chem. 1992, 96, 7404. (42) Fritzsche, S.; Haberlandt, R.; Kaerger, J.; Pfeifer, H.; Wolfsberg, M. Chem. Phys. Lett. 1990, 171, 109. (43) Demontis, P.; Suffritti, G. B.; Quartieri, S.; Fois, E. S.; Gamba, A. Zeolites 1987, 7, 522. (44) Demontis, P.; Suffritti, G. B.; Quartieri, S.; Fois, E. S.; Gamba, A. J. Phys. Chem. 1988, 92, 867. (45) Kiselev, A. V.; Du, P. Q. J. Chem. Soc., Faraday Trans. 2 1981, 77, 1. (46) Ewald, P. P. Ann. Phys. 1921, 64, 253. (47) Dufner, H.; Schlenkrich, M.; Brickmann, J. J. Comput. Chem., submitted for publication. (48) Bezus, A. G.; Kiselev, A. V.; Lopatkin, A. A.; Du, P. Q. J. Chem. Soc., Faraday Trans. 2 1978, 74, 367. (49) A Ä ngya´n, J. G.; Colonna-Cesari, F.; Tapia, O. Chem. Phys. Lett. 1990, 166, 180. (50) Northrup, S. H.; Hynes, J. T. J. Chem. Phys. 1980, 73, 2700. (51) Chandler, D. J. Chem. Phys. 1978, 68, 2959. (52) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (53) Berne, B. J.; Borkovec, M.; Straub, J. E. J. Phys. Chem. 1988, 92, 3711. (54) Voter, A. F.; Doll, J. D. J. Chem. Phys. 1985, 82, 80. (55) Carter, E. A.; Ciccotti, G.; Hynes, J. T.; Kapral, R. Chem. Phys. Lett. 1989, 156, 472. (56) Verlet, L. Phys. ReV. 1967, 159, 98. (57) Kast, S. M.; Nicklas, K.; Ba¨r, H. J.; Brickmann, J. J. Chem. Phys. 1994, 100, 566. (58) Straub, J. E.; Berne, B. J. J. Chem. Phys. 1985, 83, 1138. (59) Mosell, T. Dissertation, Technische Hochschule Darmstadt 1995.
JP952644C
