Grand Canonical Monte Carlo Simulations of Nonrigid Molecules

3910

Langmuir 2000, 16, 3910-3919

Grand Canonical Monte Carlo Simulations of Nonrigid Molecules: Siting and Segregation in Silicalite Zeolite Amit Gupta, Louis A. Clark, and Randall Q. Snurr* Department of Chemical Engineering and Center for Catalysis and Surface Science, Northwestern University, Evanston, Illinois 60208 Received June 14, 1999. In Final Form: January 28, 2000 A method is described for performing grand canonical Monte Carlo simulations of molecules with internal degrees of freedom. The approach is applied to adsorption of methanol and cyclohexane in the zeolite silicalite. Calculated adsorption isotherms and heats of adsorption compare well with experimental data from the literature. In addition, the simulations provide information on preferential siting of molecules in the different pores of silicalite. Binary adsorption results are also predicted for cyclohexane/methane, benzene/methane, and benzene/cyclohexane mixtures. The siting preferences of the single components result in clear segregation effects for these binary systems. Benzene and cyclohexane display preferred siting in the channel intersections, while methane adsorbs preferentially in the channels. In the case of benzene/cyclohexane, the cyclohexane adsorbs in the intersections and the benzene is pushed to the zigzag channels.

Introduction Adsorption in microporous materials, such as zeolites and other molecular sieves, plays an important role in many commercial separations processes.1 In many additional processes it would be desirable to replace current technology with an adsorption-based separation to achieve energy savings or improved selectivity. Adsorption thermodynamics can also play an important role in zeolite catalysis.2 To develop such processes, many experimental, simulation, and theoretical studies have been carried out to elucidate the relationship between adsorbent and adsorbate structure and the resulting thermodynamics. Adsorbent heterogeneity is expected to result in preferential adsorption, or segregation, of different species in a mixture in the different regions of an adsorbent. This idea lies at the heart of many theories of adsorption thermodynamics, yet segregation is difficult to observe by direct experiment. Molecular simulation offers an attractive means to investigate siting and segregation effects in molecular sieves because a single simulation can yield both the macroscopic adsorption thermodynamics and the microscopic siting of the molecules. Adsorption isotherms and heats of adsorption have been calculated for a variety of systems using grand canonical Monte Carlo (GCMC) simulations.3-12 Many of these studies have focused on * To whom correspondence should be addressed. Fax: 1-847467-1018. Tel: 1-847-467-2977. E-mail: [email protected]. (1) Ruthven, D. M. Principles of Adsorption and Adsorption Processes; Wiley-Interscience: New York, 1984. (2) Haag, W. O. Stud. Surf. Sci. Catal. 1994, 84, 1375-1394. (3) Soto, J. L.; Myers, A. L. Mol. Phys. 1981, 42, 971-983. (4) Snurr, R. Q.; June, R. L.; Bell, A. T.; Theodorou, D. N. Mol. Simul. 1991, 8, 73-92. (5) Snurr, R. Q.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1993, 97, 13742-13752. (6) Cracknell, R. F.; Gubbins, K. E. Langmuir 1993, 9, 824-830. (7) Smit, B. Mol. Phys. 1995, 85, 153-172. (8) Dunne, J. A.; Myers, A. L. Adsorption 1996, 2, 41-50. (9) Van Tassel, P. R.; Davis, H. T.; McCormick, A. V. Mol. Simul. 1996, 17, 239-254. (10) Jameson, C. J.; Jameson, A. K.; Lim, H. M. J. Chem. Phys. 1996, 104, 1709-1728. (11) Heuchel, M.; Snurr, R. Q.; Buss, E. Langmuir 1997, 13, 67956804.

silicalite, which has a set of straight channels intersected by a set of zigzag channels. Several groups have seen that small molecules, such as methane or xenon, prefer the channels over the intersections.13,14 June et al.15 found that at low loading linear alkanes adsorb in the channels but branched alkanes prefer the roomier channel intersections. Aromatics such as benzene and xylenes are predicted to adsorb in the intersections at low loading16 and then to begin filling the zigzag channels as loading is increased.5 Smit and Maesen17 report a step in their simulated isotherms for heptane, which they attribute to a “commensurate freezing” where chains of the right length become localized in certain channel segments of the zeolite having a matching length. Binary mixtures of methane and tetrafluoromethane (CF4) in silicalite were investigated by Heuchel et al.11 Their GCMC simulations predicted mild segregation, with methane adsorbed preferentially in the zigzag channels and CF4 preferentially in the straight channels. Binary mixtures of other spherical molecules in silicalite, mordenite, and boggsite were investigated by Clark et al.12 All three of these zeolites exhibit heterogeneity due to the presence of pores of different sizes and shapes within a single zeolite. GCMC simulations predict pronounced segregation in mordenite and boggsite. For example, adsorption of a 50/50 gas-phase mixture of methane and CF4 at 10 atm and 400 K in mordenite results in 72% of the molecules in the side pockets being methane and 72% of the molecules in the main 12-ring channel being CF4. Clark et al.12 also found that entropic contributions play an important role in determining siting in these systems. (12) Clark, L. A.; Gupta, A.; Snurr, R. Q. J. Phys. Chem. B 1998, 102, 6720-6731. (13) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1990, 94, 8232-8240. (14) Nicholas, J. B.; Trouw, F. R.; Mertz, J. E.; Iton, L. E.; Hopfinger, A. J. J. Phys. Chem. 1993, 97, 4149-4163. (15) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1990, 94, 1508-1516. (16) Snurr, R. Q.; Bell, A. T.; Theodorou, D. N. In Proceedings from the Ninth International Zeolite Conference, Montreal, 1992; von Ballmoos, R., Higgins, J. B., Treacy, M. M., Eds.; Butterworth-Heinemann: Woburn, MA, 1993; Vol. II, pp 71-78. (17) Smit, B.; Maesen, L. M. Nature 1995, 374, 42-44.

10.1021/la990756f CCC: $19.00 © 2000 American Chemical Society Published on Web 03/24/2000

GCMC Simulations of Nonrigid Molecules

Langmuir, Vol. 16, No. 8, 2000 3911

The fact that segregation has been observed for simple spherical molecules suggests that such effects might be common in zeolite adsorption. These effects would be expected to be even more pronounced for molecules of differing shape. This paper extends previous studies of siting and segregation in silicalite to nonspherical molecules, in particular, methanol, cyclohexane, and benzene. The cyclohexane/methane and benzene/methane binary systems are also interesting as model systems for studying the effects of coadsorption of larger molecules on the dynamics of smaller penetrant molecules.18 In this paper single-component adsorption results, as well as binary mixtures with methane, are reported from atomistic GCMC simulations. We start by developing a general algorithm for GCMC simulations of molecules with internal degrees of freedom in the next section. This is followed by a description of the potential models and parameters used in the simulations. We then discuss the single-component and binary results.

Z)

Simulation Algorithm GCMC simulates a system at constant chemical potential, volume, and temperature. Because adsorption equilibrium requires constant chemical potential and temperature between the adsorbed phase and the gas phase, such simulations have become a preferred method for predicting adsorption phase equilibrium. During the simulation, the number of molecules fluctuates, and one calculates the average loading. From the chemical potential, which is equal in the gas phase, one calculates the gas-phase pressure from an equation of state, thus obtaining a point on the adsorption isotherm. GCMC simulation methods have been well described in the literature for monatomic species19 and for rigid molecules.4,5 Configuration bias (CB) GCMC has also been presented for simulations of chain molecules.7,20 Here we present an algorithm for performing GCMC of nonrigid molecules for which CB GCMC may not be necessary or appropriate, e.g., methanol or cyclohexane. We start by writing the grand canonical partition function for a single-component adsorbed phase at constant chemical potential µ, volume V, and temperature T following the discussion of June et al.15 The system is described by classical mechanics. As is common in molecular simulation, some “soft” intramolecular degrees of freedom are described by appropriate potentials, while other so-called hard degrees of freedom (e.g., bond lengths) are constrained by stiff parabolic potentials. The stiffness of the potentials for the hard degrees of freedom is taken to approach infinity after performing all integrations over momentum space. The contributions of the degrees of freedom of the zeolite itself to the partition function are neglected for the moment, and the zeolite is regarded as simply providing an external potential field felt by the adsorbed molecules. With these assumptions, the grand partition function for the adsorbed phase is ∞

Ξ(µ,V,T) )

∑ N)0

N

(qtλ) N!

Z(N,V,T)

(1)

where λ ) exp(µ/kT), k is Boltzmann’s constant, and N is the number of molecules in the system. (18) Gupta, A.; Snurr, R. Q. In preparation. (19) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, U.K., 1987. (20) Macedonia, M. D.; Maginn, E. J. Mol. Phys. 1999, 96, 13751390.

∫dr exp[-V kT ]

(2)

The integration is over the center of mass position of each molecule, the molecular orientations, and all nonstiff intramolecular degrees of freedom. The potential energy of the system V includes zeolite/sorbate and sorbate/ sorbate interactions, as well as intramolecular potential energy. The variable qt contains contributions from the kinetic energy and from the degrees of freedom that are taken as infinitely stiff in the model.15,21 It is convenient to introduce a set of scaled coordinates s19,22

[-V kT ]

∫

Z ) ΩN ds exp

(3)

∫dr

(4)

ΩN )

For a spherical molecule, Ω ) V, and for a rigid nonspherical molecule, Ω ) 8π2V after integrating the center of mass position and the three Eulerian angles that define the orientation of the molecule.19 A GCMC simulation probes the statistically important regions of configuration space according to the probability density distribution for the grand canonical ensemble: N

-V 1 (qtλ) N Ω exp F) Ξ N! kT

[ ]

(5)

F depends on all of the generalized coordinates s and the number of molecules N in a given configuration. Three classes of moves are performed in a GCMC simulation. In the first class, the number of molecules in the system is not changed, but some of the molecular coordinates are perturbed. For example, a randomly chosen molecule may be translated a small distance or one of its internal degrees of freedom may be perturbed. If the move attempts to take the system from state m to state n, the move is accepted with probability Fn/Fm, which from eq 5 is min {1, exp(-∆V /kT)}, i.e., the normal Metropolis acceptance.19 The second type of move is to insert an additional molecule into the system at a randomly chosen position. Again, the move is accepted with probability Fn/Fm, where now Nn ) Nm + 1. Thus

{

pins ) min 1,

qtλΩ ∆V exp Nm + 1 kT

[

]}

(6)

Similarly, in the third type of move a molecule is randomly chosen to be removed, and the deletion is accepted with probability

{

pdel ) min 1,

Nm ∆V exp qtλΩ kT

[

]}

(7)

The acceptance probabilities for insertions and deletions just presented are not in their usual form and not convenient as written. They are, however, general to any type of adsorbed molecule: spherical, rigid nonspherical, or nonrigid. To put them in a more convenient form, we note that for a gas phase in equilibrium with the adsorbed phase23 (21) Snurr, R. Q.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1994, 98, 5111-5119. (22) Maginn, E. J.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1995, 99, 2057-2079. (23) Hill, T. L. An Introduction to Statistical Thermodynamics; Dover: New York, 1986.

3912

Langmuir, Vol. 16, No. 8, 2000

Gupta et al.

kT ln ΞG(µ,V,T) V

P)

(8)

where

(qtλ)N ZG (N,V,T) Ξ (µ,V,T) ) N! N)0 ∞

∑

G

For an ideal gas15

PIG ) )

kT V

[

ln

∞

1

∑ (qtλΩ∫e-βV N)0N!

(9)

]

intra

dφ)N

∫

intra kT q λΩ e-βV dφ V t

(10)

where the integration is over the dimensionless scaled coordinates φ (a subset of s) corresponding to the intramolecular degrees of freedom and where β ) 1/kT. Rearranging, we obtain

qtλΩ )

PIGV

∫

kT e-βV

intra

dφ

pins )

{

intra PIGV ∆V dφ)-1 exp ( e-βV kT kT(Nm + 1)

∫

pdel )

{

[

NmkT intra ∆V dφ) exp min 1, IG ( e-βV kT P V

∫

[

]} (12) ]} (13)

For molecules without internal degrees of freedom, the integral is unity and the usual acceptance rules are recovered.4,5 The above acceptance ratios can be extended to ideal multicomponent systems. The resulting equations for the acceptance criteria are

pinsi )

{

min 1,

pdeli )

PIG i V

∫

intra

-βV i

( e kT(Nmi + 1)

dφ)

-1

[

∆V exp kT

]}

(14)

{

min 1,

NmikT V iintra ∆V ( e-β dφ) exp IG kT Pi V

∫

wi(n) )

[

]} (15)

where PIG i is the partial pressure of component i in the ideal gas phase, Nmi is the number of molecules of is the intramolecular energy of component i, and V intra i one molecule of component i. Expressions for the acceptance criteria for single-component real gases can be found in a recent paper by Macedonia et al.20 It is important to perform the insertion and deletion moves for each component with equal probability to preserve microscopic reversibility. To ensure this, we followed the algorithm outlined by Karavias et al.24 (24) Karavias, F.; Myers, A. L. Mol. Simul. 1991, 8, 51-72.

exp(-βV sz i (n))

(16)

Ncubelets

∑ j)1

(11)

and we can now write eqs 6 and 7 in their more useful forms by substituting for qtλΩ.

min 1,

In the GCMC algorithm described above, the insertion and deletion moves are made randomly. However, for molecules that fit tightly in the pores of an adsorbent, random moves will be, at best, inefficient and, at worst, impossible due to overlaps of the inserted molecules with the atoms of the adsorbent or with other molecules. To overcome this, a number of biased GCMC schemes have been developed for rigid molecules by Snurr et al.5 In this paper we use the energy biasing scheme. According to this scheme, the insertions are biased so that more attempts are made in the energetically favorable regions of the zeolite. This is done by discretizing the simulation volume into Ncubelets cubelets and then selecting a cubelet according to a weighting function. The weighting function, wi(n), is assigned to each cubelet n for each species i based on the zeolite/sorbate energy felt by a test Lennard-Jones sphere at the center of the cubelet. Our weighting function is defined as

exp(-βV sz i (j))

where V sz i (j) is the sorbate/zeolite energy of the test Lennard-Jones sphere for species i at the center of cubelet j. An insertion of sorbate i is attempted by picking a cubelet n according to its weight and then placing the molecule with a random orientation and internal conformation at a randomly selected position within the cubelet. The insertion is accepted with probability min (1, pinsi) where5

pinsi ) PIG V iintra ∆V 1 Vcubelet i V dφ)-1 exp ( e-β kT wi(n) V kT(Nmi + 1)

∫

[

]

(17) where Vcubelet is the volume of a cubelet. Because the insertions are biased, we must change the acceptance probability for deletion moves to preserve microscopic reversibility. Thus, the acceptance probability for a deletion move of a molecule of sorbate i from cubelet n is given by min (1, pdeli), where

pdeli ) wi(n)

NmikT V iintra ∆V ( e-β dφ) exp Vcubelet PIGV kT i V

∫

[

] (18)

Model and Calculations The formalism developed in the last section was used to generate isotherms of methanol and cyclohexane in silicalite. Furthermore, isotherms of cyclohexane/methane, benzene/methane, and benzene/cyclohexane mixtures were generated to study siting and segregation of the molecules in the different regions of the zeolite, as well as the adsorption thermodynamics. Silicalite is the siliceous form of the commercially important zeolite ZSM-5 (three letter designation MFI). The pore structure of silicalite is composed of a zigzag channel system running along the x direction and a straight channel system running along the y direction. The channels are about 5.5 Å wide, intersecting to form large channel intersections with a diameter of about 8.7 Å. Silicalite is known to exist in three different forms: a monoclinic form, an orthorhombic form with Pnma sym-


Langmuir, Vol. 16, No. 8, 2000 3913

Table 1. Potential Parameters for Methanol in Silicalite atom

SS charge

ss/k (K)

σss (Å)

SZ charge

sz/k (K)

σsz (Å)

H(OH) O C H(CH3)

0.418 -0.683 0.145 0.040

0.0 85.578 33.224 15.102

0.0 3.120 3.500 2.500

0.400 -0.537 -0.469 0.202

40.292 95.915 59.763 40.292

2.653 2.963 3.153 2.653

parameter

value

H(OH)-O (Å) C-O (Å) H(CH3)-C (Å) ∠H(OH)-O-C (deg) ∠O-C-H(CH3) (deg) ∠H(CH3)-C-H(CH3) (deg) H-C-O-H torsion (Vtor) (kJ/mol)

0.945 1.410 1.090 108.5 109.5 109.5 1.881

metry, and an orthorhombic form with P212121 symmetry. Van Koningsveld et al.25 refer to these as MONO, ORTHO, and PARA silicalite, respectively. Transitions between the phases can be induced by varying the temperature or by adsorption of various organic molecules. Methanol was modeled using a collection of parameters from the literature. Jorgensen et al.26 have used a charged all-atom representation of bulk methanol that we took as our starting point. Charges and Lennard-Jones potentials for each atomic center are listed in Table 1. The sorbate/ sorbate (SS) parameters used by Jorgensen reproduce methanol heats of vaporization within experimental error and the liquid density at 298 K within 0.89% of experiment.27 The sorbate/sorbate interactions are handled exactly as they were originally parametrized; truncations are applied at 11.0 Å based on center of mass distances. Both Lennard-Jones and electrostatic energies were smoothed from 10.5 to 11.0 Å by multiplication with a cubic weighting function that is unity at 10.5 Å and zero at 11.0 Å.28 First derivatives of the weighting function were set to zero at the end points. Long-range tail corrections were applied only to the Lennard-Jones methanol/methanol potentials. The truncation, smoothing, and tail corrections were all applied to be consistent with the bulk parametrization of the potentials by Jorgenson. Sorbate/zeolite (SZ) Lennard-Jones parameters for methanol were obtained by mixing the sorbate/sorbate potentials with the zeolite oxygen parameters given by Snurr et al.5 Following the common assumption that the zeolite silicon atoms are well shielded by the oxygen atoms, sorbate/silicon Lennard-Jones interactions were neglected. However, simulations conducted with the LorentzBerthelot mixed set of parameters underpredicted the Henry’s constant for methanol in silicalite. Consequently, the parameter for the zeolite oxygen was increased by 20% for the final simulations of methanol to produce better agreement of simulated and experimental Henry’s constants. The final and σ values for interaction between zeolite oxygen atoms and methanol atoms are given in Table 1 as sz and σsz. A different set of sorbate charges were used to calculate the sorbate/zeolite interactions via the Ewald summation.19,29 This was done with the understanding that the Jorgensen OPLS-AA charges are empirical parameters optimized to reproduce experimental energetics and short(25) Van Koningsveld, H.; Tuinstra, F.; van Bekkum, H.; Jansen, J. C. Acta Crystallogr. 1989, B45, 423-431. (26) Jorgensen, W. L.; Maxwell, D. S.; Tirado-Rives, J. J. Am. Chem. Soc. 1996, 118, 11225-11236. (27) Kaminski, G.; Jorgensen, W. L. J. Phys. Chem. 1996, 100, 1801018013. (28) Leach, A. R. Molecular ModellingsPrinciples and Applications; Longman: Harlow, U.K., 1996. (29) de Leeuw, S. W.; Perram, J. W.; Smith, E. R. Proc. R. Soc. London, A 1980, 373, 27-56.

range interactions. Such partial charges should only be regarded as loosely related to integrated electronic densities. The methanol sorbate/zeolite charges were chosen with a different philosophy, in hopes of staying closer to the rigorous definition of partial charge. They were taken from the density functional theory (DFT) cluster calculations of Gale et al.30 Note that, in contrast to the sorbate/ sorbate charges, they specify a negative charge on the carbon atom. Unlike the Lennard-Jones potentials, Coulombic interactions were calculated using both oxygen and silicon lattice atoms. Charges for lattice silicon and oxygen atoms were 1.42 and -0.71, respectively, and were taken from the DFT cluster calculations of Kyrlidis et al.31 They agree well with the periodic Hartree-Fock calculations for silicalite by White and Hess.32 To simplify the Monte Carlo moves for methanol, we kept all bond lengths and bond angles fixed at their equilibrium geometries. With this assumption, the only remaining degree of freedom is the methyl group rotation about the carbon/oxygen axis. Rotations of the methyl are hindered by the following torsion potential:26 3

Vtor )

Vtor

(1 + cos 3φi) ∑ i)12kT

(19)

Here Vtor is the torsion angle potential parameter for each of the three hydrogen-oxygen-carbon-hydrogen torsion angles. The torsion angles are zero in the gauche conformation and are denoted by φi. Because the torsion angle potential is relatively small, i.e., 5.65 kJ/mol in a fully gauche position, we did not bias the methyl rotations. In order to calculate the acceptance probabilities for insertions and deletions, the internal potential was numerically integrated over the scaled coordinate (see eqs 12 and 13). The result is a value of 0.4341 at 300 K. Cyclohexane was modeled as six methylene united atoms confined to chair and boat conformations. The twostate model of cyclohexane simplifies the integration in eqs 12 and 13 to a summation over the intramolecular energies of the two conformations. The GCMC simulation was carried out by attempting chair and boat insertions with equal probability, and deletions were attempted by picking an adsorbed molecule with any conformation. The moves were accepted according to the criteria given in eqs 12 and 13. The contributions to the intramolecular energy were normally modeled as arising from bond-stretching, bond-bending, and torsion terms.33 We keep the bond (30) Gale, J. D.; Catlow, C. R. A.; Carruthers, J. R. Chem. Phys. Lett. 1993, 216, 155-161. (31) Kyrlidis, A.; Cook, S. J.; Chakraborty, A. K.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1995, 99, 1505-1515. (32) White, J. C.; Hess, A. C. J. Phys. Chem. 1993, 97, 8703-8706. (33) Hendrickson, J. B. J. Am. Chem. Soc. 1961, 83, 4537-4547.

3914


Gupta et al.

Table 2. Potential Parameters for Benzene, Cyclohexane, and Methanea interaction cyclohexane CH2-CH2 CH2-O methane CH4-CH4 CH4-O benzene C-C H-H C-H C-O H-O C-Si H-Si cross terms CH4-C CH4-H CH4-CH2

A

B

C

31 814 286 5 858 696

8727.95 4042.66

Ryckaert et al.36 June et al.38

35 685 432 5 386 218

13250.69 4886.62

Goodbody et al.39 Goodbody et al.39

2 923 896 27 739 384 199 1 335 511 160 330

2383.62 77.82 469.70 1807.32 510.82

10 481 158 1 178 829 33 823 020

5692.83 1105.33 10774.64

35.10 35.10 -35.10

D

-1.40 -1.40 1.40

qi

-0.15 0.15 -0.15 0.15

qj

-1.0 -1.0 2.0 2.0

ref

Shi et al.40 Shi et al.40 Shi et al.40 Snurr et al.5 Snurr et al.5 Snurr et al.5 Snurr et al.5 mixing rules mixing rules mixing rules

a V ) A/r 12 + B/r 6 + C/r 2 + D + q q /r 2. V [)] kJ/mol; r [)] Å. The charge is in units of electron charge |e|. O and Si are zeolite lattice ij ij ij ij i j ij atoms.

lengths fixed at 1.54 Å, and the bond angles are at their equilibrium value of 109.5° in both the chair and boat conformations. Hence, only the torsion term contributes to the intramolecular energy, making the chair favored by 23.4 kJ/mol over the boat conformation.34 Following Rowley et al.35 the Lennard-Jones parameters for methylene/methylene interactions were taken from Ryckaert and Bellemans.36 For the sorbate/zeolite potential energy, only sorbate/oxygen interactions were considered as in the work of Kiselev et al.37 The methylene/oxygen parameters were taken from June et al.38 The parameters are listed in Table 2. Methane was modeled as a united atom Lennard-Jones sphere, and both sorbate/sorbate and sorbate/zeolite parameters were obtained from Goodbody et al.39 LennardJones parameters for methylene/methane interactions were obtained using the Lorentz-Berthelot mixing rules. Benzene was modeled as a 12-center molecule with partial charges on the carbons and hydrogens. Benzene sorbate/sorbate interactions are described by LennardJones plus Coulomb terms (12-6-1). This has been reparametrized in 12-6-2-0 form by Shi and Bartell,40 making it more amenable to computer simulations. The sorbate/zeolite potential parameters were obtained from Snurr et al.5 To calculate the benzene/methane interactions, values of and σ for carbon and hydrogen were obtained from the 12-6 terms of the 12-6-2-0 model and used with the usual mixing rules. The parameters are listed in Table 2. Before the start of the GCMC simulation, the sorbate/ zeolite energies were tabulated on a grid dividing the zeolite unit cell into cubelets with 0.2 Å edge lengths. During the simulation the sorbate/zeolite energy was calculated using a 3-dimensional Hermite polynomial interpolation algorithm.13 This scheme results in almost a 60-fold reduction in computational time. The pretabulated values were also used to bias the insertions into (34) Magalha˜es, F. D.; Laurence, R. L.; Conner, W. C. J. Phys. Chem. B 1998, 102, 2317-2324. (35) Rowley, R. L.; Ely, J. F. Mol. Phys. 1992, 75, 713-730. (36) Ryckaert, J.-P.; Bellemans, A. Discuss. Faraday Chem. Soc. 1978, 66, 95-106. (37) Kiselev, A. V.; Lopatkin, A. A.; Shulga, A. A. Zeolites 1985, 5, 261-267. (38) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1992, 96, 1051-1060. (39) Goodbody, S. J.; Watanabe, K.; MacGowan, D.; Walton, J. P. R. B.; Quirke, N. J. Chem. Soc., Faraday Trans. 1991, 87, 1951-1958. (40) Shi, X.; Bartell, L. S. J. Phys. Chem. 1988, 92, 5667-5673.

energetically favorable regions of the zeolite. The cyclohexane and benzene insertions were biased using a spherical benzene Lennard-Jones potential with parameters taken from ref 5. The methane insertions were biased using its own sorbate/zeolite potential map. For all results reported here, the gas phase was modeled using the virial equation of state truncated after the second term. The parameters were obtained from the critical properties using the Pitzer correlation.41 Single-component simulations were done using 3 × 3 × 3 unit cells, while the binary simulations were done in 2 × 2 × 2 unit cells. Production runs were 2-6 million steps long for each point on the isotherm, where each step consists of a perturbation attempt followed by either an insertion or a deletion attempt. Configurations were written to file every 400500 steps. To ensure faster equilibration, the last configuration of each run was used as the initial configuration for the next run. The loading was obtained by averaging over all of the stored configurations after discarding the equilibration section of the run. The statistical uncertainties were calculated by dividing the stored configurations into 10 blocks and calculating the standard deviation of the block averages. Results and Discussion Single-Component Systems. A series of methanol simulations were performed to generate an isotherm at 300 K. The results are shown in Figure 1. All methanol simulations were performed for the ORTHO form of silicalite. Because the potential model was adjusted to give good agreement with the data of Dubinin et al.42 in the Henry’s law region, the simulation results match experiment very well up to about 5 kPa. Above 5 kPa the simulations agree better with the NaZSM-5 data from Hunger et al.43 Given the variability between different samples, it is difficult to say if further adjustment of the potential model would be justified. The higher loadings predicted by the simulation at higher pressures may only be due to the defect-free structure of the simulated lattice. Regardless, it is encouraging to note that good agreement (41) Perry, R. H.; Green, D., Eds. Perry’s Chemical Engineer’s Handbook; McGraw-Hill: New York, 1984. (42) Dubinin, M. M.; Rakhmatkariev, G. U.; Isirikyan, A. A. Izv. Akad. Nauk SSSR, Ser. Khim. 1990, 1950-1950 (English translation, orginally submitted in 1988). (43) Hunger, B.; Matysik, S.; Heuchel, M.; Einicke, W. Langmuir 1997, 13, 6249-6254.


Langmuir, Vol. 16, No. 8, 2000 3915

Figure 1. Comparison of simulated methanol isotherm at 300 K with experiment.42,43 Below 5 kPa, simulation results match the data from Dubinin et al.42 very well.

Figure 2. Comparison of isosteric heats for methanol in silicalite from simulation and experiment.44,45

with experiment can be had using liquid-phase potentials mixed with lattice oxygen parameters from the literature. Only a small adjustment was needed to match the isotherm in the Henry’s law region. Isosteric heats of adsorption were also calculated from the methanol simulation data using eq 2011 where R is

qst ) RT - 〈V 〉 - 〈N〉

( ) ∂〈V 〉 ∂〈N〉

(20)

T,V

the ideal gas constant, T is the temperature, N is the number of molecules, V is the potential energy per molecule, and V is the volume. Figure 2 is a comparison of the simulated isosteric heat of adsorption with experimental data from Pope44 and from Akhmedov et al.45 The high experimental heats at low loading are due to small concentrations of hydroxyl groups in the structure. Because the simulation lattice is defect-free, the simulated heats do not show this feature. All three curves show an increase in heat of adsorption with increased loading. This is typically thought to be a consequence of increasingly tight packing, giving rise to increased sorbate/sorbate interactions. Separation of the simulation energies into (44) Pope, C. G. J. Chem. Soc., Faraday Trans. 1993, 89, 11391141. (45) Akhmedov, K. S.; Rakhmatkariev, G. U.; Dubinin, M. M.; Isirikyan, A. A. Izv. Akad. Nauk SSSR, Ser. Khim. 1987, 1587-1591 (English translation, orginally submitted in 1985).

Figure 3. Comparison of cyclohexane isotherms from GCMC simulations with experimental data of Cavalcante et al.46 The open symbols connected by dashed lines represent GCMC simulations of cyclohexane fixed in the chair conformation. The open symbols connected by solid lines represent the two-state GCMC simulations of cyclohexane. The closed symbols indicate the experimental data of Cavalcante et al.

sorbate/sorbate and sorbate/zeolite contributions (not shown) supports this idea. The disagreement between simulation and experiment at higher loadings may be due to the ability of the defect-free simulation lattice to hold more molecules. Single-component cyclohexane isotherms were generated at 423, 473, and 573 K to compare against the experimental data of Cavalcante et al.46 These simulations were done in the ORTHO form of silicalite. Two sets of simulations were carried out at each temperature: first with cyclohexane in only the chair conformation and second with cyclohexane free to adopt both the boat and chair conformations. Figure 3 shows the plots of the two different sets of simulations along with the experimental data. We find a marked improvement in the agreement with experiment when we go from the single conformation simulation to the two-state simulation. The isosteric heat of adsorption for the two-state cyclohexane was calculated to be 59.5 kJ/mol at 473 K at a loading of about 1 molecule/ unit cell. This is lower than the heat reported in ref 46 by about 4 kJ/mol. Despite the lower heat of adsorption, the simulated loadings are consistently higher than the experimental loadings. This can probably be attributed to pore blockages and crystal defects which are present in the experimental zeolite samples but are absent from the simulated structures. These defects would lower the effective adsorption volume, giving a lower overall loading per unit cell. Based on the relative energies of the two conformations, the gas phase is comprised of 0.1-0.7% boat conformations depending on temperature. This ratio increases slightly in the adsorbed phase, where we find that 2-6% of the total loading is now boats. To understand why this is so, we calculated the energy of adsorption and change in entropy between the gas phase and the adsorbed phase from the Henry’s law constants at the three different temperatures. This can be done by expressing the Henry’s law constant as a ratio of the canonical partition function in the adsorbed phase, QSNVT, to that in the gas phase, 15 QG NVT, as described by June et al. (46) Cavalcante, C. L., Jr.; Ruthven, D. M. Ind. Eng. Chem. Res. 1995, 34, 177-184.

3916


Khi ) lim Pf0

〈nai〉

QSi (N)1,Vs,T) )β G PVs Q (N)1,Vs,T)

Gupta et al.

(21)

where Khi is the Henry’s constant in molecules/unit volume/ pressure, P is the pressure, and 〈nai〉 is the number of molecules of species i adsorbed in volume Vs. The gasphase and the adsorbed phase partition functions are both evaluated over the same volume Vs.22 Rearranging yields

-

( )

Khi 1 1 1 ln ) - ln(QSi ) + ln(QG) β β β β

Table 3. Calculating the Energy and Entropy of Adsorption from the Henry’s Law Constants temperature (K)

Khchair (mol/m3/kPa)

Khboat (mol/m3/kPa)

423 473 523

3040 455 37

29.3 8.87 0.774

∆U ∆S

-57.1 kJ/mol -0.058 kJ/mol/K

-35.98 kJ/mol -0.046 kJ/mol/K

(22)

Now, using A ) U - TS ) -(1/β) ln(Q), we get

-

( )

Khi 1 ln ) ∆Ui - T∆Si β β

(23)

where ∆ refers to the adsorbed phase minus the gas phase. Using eq 23, we can back out the change in entropy and internal energy from Henry’s constants at different temperatures. For our purposes here we calculate Henry’s constants from the simulation results discussed above for chair molecules and boat molecules, treating them as different “species”; i.e., we divide the number of molecules adsorbed of each conformation by the total gas-phase pressure to get the Henry’s constant for that conformation. The gas phase consists of a mixture of cyclohexane boats and chairs distributed according to their relative Boltzmann weights. The cyclohexane chair conformation is taken to be the reference energy state, giving an average gas-phase energy of about 10-3 kJ/mol. Because the energy difference for both the chair and boat conformations is calculated with respect to this reference gas phase, the internal energy of the boats includes 23.4 kJ/mol of torsion energy. In other words, ∆ for adsorbed chairs includes only the change upon going from the gas phase to the adsorbed phase. However, ∆ for most (more than 99%) adsorbed boats includes conversion from the gas-phase chair state to an adsorbed boat conformation. The Henry’s constants for the cyclohexane chair and boat conformations along with the change in energy and entropy are given in Table 3. Though the boats face a big energetic penalty in both phases because of the torsion energy, the entropic contribution upon adsorption works in their favor. The reduction in entropy upon adsorption for the boats is less than that for the chairs, making the percentage of boats slightly higher in the adsorbed state than in the gas phase. We also investigated the singlet distribution of cyclohexane to understand its siting. There is some controversy as to the preferred location of cyclohexane in silicalite. Aliev et al.47 propose that cyclohexane sits in the ZSM-5 intersections from their 2H NMR investigations of deuterated cyclohexane. In a recent paper Ashtekar et al.48 suggest that cyclohexane prefers the straight channels based on their Fourier transform Raman spectra. To obtain the preferred siting of cyclohexane, we discretized the zeolite unit cell into 0.2 × 0.2 × 0.2 Å cubelets. During the course of the simulation, we kept track of the number of times the center of mass of a molecule fell within each cubelet. We then calculated the probability of observing the molecule within each cubelet. This probability grid is then used to generate an “equal-probability” isosurface. (47) Aliev, A. E.; Harris, K. D. M. J. Phys. Chem. A 1997, 101, 45414547. (48) Ashtekar, S.; Hastings, J. J.; Gladden, L. F. J. Chem. Soc., Faraday Trans. 1998, 94, 1157-1161.

Figure 4. Singlet distribution of cyclohexane molecules in silicalite at 473 K and 10.68 kPa. The loading of cyclohexane is 3 molecules/unit cell. The light gray outline is the accessible volume of the zeolite. The darker gray volume is the region where the probability of observing the center of mass of the cyclohexane boat is 99%. Black is indicative of the region where the probability of finding a cyclohexane chair conformation is 99%.

The isosurface encloses the volume where the molecules spend 99% of the time. One such plot for cyclohexane at 473 K is given in Figure 4. We find that cyclohexane prefers the relatively open intersections over the tight channels. Because silicalite has 4 intersections/unit cell, this suggests a saturation loading of 4 molecules/unit cell, which is consistent with Figure 3. The single-component adsorption of benzene in silicalite was studied with GCMC simulations by Snurr et al.5 To provide a context for the binary systems described below, we provide a brief summary of their major results. Snurr et al. performed simulations of benzene in both ORTHO and PARA silicalite. To explain the experimentally observed steps in the adsorption isotherms, they proposed that adsorption of high loadings of benzene induces the same ORTHO to PARA phase transition that is observed experimentally upon adsorption of p-xylene.25 This transformation is thought to occur around 4 molecules/unit cell. The simulations showed that the zeolite phase change will, in turn, induce a sharp increase in the amount of benzene adsorbed, explaining the step in the experimental isotherms. Below 4 molecules/unit cell, benzene adsorbs almost exclusively in the channel intersections of ORTHO silicalite. Between 4 and 8 molecules/unit cell, benzene adsorbs in both the intersections and the zigzag channels.5 Binary Systems. Simulations of binary mixtures of cyclohexane/methane and benzene/methane were done at 420 K and a total fixed gas-phase pressure of 500 kPa. At lower pressures very little methane adsorbs when cyclo-


Figure 5. Loading of methane and benzene at 420 K as a function of the gas-phase composition. The total pressure is kept constant at 500 kPa. The open triangles represent the total loading, the filled circles the loading of benzene, and the filled squares the loading of methane.

Langmuir, Vol. 16, No. 8, 2000 3917

Figure 7. x-y diagrams of the three different binary systems. The open circles represent the mole fraction of benzene in the benzene/methane system (420 K and 500 kPa), the filled squares indicate the cyclohexane mole fraction in the cyclohexane/ methane system (420 K and 500 kPa), and the filled triangles represent cyclohexane in the cyclohexane/benzene system (328 K and 1 kPa).

Figure 6. Loading of methane and cyclohexane at 420 K as a function of the gas-phase composition. The total pressure is kept constant at 500 kPa. The open triangles represent the total loading, the filled circles the loading of cyclohexane, and the filled squares the loading of methane.

hexane or benzene is present. As described above, the adsorption of benzene is thought to induce a phase transition in silicalite from the ORTHO form to the PARA form; hence, all benzene/methane simulations were performed in PARA silicalite, because the benzene loadings were above 4 molecules/unit cell at these conditions. The ORTHO form of silicalite was used for the cyclohexane/ methane simulations. Because of the relatively small percentage of boats in the pure-component systems, the cyclohexane/methane simulations were done with the cyclohexane fixed in the chair conformation. The isotherms of the benzene/methane and cyclohexane/ methane systems are shown in Figures 5 and 6, respectively. Although the initial loading of cyclohexane increases more steeply than that of benzene, the final loading of benzene is higher. This can be attributed to the fact that cyclohexane adsorbs more strongly than the benzene, giving rise to a greater initial slope. However, because benzene is smaller in size than cyclohexane, the entropic factor becomes favorable for benzene at higher loadings. Also the PARA form of silicalite has a slightly higher volume than the ORTHO form, which further contributes to benzene having a higher loading. Both cyclohexane and benzene adsorb more strongly than methane. This results in much higher loadings of cyclohexane and benzene than methane even at very low partial pressures of cyclohexane and benzene. This suggests that silicalite could potentially

Figure 8. Site distribution of methane and benzene at 420 K and gas-phase partial pressures of 250 kPa each. The loading of benzene is 7.4 ( 0.4 molecules/unit cell, and the loading of methane is 0.36 ( 0.02 molecules/unit cell. The faint gray outline is the accessible volume of silicalite. The darker gray is the region where the probability of observing a methane molecule is 99%, and the black represents the region of 99% probability of observing benzene.

be a useful sorbent for the separation of methane from its mixtures of cyclohexane and benzene. Similar trends can be observed in the x-y plots of Figure 7. The cyclohexane increases steeply at first and then continues to rise more slowly, whereas the benzene has a smaller initial slope but reaches a plateau much more quickly. Clark et al.12 have defined four models to classify binary adsorption in zeolites having multiple types of sites. These models are based on the preferred siting of the pure components making up the binary system and the actual siting seen in the mixture. If, in the mixture, the two species do not exhibit any segregation, filling each site in proportion to its volume, then the system follows the volume-filling model. If both components have different preferred sites as single components and they retain this preference in the mixture, then the mixture is said to conform to the reconciliation model. The competition

3918


Gupta et al.

Figure 10. Loading of cyclohexane and benzene as a function of the gas-phase composition at 328 K and a constant total pressure of 1 kPa. The open triangles represent the total loading, the filled circles the loading of cyclohexane, and the filled squares the loading of benzene. Figure 9. Site distribution of methane and cyclohexane at 420 K and gas-phase partial pressures of 250 kPa each. The loading of cyclohexane is 5.1 ( 0.1 molecules/unit cell, and the loading of methane is 0.86 ( 0.02 molecules/unit cell. The faint gray outline is the accessible volume of silicalite. The darker gray is the region where the probability of observing a methane molecule is 99% and the black represents the region of 99% probability of observing cyclohexane.

classification is used to describe systems where both components prefer the same site as single components, and they compete for the same preferred sites in the mixture. Sometimes the presence of one species serves to push the second one even more into its favored site. This

system is said to obey the accommodation model. This classification scheme can be used to describe the results obtained in this work. We can get information about the siting and segregation of benzene/methane and cyclohexane/methane mixtures by looking at the site distribution plots given in Figures 8 and 9. Benzene and cyclohexane can be found in both the channels and intersections as expected from the loading. Methane, however, is excluded from the intersections. This is not surprising because we know from previous work49 with pure-component systems that methane prefers the channels to the intersections and both benzene and cyclohexane prefer the relatively open intersections. Hence, both benzene/methane and cyclo-

Figure 11. Site distribution of benzene and cyclohexane at 328 K and gas-phase partial pressures of 0.5 kPa each. The loading of cyclohexane is 4 molecules/unit cell, and the loading of benzene is 0.9 molecules/unit cell. The faint gray outline is the accessible volume of silicalite. The darker gray is the region where the probability of observing a cyclohexane molecule is 99%, and the black represents the region of 99% probability of observing benzene.


hexane/methane conform to the reconciliation model. The degree to which methane is excluded from the intersections is higher for cyclohexane than benzene because the methane molecules can come closer to the benzene because of its smaller size. Also, the molecular size affects the degree of localization of the molecules. Comparing the extent of the black region in Figures 8 and 9, we find that cyclohexane is more localized than benzene. Our simulations show that this segregation is independent of the gas-phase composition. An important point to note from the above systems is that when heterogeneity is considered, it is not sufficient to consider just the adsorbent but rather adsorbent/ adsorbate combinations. While a given adsorbent may appear homogeneous to one sorbate, it can appear quite heterogeneous to another. For example, silicalite is relatively homogeneous to small spherical molecules, and only mild segregation is observed for mixtures of methane and CF4.11 On the other hand, silicalite is not at all homogeneous for benzene adsorption, and mixtures of methane and benzene display strong segregation effects. A more interesting system perhaps is one where both species have the same preferred sites. A mixture of cyclohexane/benzene has this characteristic because both species prefer the intersection sites. This system is described by the competition model. Simulations of the cyclohexane/benzene system were performed at 328 K with variable gas-phase compositions, while the total gas-phase pressure was kept constant at 1 kPa. Note that these conditions are different from those of the binary systems above. These simulations were done at lower total pressures to facilitate better equilibration by keeping the loadings of benzene and cyclohexane lower. At these conditions the loading of benzene is less than 4 molecules/ unit cell; hence, the simulations were carried out in ORTHO silicalite. The total loading varies between 4.6 and 5.2 molecules/unit cell, as can be seen from Figure 10. Because only one molecule can occupy an intersection at (49) Snurr, R. Q.; Ka¨rger, J. J. Phys. Chem. B 1997, 101, 6469-6473.

Langmuir, Vol. 16, No. 8, 2000 3919

a given time, there is competition between the molecules to occupy the 4 coveted intersection sites. We can see from Figure 7 (filled triangles) that cyclohexane adsorbs more strongly than benzene. On the basis of this observation, we would expect cyclohexane to displace the benzene from the intersections. This can indeed be seen from the site distribution plot given in Figure 11, where cyclohexane still occupies the intersections and benzene is pushed into the zigzag and straight channels. We find that the site distribution of benzene and cyclohexane remains relatively unchanged over the entire gas-phase composition range. Summary and Conclusions We have outlined a method for performing GCMC simulations of molecules with internal degrees of freedom. We have used this method to adjust the potential parameters of methanol to match experimental data in silicalite. Furthermore, this method was extended to simulate flexible ring molecules. By just going from a system with one conformation to two conformations of cyclohexane, we saw a marked improvement in the agreement with experimental data. We also studied the siting of the cyclohexane in silicalite and found that the molecules preferred the relatively open intersections to the channels. Binary systems were studied to understand the siting and segregation of the molecules. We found that in both the cyclohexane/methane and benzene/methane systems benzene and cyclohexane prefer the intersections whereas methane favors the channels. These systems fall under the reconciliation classification of Clark et al.12 The benzene/cyclohexane system obeys the competition classification, because both species prefer the intersections. In this case the cyclohexane adsorbs more strongly, pushing benzene into the zigzag channels. Acknowledgment. This work has been supported by the National Science Foundation CAREER program and by the Petroleum Research Fund. LA990756F

Grand Canonical Monte Carlo Simulations of Nonrigid Molecules

Recommend Documents