Table Of ContentAdsorption and Diffusion in Zeolites:
A Computational Study
Adsorption and Diffusion in Zeolites:
A Computational Study
ACADEMISCH PROEFSCHRIFT
terverkrijgingvandegraadvandoctor
aandeUniversiteitvanAmsterdam,
opgezagvandeRectorMagnificus
Prof.dr. J.J.M.Fransetenoverstaanvaneen
doorhetcollegevoorpromotiesingesteldecommissie,
inhetopenbaarteverdedigenindeAuladerUniversiteit
opdonderdag12oktober2000te12.00uurdoor
ThijsJosephHenkVlugt
geborenteGeleen
Promotoren:
Prof.dr.ir. B.Smit,UniversiteitvanAmsterdam
(cid:15)
Prof.dr. R.Krishna,UniversiteitvanAmsterdam
(cid:15)
Overigeleden:
Prof.dr.ir. A.Bliek,UniversiteitvanAmsterdam
(cid:15)
Prof.dr.ir. R.A.vanSanten,TechnischeUniversiteitEindhoven
(cid:15)
Prof.dr. D.Frenkel,UniversiteitvanAmsterdam
(cid:15)
Prof.dr. F.Kapteijn,TechnischeUniversiteitDelft
(cid:15)
Prof.dr. W.J.Briels,TechnischeUniversiteitTwente
(cid:15)
dr. Th.L.M.Maesen,ZeolystInternational(PQCorp)
(cid:15)
TheresearchreportedinthisthesiswascarriedoutattheDepartmentofChemicalEngineering,
FacultyofNatuurwetenschappen,WiskundeenInformatica,UniversityofAmsterdam(Nieuwe
Achtergracht166, 1018 WV,Amsterdam,TheNetherlands)withfinancial supportbythecoun-
cil for chemical sciences of the Netherlands Organization for Scientific Research (NWO-CW).
Averylarge partof thecomputerresourceshas beengenerouslyprovidedbySARA(Stichting
Academisch Rekencentrum Amsterdam) and EPCC (Edinburgh Parallel Computing Centre).
This thesis is also available on the web: http://molsim.chem.uva.nl. The author of this thesis
canbecontactedbyemail: [email protected].
ThisdocumentwasproducedusingLATEX.
(cid:15)
Printedby: Ponsen&LooijenBV,Wageningen.
(cid:15)
CoverdesignbyThijsJ.H.Vlugt. Inspiredby“thebluebook”ofFrankT.J.Mackey.
(cid:15)
Contents
1 Introduction 1
1.1 Zeolites . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 MolecularSimulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Scopeofthisthesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 Configurational-BiasMonteCarlomethods 7
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Dualcut-offCBMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Parallel CBMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3.2 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3.4 Resultsanddiscussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4 Generationoftrialsegmentsforbranchedmolecules . . . . . . . . . . . . . . . . . 18
2.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.6 AppendixA:Modeldetails . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.7 AppendixB:Proofofequation2.21 . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.8 AppendixC:Alternativeparallelalgorithm . . . . . . . . . . . . . . . . . . . . . . 23
2.9 AppendixD:Growthofisobutane . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3 Recoilgrowthalgorithmforchainmoleculeswithcontinuousinteractions 25
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 Descriptionofthealgorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2.1 Constructionofachain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2.2 Detailedbalanceconditionandacceptanceprobability . . . . . . . . . . . 27
3.2.3 ComparisonwithCBMC. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3.1 Simulationdetails . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3.2 EfficiencyofRGcomparedtoCBMC . . . . . . . . . . . . . . . . . . . . . . 31
3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.5 AppendixA:Alternativealgorithmtocomputetheweight . . . . . . . . . . . . . 34
3.6 AppendixB:Parallelization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.7 AppendixC:Fixedendpoints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4 AdsorptionofalkanesinSilicalite 41
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.3 Simulationtechnique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.4 Linearalkanes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
vi Contents
4.4.1 HeatsofadsorptionandHenrycoefficients . . . . . . . . . . . . . . . . . . 44
4.4.2 Adsorptionisotherms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.4.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.5 Branchedalkanes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.6 Fittingofsimulatedisothermswithdual-siteLangmuirmodel . . . . . . . . . . . 62
4.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.8 AppendixA:Alkanemodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.9 AppendixB:Discussionoftheexperimentaldata . . . . . . . . . . . . . . . . . . . 64
4.9.1 Heatsofadsorption . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
4.9.2 Henrycoefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5 Adsorptionofmixtures ofalkanesinSilicalite 69
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.2 MixtureIsotherms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.3 ConsequencesforDiffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
5.3.1 TheMaxwell-Stefantheoryforzeolitediffusion . . . . . . . . . . . . . . . 74
5.3.2 DiffusionofasinglecomponentinSilicalite . . . . . . . . . . . . . . . . . 75
5.3.3 Diffusionofbinarymixtures . . . . . . . . . . . . . . . . . . . . . . . . . . 75
5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
6 DiffusionofIsobutaneinSilicalite studiedbyTransitionPathSampling 81
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.2 TransitionPathSampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
6.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
6.2.2 MonteCarlosamplingfromthedistribution . . . . . . . . . . . . 84
6.2.3 TransitionStateEnsemble . . . . . . . . . . .F.(x.0.;T.). . . . . . . . . . . . . 86
6.2.4 Integratingtheequationsofmotion . . . . . . . . . . . . . . . . . . . . . . 87
6.3 Simulationandmodeldetails . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
6.4.1 Calculating thehoppingrate . . . . . . . . . . . . . . . . . . . . . . . . . . 89
6.4.2 Transitionstateensemble . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
6.6 AppendixA:Calculation ofafreeenergyprofile . . . . . . . . . . . . . . . . . . . 93
6.7 AppendixB:Bitwisetime-reversiblemultipletime-stepalgorithm . . . . . . . . . 94
6.8 AppendixC:Parallel tempering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.8.2 Applicationtotransitionpathsampling . . . . . . . . . . . . . . . . . . . . 97
6.8.3 Modelsystem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
Bibliography 101
Summary 109
Samenvatting(SummaryinDutch) 113
CurriculumVitae 117
Publishedwork 119
Acknowledgements 121
Chapter 1
Introduction
1.1 Zeolites
Zeolitesaremicroporouscrystallinematerialswithporesthathaveaboutthesamesizeassmall
molecules like water or n-hexane (pore size is usually - ˚). The structure of a zeolite is
3 12 A
basedonacovalentlybonded tetrahedrainwhichthetetrahedralatom isusuallySilicium
or Aluminum. The very famoTuOs4Lo¨wenstine rule only allows the existencTe of zeolites with a
Silicium/Aluminum ratioofatleast . Asallcornersofatetrahedralhaveconnectionstoother
1
tetrahedra,athreedimensionalporenetworkofchannelsand/orcavitiesisformed. Currently,
these are about different zeolite structures [1], several of these can be found in nature. To
100
clarify the topology of a typical zeolite, the pore structure of the zeolite Silicalite [2] is shown
in figure1.1. This zeolite hasa threedimensional networkofstraightand zigzag channels that
crossattheintersections.
Becauseoftheirspecialstructure,thereareseveralapplicationsofzeolitesinindustrialpro-
cessessuch as(selective) adsorption,catalysis and ion-exchange[3,4]. A recentexample ofthe
useofzeolitesisthecatalyticupgradingoflubricatingoils[5]. NoblemetalloadedAEL-typesili-
coaluminophosphatemolecularsievesselectivelyabsorbthewax-like,long-chainnormalparaf-
fins from an oil feed-stock and hydro-convert them selectively into branched paraffins [5–7].
CatalystsbasedonTON-[8–11]andMTT-type[5,8,11–13]zeolitescombineastrongaffinityfor
long-chain,normalparaffinswithasignificantlyhigherselectivityforhydro-isomerizationthan
forhydro-cracking[5–14]. However,themajorityofzeolitesthatisproducedworldwideisused
asion-exchangerindetergents.
Averyimportantcharacteristicofzeolitesistheadsorptionisothermofagivensorbate[15].
Anadsorptionisothermdescribestheamountofadsorbedmaterialasafunctionofthechemical
potentialat constanttemperature. Usingtheequation of stateof thesorbateoneis able tocon-
vert this chemical potential to the pressure [16]. At very low pressures, the amount of sorbate
willbenegligible. Theamountofadsorbedmaterialalsohasamaximum (athighpressure)be-
causethespaceforguestmoleculesin a zeoliteis limited. A verypopularequationtodescribe
adsorptioniszeolitesistheLangmuirequation:
(cid:18) kp (1.1)
max =
(cid:18) 1+kp
in which is the loading of the zeolite, the maximum loading, the pressure and a
max
(cid:18) (cid:18) p k
constant. For low pressures, there is a linear relation between the pressure and the loading
(Henry’slaw):
(1.2)
max
(cid:18) = k(cid:18) p=Kp
2 Introduction
Figure1.1: Porestructureof thezeolite Silicalite (MFI typeframework). Left: projection onthe
- plane. The straight channels are perpendicular to the - plane, the zigzag channels are in
x z x z
the - plane. Right: projectiononthe - plane. Thestraightchannelsare from toptobottom,
x z x y
the zigzag channels are from left to right. The pore size of the channels is slightly larger than
˚. The dimensions of the rectangular unit cell are ˚ ˚ ˚; multiple unit cells
5A 20:1A 19:9A 13:4A
areshown. Seealsofigure4.1foraschematicrepresentationofthiszeolite.
(cid:2) (cid:2)
in which is the Henry coefficient. The value of at different temperatures can often be de-
K K
scribedbytheintegratedformofthevan’tHoffequation
exp (1.3)
K=K0 [-(cid:1)U=RT]
inwhich isthetemperatureand thegasconstant. Themeasurementofadsorptionisotherms
T R
can be quite time consuming (see, for example, ref. [17] and chapter 4 of this thesis). As the
number of zeolite structures is rapidly increasing (see, for example, refs. [1,18,19]), to design
a new zeolite-based petrochemical process one will have to perform much experimental work
to find out which zeolite will be best. Therefore, it would save much time (i.e. money) if some
experiments could be replaced by fast computer simulations. Furthermore, molecular simula-
tions are able to simulate at conditions that are difficult to realize experimentally, for example,
at high temperatures or pressures, or multicomponent systems. Another advantage of molec-
ular simulations is that one is able to localize the positions of the molecules in the pores of a
zeolitedirectly. Thiscanprovideinsightinadsorptionmechanisms,forexample,theinflections
intheisothermsof - , - ,and - inthezeoliteSilicalitethathavebeenmeasuredexperi-
mentally[20,21]. FonrCan6 enxtCen7sivereivCie4wofcomputersimulationsoftheadsorption,diffusion,
phaseequilibriaandreactionsofhydrocarbonsinzeolitesthereaderisreferredtorefs.[22,23].
1.2 Molecular Simulations
In this thesis, we will use force field based computational methods. This means that we know
exactly all interactions between the atoms of our system. Once we know these interactions,
we are able to calculate a variety of static and dynamic properties like heats of adsorption,
adsorption isotherms, and diffusion coefficients. In general, there are two methodsto obtain a
molecularforcefield:
1.2MolecularSimulations 3
1. From quantum mechanical calculations. By solving the Schro¨dingerequationusing vari-
ous approximations, we can obtain forces betweendifferent atoms and molecules. These
forces can be fitted into a force field. This usually works very well for intra-molecular
bondedinteractionslikebond-stretching,bond-bending,andtorsioninteractions,butless
well for van der Waals interactions. Note that hydrocarbon-zeolite interactions are dom-
inated by van der Waals interactions (see, for example, ref. [24] and chapter 4). Recently,
there have been several quantum-mechanical studies of water and methanol in Sodalite
[25,26] usingtheCar-Parrinello technique[27].
2. From experimental data. A force field can be fitted in such a way that experimental data
likediffusioncoefficients,heatsofadsorption,orphaseequilibriacanbereproduced. This
forcefieldcanthenbeusedtocomputeotherpropertiesofothermolecules.
Oncewehaveaforcefield,wecancalculatedynamicandstaticpropertiesofoursystem. In
general,therearetwoclassesofmethods:
MolecularDynamics(MD).ThebasicconceptofMolecularDynamicsisNewton’ssecond
law,whichstatesthatthesecondderivativeofthepositionisproportionaltotheforce:
(cid:15)
x F
d2 i i (1.4)
=
2
dt mi
inwhich isthetime, isthemassofparticle ,F istheforceonparticle ,andx isthe
positionotfparticle . Tmhei velocityv isthetimedieriivativeoftheposition: i i
i
i
x
v d i (1.5)
i =
dt
Exceptforafewtrivialcases,theseequationscanonlybesolvednumericallyforasystem
ofmorethantwoparticles. Apopularalgorithmtosolvetheseequationsofmotionisthe
socalledvelocity-Verletalgorithm[28,29]:
F
x x v i(t) 2 (1.6)
i(t+(cid:1)t)= i(t)+ i(t)(cid:1)t+ ((cid:1)t)
2mi
F F
v v i(t+(cid:1)t)+ i(t) (1.7)
i(t+(cid:1)t)= i(t)+ (cid:1)t
2mi
inwhich isthetime-stepoftheintegration. Notethatthisalgorithmistime-reversible.
(cid:1)t
Theaverageofastaticproperty canbecalculatedfromthetimeaverageof :
A A
dtA(t) (1.8)
A =
h i dt
R
Animportantdynamicquantityistheself-diffusivity ,whichcanbecomputedbyeval-
R D
uatingthemean-squaredisplacement,whichreadsinthreedimensions
x 0 x 2
t+t - (t)
1 lim (1.9)
D = 0 (cid:28)(cid:12) (cid:16) (cid:17)0 (cid:12) (cid:29)
6 t !1 (cid:12) t (cid:12)
orbyevaluatingtheintegralofthevelocit(cid:12)yautocorrelation(cid:12)function
1
1 0 v v 0 (1.10)
D = dt (t) t+t
3 0 (cid:1)
Z D (cid:16) (cid:17)E
4 Introduction
A typical time-step for MD has to be smaller than any characteristic time in the system.
For molecular systems this is in the order of . This means that we have to
-15
(cid:1)t = 10 s
integratethe equationsof motion for stepsto performa simulation of our model for
15
10
one second. In practice, we are limited to simulations of due to the limitations of
-6
10 s
modern computers. This means that using straightforward MD, we cannot obtain static
and dynamic propertiesthat have a typical time-scale of or larger. A possible way
-6
10 s
to calculate the occurrence of such infrequent events for a special class of problems is
transitionstatetheory[30].
Monte Carlo (MC). In Monte Carlo algorithms, we do not calculate time averages but
phase space averages. For example, in the canonical ( ) ensemble, the average of a
(cid:15)
NVT
staticproperty isequalto
A
x x exp x
d A( ) [-(cid:12)U( )] (1.11)
A = x exp x
h i d [-(cid:12)U( )]
R
inwhichxresemblesthepositionofallparticlesinthesystem, isthetotalenergyofthe
R
U
systemand ,inwhich istheBoltzmannconstant. Becausetheintegralsin
equation1.1(cid:12)1 a=re1i=n(tkegBrTa)lsin manykdBimensions(usually at least100) andexp x is
[-(cid:12)U( )]
nearly always zero (i.e. only for a small part of x there is a contribution to the integral),
conventional numerical integration techniques are not suited to compute . Therefore,
A
the only suitable method is MC, in which the ratio of the integrals in equation 1.11 is
h i
calculatedinsteadoftheintegralsthemselves. InaMCsimulation,wegenerateasequence
(length )ofcoordinatesx ,insuchawaythattheaverageof canbecalculatedusing
i
N A
x x
i=N i=N
lim i=1 A( i) i=1 A( i) (1.12)
A = ; N 1
h i N!1 N (cid:25) N (cid:29)
P P
by ensuring that points x in phase space are visited with a probability proportional to
exp x . Thereisaniinfinitenumberofpossibilitiestogenerateasequenceofcoor-
dina[t-es(cid:12)xU(foir)]agivensysteminsuchawaythatthisequalityholds. However,tocalculate
accuriately,somemethodswillneedanastronomicallargenumberofstates(forexam-
A
ple, ),whileothermethodsneedonlyafewstates(forexample, ). Thisis
h i 500 5
N =10 N = 10
thecharmofMCmethods,becauseonehasthefreedomtomodifythealgorithmtoobtain
anoptimalefficiency. InMDsimulationsthereusuallyisnosuchfreedom.
Asimple MCmethodistheMetropolisMC method,inwhich kisgeneratedbyaddinga
random displacement in the interval to x . When a uniform distributed random
number between and is smaller th[-a(cid:1)n;e(cid:1)x]p i k x , we choose x k,
otherwisewe choo0se x 1 x . The maximum[-d(cid:12)is(pUla(ce)m-enUt( ic)a)]n be adjusted tio+1ob=tain
acertainfractionofaccie+p1te=dtriialmoves(usuallyaround %)(cid:1). Onecanprovethatinthis
50
methodthephasespacedensityofx isproportionaltoexp x forsufficientlylarge
[29,31,32]. i [-(cid:12)U( i)]
i
Forlongchainmoleculeswithstrongintra-molecularpotentialsthisalgorithmwillnotbe
veryefficientbecauseadisplacementofasingleatomwillnotchangetheconformationof
the molecule very much. Furthermore, there might be high energy barriers (for example
torsionalbarriers)whicharenotoftencrossed;thiswillleadtopoorsamplingstatistics. A
possible solution is the use of an algorithm that regrows a chain molecule completely or
partiallyandthuschangestheconformationofthemoleculesignificantly. Suchalgorithms
arediscussedinchapters2and3ofthisthesis.
