Table Of ContentFlat histogram version of the pruned and enriched Rosenbluth method
Thomas Prellberg and Jaros law Krawczyk
Institut fu¨r Theoretische Physik, Technische Universita¨t Clausthal,
Arnold Sommerfeld Straße 6, D-38678 Clausthal-Zellerfeld, Germany
(Dated: February 2, 2008)
4
0 In this letter we present a flat histogram algorithm based on the pruned and enriched Rosen-
0 bluth method (PERM). This algorithm incorporates in a straightforward manner microcanonical
2 reweighting techniques,leading to“flat histogram” sampling in thechosen parameter space. Asan
n additional benefit, our algorithm is completely parameter free and, hence, easy to implement. We
a applythisalgorithmtointeractingself-avoidingwalks(ISAW),thegenericlatticemodelofpolymer
J collapse.
2
1
Recently, there has been revived interest in flat his- sider a set of variables mi, but for pedagogical reasons
] togramalgorithms[1],whichstrivetoevenlysamplecon- we shall in this paper restrict ourselves to the case of m
h
figuration space with respect to a chosen parameteriza- corresponding to an energy E = ǫm. Both n and m are
c
tion, e.g. microcanonical energy. These algorithms are assumed to have non-negative integer values. Moreover,
e
m particularimplementationsof“umbrellasampling”[2],in we will need the notion of “atmosphere” of a configura-
which the configuration space is sampled according to a tion, which is the number a of different waysto continue
-
t given probability distribution, the so-called “umbrella”. to growthis configuration,and is alsoa non-negativein-
a
t This umbrella distribution is generally chosen such that teger.
s
the whole configuration space of interest is accessible in While it is useful to present the algorithm in such an
.
t one simulation. One major difficulty hereby is to find a abstract setting, it may help the reader to keep the ap-
a
m suitable umbrella distribution. plication to interacting self-avoiding walks (ISAW) on a
There has also been an exciting development in regularlatticeinmind. Inthiscase,thesizeoftheconfig-
-
d stochastic growth algorithms, which are based on the uration is the number of steps of the walk, the energy is
n Rosenbluth and Rosenbluth algorithm [3]. If this algo- the number of non-consecutive nearest-neighbor bonds,
o rithm, which kinetically grows configurations, gets en- and the atmosphere is the number of non-occupied sites
c
hanced by cleverly chosen enrichment and pruning steps around the endpoint of the walk. If all the sites around
[
[4], one obtains the pruned and enriched Rosenbluth the endpoint are occupied, then the atmosphere is zero
2 method (PERM), a powerful algorithm for, e.g., simu- and the walk cannot be continued.
v lation of the polymer collapse transition. The basis of the algorithm is the Rosenbluth and
3
We present in this letter a new algorithm, flatPERM, Rosenbluth algorithm, a stochastic growth algorithm
5
which is a combination of these two types of algorithms, in which the configurations of interest are grown from
2
2 i.e. a flat histogram version of the pruned and enriched scratch. The growth is kinetic, which is to say that
1 Rosenbluth method. each growth step is selected at random from all possi-
3 As opposed to earlier work in this direction [5], in ble growth steps. Thus, if there are a possible ways of
0 which an iterative scheme similar to the multicanonical growth,oneselectsoneofthemwithprobabilityp=1/a.
/
t algorithm [6] was used, we utilize the self-tuning capa- Asthisnumbergenerallychangesduringthegrowthpro-
a
bilities of PERM directly. This leads to a considerable cess,differentconfigurationsarethus generatedwithdif-
m
simplification of the algorithm. ferent probabilities, and one needs to resort to reweight-
- While flatPERM includes umbrella sampling ideas, it ing techniques.
d
n isstriclyspeakingnotamulticanonicalmethod,asmulti- It is advantageous to view this algorithm as approxi-
o canonical sampling conventionally describes a particular matecounting,inwhichcasetheweightsoftheconfigura-
c iterativeversionofadaptiveumbrellasampling,inwhich tions serve as estimates of the number of configurations.
v: first an umbrella distribution is obtained iteratively and To understand this point of view, imagine that we were
i thenafinalsimulationisperformedwithafixedumbrella to perform a complete enumeration of the configuration
X
distribution. space. Doing this requires that at each growth step we
r Thisletterisstructuredasfollows: wefirstgiveaped- would have to choose all the possible continuations and
a
agogicalintroductiontoPERM,whichwillthenallowus countthemeachwithequalweight. Ifwenowselectfewer
to introduce flatPERM as a seemingly trivial extension. configurations,wehavetochangetheweightofthesecon-
As an application, we present simulations of interacting figurationsaccordingly,inordertocorrectformissingout
self-avoidingwalks (ISAW) onthe squarelattice andthe on some configurations. Thus, if we choose one growth
simple cubic lattice. We conclude with a description of step out of a possible ones, we replace a configurations
further applications. with equal weight by one “representative” configuration
Wewillconsideraratherabstractsettingofconfigura- witha-foldweight. Inthisway,the weightofeachgrown
tions with a certain size n, which is parameterized with configuration is a direct estimate of the total number of
an additional variable m. In general, one can even con- configurations.
2
Let the atmosphere a =a(ϕ ) be the number of dis- Define r as the ratio of weight and estimated number of
n n
tinct ways in which a configuration ϕn of size n can be configurations, r = Wn(i)/Cnest. Then we enrich if r > 1
extended. Then, the weight associated with a configu- and prune if r < 1. Moreover, the actual pruning and
ration of size n is the product of all the atmospheres ak enrichmentstepsarechosensuchthattheweightsareset
encountered while growing this configuration, i.e. as closely as possible to Cest to minimize fluctuations:
n
n−1 • r >1 → enrichment step:
W = Yak . (1) make c = min(⌊r⌋,an) distinct copies, each with
k=0 weight 1W(i) (as in nPERM [8]).
c n
After having started S growth chains, an estimator Cest
n • r <1 → pruning step:
for the total number of configurations C (the “infinite-
n continue growing with probability r and weight
temperature” partition function, in which all configura- Cest (i.e. prune with probability 1−r).
tions appear with equal weight 1) is given by the mean n
over all generated samples ϕ(i) of size n with respective Note that we perform pruning and enrichment after the
n
weights W(i), i.e. configurationhasbeenincludedinthecalculationofCnest
n
and is used for weights during the next growthstep. Ini-
1 tially, the estimates Cest canof coursebe grosslywrong,
Cnest =hWin = S XWn(i) . (2) as the algorithm knowns nothing about the system it is
i
simulating. However, even if initially “wrong” estimates
Here,themeanistakenwithrespecttothetotalnumber are used for pruning and enrichment, in all applications
ofgrowthchainsS,andnotthenumberofconfigurations considered the algorithm can be seen to converge to the
which actually reach size n. Configurations which get true values. It is, in a sense, self-tuning.
trappedbeforethey reachsizen appearinthis sumwith At this point it is now straight-forward to change
weight zero. PERM to a thermal ensemble with inverse tempera-
The problem with the Rosenbluth and Rosenbluth al- ture β = 1/kBT and energy E = ǫm by multiply-
gorithm is two-fold. Firstly, the weights can vary widely ing the weight with the appropriate Boltzmann factor
in magnitude,so that the mean canget to be dominated exp(−βE), which leads to an estimate of the partition
by very few samples with very large weight. Secondly, function Zn(β) of the form
regularly occurring trapping events, i.e. generation of
Zest(β)=hW exp(−βE)i . (3)
configurations with zero atmosphere, can lead to expo- n n
nential “attrition”, i.e. exponentially strong suppression
The pruning and enrichment proceduresare changedac-
of configurations of large sizes.
cordingly, replacing W by Wexp(−βE) and Cest by
To overcome both of these problems, enrichment and n
Zest(β), and using r = W(i)exp(−βE(i))/Zest(β). (It
pruningstepshavebeenaddedtothisalgorithm,leading n n m n
is in this setting that PERM is usually described.)
totheprunedandenrichedRosenbluthmethod(PERM)
Alternatively, however, it is also possible to consider
[4]. The basic idea is that one wishes to suppress fluc-
microcanonical estimators for the total number C of
tuations in the weights W(i), as these should on average n,m
n configurations of size n with energy m (i.e. the “density
be equal to C . Therefore, if the weight is too large,
n ofstates”). AnappropriateestimatorCest is thengiven
one enriches, i.e. one makes copies of the configuration n,m
and reduces the weight accordingly. On the other hand, by the mean over all generated samples ϕn(i,)m of size n
if the weight is too small, one prunes probabilistically, and energy m with respective weights Wn(i,m) , i.e.
and, in case the pruning is unsuccessful, continues grow-
1
ainngdwthitehrefaoprpertohperieaxtpelryesisniocnrea(2s)edfowretihgehte.stiNmoatteetCheastt,Sis, Cnes,mt =hWin,m = S XWn(i,m) . (4)
n i
not changed by either enriching or pruning steps.
We need to specify enrichment and pruning criteria Again,themeanistakenwithrespecttothetotalnumber
as well as the actual enrichment and pruning processes. ofgrowthchainsS,andnotthenumberofconfigurations
WhiletheideaofPERMitselfisstraightforward,thereis Sn,m which actually reach a configuration of size n and
now a lot of freedom in the precise choice of the pruning energy m. The pruning and enrichment procedures are
and the enrichment steps. The “art” of making PERM alsochangedaccordingly,replacingCnbyCn,mandusing
performefficientlyisbasedtoalargeextentonasuitable r =W(i) /Cest .
n,m n,m
choice of these steps. Distilling our own experience with Atthis pointit is worthpointing outthat the pruning
PERM,wepresentherewhatwebelievetobeanefficient, and enrichment criterion for PERM leads to a roughly
and, most importantly, parameter free version. constantnumber ofsamples being generatedateachsize
IncontrasttootherexpositionsofPERM(e.g. [7]),we nforPERM.Infact,onecanmotivatethegivenpruning
propose to prune or enrich constantly, to enable larger and enrichment criteria by stipulating that one wishes
exploration of the configuration space (the motivation to have a roughly constant number of samples, as this
will become clear once flatPERM is introduced below). leads to the algorithm performing an unbiased random
3
walk in the configuration size. Correspondingly, in the to overflow problems. This can easily be overcome by
flat-energy version the algorithm performs an unbiased initially restricting the maximal size of grownconfigura-
randomwalkinbothsizeandenergyoftheconfiguration, tions, e.g. by limiting the maximal size by the number
and we obtain a roughly constantnumber of samples for of started growth chains,
each size n and energy m.
It is because of the fact that the number of samples is n<cS . (6)
roughlyconstantin eachhistogramentry, thatthis algo-
The relevant number of growth chains at size n is thus
rithmcanbeseenasa“flathistogram”algorithm,which
reduced to S −n/c. For ISAW we find that a choice of
we consequently call flatPERM.In hindsight in becomes
delay of c≈10 is sufficient.
clear that PERM itself can be seen as a flat histogram
algorithm, at it creates a roughly flat histogram in size
n. Recognizing this led us to the formulation of this al-
gorithm in the first place. log10(Cnm)
Atthispointwereturntoourdiscussionofthepruning 450
400
350
and enrichment strategies. For PERM it may be more 300
250
advantageoustoallowforahighdiffusivityalongthesize 200
150
n. This is done by minimizing pruning and enrichment, 15000
0
i.e. atthecostofallowinglargerweightfluctuations. For
1000
flatPERMontheotherhandweneedtoachievediffusive 800
0 600
behavioralsowithrespecttotheenergyvariablem. Pre- 0.2 0.4 400 n
0.6 200
cisely this is achieved by allowing for much enrichment m/n 0.8 1 0
(andthusnecessarilyalsopruning),aseachsetofconfig-
urations enriched at size n will contribute to a range of
log10(Cnm)
different energies m at size n+1.
700
EventhoughweviewflatPERMasaflathistogramal- 600
500
gorithm,wehavenotyetexplicitlyconditionedtheprun- 400
300
ing and enrichment with respect to the local number of 200
100
samples Sn,m created at each size n and energy m. This 0
canbedonebymultiplyingrbyS/S ,asthisenhances 1000
n,m
800
epnlerOsi,cnhaemngdeennvteicirfeaSlvnpe,rrmsoab<.leSm,iw.ei.thwPheEnRtMheriseathreattoeonrfeicwhmsaemnt- 0 0.2 0.4 0.6 0.m8 =n1 1.2 1.4 1.6 1.8 200 400 60n0
generates large correlation between samples, so that it
FIG. 1: Logarithm of the number of configurations Cn,m
would be useful to replace the number S of samples
n,m versus internal energy m/n and length n for ISAW on the
by a number Seff of effectively independent samples,
n,m squarelattice (above) and simple cubic lattice (below).
thereby taking account of some autocorrelation time.
Thiscanbedoneheuristicallybyconsideringtheaverage
AveragesofobservablesQdefinedonthe setofconfig-
number of independent growth steps in a configuration.
urations can now be obtained by storing weighted sums
If the last enrichment has occurred at size n , then a
enr of these observables, from which one obtains
configurationofsizenandenergym hasn =n−n
ind enr
independent steps. The frequency of independent steps hQWi Q(i) W(i)
in a configuration is w(i) =n /n. We thus use Qest = n,m = Pi n,m n,m . (7)
n,m ind n,m hWin,m PiWn(i,m)
S W(i)
r= Snef,mf Cnens,,mtm with Snef,mf =Xi wn(i,)m (5) Tinhsteasenccea,nthteheenxpbeectuasteidonfovralsuuebosefqQueinnttehveacluaantoinonicsa.lFenor-
semble at a given temperature β can now be computed
for the pruning and enrichment criterion in our algo- via
rithm.
Now that one has a flat histogram method of PERM Qest(β)= PmQens,mt Cnes,mt exp(−βEm) . (8)
for the whole range of energies, one can easily modify n Cest exp(−βE )
Pm n,m m
the algorithmfurther to sample only a selected range by
dividing r by a “profile shape” f , leading to pruning We have implemented this algorithm for interacting
n,m
whenever f is close to zero. This is advantageous if self-avoiding walks on the square and simple cubic lat-
n,m
oneonlywantstoexplorearestrictedregionofparameter tice. In both two and three dimensions, we have sim-
space. ulated walks up to length n = 1024. Figures 1 and 2
Eventhoughthealgorithmdescribedthusfarisfreeof clearly show the strengthof the method. Figure 1 shows
parameters, there is a technical problem due to the fact the number of configurationsC , which vary over sev-
n,m
that initial weights are grossly wrong, which can lead eralhundredordersofmagnitude. Thisrangewouldhave
4
Snm
1.2
1e+06
1.1
1.0
100000 0.9
n1(cid:27)2(m)0.8
0.7
10000
0.6
1000
800 0.5
0 600 0.4
0.2 0.4 0.6 200 400 n 0.40 0.50 0(cid:12).60 0.70 0.80
m/n 0.8 1 0 5.0
4.0
Snm
1e+07 n1(cid:27)2(m)3.0
2.0
1e+06 1.0
100000 0.0
0.25 0.30 0.35 0.40 0.45
10000 (cid:12)
800 1000 FIG.3: Normalizedfluctuationsσ2(m)/nversusinversetem-
0 0.2 0.4 400 60n0 perature β =1/kBT for ISAW on the square lattice (above)
0.6 200 and the simple cubic lattice (below) at lengths 64, 128, 256,
Snefmf m/n 0.8 1 0 512, and 1024. The curvesfor larger lengths are more highly
peaked. Theverticallinesdenotetheexpectedtransitiontem-
1e+06
peratureat infinite length.
100000
10000
1000
100 shown for both systems in Figure 3.
1000
We have also tested our algorithm on the HP model
800
0 600 [10], which is a toy model for proteins. It consists of
0.2 0.4 400 n
m/n0.6 0.8 1 0 200 a self-avoiding walk with two types of monomers along
the sites visited by the walk, which are either hydropho-
bic (type H) or polar (type P). One considers monomer-
FIG.2: Numberofgeneratedsamplesversusinternalenergy
specificinteractions,mimickingtheinteractionwithapo-
m/nandlengthnforISAWonsquarelattice. Thetopgraph
showsthenumberofsamplesSn,m forasimulationinwhicha lar solvent such as water. The interaction strengths are
flat histogram for Sn,m was created, whereas themiddle and thus chosen such that HH-contacts are favored, e.g. by
bottom graphs showthenumberof samples Sn,m andtheef- choosing ǫHH = −1 and ǫHP = ǫPP = 0. The central
fectivenumberofsamplesSeff,respectively,forasimulation questionistodeterminethedensityofstates(andtofind
n,m
in which a flat histogram for Seff was created. thegroundstatewithlowestenergy)foragivensequence
n,m
of monomers. For various sequences taken from the lit-
eraturewe haveconfirmedpreviousdensity ofstates cal-
been inaccessible during one simulation with any canon- culations and obtained identical ground state energies.
icalmethod. As anaside,we note thatwe hadto rescale The sequences we considered had n = 58 steps (3 di-
weights during the run to avoid overflow [9]. Addition- mensions, Emin = −44) and n = 85 steps (2 dimen-
ally,wechoseadelayofc=10tostabilizethealgorithm. sions, Emin = −53) from [7], and n = 80 steps (3 di-
As can clearly be seen in Figure 2, the number of mensions, Emin = −98) from [11]. We studied also a
samples generated is roughly constant, irrespective of particularly difficult sequence with n = 88 steps (3 di-
whether one flattens withrespectto Sn,m (topgraph)or mensions, Emin =−73) from [12], but the lowest energy
Seff (middle and bottom graphs). Only results for the we obtained was E = −69. However, in this case no
n,m
squarelatticeareshown,aswefindthesamebehaviorfor lower energy has been found with any other PERM im-
the simple cubic lattice. Using Seff for the flatness cri- plementation either, see [7].
n,m
terion leads moreover to an increased sampling of walks To summarize, we have presented a new algorithm,
withveryfewandverymanyinteractions,thusovercom- flatPERM,whichisaflathistogramversionofastochas-
ing the usual difficulty of obtaining configurations with tic growthalgorithm. This algorithmcan in principle be
alargenumber ofinteractions. (The largestenergystate applied to any statistical mechanical system for which
gets repeatedly sampled in simulations in both dimen- configurations can be grown in a well-defined manner.
sions.) Next to the presented applications of linear polymers,
Once the simulations have been completed, ther- this algorithm can be applied to more complicated sys-
modynamic quantities of interest can easily be com- tems, such as lattice models of branched polymers [13].
puted. As an example, the specific heat curves C = Extensions to models with two energy parameters are in
n
k (βǫ)2σ2(m)/n near the transition temperature are preparation, e.g. to the problem of adsorbing interact-
B
5
ing polymers [14] or an extended Domb-Joyce model of Owczarek for helpful discussions, Peter Grassberger for
polymer collapse [15]. pointingoutanincorrectformula,therefereesforhelpful
The authors thank Andrew Rechnitzer and Aleks comments, and the DFG for financial support.
[1] F. Wang and D. P. Landau, Phys. Rev. Lett. 86 2050 (20) in B. A. Berg, Comp. Phys. Commun. 153 397
(2001) (2003)
[2] G.M.TorrieandJ.P.Valleau,J.Comput.Phys.23187 [10] K. A. Dill, Biochemistry 24 1501 (1985)
(1977) [11] H. Frauenkron and U. Bastolla and E. Gerstner and P.
[3] M.N.RosenbluthandA.W.Rosenbluth,J.Chem.Phys. Grassberger and W. Nadler, Phys. Rev. Lett. 80 3149
23 356 (1955) (1998)
[4] P.Grassberger, Phys. Rev E 56 3682 (1997) [12] T.C. Beutler and K.A. Dill, Protein Sci. 5147 (1996)
[5] M.BachmannandW.Janke,Phys.Rev.Lett.91208105 [13] A. Rechnitzer and E. Janse van Rensburg and T. Prell-
(2003) berg, in preparation
[6] B.A.BergandT.Neuhaus,Phys.Lett.B267249(1991) [14] A.L.Owczarekand A.RechnitzerandJ.Krawczykand
[7] H.-P. Hsu and V. Mehra and W. Nadler and P. Grass- T. Prellberg, in preparation
berger, J. Chem. Phys.118 444 (2003) [15] J. Krawczyk and T. Prellberg and A. Rechnitzerand A.
[8] H.-P. Hsu and V. Mehra and W. Nadler and P. Grass- L. Owczarek, in preparation
berger, Phys. Rev.E 68 21113 (2003)
[9] An alternative is to use logarithmic coding, see e.g. eqn.
