Efficient configurational-bias Monte-Carlo simulations of chain molecules with `swarms' of trial configurations
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] O c t Efficient configurational-bias Monte-Carlo simulations of chain moleculeswith ‘swarms’ of trial configurations
Niels Boon
Division of Physical Chemistry, Department of Chemistry, Lund University, SE-22100 Lund, Sweden
Proposed here is a dynamic Monte-Carlo algorithm that is efficient in simulating dense systems oflong flexible chain molecules. It expands on the configurational-bias Monte-Carlo method throughthe simultaneous generation of a large set of trial configurations. This process is directed by at-tempting to terminate unfinished chains with a low statistical weight, and replacing these chainswith clones (enrichments) of stronger chains. The efficiency of the resulting method is explored bysimulating dense polymer brushes. A gain in efficiency of at least three orders of magnitude is ob-served with respect to the configurational-bias approach, and almost one order of magnitude withrespect to recoil-growth Monte-Carlo. Furthermore, the inclusion of ‘waste recycling’ is observedto be a powerful method for extracting meaningful statistics from the discarded configurations.
I. INTRODUCTION
Polymer chains are challenging to model with com-puter simulations due to the vast number of possiblechain configurations that quickly increases with thenumber of monomers. Additionally, the explorationof phase space is hindered by chain entanglements.While Molecular-Dynamics (MD) methods can beapplied to follow the natural time evolution of suchsystems, Monte-Carlo(MC) approaches enable the intro-duction of unphysical ‘moves’ between configurations,which may increase the rate of generating uncorrelatedconfigurations[1, 2]. The MC method comprises twocategories of approaches to sampling phase space. A static
MC algorithm explores phase space throughsuccessive generation of uncorrelated configurationsfrom scratch[3, 4].
Dynamic
MC algorithms, on theother hand, will attempt to generate a new configurationbased on the existing state. This characterization refersto the creation of a Markov chain of configurations [5] incontrast to the physical dynamics of the system.The ‘simple sampling’ of polymer configurations by theconsecutive adding of monomers at random orientationsto form a chain is an example of a static approach. Sim-ple sampling is, however, not efficient because most gen-erated configurations will have a vanishing Boltzmannweight due to overlaps between monomers. This obstaclecan be partially avoided by considering multiple ‘probe’positions for every monomer that is added to the grow-ing chain, and select from those with a probability pro-portional to their (resulting) Boltzmann weight. Thatapproach, which largely avoids monomer positions thatlead to overlaps, is known as the Rosenbluth-Rosenbluth(RR) method. It requires keeping track of a statisticalweight W i to remove the introduced biasing in this selec-tion process of a polymer-chain configuration. Equilib-rium properties h A i follow from summing over all gener-ated configurations i h A i = P i W i A P i W i . (1) It was found, however, that for longer chains the RRmethod yields a very wide spread in the weights W i , suchthat only a few configurations dominate the weighedaverage in Equation (1). The simulation, therefore,will spend most of its time on configurations that donot contribute significantly to this weighted average[6].The pruned-enriched Rosenbluth method (PERM) wastailored to address this inefficiency[7]. This method con-trols the generation (growth) of the chain configurationsby attempting to terminate(prune) unfinished chainswith a below-average weight. On the other hand, thosewith an above-average weight are cloned(enriched) suchthat two copies of the incomplete chain continue to growwith half of the original weight each[8]. Pruning andenrichment maintains the correct statistics while morecomputational time is spent on polymer configurationsof significant weight. This leads to a much more efficientsampling of phase space. The PERM approach enabledthe sampling of very long chain molecules[7, 9, 10] andprotein structures using minimalist models[11].Static algorithms such as RR and its extension PERMare not favorable for systems with multiple polymers asthey involve finding a new configuration from scratch for all chains in the system at the same time. This approachquickly becomes ineffective as the number of chainsin the system is increased beyond one. It is, however,possible to render the RR method dynamic and theresulting approach is known as the configurational-biasMonte-Carlo (CBMC) method [3, 12–15]. Each stepin the CBMC algorithm applies the RR method togenerate a new configuration for a randomly chosenchain in the system. The existing configuration is thenalso ‘retraced’ (re-weighed) such that the acceptanceprobability of the new configuration can be determinedby comparing the weights of the new and the existingconfiguration. The CBMC method has proven to bevery effective for finding properties of chain moleculessuch as alkanes[16–22], and has been extended to widearray of other systems[23–28].This work aims at incorporating pruning and enrich-ment into the CBMC method. The approach is, nev-ertheless, quite different from earlier attempts such asDPERM[29], for which only a marginal increase in ef-ficiency w.r.t. CBMC was reported. An essential ele-ment is the simultaneous generation (growth) of a fixednumber of ‘candidate’ configurations, while at the sametime the reference configuration is retraced. This syn-chronized growth (and retracing) of candidate chains en-ables comparison of the weights of unfinished configura-tions. Pruning or enrichment can, therefore, be initiatedbased on the relative weights of the generated chains.Moreover, the size of the set can easily be kept constantduring growth: the algorithm will attempt to terminate(prune) chains that have a much lower effective statisti-cal weight than the other candidate chains in the set andreplace those with clones (enrichments) of ones with thelargest weight. The computational effort in each MC stepis, therefore, fixed by the selected number of candidatechains in the set. The performance of this approach willbe analyzed through the simulation of polymer brushes. II. ALGORITHM
Every MC step involves randomly selecting and re-moving a polymer chain from the system. The algo-rithm attempts to update the configuration of this chain,and inserts it back into the system at the end of thestep. The CBMC approach to the former is the gener-ation of a ‘probe’ configuration c n with weight W n , asdescribed in algorithm A below. Also, a weight W e forthe old configuration c e is obtained, as algorithm B de-scribes. The acceptance rate of c n as the new configu-ration is determined by a Metropolis-form probability[5] p accep = min(1 , W n /W e ). Starting at monomer ℓ = 1, thegrowth of a probe configuration c = ( r , . . . , r L ) proceedsas follows.A1 Construct a set { ζ , . . . , ζ k } consisting of k trial po-sitions to insert monomer ℓ . These positions areselected with a relative probability p ( ζ i ) = ν ℓ exp( − βu bond ℓ ( ζ i )) , (2)where ν ℓ is a normalization constant, β the usualthermodynamic beta, and u bond ℓ ( ζ ) the bonding en-ergy for monomer ℓ , which is defined w.r.t. theposition and orientation of the previous monomer.Note that the first monomer that is inserted ( ℓ =1) has a vanishing bonding energy, although forgrafted polymers u bond1 ( ζ ) is the bonding energyto the grafting surface. The Rosenbluth factor as-sociated with the set of trial positions is w ℓ = P kj =1 exp[ − βu ex ℓ ( ζ j )], where u ex ℓ ( ζ ) is the potentialenergy of monomer ℓ which takes account of inter-actions with other monomers or fields. If w ℓ > r ℓ is selected from the set with prob-ability p j = exp[ − βu ex ℓ ( ζ j )] /w ℓ and the monomeris added to the chain. A2 If a monomer was added, i.e. w ℓ >
0, then step A1is repeated until the chain is complete ( ℓ = L ).The algorithm for re-tracing an existing configuration,starting from ℓ = 1, is related to A.B1 Choose ζ = r ℓ and generate the remaining set { ζ , . . . , ζ k } of k − ℓ using Equation (2). The Rosenbluthfactor for this monomer is also given by w ℓ = P kj =1 exp[ − βu ex ℓ ( ζ j )]. The monomer is added tothe chain at position r ℓ .B2 Step B1 is repeated for every next monomer, untilthe last monomer has been weighed ( ℓ = L ).Using either algorithm, the resulting weight and en-ergy of a chain are calculated as W = Q Lℓ =1 w ℓ and U = P Lℓ =1 (cid:0) u bond ℓ ( r ℓ ) + u ex ℓ ( r ℓ ) (cid:1) , respectively. Zeroweight results for terminated chains.The proposed algorithm takes a different approachfrom CBMC through the simultaneous construction of M − M candidate configu-rations. There are no interactions between the monomersof different candidate configurations and separate Verletlists may be used for each of them(as well as one for therest of the system). Starting from monomer ℓ = 1, thegrowth of this set proceeds by repeating the followingsteps L times.C1 Monomer ℓ is added to the M − A
1. Monomer ℓ of the referenceconfiguration is weighed using algorithm B
1. The(Rosenbluth) weight of each unfinished chain W ℓ = Q ℓℓ ′ =1 w ℓ ′ and its effective weight W ∗ ℓ = 2 γ W ℓ iscalculated. Here, γ = y p − y e , with y p the number ofprune attempts, and y e the number of times cloning(enrichment) occurred(see below).C2 Chains with W ∗ ℓ = 0 are terminated. The averageeffective weight ˜ W ∗ ℓ of the other chains is calcu-lated. Chains with W ∗ ℓ < ˜ W ∗ ℓ / / γ raised by 1 otherwise. The re-traced chain, however, cannot be terminated andhas γ raised by 1 if it is marked for pruning.C3 The chain with the largest effective chain weight W ∗ max ℓ is selected. If there is more than one chainwith this effective weight then one of these is se-lected randomly . A terminated chain is now re-placed by a clone of this selected chain. Both chains To avoid comparison of floating-point numbers for equality onemay select randomly from chains with W ∗ ℓ > (1 − ǫ ) W ∗ max ℓ in-stead. Here, ǫ is a small number(e.g. 10 − ). obtain half of the effective weight W ∗ ℓ since γ is nowlowered by 1. If the retraced chain is cloned thenthis clone proceeds growing as a probe chain. Thisstep is repeated until all terminated chains havebeen replaced.Algorithm C produces a set of M configurations c m ,1 ≤ m ≤ M and corresponding (effective) weights W ∗ ( m ) ≡ W ∗ L ( m ). One of these candidates is selectedwith a relative probability P accep ( m ) = W ∗ ( m ) P Mm =1 W ∗ ( m ) . (3)This configuration updates the old polymer configurationin the system. III. JUSTIFICATION OF METHOD
For convenience, a path x will be defined as a chainconfiguration as well as the remaining trial positions thatwere generated by algorithm A or B, x = (cid:16) ( r , . . . , r L ′ ) , ( ζ k − , . . . , ζ k − L ′ ) (cid:17) , (4)where c ( x ) := ( r , . . . , r L ′ ) is the chain configurationand ζ k − ℓ := { ζ ℓ, , . . . , ζ ℓ,k } are the other k − ζ ℓ,j for each monomer ℓ . Terminated chains arecharacterized by L ′ < L and W ( x ) = 0, while successfulchains correspond to paths with L ′ = L . Every x has aunique weight W ( x ).Consider the probability of generating a path x witheither algorithm A or B. Algorithm A generates a set of k trial positions { ζ ℓ, , . . . , ζ ℓ,k } for every monomer ℓ , andselects r ℓ from this set with a probability e − βu ( r ℓ ) /w ℓ ( x ).Algorithm B, on the other hand, fixes ζ ℓ, = r ℓ andonly needs to generate the k − { ζ ℓ, , . . . , ζ ℓ,k } associated with x . Whilst in al-gorithm A the position r ℓ is selected with a probabil-ity exp[ − βu ex ℓ ( r ℓ )] /w ℓ , in algorithm B this probabilityiremarks s always 1. Therefore, the overall probability p A ( x ) of generating a path x with algorithm A is relatedto the probability p B ( x ) of generating x with algorithmB as p A ( x ) p B ( x ) = L Y ℓ =1 (cid:18) kν ℓ e − βu bond ℓ ( r ℓ ) e − βu ex ℓ ( r ℓ ) w ℓ ( x ) (cid:19) = ν e − βU ( c ( x )) W ( x ) , (5)where Equation (2) was used and the ν ≡ Q ℓ kν ℓ is aconstant. Note that the probabilities of generating theunselected trial positions have cancelled each other inthis equation.Algorithm C generates a set s of paths during ev-ery MC step, as sketched in Figure 1. Despite thedifferent coloring (red/blue) in the latter figure of the FIG. 1 Sketch of a set s of paths explored in one Monte-Carlo step. The generated configurations are colored: probeconfigurations are blue and the retraced configuration is red.The number of trial position k per monomer is set to threehere, and the k − M is fourat any time during the generation of the set. probe configurations and the retraced configuration, s itself does not hold information on which of thepaths corresponds to the retraced configuration. Thesame set s can result from any retraced configurations c ( x ) as long as x is a successful chain contained in s .Nevertheless, the probability P C ( s, x ) of generating s with algorithm C depends on the path x to which theretracing algorithm B was applied. The rationale belowwill therefore be focused on comparing P C ( s, x ) betweendifferent choices of x . For convenience, one may firstrelate P C ( s, x ) to the probability P C ∗ ( s ) of generating s with an (hypothetical) algorithm C* that does notretrace any existing configuration but applies algorithmA to all the paths. Consider selecting one of the M (initially empty) paths and following its generation,starting from ℓ = 1. Given that s is generated, theprobability that this particular path becomes x is1 / ( M · y e ( x,s ) ), with y e ( x, s ) the number of times x was cloned during its growth. The latter takes accountfor the fact that each cloning event branches off a newpath that has an equal probability of becoming x . Thisdefines P C ∗ ( s, x ) = P C ∗ ( s ) / ( M · y e ( x,s ) ), which is theprobability of generating s while at the same time theselected path yields x .In comparison with algorithm C*, algorithm C indeedselects one the M initial paths and follows(directs) thegrowth of this path through the application of algorithmB. The probability that this path results in x thereforeincreases by an additional factor ( p B ( x ) /p A ( x )) · y p ( x,s ) .Here, y p ( x, s ) is the number of times x was marked forpruning and accounts for the fact that a retraced chaincannot be terminated. Thus, the probability P C ( s, x ) ofgenerating s while having x ∈ s generated by retracingcan be expressed as P C ( s, x ) = P C ∗ ( s, x )( p B ( x ) /p A ( x )) · y p ( x,s ) , which yields P C ( s, x ) = P C ∗ ( s ) (cid:20) M γ ( x,s ) p B ( x ) p A ( x ) (cid:21) , (6)recalling that γ ( x, s ) = y p ( x, s ) − y e ( x, s ). Any move froman old polymer configuration c a to a new configuration c b must proceed through the generation of a set s thatcontains c a as well as c b , i.e. x a , x b ∈ s , where c ( x a ) = c a and c ( x b ) = c b . By defining P sc a → c b as the probability ofmoving from c a to c b via the generation of s , one finds P sc a → c b P sc b → c a = P C ( s, x a ) · P accep ( x b , s ) P C ( s, x b ) · P accep ( x a , s ) , (7)where P accep ( x, s ) is the probability of updating the poly-mer configuration to c ( x ), given the set s . By combiningEqs. (6) and (5) it is possible to rewrite Equation (7) as P sc a → c b P sc b → c a = e β ( U ( c a ) − U ( c b )) W ( x a )2 γ ( x a ,s ) W ( x b )2 γ ( x b ,s ) P accep ( x b , s ) P accep ( x a , s ) . (8)and by using the selection criterion for a candidate con-figuration 3 the existence of (a variation of) superdetailedbalance[15] is confirmed. Detailed balance is, therefore,satisfied globally as the total transition rate between c a to c b follows from all sets s that obey this balance. IV. SIMULATIONS
The algorithm that is introduced here, which will bereferred to as ‘swarm’ confrontational-bias Monte-Carlo(SCBMC), is tested on a system of polymer brushes.Polymer brushes can be used as lubricants, adhesives,or to stabilize colloidal suspensions and, consequently,have been extensively studied by simulations andtheory[30–36]. They form an interesting model systemto test the algorithm since the simulation of dense,long brushes is computationally demanding. For thesake of generality, a simple freely-jointed chain modelis considered here. The simulation box has horizontaldimensions H × H and is defined with periodic boundaryconditions in this plane. Monomers are restricted inthe vertical direction to z ≥
0, such that z = 0 definesthe grafting surface to which N = 60 polymers aregrafted. Each chain is composed of L hard sphereswith diameter d = 1, yielding a (dimensionless) graftingdensity of σ = N d /H [30]. For the initial configu-ration a random grafting of fully stretched chains is used.All simulations were run on a single core of a IntelI7-6700 CPU in a desktop computer. Calculated densityprofiles have been checked for numerical accuracy bycomparing to earlier work on similar systems[30]. Thedata in Figure 2, furthermore, confirms the Alexanderscaling of the obtained density profiles with the graftingdensity[37] for a range in σ and two different polymerlengths. Deviations from a parabolic profile[33, 38] canbe observed at extremely large grafting densities suchas σ = 0 . FIG. 2 Scaled monomer density plotted as a function of thescaled distance from the grafting plane. The blue(open) sym-bols denote data for L = 50 and the green (closed) sym-bols correspond to L = 10. Simulations ran for 10 secondsfor L = 50 and 10 seconds for L = 10. Grafting densities σ ≥ .
25 ran 10 times longer.FIG. 3 Sample chosen from simulations to demonstrate thepruning and cloning process applied on the swarm, showingthe x -coordinate for the set of candidate chains that is gener-ated in a MC step. The red curve corresponds to the retracedchain and the other completed chain configurations are darkerblue. The lighter curves show the terminated chains. This fig-ure is obtained from a simulation with L = 50 and σ = 0 . M was set to 150. Figure 3 corresponds to a set of candidate configura-tions that the SCBMC algorithm generates during oneMC step by plotting the x -coordinate of each monomeralong the chains. The red curve shows the monomerpositions of the retraced chain, while the other candidateconfigurations are drawn in darker blue. Although M was set to 150 in this simulation, the number of chain-termination events during the chain growth was muchlarger ( ≈ M enables the algorithm toexplore ‘dead-end’, as every terminated chain can easilybe replaced by a clone of a more successful chain. CBMC M = 10 M = 100 M = 1000 L =10, σ =0 .
20 19 . ± . . ± . . ± . . ± . L =50, σ =0 .
05 199 . ± . . ± . . ± . . ± . L =50, σ =0 .
20 602 ± ± . ± . ± TABLE I The effect of varying the number of candidatechains M on the calculated mean squared end-to-end distance h R z i and its standard error of the mean, σ m . The number oftrial directions per monomer, k , was set to 10 here.FIG. 4 Comparison of the standard error of the mean, SE m between the CBMC method and the SCBMC method for anextended number of candidate chains M w.r.t. to Table 1.All indicated confidence levels in SE m are estimated by boot-strapping. The efficiency of the method in exploring phasespace is analyzed by running approximately 10 -secondsimulations during which the mean-squared height ofthe end of the grafted chains h R z i was measured over10 equal time intervals. The first of those is discardedfor equilibration. For the shorter chain ( L = 10) thesimulations only ran for approximately 10 seconds intotal. Table 1 shows results for 3 different combinationsof L and σ . Most runs converged on the value of h R z i , but a discrepancy is observed for the CBMCsimulations of the longest and densest chains. Aninspection of the final configuration revealed that theCBMC simulation did not completely move away fromthe initially stretched state within the given time. Thisrendered the value of h R z i too large and underesti-mated the standard error due to correlations betweensamples. This also occurred for SCBMC simulationsusing M = 2 but was absent for other tested values of M .To get a better idea of the relative efficiency of differ-ent runs, Figure 4 plots the error of the mean in h R z i for a wider range in M than Table 1. For longer chainsthe results show that more candidates generally leads tobetter statistics, with an optimum value of M ≈
50. This is most pronounced for σ = 0 .
2. It can also beobserved that simulations of short chains do not benefitfrom SCBMC here, which may be explained by the factthat within CBMC the success rate of updating configu-rations is already high for L = 10. FIG. 5 Full chain updates during the course of 10 s sim-ulation time. The simulated brush has length L = 50 andgrafting density σ = 0 .
2. The curves result from the SCBMCalgorithm for given values of M , while varying k along thehorizontal axis. The dashed curve shows data from the recoil-growth algorithm, using a recoil length L R = 24. Based on Figure 3 it is expected that monomerslower in the chain are updated least frequently and willdetermine the longest correlation times. Consequently,the time it takes to update every monomer positionin the system may be used as an indicator for theefficiency of the algorithm. Figure 5 plots the numberof full updates during the run time of simulations withvarious algorithm parameters M and k , using L = 50and σ = 0 .
2. These results indicate the existence of anoptimal choice for the number of trial directions k forevery choice of M . The CBMC approach did not pro-duce any updates for any value of k , a full update fromthe stretched initial configuration did even not occur forsimulations that ran 20 times longer, and so it was notincluded in this plot. For the data shown, simulationsusing M ≈
100 and k ≈ k = 1 with a largenumber of candidate chains ( M = 400).The efficiency of the SCBMC method may be under-stood from the possibility of cloning chains if some ofthem encounter ‘dead-ends’. Effectively, this redirectsthese chains into a direction that could be more success-ful. Since the recoil-growth (RG) algorithm [40, 41] isbased on similar ideas it is interesting to compare the ef-ficiencies of both methods. The RG method extends onthe CBMC method by enabling the probe chain to retractup to a maximum number of L R monomers once it meetsa dead-end. This increases the chance that a successfulchain configuration is found. An implementation of thismethod was checked for consistency with other methods.Full chain updates were attempted with a probability1 /
2, and partial chain updates otherwise. Note that par-tial chain updates cannot trigger counts of full updatesdirectly, yet they were observed to be essential for equi-libration of the system in test runs. Since full updates ofall monomer positions did not occur frequently, simula-tions were ran for 5 times longer for the RG method. Theparameters L R and k were then optimized for efficiency.No full update of all monomers occurred for any k exceptfor k = 4, for which optimal performance was observedby choosing L R ≈
24. The SCBMC method, therefore,seems to be less sensitive to the choice of various simu-lation parameters than the RG method. Although RGperforms well compared to CBMC, the data in Figure 5also indicates that SCBMC outperforms RG by almostan order of magnitude in efficiency.
V. WASTE RECYCLING
FIG. 6 Probability ρ e ( z ) of finding the end monomer at ele-vation z from the grafting surface. While the colored curvesare obtained without waste-recycling Monte Carlo (WRMC),the curves consisting of (small) black dots result from Equa-tion (10). Simulations were performed for σ = 0 .
2, using M = 150, and k = 10. For short ( L = 10), longer ( L = 50),and longest chains ( L = 100) these ran for 10 , 10 , and 10 seconds respectively. Although the generation of a large number of candi-date configurations is essential for optimal performanceof the SCBMC algorithm, only one of these is recordedinto the Markov chain of the system as the simulationproceeds. Recognizing the value of the discarded con-figurations is the central idea in waste-recycling MonteCarlo (WRMC), which re-uses these configurations forcalculations of Boltzmann-weighted averages of observ-ables [4, 42]. Because of detailed balance, one may record observables A n in succession to all Monte-Carlo steps s n that attempt to update polymer n , i.e. h A n i ← P s n A n P s n . (9)It was shown in Ref. [4] that improved statistics canbe obtained if a weighted average over the accepted aswell as the rejected configurations is calculated. The rel-ative weights herein are set by the acceptance probabili-ties of the candidate configurations, so equation (9) canbe rewritten as h A n i ← P s n ! X s n P Mm =1 W ∗ ( m ) A n,m P Mm =1 W ∗ ( m ) , (10)where W ∗ ( m ) and A n,m are the effective weights andthe value of observable A n , respectively, correspondingto candidate configuration m . Equation (10) is usedto sample the probability h ρ e ( z ) i = P n h ρ en ( z ) i /N offinding the free end of a chain at a height z from thegrafting plane. Figure 6 plots results for brushes with σ = 0 . L . It canbe observed that WRMC yields a strong reductionof the statistical noise in the determined profiles.The effectiveness of this combined approach is bestillustrated by the acquired level of detail for the case L = 100, which is quite extraordinary for brushes withthis density and length. WRMC may be particularlyuseful for determining the properties of the terminalmonomer, as most of the generated candidate chainspossess a unique configuration of this particular beadin the chain. On the other hand, properties relatedto monomers earlier in the chain may benefit less astheir position is shared between multiple candidatechains. Nevertheless, waste recycling is basically ‘free’ ofcomputational effort to implement, so there is good rea-son to use it in combination with the introduced method. VI. CONCLUSION
With respect to the physics that is discussed the re-sults in this work agree with theoretical predictions forlong and dense polymer brushes. More interesting isthe efficiency of the introduced algorithm in simulatingsuch systems. The simultaneous generation of a largeset (swarm) of candidate configurations in each step ofthe dynamic MC algorithm can lead to an efficient sam-pling of phase space. An essential aspect of the intro-duced algorithm is the rigorous application of pruningand cloning(enrichment) during the generation of the setof candidate chains, which are balanced to keep the sizeof the swarm constant. It also ensures that all remain-ing configurations, including the reference configuration,have a comparable probability of being accepted. Theoptimal swarm size depends on the complexity of thesystem. This algorithm obtains accurate statistics fromsystems that do not seem to thermalize with the CBMCapproach, and a strong increase of efficiency with re-spect to the recoil-growth method is observed. The par-allelized chain generation can be distributed over mul-tiple processing elements in future work[43]. Moreover,this method is tailored for waste-recycling Monte Carlotechniques, which enables the calculation of statisticalaverages from all generated chain configurations in theset. The method is not limited to polymer systems only,and may be extended beyond chains to groups of par-ticles in general. For CBMC this has already includedmixtures of large and small particles [24, 25], phase-equilibrium calculations[26], transition path sampling[23]and the structure of water[27] or ionic solutions [28].
Acknowledgments
I would like to thank Daan Frenkel for his suggestion toconsider incorporating waste-recycling Monte Carlo intothe method. I gratefully acknowledge Alex Cumberworthfor providing numerous valuable remarks on multiple ver-sions of the manuscript.
References [1] G. S. Grest and M. Murat,
Monte Carlo and MolecularDynamics Simulations in Polymer Science (Oxford Uni-versity Press, New York, 1995).[2] M. E. J. Newman and G. T. Barkema,
Monte-CarloMethods in Statistical Physics (Clarendon, Oxford,1999).[3] D. Frenkel and B. Smit,
Understanding molecular simu-lation: from algorithms to applications (Academic Press,London, 2002).[4] D. Frenkel, Lect. Notes Phys. , 127 (2006).[5] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth,A. H. Teller, and E. Teller, J. Chem. Phys. , 1087(1953).[6] J. Batoulis and K. Kremer, J. Phys. A. Math. Gen. ,127 (1987).[7] P. Grassberger, Phys. Rev. E , 3682 (1997).[8] F. T. Wall and J. J. Erpenbeck, J. Chem. Phys. , 634(1959).[9] H. Frauenkron and P. Grassberger, J. Chem. Phys. ,22 (1997).[10] G. T. Barkema, U. Bastolla, and P. Grassberger, J. Stat.Phys. , 1311 (1998).[11] U. Bastolla, H. Frauenkron, E. Gerstner, P. Grassberger,and W. Nadler, Phys. Rev. Lett. , 3149 (1998).[12] J. I. Siepmann and D. Frenkel, Mol. Phys. , 59 (1992).[13] J. J. de Pablo, M. Laso, and U. W. Suter, J. Chem.Phys. , 2395 (1992). [14] T. J. H. Vlugt, M. G. Martin, B. Smit, J. I. Siepmann,and R. Krishna, Mol. Phys. , 727 (1998).[15] D. Frenkel, G. C. a. M. Mooij, and B. Smit, J. Phys.Condens. Matter , 3053 (1999).[16] B. Smit, S. Karaborni, and J. I. Siepmann, J. Chem.Phys. , 2126 (1995).[17] M. G. Martin and J. I. Siepmann, J. Phys. Chem. B ,2569 (1998).[18] T. J. H. Vlugt, R. Krishna, and B. Smit, J. Phys. Chem.B , 1102 (1999).[19] H. Wu, Q. Gong, D. H. Olson, and J. Li, Chem. Rev. , 836 (2012).[20] D. Dubbeldam, A. Torres-Knoop, and K. S. Walton,Mol. Simul. , 1253 (2013).[21] R. Krishna and J. M. van Baten, Phys. Chem. Chem.Phys. , 7994 (2013).[22] A. Sepehri, T. D. Loeffler, and B. Chen, J. Chem. TheoryComput. , 074102 (2014).[23] C. Dellago, P. G. Bolhuis, F. S. Csajka, and D. Chandler,J. Chem. Phys. , 1964 (1998).[24] P. Bolhuis and D. Frenkel, J. Chem. Phys. , 9869(1994).[25] T. Biben, P. Bladon, and D. Frenkel, J. Phys. Condens.Matter , 10799 (1996).[26] M. Dijkstra, D. Frenkel, and D. Waals, Phys. Rev. Lett. , 298 (1994).[27] J. C. Shelley and G. N. Patey, J. Chem. Phys. , 7656(1995).[28] J. C. Shelley and G. N. Patey, J. Chem. Phys. , 8265(1994).[29] N. Combe, T. J. H. Vlught, P. R. ten Wolde, andD. Frenkel, Mol. Phys. , 1675 (2003).[30] M. Murat and G. S. Grest, Macromolecules , 4054(1989).[31] S. Milner, Science , 905 (1991).[32] P.-Y. Lai and K. Binder, J. Chem. Phys. , 9288 (1991).[33] R. Netz and M. Schick, Macromolecules , 5105 (1998).[34] M. Patra and P. Linse, Nano Lett. , 133 (2006).[35] K. Binder and A. Milchev, J. Polym. Sci. Part B Polym.Phys. , 1515 (2012).[36] F. Lo Verso, L. Yelash, and K. Binder, Macromolecules , 4716 (2013).[37] S. Alexander, J. Phys. Fr. , 983 (1977).[38] S. T. Milner, T. a. Witten, and M. E. Cates, Europhys.Lett. , 413 (1988).[39] I. Coluzza and J.-p. Hansen, Phys. Rev. Lett. ,016104 (2008).[40] S. Consta, N. B. Wilding, D. Frenkel, andZ. Alexandrowicz, J. Chem. Phys. , 3220 (1999).[41] S. Consta, T. J. H. Vlugt, J. Wichers Hoeth, B. Smit,and D. Frenkel, Mol. Phys. , 1243 (1999).[42] D. Frenkel, Proc. Natl. Acad. Sci. , 17571 (2004).[43] K. Esselink, L. D. J. C. Loyens, and B. Smit, Phys. Rev.E51