J. Phys. Chem. 1990, 94, 1508-1516
1508
(3) Four states at 175,390,961, and 1123 K have been observed in the TD spectrum following the adsorption of CO on clean Re(0001). It has been shown that the two additional desorption features observed in previous work between 280 and 358 K correlate with a presence of surface impurities of C and 0. Adsorption of CO on Re(0001) at 115 K induces a 5 X 1 structure, whereas annealing this 5 X 1 phase to loo0 K produces a second square structure with a lattice constant 7 times that of the Re(0001) surface. This structure is incommensurate with respect to the substrate lattice and likely arises from C produced via dissociation of CO by annealing. (4)At low Cu coverages, TDS of Cu/Re(0001) following CO adsorption reveals a new desorption state that has a binding energy lying between that for CO on bulk Cu and that for CO on clean Re(0001). This state is presumed to arise from CO adsorbed on
-
a Cu overlayer perturbed by the underlying Re. At high Cu coverages (>5 ML), the TD spectra of CO are dominated by bulklike features of CO on Cu. (5) No new desorption states of nitrogen from overlayer Cu on Re(0001) compared to clean Re(0001) are observed; however, Cu is found to block Re sites for the adsorption and dissociation of nitrogen.
Acknowledgment. We acknowledge with pleasure the support of this work by the Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, and the Robert A. Welch Foundation. Registry No. Cu, 7440-50-8; Re, 7440-15-5; H, 1333-74-0; CO, N, 7727-37-9.
630-08-0
Prediction of Low Occupamy Sorption of Alkanes in SiHcalite R. Larry June, Alexis T. Bell, and Doros N. Theodorou* Department of Chemical Engineering, University of California, Berkeley, Berkeley, California 94720 (Received: April 12, 1989)
Alkane sorption in a pentasil zeolite has been studied through molecular simulations. A series of alkanes, including methane, n-butane, and three hexane isomers, were studied in silicalite by using a detailed atomistic representation that allows for torsion around skeletal bonds. Statistical mechanical principles have been employed to predict sorption equilibria at low occupancy. Henry’s constants and ismteric heats of sorption were calculated through the evaluation of configurational integrals with a Monte Carlo integration scheme. Results are in good agreement with experiment. The spatial distribution of sorbate molecules within the pore network, as well as perturbations to their conformation due to confinement in the pores, were determined via Metropolis Monte Carlo algorithm. Simulations of sorbate spatial distributions show that linear alkanes, such as n-butane and n-hexane, prefer to reside in the channels and avoid the channel intersections;on the contrary, bulky side groups in branched alkanes, such as 2- and 3-methylpentane, force these molecules toward the more spacious channel intersections. An analysis of sorbate conformations indicates that molecules are perturbed from the ideal gas dihedral angle distribution toward the more linear trans states when confined in silicalite.
Introduction Zeolites are crystalline aluminosilicates consisting of covalently bonded networks of Si04and A10, tetrahedra arranged in such a way as to form a network of connected cavities and channels with dimensions commensurate with molecular diameters (3-10 A). Because of their unique structure, zeolites have found wide application in processes such as adsorption, ion exchange, and cata1ysis.l” For the optimization of many of these applications, knowledge of the effects of zeolite structure and composition on adsorptive properties is essential. While such properties can be established experimentally, it would be highly desirable to develop theoretical methods for predicting zeolite adsorption behavior, so as to facilitate selection and optimization of zeolites for a given purpose. An attractive approach to this goal is the use of computer simulation based on an atomistic description of sorbate-zeolite interactions. Previous atomistic models of molecular sorption into narrow pores have been of two types: those dealing with idealized pore systems and those dealing with zeolites or similar cagelike structures. Simulation techniques, including molecular dynamics, grand canonical ensemble Monte Carlo, and Gibbs ensemble Monte Carlo among others, have been employed to study phenomena such as fluid structure, phase transitions, hysteresis effects, and wetting7-I4for simple Lennard-Jones fluids in idealized pores. Theoretical studies15J6 involving the application of mean field approximations have also been employed to examine the phase behavior of fluids in micropores.
* Author to whom correspondence should be addressed. 0022-3654/90/2094-1508$02.50/0
Simulation methods have also been used to investigate sorption in zeolites. Grand canonical ensemble Monte Carlo simulations have been employed to determine sorption equilibria over a full range of occupancies in model zeolite structures by Soto and ( I ) Riekert, L. In Advances in Catalysis; Eley, D. D., Pines, H., Weisz, P. B., Eds.; Academic Press: New York, 1970; Vol. 21; p 281. (2) Breck, D. W. Zeolite Molecular Sieves: Structure, Chemistry and Use; Academic Press: New York, 1978. (3) Zeolite Chemistry and Catalysis; Rabo, J. A., Ed.; ACS Monograph 171; American Chemical Society: Washington, DC, 1976. (4) Barrer, R. M. Zeolites and Clay Minerals as Sorbents and Molecular Sieues; Academic Press: New York, 1978. (5) Ruthven, D. M. Principles of Adsorption and Adsorption Processes; Wiley Interscience: New York, 1984. (6) Newsam, J. M. Science 1986, 231, 1093. (7) Snook, I. K.; van Megen, W. J . Chem. Phys. 1980, 72, 2907-2913. (8) Peterson, B. K.; Gubbins, K. E.; Heffelfinger, G. S.;Marconi, U. M. B.; van Swol, F. J . Chem. Phys. 1988.88, 6487-6500. (9) Panagiotopoulos, A. Z. Mol. Phys. 1987, 62, 701-719. (IO) Peterson, B. K.; Gubbins, K. E. Mol. Phys. 1987, 62, 215-226. (11) Tan, Z.; Van Swol, F.; Gubbins, K. E. Mol. Phys. 1987, 62, 1213-1224. (12) Schoen, M.; Cushman, J. H.; Diestler, D. J.; Rhykerd, C. L. J. Chem. P h y ~1988,88, . 1394-1406. (13) MacElroy, J. M. D.; Suh, S.-H. Presented at the 1987 Spring National Meeting of the A.I.Ch.E., Houston, TX, March 29 to April 2 1987. (14) Maada, J. J.; Tirrell, M.; Davis, H. T. J. Chem. Phys. 1985, 83, 1888-1901 .(15) Peterson, B. K.; Walton, J. P. R. B.; Gubbins, K. E. J . Chem. SO~., Faradav Trans. 2 1986.82, 1789-1800. (16fEvans, R.;Marconi, U.M. B.; Tarazona, P. J. Chem. SOC.,Faraday Trans. 2 1986, 82, 1763-1787.
0 1990 American Chemical Societv
The Journal of Physical Chemistry, Vol. 94, No. 4, 1990
Low Occupancy Sorption of Alkanes in Silicalite 5'1'1 I ' b k
____
1509
TABLE I: Bond Angles and Bond Lengths in Model Sorbate Molecules
1
a = 20 0 7 4 b = 1992A c = 1342A
a
t
c-c
Bond Lengths,
A
C-H
c-c-c
c-c-c C-C-H
MyersI7 and by Woods et a1.18 The work of Soto and Myers considered the sorption of hard spheres and Lennard-Jones spheres into a zeolite X cavity. By an adjustment of the potential parameters, Soto and Myers were able to fit the calculated sorption isotherm of krypton in zeolite X to experimental data. An approachI9 combining computer graphics and Metropolis Monte Carlo has yielded interesting results about the siting of methane at low pressures in zeolite Y. The Monte Carlo simulations showed that methane remains localized inside small regions of the zeolite at temperatures as high as 170 K. At higher temperatures, much more motion of the methane molecule was observed. Kiselev and co-workers2*22 have studied sorption thermodynamics in the Henry's law region for a number of inorganic, polar, and organic molecules in zeolites X, Y, and silicalite using configurational integration techniques. By representing all atoms of the sorbate molecules explicitly and employing a detailed atomatom potential expression, Kiselev and co-workers were able to calculate Henry's constants and isosteric heats of sorption for a large number of zeolite-sorbate systems. The present study has examined several aspects of the sorption of alkanes in silicalite. The specific alkanes considered were methane, n-butane, n-hexane, 2-methylpentane, and 3-methylpentane. Silicalite was chosen because in its ideal state this zeolite contains no aluminum and hence no cationic centers. This limits sorbate-zeolite interactions to those arising purely from dispersion and short-range repulsion potentials. Monte Carlo integration of the configurational integral was used to predict the Henry's constants and isosteric heats of sorption at very low occupancy. The distribution of sorbate rotational conformations for n-butane and n-hexane were calculated by using Metropolis Monte Carlo techniques. This technique was also used to characterize the distribution of sorbate molecules within the silicalite unit cell as a function of sorbate molecular weight and structure, as well as temperature.
Model Representation Silicalite is isostructural with ZSM-5 and contains virtually no aluminum atoms in the crystal lattice. The crystallographic coordinates of silicalite are available from detailed X-ray diffraction s t ~ d i e s . ~The ~ , lattice ~ ~ parameters of the unit cell, in (17) Soto, J. L.; Myers, A. L. Mol. Phys. 1981, 42, 971-983. (18) Woods, G. B.; Panagiotopoulos, A. Z.; Rowlinson, J. S. Mol. Phys. 1988, 63, 49-63. (19) Yashonath, S.; Thomas, J. M.; Nowak, A. K.; Cheetham, A. K. Nature 1988, 331, 601-604. (20) Beius, A. G.; Kiselev, A. V.; Lopatkin, A. A,; Quang Du, P. J . Chem. SOC.,Faraday Trans. 2 1978, 74, 367-379, (21) Kiselev, A. V.; Quang Du,P. J . Chem. Soc., Faraday Trans. 2 1981, 77, 1-15, 17-32. (22) Kiselev, A. V.; Lopatkin, A. A,; Shulga, A. A. Zeolires 1985, 5 , 261-267. (23) Olson, D. H.; Kokotailo, G. T.; Lawton, S . L.; Meier, W. M. J . Phys. Chem. 1981, 85, 2238-2243.
1.10
Bond Angle Supplements, deg linear chains (backbone) C-C-H methane H-C-H branching points in branched chains
Figure 1. Four silicalite unit cells are represented here in a 010 orientation. The straight channels of silicalite are clearly seen in this projection. The silicon atoms are only visible a t points of truncation from the rest of the lattice.
1.53
61
69
70.52 67 13.2
its orthorhombic form,23are reported as 20.07, 19.92, and 13.42
A. NMR and X-ray diffraction studies of s i l i ~ a l i t e clearly ~~-~~ indicate that the crystal structure undergoes a phase transition at approximately 350 K from a low-temperature monoclinic form to a high-temperature form that has an orthorhombic symmetry. This phase transition can also be induced at constant temperature by the presence of sorbate molecules within the pores of silicalite26 or by increasing the aluminum content of the lattice. In all simulations performed here, the phase with orthorhombic symmetry is employed. Silicalite, depicted in Figure 1, contains a channel system consisting of two sets of interconnected pores. One set of pores follows a linear path through the crystal; their cross section is elliptical, with principal axes of length roughly 5.4 and 5.6 A. A second set of pores follows a zig-zag path and has an elliptical cross sectionz3with dimensions approximately 5.1 A X 5.4 A. For the simulations performed here, the silicalite crystal is assumed to have a perfect lattice, be rigidly fixed in space, and contain no defects or aluminum atoms. As a consequence of the absence of aluminum and counterions, we neglect long-range electrostatic contributions to the sorbate-lattice potential energy. All sorbate molecules in our simulations are modeled explicitly as assemblages of carbon and hydrogen centers. Hard degrees of freedom, Le., bond lengths and bond angles which vibrate with small amplitudes around an equilibrium value, are kept fixed. In aliphatic molecules containing three or more skeletal carboncarbon bonds, rotation is allowed around all nonterminal skeletal carbon-carbon bonds. The torsion angles of terminal carboncarbon bonds are fixed with the methyl groups in a staggered conformation. This molecular representation is a significant improvement over that used by Kiselev and co-workers,22 who considered n-alkanes as rigid units with all carbon-carbon bonds placed in the trans state. The configuration of a sorbate molecule is specified by a vector. The first three elements consist of the chain start location r, which specifies the three Cartesian coordinates of the first carbon atom of a given sorbate molecule. A set of three Eulerian angles, 9 , specifies the orientation of the sorbate molecule with respect to a set of fixed coordinate axes. For molecules with internal rotational degrees of freedom (n-butane and larger), a set of n, 3 bond rotation angles, 4, must be specified, where n, is the number of carbon atoms in the backbone chain. In addition to the above degrees of freedom (6 for methane, 6 + (n, - 3) of n-butane and larger molecules), the geometry of a molecule is fixed by specifying bond lengths and angles for all carbon-carbon and carbon-hydrogen bonds. Table I lists the bond lengths and bond angles used in this study. Within a sorbate molecule, all pairs of atoms for which the distance between the two sites is not fixed by rigid bonds interact through a site-site potential. The zeolite and sorbate are also (24) Flanigen, E. M.; Bennett, J. M.; Grose, R. W.; Cohen, J. P.; Patton, R. L.; Kirchner, R. M.; Smith, J. V. Nature 1978, 271, 512-516. (25) Klinowski, J. T.; Carpenter, T. A,; Gladden, L. F. Zeolites 1987, 7 , 73-78. (26) Fyfe, C. A.; Kennedy, G. J.; Kokotailo, G. T.; Lyerla, J. R.; Flemming, W. W. J . Chem. Soc., Chem. Commun. 1985, 740-742. (27) Hay, D. G.; Jaeger, H. J . Chem. Soc.. Chem. Commun. 1984, 1433.
1510 The Journal of Physical Chemistry, Vol. 94, No. 4, 1990 TABLE II: van der Waals Radii (ro),Effective Numbers of Electrons ( n e ) ,and Atomic Polarizabilities (a)Used in Modeling All Interactions in the Silicrlite-Alkane Systems
atom
P, A
0
1.575 1.8 1.3
C
H
n. 7 .O 5.0
0.9
cy, A3
0.85
0.93 0.42
assumed to interact through an effective, pairwise-additive potential between atoms of the sorbate molecule and the zeolite lattice. In a silicalite-alkane system, this interaction is largely of a London dispersion type. All site-site interactions are modeled by a Lennard-Jones 6-12 potential
where B, is the dispersion coefficient, A , is the repulsion coefficient, and rij is the distance between interacting centers i and j . The dispersion coefficient B, is determined by the SlaterKirkwood relationshipz8
where a is an atomic polarizability, e is the charge of an electron, a. is the Bohr radius, and n, is an effective number of electrons
in an atom. The repulsion coefficient, Aij, is determined so that the potential has a minimum at distance equal to the sum of the van der Waals radii of the pair
A , = )/ZBij(Pi+ roj)6
(3)
where roi and flj are the van der Waals radii of the interacting atoms. To reduce the number of potential parameters, it was assumed, following Kiselev,z2 that the framework silicon atoms do not interact significantly with sorbate molecules. This assumption is justified because the sorbate molecule is surrounded by a latticework of framework oxygen atoms; silicon atoms, located at the centers of the Si04tetrahedra, are recessed from the outer edges of the pores. All uncertainty in the interaction potentials is absorbed into the oxygen parameters by this assumption. The parameters ne, a, and fl for carbon and hydrogen were taken from previous polymer simulations by Theodorou and S ~ t e r .The ~~ atomic polarizability of oxygen, a, was taken from Kiselev.22 An optimal value for the oxygen van der Waals radius was determined by comparing the computed Henry’s constant to the experimental data of Ma and mworkers” for methane in silicalite. This optimal oxygen radius was found to be 1.575 A, which is close to the value 1.52 8, recommended by BondL30 The potential parameters used for all the work reported here are summarized in Table 11. In modeling the intramolecular energetics of flexible sorbates, an additional term, associated with bond rotations around skeletal carbon-carbon bonds, is included. This “intrinsic torsional potential” is added to the nonbonded Lennard-Jones potential to determine the self-interaction potential energy and is modeled by Flory’s standard 3-fold symmetric e x p r e s s i ~ n ~ ~ , ~ ~ (4)
where
is the bond rotation angle measured relative to the trans
(28) Bondi, A. Physical Properties of Molecular Crystals, Liquids and Glasses; Wiley: New York,1968. (29) Theodorou, D. N.; Suter, U. W. Macromolecules 1985, 18, 14571478; 1986.19, 139-1s4,379-387. (30) Bondi, A. J . Chem. Phys. 1964, 68, 441-451. (31) Chiang, A. S.; Dixon, A . G.; Ma, Y. H. Chem. Eng. Sci. 1984, 39, 1451-1459, 1461-1468. (32) Flory, P. J . Statistical Mechanics of Chain Molecules; Wiley Interscience: New York, 1969.
June et al. position and K4 is a rotational barrier of 2.8 kcal/mol. The energetics of a dilute sorbatezeolite system can be divided into two contributions: sorbate intramolecular energy and zeolite-sorbate interactions. The overall potential energy of the zeolitesorbate system is then represented as the sum of these two components. The sorbate intramolecular, or self-interaction, potential describes energy variations arising from conformational changes in the sorbate molecule alone. The self-interaction potential is expressed as
and is only a function of the internal bond rotation angles, 4. The summation of Lennard-Jones terms is taken over all pairs of atoms in a sorbate molecule, t , for which the distances are not fixed by the structure of the molecule. The summation of intrinsic torsional potentials is taken over all nonterminal carbon-arbon backbone bonds. Note that the energy difference between trans and gauche states is included in the Lennard-Jones component of eq 5.29*32 For a rigid molecule, such as methane, @intra is zero. The potential energy arising from zeolite-sorbate interactions is obtained by a summation of pairwise contributions from all atoms, i, in a sorbate molecule, t , and all oxygen atoms in the zeolite lattice, z , which lie within a spherical cutoff radius of 13 8, from sorbate atom i:
In the evaluation of this sum for silicalite, the oxygen atoms in a symmetric block of 27 unit cells are tested to see if they fall within the potential cutoff radius of a given point in the primary unit cell. The positions of all atoms in the zeolite lattice are known exactly by virtue of the assumption of a rigid lattice. Consequently, the inner summation appearing in eq 6 can be evaluated to find the total energy experienced by the sorbate atom, i, at any position within the channels of the crystal. Such a pretabulation of the sorbate atom energies was carried out prior to the calculation of any thermodynamic properties, in the following manner. First, the Lennard-Jones potential experienced by a carbon or hydrogen center interacting with all oxygen atoms of the zeolite lattice was computed and summed on a three-dimensional grid of 0.2 8, spacing placed within the asymmetric unit of the unit cell. By virtue of the orthorhombic symmetry of silicalite, this computation needed to be performed only over one-eighth of the unit cell. Sorbate-zeolite interaction potentials needed for the evaluation of the configurational integrals were then interpolated at run time from the tabulated potential data with a three-dimensional cubic spline.33 This pretabulation scheme allows a 100-fold reduction in computational effort over the direct pairwise summation of zeolite-sorbate interactions at every step of the simulation. Calculations Our first objective was the prediction of sorption thermodynamics at very low occupancies. An important characteristic of this regime is Henry’s constant, Kh, measured in milligrams sorbate per gram zeolite per atmosphere. Considering the thermodynamic equilibrium between a pure gaseous sorbate at pressure P,and a zeolite crystal of total volume Vs at a temperature T, we define Henry’s constant as the slope of the low-pressure region of the sorption isotherm: (7)
In eq 7, na is the number of sorbate molecules present in the zeolite at equilibrium, Ma is the molecular weight of the sorbate molecule, NA is Avogadro’s number, and pz is the bulk zeolite density. By analyzing both the zeolite and the gas phase within the formalism of the grand canonical ensemble in the limit of very low sorbate (33) Sathyamurthy, N.; Raff, R. M.1 0 / 2 0 / 3 0 (QCPE 322); Department of Chemistry, Oklahoma State University: Stillwater, OK.
The Journal of Physical Chemistry, Vol. 94, No. 4, 1990 1511
Low Occupancy Sorption of Alkanes in Silicalite pressures and imposing the condition of equal sorbate chemical potentials in both phases, one can deduce the following expression for the Henry's constant (see Appendix) PzNA
-Kh(T) Ma
lexp(-P[@'(ro,f,4) =
6
+ @intra(4)]l dro d 0 d$
87r2Vs f exp(-P[ @intra(4)]) d4
where /3 = l/(kBT) and kB is the Boltzmann constant. Henry's constant emerges as the ratio of two configurational integrals. The integral in the numerator of eq 8 refers to a sorbate molecule at infinite dilution within the zeolite phase; the configurational integral in the denominator is for the same molecule in an ideal gas state where only the intramolecular potential is active. To evaluate the configurational integral in the numerator of eq 8, the potential energy experienced by the sorbate molecule must be integrated over all Eulerian angles, 9,all torsional angles, 4, of the sorbate molecule, and over all chain start positions, ro,within the crystal of total volume Vs. By virtue of the existing symmetries, the domain of integration of ro can be reduced to the asymmetric unit of the zeolite unit cell. The isosteric heat of sorption is defined as the difference between the molar enthalpy of the sorbate molecule in the gas phase and its partial molar enthalpy in the zeolite phase. The isosteric heat, under the same assumptions as eq 8 above, can be obtained (see Appendix) as follows (9)
@i"tra(4)]) dro d 0 d41-l (loa)
1aintra(4)
is an excellent method for the computation of averages with respect to a multivariate probability distribution function, such as those involved in the expression for Q,,(eq 9). Importance sampling is employed by the Metropolis algorithm to rapidly explore relevant configurations of the sorbate molecule in the zeolite unit cell. Apart from observable thermodynamic properties, one can explore equilibrium structural features of the sorbate-zeolite system with Metropolis Monte Carlo. Examples of such properties are the spatial distribution of sorbate molecules in the pores of the zeolite and the conformational changes in the sorbate molecule that occur upon sorption. To characterize the spatial distribution of sorbate molecules within the pores of the zeolite, a vector for the position of the center of mass, rcm,is defined. The expected density of centers of mass at a point r within the unit cell can be expressed as a configurational average over all sorbate degrees of freedom (see Appendix): Pcm(r;T) =
~ 6 [ r - r c m ( r o , * , 4 ) ~ ~ s ( r o , 0 ,dro 4) d 0 d4
PzNA
-S
h(T)
Ma
fis(ro,*,4)
Jf?(r0.~,4) dro d* d 4 = exp(-P[@Z(ro,*,4)
+ @i"tfa(4)11(1 1)
In eq 11, 6[r-rcm(ro,0,4)] is the Dirac delta function and rcm(r0,0,4)is a function that gives the center of mass position in terms of the microscopic degrees of freedom of the sorbate molecule. Physically, p,(r;T) is an intracrystalline concentration per unit pressure that reflects the siting preferences of the sorbate. At low gas-phase pressures, the product of pcm(r;T)and pressure tells how many molecules per unit volume of zeolite crystal we expect to find, on average, in the vicinity of position r within the unit cell, and all its crystallographically equivalent positions throughout the lattice. As a comparison of eq 8 and 11 can attest, by averaging pcm(r;nthroughout the entire crystal of volume Vs, one recovers the Henry's law constant:
exp{-P@intra(4))d 4
=
lexp{-P@i"tra(@)) d4
(lob)
Here again, the integral in eq 10a is over all degrees of freedom for a single molecule in the asymmetric unit of the unit cell. The multidimensional integrals involved in the expressions for the Henry's law constant and the isosteric heat of sorption must be evaluated numerically. Two methods of multidimensional integration, Gaussian quadrature34 and Monte Carlo3' integration, were compared. Application of both methods for methane showed that the Monte Carlo integration algorithm was approximately 20 times more efficient than Gaussian quadrature. On the basis of this result, the Monte Carlo integration scheme was used in the evaluation of all configurational integrals in this work. The standard Monte Carlo integration technique was enhanced for this problem by removing those regions of the zeolite that do not contribute accessible pore volume. The removal of unfavorable regions was accomplished as follows: Each trial was initiated with the random generation of the position and orientation of the sorbate molecule. As soon as the position of the chain start, ro, was chosen, the contribution to a2 from this first atom was evaluated by means of the spline interpolation scheme. If this interpolated contribution exceeded some preset limit, it was decided that the present trial constituted at attempt to "shoot" the molecule into a nonporous region of the unit cell: a statistical weight of zero was assigned to this outcome. Further generation of the sorbate configuration for this trial was terminated and a new trial was initiated. The isosteric heat of sorption was also evaluated by a second stochastic simulation technique, Metropolis Monte Carlo.36 This (34) Hornbeck, R. W. Numerical Analysis; Prentice Hall: Englewood Cliffs, NJ, 1975. (35) Press, W. H.; Flannery, B. P.; Teukolsky, S. A,; Vetterling, W. T. Numerical Recipes: The Art of Scientific Computing; Cambridge University Press: London, 1986.
The configurational average in eq 11 is obtained from a Metropolis Monte Carlo simulation of a single sorbate molecule carried out at infinite dilution within the zeolite unit cell, Le., with the probability density function fis(ro,0,4) in configuration space. In our Metropolis Monte Carlo simulation, three types of moves are attempted with equal probability. The first of the three moves is a simple translation; it is attempted as a random perturbation in the Cartesian coordinates of the chain start location. The second type of move involves a change in orientation of the sorbate molecule; it consists of a random perturbation in all of the three Eulerian angles. The first Eulerian angle, which corresponds to the polar angle of the first skeletal bond with respect to laboratory coordinates, is sampled from a cosine distribution so that orientationally unbiased sampling is ac~omplished.~'The third type of move consists of a random perturbation in one of the internal degrees of freedom, 4. Acceptance rates of 50% were obtained for each step by adjusting the maximum perturbation amplitude for all classes of moves. In the course of the simulation, the position of the sorbate center of mass is tracked. The number of times the molecule visits each cell of a fine three-dimensional grid running through the asymmetric unit of the silicalite unit cell is continuously updated. The ensemble average of 6[rrcm(r0,0,4)] is thus directly accumulated at each simulation temperature for each of the sorbate molecules studied. Sorbate conformation can be described by the distribution of rotation angles of the internal carbon-carbon bonds. The distribution of these bond angles is accumulated in our Metropolis Monte Carlo simulations via a method analogous to that used for (36) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.;Teller, A. H.; Teller, E. J . Chem. Phys. 1953, 21, 1087-1092. (37) Allen, M. P.; Tildesley, D. J . Computer Simulation of Liquids; Clarendon: Oxford, UK, 1987.
1512 The Journal of Physical Chemistry, Vol. 94, No. 4, 1990
-
June et al.
-3 Methylpentone
-CALCULATED A Expt: Chiong, 1984 10.0 . 0 Expt: POPP, 1984
1E8 -
e-.
E
h
1E7 -
E
0
s >
I
E
v
Y
Ost (expt) = 2 0 kJ/mol Ost (colc) = 18 kJ/mol
A
1 .o
3.0
2.5
3.5
/
-
1E5 /
1E4 -
‘
1000 2.0
4.0
/
3.0 1000/T ( K -
silicalite. The van der Waals radius of the lattice oxygen atoms was varied to achieve optimal agreement of the calculated Henry’s constant with the experimental value. The optimal value was 1.575 A. A logarithmic scale is employed for the ordinate.
3.5
4.0
)
for n-hexane, 2-methylpentane,and 3-methylpentane. Experimental data on the low-pressure sorption thermodynamics of 2-methylpentane and 3-methylpentanein silicalite are not available. 3.0 1 I
-Calculoted A Expt: Chiang, 1984
2.5 x
1E4 -
~
Figure 5. Comparison of predicted Henry’s constant and isosteric heat
1E6
-
3 Methylpentone Qst=63kJ/mol 2 Methylpentone QSt=63kJ/mol Hexane QSt=68kJ/mol
I 2.5
1000/T ( K - ’ )
Figure 2. Henry’s constant and isosteric heat of sorption for methane in
1E5
/
--Hexone
c
/ 2.0
1E6
I
/
/
2 Methylpentone
2.0
-Ideo1 Gas
t1
I
Sorbote
5 = 400 K
L Y
Ost (expt) = 47 kJ/mol Ost (colc) = 49 kJ/mol 10 2.0
2.5
3.0
3.5
4.0
-1.0
-0.5
0.0
0.5
1 .o
Dihedrol Angle (rod/rr)
Figure 6. Dihedral angle distribution for n-butane in silicalite at 400 K.
-Calculated 1 ~ .8 AExpt: Rees, 1986
/
h
E
I
5k
1E7
1E6.
E
v
II
~
1E5
-
+
Y
Ost (expt)=60-70 kJ/mol
1E4 1000 I 2.0
Ost (calc)=68 kJ/mol
I 2.5
3.0
3.5
is predicted to be greater than that of 2- or 3-methylpentane, as evidenced by the larger Henry’s constant and isosteric heat. The bulkier branched isomers are restricted in their motion within the zeolite when compared to n-hexane. No trend of deviation between calculated and experimental data is evident from the results. Predicted isosteric heats a t low occupancy for the lower molecular weight alkanes studied here can be correlated with the number of carbon units, m, in a sorbate molecule through a simple linear expression. The result, from a regression on the calculated isosteric heats is Q,, = 9Sm, 9.3 (13) where Qst is in kilojoules per mole of sorbate. Sorbate molecule structure is ignored in this correlation because too few data are available to accurately account for the effect of chain branching. Our predictions for low-pressure thermodynamic properties, obtained through configurational integration, are in good agreement with experiment. This indicates that the potential parameters used, which are identical across all systems studied, are reasonable. A direct comparison with the work of KiselevZ2 is not possible because of a different potential parametrization and the use of the monoclinic silicalite unit cell structure.23 It is noted, though, that the isosteric heats calculated by Kiselev are too large, resulting in Henry’s constants becoming too large at higher temperatures. The error in isosteric heats calculated by Kiselev increases with molecular weight. A comparison of Henry’s constants predicted by Kiselev at 303 K and experimental data shows that the Henry’s constant of methane is too small by about a factor of 3, while the predicted Henry’s constant for n-hexane is too large by about 2 orders of magnitude. The Metropolis Monte Carlo simulations reveal a variety of interesting structural information. Figures 6 and 7 show the conformational distribution for n-butane and n-hexane in silicalite, and in the ideal gas state. For butane sorbed in silicalite, a slight increase in the height of the trans peak, accompanied by a decrease
4.0
1000/T ( K - ’ )
Figure 4. Henry’s constant and isosteric heat of sorption for n-hexane
in silicalite.
the center of mass distribution. As the simulation progresses, a histogram of the Occurrence of a given bond angle within various intervals of conformation space is compiled; at the end of the simulation, the histogram is normalized to calculate frequency distributions. Results Predicted Henry’s constants and isosteric heats of sorption for methane, n-butane, and n-hexane are shown in Figures 2-4. While a comparison with experimental data is presented for methane, n-butane, and n-hexane,3’*Mexperimental data are not available for comparison with the calculated values of Kh and QSt for 2-methylpentane and 3-methylpentane. Comparison of Henry’s constants of 2-methylpentane and 3-methylpentane are made with that of n-hexane in Figure 5. The affinity of silicalite for n-hexane
1'he Journal of Physical Chemistry, Vol. 94, No. 4, 1990
Low Occupancy Sorption of Alkanes in Silicalite 3.0 2.5
1
2
1
I
1513
-Ideal Gas Sorbote T = 400 K
f I
'.oL&3J
3
D
1.5 2'o
r=
0.5 0.0
- 1 .o -1.0
-0.5
0.0
0.5
1 .o
Dihedral Angle (rod/n)
Y
Figure 7. Dihedral angle distribution for the central carbon-carbon bond of n-hexane in silicalite at 400 K. Note the qualitative difference o f hexane's conformation in the zeolite and its resemblance to butane's
conformational statistics.
Figure 9. Center of mass distribution for n-butane in silicalite at 400 K. See text for details.
Figure 8. Center of mass distribution for methane in silicalite at 400 K. The small size and low interaction energies of methane allow delocalization over much of the unit cell.
in the height of the gauche peaks, can be clearly seen relative to gas-phase butane. The constraining effects of the zeolite pore walls force the sorbate molecule into a more elongated shape than the one it assumes in the ideal gas state. In the gas phse, the trans state is populated by 53.3% of the molecules. In the zeolite, the trans peak accounts for 57.5% of the rotation angle distribution. Perturbations to the conformation of n-hexane are much more pronounced. As seen in Figure 7, the conformational distribution of the central carbon-carbon bond in the gas phase and the sorbate phase differ significantly. In the gas phase, 50% of the molecules exist with the central carbon-carbon bond in the trans state, while in the zeolitic phase 70%of the central carbon-carbon bonds are in the trans state. These calculations demonstrate, therefore, that Kiselev's assumption of all bonds being fixed in a trans conformation is not realistic. The spatial distributions for methane, n-butane, n-hexane, 2-methylpentane, and 3-methylpentane within silicalite at 400 K are shown in Figures 8-12. The siting distributions are represented by three-dimensional contour surfaces,38drawn with bold lines. Contours enclose those regions in space, in which the local (38) Wright, T. J. J . Compur. Graphics. 1979, 13, 182-189; subroutine the SCD graphics system of the National Center for Atmospheric
ISOSRF of
Research, Boulder. CO.
Figure 10. Center of mass distribution for n-hexane in silicalite at 400
K.
density of centers of mass, pcm,exceeds some preset level. These levels have been chosen in such a way that the center of mass of a sorbed molecule lies within the contoured volume with a probability of 50%. For reference, the pore structure of the silicalite unit cell is outlined in a lighter line. Siting distributions are displayed over half of the unit cell, the other half being symmetric with respect to a plane running parallel to the axes of the sinusoidal channels. The contours can be viewed as enclosing those regions in space in which the center of mass of a molecule sorbed at low occupancy spends 50% of its time. Results presented in Figures 8-10 show that methane and linear alkanes prefer to reside in the channels as opposed to the channel intersections. The diameter of the channels in silicalite is roughly 5.5 A and that of the channel intersections about 9 A. Evidently, there is a tendency for sorbate molecules to reside in the small diameter pores to maximize their attractive interactions with the zeolite. On the other hand, pore diameters that are too small will lead to conformational restriction, repulsion, and ultimately, the
June et ai.
1514 The Journal of Physical Chemistry, Vol. 94, No. 4, 1990
Figure 11. Center of mass distribution for 2-methylpentanein silicalite at 400 K. The preference for the channel intersectionsresults from the bulky structure of the branched species.
Figure 12. Center of mass distribution for 3-methylpentaneat 400 K. The molecule favors the channel intersections in a manner that is very similar to that of 2-methylpentane.
steric exclusion of the sorbate molecules from the pore. On the basis of these arguments, linear alkanes, which can assume elongated shapes and are roughly 4.5 A in diameter,39*40should prefer to reside in the channels and not in the larger intersections between two channels. It is also interesting to observe that, compared to n-butane, n-hexane prefers the somewhat narrower sinusoidal channels, over the larger straight channels. On the contrary, a singly branched alkane has a diameter of approximately 5.3 A39and, while not excluded from the interior of the channels of silicalite, is no longer comfortable in the channels. As seen (39) The approximate diameter of 4.5 A for linear alkane, calculated with the aid of the Alchemy molecular modeling package, compares well with the value of 4.3 A reported by Richards and R e e ~ . ’The ~ approximate value of 5.3 A for a singly branched alkane was also calculated with the Alchemy software package. (40) Richards, R. E.;Rees, L. V. C. Lungmuir 1987,3, 335-340.
Figure 13. Center of mass distribution for n-hexane in silicalite at 200 K. The strong localization effect due to temperature is clearly seen by
comparison to Figure 10. in Figures 11 and 12, the bulk of the center of mass distribution for 2- and 3-methylpentane is located in channel intersections. A comparison of the computed isosteric heats of 2-methylpentane and n-hexane shows that the former is 5 kJ/mol smaller; this indicates that the overall orientationally and conformationally averaged interaction with the zeolite is less favorable than that of n-hexane. The temperature dependence of the sorbate center of mass distribution is another interesting structural feature obtainable from our simulations. As the temperature decreases, the influence of attractive energetics dominates over the randomizing effect of entropy; the sorbate molecule will confine itself to the most energetically favorable regions and will not venture very far into the rest of the pore network. In the limit of zero temperature, the sorbate center of mass distribution at infinite dilution would consist of a few points, corresponding to the positions and orientations at which the total potential energy assumes a global minimum. In simpler terms, the infinitely dilute sorbate would “solidify” within the pores. To determine if n-hexane would ultimately localize itself preferentially within one of the channel systems, a series of Metropolis Monte Carlo simulations were undertaken at progressively lower temperatures. Results for n-hexane in silicalite at 200 K are presented in Figure 13. At 200 K, the center of mass distribution for n-hexane has shrunk to a set of very narrow regions which reside solely in the sinusoidal channel system. Comparison with n-hexane at 400 K shows that the spatial extent of the regions visited by n-hexane molecules is greatly restricted at 200 K. Note that the contours reflect the spatial extent, but not the density of the intracrystalline distribution of sorbate molecules. At 200 K the density of sorbate molecules is considerably larger than at 400 K. Conclusions We have presented calculations of low occupancy sorption of alkanes in silicalite with an improved sorbate model that includes internal degrees of freedom. This model, when coupled with a single set of reliable potential parameters, allows a number of alkane molecules and a range of properties to be investigated. The predicted values of Henry’s constant for methane, n-butane, and n-hexane are in reasonable agreement with experimental values, as are the predicted isosteric heats of sorption. Metropolis Monte Carlo calculations have produced very interesting results. The site at which a sorbate molecule resides within the unit cell is found to be a function of the shape of the sorbate molecule. Linear
The Journal of Physical Chemistry, Vol. 94, No. 4, 1990 1515
Low Occupancy Sorption of Alkanes in Silicalite alkanes are sorbed preferentially in the channels of silicalite, whereas branched alkanes are sorbed preferentially in the channel intersections. As a result of the constraining effects of the channels in silicalite, the conformation of n-butane is shifted slightly toward the trans configuration. This effect is much more pronounced, though, for n-hexane.
Acknowledgment. We thank the National Science Foundation for its support of this work through grant NSF-CBT-8718027. Additionally, we thank the San Diego Supercomputer Center for a generous allocation of computer time on the Cray X-MP/48. D.N.T. is grateful to the National Science Foundation for a Presidential Young Investigator Award, No. DMR-8857659. Appendix Statistical Mechanics of Low-Occupancy Sorption. The thermodynamics of sorption can be derived from microscopic structure and energetics through the formalism of the grand canonical ensemble. Previous treatments in the literature4J8 have focused on sorbates consisting of spherical, structureless molecules. Here, we present a formulation for flexible sorbate molecules that invokes the atomistic representation introduced in the second section of the paper. Consider a solid phase (S),consisting of a zeolite z and a sorbate a, at equilibrium with a gaseous phase (G) of pure sorbate a, at pressure P. The temperature of both phases is T. Let Vs be the volume of the solid phase and n, the number of sorbate molecules in the solid phase. Thermodynamic equilibrium dictates that
4
= 1: = M
(A. 1 )
from atomic Cartesian coordinates to internal coordinates. @,,htn represents the potential energy of sorbate-sorbate interactions between molecules i and j ; it is a function of the position and orientation of the two molecules and contributes to all configurational integrals beyond @(l,T,Vs). 0 ' and @intra have the meaning attributed to them in the text. In evaluating Do, sip, V ,and bond lengths and bond angles of the sorbate molecule are taken to be equal to their equilibrium values. Similarly, in evaluating W, all framework atoms are kept at the equilibrium positions about which they vibrate. Q, is the partition function of the pure zeolite solid in the absence of sorbate molecules. Note that, as a consequence of assumption 3, the zeolite component of the partition function reduces to a factor that is independent of the number and configuration of sorbate molecules. Any perturbations induced in the structure and thermal motion of the zeolite lattice by sorbate molecules are disregarded. Let us now examine the solid phase in the grand canonical ensemble. Temperature T, volume Vs, and (by the above arguments) the amount of solid present are fixed. The number of sorbate molecules is allowed to vary, subject to the requirement of a constant chemical potential, p. The grand partition function in the solid phase becomes ~ ( P , T ~ V= SQ~(T,VS)[~ ) + C$(nn~T,Vs)X".l = n,
where the fugacity, A, is defined as h = exp(0p). The number average of sorbate molecules in the solid phase under conditions T, Vs, and p is then43
where p is the chemical potential of the sorbate. The following assumptions are introduced in the statistical a In ES = mechanical formulation: ( n a )= -A ah 1 . All motion in both the solid and gaseous phases is treated classically; Le., each atom k contributes a term of the form pz/2mk to the kinetic energy of the phase of which it is a part, where P k and mk are the momentum and atomic mass, respectively. 2. The hard degrees of freedom (bond angles and bond lengths) of the sorbate molecules are constrained by stiff parabolic potentials. The force constants of these potentials are allowed to Note that the sums in both the numerator and the denominator approach infinity after all integrations over momentum space have will be polynomials of finite order in (qtA);repulsive sorbatebeen performed. (This corresponds to taking the infinitely stiff sorbate interactions will drive the configurational integrals limit of a "classical flexible" model representation, as opposed to @(n,,T,Vs) to zero at high loadings, q;the capacity of the zeolite invoking a classical rigid m 0 d e 1 . ~ ' ~ ~ ~ ) for sorbate molecules is finite. 3. All bonds and bond angles between covalently bonded atoms In the gas phase, the partition function for n molecules occuin the zeolite are governed by parabolic potentials. The stiffness pying volume V at temperature T is of these potentials is allowed to approach infinity after performing all integrations in momentum space. With these assumptions, 1 Q%,T,V) = ,l4t(T))"zC(n,T,V) (A.5) the canonical partition function for the solid phase b e c o m e ~ ~ ~ , ~ ~ 1 $(na,T, Vs) = QAT, Vs)-,lqt( na.
T)I"@(na, T,Vs) (A.2)
where
where
P ( n ,T, V) =
expI-P[C C.o,'ntcr(rgi,gi,di,roj,gj,4ji) + i>j
Z@int"(4i)]$ drio d q i ddi i
i- 1
Pressure in the gas phase can be expressed in terms of 1.1 and T using the formalism of the grand canonical ensemble43
p = - ln E G ( ~ , T , v ) PV In the above expression, B = l / ( k B T )and h is Planck's constant. The index k runs over all atoms in a sorbate molecule, each of mass mk. 14 is the determinant of an a X a matrix, the elements of which are the force constants for the hard degrees of freedom in the sorbate molecule. Do is the Jacobian of the transformation (41) GO, N.; Scheraga, H. Macromolecules 1976, 9, 535-542. (42) Honnell, K. G.; Hall, C. K.; Dickman, R. J . Chem. Phys. 1987,87, 664-614.
(A.6)
where EG(LL,T,I/)
= 1
1 + CQO(n,T,V)An= 1 + C,zG(n,T,V)(q,h)" n. n
('4.7)
In the thermodynamic limit (very large n and V), the right-hand (43) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976.
J . Phys. Chem. 1990, 94, 1516-1519
1516
side of (A.6) is idependent of V. Equations A.4 and A.6 provide a parametric representation of the equilibrium sorption isotherm
with parameter (q,X). Note that, by virtue of (A.l), this parameter is the same in both phases at equilibrium. In the limit of very low pressures (Henry's law region), the above formalism simplifies substantially. Sorbate-sorbate interactions become unimportant in both phases in this limit. Equations A.5 and A.6 reduce to their ideal gas counterparts: ZIG(n,T,V) = ( ~ T ' V ) ~lexp(-@@'"'"(+)) [ d+]"
(A.8) (A.9)
-
-
As seen from (A.9), the limit P 0 necessarily implies X 0. From eq A.4 one obtains, with NA Avogadro's number, pz the density of the zeolite crystal, Ma the molecular weight of the sorbate molecule, and Kh Henry's constant in weight sorbate per weight zeolite per pressure
lexp{-@[@z(ro,9.+)+
P
dro d* d+
8?r2Vslexp{-b@int"a(+)) d+
where f,,S(rOi.Wi&) represents the integrand appearing in the definition of P(n,,T,Vs) (see eq A.2). In the Henry's law region (X-+O), under conditions of equilibrium between gas and solid, eq A. 11 reduces to
The probability density appearing in the integral is thus fsl(ro,~,B)/zs(l,T,Vs). It is this function that is used in the Metropolis Monte Carlo calculation of molecular conformation and siting within the zeolite pore system. The isosteric heat of sorption is defined as a difference in partial molar enthalpies: Qst = -E (A.13) If the pressure is sufficiently low to warrant using the ideal gas equation of state in the gas phase, and if the partial molar volume of sorbate in the solid phase is considered negligible in comparison to the gaseous molar volume, one is led, by straightforward thermodynamic arguments, from eq A.13 to the Clausius-Clapeyron relation?
(A.lO)
For a solid phase of macroscopic dimensions, and in the absence of lattice imperfections and surface anisotropies, the ratio p (1,T, Vs)/ Vs remains the same if the integration over ro contained in the configurational integral is confined to the asymmetric unit of the zeolite unit cell. For the orthorhombic phase of silicalite, the three applicable symmetry operators are a mirror plane, a screw axis, and a point of inversion. By taking advantage of these symmetry operators, the relevant integration volume, Vs, is reduced to one-eighth of the unit cell. Equilibrium properties, such as sorbate conformation, are computed by averaging over the property in the distribution of the ensemble of interest. Let g(ro,*,+) be any function of the configuration of a molecule. The average value of the quantity g within the solid phase under given (p,T,V,) is, by definition
(A.14) By combining our grand canonical formulation eq A.4 and A.9 with eq A.14, it is possible to derive the following relation in the limit of very low sorbate coverages: (A. 15) In view of eq A.10, eq A.15 can readily be proved consistent with the following well-known thermodynamic relation, obtained from A.14 in the limit P 0:
-
(A.16) Registry No. SO2, 7631-86-9;methane, 74-82-8; butane, 106-97-8; hexane, 110-54-3; 2-methylpentane, 107-83-5;3-methylpentane,96-14-0.
Oscillatory Behavior of a Beiousov-Zhabotinskii Reverse Micelle System I. Gonda* and G . A. Rodley Department of Pharmacy, University of Sydney, Sydney, N S W 2006, Australia (Received: April 19, 1989; In Final Form: August 29, 1989) A previously reported reverse micelle system, formed by the incorporation of the components of the Belousov-Zhabotinskii (B-Z) chemical oscillator into a bis(2-ethylhexyl) sodium sulfosuccinate (A0T)-isooctane solution, has been studied in greater detail. It has been established that the AOT does not enter into chemical reaction with the B-Z components. This has enabled variations in features of the batch oscillation patterns to be directly interpreted in terms of perturbing influences. Studies have been made of the effects of varying the composition, diluting the aqueous or micelle components, and adding chloride to the system. The results show that the chemical integrity of the B-Z aqueous oscillator is retained on encapsulation into the reverse micelle system but with differences that may be attributed to sequestering of the B-Z reactants and products. Introduction We recently reported the preparation and properties of a reverse micelle-chemical oscillator system.I This displays a range of ( I ) Balasubramanian, D.; Rodley, G. A. J . Phys. Chem. 1988,92,5995.
0022-3654/90/2094-1516$02.50/0
interesting properties that may be related to the effect of combining spatial chemical ordering (micellar structure) with ordering of the temporal type (a chemical oscillator). The system consisted of the surfactant bis(2-ethylhexyl) sodium sulfosuccinate (AOT) in bulk isooctane combined with an aqueous mixture of Belou0 1990 American Chemical Society
