Table Of ContentDiffusion in Quasi-One Dimensional Channels: A Small System n,p,T, Transition
State Theory for Hopping Times.
Sheida Ahmadi1 and Richard K. Bowles1,∗
1Department of Chemistry, University of Saskatchewan, Saskatoon, Saskatchewan S7N 5C9, Canada
Particles confined to a single file, in a narrow quasi-one dimensional channel, exhibit a dynamic
crossover from single file diffusion to Fickian diffusion as the channel radius increases and the
particles can begin to pass each other. The long time diffusion coefficient for a system in the
crossover regime can be described in terms of ahoppingtime, which measures thetime it takes for
a particle to escape thecage formed byits neighbours. In this paper, we develop a transition state
7
1 theoryapproach tothecalculation of thehoppingtime,using thesmall system isobaric–isothermal
0 ensembletorigorously account forthevolumefluctuationsassociated withthesize ofthecage. We
2 alsodescribeaMonteCarlosimulation schemethatcanbeusedtocalculatethefreeenergybarrier
for particle hopping. The theory and simulation method correctly predict the hopping times for
n
a two dimensional confined ideal gas system and a system of confined hard discs over a range of
a
channelradii,butthemethodbreaksdownforwidechannelsintheharddiscscase,underestimating
J
the height of the hopping barrier due to the neglect of interactions between the small system and
4
its surroundings.
2
]
t I. INTRODUCTION or metal organicframeworks[16], and charge-carrierdif-
f
fusion in one dimensional channels [17]. At the molecu-
o
s ParticlemotioninbulksystemsisdescribedbyFickian lar level, SFD has been observed experimentally in ad-
t. (normal) diffusion, where the mean square displacement sorbate molecules confined in zeolites using pulsed field
a (MSD)ofalabeledparticleincreaseslinearlywithtimein gradient nuclear magnetic resonance [18–20]. However,
m
thelongtimelimit, andisgivenbytheEinsteinrelation, the phenomenon has also been observed in colloidal sys-
- tems [21, 22], where it is being used in nano- and micro-
d
n h(xt−x0)2i=h∆x2ti=2Dxt, (1) fluidic devices [23, 24], as well as to control flow rates in
o the delivery of drugs [25].
c where Dx is the diffusion coefficient along an arbitrary Fluids in quasi-one dimensional channels also exhibit
[ direction x, x is the position of the particle at time, t, an interesting dynamical crossover regime from SFD to
t
and x is the initial particle position. However, the fun- normal diffusion as the channel radius increases beyond
1 0
v damental nature of diffusion can change dramatically if a passing threshold that allows the particles to hop past
8 the movementoftheparticlesisrestrictedtoaquasi-one eachother [26–30]. In this crossoverregime,the MSD in
1 dimensionalnarrowchannelinwhichtheparticlescannot the long time limit is described by Eq. 1 because the
9 pass each other [1, 2]. The geometric restriction, when particles can pass each other. However, the particles
6
combined with the presence of independent stochastic are trapped between their neighbors for long periods of
0
forces [3, 4] or a Brownian background [5] acting on the time, during which they perform SFD, because passing
.
1 particles,results in a form ofanomalousdiffusion known is difficult. The resulting MSD at intermediate times is
0 as single file diffusion (SFD), because the relative order then described by Eq. 2. In the case of mixtures formed
7
of the particles is preserved and a particle must wait for from particles of different sizes, it is also possible to se-
1
its neighbors to move before being able to diffuse. This lect a channel radius that induces dual–mode diffusion,
:
v effect reduces the translational mobility of the particles where the large particle component performs SFD while
Xi in a waythat causes the MSD to increaseproportionally thesmallercomponentsexhibitnormaldiffusion[31,32].
to the square root of time in the long time limit and is The difference in particle mobilities can potentially be
r
a described by an Einstein–like relation, used for separation [33, 34].
Mon and Percus [26] used a simple phenomenological
h∆x2ti=2Fxt1/2, (2) model to show that the long time diffusion coefficient of
a system in the crossoverregime can be written,
where the diffusion coefficient has been replaced by a
D ∼1/τγ . (3)
mobility factor, Fx [6, 7]. x hop
Single file diffusion occurs in a variety of natural and where τ is the average time for a particle to escape
hop
engineered systems, including ion and water transport the cageformedbyits twonearestneighborsandthe ex-
through biologicalmembranes [8–11], diffusion in molec- ponent γ depends on the nature of the hopping event.
ularsievessuchaszeolites[12–14],carbonnanotubes[15] If the system is near the passing threshold for the par-
ticles, so that hopping events are rare and the particles
are trapped long enough to perform SFD, then γ is pre-
dicted to be 1/2, and this has been confirmed by simu-
∗ richard.bowles@usask.ca
lationfor quasi-onedimensionalharddiscs [35, 36], hard
2
spheres[26]andsoft-potentialsystems suchas Lennard– II. BACKGROUND AND THEORY
Jones [37], as well as for mixtures [31]. The hopping
time approach is useful because the important details A. The Small System n,p,T Ensemble
that influence the long time dynamics of the particles,
such as the particle–particle and particle–wall interac-
In the thermodynamic limit, fluctuations away from
tions, as well as the density or pressure, are captured in
themostprobablemicrostatesaresmallandthethermo-
a single, short time, phenomenologicalparameter that is
dynamics of the system can be described by the prop-
also accessible by simulation and theory.
erties of the states associated with the maximum term
in the probability distribution. The maximum term for
one ensemble can be shown to be degenerate with the
partition function of another so that all the ensembles
The process by which particles pass each other in a are equivalent and selecting a particular ensemble for
narrow quasi-one dimensional system can be treated in the calculation of a system property becomes a matter
terms of a free energy barrier crossing process because of mathematical convenience. However, as the system
the particle–particle and particle–wall interactions com- size decreases, fluctuations become more important and
binetoincreasetheenergyofinteraction,andrestrictthe in the nano-scale regime it is essential to consider the
availableconfigurationspacefortwoparticles,astheyat- types of fluctuations available to the system in selecting
tempttopass. Thissuggeststhattransitionstatetheory an appropriate ensemble. In the case of the isobaric-
(TST) could provide a useful description of the hopping isothermal,(N,P,T)ensemble,wherethe volumefluctu-
time. In the case of hard sphere systems, τhop diverges ates, it also becomes necessary to consider the details of
as the radius of the channel decreases towards the pass- how the ensemble itself is formulated [39–46]. The origi-
ing threshold from above because the infinite nature of nal formulationof the ensemble by Guggenhiem [47] can
theinteractionpotentialeventuallypreventstheparticles be expressed in integral form [48] as
from ever passing and TST correctly predicts the expo-
nent for the power law behavior [38]. Transition state ∆(N,p,T)=(V )−1 Q(N,V,T)e−pV/kT, (4)
o
theory also provides a qualitative description of τhop for ZV
softpotentialsystems,wherethepassingthresholdisnot
where Q(N,V,T) is the canonical partition function for
as clearly defined [37]. However, these approaches are
N particles contained in a volume V at temperature T,
not quantitative. The goal of this paper is to develop
k is Boltzmann’s constant and p is the external pressure
a rigorous approach to the calculation of the hopping
being applied to the system. The prefactor contains a
time, within the framework of transition state theory,
volumescale,V ,toensurethepartitionfunctionremains
o
thatcanprovidethebasisforthequantitativeprediction
dimensionless. In the thermodynamic limit the choice of
ofτ . Toachievethis,wedevelopthetheoryandassoci-
hop V is arbitrary and it is usually taken to be kT/p, but
o
atedsimulationmethodusingthe smallsystemisobaric–
for small systems, the value of the volume scale depends
isothermalensembletoaccountforthefluctuatingnature
on the properties of the system such as the number of
of the nearest neighbour cage that traps a particle in a
particles and the details of the boundary defining the
single file fluid.
volume.
To address these concerns Corti [45] developed a rig-
orous approach to the isobaric-isothermal ensemble for
small systems that avoids the need to specify a volume
The remainder of the paper is organized as follows: scale. Here, we will provide a brief summary of Corti’s
Section II provides a brief summary of the background derivationforasmallsystem,n,p,T ensemble,highlight-
andgeneralderivationofthe smallsystemn,p,T ensem- ingthekeyelementsofthemethodandhowtheypertain
ble in orderto highlightthe key elements ofthe method. toourparticularapplication. Theanalysisbeginsbycon-
This is followed by the description of how the small sys- sidering the division of the canonical partition function
tem n,p,T ensemble can be used to develop a transi- for N particles, in a volume V and at a temperature T,
tionstate theory for hoppingtimes inconfined singlefile intoasmallsubsystemconsistingofnparticlesinasub–
fluids. Section III outlines the Monte Carlo (MC) sim- volume v that is located at a point r and surrounded
0
ulation methods used to calculate the transition state by a system of N −n particles in the remaining volume,
freeenergybarriersassociatedwithparticlehoppingand V −v (see Fig. 1(a)). Separating the subsystem from
the direct measurement of the hopping times in a large itssurroundingsisamathematicalboundarythathasno
system. Section IV describes the application of the the- mass or momentum and is assumed to be spherical. To
ory and simulations for two simple confined systems: a avoid an over counting of configurations caused by the
twodimensionalsystemofidealgasparticlesconfinedbe- fluctuation of the volume that does not involve a change
tween hard walls and a two dimensional system of hard in particle coordinates, it is necessary to tie the location
discs confined between hard walls. Section V contains oftheboundarytoashellmolecule. Thisinvolvesrequir-
our discussion and our conclusions are summarized in ing that at least one of the n particles of the subsystem
Section VI. is located in a shell, dv, at the boundary. Now when
3
the volume fluctuates, the shell molecule moves with the in the small volume interacting amongst themselves, U
σ
boundaryto ensure the systemsamples anew configura- is the interaction between the particles in the subsys-
tion. tem and the N − n particles in the surroundings and
(dr)n−1 =dr ···dr aretheparticlevolumeelements
1,2 1,n
with positions relative to the location of the shell parti-
cle. Furthermore, the integration of the shell particle,
overthe volumeinthe shell, dv,andthe othern−1par-
ticles over v+dv, is performed with the positions of the
surrounding N −n particles held fixed so the canonical
partitionfunctionforthecompleteN,V,T systemcanbe
expressed,
V
∗
Q(N,V,T)= Q(N −n,V −v,T) Q dv, (6)
n,v 0
Z0
(cid:10) (cid:11)
where,
Q∗ e−βUN−n(dr)N−n
Q∗ = V−v n,v , (7)
n,v 0 R V−ve−βUN−n(dr)N−n
(cid:10) (cid:11)
is an average over theRconfigurations of the surrounding
N −n particles interacting amongst themselves with a
potential energy UN−n.
Now,the probabilityoffindingthe systemwith npar-
ticles in a volume v + dv and the remaining particles
outside can be written,
∗
Q(N −n,V −v,T) Q dv
P (v)dv = n,v 0
n
Q(N,V,T) (8)
(cid:10) (cid:11)
=eβnµe−βW(v) Q∗ dv,
n,v 0
where µ is the chemical potent(cid:10)ial of(cid:11)the surroundings,
W(v) is the work required to form an empty cavity in
the system and the following relation has been used:
Q(N −n,V −v,T)
=eβnµe−βW(v). (9)
Q(N,V,T)
The small system n,p,T partition function can then
be written,
∆= Q∗ e−βW(v)dv, (10)
n,v 0
Z
FIG. 1. (a) A configuration of a bulk canonical system de- (cid:10) (cid:11)
by using the normalization condition, P (v)dv = 1
composedintoasmalln,v,T-subsystem,locatedatr0,witha v n
shellmolecule, anditsN−n,V −v,T surroundings. (b)The and noting ∆=exp(−βnµ). DependingRon the size and
shape of the cavity, W(v) may contain work in addition
two dimensional, quasi-one-dimensional small n,p,T system,
to the usual volume work, such as surface work, so that
located at r0, with a shell particle located at L/2 and the
caged particlelocated at adistancealong theaxialdirection, the pressure inside the system, pˆ, does not necessarily
xc, relative to the shell particle. (c) The transition state for equal the imposed pressure. However, if the interface
thehopping process with xc =0. is planar, or the small subsystem becomes macroscopic,
then pˆ= p and W(v) = pv. If the interactions between
the subsystem and its surroundings are also negligible,
Acanonicalpartitionfunctionforthenparticlesinthe
Eq. 10 becomes,
volume v+dv can be written as,
dv ∆= Q∗ e−βpvdv. (11)
Q∗ dv = e−β(Un+Uσ)(dr)n−1, (5) n,v
n,v (n−1)!ΛDn Z
Zv+dv
Mostnotably,Eq.11containsnovolumescalebecause
where Λ is the de Broglie wavelength, D is the dimen- the uncertainty associated with identifying new configu-
sionality, U is the potential energy of the n particles rations in phase space as the volume of the container
n
4
fluctuates is avoided by the use of the shell molecule. Figure1(b) showsthe constructionofthe small n,p,T
However, considerable care needs to be taken when im- systemusedheretostudytheTSTfreeenergybarrieras-
plementing the new small system N,p,T. For example, sociated with a particle escaping its cage in a quasi-one
the integral over the shell particle position performed to dimensional channel by hopping past one of its neigh-
obtainEq.11reliesontheunderlyingsphericalsymmetry boring particles. We consider an n = 2 particle system
of the system volume. The same approach can be taken confined to a long, uniform, two dimensional (2D) chan-
forsystemswithalternativeshapessuchascubesetc,but nel that extends longitudinally along the x-axis and has
thisisnotalwayspossiblewhenconsideringnon-isotropic a radius R in the y-axis. The cell center is located at a
p
systems, as we will do in this work, and alternative ap- pointr andthe shellparticle islocatedat+L/2sothat
0
proachestoevaluatingdegreesoffreedomassociatedwith the total length of the cell is L and its two dimensional
the shell particle may be needed. volume is v = 2R L. The channel radius remains fixed
p
Corti [46] also showed that the small system isobaric- so that fluctuations in the volume are given by,
isothermalpartitionfunction is relatedto the thermody-
namic potential, dv =2RpdL. (15)
∗
G=hU i+phvi−TS, (12) The small system isobaric–isothermal partition function
∗ given by Eq. 11 then becomes,
where h···i denotes the ensemble average and U is the
internalenergyinthepresenceofashellparticle,through
the expression ∆= Q∗n,ve−βpl2RpL2RpdL, (16)
Z
G=−kT ln∆. (13)
wherep isthelongitudinalpressureactingontheendof
l
It is importantto note that G does notnecessarilyequal thechannel[52]. Particle1isthe shellparticleplacedon
the Gibbs free energy of the system, except in the ther- the boundary at one end and the 2Rp factor in dv arises
modynamiclimit,becauseitcontainsinformationregard- from the integration of the shell particle over its y-axis
ing the surroundings, as is evident from Eq. 8 where the degree of freedom.
chemical potential of the surrounding is introduced, and Thesecondparticleinoursystemrepresentsthecaged
because the small system samples a different set of vol- particle. If we define a reaction coordinate for the hop-
umefluctuationscomparedtoathermodynamicallysized ping processas the axialseparationof the cagedparticle
system. This is the case even when Uσ =0, and there is from the shell particle, xc =x1−x2, then the transition
nointeractionbetweenthe systemandthesurroundings. state occurs at xc = 0 as shown in Fig. 1(c). The con-
Consequently, a change in G, defined by Eq. 13, repre- figurationspaceassociatedwith specific points alongthe
sents the reversible work of moving the system between reactioncoordinatecanbe measuredby applyingadelta
′ ′
two states, including contributions from the surround- function, δ(xc −xc), that is one when xc −xc = 0 and
ings. A more complete discussion of these subtleties can zerootherwise,toEq.16,toobtainareactioncoordinate
be found in references [44–46]. partition function,
B. A Transition State Theory for Hopping Times ∆x′c dxc = Q∗n,v,T e−βpl2RpL δ(x′c−xc)2Rp dL, (17)
ZL
which satisfies,
Transition state theory is an equilibrium based ap-
proach to the calculation of transition rates in activated
processes,originallydevelopedbyEyring[49]todescribe ∆= ∆xc dxc. (18)
rates in chemical reactions. While the approach ne- Zxc
glects the effect of barrier recrossing, it provides an up- ′
Theprobabilityofhavingthecagedparticleatx isthen,
per bound to the rate and more sophisticated methods, c
suchasthoseproposedbyBennet[50]andChandler[51], ′ ∆x′ dxc
reducetoTSTintheregimewherethebarrierishighrel- P(xc)dxc = c∆ , (19)
ative to the thermal energy in the system and recrossing
is rare. where P(x′) is a probability density and
c
AccordingtoTST,thetimeittakesforabarriercross-
∞
ing event like hopping to occur is inversely proportional
P(x )dx =1. (20)
to the rate and is given by, c c
Z0
τ = κ =κeβ∆G∗, (14) Finally, the Gibbs free energy barrier for particle hop-
hop P∗ ∗
ping,∆G ,isdirectlyrelatedtotheprobabilityoffinding
∗
where κ is a kinetic prefactor, P is the probability of the caged particle in the transition state by,
∗
findingthe systeminthetransitionstateand∆G isthe
∗ ∗
heightofthefreeenergybarrierfortheactivatedprocess. β∆G =−lnP (21)
5
where particleoverthediameterofthechannelandtheintegral
over dL is the degrees of freedom of the shell particle
x∗
∗ c associated with the volume of the cell. If a scale factor
P = P(x )dx , (22)
c c
x =Lx¯ is introduced, then Eq 23 can be expressed as,
Z0 i i
∗
and the range of the reaction coordinate, x = 0 to x , 1
c c∗ ∆= ×
is the transition state region. Defined in this way, ∆G (N −1)!Λ2n
represents the work required to bring the caged particle
into the transition state from anywhere in the cage and, Ln−1e−βUne−βpl2RpLdLdy dx¯n−1dyn−1,
1
as such, it will always be positive because the configura- Z Z Z Z
(24)
tion space of the transition state is a restrictedsubset of
all the possible configurations available to the system.
where the U is now a function of the scaled coordinate
n
It is important to note that we have chosen to focus
system. Notethereisnoneedtorescalethey-coordinates
ona systemofjusttwoparticlesso thatwehaveasingle
because fluctuations in the volume only involve changes
transitionstatetoconsiderwhenaparticleisattempting
in the longitudinal dimension of the cell. The probabil-
to hop. It could be argued that the cage of the central
ity of finding the system in a configuration specified by
particle is actually formed by two particles, one on ei- x¯n,yn,inavolumewiththe celllengthinthe rangeLto
ther side, giving rise to two transition states. However,
L+dL, is then,
each transition state involves two particles, and both of
theseparticlesescapetheircagesforeachhoppingevent, P(x¯n,yn;L)dL∝Ln−1e−βUne−βpl2RpLdL, (25)
which suggests we only need to consider one transition
per caged particle. Explicitly including the third parti- andthe MCacceptanceprobabilityforatrialmovefrom
clewouldalsocomplicatetheanalysisasitwouldrequire an old to a new configuration becomes,
the use of two shell molecules to define the volume. The
new old
acc(old→new)=min(1,exp{ −β[U −U ]
reaction coordinate, which measures the distance of the n n
cageparticlefromthetransitionstate,wouldalsotakeon −βp 2R [Lnew−Lold]+(n−1)ln[Lnew/Lold]}).
l p
twovaluesforeachconfigurationduetothesymmetryof
(26)
the system. Furthermore, the configurations of the par-
ticles outside the cell are formally included in the anal- While our analysis is specific to the two dimensional
ysis when the full canonicalpartition function is divided channels studied here, extending Eq. 26 to the study
into the small system to formal the small system n,p,T of three dimensional, quasi-one dimensional channels,
ensemble. This include configurations where a particle where fluctuations in the volume only involvechanges in
fromthe surroundingssitsonthe boundaryofthe cellat L, simply requires the appropriate geometric factor for
−L/2, to form the other end of the cage. the channel cross section to be included in the pressure-
volume term.
The probability, P(x )dx, of finding the caged parti-
c
III. SIMULATION METHODS cle in a volume dx at a distance x from the transition
c
state can now be calculated directly using MC simula-
A. Free Energy Barrier Calculations tion. The small system n,p,T simulations are carried
out with n = 2 particles confined to a channel with ra-
Corti [46] showed that the probability for accepting a dius, Rp, and length, L. Particle one is the shell par-
trial MC move in the small system n,p,T ensemble de- ticle that defines the sub-volume and is located at L/2
pends on the shape of the small system volume, which andthe centerofthe sub-volumeislocatedattheorigin.
also constitutes the simulation cell. If we consider the Particle two is placed within the cell between −L/2 and
caseof n particles confinedto the twodimensionalchan- L/2. Since the two models being studied are the two di-
nel described in Fig. 1(b), and ignore interactions be- mensional ideal gas and two dimensional hard discs, the
tweenthe particlesinthe cellandthosein the surround- MC moves proceed as follows: A particle is selected at
ings, then the partition function for the system can be random and moved randomly by a step δx and δy, up
written, to a maximum displacement of |∆x| = |∆y| = 0.05σ.
The move is immediately rejected if the trial displace-
1 ment causes the particle to leave the channel or over-
∆= ×
(N −1)!Λ2n lap with with another particle. If particle one was se-
lected, the trial move also corresponds to a change in
e−βUne−βpl2RpLdLdy1dxn−1dyn−1, L of 2δx and the position of particle two is rescaled to
Z Z Z Z (23) ensure it remains within the simulation cell. Eq. 26 is
then used to determine the probability of accepting the
wheretheintegralsoverthedxdy coordinatesofthen−1 move. It should also be noted that the position of the
∗
particles within the cell represent Q , the integral shell particle must remain positive so that V > 0. For
n,v,T
over dy represents the degrees of freedom of the shell each state point studied, 2×107 MC moves are used to
1
6
obtainequilibriumanddataiscollectedoverthenextthe |∆x|=|∆y| is varied during the simulation in the range
next108−109moves,dependingontheheightofthefree of 0.05−0.1σ to ensure 80% of the trial moves are ac-
energy barrier. The probability is calculatedby building cepted. The shell particle approach is used. Data are
a histogram of configurations along the reaction coordi- sampled at every 10,000 time steps to obtain the aver-
nateusingbinsizesof∆x =0.005and0.03fortheideal age volume, hVi, which gives the density ρ = N/hVi.
c
gasandharddiscsrespectively,andsamplingthe system Table I shows the range densities, ρ = N/V, obtained
every1000MCmovesto ensuresconfigurationswerenot for the different radii studied in both the ideal gas and
correlated. Simulation runs were sufficiently long to en- hard disc systems. The density range for hard disc two
sure the transition state region, x < ∆x , was sampled particle system is ρ = 0.54−0.65, which is significantly
c c
at least 100 times and the probability density was ob- higher than that obtained for the large system because
tained by dividing the probability by the bin size. the lack of interactions with the surrounding allows the
system to small volume configurations that would be in-
accessible if more particles were present. When a third
B. Hopping Time Calculations particle is included in the small system N,p,T, the den-
sity range reduces to ρ = 0.47−0.53. Nevertheless, as
noted earlier, it is necessary to calculated τ using the
The hopping time, which is the average time it takes hop
twoparticlesystembutthis resulthighlightsthe needto
for a particle to pass one of its caging nearest neighbors
ultimately account for the interactions between the cell
in the single file, is calculated using MC simulation in
and the surrounding, i.e. the inclusion of U .
the canonicalensemble. The system consists of N =100 σ
particles confined to a two dimensional channel of ra-
diusR inthey–directionandlengthLalongthex–axis, Systems Rp ρ
p
where periodic boundaries are also employed. Particles 2D-ideal gas 1.1-1.2 1.00
aremovedusingthestandardMetropolisMCscheme[53] 2D-Hard discs 1.01-1.2 0.38-0.41
as a simple approximation to Brownian motion [54]. A
MC trial move involves the random selection of a parti- TABLEI.Selectedchannelradiiandsystemdensitiesforthe
idea gas and hard sphere systems at βP =1.0.
cle, followed by the random selection of a trial step in
both the x and y–directions up to a maximum step size
of |∆x| = |∆y| = 0.05σ. The trial move is rejected if it
causeseitheraparticleorwalloverlapandtheparticleis
IV. APPLICATIONS
returnedtoits originalposition,otherwiseitis accepted.
AunitoftimeisdefinedasasingleMCcycleandconsists
of N MC trial moves. The particles are initially spaced A. Two Dimensional Ideal Gas
evenly along the tube and placed randomly across the
channel such that they do not overlap then the system To illustrate our approach, we first calculate the hop-
is equilibrated over 2×107 MC cycles. At the start of pingbarrierinatwodimensionalidealgassystem,where
data collection, the cages for the particles are identified thereareno interactionsbetweenparticles,sothatU =
σ
astheirimmediaterightandleftneighbors,andtheirini- U = 0, and the hard interaction of the particles with
n
tial hopping time is set to zero. A particle hopping time wallskeepsthemconfinedbetween−R andR . There-
p p
isthenumberofMCcyclesittakesfortheparticletoes- action coordinate partition function for this system can
capeitscage. Aftereachhoppingeventthehoppingtime be written,
for the particle is recorded, then reset to zero, and the
dx +Rp +Rp ∞
new caging neighbors are identified. Data is collected ∆ dx = c e−βpl2RpL dLdy dy
over 5×107 MC cycles and the average hopping time, xc c Λ4 Z−Rp Z−Rp Zxc 1 2
τ , is calculated overallhopping times recordedfor all (27)
hop
particles. where the length of the cell is bounded in its lower limit
The hopping times for systems with different radii bythepositionofthecagedparticle,andtheintegralsy1
must be compared at the same applied external pres- and y2 can be performed independently to give,
isnurteh.eHNo,wVe,vTer,entsheemcballecubleactaiounsseovfolτuhmopeaflruecctuarartiieodnsouint ∆ dx = 2Rpe−βpl2Rpxc dx . (28)
xc c Λ4 βp c
the N,p,T ensemble cause the particle positions to be l
rescaled,leading to possible errors. To ensure τhop is ob- The partition function for the system can then be ob-
tainedatthecorrectstatepoint,thedensityforeachsys- tained by integrating over the remaining x coordinate
c
temwithadifferentradiusiscalculatedusingtheN,p,T to yield,
ensemble where βP = 1.0, N = 1000, and 108 MC cy-
∞
2R
cles are used after the system reaches equilibrium. A ∆= e−βpl2Rpxc dx
MC trial move involves the random selection of a parti- βplΛ4 Z0 c (29)
cle, followed by the random selection of a trial step in 1
= .
both the x and y–directions. The maximum step size of β2p2Λ4
l
7
To demonstrate that we have the correct partition func-
tion, we calculate the system volume,
1 ∂ln∆ 2 6
V =2R hLi=− = , (30)
p
β ∂p βp
(cid:18) l (cid:19)n,T l
4
which is the expected equation of state for an ideal )
2
x
gas [45]. P( 2
UsingEqs.28and29inEq.19yieldstheprobabilityof n
l
findingthecagedparticleatapointx alongthereaction -
c 0 R=1.1
coordinate as, R=2.0
R=4.0
P(x )dx =βp 2R e−βpl2Rpxc dx . (31) -2
c c l p c
0 1 2 3
The probability of finding the caged particle within a x
2
region dx of the transition state, where x =0, is then,
c c
P(0)dx =βp 2R dx . (32) FIG. 2. Free energy density, −lnP(xc) as a function of xc,
c l p c
for three different pore radii with fixed pressure, βpl = 1.0,
for the 2D ideal gas system. The solid lines represent results
Figure 2 shows that the free energy density, −lnP(x ),
c from the theory (see Eq. 31). The points represent the data
increases linearly as a function of xc, with a slope of obtained from simulation (see Section III for details).
βp 2R ,whichreflectsthefactthatincreasingx restricts
l p c .
thetotalnumberofconfigurationsavailabletothesystem
by increasing the lower bound on the accessible volume
fluctuations. However, it also seems to imply that the
hopping process for an ideal gas particle is “barrierless”,
∗
but this is not the case. The free energy barrier ∆G ,
4.25
givenbyEq.21,isnecessarilypositivebecauseitinvolves
anintegralofthe normalizedprobabilityovera smallre-
gion of the reaction coordinate in the transition region. op 4
This indicates that it takes work to restrict the two par- τh
ticles tothe transitionstate regionrelativetohavingthe ln
cagedparticle located anywherewithin the small system 3.75
volume. Figure 2 also shows that the probability den-
sity obtained from our simulations matches the results
3.5
obtained from our analytical analysis. While this is not
surprising for such a simple problem, it provides a proof
of principle and suggests that our simulation approach
3.25
would be applicable to systems involving more compli-
3.25 3.5 3.75 4 4.25 4.5
cated particle–particle and particle–wall interactions. β∆G*
Figure 3 shows a plot of lnτ , obtained from our
hop
hopping time simulations versus the free energy barrier
∗
to hopping obtained from the analytical results over a FIG. 3. lnτhop as a function of β∆G for the 2D ideal gas
range of channel radii, where we have used Eqs. 31 and over a range of pore radii, Rp = 1.1−4.0. The dashed line
∗ represents the best linear fit to the data and has a slope of
21, with x = 0.005. The slope is one, which sug-
c 0.98. The error bars represent the standard deviation in the
gests that our TST approach, using the small system
measured hoppingtimes.
isobaric–isothermalensemble,providesapotentiallyuse-
ful method for the calculation of hopping times in these
quasi-one dimensional systems. The intercept yields the
kinetic prefactor κ = 0.90 and essentially measures the
B. Two Dimensional Hard Discs
time it takes for the particles to move out of the transi-
tionstateregion. Equation32,withEq.14,suggeststhat
τ shouldvarylinearlyasR−1,withaslope∼κ/2βp x∗ We now calculate the hopping barrier for a system
hop p l c
(seeFig.4). Thisdiffersfromthescalinglawobservedfor of two dimensional hard discs confined to a channel by
hard disks in the narrow channel regime and provides a hard walls, where the particle-particle exclusion interac-
limitingcaseforthehoppingtimeinwidechannelswhere tion makes it more difficult for particles to escape their
the particles can easily pass. cage as the channel radius decreases. For the hard disc
8
interaction is given by,
0 if |y |≤R −σ/2,
i p
75 U (rˆ )= (34)
WHS ij
(∞ if |yi|>Rp−σ/2,
wherey isthey–coordinateoftheparticle. Tomakethe
i
solution to the small system n,p,T partition function
p
o50 tractable,weconsiderthe casewherethepressureislow.
h
τ
The interactions between the particles inside the sub-
volumeandthoseinthe surroundingsarethennegligible
and we can assume U = 0. The reaction coordinate
σ
partition function for this system is then given by,
25
0.2 0.4 0.6 0.8 1 ∆ dx = 2dxc ∞e−βpl2RpLdL Rp−σ2 dy y1−σ2 dy
1/Rp xc c Λ4 Zxc Z−Rp+σ2+σ2 1Z−Rp+σ2 2
(2R −σ−σ )2
= p 2 e−(2Rp−σ)βplxcdx ,
FIG.4. τhop,asafunctionoftheinversechannelradius1/Rp, Λ4βpl(2Rp−σ) c
for the2D idealgas system. Thedashed lineisalinear fitto (35)
the data with a slope of 83.6. The error bars represent the
standard deviation in themeasured hopping times.
where,
system the particle–particle interaction is defined as, σ = σ2−x2c if |xc|≤σ, (36)
2
(p0 if |xc|>σ,
0 if r ≥σ,
U (r )= ij (33)
HD ij
(∞ if rij <σ, and the factor of 2 appears because we need to include
the configurations in which the order of particles 1 and
where σ is the particle diameter and r =|r −r | is the 2 in the y-coordinate are reversed. Taking into account
ij i j
the distance between two particles. The particle–wall thepiecewisenatureofσ andintegratingoverx yields,
2 c
∆= 1 1+ 2e−βpl(2Rp−σ)[βplσ(2Rp−σ)+1]+β2p2lσ2(2Rp−σ)2−2 − πσ[In(z)−Ln(z)] , (37)
β2p2lΛ4 " (2Rp−σ)4 2Rp−σ #
where I (z) is the modified Bessel function of the first accessible configurations available to the system as the
n
kind,andL (z)isthemodifiedStuvefunction,bothwith transitionstate isapproached. Conversely,the pressure–
n
n=1 and z =βp σ(2R −σ) [55]. volume term increases the number of accessible config-
l p
urations because the range of allowed volume fluctua-
Figure 5 shows the free energy density, −lnP(x ), as
c tionsincreasesasx →0. Thefreeenergydensityprofile
c
a function of x for channel radii in the range R /σ =
c p obtained here, in the small system isobaric-isothermal
1.01−3. Forx >σ,thetwodiscsdonotexcludevolume
c ensemble, differs significantly from that obtained in the
with respect to each other (σ = 0) and they effectively
2 canonical ensemble where the free energy maximum is
behave as idealgas particles in a channelwith a reduced
always located at x =0 [37].
c
radius of R − σ/2, which gives rise to the linear in-
p
crease in −lnP(x ). When x < σ, the two discs begin
c c
to exclude volume, leading to a restriction in configura-
tionspace. However,whenthe channeldiameter iswide, In the context of particle motion, the transition state
the volume exclusion effect is small and the free energy forhoppingisstilllocatedatx =0,despitethepresence
c
density remainscloseto being linear. As the channelbe- of a maximum in the free energy density at values of
comes narrower, we see the appearance of a maximum x > 0, because it marks the point where the particles
c
that begins at x = 1 and moves towards x = 0 as exchange position, signifying the crossover from single–
c c
R → σ. This unusual feature results from the convo- file to normal diffusion. The probability of finding the
p
lution of two competing effects. The excluded volume cagedparticle within a regiondx of the transitionstate
c
interaction between the two discs reduces the number of can be obtained using Eqs. 35 and 37 in Eq. 19, and
9
4
10
8
0)
6 P(1
) 102 1/
c
(x 4 R /σ 0)
P p (
n /P 1 2R -σ 10
-l 2 1.01 1 p
1.10 100
1.20
0 1.70
3.00
-2 10-2
0 2 4 6 8 -2 -1 0 1
x 10 10 10 10
c R -σ
p
FIG. 5. Free energy density, −lnP(xc) as a function of xc, FIG. 6. Log–Log plot of 1/P(0) versus the Rp−σ showing
with different channel radii and a fixed pressure, βpl = 1.0, thepowerlaw behaviourof τhop for harddiscsas thepassing
forthe2D harddiscsystem. Thesolid linesrepresent results threshold is approached. The solid line shows the theoreti-
from the theory and the points represent the data obtained cal results obtained from Eq. 38, the points are the results
from simulation. from our simulation free energy calculations and the dashed
. line shows the limiting slope of −2. Insert: Log–Log plot
of 1/P(0) versus 2Rp −σ showing the wide channel scaling
law for τhop obtained from theory (solid line) and simulation
noting σ =σ (Eq. 36) when x =0, to give, (points). Thedashedlinehighlightsthelimiting slopeof −1.
2 c
4(R −σ)2dx
p c
P(0)dx = . (38)
c Λ4βp (2R −σ)∆
l p
12
The hopping time in the hard disc model is expected
to diverge as a power law, ∼ (R −σ)α, as the chan-
p
nel radius approachesthe diameter of a disc because the
10
excluded volume interaction prevents the two particles
from reaching the transition state. Once the channel ra- p
o
tdhiuesdbiseccsobmeecsosmmeaplleerrmthananentthlyepcaagsseidngbetthwreesehnotlhde,iRrpne=ighσ-, τn h 8
l
bours. AFicks–Jacobsanalysis[56]thatprojectsthedif-
fusion of the two discs onto a one-dimensional reaction
coordinate predicts [57, 58] α = −3/2, while TST pre- 6
dicts [38] α = −2. Simulation studies of the hopping
time agree with the TST result, but the scaling only be-
4 6 8 10 12
comes apparent for very narrow channels [36, 59]. Tran- β∆G*
sition state theory suggests τ ∼ 1/P(0) and Eq 38
hop
shows that the hopping time shoulddiverge with the ex-
ponent α = −2, as expected. Figure 6 also shows how FIG. 7. lnτhop as a function of β∆G∗ for hard discs for a
1/P(0) crosses over from the passing threshold limit for range of different pore radii, Rp = 1.01−1.20. The points
narrow channels to the wide channel scaling behaviour represent the data obtained from simulation with xc =0.03.
∼(2R −σ)−1, whichis consistentwith ouridealgasre- The dashed line has a slope of 1 and provides a guide to the
p
eye.
sults. Finally,Fig.7showsthatEq.14workswellforthe
harddiscsystemwhenthe channelradiusisnarrow,giv-
ingrisetohighbarriers. Thisisalsotheregimewherethe
V. DISCUSSION
hoppingtimeisdirectlyconnectedtothediffusioncoeffi-
cientthroughEq.3becausetheparticlemustbetrapped
longenoughtoperformSFDbetweenhops. However,our The phenomenological hopping time approach to un-
TST approach begins to break down for wider channels derstanding diffusion in single file fluids confined within
where we appear to underestimate the height ofthe bar- quasi–one dimensional channels is attractive because all
rier needed to predict the hopping time. theeffectsthatinfluenceparticlemotioninthelongtime
10
limit are contained within a single parameter, τ , that the channel becomes narrower. This assumption also re-
hop
describes the shorttime event of a particle hopping past structs the application of our method to low pressures,
one of its neighbours. The activated nature of the hop- where again, interactions between the small system and
ping processalsomakesthe calculationofτ amenable the surroundings are minimal.
hop
totheapplicationoftransitionstatetheory,butthechal- Thehoppingtimeisonlyconnectedtothediffusionco-
lenge is to provide a rigorous approach for the calcula- efficient through Eq. 3, with γ =1/2, when the particles
tion of the probability of being in the transition state, are caged long enough between hops that they perform
i.e. the height of the hopping free energy barrier, for single file diffusion during this time, which should occur
a nanoscale system that exhibits volume fluctuations. when the barrier is high. This suggests that our method
While the isobaric–isothermal ensemble is the appropri- should be effective in predicting the hopping time and
ate choice of ensemble, a direct application of the bulk dynamic properties of single file fluids in the crossover
thermodynamic results to a system in the nanoscale is regimes and there are a wide range of important practi-
problematic due to the system size dependence of the cal challenges in nano-fluidic systems, including the sep-
necessary volume scale. arationof mixtures, that can be addressed in the regime
where our method is valid and simple enough to provide
We have shown that this problem can be overcome by
usefulinsight. However,MonandPercus[26]alsoshowed
usingthesmallsystemisobaric–isothermalensemblepro-
thatEq.3maybeapplicableinthewidechannelregime,
posedbyCorti,wheretheuseofashellparticletodefine
but with a different hopping time exponent. When the
the volume prevents the over counting of configurations
channel is wide, the particles can pass each other eas-
whenconsideringvolumefluctuations andeliminates the
ily because the barrier is low. The hopping time is then
need for an explicit volume scale. Figures 3 and 7 high-
shortsothatthedistancetraveledbetweenhopsbecomes
light the key results of our work and show that lnτ
hop independentofthehoppingtime,andproportionaltothe
∗
varieslinearlywith∆G withthe slope ofunity forboth
average distance between nearest neighbors in the chan-
systems studied here, over a range of different radii, as
nel, which leads to an exponent of γ = 1. This has
predictedbyTST.Oneofthemainadvantagesofourap-
notbeentested,butsuggeststhehoppingtime approach
proachis thatitonlyrequiresthe simulationoftwofluid
may still be applicable even when the channel is wide.
particles and the wallinteractions,while the direct mea-
To make our TST approach work in the wide channel
surementofthediffusioncoefficientandτ involvesthe
hop regime, it would be essential to account for the interac-
simulation of hundreds of particles and requires a long
tions between the small system and the surroundings in
simulation time because τ is slow to converge. The
hop our simulation method.
small number of particles also means that the probabil-
ity density alongthe reactioncoordinatecanbe sampled
by“bruteforce”becausesimulationswithalargenumber
VI. CONCLUSIONS
of MC cycles can be performed quickly, but the method
canalsobeadaptedtoincludemoreadvancedtechniques
In this work, we have shown that the small system
such as biased umbrella sampling [53] in cases where the
isobaric–isothermal ensemble provides a rigorous frame-
barrier is large.
workfor the developmentof a transitionstate theory for
It is interestingto note that ourTST approachbreaks the calculation of the hopping time and we have applied
down for wider channels in the hard disc model, but ittotwosimplesystems. Whilethesimulationmethodis
works well over all channel radii in the ideal gas model currently restricted to the narrow channel, low pressure
even when the barrier is only a few kT. The main dif- regime,whereinteractionsbetweenthesmallsystemand
ference between the two models is the exclusion inter- its surroundings can be neglected, it is simple to imple-
action of the hard discs and the effect this has on our mentandis potentially applicable to the study ofa wide
decision to neglect the interactions between the small range of systems.
system and its surroundings. Setting U = 0 is exact
σ
for the ideal gas system. In the hard discs model, this
allowsthetwoparticlesinsmallsystemtosampleconfig- ACKNOWLEDGMENTS
urations in the transition state that would otherwise be
ruled out by the exclusion interaction with particles in WewouldliketothankNaturalSciencesandEngineer-
the surroundings sitting near the boundary. As a result, ingResearchCouncilofCanada(NSERC)andPOSBio-
our method over estimates the partition function of the Sciences for financial support. Computing resourcesand
transition state which leads to a lower predicted barrier. support were provided by Compute Canada and West-
Neglecting these interactions becomes less important as Grid.
[1] L. Gelb, K. Gubbins, R. Radhakrishnan, and systems,” Rep.Prog. Phys. 62, 1573–1659 (1999).
M. Sliwinska-Bartkowiak, “Phase separation in confined
