8866
J . Phys. Chem. 1991, 95, 8866-8878
Figure 11. The Pg starts at a lower salinity at a higher temperature and eventually the point reaches the bottom (the ternary system without salt). Thus, the phase behavior in Figure 10 is observed above the IT. Upon addition of oil to the quaternary system, middle-phase microemulsions are produced via a tricritical point or a four-phase region. This will be published in the near future.
Conclusion 1. Phase behavior in a NaCl/water/R8S03Na/hexanol system was clarified above and below the intrusion temperature (IT) of a lamellar liquid crystal at which the lamellar liquid crystal (LC) intrudes into a main miscibility gap. There is a critical point PF in a main miscibility gap in a ternary water/R8S03Na/hexano1 system above a critical end temperature (about 76 “C) which is slightly lower than the IT (78.3 “C). The critical point does not exist below the critical end temperature in the absence of salt. 2. Upon addition of salt, the ionic surfactant becomes less hydrophilic. As a result, the original critical point, P;, in a water/R8S03Na/hexanoI system continuously shifts from an
alcohol-rich to water-rich regions above the IT. A micellar solution phase (Wm) splits into two phases, W and D’, in the presence of salt, and a new critical point, P&,and a three-phase triangle, W + D’ + Om, originated. This three-phase region remains even at higher salinity. 3. Salt unstabilizes the LC phase and depresses the IT. Upon addition of salt, the LC phase retreats to a concentrate region and isotropic phases (D’ and Om) are connected at a lower temperature below the IT. As a result, the critical point, P;, originates in a water-rich region even at the lower temperature. The critical points, Pk,of D’ + Om and the W + D’ + Om are also produced at lower temperatures upon addition of salt. 4. The composition of the Om phase forming the W D’ Om region is shifted toward a water-rich region with increasing salinity.
+ +
Acknowledgmenr. We thank Messrs. H. Arai and F. Harigai for their technical assistance. Registry No. R8S03Na, 5324-84-5; NaCI, 7647-14-5; I-hexanol, 11 1-27-3.
Transition-State Studies of Xenon and SF6 Diffusion in Silfcaiite
R. Larry June, Alexis T. Bell, and Doros N.Tbeodorou* Department of Chemical Engineering, University of California, Berkeley, Berkeley, California 94720 (Received: April 25. 1991)
Selfdiffusion coefficients for infinitely dilute spherical sorbate molecules in the zeolite silicalite are computed by transition-state theory. The diffusion process is modeled as a series of uncorrelated jumps between potential minima (sites) determined from a Lennard-Jones representation of the silicalite lattice. Rate constants for jumping between different sites within the lattice are computed by transition-state theory and by the dynamically corrected transition-state theory of Voter and Doll for both a flexible and a rigid zeolite lattice model. Diffusivities are then determined from the rate constants by generating continuous-time/discrete-spaceMonte Carlo random walks. The computed diffusivities are shown to be in good agreement with molecular dynamics calculations performed on an identical model in the infinite dilution limit at low temperatures and with available experimental results. The transition-state theory and dynamically corrected transition-state theory methods afford computational savings of up to 2 orders of magnitude relative to full molecular dynamics simulations. Shortcomings in the various algorithms and zeolite models are discussed.
Introduction Transition-state theory (TST) is a well-established methodology for the calculation of the kinetics of infrequent events in numerous physical systems. Examples of commonly studied phenomena amenable to the TST approach in nonreacting systems are the conformational isomerization of molecules,’-3 surface diffusion of adatoms on metal surfaces: adsorption-desorption kinetics of adsorbed surface specie^,^^^ and the diffusion of impurity atoms within host crystals.’ While the application of TST involves a series of approximations, the exact rate constant for the assumed phenomenology can be computed by factoring the rate constant for the physical process into two portions, a TST-based rate constant computed only from a knowledge of the equilibrium-phase space distribution of the system and a dynamical correction factor computed to compensate for errors in the TST rate constant? The dynamical correction factor, a multiplicative correction to the computed TST rate constant, is often on the order of unity. Evaluation of the correction factor requires detailed dynamical information about the system from methods such as molecular dynamics, but only over relatively short times compared to the time scales associated with “reactive” barrier crossings between states. In principle, the rate constants associated with these processes could be computed directly through the application of methods such as molecular dynamics, but the computational cost Author to whom correspondence should be addressed.
0022-3654/91/2095-8866$02.50/0
is often prohibitive due to the long time scales associated with moving through bottleneck regions in phase space. An instructive example is provided by a molecule such as n-hexane in the zeolite silicalite. Experimental data suggest that the diffusivity of nhexane at 325 K is about 2 X cm2/s.* If a hexane molecule must travel a distance of 1 unit cell (20 A) to exhibit true diffusive motion, a molecular dynamics run designed to predict the diffusivity should cover a simulation time of at least 300 ns. It has been estimated that approximately 4 X IO3 CPU hours would be required on a large-scale supercomputer for a single molecular dynamics simulation run of this d ~ r a t i o n . ~ The bulk of past simulation studies of sorbate diffusion in zeolites have been of two general types. The first approach found (1) Rosenberg, R.0.;Berne, B. J.; Chandler, D. Chem. Phys. L e r r . 1980, 75, 162-168.
(2) Chandler, D. J . Chem. Phys. 1978, 68, 2959-2970. (3) Czerminski, R.; Elber, R. J . Chem. Phys. 1990, 92, 5580-5601. (4) Voter, A.; Doll, J. J . Chem. Phys. 1985, 82, 80-92. ( 5 ) Grimmelmann, E. K.; Tully, J. T.; Helfand, E. J. Chem. Phys. 1974, 74, 5300-5310. (6) Adams, J. E.; Doll, J. D. J. Chem. Phys. 1981, 74, 1467-1471. (7) Bliichl, P. E.; Van de Walle, C. G.;Pantelides, S . T.Phys. Reo. Leu. 1990,64, 1401-1404. (8) Van-Den-Begin, N.; Rees, L. V. C.; Caro, J.; BOlow, M. Zeolifes 1989, 9, 287-292. (9) June, R.L.; Bell, A. T.; Thecdorou, D. N. Unpublished results based on simulations of alkane molecules in silicalite with a locally written alkane Molecular Dynamics code.
0 1991 American Chemical Society
Xenon and SF6 Diffusion in Silicalite
The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 8867
in the literature is based on stochastic simulations of particles (sorbate molecules) on a regular lattice.'&I2 In this Monte Carlo technique, a number of particles are placed on a grid which serves as an abstract representation of the cavities within a zeolite crystal. A random walk is generated by executing a series of random moves whose frequency of success is dictated by such factors as the local concentration, temperature, and pore network topology. This method is computationally efficient and can look at the effect on the sorbate diffusivity of such factors as intracrystalline loading, pore blockage, and concentration profiles. However, the time step corresponding to the rate at which elementary move attempts are made has invariably been assumed in an ad hoc fashion; thus, the method can only determine relative changes in the diffusivity. A second approach that is more recent, but computationally more intensive, is the application of molecular dynamics to the study of sorbate diffusion in zeolites. The method requires a detailed model of the zeolite structure, sorbate molecules, and the interaction potential between all species and produces realistic dynamical trajectories for an ensemble of molecules inside the pores of the zeolite. Molecular dynamics has been used to study methane in faujasite,I3 silicalite,'"16 and sodium A.I7 Other systems investigated by molecular dynamics include water in ferrieriteI* and benzene in Na-Y.I9 Diffusivities, time correlation functions, and thermodynamic properties can also be determined from the computed trajectories. However, if significant bottlenecks to motion through phase space exist, the sorbate molecules in a molecular dynamics simulation may not be able to sample a large enough portion of phase space in a reasonable amount of computer time to obtain correct self-diffusivities. For a system whose evolution in configuration space is sluggish due to the presence of high-energy barriers, a molecular dynamics simulation would do little more than track local motion within an energetically favorable region of that space and would be incapable of accumulating statistics for long-range diffusive motion. Many zeolite-sorbate systems have diffusivities that range in magnitude from IO-" to cm2/s.20 As an order of magnitude estimate, traditional molecular dynamics simulations are limited to predicting diffusivities in systems with diffusivities greater than 10-6 cm2/s, and as a consequence, sorbate diffusion in many zeolite-sorbate systems of practical interest cannot be studied by molecular dynamics. Transition-state theory methodologies provide a bridge between simple stochastic simulations and atomistic molecular dynamics simulations. Computations utilizing the TST methodology have been attempted in zeolites in the past with variable success. Sargent and Whitford2' undertook the calculation of the diffusivity of C02in zeolite 5A by TST. Errors in the computed diffusivities as large as 104 were attributed to an assumption of a rigid zeolite lattice and, more significantly, to the neglect of any preferred molecular orientation for the linear C02molecule in the bottleneck region of configuration space. A similar computation was un-
-
(IO) Tsikoyiannis, J. G.; Wei, J. Chem. Eng. Sci. 1991, 46, 233-255.
( I I ) Theodorou, D. N.;Wei, J. J . C a r d 1983, 83, 205-224. (12) Aust, E.; Dahlke, K.; Emig, G. J . Caral. 1989, 115, 86-97. (13) Yashonath, S.; Demontis, P.; Klein, M. L. Chem. Phys. Lett. 1988, 153, 551-556. (14) Demontis, P.; Yashonath, S.; Klein, M. L.J. Phys. Chem. 1989, 93, 5016-501 9. (15) Trouw, F. R.; Iton, L. E. In Zeolites for rhe Nineties, Recenr Research Reporrs; Jansen, J. C., Moscou. L.. Post, M. F. M.,Eds.; presented
at the 1989 International Zeolites Conference, Amsterdam, The Netherlands, June 1989; pp 309-310. (16) June, R. L.; Bell, A. T.; Theodorou, D. N. J . Phys. Chem. 1990.94, 8232-8240. (17) Cohen de Lara, E.; Kahn, R.; Goulay, A. M. J . Chem. Phys. 1989, 90, 7482-7491. (18) LehertC, L.; Varcauteren, D. P.; Derouane, E. G.; Lie, G. C.; Clementi, E.: Andrt, J.-M. Zeolires: Facrs, Figures, Future; Jacobs, P. A,, van Santen, R. A., Eds.; Elsevier: Amsterdam 1989; pp 773-783. (19) Demontis, P.;Fois, E. S.;Suffritti, G. 8.;Quartieri, S . J. Phys. Chem. 1990, 94,4329-4334. (20) Karger, J.; Ruthven, D. M. Zeolites 1989, 9, 267-281. (21) Sargent, R. W.; Whitford, C. J. In Molecular Sieve Zeolites II; Gould, R. I., Ed.; Advances in Chemistry Series 102; ACS: Washington, DC, 1971.
Figure 1. Schematic outline of the pore structure of a unit cell of silkalite. Spheres representing the three types of sorption states (2 = sinusoidal-channel state, S = straight-channel state, and I = intersection state) on the zeolitesorbate potential hypersurface are shown. The thick lines denote the direction of the straight- and sinusoidal-channel systems.
dertaken by BCtemps and Jutard for argon in zeolite 5A.22 Poor agreement between the calculated and measured diffusivities in this case was attributed to uncertainties in the location of Na' ions in zeolite lattice. Ruthven and Derrah23also utilized TST to compute diffusivities for CH, and CF, in zeolite SA. Reasonable results were obtained for these spherically symmetric molecules provided that rotational degrees of freedom were taken into account. None of the above studies of zeolitesorbate systems has attempted to account for barrier recrossing events using a dynamical correction factor. In this paper we report on our calculations employing dynamically corrected transition-state theory to determine the self-diffusivity of two Lennard-Jones spheres sorbed at infinite dilution in the zeolite silicalite. The smaller of the two spheres was chosen to represent xenon, while the larger sphere roughly approximates SF6. At the low temperatures investigated here, potential energy barriers presented to the sorbate molecules by the zeolite lattice enhance the localization of sorbate molecules within discrete regions of the pore network constructed around the potential minima of the silicalite lattice. Algorithms are presented for determining the shape and extent of sorption sites, as well as the values of the fundamental rate constants, kU,for jumping between different sites within the silicalite pore network. Comparisons are then made between the computed TST diffusivities and the diffusivities derived from full molecular dynamics calculations where possible. Within the assumptions of transition-state theory, the results of the simulations are also compared at different levels of approximation. In the case of xenon, additional comparisons are made between diffusivities computed for a rigid lattice model and an animated (flexible) lattice model. The sections below present the model used to represent the sorbate-silicalite system, the relevant aspects of transition-state theory, the computational methodology, and the results of the computations. A set of concluding remarks are presented at the end of the paper.
Model Representation The crystal structure of silicalite forms a three-dimensional interconnected pore network through which guest sorbate molecules can diffuse. The structure of the silicalite lattice is wellknown from X-ray diffraction studies.24 Transformations between orthorhombic and monoclinic lattice symmetries can be induced by addition of aluminum to the lattice, sorption of large sorbate - * ~this molecules, or elevation of the lattice t e m p e r a t ~ r e . ~ ~ In (22) BCtemps, M.; Jutard, A. J. Phys. D Appl. Phys. 1980, 13, 423. (23) Ruthven, D. M.; Derrah, R. 1. J . Chem. Soc., Faraday Trans. I 1972, 68. 2332-2343. ..._.._ ~. . (24) Olson, D. H.;Kokotailo, G. T.; Lawton, S. L.; Meier, W. M. J. Chem. Phys. 1981, 85, 2238-2243. (25) Klinowski, J. T.: Carpenter, T. A.; Gladden, L. F. Zeolites 1987, 7, 73-78.
8868 The Journal of Physical Chemistry, Vof. 95, No. 22, 1991 TABLE I: Coefficients for the Attraction and Repulsion Terms of the Lenmrd-Jonw Potential for the Xewrr-oXygen and SF-ygen Systems and Panmeten for the Bonded Zeolite S i 4 Potential
interacting pair 0-Xe 04F6 interacting pair Si-0
Sorbate repulsion, kcal A’*/moI
attraction, kcal A6/moI
0.3113 X IO7 1.7776x 107
1836 4560
Lattice kb,kJ/(mol-A)
ro, A
465.0
1.605
work and in previous efforts,I6a we have considered only the phase with orthorhombic symmetry. Figure 1 contains a schematic drawing of the pore structure of silicalite which emphasizes the connectivity of the channel system. Sinusoidal channels follow a contour along the [ 1001 direction, and the straight channels trace a path along the [OlO] direction. Motion in the [OOl]direction is possible only by a tortuous series of movements between straight and sinusoidal channels. In our rigid lattice model, all silicon and oxygen atoms in the zeolite framework are considered as occupying fixed positions and perfect crystallinity is assumed. The animated zeolite lattice model allows for thermal motion of the framework atoms under the influence of forces from neighboring lattice atoms and sorbate molecules. As in our previous work, all interactions between the sorbate atoms and the oxygen atoms of the silicalite lattice are treated atomistically with a Lennard-Jones 6-1 2 potential. The silicon atoms, being recessed within the SiO, tetrahedra, are neglected in the potential summation.29 For xenon-oxygen interactions, the Lennard-Jones dispersion constant was determined with the Slater-Kirkwood equation, which requires knowledge of the effective number of electrons, ne,and atomic polarizabilities, a. The repulsion constant was determined from a force balance at the sum of the van der Waals radii for the xenon-oxygen pair, rt. Lennard-Jones e and u values for SF6 were taken from Reid et aL30 Standard geometric and arithmetic mean combining rules were used for the SF6 and 0 c and u parameters. As no experimental data are available for SF6 in silicalite for comparison against simulation results, less emphasis was placed on the accuracy of the interaction parameters for SF6. As will be shown in later sections, the large effective diameter of SF6was expected to change the probability of occupancy and connectivity of the lattice of potential minima within the silicalite network. Interaction parameters for the attraction and repulsion terms of the Lennard-Jones potential are given in Table I for both sorbateoxygen systems. The total potential energy, Yd*, of a sorbate molecule in the absence of neighboring sorbate molecules is represented by a single sum over all lattice oxygens within the cutoff distance of 15 A as
For the animated lattice, the sorbate molecule and lattice atoms interact in a similar fashion. An additional harmonic potential term arises from all bonds associated with the zeolite framework: y m l i t c =:
C Y2kb(rb- r d 2
bEz
(2)
where b is an index that runs over all bonds formed in the sim(26)Fyfe, C.A.; Kennedy, G. J.; Kokotailo, G. T.; Lyerla, J. R.; Fleming, W. W. J. Chem. Soc., Chem. Commun. 1985,740-742. (27)Hay, D.G.; Jaeger, H. J. Chem. Soc., Chem. Commun. 1984, 1433. (28) June, R. L.; Bell, A. T.; Thcodorou, D. N. J . Phys. Chem. 1990,94, 1508-1 5 16. (29)Kiselev, A. V.;Lopatkin, A. A,; Shulga, A. A. Zeolites 1985, 5 ,
June et al. ulation cell. The force constant, kb, and equilibrium bond length, ro, for the lattice S i 4 bonds were taken from previous molecular dynamics simulations in silicalite and zeolite AI4 and are reported here in Table I.
Theory The determination of the diffusivity of a sorbate molecule with the transition-state methodology begins by assuming that the diffusive motion of the sorbate molecules through the zeolite occurs by a series of uncorrelated hops between potential minima in the zeolite lattice. Around each minimum is constructed a sorption state. A first-order rate constant, ky, is then associated with the rate of transition between a given pair of neighboring states, i and j . Any neighboring pair of states, i and j , is then assumed to be separated by a dividing surface on which a saddle point of the potential energy hypersurface is located. The saddle point, s, can be viewed as the transition state between states i and j; a couple of steepest descent paths from the saddle point, s, connect the minima associated with states i and j . After rate constants, k,,, have been determined for all possible transitions in the lattice of sorption states, a continuous-time/discrete-spaceMonte Carlo calculation is used to determine the self-diffusivity of the sorbate molecules. This Monte Carlo calculation generates random walks of sorbate molecules on the lattice of potential energy minima with the hopping rates determined from TST computations. In this section, expressions for the transition-state rate constant kET, which determines the rate of escape from state i to all other states, for kFT, which determines the rate of escape from state i to a specikc state j , and for the dynamically corrected rate constant k - are presented. A continuous-time/discrete-spaceMonte Carlo afgorithm for the determination of sorbate self-diffusivities is then described. Classical transition-state theory assumes that a molecule is coupled to a heat bath, in this case, the atoms of the zeolite lattice. This assumption implies that while residing near the potential minima and also at the transition state, the molecule has distributed energy among its available degrees of freedom in accord with the canonical ensemble density function. Furthermore, TST assumes that all molecules approaching the dividing surface with sufficient kinetic energy to cross the potential barrier between states will undergo a reactive, or successful, crossing of the potential barrier.” Clearly, if the coupling of the molecule with the heat bath is weak, the molecule may fail to thermalize in the destination state and cross into either the reactant state or a second product state. Conversely, if the coupling of the molecule to the bath is strong in the region near the potential barrier, the trajectory may fail to cross the transition state. Such dynamical recrossing or multistate transition events are not addressed by conventional TST. Finally, at infinite dilution, sorbate molecules are assumed to move between sorption sites with first-order kinetics. Within the framework of these assumptions, the transition-state theory rate constant, denoted by kET, describes the rate at which atoms leave state i throu h the dividing surface surrounding state i. Note, however, that k,, contains no information about the destination state j whatsoever and only describes the total rate of egress from state i. In the formulation of Voter and Doll: any dependence of the dynamically corrected rate constant, k,, on the physical characteristics of statej is contained within a quantity known as the dynamical correction factor, defined below. The transitionstate theory estimate of the rate constant can be written as an average over the region of configuration space which corresponds to state i only:4
4
kTS7 j-.
= (2*Bm)-”2(8i(o))i
(3) where m is the atomic mass of the sorbate molecule, B is l/keT, and Si(0) is a S function which is defined to be nonzero only on the surface surrounding state i. The average represented by the brackets is taken in the canonical ensemble only over that subset of configuration space which consists of state i and its boundaries.
261-267.
(30) Reid, R. C.; Prausnitz, J. M.;Poling, B. E. The Properties of Gases and Liquids; McGraw-Hill: New York, 1987.
(31) Weiner, J. H. Statistical Mechanics ofELniciry; Wiley: New York,
1983.
Xenon and sF6Diffusion in Silicalite
The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 8869
In applying this formalism to SF6,we assume that the rotational degrees of freedom of the molecule are distributed in accord with the canonical ensemble distribution function and not hindered in any way at the transition state. For the simple systems under consideration here, the average is most easily evaluated if the 6 function is expressed by reformulating the boundary conditions of the integrals used in evaluating the ensemble average: JYP(-BY)
kFT = (2a/3m)-'I2
dr2
sY,yP(-BY) dr3
(4)
where r represents the configurational degrees of freedom in in the system. For a simple spherical sorbate molecule in a rigid silicalite lattice, only three translational degrees of freedom are present. The integral in the numerator of eq 4 is evaluated on the surface Si surrounding state i that is accessible for transitions to neighboring states, while the integral in the denominator is taken over the volume in configuration space occupied by state i. A transition-state theory estimate of the rate constant for jumping between state i and a specific state j , k r T , can be computed by evaluating the integral in the numerator of eq 4 only over the boundary surface shared by states i and j . The rate constant for the reverse process, kFT, is computed in a similar fashion, except that the limits of integration in the denominator are defined by the configurational volume of statej. In defining sorption states and transition pathways with the animated zeolite lattice, the potential hypersurface was taken as identical with the one dictated by the equilibrium X-ray coordinates used in the rigid lattice calculations. The following balance is always true:
where the index j runs over all states sharing surfaces with state i. By definition, the constants kTT satisfy the condition of microscopic reversibility: where p,the equilibrium probability of occupancy of sorption state i, is defined as
(7) A more powerful methodology for dealing with the dynamics of infrequent events is provided by dynamically corrected TST. We follow here the formalism presented by Voter and Doll4 for multistate systems with many degrees of freedom. As shown by these authors, an "exact" value of the first-order rate constant, k , incorporating dynamical recrossing and multistate jump effects, can be obtained from this formalism. The theory only assumes that an explicit separation in time scales exists between 7cwr,the time constant associated with the thermalization of a molecule after a reactive barrier crossing, and 7,. the average time between reactive barrier crossings. Postulating first-order interstate transition kinetics, the rate constant for the transition between two states can be expressed as the product of the time-de ndent portion, k p , and a time-dependent correction term,
&):
where the quantity to the left of k z T is defined as the dynamical correction factor, f i ( r ) . The term u,(O) represents the velocity of the sorbate molecule at the point of the bottleneck dividing the two states i and j , where the molecule is assumed to be at time zero, in the direction of the local normal to the bottleneck surface. The quantity O,(t) is a switching function, defined as 1 whenever a given molecule finds itself in the destination state j at time t and as 0 otherwise. The average denoted by the brackets is sampled over a phase-space distribution, w; wi = bi(O;r=O) exp(-@H)
(9) where the Hamiltonian of the zeolitesorbate system, denoted by H , is simply the sum of the kinetic and potential energies of the
system. The 6 function 6,(0;t=0) specifies that all sorbate molecules are located on the dividing surface at time zero. Conventional TST corresponds to placing @At)= e,(O+), in which case the dynamical correction factor degenerates into the Boltzmann weighted fraction of the bottleneck surface surrounding state i that leads to statej. Note that, while the dynamical correction term f i ( t ) may rapidly depart from the initial value of fi(0') corresponding to conventional TST, it reaches a plateau value between 0 andfi(0') in times on the order of 7arr. To evaluate f ( t ) ,a number of mutually "transparent" sorbate molecules are #rst placed on the dividing surface surrounding state i by random insertion. The molecules are then thermalized over the entire bottleneck surface through a number of Metropolis Monte Carlo passes which result in distributing them over that surface according to the Boltzmann distribution, exp(-/3Ymrktc). Initial momenta are assigned to all molecules according to an isotropic Maxwell-Boltzmann distribution a t the temperature of interest. Subsequent motion of the ghost molecules is tracked by microcanonical ensemble molecular dynamics. The dynamical correction factor is evaluated as a function of time by taking the trajectories generated from a series of molecular dynamics runs and computing the average indicated in eq 8. It is important to note that the total duration of the simulation run is dictated by 7", or the time required for the sorbate molecules to thermalize in a new state; beyond this time,f.(t) reaches a plateau value and experiences a slow exponential &cay characterized by a time constant of 7rxn. Thus, the corrective molecular dynamics simulations need only be run over a time scale in the neighborhood of 7arr, not a few 7,," as would be necessary to investigate the dynamical behavior of the system with a full conventional molecular dynamics simulation. Although the rate constants provide a fundamental description of the hopping process between states, they are not easily directly measurable; it is desirable to compute a self-diffusivity which can be measured e ~ p e r i m e n t a l l yand ~ ~ ,compared ~~ against full MD simulation results. For a simple cubic lattice with a constant jump rate, k+., and an equal probability of residence in any lattice site, the diffusion constant is simply D, = ki,Ax2/6
where Ax is the distance between lattice sites. In the silicalite lattice, however, this simple formalism is not applicable; the lattice of sorption sites formed by the potential minima is not simple cubic in nature, the probability of residing in any given site is dependent on the type of site, and the rate constants, kij, are, in general, different for each ij pair. While an analytical solution to the master equation describing this system is a Monte Carlo calculation based on the first-order kinetic description of the hopping process easily allows the determination of the diffusivities to an arbitrary accuracy. However, it has been shown in studies of surface diffusion that some care must be taken in the choice of transition probabilities and the interpretation of move attempts in relation to simulation time.35v36 The Monte Carlo algorithm functions by allowing sorbate molecules to move between adjacent sites with frequencies determined by the fundamental rate constants, either k y T or kip All of the allowed types of moves between states are defined by the connectivity of the states as determined in the TST calculations. The stochastic process of diffusion is assumed to be governed by Poisson statistics, which dictate that the probability of a move event from a state of type i to a state of type j occurring in a time interval between t and t + dt since the last move between those types of states is given by
(32) Heink, W.; Kirger, J.; Pfeifer, H.;Stallmach, F. J . Am. Chem. Soc.
1990, 112, 2175-2178.
(33) BOlow, M.; HBrtel, U.; Muller, U.; Unger, K. K. Ber. Bunsen-Ges. Phys. Chem. 1990, 94, 74-76. (34) Haus, J . W.; Kehr, K. W. Phys. Rep. 1987, 150, 263-416. (35) Kang, H. C.; Weinberg, W. H.J.Chem. Phys. 1989,90,2824-2830. (36) Ray, L. A.; Baetzold, R. C. J . Chem. Phys. 1990, 93, 2871-2878.
8870 The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 where pij is the rate parameter. The rate parameter is taken as equal to the first-order rate constant k, for a jump between states i and j multiplied by the number of molecules in all states of the ith type; the latter multiplication gives a rate parameter scaled for the system size. The average time between ij jump events, ( t)ij, is found to be 1 / p i j by computing the expectation value of the time between events from the distribution Fij(t). The assumption of Poisson statistics implies that the move events are uncorrelated with each other and essentially random in nature. Two properties of exponentially distributed random variables prove useful in the computation.1° The first is the overall transition rate parameter for a series of elementary events. Given n possible types of independent moves, or elementary events, with each event occurring at some time interval, t,, the interval of time after which the first move occurs is tmin = min ( t l , t 2,...,t,,)
(12)
The mean frequency of such moves is described by the overall transition rate parameter: P =
CPi
i t1
Equation 13 shows that the overall transition rate parameter is simply the sum of the individual rate parameters. The second property is that the probability Pi of event i occurring first is given by n
j= I
Finally, it can be shown that the Poisson character of the move process dictates that the system evolution is Markovian in nature.37 The Poisson process description leads to a rather simple algorithm for generating the sequence of jumps to simulate diffusion at low loadings. A random number &, between 0 and 1, is picked from a uniform distribution. The overall transition rate parameter for the system in its current state, p, is determined from eq 13. The time interval for the occurrence of the next elementary event is then chosen by integrating and inverting eq 11 for the time interval T elapsed until the occurrence of the next event:
A time increment T is then added to the global simulation time counter. Next, a particular type of move is chosen by comparing a second random number,F2, with the cumulative probability, P,, of each of the individual event types: N.
P, = 2Pi i= I
where Ne is the total number of allowed move events. The first i for which the sum defining Pq is larger than the random number l2defines the type of event from which the next move is chosen. As all molecules in a given type of state are equivalent at infinite dilution, a molecule is then chosen at random from the list of molecules in state i and moved to the new destination state j . Mean-square displacements are accumulated as the system evolves in time by sampling new states. The Einstein relation
is used along with similar expressions for D,,and D, to determine the diagonal elements of the self-diffusivity tensor. The orientationally averaged trace of the self-diffusivity tensor is determined by38
D, =
D, + Dy+ D, 3
(18)
(37) Feller, W. An Introduction to Probability Theory and its Applications. Volume II; Wiley: New York, 1957.
June et al. The theoretical estimate D,corresponds to the diffusivities determined experimentally by pulsed-field gradient NMR.
Computational Methodology Determination of State Volumes and Dividing Surfaces. Search for minima on the zeolitesorbate potential hypersurface is greatly simplified by the assumption of a rigid lattice. This assumption reduces the number of degrees of freedom to the three translational degrees of freedom associated with the sorbate molecule. A three-step process was employed to efficiently and exhaustively locate all minima on the potential energy hypersurface. In a first step, the volume of the asymmetric unit of the silicalite unit cell was discretized by a grid with a spacing of approximately 0.2 A. Due to the orthorhombic symmetry of the silicalite lattice, searching through only the asymmetric unit reduces the volume of configuration space to be searched by a factor of 8. The potential and gradient vector were then tabulated on the grid for the accessible regions of the asymmetric unit. As a matter of practice, regions of extremely high potential were excluded from the calculations to minimize computational effort because the high potential energy indicated the presence of nearby lattice atoms. In the second step, the nodes of the grid were scanned. Any cube formed by a set of nearest-neighbor grid nodes was checked for changes in sign of the three components of the gradient vector on any of the eight vertices of the cube. In the last step, if all three components of the gradient vector were found to change sign on two or more vertices of the cube, a search for the local potential energy minimum was initiated with the coordinates of the center of that cube. A quasi-Newton minimization search algorithm employing the BFGS matrix updating schemeN was used. It was later determined, however, that the aforementioned grid search algorithm placed the initial guess close enough to the minimum that iterations with a simple Newton procedure utilizing analytical second derivatives of the zeolite-sorbate potential converged quickly and accurately. Two different algorithms were evaluated for determining the location of saddle points on the zeolite-sorbate potential hypersurface. The first algorithm searches for global minima in the function gg, where g is the gradient vector of the zeolitesorbate potential hypersurface, ?Isarbate.41 As before, the hypersurface described by g g is searched for minima by using the BFGS algorithm. A disadvantage of this algorithm is that, if minimization techniques requiring second derivatives of the objective function are to be used, third derivative information, which can be expensive to compute, must be available for the zeolite-sorbate potential hypersurface. Furthermore, the function g g has a multiplicity of minima associated with the already determined minima on the zeolitesorbate potential hypersurface, ?Ihk. However, minima on the g g hypersurface are easily discerned, and after a simple examination of the eigenvalues of the Hessian matrix of ?Isarbab, an unambiguous determination of the identity of the minimum or saddle point can be made. The second algorithm, due to Baker,42searches the zeolite-sorbate potential hypersurface directly for saddle points. Baker's Newton-like algorithm accomplishes this by maximizing the objective function ?Isarbatc along the eigen direction associated with the smallest eigenvalue and by minimizing along directions associated with all other eigenvectors of the Hessian of the zeolitesorbate potential hypersurface. Some disadvantages are associated with this algorithm, too. Among the pieces of information required to start the Baker search is an eigenvector/eigenvalue pair that points, in some sense, in the direction of the nearest saddle point. The method also requires all eigenvectors and eigenvalues of the Hessian matrix at each (38) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1987. (39)Shen, D.; Rees, L. V. C. J. Chem. Soc. Faraday Trans. 1990, 86, 3687-3692. (40) Hillstrom, K.Nonlinear Optimization Roulines in AMDLIB, Technical Memorandum No. 297, Argonne National Laboratory, Applied Mathematia Division, Argonne, IL, 1976;Subroutine GQBFGS in AMDLIB. (41) Stillinger, F. H.; Weber, T. A. J . Chem. Phys. 1984.81, 5095-5103. (42)Baker, J. J. Comput. Chem. 1986,7 , 385-395.
Xenon and SF6 Diffusion in Silicalite step of the calculation; this is not a problem for a system with only three degrees of freedom, but may become a problem for larger systems. In practice, both algorithms were found to work, but Baker's algorithm was more robust and gave much more rapid convergence to the saddle points. To associate each saddle point with the transition state between two specific minima, a series of steepest descent calculations were performed. While this calculation is not strictly necessary using the formulation of Voter and Doll on a simple potential surface, it does provide insight into the connectively and topology of the network of states. The steepest descent trajectory was started by perturbing the initial placement of a test particle by 8, from the saddle point in either direction along the negative eigenvector of the Hessian. Steps of 0.01 A were then taken successively in the direction -g. After a fixed number of iterations, usually 500, the steepest descent calculation was terminated and a Newton iteration scheme was started to accurately locate the potential minima associated with the transition state. This calculation creates a table of triplets connecting each pair of states, i a n d j , with a saddle point, s. In this manner, the connectivity of the network of states present in the unit cell was determined. It should be noted that, while all saddle points were found within the asymmetric unit, the steepest descent trajectories, in some cases, terminated outside the asymmetric unit. In evaluating the rate constants according to transition-state theory, a geometric definition of each state and of the border separating it from neighboring states in configuration space is required to define the limits of integrals such as those of eq 4. It is required that these states completely tile the accessible portion of configuration space? It is also desirable from a computational viewpoint that the states be as geometrically simple as possible. Furthermore, it is required that all transitions between two states pass through the dividing surface between them. Through an optimal choice of the dividing surfaces between states, one can minimize the predicted value of the transition-state rate constant and maximize the value of any dynamical correction to the computed rate constant. Each state forms an enclosed region in configuration space, bounded either by impenetrable portions (walls) of the zeolite or by windows connecting two states. For many systems studied in the past with TST-based approaches, the boundaries between states could be defined by simple geometric constructions. Examples are found in studies of the conformational isomerization of butane, where the boundary is a point on a one-dimensional reaction coordinate,' or in the diffusion of metal adatoms on an ordered metal surface, where the dividing plane is formed by the boundary associated with the unit cell of the underlying metal ~ u r f a c e .In ~ the silicalitmrbate systems studied here no such simple geometric construction is possible, as the symmetry present in the threedimensional potential energy hypersurface is complex and is characterized by many symmetry elements. In this case it was found to be most advantageous to discretize the pore volume by using a three-dimensional grid and to assign the resulting discrete volume elements (voxels) to one of the states. The test used to determine to which state each voxel belonged was based on a steepest descent calculation similar to that used to determine the connectivity between minima and saddle points. A steepest descent trajectory was originated in the center of each voxel; membership to a specific state was assigned by the identity of the energy minimum in which this steepest descent trajectory terminated. By use of this definition, the volume of a region around a minimum (state) was determined simply as a union of the volumes of the individual voxels that comprise the state. The dividing surfaces between states were determined by looping through all voxels and investigating which state each of the six nearest neighbors of a given voxel belonged to. If all the neighbors belonged to the same state, the voxel under examination was designated as lying in the interior of the state. If one or more faces of a voxel bordered on voxels assigned to differing states, the area associated with each of those common faces was added to a list containing the window between the two states. Any face on a voxel that had no neighboring voxel associated with it was
The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 8871 considered to border an impenetrable portion of the pore wall. Calculation of Rate Constants. The calculations necessary to evaluate the transition-state theory estimates of the rate constants, k z T , were performed by Monte Carlo integration t e c h n i q ~ e s . ~ ~ The integrals in the numerator of eq 4 were performed over all facets of the volilme elements which formed the boundary between each state. For the xenon system, a total of 200 function evaluations were used for each facet in the Monte Carlo integration scheme. For the SF6system, which formed a much more tortuous potential hypersurface, 2000 integration points were required to achieve the same estimated 2% level of error in the integral. Similarly, the volume integral in the denominator of eq 4 was evaluated to an uncertainty of about 1% by utilizing 200 and 2000 function evaluations per voxel in the xenon and SF6 systems, respectively. Dynamical corrections to the TST rate constants were calculated by undertaking short molecular dynamics runs with carefully chosen initial conditions. As described in the previous section, the dynamical correction for the rate constant kET begins by distributing a number of molecules over the faces of the volume elements forming the surface of state i. In each simulation run, 1000 noninteracting molecules were randomly inserted on the facets comprising the surface of the state. Subsequently, 1000 Metropolis Monte Carlo passes were used to thermalize the molecules and distribute them according to the canonical ensemble over the surface of state i. A combination of small conservative displacements and aggressive interfacet jumps was used in this Monte Carlo procedure. The initial velocity vectors were simply assigned from a Maxwell-Boltzmann distribution a t the temperature of interest. With the initial conditions set, the equations of motion were integrated for 10 ps with a Verlet algorithm for the rigid lattice model and with a Gear fifth-order predictorcorrector algorithm for the animated lattice model. The integration of the dynamical equations proceeded serially for each sorbate molecule in the animated lattice model, and in parallel for all noninteracting molecules in the rigid zeolite model. Integration time steps for each model were set to 0.02 and 0.001 ps in the rigid and animated lattice models, respectively. Trajectories were written every 0.1 ps to disk for later analysis. The simulation cell for the animated lattice consisted of two silicalite unit cells and utilized periodic boundary conditions for all atoms, while the rigid lattice model was run without periodic boundary conditions in the spatially periodic zeolite-sorbate potential field. For the animated lattice simulations, any drift in the center of mass of the parent coordinates of the lattice atoms was zeroed before each trajectory calculation. Additionally, the temperature of the lattice was reset to the simulation temperature before each trajectory calculation to prevent heating from the thermalization of the sorbate molecules. The rigid lattice model was simple enough to allow 10 runs of 1000 molecules for each state while the computationally more intensive animated lattice simulations were limited to only three runs of 1000 serially computed trajectories for each state. According to the prescription of Voter and Doll (eq 8), correcting the rate constants in the rigid lattice model was a computationally efficient process. However, it should be noted that the process of distributing molecules over the dividing surface according to a Boltzmann distribution naturally leads to a nonuniform distribution of sorbate molecules on the dividing surface. Energetically favorable portions of the dividing surface were assigned a vast excess of molecules and produced excellent statistics for the computed dynamical correction factors associated with that region of the dividing surface. In general, transitions between states that were separated by energetically unfavorable portions of the dividing surface were assigned only a few sorbate molecules and led to correspondingly poorer statistics in the computed dynamical correction factors. To circumvent this and to preserve a detailed balance between each pair of opposite transitions, it (43) Press, W.H.;Flannery, B. P.; Teukolsky, S. A,; Vetterling, W. T. Numerical Recipes: The Art of Scientific Computing Cambridge University Press: London, 1986.
June et al.
8872 The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 was sometimes necessary to adjust the value of the statistically less significant dynamical correction factors so that the detailed balance condition, eq 6, was obeyed. The dynamically corrected rate constants, k,, calculated as discussed above, were used as inputs to the continuous-time/ discrete-space Monte Carlo code. To compute statistically significant mean-square displacements, lo00 ghost sorbate molecules were placed on a lattice of sites representing the sorption states. The initial distribution of sorbate molecules among the sites was chosen in such a fashion as to conform to the equilibrium probability distribution among the states (PIdetermined in the course of the transition-state theory calculations. Note that, in the absence of sorbatesorbate interactions, any distribution of sorbate molecules that conforms to the equilibrium distribution of states is equally likely. The lattice of sorption sites within the crystal is spatially periodic and can be considered as infinite in extent for the time scales considered here. Monte Carlo trajectories of the sorbate molecules were followed over 8 X IO5 moves. In all cases, the time duration of the simulation exceeded by far the largest time constant for a transition, l / k i j With this algorithm, there is no need for equilibration and the diagonal components of the diffusivity tensor were computed directly from the slope in a plot of the mean-squared displacements along each coordinate axis versus time according to eq 17. Two types of more traditional simulations, molecular dynamics and Metropolis Monte Carlo, were carried out for verification and comparison of the transition-state theory calculations. These simulations were carried out to investigate a number of points. First, it was desired to investigate specifically the effect of sorbate “dimerization” at low loadings and low temperature. A strong tendency to dimerize at low loadings would invalidate the concept of a potential energy hypersurface defined by the coordinates of a single sorbate molecule. Second, a recent study has shown that the thermalization of sorbate moecules at low loadings is a slow process due to the large potential field exerted by the zeolite.44 This is extremely important for sorbate molecules in a rigid lattice where thermalization only comes from interacting with neighboring sorbate molecules. Finally, diffusivities and site occupancies computed by these more traditional means could be compared against those determined by transition-state theory. Metropolis Monte Carlo simulations for the xenon-silicalite system were carried out at 150 K for sorbate loadings of 0 and 1 molecule/unit cell, while molecular dynamics simulations in a rigid lattice were carried out a t a series of temperatures between 100 and 300 K and loadings of 0.5, 1 .O,and 2.0 molecules/unit cell. These traditional simulations were not attempted in the SF6 system which exhibits laige energy barriers between states and is therefore governed by extremely long equilibration times. In all cases, the simulation model was identical with that used in the transition-state theory calculations. The Metropolis Monte Carlo simulations utilized 51 2 sorbate molecules, and statistics were gathered by making between 50K and IOOK passes over all sorbate molecules. Results were written to disk every 10 passes. In the molecular dynamics simulations, our previous code2*was used to integrate standard Newtonian equations of motion for 128 molecules for time durations of up to 6OOO ps in length; trajectories were recorded every 0.2-2 ps, depending on the simulation temperature. Uncertainties in the computed properties were estimated by block averaging the results over 10 independent blocks of data gathered during the simulation. Results As outlined in the preceding section, a number of computational
steps are necessary to produce TST estimates of the diffusivity of the sorbate molecules in silicalite. In this section, the results of each of the intermediate computations are discussed. Specifically, we discuss the location of minima and saddle points of the potential energy hypersurface, the partitioning of configuration space into states, the evaluation of the rate constants from eqs (44) Fritzsche, S.: Haberlandt, R.; Kirger, J.; Pfeifer, H. Chem. Phys. Leff. 1990, 171, 109-112.
TABLE II: Minima and SIddk Points of the Potential Energy Hvwrsurface of Xenon in SilicrliteO no. X Y 2 1 9.18 4.16 11.87 2 0.05 3.91 8.39 3 0.14 3.80 6.45 4 8.63 4.98 12.86 5 9.22 4.98 1.01 8.01 6 0.59 4.98 1 5.05 4.98 1.32 6.7 I 8 0.00 0.00 9 2.34 4.98 1S O IO 0.36 0.93 6.1 1 2.13 II 4.00 4.98 5.85 4.98 12 0.88 13 1.43 4.98 1.56 14 0.46 4.98 6.06 ~~
energy
location
-4.103 -4.238 -4.749 -4.226 -4.306 -4.140 -5.773 -5.695 -5.606 -5.876 -5.956 -5.842 -5.146 -5.1 1 1
110.1.12) (10;2;13j (14,3,10) (14.4.12) i13[5112j (1 4,6J 3)
LL LL LL
straight sinusoidal sinusoidal sinusoidal intersection
‘Entries 10-14 contain the potential energy minima found by the and Newton minimization techniques. Saddle points, located with Baker’s algorithm, are contained in entries 1-9. The coordinates, in A, of the xenon molecule in the asymmetric unit are given along with the energy, in kcal/mol, associated with each minimum or saddle point. The minima are also designated as being located in the sinusoidal-channel, straight-channel, or intersection region of the unit cell. The triplet ( i , s j ) indicates what states are connected by the saddle points; one of these states may be an image located outside the asymmetric unit. Note that low-lying saddle points are marked as LL and are ignored in further calculations. BFGS
4 and 8, and the results of the continuous-time/discrete-spaceMC
random walks from which the self-diffusivities are derived. The procedures for the two systems (Xe/silicalite, SF6/silicalite) studied are very similar. In most cases only the results for the xenon-silicalite system will be discussed in detail. Significant differences between the two systems will be highlighted whenever they occur. The first step in the calculation is the location of all saddle points and minima of the potential energy hypersurface. Table I1 contains the coordinates and energies of the minima and saddle points identified in the asymmetric unit of the unit cell for the xenon system. It should be noted that saddle points marked as LL (low lying) were neglected in further calculations; such saddle points separate two closely spaced minima by only a small potential barrier, across which rapid thermal equilibration can occur. Thus, the three minima in the sinusoidal channel and the two minima in the straight channel (one is an image of point 10 formed by symmetry) were considered as contributing a single sorption site in each of their respective channel segments. This reasonable assumption considerably simplifies the calculation of the diffusivity by reducing the number of states appearing in the master equation describing the system. It should be noted that determining the connectivity among the minima with the steepest descent algorithm was one of the computationally most demanding tasks associated with the TST calculations because of the large number of sorbatesilicalite potential energy evaluations that had to be made. Table I1 indicates that three types of sorption sites arise from the calculation: one ( S , point 10 and its images) is associated with the straight channel, one (Z, points 11-13 and their images) with the sinusoidal channel, and one (I, point 14 and its images) with the channel intersection region. In Figure 1, the sorption states are represented as spheres and superimposed on a schematic outline of the channel system of a single unit cell of silicalite. By applying symmetry operators to the locations of the minima in the asymmetric unit, one finds that only 12 (4 of each type) of the sites are unique to a unit cell. Similarly, while only six unique transition paths between states exist within the asymmetric unit, three additional transition paths can be identified by application of symmetry operators to the potential hypersurface. The three additional paths are associated with states at least partially located outside the asymmetric unit. The number of topologically distinct transition paths remains equal to 6, however. Figure 2 demonstrates the connectivity of the three types of states in the xenonsilicalite system. All transitions between adjacent pairs of states are shown in the vicinity of the channel intersection region.
Xenon and s F 6 Diffusion in Silicalite
The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 8813 Z
Intersection
Sinusoidal Intersection
Straight Channel
Straight
Figure 2. Spatial layout and connectivity of the states located around the channel intersection. The figure is a view of unit cell down the Z crystallographic axis. Connectivity between the states is indicated by the labeled arrows. The convention for labeling the arrows is identical with that used in Table 11. In other words, arrow 1 represents a reaction path going through saddle point 1, and so on. Values for the time constants associated with the transitions shown in the diagram are given in Table IV.
Figure 4. Shape of the three sorption states for xenon in silicalite. The states are constructed by the tessellation algorithm described in the text. TABLE III: Equilibrium Mole Fractions for Each of the Three States Identified in Figure 3 for SFs at Temperatures between 100 and 300 Ka temperature, K intersection sinusoidal straight 0.779 0.162 100 0.060 0.402 0.228 200 0.370 0.283 0.204 300 0.5 13
Y
A
OThe distribution of molecules among the three states is extremely temperature dependent.
Sinusoidal Intersection Straight
0
Figure 3. Connectivity between states in the SF6/silicalite system. All transitions pass through the intersection state. The connectivity between sorption states is much simpler than in the xenon-silicalite system.
It is remarkable that, in addition to the intersection statestraight-channel state and intersection state-sinusoidal-channel state transitions, there are direct transitions between channel segments (transition I, 2, and 5) which circumvent the channel intersection state. In Figure 3 the connectivity of the states in the sF6 system is represented in a similar fashion. The connectivity between states in the SF6 system is much simpler, as only one pathway between each pair of states exists. All paths between states in the SF6 system must pass through the intersection state. Each of the three types of sorption sites is surrounded by a region in the pore that contains at least one saddle point on its surface. This volume defines the state associated with the sorption site. As described in the previous section, each of these regions was defined by a large number of discrete volume elements. With a voxel volume of 0.168 A3, the computed volumes for the sinusoidal, straight, and intersection states were determined to be 87.5, 65.4, and 19.8 A3,respectively. A rough three-dimensional picture of the sorption states for the Xe/silicalite system is shown in Figure 4. The computed volumes of the channel states in the SF6system were smaller, as there exists much less accessible pore volume for the larger sorbate molecule. A volume of 24.8, 18.4, and 48.6 A3was computed for the sinusoidal,straight, and intersection states in the SF6 system, respectively. The dependence of computed properties on the size of the volume elements was checked by investigating the computed volume of the intersection state as a function of the number of grid points used in each dimension. As seen in Figure 7, the volume of the intersection state varies by only about 5% as the total number of voxels varied by a factor
of about 30. The edge length of each of the volume elements was varied between about 0.2 and 0.75 A. It should be noted that the denominator of eq 4 is a configurational integral for a given state. The probability of occupancy of a given class of states in the unit cell is determined by normalizing the configurational integral of that state by the sum of those over all accessible states in a unit cell. In the xenon-silicalite system at 150 K, the probability of occupancy of the sinusoidal, straight, and intersection states was 0.572, 0.414, and 0.014, respectively. The intersection state, as defined by our geometrical construction, is only infrequently occupied. As temperature rises, however, the effects of energetic heterogeneity among different states subside, and the distribution of molecules among states approaches the one that would be dictated by the relative volumes of the states. While the temperature dependence of the equilibrium probability distribution is fairly weak for xenon, the distribution among states for the SF6 system is strongly dependent on temperature. Table I11 contains the computed distribution of sF6 molecules among the three states at temperatures between 100 and 300 K. It can be seen that the channel intersection becomes the most favorably occupied region at higher temperatures. While the energy minima of SF6 in the straight and sinusoidal channels are slightly deeper than that associated with the channel intersection region, the large size of the molecule greatly limits the accessible pore volume in those channel segments. As a result, at higher temperatures where energetic penalties are counteracted by entropic effects, the SF, molecule is more favorably accommodated in the spacious channel intersection region. This is demonstrated in Figures 5 and 6 where the potential of mean force has been contoured for the SF6 molecule in silicalite at temperatures of 200 and 300 K. The contouring level is such that 50% of centers of mass are included within the thick-lined surfaces. The fine-lined surfaces in the figures are contours of constant xenon-silicalite potential hypersurface and serve to define the extent of the pore network. At 200 K, few molecules (none at the specified contour level) are found in the channel intersection,
June et al.
8814 The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 Silicalite/SF60 molecu1esRI.C.. 200 K
a 9 0.0
Figure 5. Potential of mean force, -kBT In [p(r)/p], where p(r) is the singlet density distribution function and p is the mean density of sorbate present in the solid phase for an infintely dilute SF6 molecule at 200 K. The potential of mean force has been contoured at such a level that the thick-lined surfaces outline the regions where the molecule would spend 50% of its time with highest probability. The thin contour surface representing the pore network of silicalite was determined by contouring the xenonsilicalite potential energy at a level of 0 kJ/mol. Note that, while the intersection region comprises a large portion of the system volume, very few molecules are found there at 200 K due to the fact that it is energetically less favorable relative to channel regions. Silicalite/SF6 0 mo1eculesRI.C.. 300 K
~
2400.0
4800.0
7200.0
9300.0
12000.0
grid points
Figure 7. Convergence of the computed volume of the intersection state as a function of the number of voxels. Total variation in the computed volume is about 5% over a factor of 30 in the total number of voxels used in the discretization.
TABLE I V Time Constants for Interstate Transitions for Xenon in Silicate at 150 K"
no.
move (i,s,j)
1 2 3 4 5 6
(10,1,12) (10,2,13) (14,3,10) (14,4,12) (13,5,12) (146,131
TST
'Ti!
2516 1179 7.64 33.8 1027 67.3
rPT 3505 1626 224 1381 1027 2722
rT,
"The time constant, presented here is defined as the inverse of the rate constant, k k ' , in picoseconds. The six topologically distinct transitions are identified by the minima that they connect. Forward rate constants represent the move between minima as defined by the triplet as written.
Figure 6. In this figure, the potential of mean force has again been contoured at a level of 50% for SF6, but at a temperature of 300 K. The pore network contour is as in Figure 5. At 300 K, the energetic penalties associated with the channel intersection are relatively less important and a large fraction of the molecules is found to reside in this region of the pore network.
while at 300 K the distribution is well formed within the channel intersection region. The inverses of evaluated at 150 K, are shown in Table IV for the xenon-silicalite system. The entries in Table IV represent the time constants characteristic of each transition. While some transitions, such as the moves between the energetically unfavored channel intersection state and the straight-channel state occur quite rapidly, most moves take place much more infrequently. The slowest process appears to be the interchange of sorbate molecules between the straight and sinusoidal channel systems through the moves denoted as 1 and 2 in Figure 2. Time
kF,
constants in the SF6system are much larger, even at temperatures as high as 300 K. Dynamical correction factors to the TST estimates of the rate constants were computed for both sorbates using the rigid lattice model, and only for xenon at one temperature with the animated lattice. The computation of trajectories in the rigid lattice was quite efficient and could be accomplished within 2 h on a workstation for all three states after implementing a fast lookup scheme for the sorbate lattice potentiaLZ8 The animated lattice model required a vastly larger computational effort, since the trajectories of the 576 lattice atoms in the model system had to be computed for each sorbate trajectory and no lookup of the sorbate lattice potential was possible. The computation of 3000 trajectories for each state required about 15 CPU hours on a Cray Y MP. The trajectory of a sorbate molecule at infinite dilution in the rigid lattice is constrained to lie on a constant-energy hypersurface and can thus never truly thermalize; nevertheless, the momentum vector is quickly decorrelated by forces from the lattice on the sorbate molecule. The molecule is effectively trapped within a state by the finite area of the exit surfaces. It was felt that trajectories computed in this way would best mimic our previous rigid lattice molecular dynamics simulations of xenon in silicalite. In the animated lattice, the xenon atoms were found to thermalize in about 5 ps by monitoring the temperature of the trajectories as a function of time. In Figure 8, the temperature of the sorbate molecules in the animated lattice simulation is plotted as a function of time elapsed since the molecule left the dividing surface surrounding the intersection state. The temperature can be seen to
Xenon and SF6 Diffusion in Silicalite
The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 8875
d
0.0
2.0
6.0
4.0
8.0
10.0
film (PI
time (PS)
Figure 8. Average temperature of xenon molecules leaving the dividing
surface that surrounds the intersection state plotted as a function of time (animated lattice simulation). The temperature rises to a maximum and decays within IO ps to a level 30 K above that of the lattice.
0.0
2.0
6.0
4.0
8.0
70.0
(Pa Figure 9. Dynamical correction factor j&t) for the transition between
the intersection state and the sinusoidal-channel state, plotted as a function of time for the xenonsilicalite system. The value of the dynamical correction factor decreases, reflecting dynamical recrossings of the ij bottleneck surface, and reaches a plateau value about 5 ps after the initiation of the calculation. return to a value that is about 30 K higher than that of the lattice within the 10-ps dynamical correction run. Two sample dynamical correction factor plots are displayed in Figures 9 and 10. The dynamical correction factor in Figure 9 starts at a value that is proportional to the total Boltzmann weighted surface area located between the intersection and sinusoidal states, Le., the bottleneck associated with the transition labeled as 4 in Figure 2. By multiplying the initial valuejj(O+) by k z T , one obtains an estimate of the TST rate constant for hopping between the two states, As trajectories recross the dividing surface, the value of the dynamical correction factor quickly drops below its initial value,fi(O+). In a period of a few picoseconds, fd,(t) reaches a plateau region which displays a characteristic exponential decay with a time constant determined by k i i ' . It is this plateau value that we use as a dynamical
kr.
Figure 10. Dynamical correction factor characterizinga multistatejump
between two straight-channel states. Two such transitions were identified for each straight-channelstate, indicating that about 20%of the sorbate molecules exiting the state do so by jumping directly past the intersection into the next straight-channelstate. correction factor in our calculations. Figure 10 represents a multistate transition not derivable from uncorrected transition-state theory calculations. While the dynamical corrections for the straight-channel state (point 10 in Table 11) were computed, it was observed that a significant fraction of the trajectories passing into the intersection state failed to thermalize and continued directly into the next straight channel. The formulation of Voter and Doll4 allows for this type of multistate jump by computing what fraction of the molecules exiting the straight-channel state make this jump. A rate constant for the direct passage of a sorbate molecule from the straightchannel segment into a straight-channel segment past the intersection is then added to the master equation that describes the dynamical evolution of the system. This type of multistate jump lowers the rate of hopping between the neighboring straightchannel state and intersection states, but can actually increase the overall diffusivity by allowing the computed fraction of jumps to have a longer displacement associated with the direct transition between two straight-channel states. In Figure 10, the dynamical correction factor for jumps between two straight-channel segments is plotted as a function of time in the xenon-silicalite system a t 150 K. As the two straight-channel segments share no common surface, there is a 2-3-ps delay before a significant fraction of the sorbate molecules begin to arrive in the neighboring straight-channel segment. The dynamical correction factor starts at an initial value of zero and rises to a plateau value in this case. Dynamical correction factors computed for xenon in the animated lattice showed no qualitative differences from those computed for the rigid lattice, but were consistently smaller in magnitude for any given transition. As is discussed below, diffusivities computed for the rigid and animated lattices at 100 K differed by about 40%, the one in the animated lattice being lower. This is in qualitative disagreement with the molecular dynamics simulations of Demontis et a1.,I9 who studied the diffusivity of methane in both rigid and animated lattices at 298 K. In their simulations, the diffusivity of methane was observed to increase by about 20% in an animated lattice when compared to simulations performed in a rigid lattice. The dynamical correction calculations performed here on the animated lattice suggest that the weak coupling of a small sorbate molecule such as xenon with the lattice enhances recrossings of the dividing surface. On the contrary, the molecular dynamics simulations of Demontis et ai. suggested that the rate constants for crossing the dividing surface actually
June et al.
88’16 The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 TABLE V Diffusivities Computed from Transition-State Tbeory and Molecular Dynamics Are Compared for Xenon and SF, in the Rigid Lattice Model” T, K system DTST DTSTlDC DmIoct DMD IO0 SF6 0.0003 0.0004 0.066 0.057 200 0.36 0.26 300 0.046 0.033 100 Xe 0.071 0.51 0.76 0.44 150 1.5 2.7 1.1 200 run time, min 20 100 1000
a
I
t
/ / /
/
I t
1 1
/
/
/
“The column labeled as DnTIDCt contains the single dynamically corrected TST estimate of the diffusivity from the animated lattice model. Results from molecular dynamics are only available at 150 and 200 K. Molecular dynamics and the dynamically corrected TST calculations agree remarkably well. Run times are for runs conducted at a temperature of 150 K on a Sun 4/330 workstation. All diffusivities are reported in units of
/
i
/ / / / /
cm2/s. /
increased in the presence of lattice motion, perhaps through cooperative motion between the sorbate molecule and the lattice. It is significant that both techniques predict small differences in sorbate mobility due to lattice vibrations for small molecules. These effects are expected to be much greater for large sorbate molecules and molecules with internal degrees of freedom, such as alkanes or aromatics, whose motion might tightly couple with those of the lattice. As previously discussed, diffusivity tensors were computed from mean-square displacements accumulated during the course of continuous-time/discrete-spacedynamic Monte Carlo simulation runs. Additionally, average time steps were computed by dividing the total simuation time by the number of executed moves. This average time step size varied with temperature; it was as large as 0.68 ps a t 100 K and as low as 0.02 ps a t 200 K for xenon. In all cases, the total simulation time considerably exceeded the longest time constant in the system, such as the times listed in Table IV. Furthermore, the continuous-time/dismete-spaceMonte Carlo technique allows probing exceedingly large displacements. At 100 K, the average root-mean-square displacement was 200 A, or a distance of roughly IO unit cell lengths. The simulation time required to traverse this distance was 1.5 ps, much beyond the capacity of traditional molecular dynamics simulations. A typical mean-square displacement plot is shown in Figure 11 for a simulation of xenon at 100 K. Differences in sorbate mobility leading to anisotropic behavior of the diffusion tensor can be clearly seen along each of the coordinate axes. The trace of the diffusivity tensor is given in Table V for a series of temperatures. Comparison is made with molecular dynamics simulations only at temperatures of 150 and 200 K, where such simulations are feasible. The diffusivities computed with transition-state theory compare quite well with those determined from molecular dynamics, which is regarded as the “exactn numerical solution to the dynamics of the model system at high temperatures. Table V includes estimates of the diffusivity at 100 K for the rigid and animated lattices. While the diffusivity for the animated lattice is somewhat smaller than that computed in the rigid lattice, the overall magnitude is quite similar and indicates that such a complicated model is unnecessary for small, simple sorbate molecules, such as xenon. Significant attention should be drawn to the required computer times listed in Table V. As the system temperature decreases, molecular dynamics simulation becomes less efficient; this is because the system evolution, as monitored by MD, slows down dramatically as potential energy barriers between states in the system become large compared to kBT. The computational cost of MD simulations increases rapidly at such temperatures because most of the simulation spends itself monitoring only the relatively uninteresting “rattling” motion of the sorbate molecules within the energetically favorable states and seldom has the chance to observe infrequent transitions across the energetically unfavorable bottleneck surfaces which pose the actual resistance to diffusive motion. The computational cost for determining diffusivities from TST is essentially independent of
________________-__--------
d
0.0
0.3
0.6
0.9
1.2
1.5
time (nlimeanlds)
Figure 11. Mean-square displacement as a function of time from a continuous-time/discrete-spaceMonte Carlo random walk simulation of xenon molecules at 100 K. Root-mean-squaredisplacements are on the order of 160 A over time scales as long as microseconds. The rate constants for hopping in this particular simulation were taken from dynamically corrected TST results.
temperature in these systems. It should be noted, however, that the cost of the TST computations grows rapidly with the total number of degrees of freedom required in the evaluation of rate constants through eq 4. For systems where a large number of degrees of freedom are involved in the reaction coordinate, only the harmonic approximation” or other approximations3to TST may be feasible. While off-diagonal elements of the diffusivity tensor are quite small and fluctuate about zero, the anisotropic nature of the silicalite pore structure is evidenced by the nonuniformity of its diagonal elements. This difference is greatest at low temperatures, where small differences in activation energy between the channel systems are amplified. For the xenonsilicalite system at 150 K, the diagonal elements of the diffusivity tensor were computed from transition-state theory to be 1.0 X 1.2 X and 0.17 X cm*/s for the Ox, Dy,and 0, components of the diffusion tensor, respectively. At the same state point, the elements of the diffusivity tensor were determined to be 0.51 X 0.73 X 10” and 0.083 X cm2/s by dynamically corrected TST. Clearly, the multistate jump implemented between adjacent straightchannel states in the Y direction in the dynamically corrected TST calculations significantly enhances the anisotropy of the diffusivity tensor. Reasonable agreement is obtained with molecular dynamics calculations a t this temperature, which gave diagonal elements of 0.43 X 1.0 X and 0.099 X 10” cm2/s. Mobility in the Z direction is greatly reduced relative to the other elments in all cases because of the tortuous path required for motion in that direction. Comparison of our simulation results can be made with a simple random walk model developed by Urger, who investigated various bounds on the ratios of elements of the diffusivity tensor on the basis of the geometry of the silicalite lattice.45 The model was designed to analyze recent NMR and simulation results in silicalite. A lower limit was identified on the ratio:
describing the mobility of a sorbate molecule in the sinusoidal and straight channels relative to that in the 2 direction, along which (45) Karger, J. J . Phys. Chem., in press.
Xenon and SF6 Diffusion in Silicalite
The Journal of Physical Chemistry, Vol. 95, No. 22, 1991 8877
TABLE VI: Experimental Sorption Data of Reed’ Are Compared with Computed Results Based on the X e n o d x y g e n Parameters Utilized throughout This Work’
0
6
*
0
experiment simulation
3.41 1.71
26.5 24.0
I
‘The computed Henry’s constants, Kh. and isosteric heats, Q,,,derived for direct evaluation of the configurational integral,16are seen to be in good agreement with experiment. no directed channel system exists. On the basis of geometrical considerations, Karger found the lower limit on this ratio to be 4.4, independent of the sorbate molecule. Our results here are in agreement with this assertion, as the ratio for xenon in silicalite determined from TST calculations decreases from a value of 8.6 a t 100 K to 5.9 at 200 K. The SF6 system shows even greater anisotropy over a wider range of temperatures; the value of the ratio f determined from TST computations decreases from of 22.2 at 100 K to 5.5 at 300 K. The computed activation energies are quite comparable between all methods employed. The activation energy for xenon in silicalite computed from simple TST is 6.0 kJ/mol; dynamically corrected TST gives a value of 5.2 kJ/mol, and molecular dynamics gives a value of 5.5 kJ/mol. Experimental values have been recently reported for xenonsilicalite as 15 and 26 kJ/mol by two different experimental method^.)^"^ These extremely large activation energies obtained experimentally for the xenon-silicalite system (larger even than that of benzene in the same zeolite) are intriguing. The value reported in ref 39 was attributed to adsorption of the xenon atoms in defect sites in the crystal at low loadings, with the caveat that this is not an entirely satisfactory explanation. One could attribute the large disparity between our model prediction of the activation energy and the experimental values to deficiencies in the potential representation we employ, which neglects the charge distribution in the lattice and any kind of lattice defects. It should be noted, however, that our potential parametrization predicts the equilibrium sorption properties reported in ref 39 quite well. To verify that the potential parameters produced a reasonable representation of experimental sorption data, we computed a Henry’s constant and isosteric heat for xenon in silicalite2* and compared them to the data of Rees at 195 and 273 K. This comparison between experimental and calculated low-occupancy sorption thermodynamics results is presented in Table VI. While the computed Henry’s constant deviates by about a factor of 2, the isosteric heat is well represented and only deviates by 2.5 kJ/mol. Furthermore, diffusivities predicted by molecular dynamics using identical potentials to those employed here compare quite well with self-diffusivities at 300 K measured by pulsed-field gradient NMR t e c h n i q ~ e s .It~ should ~ be stressed that the Lennard-Jones interaction parameters we use here are exactly the same as in our previous molecular dynamics work, with no modification.28 A number of additional comparisons can be made between the results computed from the transition-state theory calculations and the more conventional simulation techniques we used such as Metropolis Monte Carlo and molecular dynamics. These comparisons allow a check on the consistency of each of the simulation methods. The use of low, but finite, loadings in the molecular dynamics simulations avoids problems associated with thermalization of the sorbate molecules by allowing for energy exchange between colliding molecules.M The loading necessary to produce correctly thermalized trajectories both at the potential energy wells and at the transition states was investigated by directly probing the velocity distribution of the sorbate molecules in two separate regions of the unit cell during a molecular dynamics simulation. By use of the algorithm developed to compute the TST rate constants, the velocity vectors of the sorbate molecules were binned on the basis of whether the sorbate molecule was located on the dividing surface between states or within the defined states. The distribution of kinetic energies was then normalized so the results could be compared with a Maxwell-Boltzmann distribution. The
0.0
2.0
4.0
8.0
8.0
10.0
vP(AngstromP/ p s P )
Figure 12. Distribution of the velocities from a molecular dynamics simulation in the xenonsilicalite system for molecules located on or near the dividing surface between states. The plot is almost linear, indicating that the distribution of kinetic energies among sorbate molecules conforms to that prescribed by the Maxwell-Boltzmann distribution. The relatively high noise level is due to the high potential energy at the dividing surface region, as a result of which only a few sorbate molecules are sampled in that region.
t
* .-
0
0.0
2.0
4.0
8.0
8.0
10.0
vP (AngstromP/ p s P )
Figure 13. Distribution of kinetic energies for xenon molecules located
in regions not associated with the dividing surface between states. The sampled molecules are complementary to those in Figure 12. The sampled velocities closely conform to the Maxwell-Boltzmann distribution at all energies.
logarithm of the distributions accumulated in this way is plotted against the square of the magnitude of the velocity vector in Figures 12 and 13 for xenon in silicalite a t 150 K a t a loading of 1 molecule/unit cell. In Figure 12, the distribution was gathered only over the voxels associated with the dividing surface. Figure 13 is complementary and represents all molecules within the defined states. Some deviation of the distribution from linearity is noted at high velocities for the molecules located on the dividing surface in Figure 12, but in general, the distribution of velocities conforms quite closely to a Maxwell-Boltzmann distribution. The fact that both plots are linear with slope Bm/2 indicates that the velocity distribution is truly Maxwellian throughout the zeolite and that the molecules are correctly thermalized on the dividing
8878 The Journal of Physical Chemistry, Vol. 95, No. 22, 1991
June et al. TABLE VII: Equilibrium Probability of Occupancy for Each of the Three States of Xenon at 150 rad 300 K ' method T, K intersection sinusoidal straight MD 150 0.020f 0.006 0.506 & 0.051 0.475 i 0.050 TST 0.014 f 0.000 0.572 f 0.001 0.414 0.008 MC' 0.013 f 0.002 0.580f 0.035 0.407 f 0.035 MCt 0.019 f 0.004 0.572 f 0.035 0.409 f 0.033 MD 300 0.063 f 0.003 0.530 f 0.039 0.407 f 0.039 0.057 f 0.001 0.530f 0.01 1 0.413 f 0.008 TST MC' 0.059 f 0.002 0.533 f 0.017 0.408 & 0.015
"The results computed from TST (directly from the denominator of eq 4) are at infinite dilution, while the MD simulations were carried out at a loading of 1 molecule/unit cell. Metropolis Monte Carlo results are included at both 0 (ghost molecules) and 1 molecule/unit cell to demonstrate the weak occupancy dependence of the distribution at low loadings. MC' refers to 0 molecule/unit cell Monte Carlo results while MCt refers to 1 molecule/unit cell results. Molecular dynamics results are inconsistent with those determined by the other two simulation techniques at the lower temperature. 0
c:
0.0
3.0
6.0
9.0
12.0
15.0
r (Angstmms)
Figure 14. Radial distribution function for xenon in silicalite at 150 K and a loading of 1 molecule/unit cell. The first peak in the RDF is due to the formation or 'dimers" resulting from Lennard-Jones interactions
between neighboring sorbate molecules. Other peaks are due to correlations introduced by the potential energy minima in the zeolite lattice.
surface at loadings as low as 1 molecule/unit cell. Simulations performed at lower sorbate loadings displayed significant deviations from the expected velocity distribution. To validate the assumption that, at low, but finite loadings, the sorbate molecules diffuse as monomers and not dimers, the radial distribution function of the sorbate molecules was investigated by Metropolis Monte Carlo simulations. This is an important point because the potential hypersurface of a diffusing dimer could be considerably different than that of a single sorbate molecule, which was used in the computation of the rate constants. Figure 14 displays a radial distribution function for a MC simulation performed at 1 SO K and a loading of 1 molecule/unit cell. A first coordination shell can be clearly seen at a separation of about 4.5 A. Other correlations are seen at greater separations, but those peaks are due to the periodic placement of potential minima in the unit cell of the zeolite.46 By integrating the radial distribution function, g(r), to its first minimum, the number of nearest neighbors within a distance r can be determined: n(r) = 4apL'g(r)r2 dr where p is the density of sorbate molecules in the simulation. At the state point studied, this integration indicates that only one sorbate molecule in 10 exists as a dimer pair. All three simulation methods predict an equilibrium occupancy in each of three types of states. Shown in Table VI1 are the predicted occupancies by each method for each of the three states defined in the TST calculations for xenon in silicalite at a temperature of 150 K. Reported errors reflect 95% confidence levels. The occupancies have been computed at infinite dilution by TST and Metropolis Monte Carlo, but at a finite loading of 1 molecule/unit cell by molecular dynamics. The weak occupancy dependence of the distribution of molecules among states is demonstrated by the Metropolis Monte Carlo results included for loadings of 0 and 1 molecule/unit cell. The computations show that probabilities of occupancy computed by TST and Metropolis Monte Carlo are consistent with respect to the distribution of sorbate molecules between states. Those obtained with molecular (46) Snurr, R.Q.;June, Simularion, in press.
R. L.; Bell, A. T.; Theodorou, D. N . Molecular
dynamics, however, are not consistent with MC and TST results at this temperature (1 50 K) and do not exhibit a consistent trend with temperature or loading in other calculations not reported here. To further investigate this apparent anomaly in the state occupancy results derived by molecular dynamics, we carried out investigative simulations with all three methods at 300 K where molecular dynamics is expected to be more efficient than at 150 K. We chose to look at the higher temperature to investigate whether the observed anomalies at 150 K are associated with equilibration problems inherent to molecular dynamics at low temperature or to inherent differences between the molecular dynamics NVE ensemble and the Monte Carlo NVT ensemble. As seen in Table VII, the 300 K results from the molecular dynamics simulations are consistent with those derived from TST and MC, indicating that the distribution of sorbate molecules between states is a slowly evolving quantity only at the lower temperatures investigated here. Clearly, MD runs longer than 6 ns would be required to achieve agreement at temperatures lower than 150 K.
Conclusions An application of transition-state theory with and without dynamical corrections to the computation of the diffusion coefficients of sorbate molecules in zeolites without well-defined sorption sites has been presented. Results from the xenonsilicalite system show that this approach compares very favorably with the computationally more expensive molecular dynamics method and with experiment. Furthermore, it was demonstrated that the effects of thermal motion of the animated lattice on the diffusivity of small sorbate molecules are small. Unresolved, however, remains the large disagreement between the computed and experimental activation energies for the diffusion of xenon in silicalite. The SF,/silicalite system demonstrated interesting phenomena such as the temperature dependence of the anisotropy of the diffusion tensor and the strong temperature dependence in the distribution of molecules among states. The TST-based approach is expected to be even more efficient for systems that evolve more slowly, such as alkanes or aromatics in silicalite or HZSM-5. The fine detail of molecular dynamics, which produces useful short time scale information, has been complemented here with special techniques for the analysis of infrequent events to study long time processes in a reliable and practically useful manner. Acknowledgment. Support for this work was provided by the Department of Energy under Contract DE-AC03-765F00098 and by W. R. Grace and Co.-Connecticut. D.N.T. further acknowledges the National Science Foundation for a Presidential Young Investigator Award (DMR-8857659). Finally, all of the authors thank the San Diego Supercomputer Center for a generous allocation of computer time on the Cray Y-MP/864. Registry No. SF,, 2551-62-4; Xe, 7440-63-3.
