Table Of Content

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