Analytical study of tunneling times in flat histogram Monte Carlo

2 MB·English

by  Miguel D. Costa

EurophysicsLetters PREPRINT
6 Analytical study of tunneling times in flat histogram Monte Carlo
0
0
2 MIGUELD. COSTA1,2,J. VIANA LOPES1,3(∗)andJ.M.B. LOPESDOS SANTOS1
n 1 CentrodeFísicadoPorto,DepartamentodeFísica,FaculdadedeCiências,Universidadedo
a Porto,4169-007PortoPortugal
J 2 EscolaSuperiordeTecnologiaeGestão,InstitutoPolitécnicodeVianadoCastelo,Vianado

4
CasteloPortugal

3 Departamento de Física, Instituto Superior de Engenharia, Instituto Politécnico do Porto,
]
h PortoPortugal
c
e
m
PACS.75.10.Hk–Classicalspinmodels.
-
t PACS.02.70.Rr–Generalstatisticalmethods.
a
PACS.64.60.Cn–Order-disordertransformations;statisticalmechanicsofmodelsystems.
t
s
.
t
a
m
- Abstract. – Wepresentamodelforthedynamicsinenergyspaceofmulticanonicalsimulationmethods
d thatlendsitselftoarathercompleteanalyticcharacterization. Thedynamicsiscompletelydeterminedby
n thedensityofstates. Inthe±J 2Dspinglassthetransitionsbetweenthegroundstatelevelandthefirst
o
excitedonecontrolthelongtimedynamics. Weareabletocalculatethedistributionoftunnelingtimes
c
andrelateittotheequilibrationtimeofastartingprobabilitydistribution. Inthismodel,andpossiblyin
[
anymodelinwhichenteringandexitingregionswithlowdensityofstatesaretheslowestprocessesinthe

2 simulations, tunnelingtimecanbemuchlarger(byafactorofO(N))thantheequilibrationtimeofthe
v probabilitydistribution. Wefindthatthesefeaturesalsoholdfortheenergyprojectionofsinglespinflip
8 dynamics.
9
6
4
0
5 Thesimulationofmodelswithcomplicatedenergylandscapes,suchasspinglassesormodelsof
0
proteinfolding,hasalwaysprovedimpracticalforconventionalcanonicalMonteCarlomethods[1].
/
t InaMonteCarlosimulationofasystemofN degreesoffreedom,atfixedtemperatureT,energy
a
m fluctuationsabouta mean energyE(T) are of order kBT√Ncv. As a result, only states in a range
of energies E δE < E < E +δE with δE k T√Nc are accessible in the simulation.
- h i− h i ∼ B v
d Insystemswith complexenergylandscapes,phasespacein thisnarrowenergyrangebreaksupinto
n many regions, connected only by states requiring energy fluctuations ∆E δE; tunneling times
o ≫
betweentheseregionsbecometoolongtoretainanyhopeofachievingtheasymptoticdistributionin
c
: areasonablesimulationtime.
v
To overcomethis problem, one would like to broadenthe energyrange of the states sampled in
i
X a simulation, forsaking the canonical ensemble at a fixed temperature. A variety of methods have
r been proposed to implement this idea, such as entropic sampling [2], multicanonical Monte Carlo
a [3,4], simulated and paralleltempering[5–7], Wang-Landausampling[8–10], broadhistogramand
transitionmatrixmethods[11–14].Insomecases,likemulticanonicalMonteCarloandWang-Landau
sampling, the aim is to sample all energy levels of the spectrum with equal probability, producing
(∗)E.mail:[email protected]
(cid:13)c EDPSciences
2 EUROPHYSICSLETTERS
flathistogramsin energyspace. Whilethisallowsovercominglargeenergybarriersinsystemswith
ruggedenergylandscapes,criticalslowingdowncanremainaproblem[14,15].
Given the practical impossibility of measuring equilibration time of an initial distribution, the
performanceofthesemethodsisgenerallyassessedbymeasuringthesocalledtunnelingtime[15–17],
thetimerequiredtocrosstheentireenergyspectrum,fromgroundstates,tostateswiththeenergyof
maximumdensityofstates, andback,orfromgroundstatetoanti-groundstates, statesofmaximum
energy.
Inthisletterwefocusonthedynamicsofthisrandomwalkinenergyspace. Weproposeamodi-
ficationofthedirectednetworkofnon-zerotransitionprobabilitieswhichunderliestheMarkovchain
ofaMonteCarlosimulation,which,whileretainingthelongtimebehaviouroftheenergyprojection
ofconventionalsimulations,lendsitselftoa verycompleteanalyticdescription. Ourthreemostim-
portantresultsconcernthe J 2Dspinglass: (a)thetransitionsbetweengroundstate levelandfirst
±
excitedonecompletelycontrolthelongtimedynamicsofsimulationsinthismodel;(b)theequilibra-
tiontimeofasimulationcanbeconsiderablyshorterthantunnelingtime; (c)theseconclusionsalso
holdforthecaseofconventionalsinglespinflip(SSF)dynamics.
To calculate thermodynamic averages in any system, we need only specify the set of points of
phasespaceandthecorrespondingprobabilities.TosolvethisproblemusingMonteCarlo,weimpose
onthisstructureadirectednetworkofnon-zerotransitionrateswhichdefinestheMarkovchainused
to generate the asymptotic distribution. The traditional SSF connects each point of phase space to
N first neighbours. The SSF algorithm is, usually, local in energy, i.e., energydifferencesbetween
connected phase space points being of O(1), independentof system size. It has been claimed that
differencesinthescalingbehaviourofthe J 2Dspinglassmodelandafullyfrustratedmodel,both
±
withextensiveground-stateentropy,resideintherestrictionsimposedbythenetworkdefinedbySSF
dynamics.Still,thescalinglawsofrelevanttimescaleswithsystemsizefoundwithourdynamics[16]
arealsosimilartothoseobtainedwithSSFdynamics[15].
ThesimplestchoicethatcircumventstheproblemoflocalminimaistoconstructaMarkovchain
that connects each point of phase space to all the remaining points that have the same or adjacent
energies. Inthisway,wepreservethelocalityinenergyspaceofthestandardSSFnetworkandavoid
theneedforlongpathstoconnectstatesthatareadjacentinenergy.
Theprocedureinbestexplainedinreferencetoaspecificmodellikethe J 2DIsingspinglass.
±
The correspondingenergylevels, E , and degeneracies, N(E ), can be calculated exactly using the
i i
programofSaulandKardar[18]. Non-zerotransitionratesconnectstateswithenergyE toallstates
i
withenergyEi orEi±1. Thesimplestchoiceoftransitionratesistoproposethefinalstateuniformly
fromthesetofallowedstates,andaccepttheproposedstatewithaprobabilitywhichensuresthatthe
asymptotic probability of each state equals 1/N(E). Whereas in a SSF simulation most transition
probabilitiestostatesofthesameornearbyenergyarezero,inourmodeltheyareallnonzero,with
valuesadjustedtoensurethesameequilibriumdistribution. Thisprocedureisakintoanaveragingof
transitionrates.Whilethiscouldleadtosignificantchangesinthedynamics,weneverthelessfindthat
themostrelevantfeaturesofthelongtimeSSFdynamicsinenergyspacestillholdinourmodel.The
advantage of the current model is that the projection of the multicanonicalMarkov chain in energy
spaceremainsaMarkovprocess,allowingustheuseofasetofanalytictoolstostudythedynamics.
Wedenotebyp (t)theprobabilityofbeinginastateα;N isthedegeneracyofenergylevelE
α i i
andMi Ni−1+Ni +Ni+1 1isthenumberofstatesaccessibleinadirecttransitionfromany
≡ −
stateofenergyE (thechoiceoftentativefinalstatesexcludesthecurrentstate). Inequilibrium,the
i
flatenergyhistogramconditionrequires,forastateαofenergyE ,
i
1 1
p = (1)
α
N(E ) ≡ N
α i
Forreasonsofefficiency,weexcludethenegativetemperatureregion: forenergiesgreaterthanE ,
m
MIGUELD.COSTAetal.:ANALYTICALSTUDY OF TUNNELINGTIMES IN FLATHISTOGRAM... 3
theenergyatthemaximumofdensityofstates,wereplaceN(E)byN(E )ineq.1. Thetransition
m
ratefromastateαtoanaccessiblestateβ,ω ,satisfiesdetailedbalance,withtheequilibriumdistri-
βα
butiongivenineq.1;asusual,itiswrittenasaproductofaproposalprobability,whichwetaketobe
uniformamongallpossible finalstates, and anacceptanceratio, determinedbythe detailed balance
condition:
1 M(E )N(E )
ω = min 1, α α . (2)
βα
M(E ) (cid:18) M(E )N(E )(cid:19)
α β β
Inthemasterequation,thetransitionprobabilitiesω donotchangeifwevaryβ (orα) withinthe
βα
sameenergylevel.Thismakesitpossibletowritethefollowingmasterequationfortheprobabilityof
havingenergyE ,P = p δ
i i α α Eα,Ei
P
Pi(t+1) = ΓiPi+1(t)+Γi−1Pi−1(t)
+ (1 Γi Γi−1)Pi(t) (3)
− −
= Ω P (t) (4)
ij j
Xj
with
Γ = Ni min 1,Ni+1Mi+1 . (5)
i
Mi+1 (cid:18) MiNi (cid:19)
This equation defines the projection of the original Markov chain onto to the energy variable; the
random walk in energy space remains a Markov chain with transition probabilities defined by the
density of states of the original spin model. This is,of course, an important simplification afforded
by our choice of transition rates. With the exception of special models, like the infinite range Ising
model[19],theprojectionofconventionaldynamicsinenergyspace isnon-markovianandmemory
effectscanhaveasignificantimpactonthedynamics[14,15].Nevertheless,inthemodelunderstudy,
wefindbelowthatimportantaspectsofthelong-termSSFdynamicsarepreserved.
Assumethesystemstartsfromanenergylevelrattimet=0:
P (0)=δ
i i,r
andletQ (t)nowrefertotheprobabilityofhavingenergyE attimet,giventhatthesystemhasnever
i i
visited energylevels (s > r). Q (t) satisfies the master equationin eq. 3 fori < s, with the same
i
initialconditionasP (t),namelyQ (0)=δ ,andanadditionalconditionQ (t)=0replacingeq. 3
i i i,r s
fori=s. Theprobabilityoffirstpassageinsbecomes
Hsr(τ)=Γs−1Qs−1(τ). (6)
The master equationcan be solved using a normalmode expansion,Q (t) = a fγλt. The di-
i γ γ i γ
mensionofthetransitionmatrixinenergyspace,Ω,scaleslinearlywithsystemsPize,notexponentially
as in thecase of thetransitionmatrixof theMarkovprocessin phasespace. Thediagonalizationof
Ω,becomesamanageableproblemallowingthecalculationoftheeigenvalues,λ ,theleftandright
γ
eigenvectors,gγandfγ,andthecoefficientsa = gγ/ gγfγ. Itisalsopossible,ofcourse,todi-
i i γ r i i i
agonalizethetransitionmatrixfortheactualprobabilityPdistributionPi(t)usingthesameformalism,
andaccessthedecaytimeofthevariouseigenmodesofthemasterequation.Thelargestfiniteone(the
equilibriumdistribution,P ,doesnotdecay)istheequilibrationtimeofP (t),τ .
eq i eq
Thedistributionoftunnelingtimemeasuredinoursimulationscanbewrittenintheform
τ
′ ′
H(τ)= H0m(τ τ )Hm0(τ ) (7)
−
τX′=0
4 EUROPHYSICSLETTERS
P(τ/N2) L=6 P(τ/N2) L=12
EExxaacctt PPddff
0.004 MMCC HHiissttooggrraamm
0.03
0.003
2N)) 2N))
0.02 τlog(P(/ 0.002 τlog(P(/
0.01
τ/N2 0.001 τ/N2
0 0
0 20 40 60 80 100120 0 200 400 600 800 1000
τ/N2 τ/N2
Figure1–ComparisonbetweenexactresultgivenbymasterequationapproachandtheMonteCarloresultsfor
twospinglasssampleswithL=6andL=12.
where E is the energy level corresponding to the maximum of N(E). In Fig. 1 we superpose a
m
distributionmeasuredinaMonteCarlosimulationwiththeratesgivenbyeq. 2withacalculatedone
fortwo spin glasssamplesof linear dimensionL = 6, 12; the simulationreproducesthe calculated
distributionveryaccurately.
The distribution of tunneling time between any two energy levels decays exponentially with a
time constant given by τ = 1/log(λ ), where λ < 1 is the largest eigenvalue of the
max max max
correspondingΩ(eq.6).Foragiv−enspinglasssample,wedenotebyτ andτ thelongestdecaytime
u d
fortunnelingE0 Em andEm E0,respectively;H(τ)(eq.7)willdecayexponentiallywiththe
→ →
largestof τ andτ forthe givensample. The powerlaw decayseen in [15]appearsonlywhen we
u d
aggregatetunnelingtimesfromdifferentrealizationsoftherandominteractions J.Ifthedistribution
ofthelargestdecaytime,intheensembleofspinglasssamples,hasapowerlaw±tail,ρ(τ ) τ−ν
forτ ,weobtain,H(τ) dxρ(x)e−τ/x t1−ν ast . max ∼ max
max
→∞ ∼ ∼ →∞
This asymptotic behaviour is dirRectly related to features of the density of states of the J spin
±
glass. Theinversetransitionratefromthegroundstateleveltothefirstexcitedlevelisgivenby(eq.
5),
τ0 =Γ−01 =1+ N1N+0N2.
It was noted in [15] that N1/N0 has a Fréchet probability distribution function, in an ensemble of
spinglasssamples;asimilarresultoccursforτ0 (Fig.2). Thisfactcompletelycontrolsthelongtime
dynamicsoftheMonteCarlosimulations.
InFig.3weplotτ ,theequilibrationtime,andτ andτ ,thedecaytimeconstantsfortunneling
eq u d
fromE0toEmandfromEmtoE0,againstτ0. Inthesampleswithlargerτ0,weobservethefollowing
relations:
τeq τ0 (8)
≈
τu τ0. (9)
≈
Notethatτd canbetwoordersofmagnitudelargerthanτu;tunnelingfromEm toE0 ismuchslower
than from E0 to Em. In fact, as can be seen in the inset of Fig. 3 (b), the following relation holds
asymptotically:
τ N τ (10)
d b u
≈
MIGUELD.COSTAetal.:ANALYTICALSTUDY OF TUNNELINGTIMES IN FLATHISTOGRAM... 5
8×10-5
6×10-5 τ))0
P(
τ)04×10-5 log(
(
P
τ
0
2×10-5
Histogram
Fit
0
0 2×104 4×104 6×104 8×104 1×105
τ
0
Figure2–Comparison betweenthehistogramof τ0 from10000 uniformlygenerated samplesof ±J 2D spin
glasswithL = 10andthemaximumlikelihoodfittoaFrechetdistributionfunctionobtainedfromtheoriginal
data.
whereNb O(N)isthenumberofenergylevelsbetweenE0 andEm (recallthatinoursimulation
∼
theenergyhistogramislimitedtothisenergyrange). Theseresultsarequitesimpleconsequencesof
theexistenceofabottleneckinthedynamicsduetotransitionsbetweenE0andE1(largeτ0).
TheequationforP0is(forlongtimewereplacePi(t+1) Pi(t)bydPi/dt)
−
dP0(t)
= Γ0(P0(t) P1(t)). (11)
dt − −
We denote the tunnelingtime scale between E1 and Em by τ1 ; for sampleswith verylarge τ0, we
assumeτ0 τ1. Fort τ1,wewillhavefori=0
≫ ≫ 6
1 P0(t) Nb
Pi(t) − = Peq(1 P0(t))
≈ N 1 N 1 −
b b
− −
whereP = 1/N istheequilibriumdistribution. Replacingineq.11weobtainatimeconstantfor
eq b
thedecayofP0(t) Peq givenby
−
−1
N
τeq =(cid:18)Γ0N b 1(cid:19) ≈Γ−01.
b
−
TocalculatethedistributionoftunnelingtimefromE0 toEm,H(τ),weusetheinitialcondition
Pi(t) = δi,0. DefiningPi(t) = Qi(t)+Ri(t), where Qi(t) is the probabilityof having energyEi
attimet giventhatthesystem hasnotreachedE , Q(t) = Nb−1Q (t)istheprobabilitythatthe
∞ m i=0 i
tunnelingtimeisgreaterthant,Q(t)= dτH(τ). Fort Pτ1wehave
t ≫
R
P0(t) Q0(t)
≈
P (t) R (t) fori=0
i i
≈ 6
This expresses the fact that the system either has energy E0, and has not tunneled, or has left the
groundstateandhasalmostcertainlytunneledtoE . Therefore
m
dQ(t) dQ0 dP0
dt ≈ dt ≈ dt
6 EUROPHYSICSLETTERS
τ
eq
109 τeq=τ0 τd/Nb
106 L = 20 (a)
103 L = 6
τ 103 106 109 τ0 106
d,u τu=τ0 τd L=20
108
L = 20
τ
L = 6 u (b)
104
104
100 τ/N=τ
103 106 109 τ d b u
0
τ /N
d b
106 (c) 102 L=6
103
103 106 τ 102 104 106 τ
u u
Figure3 Figure4
Figure3–(a)Correlationbetweenequilibrationtimeτeq andτ0 forlatticesizeL = 6to20. Thedashedline
showstheasymptoticbehaviourwhereτeq = τ0. (b)Correlationbetweenbothτd andτu withτ0. Thedashed
lineshowstheasymptoticbehaviour whereτu = τ0. (c)Correlationbetweenτd/Nb andτu. Thedashedline
showstheasymptoticbehaviourwhereτd/Nb =τu.
Figure4–Correlationbetweenτd/Nb andτu obtainedfromMonteCarlosimulationsofthe±J 2Dspinglass
withSSFdynamics.Thedashedlineshowstheasymptoticbehaviourwhereτd/Nb =τu.
ThedecaytimeofQ(t)(andH(t))isτu τ0.
≈
Theasymmetryofτd andτu shouldbynowbeobvious. FortunnelingfromEm toE0 theinitial
conditionisPi(0)=δi,m andQ(t)= mi=1Qi(t). Fortimesmuchlargerthanτ1 weexpect,now,
P
Q(t)
Q (t) . fori=0
i
≈ N 1 6
b
−
SincedQ(t)/dt= Γ0Q1(t)(Q(t)onlychangesduetotransitionsbetweenE1andE0),weobtain
−
dQ(t) Γ0
= Q1(t),
dt −N 1
b
−
i.e.
τ =(N 1)τ N τ .
d b u b u
− ≈
In simple language, these results can be understood as follows. Tunneling from E0 to Em is
controlledbytheprocessofexitingthegroundstatelevel:thesystemcannothavetunneledifitisstill
inagroundstate. Ontheotherhand,asystemcanonlyenterthegroundstatelevelifitisinenergy
levelE1; thereforethetunnelingratefromEm toE0 hasanextrafactorPeq = 1/Nb corresponding
totheprobabilityofhavingenergyE1.
These results provide a very clear explanation of the correlation between tunneling time and
N1/N0 found in [15]. The controlling time scale τ0 = 1 + (N1+N2)/N0 is clearly related to
theratioN1/N0. Thisdifferenceintimescalesτu andτd shouldbepresentinothermodels;namely,
MIGUELD.COSTAetal.:ANALYTICALSTUDY OF TUNNELINGTIMES IN FLATHISTOGRAM... 7
thoseinwhichtheslowestprocessesinthesimulationsinvolvegettingacrosssteepentropychanges.
The2D J spinglassisanextremecaseofthisbehavior,thedynamicsbeingcontrolledbythetran-
±
sitionbetweenthetwolowestenergylevels. Onecouldaskifthisresultisanartifactofoursimplified
dynamics. Infact,weverifiedthesamebehaviorwithSSFdynamics. Fig.4showsthesametypeof
plotasinFig.3withaveragedtunnelingtime(averagestakenwithinasample)inSSFdynamics: for
longtimestherelationofeq. 10is stillverified. Noticethatthe equilibrationtime scale ofaproba-
bilitydistribution,τ ,maybemuchsmallerthantunnelingtimescale,which,asusuallymeasured,is
eq
dominatedbyτ . Thiscanleadtoaverypessimisticestimateofthetimerequiredtoreachequilibrium.
d
Insummary,modellingthedynamicsofamulticanonicalsimulationinasortofmean-fieldway,
byaveragingtransitionratestostatesofthesameornearbyenergies,wewereabletodefineaMarkov
process in the energy variable, reducing the dimension of the Markov matrix to a manageable size.
Nevertheless,thisprocedurepreservesthemainfeaturesofthelongtimedynamicsofaconventional
simulation. A very complete and physically transparent description of the more salient features of
thedynamicsofmulticanonicalsimulationsbecomespossible. Inparticular,weclarifiedtherelation
betweenequilibrationandtunnelingtime,andthedifferencebetweentunnelingawayfromregionsof
lowdensityofstatesfromtunnelingintosuchregions.
∗∗∗
The authors would like to thank E.J.S. Lage and J. Penedones for very helpful discussions and
S.SabhapanditforprovidinganimplementationofSaulandKardar’salgorithm. Thisworkwassup-
portedbyFCT(Portugal)andtheEuropeanUnion,throughPOCTI(QCAIII).Twooftheauthors,JVL
andMDC,weresupportedbyFCT grantsnumbersSFRH/BD/1261/2000andSFRH/BD/7003/2001,
respectively.
References
[1] NEWMANM. E.J.andBARKEMAG. T.,MonteCarloMethodsinStatisticalPhysics(OxfordUniversity
Press)1999.
[2] LEEJ.Y.,Phys.Rev.Lett.,71(1993)2353.
[3] BERGB.A.andNEUHAUST.,Phys.Rev.Lett.,68(1992)9.
[4] BERGB.A.andNEUHAUST.,Phys.Lett.B,267(1991)249.
[5] MARINARIE.andPARISIG.,Europhys.Lett.,19(1992)451.
[6] HUKUSHIMAK.andNEMOTOK.,J.Phys.Soc.Jpn,65(1996)1604.
[7] LYUBARTSEVA.P.,MARTSINOVSKIA.A.,SHEVKUNOVS.V.andVORONTSOVVELYAMINOVP.N.,J.
Chem.Phys.,96(1992)1776.
[8] WANGF.G.andLANDAUD.P.,Phys.Rev.E,64(2001)056101.
[9] WANGF.G.andLANDAUD.P.,Phys.Rev.Lett.,86(2001)2050.
[10] TROYERM.,WESSELS.andALETF.,Phys.Rev.Lett.,90(2003)120201.
[11] OLIVEIRAP.M.C.,PENNAT.J.P.andHERRMANNH.J.,Braz.J.Phys.,26(1996)677.
[12] WANGJ.S.andSWENDSENR.H.,J.Stat.Phys.,106(2002)245.
[13] WANGJ.S.,TAYT.K.andSWENDSENR.H.,Phys.Rev.Lett.,82(1999)476.
[14] TREBSTS.,HUSED.A.andTROYERM.,Phys.Rev.E,70(2004)046701.
[15] DAYALP.,TREBSTS.,WESSELS.,WURTZD.,TROYERM.,SABHAPANDITS.andCOPPERSMITHS.
N.,Phys.Rev.Lett.,92(2004)097201.
[16] LOPESJ.V.,COSTAM.D.andLOPESDOSSANTOSJ.M.B.,inpreparation,.
[17] BERGB.A.andCELIKT.,Int.J.Mod.Phys.C,3(1992)1521.
[18] SAULL.andKARDARM.,Nucl.Phys.,432(1994)641.
[19] WUY,KOERNER,M.,COLONNA-ROMANOL.,TREBSTS.,GOULDH.,MACHTAJ.andTROYERM.,
condmat/0412076,(2004).