J. Phys. Chem. 1993, 97, 13742-13752
13742
Prediction of Adsorption of Aromatic Hydrocarbons in Silicalite from Grand Canonical Monte Carlo Simulations with Biased Insertions Randall Q.Snurr, Alexis T. Bell, and Doros N. Theodorou' Department of Chemical Engineering, University of California and Center for Advanced Materials, Lawrence Berkeley Laboratory, Berkeley, California 94720 Received: July 1, 1993"
Adsorption isotherms and isosteric heats of adsorption for benzene andp-xylene in silicalite have been calculated from molecular simulations. The simulations were performed using newly developed grand canonical ensemble Monte Carlo (GCMC) techniques in which insertion attempts are biased toward the most favorable regions of the zeolite pore space. The new techniques result in a substantial improvement in the efficiency of the simulations compared to traditional GCMC. The adsorption thermodynamics and molecular-level structure were studied for benzene and p-xylene in silicalite with Pnma symmetry (ORTHO) and P2,2121 symmetry (PARA). The subtle differences between ORTHO and PARA silicalite result in qualitatively different sorption behavior. An explanation of the experimentally observed step in the adsorption isotherm is presented, based on the results of the simulations and the ORTHO to PARA framework transformation that is observed experimentally. Predictions of the adsorption isotherms, isosteric heats, and siting locations of the adsorbates are in good agreement with experiment.
Introduction The adsorption and diffusion of aromatic sorbates in ZSM-5 zeolites are thought to be crucial to the para-selectivity displayed in processes such as xylene isomerization and methanol-to-gasoline conversion.l The importance of these systemshas prompted much experimental work on the fundamental adsorption properties of aromaticsorbates in zeolites of the ZSM-5 family. These studies have revealed certain unexpected features which are poorly understood. For example, adsorption isotherms for p-xylene in silicalite display a step at certain temperatures,*J and the heat of adsorption for benzene in silicalite exhibits an unusual dependence on sorbate 10ading.~~~ This paper discussesour efforts to use molecular simulations as a means for understanding these phenomena better and for predicting the adsorption properties of benzene and p-xylene in silicalite. Energy minimization studies and Monte Carlo simulations of benzene in silicalite at infinite dilution have been used by several authors to calculate sites of minimum energy within the zeolite lattice, the potential energies at those sites, and the heats of The studiesalso provide some insight into the diffusion mechanism and estimates of the activation energy for diffusion. To investigate adsorption at higher loadings, VignC-Maeder and Jobic' have computed interaction energies for 8 and 12 benzene molecules per unit cell. Minimum-energy locations and energetics for the closely related system of p-xylene in silicalite have been computed by Talu9and Pickett et a1.I0at infinite dilution. Higher loadings have been explored by Reischman and co-workers.lI de Vos Burchart et ale1* have used energy minimizations of a flexible silicalite lattice to investigate transformations in the zeolite structure upon adsorption ofp-xylene, and SchrliderI3has studied the differences between p-xylene adsorption in the monoclinic form of silicalite and in the high-loading orthorhombic form. Cheetham and Bull14recently reviewed the literature on computer simulations of p-xylene in silicalite. Molecular simulations of the adsorption equilibrium between the gas phase and the adsorbed zeolite phase for aromatics in silicalitehave been limited to the region of infinitedilution;Henry's constants for several aromatics in silicalite have been predicted by evaluation of configurational i n t e g r a l ~ . ~ J ~Isosteric -l~ heats
* To whom correspondence should be addressed.
Abstract published in Advance ACS Abstracrs, December 1, 1993.
0022-3654/93/2097- 13742304.00/0
of adsorptionand the equilibrium probability density of molecules at the various sorption sites have also been calculated. Energy minimization calculations have identified sorption sites in the channel intersections, the straight channels, and the sinusoidal channels but do not provide information concerning which of these sites would actually be occupied at finite temperatures due to the inherent neglect of entropic effects in the energy minimization technique. Evaluation of the configurational integrals predicts that, at low loading, benzene, toluene, and the three isomers of xylene all occupy the channel intersections preferentially rather than the straight or sinusoidal channels.l6 The conclusionsreached in the simulation studies of aromatics in silicalite thus far have been useful in understanding experimental data on these systems and provide some confidence in the ability of molecular modeling to predict adsorption properties a priori. However, the important problem of predicting the phase equilibrium behavior at higher loading is not addressed in the investigations published to date. Our objective here is to predict the adsorption isotherms and the isosteric heat of adsorption as a function of loading. We can then explore how sorbate siting within the zeolite varies with loading and how this affects the isotherms and isosteric heats. Because aromatic species fit very tightly in the pores of silicalite,some of the conventionalsimulation techniques for the prediction of sorption equilibria break down, and so an additional goal of this work has been to develop computational techniques for attacking these more difficult systems.
Model Silicalite is a completely siliceous MFI zeolite, having two sets of intersecting channels, each about 5.5 A in diameter. A set of straight channels run along they direction and are intersected by a set of sinusoidal channels which run along the x direction. One unit cell contains four straight channel segments, four sinusoidal channels, and four channel intersections. Silicalite is known to exist in three distinct forms: a monoclinic form, an orthorhombic form with Pnma symmetry, and an orthorhombic form with P212121 symmetry, which van Koningsveld et al. refer to as MONO, ORTHO, and PARA, respectively.ls-20 A reversible phase transition from the low-temperature monoclinic structure to ORTHO occurs at about 350 K. This transition can also be brought about by the adsorption of various organic 0 1993 American Chemical Society
Adsorption of Aromatic Hydrocarbons in Silicalite
The Journal of Physical Chemistry, Vol. 97, No. 51, 1993 13743
TABLE I: Parameters Used in Calculating the Zeolite/ Sorbate Lennard-Jones Interaction Potentials: Polarizabilities, a;van der Waals Radii, r0; Numbers of Valence Electrons,n~ Lennard-Jones c and tP atom C$,+ H
a(A3)
1.34[15] 0.40 [15]
P(A) 1.80 [15] 1.35 [IS]
nu
a
0.85 [24]
1.575 [24]
.(A)
c/k(K)
CH3 C
4 1
CH3 0
TABLE 11: Partial Charges for pxylene Used in Calculating tbe hlite/Sorbate Interaction Energy, Given in Units of the Electron Charge C -0,2835 C -0.2835
6
H 72 [25] 89.6
C H
3.923 [25] 2.806
The literature sources for all parameters are indicated.
molecules. the as-synthesized structure (containingthe template molecule) is also known to be ORTHO. The third form, PARA, is observed for silicalite with high loadings of p-xylene. In this paper we report calculations for both the ORTHO and PARA silicalite lattices. In the simulations the silicon and oxygen atoms of the zeolite framework are assumed to be fixed at the crystallographic coordinates determined from X-ray diffraction studies.1s+21Benzene and p-xylene are modeled as rigid, planar molecules, consisting of carbon atoms, hydrogen atoms, and methyl groups. Carbonxarbon, carbon-hydrogen, and carbonmethyl bond lengths are fixed at 1.40, 1.08, and 1.51 A, respectively. The zeolite and the sorbates are assumed to interact through a pairwise-additive potential between atoms of the adsorbed aromatic molecules and atoms of the zeolite lattice. The sitesite interactions are modeled with a Lennard-Jones plus pointcharge potential:
a
0.2025 -0.1161 0.1566 -0.1161 0.1566
CH3
0.2025 -0.1161 0.1566 -0.1161 0.1566
C
H C H
The resulting quadrupole moment equals -1.9 X
w - -
L -
m
w
esu cm2.
1
-
Figure 1. Contour plot of the calculated electric field along the center plane of a straight channel. Contour levels shown are at 0.1,0.2,0.3, 0.4,0.5,and 0.6 V A-I. Channel intersections can be seen at 1/4 and 3/4 of the distance along the y axis.
and simulations of benzene in zeolites NaXZ9and NaY .30 Partial charges off 0.15 are also in reasonableagreement where i a n d j are atoms of a sorbate molecule and of the silicalite with charges of fO.17 obtained from ab initio calculation^.^^ The lattice, respectively, and rijis the distance between them. Ai, and partial charges used for p-xylene are reported in Table 11. The Bi, are the Lennard-Jones constants, and qi and q, are the partial choice of the p-xylene values has been discussed earlier.I6 charges of the sites. A cutoff radius of 13 A is applied to the Partial charges are also placed on the atoms of the silicalite Lennard-Jones interactions, and the long-range electrostatic lattice: +2 on silicon and -1 on 0xygen.3~ These values were interactions are calculated using the Ewald summation techobtained from MNDO cluster calculations and are in close nique.Z2 The potentials between a sorbate atom and the zeolite agreement with more recent embedded cluster calculations using are tabulated on a three-dimensional grid throughout the density functional theory.33 A contour plot of the electric field asymmetric unit of the zeolite, and then during the simulations inside the zeolite resulting from these charges is shown in Figure the potential energy between a sorbate atom and the zeolite lattice 1. The electric fields are quite large, and a qualitative comparison is determined from the tabulated information using a threecan be made with the experimental value of 0.24 V A-I, as dimensional cubic hermite polynomial interpolation a l g ~ r i t h m . ~ ~ determined by infrared spectroscopyof adsorbed probe molecules The Lennard-Jones potential accounts for dispersive and in H-ZSM-5.34 repulsive interactions. The dispersion constant, Bij, for carbonThe zeolite/sorbate interaction model described here has been oxygen and hydrogen-xygen interactions is calculated from the used in earlier simulations of low-occupancy adsorption.I6 The Slater-Kirkwood equation using parameters from the literamodel yields Henry's constants and isosteric heats of adsorption ture,15J4 listed in Table I. The oxygen van der Waals radius was for benzene, toluene, and xylene that are in very good agreement adjusted by June et al.Z4to match the calculated Henry's constant with experimental data. For simulations at higher loadings, the for methane in silicalite with experimental data at a single interactions between sorbate molecules must also be included. temperature. The repulsiveconstant,A p is computed by requiring For these interactions we again use a Lennard-Jonesplus Coulomb the Lennard-Jones potential to have its minimum at a distance potential. The benzene parameters are taken from the work of equal to the sum of the van der Waals radii of the interacting Shi and Bartell,35who reparametrized a 12-6-1 Lennard-Jones atoms. For the methyl-oxygen interactions, Lennard-Jones e plus Coulomb potential to a 12-6-2-0 form, which is more and uvalues were taken from the l i t e r a t ~ r e .Following ~ ~ ~ ~ ~earlier convenient for large-scale computer simulations. The p-xylene work,l5qZ4 the dispersiveand repulsive interactions with the lattice interactions are modeled with the OPLS potential parameters of silicon atoms are neglected due to the very small polarizability Jorgensen et It should be noted that theCoulombinteractions of silicon, its small radius, and the recessed location of the silicon in the benzene/benzene potential are consistent with the values atoms away from the pore walls. of the partial charges used for the benzene/silicalite potential, Electrostatic interactions are included in an effort to account but the OPLS partial charges for the p-xylenelp-xylene interfor the quadrupole moments of the sorbates. For benzene, actionsare different from the partial chargesused for thep-xylene/ electrostatics are thought to have a significant effect on the liquid silicalite potential. and crystalline phases. Charges of -0.15 le1 and +0.15 le1 are placed on the carbons and hydrogens of benzene, giving a Metbodology quadrupole moment of -9.1 X l e z 6 esu cm2 for the modeled benzene, in good agreement with the range of experimental The most common technique for the prediction of zeolite values.26 These charges have been used previously in calculations adsorption phase equilibria from molecular simulations has been of benzene crystal structure^,^^ Monte Carlo studies of liquid grand canonical ensemble Monte Carlo (GCMC),374Salthough
Snurr et al.
13744 The Journal of Physical Chemistry, Vol. 97, No. 51, 1993
other techniqueshave been used, such as Monte Carlo simulations with a coupling parameter i n t e g r a t i ~ n ~or~molecular ,~’ dynamics simulationswith test particle insertions.43In a GCMC simulation three types of moves are performed. The first is a displacement and/or rotation step, handled by thenormal Metropolis method.22 In the second type of move, a new molecule is inserted into the system at a randomly chosen position. The new configuration is accepted with a probability
wherefis the gas-phase fugacity, Vis thevolume of the simulation box, N is the number of molecules present before the attempted insertion, k is Boltzmann’s constant, Tis the temperature, and A T is the change in the potential energy of the system due to the insertion. In the third type of move, a molecule is randomly chosen to be removed, and the removal is accepted with probability
p = min (l,N+exp[-z]J A T
(3)
GCMC simulations have- proven successful in predicting sorption thermodynamicsof a variety of simple sorbates in zeolites X,’7-39.42and ~ i l i c a l i t e . ~Attempts ~ , ~ ~ . ~to apply the normal GCMC technique to larger, more complex sorbates or to systems in which the sorbates fit very tightly in the pores are bound to be frustrated by the low acceptance rates of the insertion and deletion steps needed to sample the grand canonical ensemble correctly. Our initial attempts to perform GCMC simulations for benzene in silicalite illustrate this. A simulation at 70 OC and lo“ atm using 27 unit cells yielded an average loading of 0.27 molecules of benzene per unit cell, in agreement with the Henry’s constant calculated by direct evaluation of the configurational integralI6 and in agreement with the experimental isotherm.2 However, only0.3%of theattemptedinsertion and deletion moves were accepted during the simulation. This acceptance rate is high enough to yield proper ensemble averages at these conditions,48 but for simulations at higher loadings it will likely drop to an unacceptably low value, making GCMC infeasible. In normal GCMC, insertions are attempted uniformly throughout the volume of the simulation box. For zeolite systems, however, we know a priori that much of this volume is filled by atoms of the zeolite framework and is inaccessible to sorbate molecules. Furthermore, not all portions of the accessible space are equally favorable; there often exist preferred regions in which sorbate molecules are localized. This information can be incorporated into a GCMC simulation if insertions are not attempted randomly throughout the volume, but rather are performed such that sorbates are inserted preferentially into the most favorable locations in the zeolite according to some predetermined probability density p(r). Additionally, we can avoid insertion overlap with existing sorbate molecules in the system by searching for regions in which an additional molecule would fit and attempting insertions only in those regions. If the Monte Carlo insertion moves are biased in some manner, the acceptance rules for insertions and deletions must be altered to ensure that microscopicreversibility is satisfiedand that the grand canonical ensemble is still correctly sampled. Four biased GCMC schemes will now be described. The first scheme is the cavity-bias GCMC (CB GCMC) of In this method, insertions are only attempted in regions where the sorbate molecules leave “cavities” which can accommodate an insertion. This technique has been shown to increase the acceptance rates of insertion and deletion moves compared to normal GCMC. Asdeveloped by Mezei, however, it isappropriate only for liquid phases and does not incorporate any information about an adsorbent. The second scheme, energy-bias GCMC (EB GCMC), is a new bias scheme applied for the first time in this work. It exploits the knowledge of the structure of the zeolite and the known energetics experienced by sorbates at different
locations within the zeolite lattice. In this method, insertions are biased so that they are attempted more often in the energetically more favorable regionsof the zeolite. The tbird schemediscussed is a combined energy/cavity-bias GCMC (ECB GCMC) which biases the insertions according to the energetics of the zeolite but also tries to avoid overlaps with existing sorbate molecules. The fourth scheme, orientation-bias GCMC (OB GCMC), includes orientational biasing of insertions for systems with strong orientational ordering such as p-xylene in silicalite. Cavity-Bm GCMC. In Mezei’s cavity-bias GCMC method, translationand rotation moves are handled by the usual Metropolis rules, as in regular GCMC. An insertion is performed as follows. A number, Nlrial,of test points are randomly chosen throughout the simulation volume, V, and each point is tested to see whether it is at the center of a cavity of suitable radius, r,, with respect to the already existing sorbates in the system. A running average, Pc(N),is kept of the fraction of test points found to lie in suitable cavities; this provides a steadily improving estimate of the probability of finding a cavity of radius r, or larger in a configuration of N molecules. One of the suitable test points is chosen at random and the insertion is attempted there. The move is accepted with probability
A deletion is performed by randomly choosing one of the existing molecules and accepting its deletion with probability
For adsorbed phases, the Nlrialtest points can be chosen randomly from a list of points in the zeolite pores so that they are never in the impenetrable ’walls”, in which case the simulation volume, V, is replaced by the pore volume, V,, in eqs 4 and 5. Eaergy-BlrsCCMC. Allen and TildesleyZ2 discuss Monte Carlo schemes with preferential sampling; the discussion of energybias GCMC here will closely follow their ideas and notation. Let II,,,,,be the one-step transition probability of going from state m tostateninaMonteCarloalgorithm,andlet a,,,,be the probability of attempting to move the system from state m to state n. In a normal Metropolis scheme (such as GCMC), amn = a,,,,,and we have
where pn is the ensemble probability distribution for the system in state n. In a preferential sampling scheme such as EB GCMC, a,,,,,# an,,,and we have nmn
= a,,,,, anmpn 2 a,,,,,~,,,,m Z n
n#m
The Markov chain still obeys microscopic reversibility with am,, # a n m and is generated by making random trial moves from state m to state n according to the attempt probability am,,.The trial move is accepted with a probability given by min (1, anmp,,/ amnpm).
In the energy-bias GCMC scheme, translation and rotation moves are handled by the usual Metropolis rules as in regular
Adsorption of Aromatic Hydrocarbons in Silicalite
The Journal of Physical Chemistry, Vol. 97, No. 51. 1993 13745
GCMC. For insertions and deletions, we discretize the simulation volume into small cubic regions before the simulation begins, and to each cubelet we assign a weight, pi, normalized so that Xipi = 1. An insertion is performed by choosing a cubelet i according to the weighting pi and placing a molecule with its center at a randomly selected point within cublet i and with a random orientation. Themove is accepted with probability min (1, anmpn/ amnpm), where now
A deletion is performed by choosing one of the existing molecules at random, determining in which cubelet i its center resides, and removing it with probability min (1, anmPn/amnPm), where f f n m P n - pi-amnpm
V 'cubelet
NkT -
fV
(9)
Even though the deletions are not biased, the deletion acceptance probability must reflect the fact that thereversemoves (insertions) are biased. This is to ensure that microscopic reversibility is satisfied. The weightingpi should be chosen so as to increasethe efficiency of the Monte Carlo calculation. Note that if p i is chosen such that molecules are inserted uniformly throughout the volume V, i.e. pi = l/(number of cubelets), then eqs 8 and 9 revert to eqs 2 and 3, respectively,and theusual GCMC technique is recovered. If pi is chosen such that insertions are attempted uniformly but only throughout a portion of the zeolite defined as accessible (i.e., the pore space), then the method resembles the cavity-bias GCMC technique or the excluded-volume map sampling of Deitrick et except that the cavities are now defined with respect to a fixed adsorbent instead of mobile fluid molecules. To realize the full advantage of the EB GCMC technique, pi should contain information about the energetics of placing a molecule in cubelet i. In this work, we have defined pi as
microscopic reversibility, the move is accepted with probability
Orientation-Bias CCMC. From simple geometric considerations, it is clear that a long molecule such as p-xylene can enter the straight or sinusoidal channels only in certain orientations. Xylene molecules in the channel intersections might be expected to have more flexibility but also to be preferentially oriented along the x or y axis, corresponding to the sinusoidal or straight channels. It is desirable, therefore, to choose the orientation of an inserted p-xylene molecule such that it is aligned in an appropriate direction. The orientation of a rigid sorbate molecule with respect to the zeolite coordinate axes can be specified by three Eulerian angles, #,e, and t,b.22 and the desired orientations can be achieved by biasing the choice of 4 and t,b according to predetermined probability densities, p+ and ps. Note that p ) and ps are functions of position within the zeolite, so that molecules inserted in the channels tend to be oriented along those channels while molecules inserted in the channel intersections are preferentially oriented along either the x or y axis. In conjunction with ECB GCMC, orientation-bias GCMC can be implemented as follows. Insertions are handled as in ECB GCMC except that the orientations are not chosen randomly. Rather, 4 and are chosen from the probability densities p4 and ps, respectively, and the move is accepted with probability
+
(13) Deletions are performed by choosingone of the existing molecules at random and removing it with probability
Calculations where V i is the potential energy a test Lennard-Jones sphere would feel at the center of cubelet i due to the zeolite and j3 is I/kT. The Lennard-Jones sphere is only used to bias the insertions; in calculating AV in eqs 8 and 9, the explicit atomistic model for benzene or p-xylene described previously is used. The three-dimensional probability density at infinite dilution from the full atomistic model16 could also be used for the biasing. Energy/Cavity-BiasGCMC. The energy/cavity-bias GCMC scheme combines the ideas of Mezei's cavity-bias GCMC and the energy-bias GCMC introduced above. As in EB GCMC, the simulation volume is discretized before the simulation begins and each cubelet is assigned a weight pi. Translations and rotations are handled as before. An insertion move is begun by choosing Ntrial cubelets according to the probability distribution (pi) (rather than at random as in CB GCMC). A random point is chosen within each of the Ntriaj cubelets and tested to see whether it lies in a cavity of a suitable radius, rc, with respect to the sorbate molecules. A running average, P,'(N), is kept of the fraction of test points found to be in cavities at each occupancy, N. One of the suitable test points is chosen at random and the insertion is attempted there, giving the inserted molecule a random orientation. The move is accepted with probability
Deletion moves are not biased; one of the existing molecules is randomly chosen as the deletion candidate, and to preserve
Test Case: Argon in Sicalite. Before implementing any of the biased GCMC schemes for aromatics in silicalite, we first tested our algorithms with simulationsof thesimpler systemargon in silicalite. This system is simple enough that normal GCMC can predict the adsorption isotherm over the full range of sorbate occupancies. Three test conditions were chosen, corresponding to low, moderate, and high loadings. These calculations were performed to verify that all of the GCMC variations produced the correct results and also to assess the relative efficiencies of the various algorithms. The simulations were all performed for a temperature of 77 K with a system consisting of 27 unit cells. Each iteration consisted of a translation attempt followed by either an insertion or a deletion attempt (with 50% probability of each). The Lennard-Jones test sphere used to calculate Vi in eq 10 for biasing the insertions was simply the Lennard-Jones argon used in calculating the potential. The statistical uncertainties in the results were estimated by dividing each simulation run into 10 blocks and calculating the standard deviation of the block averages. Results are presented in Table 111. The first observation is that all of the GCMC algorithms yield the same results (within the statistical uncertainty of the calculations) for the number of atoms adsorbed and the average energy per atom. Although not shown, the fluctuations in N are also well reproduced by the biased sampling schemes. It can also be seen that, for a small simple sorbate like argon, the acceptance ratio for insertions and deletionscan be improved significantlyfor normal GCMC simply by restricting insertion attempts to the pore volume, V,,and then using V, instead of V in eqs 2 and 3. Note that there is some
Snurr et al.
13746 The Journal of Physical Chemistry, Vol. 97, No. 51, 1993
TABLE III: Comparison of CCMC Algorithm Efficiency for Test Case of Argon in Silicalite method
(atoms/uc)
WIN) (kJ 1mol)
GCMC GCMC GCMC CB GCMC EB GCMC ECB GCMC GCMC GCMC GCMC CB GCMC EB GCMC ECB GCMC GCMC GCMC GCMC CB GCMC EB GCMC ECB GCMC
2.25 f 0.04 2.26 f 0.04 2.25 f 0.02 2.26 f 0.03 2.26 f 0.04 2.26 & 0.04
-1 1.6 -11.6 -1 1.6 -1 1.6 -1 1.6 -11.6
7.86 & 0.08 7.87 & 0.06 7.88 f 0.03 7.8 0.1 7.92 f 0.08 7.92 & 0.07
-1 1.7 -11.7 -1 1.7 -1 1.7 -1 1.7 -1 1.7
23.61 f 0.09 23.5 f 0.1 23.57 f 0.07 23.57 f 0.07 23.6 f 0.1 23.3 f 0.1
-12.1 -12.1 -12.1 -12.1 -12.1 -12.1
(N)
VPl 0.23 0.12 0.12
Nrria~
5% of insertions and deletions accepted
relative efficiencyb
5.2 16.1 25.0 27.7 67.6 75.1
32.4 86.0 132.4 13.5 295.6 12.7
3.5 11.3 17.6 26.3 35.0 59.0
9.7 27.4 42.0 5.7 74.9 6.8
