Adsorption 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.