Molecular dynamics simulation of propane and methane in silicalite

J. Phys. Chem. 1993,97, 4149-4163

4149

Molecular Dynamics Simulation of Propane and Methane in Silicalite John B. Nicholas,'*+Frans R. Trouw,* John E. Mertz,B Lemox E. Itondl and Anton J. Hopfinger' Department of Chemistry, University of Illinois at Chicago, Chicago, Illinois 60680, Intense Pulsed Neutron Source Division and Materials Science Division, Argonne National Laboratory, Argonne, Illinois 60439, and Cray Research, Inc.. 655 East Lone Oak Drive, Eagan, Minnesota 55121 Received: September 16, 1992; In Final Form: January 4, 1993

We have investigated the diffusion of propane and methane in the molecular sieve silicalite by computer simulation using energy minimization and molecular dynamics techniques. We present heats of adsorption and selfdiffusion constants, calculated using four sets of nonbonded interactions, and compare them to experimental values. Extensive simulation results for large ensembles of methanes in a rigid-molecule approximation are presented. Methane is studied a t infinite dilution and a t loadings of 2 , 4 , 8, 12, and 16 molecules per unit cell. Theoretical self-diffusion constants range from 11.5 X to 2.0 X 10-5 cm*/s a t 300 K, in excellent agreement with pulsed field-gradient spin-echo nuclear magnetic resonance measurements. Simulations of propane allowed free movement of all the internal coordinates of the molecule and incorporated large ensembles to achieve accurate representations of bulk properties. Propane is studied a t infinite dilution and loadings of 4 and 12 molecules per unit cell. The corresponding theoretical self-diffusion constants are 2.3 X and 6.0 X lo-' cm2/s a t 300 K. These simulated diffusion rates are also in excellent agreement with NMR measurements. Center-of-mass time distributions were calculated and energy minimizations of the molecules within the zeolite lattice were done. These analyses show that the zigzag channels are the favored residence sites for both methane and propane. The calculated isosteric heats of adsorption of methane and propane are -5.8 and -10.3 kcal/mol, respectively, in good agreement with experimental values. Both adsorbate-silicalite and adsorbate-adsorbate interactions are shown to have an effect on the packing of the molecules in the zeolite. In addition, dynamic effects related to the anisotropy of diffusion in the silicalite lattice also have a role in packing. The microscopic diffusive behavior of the molecules is dicussed and compared to the classical jump diffusion model.

I. Introduction Molecular shape-selectivity and crystalline regularity are the distinguishing features that set apart zeolites and related molecular sieve materials from heterogeneous catalysts that have nonporous or only mesoporous internal structure. Although some zeolites exhibit external surface reactivity, the primary reactive sites are inside the cages and channels of the micropores defined by the crystalline framework. Because the channels have sizes of "molecular dimension", small organic molecules are selectively adsorbed based on their size and shape. Examples of important industrial processes that exploit shape-selectivity are the alkylation and isomerization of aromatic hydrocarbons' and the conversion of methanol into gasoline-range hydrocarbons.* The mediumpore zeolite ZSM-5 (structure type MFI) is used in both of these applications. Silicalite is the aluminum-free form of this zeolite and is essentially devoid of the charge-balancing cations and/or Bransted acid sites necessary to catalyze reactions and coordinate polar molecules. Silicalite is thus hydrophobic and organophilic and is considered to offer an energetically homogeneous surface for the adsorption of small, nonpolar organic molecules. It finds use in the removal of organics from wastewater and in the separation of mixtures of hydrocarbons. Two modes of molecular shape-selectivity have been distinguished: mass transport molecular shape-selectivity and transition state molecular shape-~electivity.~The latter depends on the steric accommodation of the reactive transition state complex at the catalytically active site in the channel or cage of the zeolite. The former relates directly to the intracrystalline Author to whom correspondence should be addressed. Current address: Molecular Science Research Institute, Battelle, Pacific Northwest Laboratories, Richland, W A 99352. University of Illinois at Chicago. t Intense Pulsed Neutron Source Division, Argonne National Laboratory. Cray Research, Inc. 11 Materials Science Division, Argonne National Laboratory. +

translational diffusion of the molecules, both for reactants and products, and how the diffusion is influenced by the strength of the adsorptive interaction with the framework and by the molecular dimensions relative to the pore geometry. Mass transport shape-selectivity is the only mode that is pertinent to silicalite because it is chemically unreactive. We have applied various methods of computational chemistry to examine the dynamics of the silicalite framework and the adsorption and translational dynamics of small hydrocarbon adsorbates on silicalite. In this paper we present the results of the studies of adsorbed propane and methane. Applicationsof computational chemistry encompass computer simulation and modeling of chemical phenomena based on the concepts of theoreticalchemistry. Molecular modeling techniques are excellent tools for the study of zeolite structure and the interaction of charge-balancing cations and adsorbed molecules with the zeolite framework. The theoretical methods can be divided into those that are based on classical mechanics (CM) and those that are based on quantum mechanics (QM). Although many QM studies have been performed on molecular fragments representing zeolites, calculations on the scale of the unit cells are only now becoming feasiblec~mputationally.~ This is because the size of QM calculations on a system of n atoms scales approximately as n3-n6, making these calculations on large systems very computation intensive. QM methods have thus been limited generally to the investigation of the localized structural and electronic properties of zeolites, e&, proton affinity of the framework5 and complexation of adsorbate at active sites: which are some of the crucial considerations in problems of catalysis and reactivity. However, the simulation of molecular and framework dynamics, which is the essence of mass transport molecular shape-selectivity, presents too large a problem for QM methods at the current state of the art. In contrast to QM, CM methods scale approximately as n2, making the simulation of macroscopic properties more tractable.

0022-365419312097-4149%04.00/0 0 1993 American Chemical Society

4150 The Journal of Physical Chemistry, Vol. 97, No. 16, 1993

For example, molecular mechanics energy minimizations (EM) can be used to explore the energetics of the zeolite framework,’ potential surfaces within the zeolite pores, and the preferred sites of interaction of counterions and adsorbed molecules with the zeolite framework. The C M calculations also allow the dynamics of these systems to be studied. These molecular dynamics (MD) simulations can be used to explore all the time-dependent properties of the system, e.g., transport coefficients, vibrational spectra, and y-ray and neutron scattering functions. Both Q M and CM are powerful techniques for providing insight into experimental results; they become particularly valuable for making predictions about systems that are difficult or impossible to investigate experimentally, e.g., determining diffusivities of adsorbates at low loading or high temperature. The results obtained with all molecular modeling methods that are based on C M depend heavily on the force field used to represent inter- and intramolecular interactions. Generally, the force field contains “valence” terms that represent bond stretches, angle bends, and rotations about dihedral angles. Also included are “nonbonded” interaction terms that represent nuclear repulsion, dispersion forces, and electrostatic interactions. Various welltested force fields for organic molecules are widely used in molecular modeling. Unfortunately, the available force field data for inorganic atoms are sparse. Because of the importance of the force field to our simulations, part of the focus of this work was to assess the quality of four different force field representations of the nonbonded interactions pertinent to organic adsorbates in zeolites. This was done by comparing with experimental data the calculated heats of adsorption and calculated diffusion constants for themethantwilicalite system. Thechoice was made on the basis of the best agreement for both of these measures, one static and the other dynamic. The cardinal assumption thereafter is that this force field is transferable to the systems involving larger hydrocarbons. After establishing the preferred force field, we proceeded to investigate some properties of adsorbats-silicalite system for which experiments are inconclusive or in disagreement. One issue we have sought to address is that of the most probable distribution of the adsorbate molecules in the silicalite pores. There have been divergent hypotheses on this issue offered in the literature; at lower loadings, n-alkanes have been assumed by some authors to preferentially occupy sites in the channel intersection^,^^^ while preferential occupation of sites in the zigzag channel segments has been assumed by others.l03l1 The difference between the two interpretations hinges on the relative importance of adsorbateadsorbate and adsorbatezeolite interactions. The MD simulation is the only plausible approach to answering this question for the system at room temperature where the adsorbate is diffusing rapidly. A second important issue is the calculation of the translational diffusion constant of the adsorbed molecules, its dependence on adsorbate loading and temperature, and the relative magnitude of its anisotropic components. The intracrystalline self-diffusion constants of methane and propane in silicalite have been measured by the microscopic technique of pulsed field-gradient spin-echo nuclear magnetic resonance (PFGSENMR)8 and by several macroscopic uptake measurements, using gravametric, volumetric, and gas chromatographic approaches. Comparisons of the microscopic and macroscopic results have been made by Karger and Ruthven.I2 It is remarkable that the macroscopic results obtained for propane by zero length column (ZLC),13frequencyresponse (FR),I4 and membrane methods,’5 are in excellent agreement with each other, but diverge widely from the microscopic value determined by the PFGSENMR method.8 For example, Eic and Ruthven have recently measured the diffusion of propane in silicalite using the ZLC method.13 For a concentration of 4 molecules per unit cell (mol/u.c.), they report a value of (7.4-8.2) X cm2/s. Their results are more than 2 orders

Nicholas et al. of magnitude slower than the PFGSENMR measurements of the same system done by Car0 and co-workers.8 In contrast, the new square-wave modulation variation of the FR macroscopic method yielded a value of 3 X cm2/s at 334 K, in good agreement with PFGSENMR results. The FR result requires that the sample be well-dispersed so that the silicalite crystallites were separated from each other.I3 The implication of the new results obtained by the squarewave modulation method is that intercrystalline diffusion effects may yet dominate macroscopic diffusion measurements, even when the methods used are ones refined so as to minimize these errors, specifically, the ZLC and FR methods. These heat- and/ or mass-transfer resistance effects are expected to be most pronounced when the true intracrystalline mobility is high, as it is for methane and propane in silicalite. One of our primary goals of the propane/silicalite simulation was to resolve this discrepancy between the results of the different experimental procedures. Given a reliable transferable force field, the best test of the accuracy of the PFGSENMR measurements is provided by the MD simulation. The PFGSENMR data show, in addition, a very strong dependence of the diffusivity on adsorbate loading,8 and this is also conveniently explored in the MD simulations. Loading is also expected to influence diffusion measurements made by the membrane method. In this method, the gas is allowed to diffuse through a single crystal membrane under the influence of a pressure gradient. It is significant that the diffusion in these membrane method experiments is parallel to the c-axis of the crystal.16 The c-axis is shown in the MD simulations to have a diffusivity component substantially smaller than those in the aand b-axis directions. In contrast, the N M R pulsed field-gradient method is expected to give an average value for the diffusion constant due to the random orientation of the silicalite crystals in the powder sample.8 It will be shown from simulations of adsorbed methane and propane that the force fields we have used for the hydrocarbons in silicalite yield theoretical diffusivities that agree very well with PFGSENMR results. There have been several previous computer simulation studies of adsorbates in ze01ites.l~~~~ The adsorption of alkanes in silicalite was studied by June et al.I9 using Monte Carlo techniques. Cohen de Lara and co-workers performed a MD study of methane in zeolite Na-Ae20 Recently, MD simulation studies have been published on the xenon/silicalite s y ~ t e m , I ~methane/silicalite,l9 .~l and water/ferrierite.22 The potential energy surface for the adsorption of methane in silicalite has been described by VigneMaeder and A ~ r o u x Since . ~ ~ the completion of the simulations reported here, Nowak et al. have published MD simulations on the mobility of methane in various zeolite structures, as well as results on the mobility of ethane and propane in silicalite.24 In this paper, we present the first complete MD simulations on hydrocarbon diffusion in zeolites for which all internal degrees of freedom of the adsorbate, propane, are treated. The importance of the internal degrees of freedom on diffusion through the micropores increases with the chain length of the linear alkanes. The assumption of rigid-molecule representations becomes untenable for molecules longer than ethane. We have also chosen to include Lennard-Jones and electrostatic interactions between all atoms in the systems, foregoing the assumption made in previously published work that interactions involving silicon atoms can be ignored. Furthermore, large ensembles of molecules (54432) have been employed, so that the estimated properties have good statistical reliability. Hence, we are presenting what should serve as the prototype for even more ambitious simulations that, in addition, will incorporate the dynamics of the zeolite framework. In anticipation of the detailed treatment of the propane adsorbate, a large body of work was done in developing the force field for the simpler problem of methane adsorbed on silicalite. Methane can be satisfactorily treated in a rigid-molecule approximation, and the relative computational tractability of the simulations

Diffusion of Propane and Methane in Silicalite

The Journal of Physical Chemistry, Vol. 97, No. 16, 1993 4151

h

.

..

a : Zig-Zag Channel

b : Straight Channel c : Combination of a and h Figure 1. Schematic diagram of the silicalite pore structure illustrating the channels and how diffusion is effected in the a, b, and c directions.

a

I

0 Silicon

0

Oxysen

Figure 3. Silicalite unit cell looking down zigzag channels (a-direction).

*OS, c = 0. Assuming an oxygen van der Waals radius of 1.35

A, the straight channels are elliptical with major and minor axes of 5.8 and 5.2 A. Figure 3 is a similar plot of silicalite viewed a

I

b C

0 Silicon

Oxygen

Figure 2. Silicalite unit cell looking down the straight channels (b-direction).

permitted their extensive use in testing the interaction potentials. Hence, we present some of this work as a significant component of the paper. We comment, as appropriate, on the important departures of our simulations from those published by others. 11. Theoretical Methods

A. Silicalite. The silicalite unit cell is orthorhombic with a = 20.044 A, b = 19.918 A, and c = 13.395 A and contains 96 silicon and 192 oxygen atoms (Si02). Atomic coordinates for the silicalite lattice are taken from Flanigen et al.2s The silicalite structurecontains a three-dimensionalnetwork of interconnecting channels (Figure 1). The channels are defined by 10-rings that allow the adsorption of molecules with a kinetic diameter of up to =6.0 A. Thus, selectivity of adsorption is based in large part on the size of the adsorbate molecule. For example, benzene, with a kinetic diameter of 4 . 8 5 A, can enter the silicalite pores. "pentane, with a larger kinetic diameter of 6.2 A, is largely excluded at room temperature. In contrast to the isotropic diffusion in gases and liquids, the nature of the silicalite channel system forces diffusion to be anisotropic. Figure 2 shows the atomic structure of the silicaliteunit cell looking along the 6-axis. Half of the straight channels are shown on each edge of the unit cell. The straight channels run parallel to the b-direction and are centered at fractional coordinates of a = 0, b = &OS, and c =

along the a-axis. The two zigzag channels run parallel to a and are at fractional coordinates of b = 0.25 and b = 0.75. Although it is not apparent from this view of the crystal structure, the zigzag channels have a morecircular cross section than the straight channels and are slightly smaller, with a radius of 5.3 A. The zigzag channels are usually depicted as very cylindrical in shape. However, plots presented later will show that the path of the molecules through the zigzag channels is more irregular than implied by the schematic in Figure 1. From a strictly geometric point of view, diffusion is expected to be faster in the larger, straight channels (6-direction) than in the zigzag channels (adirection). Long-range diffusion along the c-direction requires a "two-story" process that involves movement along a zigzag channel to an intersection, movement in a straight channel to another intersection,and, finally, further movement in the zigzag channel. Diffusion along c is therefore expected to be slowest because of the convoluted path the molecules must take and because traversing the length of any single channel segment produces a larger displacement along either a or b than along c. There are no water molecules or charge-balancing cations necessary in the present calculationsbecause the simulated lattice is neutral. This simplifies the calculations since only the adsorbate-adsorbate and adsorbate-zeolite nonbonded interactions are needed. The similarity in the van der Waals diameter of the silicalitechannels and the methane and propane molecules gives negligible short-range repulsions, leading to an overall favorableinteraction between the organic adsorbatesand silicalite. The simulated silicalite lattice is a 3 X 3 X 3 arrangement of unit cells. The total size of the periodic simulation box is thus 60.132 X 59.754 X 40.1 85 A, containing 7776 framework atoms. The large box size is needed to minimize the effects of the periodic boundary conditions, used to simulate the infinite lattice. The large box also allows the inclusion of a sufficiently large enough number of molecules in thesimulation toobtain accurate statistical averages of their calculated physical properties. The zeolite framework is the same for both the methane and propane simulations. The framework atoms are held fixed; thus the total momentum is not conserved in the MD. The framework atoms interact with the organic molecules via Lennard-Jones and electrostatic potentials (see below). B. Methane. The atomiccoordinates of methane are obtained from the tetrahedral H-C-H angle and a C-H bond length of


1.09 A. The M D runs are done within the microcanonical ensemble (NVE) at 300 K with the methane molecules held rigid. The rotational motion of the molecules is treated with the quaternion algorithm.26 To propagate the molecules through time, the translational equations of motion are integrated by a fifth order predictor-corrector scheme,z7 while the rotational motion is handled by a fourth order predictorxorrector. An integration time step of 1.5 fs was found to maintain energy conservation. Each system is equilibrated for 40 ps or more by periodically scaling the rotational and translational velocities toward the desired temperature. The average rotational and translational temperatures agreed to within =l-5 K, suggesting that equilibration was achieved. The 60-ps simulations require 40 000 timesteps. Trajectory information is saved every 10 steps. The nonbonded interactions between methane and silicalite atoms are represented by a Lennard-Jones and Coulomb potential and are evaluated for all atoms within a cutoff radius of 10.0 A. Interactions between methane molecules are treated by an isotropic, center of mass potential. Electrostatic interactions at the cutoff distance (8-10 A) can be significant, resulting in a discontinuity in the potential. This in turn gives a disastrous lack of energy conservation in MD and also causes problems with convergence in energy minimizations. There are several techniques that can be used to smooth the potential and forces to zero at the cutoff distance. We have used the shifted force potential.28 This method adjusts the potential so that both the energy and force go smoothly to zero at the cutoff distance. With the shifted force potential, the total energy of the systems in the MD simulations showed no drift and had average RMS fluctuations of ~ 0 . 0 5 %over the entire run. C. The Force Field. For the study of nonpolar organics in silicalite, the organic-zeolite interactions are primarily of a van der Waals nature. Therefore, in the calculation of heats of adsorption and/or diffusion of molecules in zeolites, the potentials that describe adsorbate-adsorbate and adsorbate-zeolite dispersion forces are crucial. The electrostatic potentials are expected to play a smaller role in the modeling of nonpolar organics unless aluminum or other metals and charge-balancing cations are present. There are several well-tested molecular mechanics force fields available for C M energy minimization and MD simulations. The most popular, such as AMBER,29 CHARMM,30 and GROMOS,)' are specialized for proteins and organic molecules and cannot be applied to modeling of zeolites without extension or modification. Since the modeling of zeolites is very recent, there is still much to be learned about the proper methods of simulating these systems and the force fields required to represent intra- and intermolecular interactions. While both bonded and nonbonded interactions involving carbon, hydrogen, and oxygen have been parametrized in the protein force fields, little is known about the correct potentials for silicon and other metal atoms found in zeolites. In addition, it is not clear that existing parameters for oxygen are valid when it is in the silyl ether type environment (Si-&Si) found in zeolites. The nonbonded potentials can be derived in several ways. For example, potentials are often obtained by utilizing experimental crystal structures and heats of formation and sublimation. Several worker^^^.^^ have suggested that Lennard-Jones terms for zeolites can be derived via the Slater-Kirkwood equation,32 in a slightly modified form first proposed by Scott and Scheraga.33~~~ This method requires values for the van der Waals radii and polarizabilities of the atoms in question. Unfortunately, these values are considerably in doubt. For example, the reported van der Waals radius of silicon ranges from 0.535 to 2.7 A.36 LennardJones potentials can also be parametrized by ab initio methods, although this is difficult to do correctly. Hartree-Fock calculations by themselves give only a representation of short-range repulsion and electrostatic interactions between molecules.

Nicholas et al. Higher levels of theory, including somelevel of electron correlation, are required to study dispersion forces. Readers interested in the amount of effort needed to obtain nonbonded potentials via quantum mechanical methods are referred to Clementi et Although force fields have been proposed for zeolite work,19,2',38-39 no definitive potentials are currently available that contain all the necessary valence and nonbonded terms for zeolites and adsorbates. Thus, the correct form of the nonbonded interactions is still open to considerable debate. An important part of our initial work was assessing the accuracy of four force fields that contain the nonbonded interactions needed to model the adsorbate-silicalite system. The force field is represented by a sum of potential energy terms

where the first three *valence" terms represent harmonic bond stretches, harmonic bond angle bends, and periodic torsional potentials. The remaining terms represent Lennard-Jones and electrostatic nonbonded interactions. The dielectric constant c is chosen to be 1.O. Because the geometry of themethane molecule is held fixed during the simulation, the valence terms are not used. Thus, the methane interacted with the zeolite only through the Lennard-Jones and electrostatic potentials. The exception to the use of eq 1 was that the methanemethane interactions are not treated on an atom-atom basis but are instead represented by a Lennard-Jones potential based on the centers of mass of the molecules.4o This potential is already well established for the representation of methanemethane interactions, thus assuring that any deviations of our results from experimental values are due only to the methane-silicalite potentials. The partial charges for both the adsorbate and zeolite atoms can be obtained by several methods. We feel that the best partial charges are those that duplicate the electric field around the molecule. The possible point charges for the zeolite lattice have been discussed in detail in our previous study of silica sodalite.' For this work, the zeolite charges are q(Si) = 1.2 and q(0)= -0.6. For methane, partial charges of q(H) = 0.143 and q(C) = -0.572 are used, which duplicate the experimental octapole moment of the molecule. D. Propane. After the accuracy of the potentials is verified with an extensive set of methanesilicalite simulations, a more ambitious approach is taken for the propane simulations. Here the internal coordinates of propane are allowed complete flexibility. The complete force field as presented in eq 1 is used for the simulations of propane. Parameters for the valence and Lennard-Jones part of the potential are from the MM2 force field.41 An energy minimization of the propane molecule in the vacuum gives internal coordinates of r(C-H) = 1.1 A, r(C-C) = 1.5 A, and e(C-C-C) = 11 1.7'. The nonbonded interatomic interactions are evaluated for all atoms within a cutoff radius of 8.0 A. Intermolecular nonbonded and electrostatic interactions are calculated on an atom-atom basis for all the possible zeolite propane and propane-propane interactions. Intramolecular nonbonded interactions in the propane molecules are calculated for all atoms except those that are separated by one bond (1-2 interactions) or by two bonds (1-3 interactions). The 1-2 and 1-3 interactions are assumed to be entirely covalent and accounted for by the valence potentials. For the M D simulations the equations of motion for the system are integrated using the leapfrog algorithm42 in a modified version of the program M o l ~ i m .A~ ~time step of 0.5 fs is needed to achieve energy

Diffusion of Propane and Methane in Silicalite conservation with C-H bond vibrationsincludedin thesimulation. Each system is equilibrated at 300 K for 40 ps or more by loosely coupling the system to a temperature bath44with a temperature relaxation time of 0.1 ps. Following equilibration, particle trajectories are computed for an additional 40 ps, during which the atomic positions and velocities are saved for later analysis. Note that both the equilibration periods and final MD runs for propane are much longer in this work than the 5- and 3-ps runs used in previously published This approach greatly increases the statistical accuracy in time-averaged properties. Determination of acceptable partial charges for propane is considerably more difficult than for methane. Partial atomic charges derived by the distributed multipole analysis (DMA) have been found to give a representation of the electric field around a molecule similar to charges that are fit to the electric fieldcalculated by ab initio methods.45 In keeping with our desire to duplicate the electron field of the adsorbates, the charges for propanearederivedusing the DMA. Thegeometryof the propane molecule was first optimized at the Hartree-Fock level in C2” symmetry using a TZP (triple zeta plus polarization functions)46 basis set. A force calculation at the optimized geometry gave only positive eigenvalues, indicating that the molecule was not at a saddle point and a local minimum in the potential surface had been found. The DMA was then calculated from the TZP wave function. Because the methyl hydrogens are not equivalent, the DMA assigns different charges to them. In the MD, the methyl groups are allowed full rotation, which requires that the methyl hydrogens have equivalent charges. Therefore, each methyl hydrogen is assigned an average charge. The charges for the terminal carbon and hydrogen charges are q(C) = 0.099 and q(H) = -0.035. The central carbon and hydrogen have charges of q(C) = 0.15 and q(H) = -0.056. All the ab initio calculations were done with the program GAMESS-UK.47 E. General Methods. The number of adsorbate molecules used in the simulations is determined by the fixed volume of the zeolite lattice and thedesired concentrations. Previous simulations of methane and propane used only a 1 X 1 X 2 grouping of silicalite unit cells.24 This approach requires that a relatively small number of molecules be used, even for high loadings. By using a larger simulation cell in this work, we increase the number of molecules substantially. For example, we use 216 molecules to create a loading of 8 mol/u.c., while previous workused only 16molecules. Increasing the number of molecules provides a much better statistical estimation of average properties of the system, such as diffusion rates. Simulation of infinite dilution is done by turning off all methane-methane interactions so the molecules interact only with the zeolite. Infinite dilution calculations are done with 108 molecules to achieve adequate statistical averages. For all the simulations, the molecules are initially distributed randomly in the lattice with an equal number in the intersections, straight channels, and zigzag channels. Energy minimizations, using steepest descents and conjugate gradients,4s are then done to relieve any large repulsions between the molecules and the zeolite. To begin the MD equilibration, each molecule is given a random translational and rotational velocity corresponding to the desired temperature. In previous studies of adsorbates in silicalite the interactions of the adsorbed methane with the silicon in the lattice have been ign0red.2~ In the case of methane, diffusion is largely controlled by the size and shape of the channels, which is in turn is determined primarily by the interactions of adsorbate hydrogens and zeolite oxygens. Indeed, experimentally, the diffusion of methane in ZSM-5 is not affected by the Si/Al ratio.s This highlights the weak influence of the electrostatic interactions in methane diffusion. This is not, however, the general case. For polar molecules the roleofelectrostaticswill be important. For example, the diffusionof water becomes slower as the amount of aluminum in the lattice is increased, which is indicativeof significant cation-

The Journal of Physical Chemistry, Vol. 97, No. 16, I993 4153 water interactions.* By including silicon in the simulation at this stage, we show that the nonbonded interactions between adsorbate and silicon are reasonable, paving the way for further simulations that include aluminum or other framework metal atoms and charge-balancing cations. We have also chosen to include electrostatic interactions between all atoms. Including electrostatic forces makes the simulation more realistic and will be required in future simulations that involve charges species. Finally, for the simulations of propane, we included all internal degrees of freedom in the adsorbate molecule. The effects of molecular flexibility are vital in modeling the true diffusive behavior of propane, where the molecule might flex as it negotiates the movement between channels and intersections. Indeed, molecular migrations in previous simulations of propane may have been hindered by using a rigid-molecule as~umption.2~ The simulations of both methane and propane are computationally intensive because of the large number of atoms (over 10 000) and atom-atom interactions (over 200 000) involved. The calculations required the use of supercomputers and were performed on a CRAY-2 at the Magnetic Fusion Energy Computer Center (MFECC) and a CRAY-2 at Cray Research, Inc., in Minnesota. Data analysis was done on an Ardent Titan at the Theoretical Chemistry Group, Argonne National Laboratory, Argonne, IL. 111. Results and Discussion

A. Comparisonof Heats of Adsorption and DiffusionConstants for the Different Force Fields. The Lennard-Jones terms for carbon, hydrogen, and oxygen are well-knonw from the large number of studies that have been done on organic molecules and proteins. Still, the values for these terms are not in complete agreement in different force fields. Considering that we desired to include all atom-atom interactions, Lennard-Jones terms are also required for silicon, for which much less work has been done. Nonbonded potentials were taken from the following sources that contain the necessary Lennard-Jones terms to model all the atomatom interactions of hydrocarbons with the silicalite lattice: 1. Allir~ger.~~ The MM2 force field of Allinger has been shown to accurately describe the structure and thermodynamicproperties of a wide range of organic molecules. The MM2 nonbonded potentials have been derived primarily by calibration to experimental properties. The nonbonded parameters for Si were successfully used by O ~ e l l e t t in e ~molecular ~ mechanics studies of organometallic compounds. We have also found these potentials to give very good agreement with experimental structural and dynamic data in our simulation of silica sodalite.’ 2. G r i g o r a ~ .These ~ ~ nonbonded potentials were derived from ab initio calculations at the Hartree-Fock level with a 3-21G* (modified) basis set. The Lennard-Jones potentials, used in conjunction with charges from extended Huckel theory (EHT) calculations, were shown to model accurately the conformations of chain and ring siloxanes. We tested the Grigoras potentials in this work using charges that are similar to the EHT values. 3. Trouw. Nonbonded potentials were derived by extrapolation of values for other at0ms.2~ The potential parameters are given in Table I. 4. S p a ~ k m a n .This ~ ~ work presents a set of internallyconsistent Lennard-Jones potentials calculated using Gordon-Kim electron gas theory. These potentials are notably more repulsive than potentials derived by the other methods. A comparison of the nonbonded interactions for a given atom shows that the potentials differ greatly in both the depth of the potential well and the position of the minimum. In addition, the differences between the four force fields are not consistent across the different atom types. This makes it very difficult to assess the importance of the energetics (potential well depth) or atomic volume (van der Waals radius) on the properties of the system that we are interested in. To illustrate the differences between

The Journal of Physical Chemistry, Vol. 97, No. 16, 1993

4154

TABLE I: Valence Force Parameters and Leonard-Jones Parameters Bond Stretch Parameters C-H

4.6 4.4

c-c

1.130 1.523

0.32 0.36 0.45

c-c-c

TABLE II: Comparison of Experimental Heats of Adsorption (( v) - R n and Self-Diffusion Constants for Methane at 300 K with Theoretical Values Calculated from Four Mfferent Sets of Nonbonded Potentials potentials

Bond Angle Bend Parameters H-C-H H-C-C

Nicholas et al.

109.40 109.39 109.50

enthalpy of adsorption (kcal/mol)

self-diffusion constant (XIO-s cm2/s)

Spackman Grigoras Allinger Trouw

-3.6 -5.2 -5.8 -7.6

5.9 0.55 7.9 11.9

experiment

-6.7 f 1.3" -6.lb -4%

11.0 & 5Sd 0.11e

Torsional Parameters R T a t 300 K = 0.6 kcal/mol

0.0 0.0

H-C-C-H H-C-C-C

0.0 0.0

0.2370 0.2670

MM2 Lennard-Jones Parameters H-H

c-c H-C H-0 H-Si

c-0 C-Si

Ai/

BiF

74.3949 287.655 1 152.5339 121.7646 489.7994 236.8849 870.5554

25 177.5410 402 083.4063 109 391.2188 65 393.371 1 632 337.0625 255 792.8750 2 064 560.1250

Trouw Lennard-Jones Parameters AJ

B. k

C-Si

101.1940 723.7 142 435.5883 2774.7623

51 144.3989 1 874 992.8350 425 539.1940 10 962 359.7664

meth-meth

3388.5277

10 562 141.0000

H-0 H-Si

c-0

(mdyn.A)/(mol.A2). A. mdyn/(mol.deg2). deg. kcal/mol. (kcal.A6)/mol. (kcal.AI2)/mol.

I\ .-G

-0.05

I

I\

I

\\

1

\ \

I

\ Allinger

-0.1

-,

r

the force fields, the four Lennard-Jones potentials for 0-H interactions are plotted in Figure 4. The considerable differences in the potentials can be readily seen. For this particular interaction, the MM2 parameters are "softest", with the smallest van der Waals radius and the most negative well depth. The potentials of Trouw are very similar with a slightly shallower well. In contrast, the Grigoras and Spackman potentials are

Adsorption isotherm method of Papp, et ale5' Calorimetric method ofPapp,etal.sl Gaschromatographi~methodofChiangetal.~~ Pulsedfield gradient spin-echo N M R method of Caro et al.* e Membrane flow technique of Hayhurst and Paravar.Is

much more repulsive, having a van der Waals radius -0.7 A greater than the MM2 potentials. The potentials in each force field were derived in a consistent fashion using accepted theoretical techniques, but Figure 4 plainly illustrates the wide differences of opinion currently existing about the form of Lennard-Jones interactions for the atoms involved in zeolite studies. To assess the accuracy of the potentials, we calculated isosteric heats of adsorption and diffusion constants for methane in silicalite using the four nonbonded potential sets. The theoretical methods used for determining these properties from the simulation data are described in a later section. The results of the calculations are presented in Table 11. Samples of experimental values are included both to compare to the theoretical methods and to illustrate thevariance in values obtained by different experimental techniques. Papp et aLsl have reported two values for the heat of adsorption of methane on H-ZSM-5, a zeolite essentially the same as silicalite except for the Si/Al ratio. Using an adsorption isotherm method, Papp determines the heat of adsorption to be 6.74 i 0.12 kcal/mol. This compares to a value of 6.1 kcal/mol using a calorimetric method. In these experiments, there is evidence that the methane is weakly chemisorbed on Brmsted acid sites, giving a heat of adsorption that is probably more favorable than would be expected on silicalite, where there are, in principle, no acid sites. Chiang et report a value of 4.8 kcal/mol for methane on silicalite, substantially lower than the Papp measurements. Stachs3 has reported heats of adsorption for several n-alkanes in silicalite, in which the relationship between the number of carbons and the heat of adsorption is highly linear. Although methane is not included in this study, extrapolation of the linear relationship to the case of one carbon gives a value of 4.8 kcal/mol, which is identical with Chiang's. The calculated heats of adsorption using the four nonbonded potential sets vary over a large range of 4 kcal/mol, on the order of the heat of adsorption itself. Considering the appropriate experimental heat of adsorption to compare to the simulation to be Chiang's value of 4.8 kcal/mol, the Grigoras and MM2 potentials give the most accurate duplication of the experimental measurement. The diffusion data are reasonably well represented by the various potentials; three of the potential sets give diffusion constants that are in good agreement with the PFGSENMR data. However, the relationship between experimental and theoretical values does not parallel that found for the heat of adsorption. For example, while Grigoras' potentials give an intermediate value for the heat of adsorption, the diffusion is much slower than that calculated with any of the other potential sets. Obviously, potentials that give accurate theoretical heats of adsorption do not necessarily give accurate diffusion rates, or vice versa. From the data presented in Table 11, it is clear that both the MM2 and Trouw potentials provide a good representation of the experimentalquantities. However, besides providing a reasonable



a

h

C

Figure 5. Contours representing the time-averaged distribution of the centers of mass of the methanes in the straight channels at intinitedilution.

agreement between theory and experiment, the MM2 potentials are advantageous in that they contain the valence terms needed to model flexible organic adsorbates, which the other force fields do not. We therefore choose to use the MM2 potentials for the remaining simulations of methane and propane. The complete force field is given in Table I. We know of no attempts to utilize the MM2 potentials in molecular dynamics adsorbatezeolite simulations other than our own. The results presented heresuggest that these potentials might be widely used for these types of systems. Because of the tremendous amount of work needed to develop and test a force field, verification of the utility of the MM2 force field for simulations of zeolites is highly desirable. We have recently reported additional valence parameters that, used with the MM2 nonbonded potentials, allow the accurate molecular mechanics treatment of the zeolite lattice.’ Thus, the potentials are now available for simulations of flexible organic adsorbates in a flexiblesiliceous zeolite. Of course, theduplication of experimental heats of adsorption and diffusion constants does not guaranteethat other properties calculated with these potentials will be correct. Still, the duplication of two very different experimental measurements by the MM2 parameters gives us confidencethey are a reasonable description of the van der Waals interactions between the adsorbate and the zeolite. B. Detailed Results for Methane. (i) Heat of Adsorption and Potential Energy Surface for Adsorption. The theoretical heat of adsorption is defined as

H = (U)- R T (2) where ( U)is the ensemble average nonbonded interaction energy between adsorbate and silicalite from the MD simulation at the desired temperature, R is the gas constant, and the temperature T equals 300 K. For comparison to experimental heats of adsorption, (U)is taken from the simulation at infinite dilution. The calculated heat of adsorption agrees with the experimental values that were presented in the previous section (Table 11). At higher loading, favorable adsorbate-adsorbate interactions are present, ranging from -0.1 kcal/mol at 2 mol/u.c to -0.8 kcal/ mol at 12 mol/u.c. However, the average interaction energy between the methanes and silicalitealsochanges at higher loading. Thus, when adsorbate-adsorbate interactions are included in ( U), the heat of adsorption remains essentially unchanged, which is consistent with experimental observations. Figure 5 and 6 are contour plots that represent the proportion of time the center of mass of the membranes spend at various positions in the channels taken from the MD simulation at infinite

c Figure 6. Contours representing the time-averaged distribution of the centers of mass of the methanes in the zigzag channels at infinite dilution.

0 Silicon Oxygen Carbon @ Hydrogen Figure7. Methanemolecules in the twoenergy minima in the intersections. dilution. From Figures 5 and 6 it is apparent that the molecules have sites of preferential residence and that these sites are likely energy minima in regard to methane-silicalite interactions. To obtain the interaction energy, individual molecules were placed at each apparent minima and then subjected to energy minimization. Each of the four intersections in the unit cell has two minimum energy sites of methane-silicalite interaction (Figure 7). The two positions are related by a displacement in 6 of 0.64 A. The repulsion energy between methane molecules prevents the simultaneous occupation of the two sites: however, an individual molecule can rapidly interchange between them. The methane molecule is displaced away from the center toward the walls of the intersection and away from the openings of the zigzag channels. The methane-silicalite interaction energy at these sites is -5.92 kcal/mol. In these positions, the closest 0-H distance is =2.9 A and the closest Si-H distance is a3.4 A. The straight channels also provide eight minimum energy sites per unit cell, with two symmetry-related sites in each straight channel segment (Figure 8). The neighboring sites are a2.57 A apart, which also prevents simultaneous occupation. The interaction energy in the straight channels is more favorable than in the intersections at -6.7 1 kcal/mol. The zigzag channels provide four possible sites

4156 The Journal of Physical Chemistry, Vol. 97,No. 16, 1993

Nicholas et al.

....... ......

^_

t

0 Silicon

0 Hydrogen

Carbon

Oxygen

t

Figure 8. Methane molecules in the two energy minima in the straight channels.

I-. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .t . "

1x10'

0

2

4

6

8

10

12

14

16

Molecules per Unit Cell

CI)

4

:.e

Figure 10. Diffusion constants for methane and propane at various loadings compared to experimental NMR measurements.

"1

to the total MSD. The data shown are for the loading of 4 mol/ The overall appearance of the MSD and relative rate of movement in a, b, and c are very similar for all the loadings. Figure 9 illustrates the anisotropic nature of methane diffusion in silicalite. As expected from thegeometry of thezeolitechannels, diffusion is fastest in thestraight channels (b-direction). Although the sinusoidal channels (a-direction) are only slightly smaller than the straight channels, diffusion in the sinusoidal channels is considerably slower. The calculated rate of diffusion in the b-direction is on averagetwice as fast as diffusion in the adirection. NMR measurements indicate that the rates of diffusion in a and b differ by less than a factor of 5,s4 which agrees with the our findings. Also, as predicted, diffusion in c is very slow. There is an initial sharp rise in the MSD values as the molecules sample the available free volume. This nondiffusive motion continues until collisions with the channel walls or other molecules take place. After this initial movement, the MSDs should reflect the actual translational diffusion of the adsorbate molecules, and the Einstein relation

U.C.

204

0

10

30

20

40

50

60

Time in ps

Figure 9. Raw mean squared displacements of the methane molecules at 4 mol/u.c. showing the anisotropic nature of diffusion.

with energy -6.74 kcal/mol, essentially the same as the interaction energy in the straight channel. Thus, on a strictly energetic basis, there is essentially no difference in preference for the methanes to reside in the straight or zigzag channels. Although the interactions at intersections are defintely less favorable, the k0.8 kcal/mol energy difference and the small energy barriers between sites substantiate the assumption that silicaliteis an energetically homogeneous surface. Although the time distributions change somewhat at higher loading, the positions of the energy minima remain the same until significant crowding of the molecules takes place. (ii) Anisotropy of Diffusion and the Effects of Loading on Diffusion. Figure 9 illustratesthe time evolution of the ensemble average mean square displacements (MSD) of the centers of mass of the molecules in the lattice as defined by I

where rj represents the coordinates of the molecule i and the ensemble average is over all the molecules in the simulation and multiple time origins. The raw MSD are shown for the components of movement in the a, b, and c directions, in addition

(4)

can be used to obtain the self-diffusion constant from the MSD. To investigate the effects of concentration on the rate of diffusion, simulations were performed at six different loadings of methane in silicalite. In simplest terms, increased sorbate concentration is expected to lead to a reduced volume for the individual molecules to occupy and, consequently, reduced translational mobility. A plot of thecalculated diffusion constants as a function of loading is given in Figure 10 with experimental data from the PFGSENMR method. For both methane and propane, the results are in excellent agreement with experiment. We must note that the experimental measurements have an estimated error of 50% for methane and 100%for propane.* In all cases where comparison with experiment is possible, the theoretical diffusion rates are slightly slower. (iii) Energy of Activation. The force field has been shown to duplicate the NMR diffusion measurements at 300 K. Another test of the potentials is their ability to duplicatediffusive behavior at different temperatures. Although nonbonded potentials are generally assumed to be independent of temperature, there is evidence that this is not necessarily true. For example, Momany et al.ss found large changes in Lennard-Jones parameters derived

Diffusion of Propane and Methane in Silicalite from crystal structures of benzene a t various temperatures due to the effects of thermal vibration. The MM2 potentials we use were derived from room temperature data and wil therefore also contain thermal vibrational effects. The MD then adds an additional temperature contribution. Thus, the correct temperature dependence of the potentials cannot be assumed. N M R and quasi-elastic neutron scattering (QENS) diffusion measurements of methane in Na-ZSM-5 at 200 K have been reported by Jobic et al.56 For a loading of 4 mol/u.c., they find a diffusion constant of (2.5-3.7) X c m 2 / ~ .The ~ ~calculated value from our simulations a t 200 K is 3.38 X cm2/s, which is in good agreement with the experimental value. It should be noted that the determination of the diffusion constant from the QENS measurements in these anisotropic systems cannot be made in a straightforward manner. We do not assign the same merit to these results that they would warrant in the case of an isotropic zeolite -57 The diffusion constants at different temperature can be used to determine the activation energy for diffusion58 from

The Journal of Physical Chemistry, Vol. 97, No. 16, 1993 4157 I

I,

16 mol/u.c.

i 4

From the calculated diffusion constants, the energy of activation is 1.01 kcal/mol. This compares to experimental values ranging from 0.9 to 1.2 k c a l / m 0 1 . ~ ~ ~ ~ (iv) The Effects of Volume, Energetics, Anisotropic Diffusion Rates, and Intermolecular Interactions on Packing. The manner in which the adsorbate molecules pack in the framework has been previously proposed by several i These models were based primarily on the size of the adsorbates, adsorbateadsorbate interactions, and the relative energy preferences for the possible residence sites. Although these are clearly important factors that influence adsorbate packing, the simulations show the importance of dynamic effects that have not been previously considered. To study the effects of loading on packing, the time distribution of the centers of mass of the methanes was calculated for thevarious possible locationsof the molecules within thelattice. From the study of the microscopic behavior of the molecules in the zeolite, it is apparent that much time is spent in transit between the various sites. This leads to some ambiguity in the interpretation of the time distribution. For example, a molecule diffusing along the straight channel is classified as being in an intersection as it passes through, even if it continues in the straight channel. We have chosen to use the simple classification scheme, rather than trying to further classify molecules as "in transit". Despite this simplification, this scheme does reveal useful information about the distribution of molecules and the effects of packing. The nature of the packing is thus due to a complex mixture of effects that we describe in the following sections. (a) Volume. In discussion of the packing of adsorbates in silicalite, the total available frameworkvolume can be conveniently divided into three areas, the intersections, the straight channels, and the zigzag channels. These areas are generally considered to be the residence sites for small organic molecules up to the size of propane. For longer chain alkanes, the molecule can no longer be contained entirely within a channel segment. In each unit cell there are four intersections, four straight channel segments, and four zigzag channel segments, for a total of 12 residence sites. For purposes of calculation we have based the division of the volume on a length of 4.5 A in b for the straight channel and 5.4 A in a and b for the intersections. If the placement of the molecules were strictly random and based only on the available volume, we would expect the molecules to be found 33% of the time in the intersections, 27%in the straight channels, and 40% in the zigzag channels. Thus, the zigzag channels have a considerably greater volume than the straight channels or the intersections. The molecules are diffusing relatively rapidly over an energetically homogeneous surface. Thus, in simplest terms, the molecules

I

6

8

10

12

14

16

Distance in A

Figure 11. Radial distribution function (RDF) of the centers of mass of the methanes at various loadings.

might bedistributed according to the availablevolume. Deviations from this volume partitioning highlight the manner in which other factors influence packing. (b) Energetics. From the study of the potential energy surface for adsorption (discussed above), the methanes have a more favorable energy of interaction in the zigzag or straight channels than in the intersections. If we assume that there are an equal number of occupation sites in the intersections, zigzag, and straight channels, we can make a simple approximation of the Boltzmann probability for residence in the various sites by

where P is the probability and the partition function qo is a sum over the various sites of interaction energy Ei. This analysis predicts the molecules will spend 23% of the time in the intersections and 38.5% of the time in the straight and zigzag channels. In contrast to the greater percentage of molecules that would be expected to reside in the zigzag channels due to channel volume, the energetics predict equal populations in both the straight and zigzag channels. (c) Anisotropic Diffusion Rates. We have previously shown that diffusion is slower in the a-direction than in they direction. We have also just shown that this slower rate of diffusion in a is not due to a more favorable interaction energy in the zigzag channels. Motion in the a-direction is not actually parallel to the a-axis due to the twisted geometry of the zigzag channel. Thus, for a given value of the diffusion constant, a larger absolute displacement is needed in a than in b. This does not fully account for the slower diffusion in a. Therefore, the difference in diffusion rates must involve the slightly smaller radius and the difficulty the molecules have in navigating the contorted geometry of the zigzag channels. The slower rate of diffusion in a also suggests that molecules might be expected to spend more time in the zigzag channels than in the straight channels. (d) Adsorbate-Adsorbate Interactions. A more subtle effect in packing is the existence of favorable adsorbate-adsorbate interactions. Figure 1 1 shows the center-of-mass radial distri-


Nicholas et al.

TABLE III: Percentage of Time Molecules Spent in the Intersections and Channels as a Function of Loading

1.6

,A

loading

0'4: 0.2

01 2

1

1

1

1 , 4

6

8

10

12

14

16

I 18

Distance in 8,

Figure 12. Radial distribution function (RDF) of the centers of mass of the methanes at various loadings with attractive potential removed.

bution functions (RDFs) for methane in silicalite over the range of loadings simulated. The integral of the RDFs gives the probability of finding another molecule at a distance r from any given reference molecule. Clearly the RDFs differ greatly from what would be observed for methane gas and show the effects of the geometric restrictions of the zeolite channels. All the RDFs haveapeakat -12.0A thatvariesonlyslightlywithconcentration. This distance corresponds to the distance between favorable sites in the zigzag channels. The time distributions have shown that the zigzag channels are preferentially occupied. The peak at 12 A suggests that if there is a molecule in a zigzag channel, it is likely that the adjacent zigzag channel is also occupied. In contrast to this effect, there is a peak in the RDFs at -4.0 A that is very concentration dependent. This intermolecular distance can be achieved by having one molecule in a favorable site in an intersection and another in a straight channel. Thus, its presence in the RDFs could also be due to the geometry of the lattice. However, the peak is present even a t a very dilute concentration of 2 mol/u.c., in which the molecules would not be expected to spend appreciable time in contact. The m e t h a n e methane Lennard-Jones potential has a minimum of 4 . 0 5 kcal/ mol a t 3.95 A. Therefore, the possibility exists that themolecules maintain close distance because of the favorable Lennard-Jones interaction. To test this case, simulations were done a t 4 and 8 mol/u.c. with the attractive part of the methanemethane Lennard-Jones potential removed. The RDFs for these simulations are shown in Figure 12. There is a pronounced reduction in the peak at -4.0 A for both concentrations. This clearly demonstrates that the methane-methane interactions have an effect on packing even at low concentrations. Noticeably absent from theRDFisapeakat=lO.OAthatwouldappearifmolecules tended to occupy sites in adjacent intersections or straight channels. (e) Discussion. In Table I11 the percentage of time the molecules spend in the various locations as a function of loading is presented. Also given is the average number of molecules per unit cell that the percentages represent. Although the number of molecules per unit cell has no physical meaning at infinite dilution, these numbers are given for comparison. There is a clear tendency for the molecules to preferentially reside in the channelsand not in the intersections, compared to thedistribution predicted by volume or energetics. To illustrate the nature of the

methane volume Boltzmann population infinite dilution 2 mol/u.c. 4 mol/u.c. 8 mol/u.c. 12 mol/u.c. 16 mol/u.c. propane Boltzmann population infinite dilution 4 mol/u.c. 12 mol/u.c.

locations of the molecules straight zigzag intersections channels channels 33.0 23.0 12.7 13.5 (0.25)a 14.9 (0.60) 18.2 (1.46) 20.8 (2.50) 23.6 (3.77)

27.0 38.5 34.1 30.2 (0.60) 29.0 (1.16) 32.0 (2.56) 29.9 (3.59) 28.2 (4.52)

40.0 38.5 53.2 56.3 (1.13) 56.1 (2.24) 49.8 (3.98) 49.4 (5.93) 48.2 (7.71)

0.9 25.0 34.9 (1.40) 26.8 (3.22)

33.6 17.4 24.7 (0.99) 33.5 (4.02)

65.5 57.6 40.4 (1.62) 39.7 (4.74)

"Average number of molecules per unit cell represented by the percentage.

packing, consider the distribution of molecules at a concentration of 12 mol/u.c. It could be naively assumed that for this loading, each of the occupation sites could be occupied, giving a total of four molecules in each area. However, the molecules clearly favor the zigzag channels, such that there are on average six molecules per unit cell in the zigzag channels. An analysis of the center-of-mass time distribution a t this loading reveals that in addition to the primary maximum found in the zigzag channels a t lower loading, two new sites in the zigzag channels are now occupied, on opposite sides of the site favored a t infinite dilution, -4.0 A apart. The lack of molecules in the intersections is reasonable considering either the volume differences or that the interaction energy is less favorable by -1 .Okcal/mol. However, the favoring of the zigzag channel over the straight channel cannot be explained on strictly energetic grounds considering that the interaction energy isessentially thesamein thestraight andzigzag channels. The difference apparently lies in the geometry and smaller size of the zigzag channel and the resulting slower diffusion through this channel. As the loading is increased, the distribution of molecules lends toward that which would be predicted on a volume basis, but still maintains a preference for the zigzag channels over the intersections. Although theories that predict the preferred distributions of molecules in the framework have taken into account energetic effects, the simulation demonstrates that dynamic factors previously overlooked are also important. (v) Jump Diffusion Model. The manner in which molecules diffuse in zeolites has often been described in terms of a jump diffusion model (JDM). The JDM assumes that molecules occupy particular sites in the zeolite lattice for a mean residence time of T and periodically make jumps between sites of an average length ( 12) li2. The transit time between the localization sites is assumed to be much less than the residence time. The diffusion constant is related to these values by

D = ( l2)'I2/6r For example, Caroet a1.8have proposed that diffusionof methane and propane in silicalite occurs as activated jumps between channel intersections. Using the JDM model Car0 et al. obtain jump lengths that are similar to the distance between intersections. For propane they find an average jump length of -11.7 A a t concentrations less than 12 mol/u.c. that abruptly decreases to -7.7 A a t higher loading. The distance between intersections is -10 A in the straight channels and -12 A in the zigzag channels. Note that the distance between the most favored occupation sites in the zigzag channels is also -12.0 A. To illustrate the nature of diffusion at low loading in the simulation, and to contrast it to the type of behavior predicted by the JDM, Figure 13 presents the actual 60-ps MD trajectory



t

I

I

I

Trans + Rot 50

\

40

20

10

Total - (Trans + Rot) 0

10

20

30

40

50

80

Time in ps

Figure 13. Mean squaredisplacementsof an individualmethanemolecules indicating the detailed nature of molecular motion.

0

50

100

150

200

250

cm"

Figure 14. Spectral density for translational, rotational, translationalrotational coupling, and total velocity.

of two representative molecules from the 4 mol/u.c. simulation. The coordinates are given in A. From these, the molecules can be classified according to their location in intersections or in straight or zigzag channels. In Figure 13, molecule A begins in an intersection, moves to a zigzag channel at -25 ps, and rattles back and forth between the intersection and zigzag channel between 40 and 50 ps before finally moving to a straight channel. In contrast, molecule B moves back and forth between various favored sites before making a long excursion in the straight channel between 25 and 40 ps. This extensive excursion of a single molecule is reflected in an increase in the ensemble average mean square displacements (Figure 9)during this time period. We do see that methane diffusion is somewhat jumplike, with principal jumps occurring along the straight channel (b-direction). Residence time is longer in the zigzag channels. For a mean jump length of 11 A between intersectionsand our theoretical diffusion constant of 7.9 X 10-5 cm2/s at a concentration of 4 mol/u.c., the jump diffusion model predicts a mean residence time between jumps of 25.5 ps. Thus, the JDM suggests that we should witness many diffusivejumps during the 60-ps simulations, as we illustrate. However, our data also show that the transient diffusive motions of individual molecules can bevery different and that their motion is much more complexthan the JDM assumes. At higher loadings the crowding of the channels prevents long excursions between sites and diffusion follows the JDM more closely. We must mention that the smooth methane-silicalite potential energy surfacewould tend to mitigate against the type of motion predicted by the JDM. As the potential energy barriers between sites is increased, the JDM becomes more viable. The confinement of the molecules in the zeolite channels not only causes diffusion to be highly anisotropic but also leads to a very inhomogeneous type of diffusive behavior in the individual molecules. An important practical consequence of this type of diffusion pertains to the number of molecules and the amount of time used in the MD simulation. We have shown that the movement of one molecule out of 108 can substantially affect the ensemble average diffusion rate over a period of approximately 15 ps. This suggests that the use of small numbers of molecules and short MD runs gives a large possiblity for error in the calculated properties of the system. Although it would be desirable to use smaller simulations systems, calculations of the size presented here seem close to the lower limit for obtaining proper statistical averages.

(vi) Velocity Autocorrelation Functions. One major strength of theoretical simulations is their ability to explore properties of the system that are difficult or impossible to obtain experimentally. For example, the velocity autocorrelation function (VACF) can be used to investigate the possibility of coupling between the rotational and translational motion of the molecules in the lattice. The VACF is obtained by calculating the dot product of the velocity of a molecule at time zero with later times t . The VACF thus contains information about any periodic fluctuations in the velocity of the molecules. To minimize the effects of spurious fluctuations and noise, the VACF is calculated over the ensemble of molecules and many time origins. The frequency spectrum for the motion of the molecules, which correspondsapproximately to an IR spectrum, can be obtained by the Fourier transform of the total velocity autocorrelation function

Following Impey et al.,59the total velocity of the rigid molecules can be decomposed into the translational velocity of the center of mass of the molecule and a rotation about the center of mass (8) Vtotal = Vtrans + 'rot By expanding the total velocity into the translational and rotational parts (eq 8), the frequency spectrum for total velocity can be shown to also equal the sum of the Fourier transforms of the translational, rotational, and translation-rotational coupling autocorrelation functions

F(U) = So(vtra,,(o)'vtnns(t)) + (vrot(o)*vrot(t)) + 2(v,rans(O)*vrot(t)) e x ~ ( i ~dtt ) ( 9 ) The frequency spectrum for total, translational, rotational, and coupling velocity can each be obtained separately from the simulation by the Fourier transform of the appropriate autocorrelation function. The frequency spectrum due to coupling can also be obtained by subtraction of the translational and rotationalparts from the total. Figure 14 illustrates the spectrum for the various types of motion obtained for a concentration of 4 mol/u.c. The translational velocity spectrum exhibits weak peaks at 20 and 60 cm-I while the rotational velocity spectrum has a more intense peak at 100 cm-I and some suggestion of a

4160


shoulder at 80 cm-1. The sum of the translational and rotational parts closely matches the spectrum calculated for the total velocity. The difference between the total and the sum of the rotational and translational parts is negligible, indicating little or no translational-rotational coupling. (C) Detailed Results for Propane. The results obtained in the simulations of rigid methane suggest that the MM2 nonbonded potentials can reasonably duplicate available experimental measurements and therefore give an accurate and predictive description of the real system. These preliminary simulations set the groundwork for the more ambitious simulations of propane. In these calculations we were dealing with a larger number of atoms (over 11 000) and a much large number of atom-atoms interactions (over 1 600 000) than in the methane/silicalite systems, requiring considerably more computational effort. In addition, the propane molecules were allowed full internal flexibility, which also requires some additional computation. These simulations offered a test of the MM2 valence potentials and represent the first attempt to model a totally flexible adsorbate in a zeolite. (i) Heat of Adsorption and Potential Energy Minima. At infinite dilution, the ensemble average interaction energy between the propane molecules and the zeolite gives a theoretical heat of adsorption of-10.3 kcal/mol. Thisvalue is in excellent agreement with an extrapolated value of -9.7 kcal/mol from the data of Stach53 and -9.6 kcal/mol found by Richards and Rees." As a caution to other researchers we must mention that at the very low diffusion rates of propane, the equilibration of the distribution of the molecules in the infinite dilution simulations is very slow. The monitoring of the conservation of energy would suggest that equilibration is achieved in 40 ps or less. However, the time distributions of the molecules in the various locations show that 180ps is needed before the configurations stabilize. The relatively faster diffusive motion of methane and the presence of moleculemolecule interactions in the simulations at higher loading allow equilibration to take place in much less time. To provide some measure of the accuracy of the calculated properties, five 40-ps simulations of propane at infinite dilution were done. Unfortunately, the computational demands of the simulations prevented us from obtaining multiple runs of the other loadings. From the five infinite dilution results the root mean square error in the heat of adsorption is f0.08 kcal/mol. The error in the diffusion constant is f0.15 X 1Ck5cm2/s. Because the other propane simulations used an equal or greater number of molecules, we expect similar accuracy in corresponding diffusion constants. Figures 15 and 16 show the time distributions of the centers of mass of the propane molecules at infinite dilution in the straight and zigzag channels, respectively. The positions of maximum density are grossly similar to those of methane; however, notable differences are present. In Figure 15, the most favored location for propane is a broad area in the center of a straight channel, unlike the two symmetry-related maxima of methane. In contrast, the favored location in the zigzag channel is very close to that of methane; however, the propane site is in the middle of the channel rather than toward the wall. Notice that the volume available for occupation by the molecules is much less than that of methane, as expected due to the larger bulk of propane. Increasing the loading of 4 mol/u.c. does not appreciably change the time distributions; however, a t 12 mol/u.c., marked changes take place. In the straight channels, the molecules are largely isolated in the middleof thechannel. In contrast to the continuous distributions found for methane, there are areas where the molecules are not found during the time scale of the simulation. The trajectories were sampled at 100-step intervals during the simulation, so some information is possibly lost about transitions between sites. However, isolation of the molecules in preferred sites is primarily due to the very slow rate of diffusion.

Nicholas et al.

h

a Figure 15. Time distribution of the center of mass of propane in the straight channels at infinite dilution.

a

c Figure 16. Time distribution of the center of mass of propane in the

zigzag channels at infinite dilution. The positions of the energy minima were obtained from the time distributions in the same way as was done for methane. However, the potential energy surface for motion around the minima appears to be very flat. In addition, the conformational flexibility of propane makes the determination of the true minima considerably more difficult since many conformationsof propane have nearly the same interaction energy with the silicalite framework. Without a deep potential well, the molecules can assume many conformations that are very close in energy. We have made a concerted effort to locate the global minima, but we are not certain they have been found. Similar to the case of methane, each of the four intersections in the unit cell offers a minimum energy site of propane-silicalite interaction. The propane molecule is slightly displaced away from the center of the intersection toward the walls of the intersection and away from the openings of the zigzag channels. The propane-silicalite interaction energy a t these sites is -9.7 kcal/mol. The straight channels also provide four minimum energy sites with themoleculesituated in the middleof thechannel. The interaction energy in the straight channels is more favorable than in the intersections at -1 1.9 kcal/mol. The zigzag channels


Diffusion of Propane and Methane in Silicalite 6Oi

B

6 B

'3

l

t

3 ;

I

2-

k E

'=

b !i

1.5

a

0'5] 0

0 0

5

10

15

20

25

30

35

40

Time in ps

Figure 17. Mean squaredisplacements of propane at 4 mol/u.c. illustrating the anisotropic nature of diffusion.

provide four possible sites with energy -12.3 kcal/mol, which is slightly more favorable than the interaction energy in the straight channel. Compared tomethane, there is some energetic preference for the zigzag channel compared to the straight channels, indicating that the larger bulk of the propane gives a better fit to the smaller diameter of the zigzag channels, thereby enhancing the van der Waals interactions. The interaction energy at the intersections is less favorable by =2.0 kcal/mol, indicating that silicalite offers a less energetically homogeneous surface to propane than methane. (ii) Anisotropy of Diffusion and the Effects of Loading on Diffusion. The mean square displacements for the propane molecules (Figure 17) demonstrate an isotropic behavior very similar to the methane molecules. Diffusion is dominated by motion along the straight channels (b)with diffusion in the zigzag channels (a) being about a factor of 2 slower. The calculated diffusion constants for propane are shown in Figure 10. As expected, the overall diffusion rates are much slower than for methane, due to the larger size of the propane. The diffusion constants are also in excellent agreement with the PFGSENMR measurements. Interestingly, the rate of diffusion at infinite dilution is slightly slower than a t a loading of 1 mol/ U.C. The attraction of the propane molecules for the channels walls is considerably stronger than methane. The kinetic energy of the system must overcome this attraction for diffusion to take place. Apparently the propane-propane repulsions and collisions that take place at higher loading actually facilitate diffusion. Thus, it appears that propane-propane repulsions and/or collisions actually increase diffusion by overcoming the attraction of the propane for the silicalite channels. This leads to a faster diffusion rate than that at infinite dilution. Experimental data are not available to verify this behavior; however, June et al.19 found a similar behavior in their simulation of methane at 200 K, and this was also observed experimentally.60 (iii) The Effects of Volume, Energetics, Adsorbate-Adsorbate Interactions, and Anisotropic Diffusion Rates on Packing. The proposed models of packing were discussed in the previous section on methane. Some of the factors that determine packing are similar to those that effect methane. For example, the division of the channel volume does not change, and the relative rates of

2

;/

,

,

,

,

,

,

4

6

8

10

12

14

16

18

Distance in A

Figure 18. Radial distribution function for propane in silicalite at 4 and 12 mol/u.c.

diffusion in a and b have been found to be the same. The factors that do change in the propane case are described in the following sections. (a) Energetics. Compared to methane, the relative interaction energies between the propane and the silicalite surface vary over a wider range. From eq 6 the energetics predict the molecules to be distributed as 0.9% in the intersections, 33.6% in the straight channels, and 65.5% in the zigzag channels. As we had stressed before, this partitioning is quite approximate. Still, it shows the occupation sites in the intersections are particularly disfavored for propane. (b) Adsorbate-Adsorbate Interactions. The radial distribution functions for the centers of mass of the two loadings of propane are shown in Figure 18. Comparison of the RDF shows that peaks are present at the same intermolecular distances for the two loadings. The first peak is again due to propanepropane interactions, while the origin of the other peaks can be attributed to the different distances between favored sites in the lattice. The peak a t -8.0 A is the distance between the preferred sites in the straight and zigzag channels. The distance between sites in the straight channels shouldgivea peak at -10.0 A, while the distance in the zigzag channels is -12.0 A. Also present is a peak at .=13.5 A, which represents the translational symmetry of the silicalite in the c direction. In a similar fashion, the peak at -16.0 A could be related to the peak at -8.0 A for molecular neighboring pairs that span two unit cells. The most notable difference between the two loadings is the peakcorresponding to propanepropane interactions. At thelower loading of 4 mol/u.c., this distance is =6.0 A. As the molecules are crowded together at 12 mol/u.c., the intermolecular distance is reduced to 4 . 0 A. This behavior was not seen for methane, in which the first peak at -4.0 A did not shift as loading was increased. In contrast to the isotropic methane molecules, the propanes have orientational and internal flexibility. This could allow for different conformations of the molecule that may tend to favor different intermolecular distances. At the lower loading of propane, the intermolecular repulsions are not apparently high enough to induce orientational or conformational changes in the propane molecules. This might be termed a loose packing arrangement. At higher loading, the molecules are induced to change orientation or conformation to achieve an optimal packing. Similar to the methanes, the propane-propane peak height is much more pronounced at the higher loading. It is obvious the

4162


Nicholas et al. 19. Also shown is the distribution calculated for the free molecule. These data are taken from a simulation of 108propanes in vacuum (no propane-silicalite or propanepropane interactions were included). From the distribution it can be seen that in both cases the C-C-C angle varies over approximately 100-12S0. The distribution is shifted very slightly to larger values when the propane is in the lattice and the average angle increases from 113.1 to 113.6O. Thus, the zeolite has a very minor effect on the internal motion of propane, but the internal flexibility of the molecule is much as it would be in a pure gas. IV. Conclusion

0 100

105

110

115

120

125

C-C-C Angle in Degrees

Figure 19. Distribution of C-C-C bond angles for the free molecule and the molecule in silicalite.

propane is moving so slowly at the 12 mol/u.c. loading that it appears to be similar to a solid, with very sharply defined peaks in the RDF. (c) Discussion. The relative percentages of time the propanes are found in the various channel locations are given in Table 111. The preference of the propanes to reside in the various channel segments is markedly different from that of methane at lower loadings. A substantially larger number of molecules are found in the intersections than the straight channel segments. This is despite the calculated interaction energy being less favorable in the intersections. Contrary to the distribution expected on energetic grounds, the propane molecules might spend more time in the intersections because there are orientational restrictions that make it difficult to re-enter a straight or zigzag channel segment once they are in an intersection. Thus, the simulations show that preferred distribution of molecules in the zeolite cannot be predicted on a strictly energetic basis. At higher loading, the propane is moving slowly enough that the dynamic effects we found to be important in methane are diminished. At 12 mol/u.c. the molecules are close to being evenly divided between the three residence sites, which is in agreement with naive assumptions of packing. Note that the volume of the propane molecules prevents more than one molecule from occupying a residence site, so that the Boltzmann population distributions cannot be achieved at this loading. (iv) Internal Flexibility of Propane. In our study of methane diffusion, the molecule was held rigid because we assumed the internal flexibility of methane would make a negligible contribution to diffusive behavior. This assumption would certainly not hold for long-chain hydrocarbons and might not be valid for propane, because the actual molecule has considerable flexibility in the C-C-C bond angle. Because the assumption of a rigid propane seemed untenable, the propane simulation allowed complete internal flexibility of the adsorbate. Although the H-C-H bonds were allowed to flex and the end-methyl groups could rotate around the C-C torsions, these motions do not substantially effect the shape or size of the molecule and are also expected to have only slight effects on the energy of propanesilicalite interactions. Thus, we will consider only the movement of the C-C-C bond angle. To illustrate the internal flexibility of propane, the distribution of C-C-C angle for the 4 mol/u.c. simulation is shown in Figure

The important contribution of this work is the development of a reasonable force field for the treatment of adsorbate diffusion and dynamics in molecular sieve silicates that explicitly takes into account the interactions with all the framework atoms and which incorporatesall internaldegreesof freedomoftheadsorbate. The results represent the largest and, hence, the most statistically reliable simulations yet made in this area and reveal the importance of the diffusive dynamics of molecules in determining the time averaged spatialdistribution of theadsorbatein the pores. Simple concepts based on the assumed energetics of adsorbate interactions give a misleading picture for both methane and propane in silicalite, and simplified models of the diffusive motions for weakly interacting adsorbates are also misleading. The results also reinforce the necessity for using large ensembles and long run times in these simulations and illustrate the difficulty faced in tackling simulations wherein the adsorbate diffusivity isvery slow. Specifically, the results presented here clearly demonstrate that the diffusive and energetic behavior of small hydrocarbons in silicalite can be modeled by energy minimization and MD methods. The results also indicate that the MM2 force field provides a reasonable representation of the inter- and intramolecular forces that are needed for zeolite modeling. Another important result obtained by this work is the demonstration that various nonbonded parameter sets all derived in a reasonable and consistent manner give widely differing values for the heat of adsorption and diffusion constant of methane in silicalite. Clearly, this illustrates the danger in accepting theoretical properties calculated with the force field concept without prior testing or calibration of the potentials against experimental data. This is especially true when working with materials, such as zeolites, that have not been subject to a large amount of theoretical study. We feel there are still improvements that need to be made to the MM2 potentials to better suit them for adsorbate-zeolite simulations. However, the considerable effort already invested in their derivation and calibration substantially lessens the amount of additional work needed, particularly compared to developing completely new potentials. The calculated diffusion constants agree very well with the experimental values obtained with N M R methods. This suggests that these microscopic experimental methods provide the most accurate data in adsorbate diffusion. The large deviations in diffusion rates obtained by macroscopic techniques could be attributed to heat and mass transfer effects, the methods of preparation of the samples, and real differences in diffusion in different zeolite samples. However, the simulation results suggest that the macroscopic methods that have purported to minimize these effects have done so. An important aspect of the simulations to consider is the assumption that the framework can be held rigid. Although silicalite at room temperature can be expected to exhibit only modest displacement amplitudes for the atoms in the framework, it is reported to undergo a slight change in symmetry upon the adsorption of methane. The flexibility of the lattice could influence thediffusion of adsorbates in several ways. For example, although the displacements of individual atoms in the lattice would be small, slow breathing motions of the lattice occur that might

Diffusion of Propane and Methane in Silicalite facilitate diffusion of adsorbates. A more fundamental aspect of this concerns energy transfer between adsorbate and framework. This energy transfer is not permitted in the rigid framework approximation, and energy exchange in these simulations occurs only via intermolecular collisions. The significance of this is greater in situations where the temperature is low or the molecular diffusivity is intrinsically slow so as to make encounters infrequent in the duration of the simulation. It does not appear to cause significant inaccuracy in the results presented.

V. Areas of Future Work The verification of the MM2 potentials sets the groundwork for simulationof many interestingorganic adsorbates in the zeolite lattice. Combined with our previous work in deriving accurate force field parameters for simulation of the zeolite lattice, it is possible to perform energy minimization and MD simulations of flexible adsorbates in a flexible lattice. More importantly, they will enable the reliable calculation of systems at elevated temperatures and low loading where there is a dearth of experimental data. Simulations of this type can be used to investigate phase changes in zeolites and changes in zeolite structure that accompany adsorption.

Acknowledgmeat. John Nicholas would like to personally thank Dr. Thom Dunning, Jr., Material Science Research Center, Battelle, Pacific Northwest Laboratories, and Dr. Albert F. Wagner, Dr. Robert J. Harrison, Dr. Lawrence B. Harding, and Dr. Ron L. Shepard of the Theoretical Chemistry Group at Argonne National Laboratory for their generous discussion, assistance, and encouragement. This work is supported in part by the U.S.Department of Energy, BES-Material Sciences, under Contract No. W-31-109-ENG-38. References and Notes (1) Young, L. B.; Butter, S. A.; Kaeding, W. W. J . Catal. 1982,76,418. (2) Meisel, S. L.; McCullough, J. P.; Lechthaler, C. H.; Weisz, P. B. CHEMTECH. 1976,6, 86. (3) Haag, W. 0.;Lago, R. M.; Weisz, P. B. Faraday Discuss. 1982,72, 317-330. (4) White, J. C.; Hess, A. C. To be submitted for publication in J .

Phys. Chem. (5) Curtiss, L. A.; Brand, H.; Nicholas, J. B.; Iton, L. E. Chem. Phys. Lett. 1991, 184, 215-220. (6) Vetrivel, R.; Catlow, C. R.; Colbourn, E. A. J . Chem. SOC.,Faraday Trans. 1989, 85, 497-503. (7) Nicholas, J. B.; Hopfinger, A. J.; Iton, L. E.; Trouw, F. R. J. Am. Chem. Soc. 1991, 113, 4792. (8) Caro, J.; Bulow, M.; Schirmer, W.; Karger, J.; Heink, W.; Pfeifer, H.; Zdanov, S. P. J. Chem. SOC.,Faraday Trans. I1985, 81, 2541-2550. (9) Jacobs, P. A.; Beyer, H. K.; Valyon, J. Zeolites 1981, I , 161. (IO) Thamm, H. Zeolites 1987, 7, 341-346. (1 I ) Richards, R. E.; Rees, L. V. C. Langmuir 1987, 3, 335-340. (12) Karger, J.; Ruthven, D. M. Zeolites 1989, 9, 267. (13) Eic, M.; Ruthven, D. M. In Zeolites: Facts, Figures, Future; Elsevier: Amsterdam, 1989. (14) Bulow, M.; Schlodder, H.; Rees, L. V. C.; Richards, R. E. in Proceedings of the 7th InternationalZeolite Conference;Elsevier: Amsterdam,

The Journal of Physical Chemistry, Vol. 97, No. 16, 1993 4163 (18) Demontis, D.; Suffritti, G. B.; Alberti, A,; Quartieri, S.; Fois, S.; Gamba, A. 1986, 116,459-466. (19) June, R. L.; Bell, A. T.; Theodorou, D. N. J . Phys. Chem. 1990,94, 8232.

(20) Cohen De Lara, E.; Kahn, R.; Goulay, A. M. J. Chem. Phys. 1989, 90, 1485. (21) Pickett, S. D.; Nowak, A. K.; Thomas, J. M.; Peterson, B. K.; Swift, J. F. P.; Cheetham, A. K.; Ouden, C. J. J. d.; Smit, B.; Post, M. F. M. J. Phys. Chem. 1990, 94, 1233. (22) Leherte,L.; Vercauteren, D. P.;Derouane, E. G.; Lie,G. C.;Clementi, E.; Andre, J.-M. In Zeolites: Facts, Figures, Fufure; Elsevier Science Publishers: Amsterdam, 1989. (23) Vigne-Maeder, F.; Auroux, A. J . Phys. Chem. 1990, 94, 316. (24) Nowak, A. K.; et al. J. Phys. Chem. 1991, 95, 848. (25) Flanigen, E. M.; Bennett, J. M.; Grose, R. W.; Cohen, J. P.; Patton, R. L.; Kircher, R. M.; Smith, J. V. Science 1978, 271, 512. (26) Evans, D. J.; Murad, S. Mol. Phys. 1977, 34, 327. (27) Allen, M. P.; Tildesley, D. J. ComputerSimulation of Liquids;Oxford University Press: Oxford, 1987. (28) Streett, W. B.; Tildesley, D. J.; Saville, G. in Computer Modeling of Matter, American Chemical Society: Washington, DC, 1978. (29) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; ALagona, G.; Salvatore Profeta, J.; Weiner, P. J.Am. Chem. SOC.1984,106, 765-784. (30) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J . Comput. Chem. 1983, 4, 187-217. (31) van Gunsteren, W. F.; Berendsen, H. J. C. Groningen Molecular Simulation (GROMOS)Library Manual; Biomos: Nijenborgh 16, Gronigen, The Netherlands, 1987. (32) Slater, J. C.; Kirkwood, J. G. Phys. Reu. 1931, 37, 682. (33) Scott, R. A.; Scheraga, H. A. J . Chem. Phys. 1966,45, 2091. (34) Hopfinger, A. J. Conformational Properties of Macromolecules; Academic Press: New York, 1973. (35) Leherte, L.; Lie, G. C.; Swamy, K. N.; Clementi, E.; Derouane, E. G.; Andre, J. M. Chem. Phys. Lett. 1988, 145, 237. (36) Grigoras, S.; Lane, T. H. J . Comput. Chem. 1988, 9, 25-39. (37) Clementi, E. Modern Techniques in Computational Chemistry: MOTECC-90; ESCOM: Leiden, 1990. (38) Kiselev, A. V.; Du, P. Q. J . Chem. SOC.,Faraday Trans. 2 1981,77, 1. (39) Kiselev, A. V.;Du, P. Q. J. Chem. Soc., Faraday Trans. 2 1981,77, 17. (40) Speedy, R. J. Mol. Phys. 1989, 66, 580. (41) Burkert, U.; Allinger, N. L. Molecular Mechanics; American Chemical Society: Washington, DC, 1982. (42) Hockney, R. W. Methods Comput. Phys. 1970, 9, 136. (43) Doherty, D.; Hopfinger, A. J. Macromol. 1990, 23, 676-678. (44) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNole, A.; Haak, J. R. J . Chem. Phys. 1984,81 (8), 3684-3690. (45) Price, S. L.; Harrison, R. J.; Guest, M. F. J . Comput. Chem. 1989, 10, 552. (46) Dunning, T. H.; Hay, P. J. In Modern Theoretical Chemistry; Schaefer, H. F., Ed., Plenum: New York, 1977. (47) Guest, M. F.; Harrison, R. J.; van Lenthe, J. H.; van Coder, L. C. H. Theor. Chim. Acta 1987, 71, 117. (48) van Gunsteren, W. F.; Karplus, M. J . Comput. Chem. 1980, 1 1 , 266-274. (49) Ouellette, R. J. J . Am. Chem. SOC.1972, 94, 7674. (50) Spackman, M. A. J . Chem. Phys. 1986,85,6579. (51) Papp, H.; Hinsen, W.; Do, N. T.; Bearns, M. Thermochim. Acta 1984, 82, 137-148. (52) Chiang, A. S.; Dixon, A. G.; Ma, Y. H. Chem. Eng. Sci. 1984, 39, 1461-1 468. (53) Stach, H.; Thamm, H.; JLnchen, J.; Fiedler, K.; Schirmer, W. In

1986. (15) Hayhurst, D. T.; Paravar, A. R. Zeolites 1988, 8, 27-29. ( 1 6) Paravar, A.; Hayhurst, D. T. In Proceedings of the 6th International

Proceedings of the 6th International Zeolite Conference, 1983. (54) Zibrowius, B.; Caro, J.; Krager, J. Z . Phys. Chem. 1988,269,1101. (55) Momany, F. A.; VanderKooi, G.; Scheraga, H. A. Proc. Natl. Acad. Sci. U.S.A. 1968, 61, 429. (56) Jobic. H.: BC. M.: Kearlev. G. J. Zeolites 1989.. 9., 312. (57) Trouw, F: R. Specfrochim: Acta, in press. (58) Wu, P.; Debebe, A,; Ma, Y . H. Zeolites 1983, 3, 118. (59) Impey, R. W.; Madden, P. A,; McDonald, I. R. Mol. Phys. 1982.46,

Zeolite Conference 1983. (17) Demontis, P.; Fois, E. S.; Gamba, A.; Manunza, B.; Suffritti, G. B. J . Mol. Struc. 1983, 93, 245-254.

(60) Van den Begin, N. G.; Rees, L. V. C. In Zeolites: Facts, Figures, Future; Elsevier: Amsterdam, 1989.

-4 1_7- .

Molecular dynamics simulation of propane and methane in silicalite

Recommend Documents