Xenon Diffusion in Zeolite NaY: Transition-State Theory with

Zeolite Adsorption Free Energies from ab Initio Potentials of Mean Force. Hui Li , Christopher ... High-Accuracy Estimation of 'Slow' Molecular Diffus...
6 downloads 6 Views 549KB Size
Download PDF
4582

J. Phys. Chem. 1996, 100, 4582-4590

Xenon Diffusion in Zeolite NaY: Transition-State Theory with Dynamical Corrections Thomas Mosell,† Gerhard Schrimpf,†,§ 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 coefficient for xenon in zeolite NaY at infinite dilution is calculated at low temperatures within a hopping model based on cage-to-cage migration only. Diffusion is modeled as a sequence of jump events that may consist of several barrier passages. The number of jump events is calculated from transition-state theory using the potential of mean force as a basis. The potential of mean force is more or less independent of temperature. In the conversion of the jump rate to the diffusion coefficient, dynamical corrections are taken into account. The mean number of barrier passages per jump event increases significantly when the temperature is raised. However, the contribution of the dynamical corrections to the activation energy remains small. In the range from 140-210 K, the diffusion coefficients obtained from the hopping-model are in excellent agreement with corresponding data from conventional MD simulations.

I. Introduction The existence of specific adsorption sites for a guest molecule in a zeolite is related to the occurrence of a possibly considerable activation energy for the diffusion of that guest molecule within the zeolite. Diffusion in the regular micropores can often be described as a sequence of jumps between adsorption sites. At high activation energies, jump events can become so rare that diffusion cannot be observed on the time scale of conventional molecular dynamics (MD) simulations, which are presently limited to a few nanoseconds for typical simulation systems. In recent years, the combination of transition-state theory and the calculation of dynamical corrections has evolved as a promising approach for the simulation of systems governed by rare events.1-3 The crucial rare event is investigated in detail, and conclusions can be based on the dynamics of this process, whereas most parts of the system are considered in the form of statistical data only. In dynamically corrected transition-state theory, the rate constant k can always be written in the form

k ) κkTST

(1)

where kTST is the rate constant from transition-state theory and κ is the transmission coefficient for trajectories that pass through the transition state. The idea to calculate absolute rate constants from probability densities by transition-state theory4,5 has been known and used since the 1930s. Keck6 introduced MD simulations of trajectories that originate in the transition state in order to obtain dynamic corrections to statistical rate constants. The transmission coefficient can be determined either from an analysis of the final configurations of these trajectories6,7 or from the reactive-flux correlation function.8 Voter and Doll9 generalized this approach to free-energy hypersurfaces with multiple minima. They introduced dynamic correction factors that have to be calculated for the different combinations of minima, resulting in a multistate analogue of the transmission coefficient to that of the two-state model. †

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. * To whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, February 15, 1996. ‡ §

0022-3654/96/20100-4582$12.00/0

Transition-state theory has repeatedly been applied to the diffusion of guest molecules in zeolites.10,11 On the basis of these theories, the macroscopic diffusion coefficients can be obtained from a hopping model. Dynamical corrections were first taken into account by June et al.,12 who based their calculations on the multi-state approach of Voter and Doll.9 These authors obtained macroscopic diffusion coefficients at finite loadings by performing Monte Carlo simulations on a lattice model with jump rates obtained for the different combinations of minima. At infinite dilution, a different approach can be used. Diffusion at infinite dilution in a zeolite is dynamically equivalent to the diffusion of one particle on a regular lattice, because no sorbate-sorbate interactions have to be considered. The diffusion coefficient can be calculated with a hopping model that is based on the geometric arrangement of the lattice knots and the transition probability for a jump.13 In this work the diffusion will be described in terms of jump events that consist of a sequence of barrier passings without an intermediate energetic deactivation of the trajectory. The multistate approach of Voter and Doll is not necessary, if all minima are equivalent and the different jump lengths are taken into account by an additional correction factor. In a preceding publication14 (referred to as paper I hereafter), we introduced an absorbing boundary that surrounds the region of the two minima on either side of the central transition state and defined an absorbingboundary factor χ as the percentage of trajectories that is absorbed by this boundary during a single jump event. This absorbing-boundary factor can be used to determine the mean number nBP of barrier passages per jump event. In the present work, a hopping model is presented that takes into account both a transmission coefficient and an absorbingboundary factor. This hopping model considers jumps between different supercages of zeolites with faujasite structure. The hopping model is applied to the diffusion of xenon in zeolite NaY at low temperatures at which the diffusion can be described as a sequence of jump events. Transmission coefficients for cage-to-cage jumps of xenon in zeolite NaY and the corresponding absorbing-boundary factors have been presented in paper I14 for several temperatures from 77 to 210 K. In this work, an approach of Ciccotti et al.15 is used to obtain the potential of mean force along the coordinate perpendicular to the transition state. The probability density for the transition © 1996 American Chemical Society

Xenon Diffusion in Zeolite NaY

J. Phys. Chem., Vol. 100, No. 11, 1996 4583

Figure 2. Unit cell of a diamond lattice. The length ∆r is the distance between two lattice points, φi the angle between the direction of the jump and one of the axes of the unit cell, and a is the length of the unit cell.

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.

state and, consecutively, the transition-state theory value for the jump rate are determined from the potential of mean force. The diffusion coefficient for xenon in zeolite NaY is calculated by applying the hopping model. The results are compared to those of conventional MD simulations with the same model system.16 The paper is organized as follows: The structure and the force-field representation of zeolite NaY are shortly reviewed in section II. In section III, the hopping model for the calculation of the diffusion coefficient is presented. Section IV is concerned with the determination of the potential of mean force. The results of the simulations are reported in section V. Section VI contains the results from calculation of the diffusion coefficients and a discussion of these values. Concluding remarks are summarized in section VII.

to the polarization of the xenon atom is taken into acount. For the calculation of the induction energy, charges are assumed on the sodium and oxygen atoms. The Lennard-Jones parameters are those of Kiselev and Du20 and the polarizability of xenon was taken from Bezus et al.21 III. Hopping Model for the Diffusion Coefficient In previous simulations,16,22 only one specific adsorption site, which is in front of the 4-T ring of the sodalite cage, has been observed for xenon in zeolite NaY. Cage-to-cage jumps were found to be the rate-determining step of the diffusion process in these simulations. If cage-to-cage jumps are the ratedetermining event, the supercage area can be considered as one joint valley in the hopping model. The hopping model can, consequently, be based on the diamond lattice formed by the centers of the supercages. The diffusion of a single particle on a void lattice is equivalent to the diffusion of a vacancy in a crystal.13 For a cubic system, the diffusion coefficient D is given by23 z

D ) 1/2Γ0∆r2∑cos2 φi

(2)

i)1

II. Model Representation of Zeolite NaY Zeolite NaY (see Figure 1) has a cubic unit cell with a cell parameter a ) 24.85 Å.17 The zeolite can be synthesized with varying contents of aluminum and counterions for charge compensation. In our model, the number of Na atoms corresponds to a Si/Al ratio of 3:1. The force field does, however, not distinguish between silicon and aluminum atoms in the tetrahedral TO4 sites. Zeolite NaY contains two kinds of cages. The cubooctaedral sodalite cages are inaccessible to xenon. Linkage of the sodalite cages by hexagonal prisms leads to the formation of a three-dimensional net of large supercages, which are accessible for xenon through windows framed by rings of 12 TO4 units. The centers of the supercages form a diamond lattice. In the model system, the cation sites SI und SII are completely occupied by sodium atoms. The force field for xenon and zeolite NaY is described in detail in paper I.14 Only a short summary is given here. The force field for the fully flexible zeolite lattice is based on the force-field representation for zeolite A by Demontis et al.18 Only the Na-Na and the Na-O interactions had to be modified by Schrimpf et al.19 in order to adapt the force field to zeolite NaY. The interactions between xenon and the lattice atoms are as well taken from Schrimpf et al.16,19 Not only the contributions from Lennard-Jones potentials but also the induction energy due

where Γ0 is the jump rate, ∆r the distance between two lattice points, φi the angle between the direction of the jump and one of the axes of the unit cell, and z the number of next neighbors on the lattice. The diffusion coefficient for the diamond lattice23

D ) 1/8Γ0a2

(3)

results for z ) 4, φi ) 54,7°, and ∆r ) 0,433 a, where a is the length of the unit cell (see Figure 2). III.1. Dynamical Corrections. The hopping model is based on the assumption that the probability γ for the continuation of a jump event is independent from the previous number of barrier passages in the jump event. It will be assumed that the guest particle leaves the supercage through each of the windows with the same probability. As a result, a jump event that consists of n barrier passages is equal to n jumps over one barrier. Thus, the right-hand side of eq 3 has to be multiplied by the mean number nBP of barrier passages per jump event in order to obtain the dynamically corrected diffusion coefficient:

D ) 1/8nBPΓ0a2

(4)

The dynamic corrections for the hopping model are suitably expressed in terms of two factors that take into account the

4584 J. Phys. Chem., Vol. 100, No. 11, 1996

Mosell et al.

effects during a single barrier passage and the probability γ for the continuation of the jump event separately. Since the continuation of the jump event may involve a second passage over the same barrier, the transmission coefficient κ determined from trajectories started in the transition state cannot be used directly to describe the former of the two effects. A “local” transmission coefficient κ′ that considers recrossings during a single barrier passage only has to be introduced to convert the jump rate kTST obtained by transition-state theory to the dynamically corrected jump rate Γ0:

Γ0 ) κ′kTST

(5)

different barrier passages results in: n

Ωn(II) ) ∑(4n-i - 1)(4i-1 - 1) ) n(4n-1 + 1) i)1

2

/3(4n - 1) (10)

n

Ωn(I) ) ∑(4n-i - 1 + 4i-1 - 1) ) 2/3(4n - 1) - 2n (11) i)1

The total number of differents path with n barrier passages is

Ωn ) n4n-1

resulting in

D ) 1/8κ′nBPkTSTa2

(6)

Thus, the deviations from the transition-state theory value of the diffusion coefficient are included in the product of the mean number nBP of barrier passages per jump event multiplied by the local transmission coefficient κ′. The difference between the local transmission coefficient κ′ and the transmission coefficient κ of the MD simulations results from the trajectories that leave the supercage through the window of entry. The probability for this to happen is γ/4 at each stage of the jump event: ∞

( )

κ ) κ′ + ∑ -

γ

n)1

4

) κ′ -

γ 4+γ

Pn ) nγn-1(1 - γ)2

(7)

)

∞ Ωn(II) 1 Ωn(I) + ‚ Pn κ′χ ) ∑ 1‚ Ωn 2 Ωn n)1

(9)

where Pn is the probability for a jump event with n barrier passages, Ωn is the total number of jump events with n barrier passages, and Ωn(I) and Ωn(II) are the numbers of jump events where the two supercages are left at one and two ends of the jump event, respectively. For barrier passage i, there are 4i-1 different paths in one direction and 4n-i in the other one. Only one path in each direction remains within the two supercages, whereas (4i-1 1) (4n-i - 1) different paths contribute to Ωn(II) and (4i-1 - 1) + (4n-i - 1) different paths to Ωn(I). The sum over the n

(14)

and consequently

4 4 (1 - γ) 3γ κ′χ ) 1 - (1 - γ) + ) 3 3 4-γ 4-γ 2

(15)

Finally, the local transmission coefficient can be substituted by eq 7:

γ 3γ ) κ′χ ) χ κ + 4-γ 4+γ

(

(8)

It remains to determine γ as a function of the transmission coefficient κ and the absorbing-boundary factor χ. III.2. Probability for the Continuation of a Jump Event. The part of jump events that extend further than the two supercages on either side of the transition state is equal to the product of the absorbing-boundary factor χ multiplied by the local transmission coefficient κ′. These jump events can extend into other supercages either at one end or at both ends. The product κ′χ may be determined from

)

(

n)1

(

(13)

∞ 4 1 κ′χ ) ∑ n - + γn-1(1 - γ)2 n-1 3 3‚4 n)1

∞

˜ n ) ∑nγn-1(1 - γ) ) 1/(1 - γ) nBP ) ∑nP n)1

The probability Pn to observe a jump event with n barrier passages in the simulation procedure does not correspond to a standard geometric series. Instead, it has to be taken into account that the criterion for the breakoff of the jump event has to be fulfilled in both directions. Moreover, a jump event with n barrier passages is accessible from n phase points in the transition-state plane and will therefore be sampled with an n-fold probability. Therefore, Pn has to be substituted by

and eq 9 rearranges to n

Jump events that intermediately extend into other supercages are not to be considered here, because they will be taken into account by the absorbing-boundary factor. The mean number nBP of barrier passages per jump event can also be expressed in terms of γ. The probability P ˜ n for a jump event with n barrier passages is given by the geometric distribution γn-1(1 - γ), because the jump event is continued n - 1 times, before it is ended with the probability (1 - γ): ∞

(12)

)

(16)

which can be rearranged to a quadratic equation for γ. Only the positive solution is useful, because γ g 0 must be valid:

γ)2

x(3 - χ)2 + 4κχ(3 + (1 + κ)χ) - (3 - χ) 3 + (1 + κ)χ

(17)

The probability γ for the continuation of the jump event can now be used to calculate the mean number nBP of barrier passages per jump event by eq 8 and the local transmission coefficient κ′ by eq 7. IV. Simulation Procedure IV.1. Probability Density in the Transition State. The aim of the simulations is the calculation of the transition-state theory value of the rate constant kTST, which can be determined from the probability density W(rq) for the transition state at rq:24

kTST )

x

kBT W(rq) 2πµ

(18)

where µ is the reduced mass for a vibration along the reaction coordinate r.

Xenon Diffusion in Zeolite NaY

J. Phys. Chem., Vol. 100, No. 11, 1996 4585

The probability density in the transition state can be determined from MD simulations, if a special approach that makes possible a sampling of the energetically unfavorable parts of configurational space is applied in the simulations. The umbrella-sampling technique25 uses biased potentials to obtain a modified probability distribution. Alternatively, free energy simulations26 where a suitably defined coordinate is varied systematically to reach unfavorable regions of configurational space can be performed. The present work is based on an approach by Ciccotti et al.,15 who constrained their system at different points on the reaction coordinate. From these simulations, which have been termed constrained reaction coordinate dynamics (CRCD), the mean value 〈FRC(r)〉 for the force exerted on the reaction coordinate is determined for each point. The potential of mean force ∆A(r) along the reaction coordinate is then obtained by integration of these forces. In the integration, a Jacobi determinant has to be taken into account due to a coordinate transformation between the general coordinates of the reaction coordinate and the Cartesian coordinate system.24 However, if the reaction coordinate is defined in a way that the contributions of the Jacobi determinant is equal to 1, the integration can be performed by

∆A(r) ) -∫0 〈FRC(r′)〉 dr′ r

(19)

The probability density is subsequently determined from the potential of mean force ∆A(r) by

W(ri) ) W0 e-∆A(ri)/kBT

(20)

where W0 is a normalization constant. IV.2. Definition of the Reaction Coordinate. The reaction coordinate and the transition state for cage-to-cage diffusion are the same as in paper I.14 The reaction coordinate d is parallel to the body diagonal [1 1 1] of the unit cell. Note that this reaction coordinate is not defined by a specific minimum-energy path. Instead the definition takes advantage of the symmetry of the system. Moreover, this definition in Cartesian coordinates forces that the Jacobi determinant to be equal to 1, which makes the calculation of the potential of mean force much easier. Each point on the reaction coordinate corresponds to a plane in three-dimensional space:

n‚(p - a) ) d

(21)

where p is the position of the guest particle, n the normal vector of the plane, and a a point in the window that defines the zeropoint of the reaction coordinate. Thus, d equals the distance between the xenon atom and the window plane. The results for coordinate point d can be transferred to -d due to the symmetry of the system. IV.3. Details of the Calculations. During the CRCD simulations the xenon atom was confined to a plane by a constraint. Simulations were performed in the range from d ) 0 Å to d ) 9.75 Å for points on the reaction coordinate distanced by 0.25 Å, covering the space contained in one supercage. Simulations for different points on the reaction coordinate were started consecutively. At the beginning of each simulation, the xenon atom was shifted to the new plane by a move parallel to the normal vector n. The results for d ) 0 Å were obtained from the 1.1 ns simulation described in paper I.14 For all other points, a 20 ps period was used for the equilibration of the local distortions introduced by the shift, before a 40 ps simulation was performed in order to sample the mean force along the reaction coordinate.

In some of the planes the guest particle can in principle leave the supercage through one of the other windows. An additional potential VADD(t) was introduced to prevent the particles from permanently leaving the region of space that has to be sampled in the simulations:

VADD(t) ) V0,ADD(rWIN(t) - rADD)2

rWIN(t) > rADD (22)

where rWIN is the appropriate coordinate for the passage through the window in question. The definition of rWIN is analogous to that for the reaction coordinate. The additional potential does not influence the probability densities within the supercage, because it is applied only to configurations where the xenon atom has left the supercage by more than rADD ) 0.2 Å. A value of 0.5 J/m2 was used for V0,ADD. Only configurations where the xenon atom is located within the supercage are taken into account in the calculation of the mean forces. The mean forces were integrated to the potential of mean force by interpolation with a cubic spline. In the calculation of kTST the reduced mass µ was approximated by the mass of the guest particle.12 The constant W0 of eq 20 normalizes the probability to find a particle within a single supercage to 1. The simulations were performed with periodic boundary conditions at constant volume. The Verlet algorithm27 with a time step of ∆t ) 1 fs was used for the integration of the motion. The stochastic-thermalization method of Kast et al.28 was applied to obtain the canonical ensemble. The temperatures studied are 77, 100, 140, 160, 180, and 210 K. All simulations were performed with a MD simulation code developed in our group. V. Results of the CRCD Simulations V.1. Potential of Mean Force. In Figure 3, the mean force 〈FRC(d)〉 along the reaction coordinate d is depicted as a function of this coordinate for the different temperatures. A positive value of 〈FRC(d)〉 indicates that the force vector is predominantly directed towards higher values on the reaction coordinate. In Figure 3a, only the results for 210 K are displayed, because the temperature dependence of the results is practically invisible on the scale of the plot. In the simulations for d g 9 Å, the wall of the supercage has been reached. The results for d < 9 Å are displayed separately in Figure 3b to show the differences between the results at different temperatures. Differences between the mean forces at different temperatures are observed mainly in some distinct regions. The potential of mean force ∆A(d), displayed in Figure 4, is similar for all temperatures. The first minimum at d ≈ 2.7 Å is chosen as the zero point of the potential of mean force to allow for a better comparison of the results. This minimum, as well as the second one at d ≈ 8.2 Å, is caused by xenon atoms adsorbed in front of the 4-T rings of the sodalite cages. Each supercage contains six adsorption sites, but they occur at only two distinct positions on the reaction coordinate due to the 3-fold symmetry of the system. The location of the adsorption site for xenon is in accordance with the results of the simulations for the transmission coefficient14 and with those of other simulations.16,22 The second minimum exhibits a shoulder around 6.5 Å, which is created by xenon atoms close to a 4-T ring of the hexagonal prism. The high mobility along the sequence of 4-T rings has also been reported before.16 A similar shoulder remains invisible for the first minimum, because the movements in front of the sequence of 4-T rings are nearly perpendicular to the reaction coordinate for that adsorption site. The probability densities for the whole sequence of 4-T rings contribute to the value of

4586 J. Phys. Chem., Vol. 100, No. 11, 1996

Figure 3. Component of the mean force 〈FRC(d)〉 along the reaction coordinate d as a function of this coordinate. (a, top) The results for 210 K are displayed in the range from 0-9.75 Å. (b, bottom) The results at different temperatures are displayed for d < 9 Å. A positive value of 〈FRC(d)〉 indicates that the force vector is predominantly directed toward higher values on the reaction coordinate.

the potential of mean force for the same points on the reaction coordinate. As a result, the first minimum is slightly lower than the second.

Mosell et al. A maximum of the potential of mean force is located in the window at d ) 0 Å. This maximum corresponds to a minimum of the probability density. In transition-state theory with dynamical corrections, the product of kTST multiplied by the transmission coefficient κ is constant.29 Therefore, a low probability density in the transition state is connected with a high transmission coefficient. As a result, only a small number of trajectories is necessary to get a good estimate of the transmission coefficient.29 Therefore, the choice d ) 0 Å for the transition state is not only preferable due to the symmetry of the system but also the best choice for a transition-state plane by energetic reasons. At very low temperatures, the barrier of the potential of mean force at d ) 0 Å increases slightly with decreasing temperature. The barrier at 77 K is 9.7 kJ/mol, in comparison to results between 8.7 and 9.3 kJ/mol at higher temperatures. Moreover, the barrier 0 K must be equal to the potential-energy difference of 10.5 kJ/mol16 between the most favorable configuration at the adsorption site and in the window plane. At temperatures above 100 K the differences in the barrier height seem to be more or less statistical and the dependence of the potential of mean force on temperature must at least be significantly less strong. On the basis of the present results, it is impossible to decide whether there is any systematic temperature dependence. A tendency for a decrease in the barrier of the potential of mean force with increasing temperatures is observed for the second maximum at 5.3 Å. The difference in potential energy between the position at an adsorption site and the maximum of the barrier between two adsorption sites in the same supercage is 7.1 kJ/mol.16 When the temperature is raised from 77 to 210 K, the barrier of the potential of mean force is reduced from 6.3 to 5.2 kJ/mol. The barrier for jumps within a supercage is significantly smaller than that for jumps between different supercages. Therefore, changes between adsorption sites within a supercage should occur more frequently than cage-to-cage jumps. This supports the applicability of a hopping model that treats the supercage as a single area in reaction space. Probability densities for the plane at d ) 0 Å were calculated from the potential of mean force using eq 20. The probability densities for the different temperatures are listed in Table 1

Figure 4. Potential of mean force A(d) for the reaction coordinate d at different temperatures.

Xenon Diffusion in Zeolite NaY

J. Phys. Chem., Vol. 100, No. 11, 1996 4587

Figure 5. Sum function of the probabilities for the reaction coordinate d at different temperatures.

TABLE 1: Probability Density W(0 Å) in the Transition State and Jump Rate kTST Obtained from Transition-State Theory T [K]

W(0) [Å-1]

kTST [s-1]

77 100 140 160 180 210

2.1 × 1.6 × 10-5 3.7 × 10-4 5.3 × 10-4 1.4 × 10-3 2.2 × 10-3

5.8 × 104 4.9 × 106 1.4 × 108 2.1 × 108 6.0 × 108 1.0 × 109

10-7

together with the jump rates kTST calculated from these probability densities. A discussion of these values will be presented below. V.2. Sum of Probabilities. In a region from 5.5 to 6.25 Å, there are differences in the mean forces 〈FRC(d)〉 at different temperatures that cannot be attributed to an unambiguous temperature in this region. There is no preferred site for xenon at these coordinate values. If the component of the force along the reaction coordinate depends strongly on the position in the plane, the mean force will not properly be determined from simulations that do not reach an equilibrium distribution between the different local parts of configurational space. Moreover, the three other windows of the supercage are located in the region from 5.5 to 8.0 Å. Xenon atoms have the possibility to leave the supercage through one of these windows. Despite the introduction of the additional potential VADD, the time during which the xenon atom was present in the supercage was reduced in some of the simulations from 40 to up to 28 ps. The deviations in 〈FRC(d)〉 may, therefore, be attributed to a greater statistical uncertainty in this region. In principle, it would be possible to normalize the probability densities to the first minimum in the supercage, because all adsorption sites are equivalent. The sum function of the probabilities as a function of the reaction coordinate in Figure 5 shows, however, that there is no unambiguous boundary between the two minima of the potential of mean force at the highest temperatures investigated in this work. Therefore, the probability density in the transition state was normalized on the probabilities for the whole supercage. The effect of the statistical uncertainty on the barrier at 5.3 Å may be estimated from the sum of the probabilities. Since

the two minima have to be equally populated, the sum functions should intersect ideally at a value of 0.5 at one point that forms the boundary between the two regions. It is seen from Figure 5 that the deviations from the ideal result are within a range of (0.1. In a series of six independent numerical experiments on the potential of mean force at 180 K, the deviations from 0.5 are never larger than (0.15 in any case. Two of these simulation series at 180 K were started at d ) 10 Å and followed in the reverse direction. There are no indications that there might be a systematic dependence of the mean forces on the direction of the shifts on the reaction coordinate. The order of magnitude of the deviations for the barrier at 5.3 Å should be transferable to the barrier at d ) 0 Å. Considering the simplicity of the force field applied in the present calculations, we believe that a prolongation of the simulation to reduce the statistical uncertainty would be inapproriate. V.3. Locating the Xenon Atoms. Contour maps for the probabilities to find a xenon atom in different locations on the plane are depicted in Figure 6 for the planes d ) 0 Å and d ) 2.75 Å at 180 K. The perpendicular coordinates q1 and q2 on the plane are pointing along [1 1-2] and [1-1 0], respectively. The origin of the coordinate system q1q2 is located in the point of intersection between the plane and the body diagonal [1 1 1]. For d ) 0 Å, Coordinate q2 is pointing from the origin to one of the six oxygen atoms of the 12-T ring that are located in the window plane. The most probable position in the window plane (see Figure 6a) is located in front of a T atom site. In this position, there are strong attractive interactions between the xenon atom and three oxygen atoms, which are all part of the 12-T ring. However, the position in front of one of the six oxygen atoms in the window plane is only slightly less favorable. The difference in the probabilities corresponds to a (free) energy difference of 0.5 kJ/mol. In contrast, the probability of finding a xenon atom near the center of the window is very low. The xenon atom, therefore, remains near the walls of the zeolite during the passage of the transition state. The plane at d ) 2.75 Å is typical for the planes close to the adsorption site. High probability values are restricted to a very small space, although the total void space is larger than that in

4588 J. Phys. Chem., Vol. 100, No. 11, 1996

Mosell et al. TABLE 2: Dynamic Corrections for the Cage-to-Cage Diffusion of Xenon in Zeolite NaYa T[K]

k

c

g

κ′

nBP

77 100 140 160 180 210

0.91 0.89 0.87 0.89 0.82 0.77

0.04 0.10 0.15 0.19 0.22 0.34

0.04 0.11 0.17 0.23 0.24 0.35

0.92 0.92 0.91 0.94 0.88 0.85

1.04 1.13 1.21 1.29 1.32 1.53

a Transmission coefficient κ, absorbing-boundary factor χ, probability γ for the continuation of a jump event, local transmission coefficient κ′ and mean number nBP of barrier passages per jump event.

TABLE 3: Hopping-Model Results for the Diffusion Coefficient for Xenon in Zeolite NaYa T [K]

D [cm2 s-1]

κ′χb

77 100 140 160 180 210

4.3 × 10-10 3.9 × 10-8 1.2 × 10-6 2.0 × 10-6 5.3 × 10-6 1.0 × 10-5

0.96 1.03 1.10 1.22 1.16 1.30

a Diffusion coefficient D at infinite dilution. b Deviation from transition-state theory.

Figure 6. Contour maps of the relative probabilities for the location of xenon atoms on the planes d ) 0 Å in the window (a) and d ) 2.75 Å near the adsorption site (b) at 180 K. The definition of the coordinate system q1q2 is explained in the text. Both plots have been drawn on the same scale.

the window. In the potential of mean force ∆A(d), the differences in the width of the area of high probability represent an entropic effect, e.g., the window appears to be a region of high entropy. The temperature dependence of the entropic contributions will usually be different for different points on the reaction coordinate. These differences can explain the decrease in the temperature dependence of the barriers of potential of mean force ∆A(d). VI. Diffusion Coefficient and Activation Energy VI.1. Diffusion Coefficient. The dynamical corrections for the diffusion coefficient are listed in Table 2. The values of the transmission coefficient κ and the absorbing-boundary factor χ are taken from paper I.14 The probability γ for the continu-

ation of the jump event increases signifcantly when the temperature is raised. The differences between γ and the absorbing-boundary factor χ are very small, although a ratio of about γ/χ ≈ 4/3 would be expected at a first guess, because only three of the four windows of the supercage contribute to χ. The similarity of both parameters is a result of the reduction of this ratio due to the consideration of the local transmission coefficient κ′ and the influence of jump events with several barrier passages. The increase in γ induces a significant increase in the mean number nBP of barrier passages per jump event. Whereas at 77 K nearly all jump events consist of one barrier passage, the mean number nBP of barrier passages increases to 1.53 passages at 210 K. Therefore, the slow energetic deactivation of the trajectories at higher temperatures leads to a significant increase in the diffusion coefficient. There is only a small temperature dependence of the local transmission coefficient κ′. In the simulations at up to 160 K the results are identical within the error boundaries. At higher temperatures, the local transmission coefficient κ′ decreases moderately. The temperature dependence is smaller than that of the transmission coefficient κ, because κ is reduced by the trajectories with several barrier passages within a jump event. However, the remaining temperature dependence of κ′ indicates that the temperature depencence of κ may not simply be explained by the occurrence of jump events with several barrier passages. The results for the diffusion coefficient D determined by the hopping model are listed in Table 3. The diffusion coefficient increases from 4.3 × 10-10 cm2/s at 77 K to 1.0 × 10-5 cm2/s at 210 K. For xenon in zeolite NaY, neither κ′ nor nBP deviate strongly from 1. Therefore, the classic transition-state theory without dynamical corrections is a good approximation for the diffusion coefficient of xenon in zeolite NaY at very low temperatures. However, jump events with several barrier passages become important at higher temperatures and transition-state theory underestimates the diffusion coefficient by ca. 30% at 210 K. VI.2. Activation Energy. The temperature dependence of the diffusion coefficient can often be described by an Arrhenius activation energy EA:

Xenon Diffusion in Zeolite NaY

J. Phys. Chem., Vol. 100, No. 11, 1996 4589

Figure 7. Arrhenius plot for the diffusion coefficient of xenon in zeolite NaY. The horizontal bars indicate results obtained with the hopping model, whereas the diamonds indicate the results of conventional MD simulations of ref 16 at a loading of 1 xenon/supercage. The straight line was obtained from the latter results in the temperature range 150700 K.

∂ ln D ) -EA/R ∂(1/T)

(23)

The Arrhenius plot for the diffusion coefficients of xenon in zeolite NaY is depicted in Figure 7. The diffusion coefficients obtained by Schrimpf und Brickmann16 in conventional MD simulations with a loading of 1 xenon atom/supercage are also displayed in the plot. In addition, the results of these conventional MD simulations in the range 150-700 K are indicated by a straight line. In the temperature range 140-210 K, the agreement between the hopping-model results and the conventional MD simulations is very good. An activation energy of 7.84 ( 0.85 kJ/mol for the diffusion coefficients of the hopping model was obtained by linear regression, whereas the conventional MD simulations in the range 150-700 K resulted in an activation energy of 7.81 ( 0.40 kJ/mol.16 The prefactor D0 was (9 ( 4) × 10-4 cm2/s in the hopping model in comparison to D0 ) (6.6 ( 1.6) × 10-4 cm2/s in the conventional MD simluations. The difference in D0 may result from the different loadings. The hopping model assumes infinite dilution, whereas a part of the attempted barrier passages will be unsuccessful in the conventional MD simulation at a finite loading due to collisions between xenon atoms. The difference between the two values for D0 is, however, statistically not significant. The dependence of the diffusion coefficient on the loading becomes stronger at higher loadings.16 Therefore, the experimental results of Ka¨rger et al.30 for xenon in zeolite NaX at a loading of 3-4 xenon atoms/supercage cannot be compared directly to the hopping-model results. However, the activation energy of (6 ( 1.5) kJ/mol obtained by 129Xe PFG NMR is in good agreement with the results of the conventional MD simulations.16 The results of the hopping model at 100 and 77 K are lower than those expected from an extrapolation of the diffusion coefficients at higher temperatures. Thus, the diffusion of xenon in zeolite NaY at very low temperatures cannot be described by the same temperature-independent activation energy as at higher temperatures. The height of the barrier for the potential of mean force (see Figure 4) at 0 Å depends on temperature at very low temperatures, whereas it is nearly constant at higher temperatures. This microscopic difference must result in a temperature dependence of the activation energy. A tempera-

ture-dependent activation energy for the diffusion of xenon in zeolite NaY has already been reported by Yashonath and Santikary.31 The influence of the dynamical corrections on the activation energy is relatively small. A separate Arrhenius plot for the temperature dependence of the product κ′ nBP results in a contribution of 0.3 kJ/mol to the activation energy in the range 77-210 K. The potential of mean force ∆A(d) is of much greater importance for the Arrhenius activation energy for the diffusion of xenon in zeolite NaY. VI.3. Discussion of the Correlation Effects. The omission of correlation effects between different jump events can be justified on the basis of the adsorption-site changes within a supercage. These changes are more frequent than the cage-tocage jumps, and no memory effect should remain until the next jump event occurs. On the contrary, barrier passages within a jump event may be correlated. For example, the local conditions in the transition state can be correlated, if xenon leaves the supercage through the window of entry. An investigation of the dependence of recrossings on the local conditions in the transition state has already been discussed in paper I.14 In that investigation, a dependence on the initial conditions was only observed for trajectories that remain in the window region before the recrossing occurs. These trajectories are taken into account by the local transmission coefficient κ′ and do not count as successive barrier passages in the hopping model. Nevertheless, a correlation effect will occur, if xenon preferentially leaves the supercage in a successive barrier passage either through the window of entry or through one of the other windows. Xenon atoms that cross the supercage to leave through one of the other windows directly should induce a decrease of the probability to leave the supercage through the window of entry below γ/4. As a result, the local transmission coefficient κ′ determined by eq 7 would be overestimated. At higher temperatures, the absorbing-boundary factor χ becomes more important, and accordingly, the overestimation of the local transmission coefficient κ′ should increase. Since the importance of trajectories that remain in the window region before the recrossing occurs is negligible at higher temperatures, an artificial increase of the local transmission coefficient κ′ should result. However, κ′ decreases with increasing temperatures. Therefore, a return through the window of entry is slightly favored. This must be attributed to the effects of collisions with the walls. Due to the weak correlations between successive barrier passages, the “local transmission coefficient” κ′ can no more be exclusively linked to local effects. However, the fact that κ′ is influenced by the correlation effects avoids a full repercussion of the correlation effects on the diffusion coefficient. VII. Summary and Conclusions The diffusion coefficient for xenon in zeolite NaY at infinite dilution has been treated at low temperatures by a hopping model in which the supercage is a joint minimum and where only cage-to-cage migration is considered. Diffusion is modeled as a sequence of jump events that may consist of several barrier passages in a row. It is assumed that the probability for the continuation of a jump event is independent from the previous number of barrier passages. The number of jump events was calculated from transition-state theory. In the conversion of the jump rate to the diffusion coefficient, dynamical corrections are taken into account. The transition-state theory value of the jump rate was obtained from the potential of mean force for a reaction coordinate that

4590 J. Phys. Chem., Vol. 100, No. 11, 1996 is perpendicular to the plane of the window between two supercages. The potential of mean force was determined from the mean forces in a series of MD simulations where the system was constrained to remain at the same value on the reaction coordinate. The results for the potential of mean force are similar at all temperatures. A temperature dependence of the height of the barrier for cage-to-cage migration was observed for very low temperatures only. The transmission coefficient and the absorbing-boundary factor defined in paper I14 were used to determine a local transmission coefficient and the probability for the continuation of the jump event. The mean number of barrier passages per jump event increases significantly when the temperature is raised. Whereas transition-state theory is a good approximation at very low temperatures, the consideration of the dynamical corrections increases the diffusion coefficient by 30% at 210 K. However, the contribution of the dynamical corrections to the Arrhenius activation energy is small. In the temperature range 77-210 K it does not exceed 0.3 kJ/mol. In the range from 140-210 K, the hopping-model results are in accordance with those of conventional MD simulations. Both the absolute value of the diffusion coefficient and the activation energy are in excellent agreement. Therefore, it seems appropriate to apply the hopping-model approach to the diffusion of other guest particles in zeolites. The application to diffusion processes with higher activation barriers will allow simulations to be used to investigate diffusion at room temperature for guest molecules that practically remain at the same adsorption site on the time scale of conventional MD simulations. 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) Laria, D.; Ciccotti, G.; Ferrario, M.; Kapral, R. J. Chem. Phys. 1992, 97, 378. (2) Verkhivker, G.; Elber, R.; Gibson, Q. H. J. Am. Chem. Soc. 1992, 114, 7866.

Mosell et al. (3) Hayoun, M.; Meyer, M.; Turq, P. J. Phys. Chem. 1994, 98, 6626. (4) Eyring, H. J. Chem. Phys. 1935, 3, 107. (5) Wigner, E. Trans. Faraday Soc. 1938, 34, 29. (6) Keck, J. Discuss. Faraday Soc. 33, 1962 , 173. (7) Bergsma, J. P.; Gertner, B. J.; Wilson, K. R.; Hynes, J. T. J. Chem. Phys. 1987, 86, 1356. (8) Chandler, D. J. Chem. Phys. 1978, 68, 2959. (9) Voter, A. F.; Doll, J. D. J. Chem. Phys. 1985, 82, 80. (10) Ruthven, D. M.; Derrah, I. H. J. Chem. Soc., Faraday Trans. 1 1972, 68, 2332. (11) Ka¨rger, J.; Pfeifer, H.; Haberlandt, R. J. Chem. Soc., Faraday Trans. 1 1980, 76, 1569. (12) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1991, 95, 8866. (13) Ka¨rger, J.; Ruthven, D. M. Diffusion in Zeolites and Other Microporous Solids; John Wiley: New York, 1992. (14) Mosell, T.; Schrimpf, G.; Hahn, C.; Brickmann, J. J. Phys. Chem. 1996, 100, 4571. (15) Ciccotti, G.; Ferrario, M.; Hynes, J. T.; Kapral, R. Chem. Phys. 1989, 129, 241. (16) Schrimpf, G.; Brickmann, J., J. Comput.-Aided Mater. Design 1995, 2, 49. (17) Fitch, A. N.; Jobic, H.; Renouprez, A. J. Phys. Chem. 1986, 90, 1311. (18) Demontis, P.; Suffritti, G. B.; Quartieri, S.; Fois, E. S.; Gamba, A. J. Phys. Chem. 1988, 92, 867. (19) Schrimpf, G.; Schlenkrich, M.; Brickmann, J.; Bopp, P. J. Phys. Chem. 1992, 96, 7404. (20) Kiselev, A. V.; Du, P. Q. J. Chem. Soc., Faraday Trans. 2 1981, 77, 1. (21) Bezus, A. G.; Kiselev, A. V.; Lopatkin, A. A.; Du, P. Q. J. Chem. Soc., Faraday Trans. 2 1978, 74, 367. (22) Santikary, P.; Yashonath, S.; Ananthakrishna, G. J. Phys. Chem. 1992, 96, 10469. (23) Manning, J. R. Diffusion Kinetics for Atoms in Crystals; van Nostrand: Princeton, 1968. (24) Carter, E. A.; Ciccotti, G.; Hynes, J. T.; Kapral, R. Chem. Phys. Lett. 1989, 156, 472. (25) Torrie, G. M.; Valleau, J. P. J. Comput. Phys. 1977, 23, 187. (26) Mezei, M.; Beveridge, D. I. Ann. N.Y. Acad. Sci. 1986, 482, 1. (27) Verlet, L. Phys. ReV. 1967, 159, 98. (28) Kast, S. M.; Nicklas, K.; Ba¨r, H. J.; Brickmann, J. J. Chem. Phys. 1994, 100, 566. (29) Miller, W. H. Acc. Chem. Res. 1976, 9, 306. (30) Ka¨rger, J.; Pfeifer, H.; Stallmach, F.; Feoktistova, N. N.; Zhdanov, S. P. Zeolites 1993, 13, 50. (31) Yashonath, S.; Santikary, P. J. Phys. Chem. 1993, 97, 3849.

JP9526455