J. Phys. Chem. 1992, 96, 1051-1060 temperature dependence of vibrational dephasing for the v1 mode below 100 K excludes phonon emission as a reasonable mechanism9J0and implies that the low-temperature dephasing rate at ambient pressure is due to isotopic scattering or inhomogeneous strain. Realistic intermolecular potential functions evaluated in the full crystal symmetry have been used to model the dynamics of molecular crystals with some success.16 A comparison of density-dependent experimentalmeasurements of vibrational relaxation with such models should allow the refinement of atomatom
IO51
potential functions to evaluate the importance of higher order electrostatic interactions and orientational interactions.
Acknowledgment. We thank Dr.David S.Moore for teaching us the basics about diamond anvil cell spectroscopic techniques. We also acknowledge the National Science Foundation (Grant CHE-9008551), the Los Alamos Center for Nonlinear Studies INCOR Program, and the donors of the Petroleum Research Fund, administered by the American Chemical Society, for financial support of this research.
Molecular Dynamics Studies of Butane and Hexane in Siiicalite R. Larry June, Alexis T.Bell, and Doros N.Theodorou* Department of Chemical Engineering, University of California, Berkeley, California 94720-9989 (Received: June 20, 1991)
Molecular dynamics simulations are performed for n-butane and n-hexane in the zeolite silicalite over a range of sorbate loadings at 300 K. The simulations show that the molecular conformation of butane is sensitive to its local environment, while hexane is found to remain in a relatively linear conformation at all times. Furthermore, the molecules are found to reside in the two channel segments with 95% or greater probability while avoiding the channel intersection region. The rate of isomerization between conformers is largely independent of position within the zeolite pore structure. Diffusivities of the sorbate molecules are computed to be monotonically decreasing functions of sorbate loading. The computed self-diffusivities are near the largest experimental values, which scatter over many orders of magnitude. The long time constants governing the decay of orientational correlation functions in the hexane system suggest the need for simulations in excess of lo00 p to properly probe the long time dynamics of the sorbate molecules.
Introduction The study of sorbate molecules by molecular dynamics in confining pores, such as those of zeolites, holds the promise to yield interesting information concerning the equilibrium and dynamical properties of molecules in constrained environments. These properties are of interest because zeolites, such as ZSM-5, are used for a number of industrially important processes involving the processing and upgrading of hydrocarbon streams. Understanding the behavior of the sorbate molecules in these environments should aid in elucidating the origin of shape selective catalytic reactions. Previous simulations of sorbate molecules in zeolites have dealt principally with small molecules such as methane and xenon with only a few notable exceptions. The reasons for this are threefold. First, there exists a wealth of complementary experimental data amcaningboth the dynamical and equilibrium properties for small molecules such as methane and xenon in zeolites from techniques such as pulsed field gradient nuclear magnetic resonance (PFG NMR) spectracopy and gravimetric sorption uptake experiments. Second, simulations involving large molecules with many internal degrees of freedom are inherently more complex computationally and require the use of smaller time steps in molecular dynamics or smaller random displacements in Monte Carlo. With the availability of faster computational resources, simulations of such molecules are only now becoming commonplace in the liquid phase. Finally, the temporal evolution of large sorbate molecules in a zeolite such as silicalite is very slow due to the relatively large potential energy barriers that exist between adjacent channel segments. Simulations in such systems are likely to be nonergodic, unless care is taken to investigate the properties of interest over sufficiently long time scales to ensure that the relevant portions of phase space have been thoroughly explored by the simulation algorithm. The sorption and diffusion of sorbate molecules in silicalite has been studied extensively because of the zeolite’s interesting Author to whom correspondence should be addressed.
0022-3654/92/2096-1051$03.00/0
three-dimensional pore network. Moreover, since the charge distribution in silicalite is largely delocalized, electrostatic interactions between sorbate molecules and the host lattice can be neglected in many cases. In past studies, methane sorbed in silicalite has been investigated using a variety of models to repesent both the sorbate and sorbent by methods such as codigurational integral evaluationI1JMonte Carlo,3s and molecular dynamics.511 Other studies of small molecules in silicalite include an investigation of xenon and ethylene diffusion by molecular dynamia.llJz Simulation studies of large molecules in silicalite are practically nonexistent, but aromatics have been investigated by energy minimization techniques to determine activation energies for diffusion and to locate sorption sites in the pore structure of the zeolite.13J4 Preliminary results from a molecular dynamics study (1) Kiselev, A. V.; Lopatkin, A. A,; Shulga, A. A. Zeolffes 1985, 5, 261-267. (2) June, R. L.; Bell, A. T.; Theodorou, D. N . J. Phys. Chem. 1990,94, 1508-1516. (3) Smit, B.; den Ouden, C. J. J. J . Phys. Chem. 1988, 92,7169-7171. (4) Snurr, R. Q.;June, R. L.; Bell, A. T.; Theodorou, D. N. Mol. Simul., in press. (5) Goodbody, S. J.; Watanabe, K.; MacGowan, D.; Walton, J. P. R. B.; Quirke, N. J . Chem. Soc., Faraday Trans. 1991,87, 1951-1958. (6) Trouw, F. R.; Iton, L. E. Zeolitesfor the Nineties, Recenf Research Reports, 1989 International Zeolites Conference, Amsterdam, The Netherlands, June 1989; Jansen, J. C., Mcwcou, L., Post, M. F. M., Eds.;p 309-310. (7) Nowak, A. J.; den Ouden, C. J. J.; PicLett, S. D.; Smit, B.; Chcetham, A. K.; Post, M.F. M.; Thomas, J. M. J. Phys. Chem. 1991, 95, 848-854. (8) Demontis, P.; Fois, E. S.; Suffritti, G. B.; Quartieri, S. J. Phys. Chem. 1990, 94,4329-4334. (9) den Ouden, C. J. J.; Smit, B.; Weilers, A. F. H.; Jackson, R. A.; Nowak, A. K. Mol. Simul. 1989,4, 121. (10) June, R. L.; Bell, A. T.; Theodorou, D. N . J . Phys. Chem. 1990.94, 8242-8240. (11) Catlow, C. R. A.; Freeman, C. M.; Vessal, B.; Tomlinson, S. M.; Leslie, M. J. Chem. Soc., Faraday Trans. 1991,87, 1947-1950. (12) Pickett, S. D.; Nowak, A. K.; Thomas,J. M.; Peterson, B. K.; Swift, J. F. P.; Chatham, A. K.; Den Oudcn, C. J.; Smit, B.; Post, M. F. M. J. Phys. Chem. 1990. 94. 1233-1236. (13) Nowak,’AyK.; Cheetham, A. K.; Pickett, S. D.; Ramdas, S. Mol. Simul. 1987, 1, 61-77. (14) Pickett, S. D.; Nowak, A. K.; Cheetham, A. K.; Thomas,J. M. Mol. Simul. 1989, 3, 353-360.
0 1992 American Chemical Society
1052 The Journal of Physical Chemistry, Vol. 96, No. 3, 1992
of butane in silicalite have been presented in ref 5 , and results from short molecular dynamics simulations have been presented for alkanes as large as propane in ref 12. In this study, united atom models of butane and hexane are investigated in the zeolite silicalite by molecular dynamics. The dynamical properties determined include the self-diffusivity of the sorbate molecules, time constants associated with conformational isomerization about nontenninal carbonarbon bonds, and the rotational deoomlation of the sorbate molecules. Equilibrium properties investigated include the spatial distribution of sorbate molecules among the channel segments and perturbations to the conformation of the sorbate molecules due to the presence of the zeolite. All properties have beem investigated over sorbate loadings of up to 8 molecules per unit cell. As the temporal evolution of these alkanesilicalite systems is very slow, the equilibration of the simulations has been checked where possible.
Model The zeolite lattice model used in this study is identical to that used in our previous simulation^.^*^*^^ The silicalite lattice is isostructural with ZSM-5, but contains no aluminum, and, in its orthorhombic form (Pnma), has unit cell parameters of 20.07, 19.92, and 13.42 A.1s The lattice contains two distinct channel systems; the first, known as the sinusoidal or zig-zag channel, is directed along the [1001 direction while the second channel system, known as the straight channel, is directed along the [OlO] crystal axis. An overall three-dimensional connectivity of the channel system results from periodic intersections between the two channel systems, allowing for motion of the sorbate molecules in the [001] direction. A phase transition which takes the crystal lattice from monoclinic to orthorhombic symmetry is known to exist.1618 For computational simplicity, the model lattice is assumed to be rigid, in the orthorhombic phase, and to have perfect crystallinity. The model sorbate molecules, n-butane and n-hexane, consist of the appropriate number of united atoms representing both the methyl and methylene groups located along the backbone of the alkane molecule. The same potential parameters represent both the methyl and methylene groups, although each of the groups is assigned the correct mass, 15 and 14 g/mol, respectively. The united atom representation considerably simplifies the simulation by reducing the total number of atoms and, hence, the total number of potential evaluations which must be made at any given time step. Furthermore, the use of the united atom representation removes the high-frequency motion associated with hydrogen atoms from the simulation model and allows a longer time step to be used. However, by assuming that all mass of the methyl and methylene groups is located at the center of force on the backbone of the chain, the moment of inertia of the molecule about the longest molecular axis is considerably reduced relative to an explicit representation of all atoms.19 Bond lengths were assumed rigid and constrained to a length of 1.53 A. Bond angles between three adjacent united atoms were allowed to fluctuate under the influence of a harmonic potential expression around an equilibrium value, e,, of 112'.
June et al. TABLEk I . t c r a d s a h r d I . b . m o k e p L r P ~ t i d P u u " (Refer to Text far DcAdtku) k,/(K rad-*) 62608 es/K 72 CO/K
u~/A
1116 1462 -1578 -368
c,/K C2/K
c3/K C4lK
3156
c5/K
-3788
3.923 83.8 3.364
csz/K QZ/A
TABLE Ik Compuiroll of Heary'e
Coaetmts d hteric Heats of Butane at 393 I( Preamted as a cbeclr of the POteatirl Re~rew~tatioa Uaed in Theae d MOIB Simulrtiws2" Kh,mg/g/atom Q",,,kJ/mol experiment 1543 47 CIEt 1440 49 CIES 1612 48
Results from c o ~ @ ~ t i o n aintegral l evaluation simulations (CIE)l are subject to an uncertainty of about 2.5%. CIEt refers to Henry's constants computed with an explicit molecular representationand CIES refers to Henry's constants computed with the lumped parameter rep resentation used here. Deviations from the experimental data of Chiang and co-worker@ are approximately equal for the two model reprtscntations. (I
BellemansZ1and consists of a six-term expansion in the cosine of the dihedral angle, 4:
In the gas phase, this potential governs the relative populations of trans and gauche conformers and the rate of isomerization between them. For hexane an additional intramolecular potential term, based on the Lennard-Jones potential, was u t i l i to r e p resent the interaction of atoms separated by more than three carbon-carbon bonds:
This potential prevents the backbone of the alkane molecule from crossing over itself. A pairwise additive Lennard-Jones intermolecular potential term was computed over all atoms 011 different sorbate molecules in the simulation cell separated by a distance of less than 13 A: Npa-1 NIp. N-
YS =
Nm
c c c /=Ie,[
in1
j = l k>i
(-)"I
(q2 rijkl
rijkl
(4)
A final potential term, arising from zeolittsorbate interactions, is computed by summing a Lennard-Jones potential expression over all sorbate atoms with the oxygen atoms of the lattice that are less than 13 A away:
N p a N--2
YB= C i=I
j=l
Y2ke(Bij- e,)2
(1)
where Nam,specifies the number of segments in the sorbate specifies the number of molecules in the molecule and Nmolc simulation cell. A force constant, ke, of 520 kJ mol-' rad-2 was assumed.*O Rotation around nonterminal carbonarbon bonds was allowed under the influence of a four-body potential expression. The potential expression used was that of Ryckaert and (15) Olson, D. H.; Kokotailo, G. T.; Lawton, S.L.; Meier, W. M. J. Chem. Phys. 1981.85, 2238-2243. (16) Klinowski, J. T.; Carpenter, T. A,; Gladden, L. F. Zeolites 1987, 7 , 73-78. (17) Fyfe, C. A.; Kennedy, G. J.; Kokotailo, G. T.; Lyerla, J. R.; Fleming, W.W.J . Chem. Soc., Chem. Commun. 1985, 740-742. (18) Hay, D. G.; Jaeger, H. J . Chem. Soc., Chem. Commun. 1984,1433. (19) Leggetter, S.;Tildesley, D. J. Mol. Phys. 1989, 68, 519-546. (20) Ryckaert, J. P.; McDonald, I. R.; Klein, M. L. Mol. Phys. 1989,67, 957-979.
where N,is the total number of oxygen atoms in the simulation cell. Following previous work,l the silicon atoms of the lattice are neglected in the summations of eq 5 . The silicon atoms, located near the center of the S O , tetrahedra, are recessed from the surface of the pore and do not define the "shape" of the pore to as great an extent as the oxygen atoms do. However, the parameters describing the oxygen atoms are to be considered as only effective parameters describing the total zeolittsorbate potential hypersurface. The Lennard-Jones e and u parameters for all interacting species are given in Table I. Interaction parameters (21) Ryckaert, J. P.; Bellemans, A. Faraday Discuss. Chem. Soc. 1978, 66, 9 5- 106.
The Journal of Physical Chemistry, Vol. 96, No. 3, 1992 1053
The Molecular Dynamics of Butane and Hexane in Silicalite
TABLE EL
Simuhtion Results for Butane d H e w loading, molecules/unit cell ( Ta), K (Ysz), kl/mol 1 4 8
289.1 311.9 311.9
1 4
308.1 299.8
** 0.2 * 0.2 0.1 -66.3 * 0.2 -67.6 i 0.1
-45.6 -45.4 -45.4
(Y% W/mol
(9% kl/mol
Butane -0.48 f 0.08 -1.55 i 0.08 -3.06 0.05
2.45 2.57 2.52
* 0.09 0.16 * 0.06
(‘VT),kJ/mol
(V), kJ/mol
**
2.1 0.2 2.1 0.1 2.4 & 0.1
Hexane -0.88 t 0.2 -3.35 0.1
5.1 0.2 4.9 i 0.1
6.5 0.2 6.3 & 0.2
-1.12 f 0.02 -1.13 0.01
“Ensemble average potential energies and temperatures are reported as functions of loading for both butane and hexane. for the alkane molecules were taken from Ryckaert,21and the zeolite parameters were taken from our earlier simulation work in zeolites.2 The zeolitesorbate c and u parameters were determined by geometric and arithmetic mean combining rules, respectively. The accuracy of the new potential parameterization was verified by computing a Henry’s constant and isosteric heat2 for butane at 393 K and comparing the computed values with experimental data. Results of the computation, preaented in Table 11, show that the new potential representation is as valid as the explicit representation used in the past, and computationally more efficient. It should be stressed that it is a necessary but not sufficient condition that the potential model correctly represent the thermodynamics of a system for a realistic representation of the dynamics of the system. In this case, the neglect of the hydrogen atoms has been accounted for thermodynamically, but it may underestimate steric hindrances between hydrogen atoms and the lattice which could significantly retard the temporal evolution of the system. As in our previous work with molecular dynamics simulations,I0 the inner summation of eq 5 was pretabulated on a grid within the unit cell and interpolated at run time with the aid of a threedimensional Hermite interpolation algorithm?2 The overall accuracy of the interpolated potential is within approximately 0.1% of that of the explicit summation (itself describing an approximate and empirical potential expression), and its computation is 1-2 orders of magnitude faster than the full summation. The Hermite polynomials provide a continuous description of the first and second derivatives of Ysz across the asymmetric unit and its boundaries. Analytical derivatives of the Lennard-Jonespotential are employed in fitting the potential expression to a threedimensional mesh that requiresa resolution of about 0.2 A for the stated level of accuracy in representing the zeolitesorbate potential.
calculatiom As in our previous in~estigations,’~ the simulations were carried out in the microcanonical (NVE) ensemble by integrating Newton’s equations of motion for the constrained molecules with Gear’s fifth-order predictorcorrector algorithm.23 Constraint forces for the bond lengths were computed by the algorithm of Edberg, . ~ ~prescribed by Edberg and co-workers, Evans, and M O K ~ S SAs deviations in the bond lengths from the proper constrained values were corrected by minimizing a function formed from the sum of the squares of the deviations of the individual bond lengths with a quasi-Newton minimization algorithm.2J A similar function involving the velocities of the constrained atoms was also utilized to correct for any drift in the constrained distances. As a result of the constraints, the molecules have 2Nato, + 1 degrees of freedom. Corrections to equilibrium properties due to constraints on the phase space distribution of trajectories arising from the constant bond lengths are expected to be less than 1% so a Fixman potential was not emp10yed.I~Temperatures based on the center (22) Schultz, M. Spline Analysis; Prentice Hall: Englewood Cliffs, NJ, 1973. (23) Gear, C. W . Numerical Initio1 Value Problems in Ordinary Dqferentiol Equotionr; Prentice-Hall: Englewood Cliffs, NJ, 197 1 . (24) Edberg, R.; Evans, D. J.; Morriss, G. P. J. Chem. Phys. 1986,84, 69334939. (25) Hillstrom, K . Nonlineor Oprimizorion Routines in AMDLIB. Technical Memorandum 297; Argonne National Laboratory, Applied Mathematics Division: Argonne, IL, 1976; Subroutine o ~ e m in s AMDLIB.
of mass velocity and atomic velocity were monitored as one measure of equilibrium in the system. In addition to the pretabulation of the lattice forces described above, a Verlet neighbor list was utilized to efficientlyevaluate intermolecular interacti~ns.~ The simulations for both systems were performed with a total of 128 molecules. The simulation cell size was varied from a minimum dimension of 40.14 X 39.84 X 53.68 A for a loading of 8 molecules per unit cell to a maximum of 80.28 X 79.68 X 107.36 A at a loading of 1 molecule per unit cell. Previous investigations have shown that much smaller simulation cells allow relevant ensemble fluctuations to Butane simulationswere carried out at loadings of 1,4, and 8 molecules per unit cell, while hexane simulations were performed at loadings of 1 and 4 molecules per unit cell. Experimental isotherms suggest that the saturation loading of hexane in silicalite is about 7.9 molecules per unit 1.~11.~’ Initial configurations were generated by random insertion of the sorbate molecules under the influence of the full model potential. Up to 20 Monte Carlo passes at 300 K were utilized to relax any initial excessive steric repulsions between the sorbate molecules and the lattice before the integration of the equations of motion was started. Initial velocities were assigned to the centers of mass of the sorbate molecules from a Maxwell-Bolt” distribution. Once the integrationof the equations of motion was started, no data were recorded over an equilibration period of 40 ps, during which the velocities were scaled to achieve the desired temperature of 300 K. Excellent energy conservation, about 0.012%in 104 steps, was maintained by using an integration time step of 2 fs. During each 500 ps long production run, data were recorded to magnetic storage every 0.05 ps for later analysis. For butane, two simulation runs were performed at each loading, the second utilizing the positions and velocities from the last time step of the first simulation as the initial configuration. Overall CPU requirements were about 32 CPU seconds for every picosecond of simulation time on a Cray YMP. ReSdtS
In this section, results for equilibrium properties, such as average sorbate energies, will be prwnted first. These properties converge rather quickly if the initial conditions of the simulation are chosen C O I T ~ Y . In a second section, dynamical properties of the sorbate molecules will be investigated. The characteristic times of dynamical processes, such as the rotational decorrelation or the translational diffusion of the sorbate molecules, can all be obtained from the decay of various correlation functions describing these phenomena. To properly probe these dynamical phenomena, the molecular dynamia simulation must span time scales several times longer than the largest relaxation time in the system. We will show for the hexane system that, while equilibrium properties have converged to apparently stable values, there exist processes whose time constants are on the order of magnitude of the duration of the simulations performed here. Equilibrium. Average values of the energy quantities defined in eqs 1-5 are summarized in Table 111. The reported uncertainties are 95%confidence intervals, computed from averaging the simulation data over 10 bl0cks.3~The temperatures reported in Table 111, T,, are computed from the velocities of all atoms (26) Verlet, L.Phys. Reo. 1967, 259, 98-103. (27) Richards, R.E.; Rees, L. V. C. Langmuir 1987, 3, 335-340.
1054 The Journal of Physical Chemistry, Vol. 96, No. 3, 1992 TABLE W. Percentage of Mokcsles in Each of the Three States D e f d by the Steepest Merit Algorithm, Reported eo a F~laetioaof sorbate Lorrdiag" loading, molecules/ straight sinusoidal intersection channel channel unit cell Butane 1 4.6 f 1.5 44.9 f 3.1 50.5 i 4.2 4 4.3 f 1.3 42.4 f 3.9 53.3 i 3.5 8 4.1 f 1.1 45.6 f 1.1 50.3 f 1.0 1 4
*
Hexane
2.4 0.5 2.0 f 0.5
48.6 f 5.3 44.0 f 4.7
49.0 i 5.2 54.0 i 4.9
"Occupancy of a state is defined by the presence of the center of mass of the sorbate molecule within the bounds determined for the state. Molecules are seen to avoid the channel intersection region of the pore network and to fill the two states located in the channel segments with comparable probabilities.
June et al.
are assumed to preferentially pack into sinusoidal channel acgmcnts at loadings less than or equal to 4 molecules per unit cell and to fill both the sinusoidal and straight channel segments only at loadings between 4 and 8 molecules per unit cell. To explore this discrepancy between the packing model of Richards and Rees and the observed simulation results, the potential energy of an isolated butane molecule sorbed in the straight and sinusoidal channels was estimated by energy minimization techniq~es.2~ With the exception of rigid bond angle constraints and a latticesorbate potential cutoff of 15 A, all other potential parameters were identical to those used in the molacular dynamics simulations. Energy minima of -52.2 and -53.0 kJ/mol were determined for molecules located in the straight and sinusoidal c h a ~ & ,respectively. The observed difference in energy between the two channel systems is sufficiently small that, at 300 K, entropic factors associated with filling the two channel systems offset any penalties associated with the energetically less favorable straight channel segments. A lower system free energy is thus obtained with the random filling of the straight and sinusoidal channel systems. As can be seen in Table 111, m l i w b a t e energetics dominate the total potential energy experienced by the sorbate molecules. The highest h t e loading of 8 molecules per unit cellcorresponds to the total number of straight and sinusoidal channel segments contained within the unit cell of silicalite. The relative constancy of (V z across ) the various sorbate loadings studied indicates that the alkane molecules fill the channel segments but are not forced into the energeticallyless favorable ChaMel intersection regions as the sorbate loading is increased. Interactions between neighboring sorbate molecules are largely limited to those occurring between the terminal methyl groups, as the molecular centers of mass are located on a roughly square array of sorption sites with a spacing of 10 A formed by centers of the channel segments. Estimates of the isosteric heat of sorption at infinite dilution can be made from the molecular dynamics (MD)simulations at the lowest loadings studied here and compared with available experimental data. The isosteric heat, defined as the difference in partial molar enthalpy of the sorbate molecules in the zeolite and the gas phases, can be approximated at low pressures as2
in the sorbate molecules rather than from center of mass motion. In the course of the short equilibration runs, temperatures were rescaled in an attempt to reach the desired target temperature of 300 K. Some temperature drift occurred during the NVE production runs, however, resulting in the reported averages. A comparison between the temperature of the sorbate molecules located in different regions of the zeolite to ensure that the molecules properly sampled a Maxwell-Boltzmann velocity distribution, even in the absence of strong thermalizing effects at low loadings and in the absence of a vibrating zeolite lattice. To perform this calculation, the pore volume of the zeolite was discretized into regions, or states, associated with the straight channel segments, the sinusoidal channel segments, and the intersection region formed between the two channel systems. Temperatures were measured for molecules whose centers of mass were within each of the defined states for molecules whose centers of mass were near the dividing surfaces between the states. The algorithm used to divide the pore volume into discrete states was based on series of steepest descent paths originating from a regular array of points placed within the pore system of the silicalite unit ce11Q2*The steepest descent paths terminated in a number of minima associated with the potential energy hypersurface experienced by a single Lennard-Jones interaction site with the zeolite lattice. This constitutes a crude but practically useful discretization of the pore volume accessed by the alkane molecules. The volume element surrounding each of the original grid points was then assigned to one of the states based on the identity of the minimum in which the steepest descent trajectory terminated. The identities of the minima were determined by their location within the straight channel, sinusoidal channel, or intersection regions. Within the unit cell of silicalite, there are four states of each type. The state volumes are 87.5,65.4, and 19.8 A3for the sinusoidal, straight, and intersection states, respectively. The temperatures of the molecules within the defined states and those near the dividing surfaces between the states were found to be always within a few degrees of each other, indicating that sufficient energy transfer between sorbate molecules was oocurring even at the lowest loadings studied. The distribution of sorbate molecules among the three types of states defined by the algorithm described above is reported in Table IV for all simulations performed in this study. The molecules reside in the channel segments and largely avoid the channel intersection region which is energetically less favorable. Furthermore, the moleculas tend to reside in the sinusoidal and straight channels with about equal probability at this temperature. The observed packing of sorbate molecules at loadings of less than 4 molecules per unit cell is seen to differ from the model proposed by Richards and which was used to explain sorption isotherms and isosteric heats for a series of alkane molecules in silicalite in the vicinity of 300 K. In that model, butane and hexane
where (U)z and ( U)Grepresent an ensemble average of the total internal energy of a mole of molecules in the sorbate and gas phases, respectively. The gas phase in equilibrium with the mlite in the low occupancy region of the sorption isotherm is ideal, so that only energy contributions from YB, YT,and Ysare included in the (U)G, while Ysz is an additional and dominant term Estimates of (U),for a low-density gas contributing to (& phase were obtained by performing molecular dynamics simulations in the NVE ensemble at 300 K. The simulations, 500 ps in duration, were performed at a gas-phase density of 0.036 g/cm3 to ensure that the molecules experienced the periodic collisions necessary for the thermalization of all degrees of freedom. At thisdensity, the weak solvation effects associated with neighboring sorbate mdeculca should not a f k t the conformational or energetic properties. The computed isosteric heats for butane and hexane are 48.3 and 70.2 kJ/mol, respectively. Experimental isosteric heats for butane and hexane at 300 K are 50 and 60 kJ/mol, respectively?' While agreement with the data of Richards and Rees is reasonable, it should be noted that isosteric heats derived from other, less extensive sorption measurements for the hexane/silicalite system are significantly larger, approaching values of 75 kJ/mol.B The homogeneity of the ptential energy surface of silicalite is evident from the fact that the predicted isosteric heat of hexane is only a few W/mol less than would be estimated from that of butane by assuming that the isosteric heat is proportional to the number of skeletal carbons. While the flexibility of bond angles is known to be important to the equilibrium distribution of torsional angles,3O excursions
(28) June, R. L.; Bell, A. T.;Theodorou, D.N.J . Phys. Chem. 1991,95, 8866-8878.
(29) Flanisen, E. M.;Rennett, J. M.;G r w , R. W.;Cohen, J. P.; Patten, R. L.; Kircbner, R. M.;Smith, J. V. Nature 1978, 271, 512-516.
at = -((U)Z- (VG) + RT
(6)
The Molecular Dynamics of Butane and Hexane in Silicalite
The Journal of Physical Chemistry, Vol. 96, No. 3, 1992 1055
? ! d
2
A
d
8 d
8
d
ad s
-180.0
-lntemsction - - Straight
-108.0
Dihedral Angle (Degrees)
-Intersection - - Straight
.......".....SlnusoMel
-36.0
36.0
108.0
180.0
Dihedral Angle (Degrees)
.............. Sinusoidal
Figure 1. Torsional angle distribution for butane in silicalite at 300 K and a loading of 4 molecules per unit cell. The distribution has been subdivided into sorbate molecules whosc center of mass are located in the straight, sinusoidal, and intersection states. The larger diameter of the channel intersection region allows the sorbate molecule to populate the gauche state substantially. Molecules in the straight and sinusoidal channels are confined by the 5.5 A wide pores. TABLE V Percmtaga of Tnss/C.pcbe Conformers for Butane and Hexme In suialite at 300 IC, at a Lording of 4 Molecuks/Unit Cell# straight sinusoidal butane intersection channel channel liquid butane 56.2J43.7f 13.2 79.0/21.0& 6.1 82.1/17.9f 4.6 60.9/38.1 hexane 86.9113.1 f 7.8 85.6114.4f 5.1 76.9123.1 5.8 70130
*
'Listed values are for the central carbon-carbon bond in the three regions of the unit cell. Estimates of the trans/gauche conformer distribution fnr hiitanr s n d hrxanr in the liniiid nhaar Pnmr frnm aimiilatinn etwiirc
performed e l ~ e w h e r ewith ~ ~ *similar, ~~ but not identical, potential models.
from the equilibrium bond angle, e, during the simulations were limited to about 4'. Such small perturbations justify neglecting the metric tensor correction to equilibrium ~r0perties.I~ The distribution of torsional angles is drastically affected by the zeolite potential field, however. In a liquid-phase simulation of butane at 291 K with a similar potential representation (differing only in the use of an equilibrium bond angle of 1 0 9 . 5 O and uniform segment masses), Leggetter and Tildesley found that an equilibrium distribution contained about 62.5% of the molecules in the trans state defined by a central torsional angle, 9, in the range -r/3 to a/3. In Figure 1, the distribution of torsional angles for butane at a loading of 4 molecules per unit cell is presented. The distribution has been further subdivided into molecules whose centers of mass are located in the intersection, straight, and sinusoidal regions of the channel system, respectively. The percentages of trans and gauche isomers, integrated for each of the three pore regions, are presented in Table V. In the straight and sinusoidal regions of the channel system, the distribution of torsional angh of butane deviates significantly from that expected in the gas phase. The n a m w pom, about 5.5 A in diameter, favor the more linear trans conformer. In the channel intersection region, whose diameter is roughly 9 A, the butane molecule explores a larger conformational space. As far fewer butane molecules reside in the channel intersedon region, the uncertainty
F m 2. Torsional angle distribution for the central carbon-carbon bond in hexane at a loading of 4 molecules per unit cell. Very few molecules are found in the gauche in any region of the channel system. TABLE M: Computed Self-Diffusivities,in Units of lo4 cm2/s, Compared with the Range of Reported Experimental Results as a FuncHoa of Sorbate Loading at 300 K" 1 4 8 molecule/ molecules/ molecules/ unit cell unit cell unit cell experiment butane 3.2 i 1.3 2.4 f 0.79 0.51 f 0.49 1 X 10-'-3.7 X lo-' hexane 2.2 f 0.53 1.0 f 0.78 7 X 10-'-7.5 X lo-' 'The experimental data are at various temperatures in the vicinity of 300 K and are discussed in the text.
in the trans/gauche distribution is much larger there. The computed averages for butane in the intersection region are statistically indistinguishable from those expected in a liquid or gas phase. Hexane shows a behavior different than the more compact butane molecule. Figure 2 shows the distribution of torsional angles about the central carbon-carbon bond for hexane molecules at a loading of 4 molecules per unit cell with the same subdivision between sorption states as in Figure 1. As with butane, the torsional angle distribution greatly favors the trans state within the channels. In this case, however, the trans state is enhanced in the intersection state as well. This result can be rationalized by the fact that the hexane molecule is less well characterized by its center of mass than the butane molecule. With a length of 10.3 A, hexane is sufficiently long to span entire channel ~egments.'~Thus, while the center of mass of the molecule may be contained within the channel intersection, the terminal methyl groups tend to reside within the energetically more favorable channel segments on either side of the channel intersection with the molecule in a straight, all-trans conformation. Interestingly, the distribution of trans/gauche conformers is statistically different between the straight and sinusoidal channels for hexane. Evidently, the difference in channel architecture is large enough to affect the distribution of torsional angles for the longer hexane molecule but not for the more compact butane molecule. (32)The length of alkane molecules in the all-trans conformation was estimated as the backbone chain length plus the Lennard-Jones diameter of the terminal methyl group. The following expreuion results for the molecular length: (N,- - 1)1, sin (4,/2) + u s , where 1, is the constrained carbon-carbon bond length. a CH2 segment to account for excluded volume interactions of
(30)Toxvaerd, S.J. Chem. Phys. 1987,87,6140-6143. (31) Clarke, J. H. R.; Brown, D.Mol.Phys. 1986, 58,815-825.
1056 The Journal of Physical Chemistry, Vol. 96, No. 3, 1992
Dynamk The dynamical property of mast inkrest in the study of zeolitesorbate systems is the diffusion coefficient of the sorbate molecules. Elements of the self-diffusivity tensor were computed over time scales from mean-squared center of mass dis~lacements~~ of up to 125 ps and over time origins spaced at 5-ps intervals. Computed selfdiffusivities for the two sorbate molecules are given in Table VI as functions of occupancy. For both sorbate molecules, the diffusivities are seen to be monotonically decreasing functions of occupancy. Also included in Table VI are the ranges of values spanned by experimentaldiffusivity data determined by various techniques. For butane, the range of experimental values shown in Table VI comes from frequency response techniques at 323 K (1 X 10” cmz s) and from membrane permeation rates at 334 K (3.7 X 10- cm2/s). The frequency response methods are believed to work well with fast diffusing systems and to be free of experimental artifact^.^^.^^ Measurements by other techniques, such as membrane permeation rates,36are as much as 2 orders of magnitude lower. Similar discrepancies also exist in the experimental literature for hexane, where estimates of the uptake (303 K) diffusivity between 7 X 10” (298 K) and 7.5 X cm2/s have been r e p ~ r t e d . ~ ’ ?Clearly, ~~ there is great variability in the zeolite samples and experimental techniques used. With smaller hydrocarbon molecules such as methane or propane, PFG NMR and inelastic neutron scattering experiments give consistent estimates of the self-diffusivity of the sorbate molecule and are a direct measure of the microscopic motion of the sorbate molecules. These techniques are in favor of the upper limit of the reported data and, thus, closer to the self-diffusivitiescomputed by MD. However, the reason for these discrepanciescould fall into either of two broad categories. First, other physical processes may mask the dynamics of the sorbate diffusion during the experimental uptake experiments. Unfortunately, the diffusivities for the alkane molecules studied in this work are outside the range of values accessible by PFG NMR experiments. Sorption uptake experiments occur in the presence of a number of physical processes, some of which may be slower than sorbate self-diffusion. Examples of such processes are diffusion of gas-phase molecules through a small bed of Zeolite crystals in the experimental sorption apparatus and the diffusion of the molecules through an amorp hous layer of silicate material on the external surface of the zeolite crystals. Second, the discrepancy between our results and experimental values may stem from deficiencies in the simulation technique and simulation method. Possible candidates are the presence of defect sites in the experimentalcrystal lattice which are not represented in the model lattice, the assumption of a rigid zeolite lattice which ignores energy transfer between host and guest, an inappropriate potential parameterization in the model system, and simulations that are too short to observe the true diffusive behavior of the sorbate molecules. While it is impossible to quantify the effect of such factors as defects in the silicalite lattice, other aspects of the model assumptions can be checked. Previous molecular dynamics simulations of methane in silicalite have investigated the differences between simulations performed with rigid and vibrating lattices? Computed diffusivities differed by only 20% for the two lattice representations. Clearly, the effect of the vibrating lattice is small for a molecule as small as methane, and while it may be larger for molecules with internal degrees of freedom, it is not expected to account for the order of magnitude differences observed here. The potential parameterization employed in these simulations accurately predicts measured Henry’s constants and isosteric heats and therefore can be assumed to give a reasonable representation of the zeolitesorbate energetics and
June et al.
I ?I / /
/
/
I
(33) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, U.K., 1987. (34) Van-Den-Benin. N.: Rees. L. V. C.: Caro. J.: Biilow. M. Zeolites 1989, 9, 312-317. (35) Shen. D.; Rees. L. V. C.: Caro, J.: Biilow, M.: Zibrowius. B.: Jobic. H. J. Chem. Soe., Faraday Trans. 1990.86, 3943-3498. (36) Hayhurst, D. T.; Paravar, A. R. Zeolites 1988.8.27-29. (37) Beschmann, K.; Fuchs, S.:Riekert, L. Zeolites 1990, ZO, 798-801. (38) Wu. P.; Ma, Y. H. In Proceedings of the 6th International Zeolite Conference, Reno, NV, July 10-15, 1983; Olson, D., Bisio, A., Eds.; Butterworths: Guildford, U.K., 1984; p 251-260.
0.0
25.0
50.0
75.0
100.0
125.0
time (pkxasmnok)
-x - - Y
----z Filpw 3. Mean square displacement for butane as a function of time at 300 K and a loading of 4 molecules per unit cell. Displacements in the X and Y directions are greatest. Displacements in the Z direction, along which there is no directed channel system, are limited to about 3.4 A in periods of 125 ps.
thermodynamic properties. Finally, it should be noted that molecular dynamics simulations performed on currently available computing resources are extremely inefficient as a means of estimating diffusivities below lod cm2/s. This statement is based on the argument that sorbate molecules, on average, must traverse a distance of about 1 unit cell for long time diffusive behavior to be realized. By the Einstein relation, the CPU time required to safely estimate a diffusivity could be expressed as (7)
where Q describes the ratio of required CPU time to simulation is the required mean square displacement so that all time, (9) accessible portions of the potential surface have been sufficiently explored by the diffusing sorbate molecule, and Osis the selfdiffusivity of the sorbate molecule. A reasonable atimate for Q on current vector processors (see Model section of this work) is 32 CPU s/simulation ps. With (s)’/2 equal to 20 A, the simulation time required to safely estimate a diffusivity of 10-* cmz/s would ex& 6000 CPU hours. The simulations performed here largely satisfy this criterion. However, as discussed below, the 2 component of the diffusivity tensor is extremely small in magnitude and, in general, does not satisfy eq 7 with the above
( 3 )‘/Z. In Figure 3, the mean square displacement for butane at a loading of 4 molecules per unit cell is presented as a function of time. The relationship between the mean-squared displacement and time is linear to a very good approximation. The anisotropy in the diffusion tensor is evident from the fact that the s l o p are very different along the three Cartesian directions. The high mobility of the model system is also evident from the fact that the displacements in the X and Y directions are on the order of magnitude of the unit cell edge lengths when examined over the entire length of the 500-ps simulation. Mobility in the Z direction perpendicular to the straight and sinusoidal channels, however, is much smaller. As can be seen in Figure 3, the average root mean square displacement for a molecule in the Z direction during a 125-pstime period is 3.4 A, or less than the diameter of the channels themselves. As there is no continuous channel system directed along the [Ool]direction of the silicalite structure, motion
The Molecular Dynamics of Butane and Hexane in Silicalite TABLE W: Adsotropy in the Diffusion Teasor Characterized through tbe Ratio 0.5(Dx, + D,)/D, for Butane. Tbe Ratio Appears to Reach a Maximum in tbe Vicinity of 4 Molecules/Unit Cell 1
butane hexane
4
8
unit cell
molecule/ unit cell
molecules/
molecules/ unit cell
1.9 8.5
9.8 11
3.8
The Journal of Physical Chemistry, Vol. 96, No. 3, 1992 1057 a
-*! I
of the molecules in the Z direction can only occur through a series of passages between straight and sinusoidal channel segments. The anisotropy in the diffusion tensor is described by a ratio
f:
where D,, represents the diffusivity of the sorbate molecules in the sinusoidal channels of the crystal. A lower limit of 4.4 is p l a d on the value off by models based on a stochasticjump description of the diffusion process.3g Contained in Table VI1 are the predicted values off for three different sorbate loadings. At the highest loading studied for butane, the value offfalls below the lower limit of 4.4 for a jump model; motion of the sorbate molecules in this regime is highly correlated at the near-saturation loading and is not expected to be well described by a simple jump phenomenology. A second measure of anisotropy in the diffusion tensor is the ratio Dxx:Dy,,.This ratio compares the mobility of the sorbate molecules in the straight and sinusoidal channels. At a loading of 4 molecules per unit cell, the value of this ratio for butane and hexane is 0.32 and 0.19,respectively. The ratio is expected to decrease with increasing molecular weight as larger sorbate molecules must distort their torsional angles away from the energeticallyfavorable trans potential well to move through the tortuous sinusoidal channels. The detailed microscopic motion of the sorbate molecules can be studied through the velocity autocorrelation function (VACF) and its Fourier transform. By examination of the frequency components of the VACF, the time and length scales that make up the motions can be identified. Figure 4 presents the Fourier transform of the VACF for the center of mass motion of butane at a loading of 4 molecules per unit cell. Two distinct peaks are seen in the spectrum for the X and Ycomponent of the VACF, while only one frequency maximum is prominent in the Z component. As determined in earlier molecular dynamics studied0 of methane and xenon, the low-frequency component (0.4 ps-’) corresponds to the motion of sorbate molecules within the channel segments between the potential barriers presented by the channel intersection. The high-frequen component (1.5p-’)corresponds to a length scale of about 1-2 at 300 K and is associated with a “rattling” motion of the molecules perpendicular to the channel axis. It is interesting to note that this separation of motion into directions parallel and perpendicular to the pore axes can be directly observed by animating the molecular dynamics data on a graphics workstation. The time scales associated with motion of the internal degrees of freedom for both sorbate molecules were investigated through a series of short simulations (20 ps) at loadings of 4 molecules per unit cell and 300 K. During these simulations, data were written to disk every 0.02 ps to avoid limitations imposed on the Fourier transform process by the Nyquist criterion.@ The Fourier transforms of velocity autocorrelation functions of individual segments are presented in Figure 5a,b as a function of frequency. Distinct peaks in Figure 5a,b at frequencies greater than 3.0 ps-’ correspond to the motion of internal degrees of freedom. For both sorbate molecules, peaks are Seen at frequencies in the range 3-16 ps-I. To determine which frequency peak corresponded to a given
1.0
0.0
2.0
3.0
4.0
5.0
-x - - Y
----z Figure 4. The Fourier transform of the center of mass velocity autocorrelation function of butane at 300 K and a loading of 4 molecules per unit cell. Two distinct peaks are seen in the X and Y directions. The high-frequency peak,common to all components of the VACF, originates in a rattling motion normal to the pore axis. The low-frequency peak, present only in the X and Y directions, indicates that the molecules shuttle between channel intersections with a path length of about 7 A. a: Butane
a
51
A ’
I
frequency(invarsa picasecon&)
b: Hexane
a
1
1
(39) Kirgcr, J. J . Phys. Chem. 1991, 95, 5558-5560. (40)Press, W.H.;Flannery, B. P.; Teukolsky, S.A.; Vetterling, W . T. Numerical Recipes: The Art of Scientific Computing; Cambridge University Press: London, 1986.
0.0
4.0
6.0
12.0
16.0
20.0
traquamy (invam ~ a w n d s )
Figure 5. The Fourier transform of the segmental VACF for butane (a) and hexane (b) at 4 molecules per unit cell and 300 K. Low-frequency portions of the spectra are similar to the center of mass VACF in Figure 4. High-frequency components repre-sent torsional angle motion (3-7 ps-’) and bond angle motion (8-16 ps-I). The hexane spectrum (b) has
a much richer structure due to the increased number of degrees of freedom in the sorbate molecule. Note that 0.03 ps-’ is the equivalent to 1 cm-’.
degree of freedom, the potential energy parameters in Table I governing the torsional and bond angles were modified in two separate simulations by factors of 2. Resulting shifts in the peaks could then be assigned unequivocally to particular degrees of freedom. In this manner, peaks in the range 3-7 p-’were assigned
1058 The Journal of Physical Chemistry, Vol. 96, No. 3, 1992
June et al. TABLE MI: Time C c w t ~ t sfor tbe Decay of the End to End Vector Time Autworrehtiw~Fsnction, Pl(t)@
butane hexane Y
d
,
0.0
I
2.0
4.0
6.0
8.0
10.0
1
4
8
molecule/
molecules/ unit cell 123163.7 5311353
molecules/ unit cell 222/78.1
unit cell 115/70.9 182/155
'Relaxation times for the X and Y components of the vector are presented here in pairs (TJTJ, in units of picoseconds. The correlation function decays on time scales of hundreds of picoseconds, indicating the need for very long molecular dynamics simulations.
time (piwseconds)
-
b: long time
9 ,
Y
0.0
25.0
50.0
75.0
1
100.0
125.0
time (picaseconds)
Figure 6. hrrelation of the product (u(t).u(O)) as a function of time for butane at 300 K and a loading of 4 molecules per unit cell. The X and Y components d a y with time constants longer than 50 p, as clavly seen in Figure 6b. The slow decay in the X and Y directions is associated with the interchange of molecules between the straight and sinusoidal channels. The Z component decays with a relaxation time of lcss than 5 p. A rocking motion is exhibited by the sorbate molecules in Figure 6a at short times. The rocking motion has a characteristic frequency of
about 1.0 ps-' and is quickly damped out.
to torsional angle oscillations and those in the range 8-16 ps-' were assigned to bond angle vibrations that would be observable by IR or Raman s m m p y . The spectrum for hexane in Figure 5 is much richer given the greater number of internal degrees of freedom. For both molecules, the characteristic frequencies of center of mass translational motion are much lower than those associated with the internal degrees of freedom, suggesting that coupling between the two modes of motion would be inefficient. The observed long time scale phenomena in the sorbatezeolite systems studied here are best exemplified by the decorrelation of a vector that describes the orientation of the alkane molecule in the pore. The vector we chose to study this behavior is the normalized end to end vector, u(t), of the sorbate molecule. The decorrelation of this vector as a function of time is studied through the first Legendre polynomial, P l ( t ) : Pdt) = (U(t).U(O)) = (u,(t)-u,(O)) + (u,W.u,(O)) + (u,(?).u,(o)) (9) where u(t) is a unit vector in the direction connectingthe two ends of the molecule. As the two channel systems in the silicalite lattice are orthogonal and occupied with nearly equal probability, the average value of the correlation function, P , ( t ) , is expected to approach zero at long times. Since the alkane molecules studied here are large, they experience significant changes in orientation only when moving between channel systems. For hexane, significant decorrelation of n(?) is thought to occur when molecules pass from the sinusoidal channels to the straight channels and vice versa. It should be noted that butane is small enough that it is expected to tumble in the channel intersection region, which would accelerate the decorrelation of u(t). Thus, the relaxation times extracted from the decay of P l ( t ) provide a reasonable, but approximate, measure of the rate of interchange of sorbate molecules between the two channel systems. Figure 6 shows the three Cartesian components which make up P l ( t ) for butane in silicalite at a loading of 4 molecules per unit cell. The vector u(t) was correlated for times up to 125 ps with a gap of 5 p between data windows. The short time behavior is demonstrated in Figure 6a where a rocking motion of u(t) is
evidenced by small oscillations in the time correlation function P l ( t ) . This motion has a frequency of about 1 ps-" and damps out quickly. Significant variation is seen in the components of P l ( t ) along different directions at long times in Figure 6b. The long-time decay in the X component is mainly associated with molecules leaving the sinusoidal channels and moving into the straight channels. Similarly, the longtime decay in the Y component is mostly associated with molecules initially in the straight channels moving into the sinusoidal channels. As there is no channel system in the [001] direction, the Z component of u(t) is small and decorrelates quickly through a small amplitude rocking of the sorbate molecule. The decay in this component does not reflect any long-range motion in the pore system. The time constants listed in Table VI11 for all simulations performed in this study were derived from a least-squares fit to the logarithm of P l ( t ) for times between 25 and 125 ps. The large values of these time constants demonstrate the extremely slow rate of interchange between the two channel systems. Assuming that the zeolite-sorbate system should be sampled for at least 2 or 3 relaxation times to gather adequate statistics on interchange between the two channel systems, molecular dynamics simulations of at least 1000 ps in duration should be utilized to fully probe the long time behavior of the hexane system. This is in stark contrast to what happens in bulk monoatomic liquids often studied by molecular dynamics where all relaxation phenomena occur on time scales of picoseconds. The process of conformational isomerization is another long time scale phenomenon that is intimately linked to the interchange of sorbate molecules between the straight and sinusoidal channel systems. For hexane, whose length is greater than the diameter of the channel intersection, movement between the straight and sinusoidal channels requires that the molecule undergo conformational isomerization. On the basis of this argument, the rate of interchange between the two channel systems should not be greater than the rate of conformational isomerization. It is possible to quantify the rate of conformational isomerization processes by defining a rate constant, k. In previous investigations of conformational isomerization in the liquid state, the form of the rate expression governing the isomerization process has invariably been assumed to be first-order reaction k i n e t i c ~ . 4 ~To * ~verify ~ these assumptions, the rate of conformational isomerization of the sorbate molecules was monitored by two methods. The first method directly counts the number of trans/gauche potential barrier crossing events in a given period of time but gives no other information on the dynamics of the isomerization process.4' A barrier crossing event is defined as the movement of the torsional angle past either of the two maxima in the torsional potential which are located at fn/3. By such a counting of the number of barrier crossings, a transition-state theory (TST) estimate of the rate constant, k", is obtained.45 The overall rate constant contains contributions from the four likely isomerization events (gauchetrans, gauche+ trans, trans gauche+, and trans gauche-) between the three states. Direct transitions between the two gauche states via an eclipsed conformation are extremely unlikely, and although others4I have included an adiabatic gauche to gauche transition via the trans state in previous studies, this
-
-
-
+
(41) Brown, D.;Clarke, J. H. R. J . Chem. Phys. 1990, 92, 3062-3073. (42) Romberg, R.0.; Beme, B.J.; Chandler, D. Chem. Phys. Lerr. 1980, 75, 162-168.
The Molecular Dynamics of Butane and Hexane in Silicalite pathway is not considered here. Estimates of kTSTwere obtained as
where N E is the total number of tabulated barrier crossings, N,,, - 3 the number of nonterminal carbon-carbon bonds which may conformationallyisomerize, and t the length of the simulation run. The data were examined for barrier crossings every 0.05 ps. The TST rate constant overestimates the true rate of conformational isomerization because many bonds fail to thermalize in the destination state following a barrier crossing and quickly return to the original state. Furthermore, the formulation of eq 10 gives no information on the dynamics of the isomerization process other than in the context of the assumed first-order kinetics; Le. no information is provided to determine the true rate law of the isomerization process. It is possible to account for the short time dynamical recrossings of the potential barrier and to obtain insight into the dynamics of conformationalisomerization through a more detailed analysis of the process of conformational isomerization from long molecular dynamics simulations?' In the following discussion, we briefly review the necessary formalism presented by Brown and Clarke4' and then implement it to discuss the conformational isomerization process of alkane molecules in silicalite. As a consequence of the Onsager fluctuation-dissipation theorem, the relaxation of any spontaneous fluctuation from an equilibrium distribution is governed by the same characteristic times as the relaxation from a nonequilibrium distribution to an equilibrium di~tribution.~~ In this context, we can follow the relaxation of a population of sorbate molecules in a single conformational state to an equilibrium distribution of conformers as a function of time during a molecular dynamics simulation. The subgroup of sorbate molecules characterized by a single conformational state taken by themselves constitutes a nonequilibrium distribution which will relax into an equilibrium distribution over time. For example, the population of interest could be chosen as the subgroup of all molecules that find themselves in the trans state at time t during a molecular dynamics simulation. By the fluctuation-dissipation theorem, the dynamics by which this subgroup decays to an equilibrium distribution is representative of the rate of conformational isomerization for all molecules in an equilibrium simulation. To follow the relaxation of the molecules in the trans state, we Can define a function &[f#q(t)]: 1 if
HT[di(t)l HT[&(f)]
= 0 if
I'$i(t)l
5
(11)
where c$i(t) represents the torsional angle of molecule i at time t." Using &[&(?)I, it is then possible to define a relaxation function, Rm(t), by which the dynamics of the isomerization process can be followed: 1 "rl. RTT(t> = -( HT[h(fo)lHT[#i(fo + t)l> (12) NmoIe i=1
where the average is taken over all time origins, to. The relaxation function Rr,(t) decays from an initial value equal to the equilibrium mole fraction of trans molecules, (&), to a long time value of (XT)2.Using the short and long time vaiues of the relaxation function, a normalized relaxation function, Rm(t), can be defined
Excluding any short time effects not described by simple first-order (43) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: New York, 1987. (44)Chandler, D. J. Chem. Phys. 1978,68,2959-2970. (45) Brown, D.; Clarke, J. H. R. J. Chem. Phys. 1990, 93, 4117-4122. (46) Chiang, A. S.; Dixon, A. G.; Ma, Y. H.Chem. Eng. Sci. 1984, 39, 1451-1459, 1461-1 468.
The Journal of Physical Chemistry, Vol. 96, No. 3, 1992 1059 TABLE Ix: Rate Constant, itm, for the C O ~ O ~ ~ ~ O I U ~ Isomerization of the Central CIrbon-Carbon Bond of the Sorbate Mdecdea Studied, Presented in Unita of ns-l a
butane hexane
1
4
8
molecule/ unit cell 70* 10
molecules/ unit cell 92 10
molecules/ unit cell
75
* 15
67
** 13
49
*8
'The TST rate constant provides only an estimate of the true sorbate conformational isomerization kinetics, which are more compiicated than first order over the timescales of the simulation.
-1 p io'
