Table Of ContentProspects of Transition Interface Sampling simulations for the theoretical study of
zeolite synthesis
∗
Titus S. van Erp, Tom P. Caremans, Christine E. A. Kirschhock, and Johan A. Martens
Centre for Surface Chemistry and Catalysis, K.U. Leuven, Kasteelpark Arenberg 23, 3001 Leuven, Belgium
7 (Dated: February 6, 2008)
0
The transition interface sampling (TIS) technique allows to overcome large free energy barriers
0
within reasonable simulation time, which is impossible for straightforward molecular dynamics.
2
Still,themethoddoesnotimposeanartificialdrivingforce,butitsurmountsthetimescaleproblem
n byanimportancesamplingoftruedynamicalpathways. Recently,itwasshownthattheefficiencyof
a TIStocalculatereactionratesislesssensitivetothechoiceofreaction coordinatethanthoseofthe
J
standard free energy based techniques. This could be an important advantage in complex systems
5 for which a good reaction coordinate is usually very difficult to find. We explain the principles of
1 thismethod and discuss some of the promising applications related tozeolite formation.
]
h
I. INTRODUCTION ies4,5 to precise framework fragments6.
c
e The formation of crystalline zeolite particles is initi-
m
ated when this suspension is heated upto temperatures
Gaining insight in the zeolite formation has not only
- of 350 K.Light scattering experiments show that the in-
t fundamental scientific importance, but could also accel-
a tensity scatteredby the suspension increasesonly slowly
t eratemomentoustechnologicaldevelopments. Theappli- in time during the first period of the synthesis. This is
s
cationsofzeolitesareuncountablerangingfromcracking
. thenfollowedbyasharpincrease,indicatingthestarting
t
a catalysis of crude oil, gas separation, detergent builders, point of growth of what will become the final crystals.
m and sensors for pharmaceutical formulations. The spe-
The first period can be associated to a nucleation pro-
cific catalytic properties of zeolites lie in their unique
- cess, in which a particle has to be formed with a size
d open crystalline structure that incorporates cages or
beyond its critical nucleation radius. The formation of
n channelswithtypicallynanoscalediameters. Thegrowth
zeolites consist hence of several stages. First a polymer-
o
of the open structure silicon dioxide polymorphs is me-
c ization process which eventually leads to the formation
diated by so-called template molecules that can be re-
[ of sub-colloidal particles, second the nucleation process,
moved out of the zeolite pores after the crystallization
and finally the crystal growth.
1 process. Besides template molecules, solvent, Si/Al ra-
v One of the difficulties in the investigation of the ze-
tio, temperature,pH, and manyother factors playa role
5 olite formation process is that the relevant lengthscales
indeterminingwhichzeolitetopologyisformed. Aseach
3
of the zeolite formationlie just in between the accessible
structure and composition has its unique catalytic prop-
3 lengthscale of NMR and diffraction techniques7. More-
1 erties, the synthesis of new zeolite materials has been
over,itisuncleariffreezed-driedextractionsareidentical
0 an important branch of chemical research. This devel-
7 opment has progressed mainly on the basis of trial-and- to the silicate particles existing in solution. Since many
0 error and ’chemical intuition’ as a fundamental under- experiments do not allow unequivocal interpretation, it
/ is not a surprise that several crystallization mechanisms
t standing of zeolite formation is lacking. The clear so-
a have been proposed. These theories concentrate on the
lution synthesis studies of silicalite-1 zeolites initiated
m
by Schoeman et al.1 were an important step forward for structure and shape of the colloidal particles, how these
- the experimental analysis. The use of clear solutions in- particles are formed and how these particles finally con-
d
tribute in the formation of the zeolite crystal.
n stead of gels made the analysis of zeolite synthesis much
o more accessible by experimental techniques. Since then, An example is the nanoslab hypothesis that was pos-
c this model system has been subject of many studies in- tulated by some of us. It was inspired by several ex-
v: cluding x-rayandneutronscattering,infrared(IR) spec- perimental observations6,8,9. This theory assumes that
i troscopy, Nuclear magnetic resonance (NMR), and dy- at an initial stage precursor particles are formed that
X
namiclightscattering(DLS).Thesestudiesrevealedthat consist of 30 to 33 Si atoms enclosing a single template
r uponmixingtetraethylorthosilicate(TEOS),tetrapropy- molecule. These precursors stick together in a block
a
lammonium hydroxide (TPAOH) and water at a cer- shaped particle, the nanoslab, that has already the cor-
tainratioatroomtemperaturesub-colloidalparticlesare rect crystalline structure. These particles finally form
formed of several nanometers. Using the freeze drying the zeolite by a ’clicking-mechanism’ when the solution
technique2,3,theseparticleshavebeenextractedfromthe is heated up. Others have claimed that the apparent
solutionandexaminedbyvarioustechniquessuchassolid evidence of the Si-30/33 precursor particle should be
state NMR, Fourier Transform IR (FTIR), transmission attributed to other Si containing species10 or that the
electronmicroscopy(TEM),andatomicforcemicroscopy nanoshaped particle is actually an amorphous identity
(AFM). Various models for the structure of these parti- with a layer of template molecules around it4,5. Also
cles have been proposed ranging from amorphous bod- the role of the nanosizedparticles for the nucleationand
2
crystal growth has been subject of debate. According gationwere important mechanisms for diluted solutions,
to some groups, the nanoparticles add one by one to while the first mechanism was dominant in the concen-
the growing crystals11,12. Others regard the particles as trated systems. Using a implicit solvent model not in-
monomerreservoirs: monomerdissolvesintothesolution cluding template molecules, Wu and Deem analyzed the
and attaches to the growing nucleus13,14. Recent publi- free energy barriers and critical cluster sizes as function
cations14,15 state that an aggregativegrowthmechanism of pH and Si-monomer concentration at ambient condi-
of discrete nanoparticles may dominate the early stage tions using a series ofadvancedMonte Carlo(MC) tech-
of the growth process. However, after a certain size is niques28. They found that the critical clusters for the
attained, the growth mechanism seems to switch to ad- polymerization contained relatively few (≈ 30−40) Si
dition of low molecular species, probably monomers. atoms. No attempt was made to derive reaction rates
In conclusion, despite many years of abundant exper- by calculating transmissioncoefficients. Evenlargersys-
imental research, zeolite synthesis still contains many tems and timescales have been simulated using lattice
mysteries. Therefore,this fieldofresearchis aprototype models29,30 and kinetic MC (KMC)31. Relative rates for
example where computer simulations could give invalu- different crystal growth mechanisms via kink and edge
able information. However, before truly realistic simu- sites can be derived by mimicking atomic force micro-
lations of all stages in the zeolite synthesis can be per- graphsviaatomisticsimulations32,33,34. Still,evenKMC
formed, a long way has to be gone. Reason for the dif- simulations are usually restricted to growth35,36. The
ficulty is that the typical system sizes and timescales at time, before the critical nucleus of a zeolite is formed,
which the zeolite formation takes place are generally be- is still too long even for this ultrafast type of dynamical
yond the capabilities of present computer resources. For simulations.
a correct modeling of the nucleation process, the sim-
ulation box should at least be larger than the critical
nucleus. Arequirementthat isoutoflimits forquantum It is clear that the simulation methods have made sig-
mechanicalcalculationsanddemandsthedevelopmentof nificant progress in recent years. At the early stages of
accurate reactive forcefields. Si polymerization, fully quantum mechanical MD stud-
So far fully quantum mechanical calculations using ies are in our reach using Born-Oppenheimer37 or Car-
Density Functional Theory (DFT) have been applied to Parrinello37,38 simulations. Thanks to newly developed
silica polymerization clusters16,17. These studies showed potentials and coarse grained systems, simulations ap-
a stronger stability of silicate 6 rings and linear poly- proach the system sizes that are needed to describe the
mers compared to smaller rings and branched polymers. template directed zeolite synthesis in solution. Nonethe-
Ph effects were considered in [18,19] by analyzing nega- less, each stage in the zeolite synthesis involves signif-
tivelychargedsilicaclustersthatarefavorabletoneutral icant reaction barriers. This makes the chance to ob-
ones in an alkaline environment. This study revealed serve important reactive events at experimental condi-
that internal cyclization if preferred over further linear tionswithin the durationofthe simulationperiodhighly
growth18. Barriers for oligomerization were significantly unlikely. Thereactionitselfisusuallyveryfastandcould
reduced for single charged cluster compared to neutral fit perfectly within the window oftimescales that are at-
ones19. tainablebythe simulationmethod. However,thesystem
will likely spend extensively long periods within the well
Based on ab initio calculations or experimental
of the reactant state without any reactive event taking
data, several classical forcefields have been devel-
oped20,21,22,23,24. These potentials allow the study of place. It is, therefore, important to have a method that
focuses the costly simulation time on the important but
largersystemsincluding solventandtemplate molecules.
rarereactiveevents,whilelimitingthesuperfluentexplo-
However, the existing potentials are not yet very accu-
rationofthereactantwell. Inthisarticle,wereviewsuch
rately describing the breaking and making of chemical
a method, the transition interface sampling (TIS)39,40,41
bonds, which presumably requires complex many-body
method, that allows to concentrate only on those trajec-
terms and polarizable forcefields. Still, studies using
toriesthatareimportantforthechemicalprocess. More-
these approximate potentials can give valuable insights.
over, the TIS technique can calculate the frequency for
For instance, classical molecular dynamics (MD) simu-
occurrence of these successful trajectories within an in-
lations have shed some light to the role of solvent and
template molecules25,26. These simulations showedthat, finitely long straightforward simulation. Hence, TIS al-
lows the determination of the rate of the rare event.
contrarytofully formedcagesandrings,openstructures
collapse in the presence of solvent, unless it contained
strongly bonded template molecules. The early stages
of silica polymerization dynamics were studied by Rao The aim of this article is not to give a fully detailed
and Gelb27 at high temperatures & 1500 K. These al- theoretical derivation of the method. This has already
leviated temperatures were required as upto 600 K, no been published elsewhere39,40,41. The goal of this article
polymerization reaction could be observed within the is to give an educative overview of the practical algo-
nanoseconds simulation periods. They found that both rithms and their possible applications related to zeolite
themonomerincorporationandthecluster-clusteraggre- studies rather than on mathematical aspects.
3
II. TRANSITION INTERFACE SAMPLING canbecollectedthatsatisfysomepredeterminedcriteria.
For instance, one can constrain the start- and end-point
of the path in such a way that each trajectory connects
A. historic perspectives of rare event simulations
the reactant and product state. An important point is
that this sampling of successful reactive events does not
The first theories for treating rare events from a
require a RC that captures the reactive mechanism, but
microscopic perspective where pioneered by Eyring42,
only needs an order parameter that can distinguish be-
Wigner43, and Horiutu44 about 20 years before the first
tween reactant and product state. In addition to this,
MD simulation was performed45. They introduced the
the first series of TPS papers51,52,53,54 also provided a
concept of Transition State (TS) and the so-called TS
routetocalculatereactionrates. However,thisapproach
Theory(TST)approximation. LaterKeckdemonstrated
hasseldombeenusedduetoitshighcomputationalcost.
howthe TST approximationcanbe made exactbya dy-
Moreover, within the context of the reaction rate calcu-
namical correction, the transmission coefficient46. The
lation, it is not so obvious to state that the TPS order
actualapplicationofthese theoriesfor molecularsimula-
parameter is actually very different from a RC. In this
tion was directed by the works of Bennett47 and Chan-
approach,the end-point of the path is forced to progress
dler48, which have made this a standard approach in
insuccessivestepsfromreactanttoproductstate. Hence,
molecular simulation. A crucial point in this reactive
the TPS order parameter needs to describe the interme-
flux (RF) method is the definition of a suitable reaction
diatestatesaswelljustasaRCinthestandardmethods.
coordinate (RC). As a first step, the free energy needs
Luckily, the algorithmic procedure to calculate reac-
to be determined along this RC using importance sam-
tion rates using the same path sampling framework im-
pling techniques such as Umbrella Sampling (US)49 or
proved considerably when the transition interface sam-
Thermodynamic Integration (TI)50. This result alone is
pling (TIS) technique was devised39. TIS uses a flexible
sufficient to obtain the TST approximation of the rate,
pathlength which reduces the number of required MD
whichisanupperlimitfortheactualrate. Inthesecond
steps significantly. Moreover,the TIS method also elim-
step, the correction to this approximation can be calcu-
inates the need of so-called MC shifting moves that re-
lated by releasingdynamical trajectoriesfrom the top of
quired a considerable percentage of the simulation time
the free energy barrier. Only when both steps are com-
in the TPS scheme. In addition, one can show that the
pleted, the exact reaction rate can be calculated. Both
newmathematicalformulationofthereactionrateisless
the free energy barrier and the transmission coefficient
sensitive to recrossing events which guarantees a faster
depend on which RC is taken, but the final result that
convergence.
combinesthetwooutcomesisindependentofthischoice.
While TPS imposes conditions to the start- and end-
The RF method as provedits value for many systems,
point of the path to be within certain intervals of the
but also has its drawbacks. Although its result is inde-
RC,TISimposesaninterfacecrossingcondition. Except
pendent of the chosen RC, its efficiency does and sen-
for the technique proposed in Ref. [41], TIS needs a RC
sitively determines its success or failure. A non-suitable
just like the original TPS scheme. The RC is required
choiceofRCcanresultinhysteresiseffectsinthefreeen-
to define a set of interfaces between the stable reactant
ergycalculation,whichfrustratesanaccurateestimation
and product states. However, unlike the standard RF
of the barrier. Besides, even if an accurate value for the
methods, the TIS efficiency is relatively insensitive to
free energy barrier can be obtained, the corresponding
the choice of RC as was first proven in [55]. This point
transmissioncoefficientwillbeverysmallanditsevalua-
is a strong advantage in complex systems where a ’good
tionwillrequireanextremelylargenumberofpathways.
RC’ can be extremely difficult to find.
In practice, it has been experienced that finding a good
RC can be extremely difficult in high dimensional com-
plex systems. Notable examples are chemical reactions
insolution,wherethe reactionmechanismoftendepends B. the TIS algorithm
onhighly non-trivialsolventrearrangements. Also,com-
putersimulationsofnucleationprocessesuseverycompli- The TIS algorithm works as follows. First step is to
cated order parameters to distinguish between particles define a RC and a set of related values λ0,λ1,...,λn
belonging to the liquid and solid phase. This makes it with λi < λi+1. The subsets of phase- or configuration
unfeasible to construct a single RC that accurately de- points for which the RC is exactly equal to λi basically
scribes the exact place of cross-over transitions. As re- define multidimensional surfacesor interfaces which give
sult, hysteresis effects and low transmission coefficients thenametothismethod. Thesevalues/interfacesshould
are almost unavoidable. obey the following requirements: if the RC is lower than
TheproblemoffindingsuitableRCs,hasurgedthede- λ0 =λA,thesystemshouldbe inthe reactantstateA; if
velopmentofalternativemethods. In1998,Dellagoetal. the RC is higher than λn = λB the system should be in
cameupwithsuchanalternativemethodthattheycalled theproductstateB;nandthepositionsfortheinterfaces
transitionpathsampling(TPS)51,52,53,54. Thisapproach in between, λi with 1 = 1,2,...,n−1, should be set to
can be described as a MC sampling of MD pathways. optimize the efficiency. Further, the surface λA should
Using a detailed balance technique, a set of trajectories be set in such a way that whenever a MD simulation is
4
a) λ4=λB
Rejecttrialmove.Keeptheoldpath.Countitagain λ3
andrestartprocedureusingthesamepath.
b)
λ2=λi
λ1
o()KIN max λ0=λA
≫ L
hasnocrossingwithλ.i wrongbackwardending cn())KIN trialL> IorIoIm[aa.IlfIdsin.IPnad.iAlp{nilocCa1dkmpdt,odahe)rnrcxas]tt[h.mnLipincadO[ialuβnkeolteegs(lmh.KewGteslh[IyriLaNdwetuihapk(irsos’eetsa)sedii−imatpRnondreeKoocstdjtbl)ehIiicNacseiotnteb(rrmni[ioblab)diof)ts)u]m)yt]}ti].hneeednct)a]. d)
y e) f)
max>L trajector atIhVne.umnTiaafxokriemmaudmraisntadrlilobomwuteindounpma[0bthe:rl1eαn].gfDtrhoemfine g)
trialL Trial Lmax=int[L(o)/α].Continue.
V.IntegratetheNewtonianequations
ofmotionbackwardintimeuntil
reachingλ0[seeg)],λi+1[seee)]or
untilthepathlengthexceedsLmax[f)].
Continueincaseg).Otherwisereject.
h) i) j) k)
VI.Integratetheequationsofmotion
forwardintimeuntilreachingλ0
[h)andj)],λi+1[k)],orwhenthe
wholepathexceedsLmax[i)].
VII.Rejectincasei)orwhenthe
pathhasnocrossingwithλi[caseh)].
Otherwise,acceptthenewpath.
Replacetheoldpathbythenew
one[Sayk)].Updatestatistics.
l)
Acceptthenewpath.
RepeatprocedurefromstepIusingthisnewpath.
FIG.1: Illustration oftheshootingalgorithm inTIS.Thepanelsa)-l)depicttrajectories/trajectory-segments onafreeenergy
surface. The dashed horizontal lines are the TIS interfaces (n = 4 and i = 2 in this case). The algorithm requires an initial
path a) to start the loop. The length of this particular path L(o) is eighteen timeslices (endpoints are not included). At step
I, a random point is picked from this old path and some small randomized changes are applied to the velocities of all the
particles(II), followed byaMetropolis acceptance/rejection step(III). Inc)thenewvelocities haveresultedin amuchlarger
kinetic energy KIN(n) and therefore this trial move is most likely rejected. Step IV is requir ed to maintain detailed balance
between pathways of different lengths. For example, if the random numbergenerator assigns α=.59 then Lmax =30 and we
canreject whenthepathisunfinished,butalready contains31timeslices asinpanelf)andi). AtV,theequationsofmotion
areintegrated backwardsintimebyaMDalgorithm usingtheshootingpointwith reversevelocitiesasstartingpoint. AtVI)
the equations of motion are integrated forward in time starting from the same shooting point (without reversed velocities).
After a rejection the old path is kept and counted again. If accepted, the new path will automatically start at λA and cross
λi. The path can end at eitherλA as in j) or at λi+1 asin k). The fraction of sampled pathwaysthat end at λi+1 determines
PA(λi+1|λi)
releasedfromwithinthereactantwell,thissurfaceshould Here, fA is the escape flux through the first interface.
be frequently crossed. The TIS rate expression can then In a long MD simulation, this simply corresponds to the
be formulated as number of detected crossings through the surface λA di-
videdbythe totalsimulationtime (here,weassumethat
kAB =fAPA(λB|λA) (1) we will not observe a spontaneous transition to state B
5
during the simulation. For a more formal definition see loop,one needs to havea single paththat obeys the cor-
[39]). The other term PA(λB|λA) is the overall crossing rectrequirements. Thiscanalreadybequitedifficultand
probability. This is the probability that whenever the several techniques to get such a first initial path have
system crosses λA, it will cross λB before it crosses λA been suggested57. However, in TIS these initial paths
again. As λB is a surface at the other side of the bar- are generated automatically when the different types of
rier, this probability will be very small and can not be simulations are consecutively performed (See Fig. 2).
calculateddirectly. Thisprobabilitycan,however,bede- First the MD simulation is performed to calculate fA.
termined by a series of path sampling simulations using Then,aseriesofpath-samplingsimulationsfollowstocal-
the following factorization: culate PA(λ1|λ0),PA(λ2|λ1),...,PA(λn|λn−1). When
these simulationsareperformedinthis order,eachpath-
n−1 samplingsimulationcanobtainthenecessaryinitialpath
PA(λB|λA)=PA(λn|λ0)= YPA(λi+1|λi). (2) from the previous simulation (See Fig. 2).
i=0 It is important to note that the final result, the reac-
ThetermsPA(λi+1|λi)arehistorydependentconditional tion rate kAB, does not sensitively depend on the posi-
crossing probabilities which are much higher and can be tions of the outer interfaces λA and λB as long as they
are reasonable. The number of interfaces n and their
computed. In words, PA(λi+1|λi) is the probability that
positions only influence the efficiency of the method. It
λi+1 will be crossed before λA under the twofold condi-
was found that the total efficiency is optimized when for
tionthatthe systemis atthe pointto crossthe interface
each path-simulation one out of five trajectories reaches
λi in one timestep while λA was more recently crossed the next interface41,55. Hence, using some initial trial
than λi in the past. It is due to this history dependence
simulations,one canadjust the number ofinterfacesand
thatEq.(2)isexactandshouldnotbe misinterpretedas
their position to satisfy this condition. The easiest way
aMarkovianapproximation. PA(λi+1|λi)isalsoequalto
toachievethisistouseaslightvariationofthealgorithm
thenumberofallpossiblepathsthatstartatλA andend
that is shown in Fig. 1. Instead of stopping the integra-
soatntaeλrtic+rao1tsdsλiinAvgi,dewendidthbayλttie.hietHhneeunrmcλeb,i+etr1hioosfrtaeλlrlAmp,oacsansndibblheeavpceaatlachtusllatehtaeasdtt tki)o,nonwehceannthcoenttrianujeecttohreytcrarojesscetosrλyi+u1ntailsiitnrepaacnheelseλ)Aanodr
ifwecangeneratetheappropriatetrajectorieswiththeir λB. Thisalgorithmonlyrequiresknowingthepositionof
correct statistical weight. It is however not so obvious λA, λB, and λi. By examining the progress of the paths
to generate these pathways especially when λi is in the alongtheRCbeyondλi,onecandefinethenextinterface
reaction barrier region. This difficulty can be overcome λi+1 exactly at the point where 80 % of the paths have
by a MC algorithm that employs a variation of the TPS returned to λA.
shooting move. The algorithm is explained in Fig. 1.
The algorithm requires to have one path fulfilling the
correct condition. That is starting at λA and crossing C. analysis of the reaction mechanism
λi at least once before ending at either λi+1 or λA. A
crucial point is step II. After picking a random timeslice Whenthecompleteseriesofsimulationsisfinished,the
(a point that constitutes all the particle positions and reaction rate follows simply from Eqs. (1,2). In addition
momentaatacertaintimestepalongthepath),oneadds to this, the ensemble ofpathwayscanbe analyzedwhich
randomvaluestoallthemomenta. Inpractice,theseran- can yield valuable information about the reaction mech-
dom values are taken from a Gaussian distribution with anism. In this respect, the TIS path-ensembles might
a certain width σ, that should be adapted to obtain the actually prove to be more useful than the ones obtained
optimumefficiency. Ifσ issmall,therandommomentum by the original TPS method. As each ensemble contains
changes will be small as well and the new path will lie thecorrectratioofpathsprogressinguptoacertainlevel,
closelytotheoldone(ifweassumedeterministicdynam- but then either return or make a little step further, one
ics). This small deviation results in a significant chance can try to understand the characteristic differences be-
that the trial path will satisfy the requiredconditions as tween the ’successful’ and unsuccessful’ pathways. In
well which yields a good acceptance rate. However, a contrast, the TPS method aims to generate successful
too small value of σ will result in too strong correlations trajectoriesonly. Oneofthe propertiesthatcanimprove
between the accepted moves (In the extreme case when understandingofmechanismsistheoverallcrossingprob-
σ =0, one regenerates exclusively the same path). Usu- ability function. This represents a sort of path survival
ally, one tries severalvalues for σ in a seriesofshorttest probability along the RC. This function equals 1 at λ
0
simulations. It is generally assumed that the value that and PA(λB|λA) at λB. In transition, this function is
yields anacceptance rate of50 % is close to anoptimum monotonically decreasing and terminates in a horizon-
value for σ. If one wants to simulate at constant energy tal plateau when the barrier ridge is crossed completely.
insteadofconstanttemperature,stepIIIcanbereplaced Thisfunctioncouldbeconsideredasadynamicalequiva-
by a proper velocity rescaling procedure that, if needed, lentofthe free energyprofilealongthe RC.InFig.3 the
can also preserve linear and angular momentum56. overallcrossing probability function is depicted together
Another important point is that, in order to enter the with the free energy profile obtained from a nucleation
6
process of Lennard-Jones particles.
1 30
0.1 βF 25
bility 0.01 20
a 15
ob 0.001
g pr 1e-04 10Fβ
n
si 5
s
o 1e-05
cr PPTIS 0
1e-06 -5
TIS
1e-07 -10
0 10 0 20 0 30 0 40 0 50 0 60 0 70 0 80 0
largest solid cluster
FIG. 3: The crossing probability function obtained by TIS
and partial path TIS (PPTIS)59 for nucleation process of LJ
particles. The simulation data are obtained from [58]. The
free energy profile, that was calculated simultaneously using
the technique of [72], is also shown. The simulated system
contained10648particlesintotal. TheRCwasdefinedasthe
largest solid clusterinthesystem. Moredetailscanbefound
in [58]. The final crossing probability was PA(λB|λA) =
(7.8±5.5)·10−7 (TIS) and (14.1±8.7)·10−7 (PPTIS). The
crossing probability shows a plateau at 410.5 indicating that
the reaction barrier ridge is crossed. The free energy barrier
hasamaximumat243. Afterthispointabout50%(PPTIS)
till 75 % (TIS) of the paths still fail to reach the reactant
state. Although, the PPTIS and TIS results are within each
others error-bars, it is likely that PPTIS overestimates the
crossing probability dueto theimprecise choice of RC55.
The results were obtained from from Ref. [58]. In-
terestingly, after the maximum of free energy barrier is
crossed(clustersize243),themajorityofthetrajectories
(∼75%) still fail to reach the reactant state. This effect
mightbepartlyduetodiffusivemotion,butismostlikely
an effect of an improperly chosen RC. This shows that
the projection on a single RC using static free energy
calculationscanbe misleading. Neither the heightofthe
FIG. 2: This pictures illustrates some typical trajectories on
barrier nor the position of the TS dividing surface have
a free energy surface for a series of TIS simulations. The
to reflect the actual height and position of the reaction
glassy plates represent the TIS interfaces. The first type of
barrier. Indeed, Moroni et al. found that a ’good RC’
simulation(top),isastraightforwardMDsimulationwhichis
shouldat leastincorporateone moreimportant quantity
required to calculate the flux fA through the first interface.
The nextstep is apath-sampling simulation which generates which is the crystallinity of the cluster. Small clusters
pathways that start at λA and end at either λA or λ1. This < 243 with a high crystallinity were found to grow fur-
simulation yields the result of PA(λ1|λ0) and is illustrated ther easily, while large clusters with less structure were
in the middle panel. The initial path to start this simula- unstable and broke-up into smaller pieces.
tion can be obtained by taking a trajectory segment from
theinitial MD simulation. Thebottom panelshows thenext
path-sampling simulation to calculate PA(λ2|λ1). It gener-
D. recent and future developments and
ates pathways that start at λA and cross λ1 at least once.
applications
Also here, the initial path can be obtained from the previ-
ous simulation by taking one of the paths that successfully
reached λ1. Besides anefficientalgorithmfor the calculationofre-
action rates, the TIS method has provided a new math-
ematical framework to describe rare events. Recent new
simulation techniques have exploited this TIS theory.
7
Forexample,fordiffusivebarriercrossings,wheretransi- achieved by KMC, the expectation time to form a criti-
tion paths become very long, the partial path TIS (PP- cal nucleus starting from a disorder solution is generally
TIS) method was devised59. Here, using the assump- still out of reach. Therefore, most KMC studies have
tion of memory loss, much shorter paths are generated concentrated on growth rather than nucleation. Hence,
afterwhichthe overallcrossingprobabilitycanbe recon- thestudyofzeolitenucleationmightbenefitsignificantly
structed by a recursive formulation. The forward flux using combined KMC and path sampling techniques.
samplingtechnique60,61 isbasicallythe sameasTIS,but
the way to generate pathways is different. In this ap-
proach, the endpoints of all the pathways successfully
reaching the next interface are stored and starting from III. SUMMARY
eachpointasetofnewpathwaysisgeneratedinthenext
simulation. The main advantageis that the FFS scheme We have reviewed the TIS method, its variations and
supplies a route to handle non-equilibrium systems. A its possible applications for the theoretical study of zeo-
disadvantage is that this only works for stochastic dy- lite synthesis. The TIS method is a an elegant approach
namics and will always yield much stronger correlations circumventingthe timescale problemnot by speeding up
between the generated pathways and the different path the dynamics of the system itself, but by concentrating
ensembles even for the pure Brownian dynamics case. on the short time trajectories which are of interest with-
Another drawback of PPTIS and FFS methods is that out using any approximation. TIS allows to overcome
they do not possess the same RC insensitivity as TIS55. reaction barriers by a sequence of simulation series. It
If necessary, even a fully RC free approach is possible is important to realize that the barrier crossing event is
as was suggested in [41] using the TIS pathlength as a not enhanced due to some artifical force but only due
transition parameter. A nice feature of this approach to the MC acception/rejection steps that include the in-
is that it does not require to specify a specific product terface crossing condition. Hence, each trajectory in the
state. Combinations with Configurational Bias MC41,61 TIS path ensembles satisfy the correct dynamics on the
andpathswappingtechniques41,62mayalsoyieldpromis- true potential energy surface. This makes the method
ing advances for the computational efficiency. fundamentally different from, for instance, the metady-
The TIS and its variations have been applied to namics70 approach. The TIS method makes use of the
various systems ranging from simple test-systems39,59, fact that the time needed to actually cross the barrier,
nucleation58,63,64, protein folding65,66, biochemical net- the transitiontime,is muchshorterthanrelaxationtime
works60, driven polymer translocation through pores61, k−1,whichisthetimewhereinonecanexpectareactive
AB
micelle formation67, ab initio simulation of chemical re- event from an arbitrary point within the reactant well.
actions40 and DNA denaturation62. TIS can be combined with any type of dynamics such
The TIS technique can open many possible avenues asabinitio MD,Langevin,pureBrownianmotion,classi-
in the field of zeolite formation simulations at several cal MD and KMC. A requirement for application of this
stages of the process. For instance, the first elementary method is that the simulation of short trajectories can
step to Si polymerization is the condensation reaction occursufficiently fast. This limitsthe sizeofthe systems
2Si(OH) →Si O(OH) +H O.Thishasbeenstudiedby whichcanbe studied, rangingfromseveralmolecules for
4 2 6 2
abinitio staticanalysisincludingimplicitsolvent68. This abinitio dynamicstoseveralthousandmoleculesforMD,
study has revealed two possible reaction mechanisms. and even larger assemblies for KMC simulations.
Such a system is small enough to be treated by ab initio Still, substantial work has to be done before fully re-
MD37,38 including explicit solvent molecules. The TIS alistic modeling of zeolite synthesis is our reach. An
method could hence give valuable insight which reaction important requirement is the development of more ac-
mechanism dominates when dynamics and explicit sol- curate reactive force fields that can describe chemical
vent is taken into account. Classical reactive forcefields events within the environment of solvent and template
andrare eventmethods, such asTIS and PPTIS,should molecules. Recently,amoresystematicapproachforthis
makeitpossibleto simulatethe dynamics ofSipolymer- development was suggested71. Even though, the lattice
ization at much lower temperatures than hitherto was and KMC models are making substantial steps forward,
possible69. This allows to study this process at condi- inclusion of solvent effects in a lattice-type models has
tions that aremuch closer to the experimental situation. proven to be a difficult problem that has not yet been
Moreover, by a right construction of the interfaces, the solved. To conclude, the simulation methods have made
TIS method allowsto focus onreactionmechanisms and prodigious advancements in recent years and might ul-
rates of some very specific polymerization reactions, for timately give answers to important questions regarding
instance the formation of the Si-30/33 precursor8. zeolitesynthesis,thatcannotbeunambiguouslyaccessed
WiththedevelopmentoflatticeandKMCmodels,the by experimental techniques. The TIS methods can help
study of the next stagesofnucleationandzeolite growth in obtaining dynamical information for the crucial but
alsocomeintoreach. ThecombinationofKMCandpath rare reaction steps in the zeolite process. In the near fu-
sampling is a promising yet unexplored territory. De- ture, we are going to explore the application of TIS for
spite the enormous long simulation periods that can be the study of zeolite genesis.
8
Acknowledgments Government through the IAP-PAI network. T.C. ac-
knowledges the Institute for the Promotion of Innova-
This work was sponsored by the Flemish Government tion through Science and Technology in Flanders (IWT-
via a concerted research action (GOA), and the Belgian Vlaanderen) for a Ph.D. scholarship.
∗ E-mail: [email protected] 22 P.Vashista,R.Kalia,andJ.Rino,Phys.Rev.B,1990,41,
1 B. J. Schoeman, , J. Sterte, and J. E., Zeolites, 1994, 14, 12197.
110–116. 23 B. W.H.vanBeest, G.J. Kramer, andR.A.van Santen,
2 J.Hedlund,B.J.Schoeman,andJ.SterteIned.H.Chon, Phys. Rev. Lett., 1990, 64, 1955.
S.-K. Ihm, and Y. S. Uh, Progress in zeolites and micro- 24 B. P. Feuston and S.H.Garofalini, J. Phys. Chem.,1990,
porous materials, pp. 2203–2210, Amsterdam, 1996. Else- 94, 5351–5356.
vier. 25 C.R.A.Catlow, D.S.Coombes, D.W.Lewis, J. Pereira,
3 R. Ravishankar, C. E. A. Kirschhock, B. J. Schoeman, and G. Carlos, Chem. Mater, 1998, 10, 3249–3269.
P. Vanoppen, P. J. Grobet, S. Storck, W. F. Maier, J. A. 26 D.W.Lewis, C.R.A.Catlow, andJ. M.Thomson, Fara-
Martens, F. C. D. Schrijver, and P. A. Jacobs, J. Phys. day Discuss, 1997, 106, 451–471.
Chem. B, 1998, 102, 2633–2639. 27 N. Z. Rao and L. D. Gelb, J. Phys. Chem. B, 2004, 108,
4 D.D.Kragten,J.M.Fedeyko,K.R.Sawant,J.D.Rimer, 12418–12428.
D. G. Vlachos, R. F. Lobo, and M. Tsapatsis, J. Phys. 28 M. G. Wu and M. W. Deem, J. Chem. Phys., 2002, 116,
Chem. B, 2003, 107, 10006. 2125.
5 J.D.Rimer,J.M.Fedeyko,D.G.Vlachos,andR.L.Lobo, 29 M. Jorge, S. M. Auerbach, and P. A. Monson, J. Am.
Chem. Eur. J., 2006, 12, 2926–2934. Chem. Soc., 2005, 127, 2005.
6 R.Ravishankar,C.E.A.Kirschhock,P.P.Knops-Gerrits, 30 F.R.SipersteinandK.E.Gubbins,Stud.Surf.Sci.Catal.,
E.J.P.Feijen,P.J.Grobet,P.Vanoppen,F.C.D.Schri- 2002, 144, 647–654.
jver, G. Miehe, H. Fuess, B. J. Schoeman, P. A. Jacobs, 31 D. T. Gillespie, J. Comput. Phys., 1978, 28, 395.
and J. A. Martens, J. Phys. Chem. B, 1999, 103, 4960– 32 M. W. Anderson, J. R. Agger, N. Hanif, O. Terasaki, and
4964. T. Ohsuna, Sol. State. Sci., 2001, 3, 809–819.
7 S. M. Auerbach, M. H. Ford, and P. A. Monson, Curr. 33 M. W. Anderson, J. R. Agger, N. Hanif, and O. Terasaki,
Opin. Colloid Interface Sci, 2005, 10, 220–225. Microporous Mesoporous Mater., 2001, 48, 1–9.
8 C. E. A. Kirschhock, R. Ravishankar, F. Verspeurt, P. J. 34 J. R. Agger, N. Hanif, and M. W. Anderson, Angew.
Grobet, P. A. Jacobs, and J. A. Martens, J. Phys. Chem. Chem., 2001, 40, 4065–4067.
B, 1999, 103, 4965–4971. 35 S. Piana and J. D. Gale, J. Am. Chem. Soc., 2005, 127,
9 C. E. A. Kirschhock, R. Ravishankar, L. van Looveren, 1975–1982.
P.A.Jacobs, andJ. A.Martens, J. Phys. Chem. B,1999, 36 S. Piana, M. Reyhani, and J. D. Gale, Nature, 2005, 438,
103, 4972–4978. 70–73.
10 C.T.G.Knight,J.Wang,andS.D.Kinrade,Phys.Chem. 37 D. Marx and J. Hutter In ed. J. Grotendorst, Modern
Chem. Phys., 2006, 8, 3099–3103. Methods and Algorithms of Quantum Chemistry, pp.301–
11 V. Nikolakis, E. Kokkoli, M. Tirrell, M. Tsapatsis, and 449, FZ Ju¨lich, 2000. NIC.
D. Vlachos, Chem. Mater., 2000, 12, 845–853. 38 R. Carr and M. parrinello, Phys. Rev. Lett., 1985, 55,
12 W. H. Dokter, H.F. Vangarderen, T. P. M. Beelen, R.A. 2471–2474.
vanSanten,andW.Bras,Angew.Chem.,1995,34,73–75. 39 T. S. van Erp, D. Moroni, and P. G. Bolhuis, J. Chem.
13 B.Schoeman,Micropor. Mesopor. Mater.,1998, 22,9–22. Phys., 2003, 118, 7762–7774.
14 C.S.CundyandP.A.Cox,Micropor.Mesopor.,2005,82, 40 T. S. van Erp Solvent effects on chemistry with alcohols
1. PhD thesis, Universiteit van Amsterdam, 2003.
15 T.M.Davis,T.O.Drews,H.Ramanan,C.He,J.S.Dong, 41 T. S. van Erp and P. G. Bolhuis, J. Comput. Phys., 2005,
H. Schnablegger, M. A. Katsoulakis, E. Kokkoli, A. V. 205, 157–181.
McCormick,R.L.Penn,andM.Tsapatsis,NatureMater., 42 H. Eyring, J. Chem. Phys., 1935, 3, 107.
2006, 5, 400–408. 43 E. Wigner, Trans. Faraday Soc., 1938, 34, 29.
16 J. C. G. Pereira, C. R. A. Catlow, and G. D. Price, J. 44 J. Horiuti, Bull. Chem. Soc. Jpn, 1938, 13, 210.
Phys. Chem. A, 1999, 103, 3252–3267. 45 B. J. Alder and T. E. Wainwright, J. Chem. Phys., 1957,
17 J. C. G. Pereira, C. R. A. Catlow, and G. D. Price, J. 27, 1208.
Phys. Chem. A, 1999, 103, 3268–3284. 46 J. C. Keck,Adv. Chem. Phys., 1967, 13, 85.
18 M. J. Mora-Fonz, C. R. A. Catlow, and D. W. Lewis, 47 C. H. Bennett In ed. R. Christofferson, Algorithms for
Angew. Chem.,2005, 44, 3082–3086. Chemical Computations, ACS Symposium Series No. 46,
19 T. T. Trinh, A. P. J. Jansen, and R. A. van Santen, J. Washington, D.C., 1977. American Chemical Society.
Phys. Chem. B, 2007, 110, 23099–23106. 48 D. Chandler, J. Chem. Phys., 1978, 68, 2959.
20 C.R.A.Catlow, C.M.Freeman,M.S.Islam,R.A.Jack- 49 G. M. Torrie and J. P. Valleau, Chem. Phys. Lett., 1974,
son,M.Leslie,andS.M.Tomlinson,Philos.Mag.A,1988, 28, 578.
58, 123–141. 50 E. A. Carter, G. Ciccotti, J. T. Hynes, and R. Kapral,
21 S.Tsuneyuki,M.Tsukada,H.Aoki,andY.Matsui, Phys. Chem. Phys. Lett., 1989, 156, 472.
Rev. Lett., 1988, 61, 869. 51 C. Dellago, P. G. Bolhuis, F. S. Csajka, and D. Chandler,
9
J. Chem. Phys., 1998, 108, 1964. 63 C. Valeriani, E. Sanz, and D. Frenkel, J. Chem. Phys.,
52 C. Dellago, P. G. Bolhuis, and D. Chandler, J. Chem. 2005, 122, 194501.
Phys., 1998, 108, 9236. 64 A. J. Page and R. P. Sear, Phys. Rev. Lett., 2006, 97,
53 P. G. Bolhuis, C. Dellago, and D. Chandler, Faraday Dis- 065701.
cuss., 1998, 110, 421. 65 P. G. Bolhuis, Proc. Nat. Acad. Sci. USA, 2003, 100,
54 C. Dellago, P. G. Bolhuis, and D. Chandler, J. Chem. 12129–12134.
Phys., 1999, 110, 6617. 66 P. G. Bolhuis, Biophys. J.,2005, 88, 50–61.
55 T. S.van Erp, J. Chem. Phys., 2006, 125, 174106. 67 R.PoolCoarse grained soap A molecular simulation study
56 P. L. Geissler, C. Dellago, and D. Chandler, Phys. Chem. on Micelle formation in dilute surfactant solutions PhD
Chem. Phys., 1999, 1, 1317. thesis, Universiteit van Amsterdam, 2006.
57 P.G.Bolhuis,D.Chandler,C.Dellago,andP.L.Geissler, 68 J.C.G.Pereira,C.R.A.Catlow, andG.D.Price,Chem.
Ann. Rev. Phys. Chem.,2002, 53, 291–318. Commun.,1999, 13, 1387–1388.
58 D. Moroni, P. R.ten Wolde PR, and P. G. Bolhuis, Phys. 69 N. Z. Rao and L. D. Gelb, J. Phys. Chem. B, 2004, 108,
Rev. Lett., 2005, 94, 235703. 12418–12428.
59 D. Moroni, P. G. Bolhuis, and T. S. van Erp, J. Chem. 70 A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. USA,
Phys., 2004, 120, 4055–4065. 2002, 99, 12562.
60 R.J.Allen,P.B.Warren,andP.R.tenWolde,Phys.Rev. 71 t. verstraelen, d.van neck,p. w. ayers, v.van speybroeck,
Lett., 2005, 94, 018104. and m. waroquier, submitted.
61 R. J. Allen, D. Frenkel, and P. R. ten Wolde, J. Chem. 72 D. Moroni, T. S. van Erp, and P. G. Bolhuis, Phys. Rev.
Phys., 2006, 124, 024102. E, 2005, 71, 056709.
62 T. S.van Erp, In Progress.
