Designed-walk replica-exchange method for simulations of complex systems Ryo Urano1 and Yuko Okamoto1,2,3,4 1Department of Physics, Graduate School of Science, Nagoya University, Nagoya, Aichi 464-8602, Japan 2Structural Biology Research Center, Graduate School of Science, Nagoya University, Nagoya, Aichi 464-8602, Japan 3Center for Computational Science, Graduate School of Engineering, 5 Nagoya University, Nagoya, Aichi 464-8603, Japan and 1 4Information Technology Center, Nagoya University, Nagoya, Aichi 464-8601, Japan 0 2 Weproposeanewimplementationofthereplica-exchangemethod(REM)inwhichreplicasfollow p a pre-planned route in temperature space instead of a random walk. Our method satisfies the e detailed balance condition in the proposed route. The method forces tunneling events between S the highest and lowest temperatures to happen with an almost constant period. The number of tunnelingcounts is proportional to that of therandom-walk REM multiplied by thesquare root of 0 moving distance in temperature space. We applied this new implementation to two kinds of REM 2 and compared the results with those of the conventional random-walk REM. The test system was a two-dimensional Ising model, and our new method reproduced the results of the conventional ] h random-walk REM and could adjust the tunneling counts by two times or more than that of the c random-walk REMby replica-exchange attempt frequency. e m Keywords: replica-exchange method(REM),randomwalk,efficiencyofREM,tunnelingevents - t a t s . t a INTRODUCTION m - d Theenhancementofconfigurationalsamplingensuresaccuracyandefficiencyofmoleculardynamics(MD)orMonte n Carlo (MC) simulations. Replica-exchange method (REM)[1–4] (or parallel tempering) is a popular way to improve o c efficiency of configurational sampling (for reviews, see, e.g., Refs.[5, 6]). This method achieves a random walk in [ temperaturespaceandallowsthesystemto overcomeenergybarriersbetweenlocal-minimumfreeenergystates. The number of tunneling events[7, 8] is defined to be the number of times where the simulation visits from the lowest 3 v temperature throughthe highesttemperature andback to the lowesttemperature (hereafter, we referto this number 2 as tunneling counts). The tunneling count is a measure of configurational sampling efficiency: the more tunneling 7 counts are obtained with a fixed number of simulation steps, the better the sampling is. However, as the number of 7 replicas increases, tunneling events require more steps by random walk. This is a difficulty for large systems because 0 they require many replicas in general. 0 1. In the Metropolis criterion with pseudo random numbers, randomness prevents the design of a trajectory in tem- 0 perature space. Moreover, the moving distance by a random walk is proportional to the square root of the number 5 of trials. Therefore, previous improvements of tunneling count and efficiency in the REM was fulfilled by approaches 1 such as the temperature selection[9–12], attempts of non-neighboring pairs, and increase of acceptance probabilities. : v Anexampleisall-paringmethod,whichtriesallcombinationsoftemperaturepairs[13,14]. AnotherexampleisGibbs i sampling heat-bath replica-exchange method[15]. Other examples are to employ the global balance condition such X as replica permutation method[16] with Suwa-Todo algorithm[17] and all possible pair exchange[18]. These methods r a reportedhigher tunneling counts thanthe conventionalREM.However,a productionofan equilibriumstate without random walk leads to another possibility of an improvedREM, which enables a simulation to determine a trajectory ofreplicasinatemperaturespace. Infact,someofspecialmanipulationofstatechangeswithoutthedetailedbalance condition shows a rapid convergence for thermal equilibrium[19, 20]. Based on Chaotic Boltzmann machine[21, 22], we have also proposed the deterministic replica-exchange method (DETREM)[23], which performs replica exchange based on a differential equation without pseudo random numbers. This method produced the same efficiency as the conventional REM. Recently, Spill et al. showed that a planned route trip in only one randomly selected replica during a simulation showed an improvement compared with the conventional REM[24]. We can generalize this idea so that planned route forallreplicaswill realizethe maximumefficiency withwide configurationalsampling,andwe here proposethe designed-walk replica-exchange method (DEWREM) by even-odd sequential replica exchange. 2 METHODS We now give the details of our methods. We prepare M non-interacting replicas at M different temperatures. Let thelabeli(=1, ,M)standforthereplicaindexandlabelm(=1, ,M)forthetemperatureindex. Werepresent ··· ··· the state of the entire system of M replicas by X = x[1] , ,x[M] , where x[i] = q[i],p[i] are the set of coordinatesq[i] andmomentap[i] ofparticlesinreplicain(amt(t1e)m·p··eratmur(MeT)o). Theprmobabi(cid:8)lityweig(cid:9)hmtfactorforstate m X is given by a product of Boltzmann factors: M W (X)= exp[ β H(q[i],p[i])], (1) REM m(i) − Yi=1 whereβ (=1/k T )istheinversetemperatureandH(q,p)istheHamiltonianofthesystem. Weconsiderexchanging m B m a pair of replicas i and j corresponding to temperatures T and T , respectively: m n X = ,x[i], ,x[j], X′ = ,x[j]′, ,x[i]′, , (2) n··· m ··· n ···o→ n··· m ··· n ···o where x[i]′ q[i],p[i]′ ,x[j]′ q[j],p[j]′ , and p[j]′ = Tmp[j],p[i]′ = Tnp[i] [4]. n ≡n on m ≡n om′ qTn qTm Here, the transition probability ω(X X ) of Metropolis criterion for replica exchange is given by → ′ W (X ) ′ REM ω(X X )=min 1, =min(1,exp( ∆)), (3) → (cid:18) W (X)(cid:19) − REM where ∆=∆ =(β β )(E(q[i]) E(q[j])). (4) m,n n m − − Becauseeachreplicavisits varioustemperatures followedby the transitionprobabilityofMetropolis algorithm,REM performs a random walk in temperature space. We now review two REMs, which are based on random walks in temperature space. Without loss of generality, we can assume that M is an even integer and that T < T < < T . The conventional REM[1–4] is performed by 1 2 M ··· repeating the following two steps: 1. We performa conventionalMD orMC simulationofreplica i (=1, ,M)attemperature T (m=1, ,M) m ··· ··· simultaneously and independently for short steps. 2. Pairsofexchangeattemptsareselectedinreplicapairswithneighboringtemperatures,forexample,fortheodd pairs (T1,T2), (T3,T4), , (TM−1,TM) or even pairs (T2,T3), (T4,T5), , (TM−2,TM−1). ··· ··· All the replica pairs thus selected are attempted to be exchanged according to the Metropolis transition probability in Eqs. (3) and (4) with n=m+1. WerepeatSteps1and2untiltheendofthesimulation. Thecanonicalensembleatanytemperatureisreconstructed by reweighting techniques[25–27]. We next present the deterministic replica-exchange method (DETREM)[23]. Only Step 2 is different from the conventionalREM.Atfirst,weintroduceaninternalstatey asanindexofapairofreplicasiandj attemperatures m,n T and T , and consider the following differential equation: m n dy 1 m,m+1 =σ , (5) m dt 1+exp(∆ ) m,m+1 where t is a virtual time, ∆ is the same as in Eq. (4) with n = m+1, and the signature σ of the pair of m,m+1 m (T ,T )changesto1or 1tocontrolthesignatureofthechangeofy whichmonotonicallyincreasesordecreases. m m+1 m − In Step 2, instead of applying the Metropolis criterion in Eqs. (3) and (4), we solve the differential equation in Eq. (5)forthe internalstatesy 1,1 for(T ,T ),wherethe totalnumberofinternalstatesisM-1withthe m,m+1 m m+1 ∈{− } following pairs: (1,2), (2,3), , (M-1,M) for the random-walk DETREM and the pairs: (1,2), (3,4), , (M-1,M) ··· ··· and (2,3), (4,5), , (M-2,M-1) for designed-walk REM. The replica exchange is done as follows[23]: ··· if updated y ≷ 1, then (T ,T ) (T ,T ), m,m+1 m m+1 m+1 m ± → y y 1,σ 1. m,m+1 m,m+1 m ← ∓ ←∓ 3 Fortherandom-walkDETREM,ify performsexchanges,y isnottimeevolvedandy isevolved m,m+1 m+1,m+2 m+2,m+3 to avoid the leap exchange of temperature such as from T to T . m m+2 Finally, the designed temperature walk can be implemented to both conventional REM and DETREM (and other REMs) as follows. Namely, the designed-walk replica-exchange method (DEWREM) is performed by repeating the following steps. 1. We perform a conventionalMD or MC simulation of replica i (=1, ,M) at temperature T (m=1, ,M) m ··· ··· simultaneously and independently for short steps. 2. Replica exchange is attempted for all the odd pairs (T1,T2), (T3,T4), , (TM−1,TM). ··· 3. RepeatSteps1and2untilalloddpairsperformreplicaexchangeexactlyonce. Namely,onceapairisexchanged, theexchangedpairstopsexchangeattemptsandkeepperformingthesimulationinStep1withthenewtemperatures. Replica exchange attempt in Step 2 is repeated until all the other odd pairs finish exchanges. 4–6. Repeat Steps 1–3 where the odd pairs in Steps 2 and 3 are now replaced by the even pairs (T ,T ), (T ,T ), 2 3 4 5 , (TM−2,TM−1). ··· 7. The cycle of Steps 1 to 6 is repeated until the number of cycles is M, which is equal to the tunneling count and all replicas have the initial temperatures. 8. Begin the above cycle of Steps 1–7 with Steps 1 to 3 and Steps 4 to 6 interchanged. These eight steps are repeated until the end of the simulation. The schematic picture of this procedure is shown in Fig. 1. We remark that Step 8, namely, reversing the cycle of Steps 1–3and4–6,is necessaryforthe detailedbalancecondition, becausethe entering statesarethe sameasleaving states. Forexample,thestate(x[1],x[2],x[3],x[4],x[5],x[6])isreachedfromonlytwostates(x[1],x[2],x[3],x[4],x[5],x[6]), 1 2 3 4 5 6 2 1 4 3 6 5 [1] [2] [3] [4] [5] [6] [i] (x ,x ,x ,x ,x ,x ) only and makes transition to the two states as shown in Fig. 1, where x is the state 1 3 2 5 4 6 m of replica i at temperature T . This exchange procedure satisfies the detailed balance condition for replica and m temperature pair because the trials of exchange pair γ i(m) i(m+1) ω (xi(m),xi(m+1)) (xi(m),xi(m+1)) (cid:16) → (cid:17) (cid:16) m m+1 → m+1 m (cid:17) =γ i(m+1) i(m) ω (xi(m),xi(m+1)) (xi(m),xi(m+1)) (6) (cid:16) → (cid:17) (cid:16) m+1 m → m m+1 (cid:17) isequalintherouteasisshowninFig. 1,whereγ(i(m) i(m+1))istheselectedprobabilityoftheexchangeattempt. → Notethatforexactlythesameconditionalprobabilityofoddpairreplicaexchange,(T1,T2),(T3,T4), ,(TM−1,TM), ··· ∆is Mk=/12(β2k−1−β2k)(E2k−E2k−1), wherek =1,··· ,M/2. However,becauseeachreplicais non-interactingwith otherPreplicas, waiting for other exchanges does not influence the transitions of others. This sequential exchange achieves one tunneling count when M cycles for each replica are finished. In theory, the estimated ratio of tunneling count between the odd-even sequential exchange and the conventional random walk is given by N PDEW trial× correction TCsequential 2M = N , (7) TC √N PRW ∝ trial random walk trial× correction p 2M where N is the number of exchange attempts, PDEW is the correction for waiting for all the replica exchanges trial correction in Steps 3 and 6, and PRW is the correction for the deviation of random-walk probability from the value 1/2. correction RESULTS In order to test the effectiveness of the present methods, we applied them to the 2-dimensional Ising model. The lattice size in a square lattice was 128 (hence, the number of spins was N = 1282 = 16384). We have performed conventional random-walk simulation and designed-walk (DEWREM) simulation of both conventional REM and DETREM. We have also performed a mixed random-walk and designed-walk simulation of DETREM, where we repeated the two walks alternately. The total number of replicas M was 40 and the temperatures were 1.50, 1.55, 1.60, 1.65, 1.70, 1.75, 1.80, 1.85, 1.90, 1.94, 1.98, 2.01, 2.04, 2.07, 2.10, 2.13, 2.16, 2.19, 2.22, 2.25, 2.28, 2.31, 2.34, 2.358, 2.368, 2.38, 2.40, 2.42, 2.44, 2.47, 2.51, 2.57, 2.63 ,2.69, 2.75, 2.82 ,2.90, 3.00, 3.10, and 3.15. Boltzmann ∗ constant k and coupling constant J were set to 1. Thus, β =1/k T =1/T =β , and the potential energy is given B B byE(s)= s s ,wheres = 1,andthesummationistakenoverallthenearest-neighborpairsinthesquare − <i,j> i j i ± lattice. P 4 Fortheconventionalrandom-walkREMandDETREM,replica-exchangeattemptwasmadeevery1MCstep. One MC step here consists of one Metropolis update of spins. The total number of MC steps for all the simulations was 100,000,000. To integrate Eq. (5), we used the fourth-order Runge-Kutta method with virtual time step dt = 1. For DEWREM, replica-exchange attempt was made every 10, 20, 50, 100, and 150 MC steps in the conventional REMsimulationsandevery10,20,50,and100MCsteps inthe DETREMsimulations(seeTableI).Themixed-walk simulationwasperformedinwhichafter4M(=160)even-oddorodd-evencyclesofdesigned-walksimulations(replica- exchange attempt was made at every 20 MC steps) were performed, 200,000 MC steps (which roughly corresponds to 2M cycles) ofrandom-walksimulations (replica-exchangeattempt was made at every 1 MC step) were performed, andthenthis procedurewasrepeated. For reweightinganalyses[25–28], the totalof10,000spinstate data weretaken with a fixed interval of 1,000 MC steps at each temperature from the REM simulations. Table I lists the mean tunneling counts in temperature space and energy space per replica for each method, which isthenumberoftimeswherethereplicasvisitfromthelowesttemperaturethroughthehighesttemperatureandback to the lowest, and the mean energy at the lowest temperature through the mean energy at the highest temperature and back to the mean energy at the lowest temperature during the simulation, respectively. The mean tunneling counts in temperature space per replica of the designed-walk simulations at every 100 MC attempts were about four times larger than random-walk Metropolis REM and twice larger than random-walk DETEM. On the other hand, the mean tunneling counts in energy space per replica of the designed-walk simulations at every 100 MC attempts were about six times larger than random-walk Metropolis REM and twice larger than random-walk DETEM. These large tunneling counts imply that in designed-walk method all replicas traversed more efficiently in temperature space,and ourdesignto maximize the tunneling counts for all replicaswithout randomwalks wassuccessful. For the mixed-walk simulation, the maximum tunneling count was about twice larger than that of random-walk DETREM. The mean tunneling count was almost the same as that of designed-walk DETREM. We next examine physical quantities obtained fromthe designed-walksimulations with various replica-exchangeattempt frequencies and mixed walksimulationsandcomparethemtothosefromtheconventionalrandom-walksimulations. Fig. 2(a)andFig. 2(b) comparethecanonicalprobabilitydistributionsofenergydensityǫ=E/N atfourtemperaturesasafunctionofT and the average total energy density ǫ as a function of T obtained from the random-walk and designed-walk simulations of REM and DETREM, for DETREM including the mixed-walk simulation. Most of the probability distributions of energy density are the same. However, in the DEWREM simulations, the probability distributions with high frequencyreplica-exchangeattempts atT =2.25nearthe exactcriticaltemperature T =2.269aredeformedslightly c compared to the results of the random-walk simulation. This can also confirm the results of heat capacity as shown in Fig. 3. Fig. 2(c) and Fig. 2(d) show the average total energy density ǫ =E/N as a function of T from the same simulations. They were obtained by the reweighting techniques[25–28]. The average total energy densities in all the simulations are the same. Fig. 3(a)andFig. 3(b)showthespecificheatC asafunctionofT duringtheconventionalREMsimulationsandthe DETREM simulations, respectively. This shows that designed-simulationwith shorter replica-exchangeinterval such as every 10 and 20 MC steps underestimated the heat capacity near the critical temperature although the transition point is sufficiently similar to the exact critical temperature at T = 2.269. As the intervals of replica-exchange c attempts are longer, the accuracy of heat capacity is higher. Moreover, the combination of the random-walk and designed walk also increased the accuracy. This suggests that the designed-walk replica-exchange attempts caused correlation between replicas. The correlation seems to be very strong near the critical temperature. As a result, the heatcapacityis underestimatedslightly. However,these physicalquantities showthatmixed-walksimulationcan increasetheaccuracyofresultsandthenumberoftunnelingcountsandDEWREMsimulationissuitedforsimulations with longer replica-exchange time intervals between replica-exchange attempts. CONCLUSIONS Inthisarticle,weproposedanalgorithmbyadesignedtemperaturerouteforreplica-exchangemethods. Theaimfor developing this method was to maximize the tunneling counts. We reproducedthe results of conventionalREM. Our designed-walk simulations showed that the tunneling counts increased by at least twice more than the random-walk simulations. This method also showedthat in REM replicas have no interactionamong replicas but replica-exchange trialswithartificialorderscauseacorrelationbetweenreplicas. Bydecreasingthecorrelationsintroducingmixed-walk simulation,usinglonger-intervalreplica-exchangetrials,orotherparingofreplicas,thismethodwillbemoreefficient. Inafuturework,wewillgivetheformulationfordesignedroutesformultidimensionalreplica-exchangemethod[29], which is a power tool in massive parallel computing[30–32]. This will be useful for creating another design route for non-neighboring update for exchanges. 5 FIG.1. AnschematicpictureoftimeseriesoftemperatureindicesinDEWREMwith6replicas. Theleftcyclebeginswiththe temperature exchange of odd index pairs (T1,T2), (T3,T4), and (T5,T6), then tries with even pairs (T2,T3) and (T4,T5). The right cycle begins with even pairs and next tries odd pairs. They are the reverse cycles of each other and their combination satisfies thedetailed balance condition of replica exchange. TABLE I. Themean numberof tunneling countsper replica. TC Random walk Designed walk Mixed T TC Met DETREM Met DETREM DETREM Interval 1 100 1 100 20 50 100 150 20 50 100 1 & 20 Mean 173 37 178 58 292 197 131 99 231 144 93 293 ± SD 10 3 9 5 56 41 27 21 48 30 20 6 E TC Met DETREM Met DETREM DETREM Interval 1 100 1 100 20 50 100 150 20 50 100 1 & 20 Mean 55 13 55 34 80 79 77 70 81 79 69 75 ± SD 5 2 5 4 8 6 7 5 6 6 4 5 T TC, E TC, Interval,SD,and Met stand for tunnelingcountsin temperature space, tunneling count in energy space, the numberof MC stepsbetween replica-exchange attempts, standard deviation, and REM based on Metropolis criterion, respectively. The frequency (1 & 20) of Mixed means that it was 1 MC step for random walk REM and 20 MC steps for designed-walk REM. ACKNOWLEDGMENTS Some of the computations were performed on the supercomputers at the Institute for Molecular Science, at the Supercomputer Center, Institute for Solid State Physics, University of Tokyo, and Center for Computational Sci- ences, University of Tsukuba. This work was supported, in part, Grants-in-Aid for Scientific Research (A) (No. 25247071),for Scientific Researchon Innovative Areas (“Dynamical Ordering & Integrated Functions”), Programfor Leading Graduate Schools “Integrative Graduate Education and Research in Green Natural Sciences”, and for the Computational Materials Science Initiative, and for High Performance Computing Infrastructure from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. [1] K.Hukushimaand K. Nemoto, J. Phys.Soc. Jpn.65, 1604 (1996). 6 FIG. 2. Probability distributions of energy density at four temperatures (from left to right, 1.50, 2.10, 2.25, 2.47, and 3.15) obtained from (a) REM and (b) DETREM simulations including the mixed-walk simulation, and average energy density as a functionofT from(c)REMand(d)DETREMsimulationsincludingthemixed-walksimulationandtheerrorbarsaresmaller than thesymbols. [2] R.H. Swendsen and J.-S. Wang, Phys. Rev.Lett. 57, 2607 (1986). [3] C. J. Geyer, Comput. Sci. Stat.: Proc. 23rd Symp.Interface, Interface Foundation, Fairfax Station, VA , 156 (1991). [4] Y.Sugita and Y. Okamoto, Chem. Phys. Lett.314, 141 (1999). [5] A.Mitsutake, Y. Sugita, and Y. Okamoto, Biopolymers 60, 96 (2001). [6] Y.Iba, Int.J. Mod. Phys.C 12, 623 (2001). [7] B. A. Berg and T. Neuhaus,Phys. Rev.Lett. 68, 9 (1992). [8] A.Mitsutake, Y. Sugita, and Y. Okamoto, J. Chem. Phys. 118, 6676 (2003). [9] S.Trebst, D. A.Huse, and M. Troyer, Phys. Rev.E 70, 046701 (2004). [10] H.G. Katzgraber, S. Trebst, D.A. Huse, and M. Troyer, J. Stat. Mech.: Theory. Exp.2006, P03018 (2006). [11] S.Trebst, M. Troyer, and U. H.E. Hansmann,J. Chem. Phys.124, 174903 (2006). [12] W. Nadlerand U.H. E. Hansmann,Phys. Rev.E 76, 065701 (2007). [13] F. Calvo, J. Chem. Phys. 123, 124106 (2005). [14] P.Brenner, C. R.Sweet, D. V. Handorf, and J. A. Izaguirre, J. Chem. Phys.126, 074103 (2007). [15] J. D. Chodera and M. R.Shirts, J. Chem. Phys.135, 194110 (2011). [16] S.G. Itoh and H. Okumura,J. Chem. Theory Comput. 9, 570 (2012). [17] H.Suwa and S. Todo, Phys. Rev.Lett. 105, 120603 (2010). 7 FIG.3. Specificheat C asafunction of T from the(a) REM,(b)DETREM simulations includingthemixed-walk simulation. The error bars are smaller than the symbols. The exact results for L =128 (black curves) were obtained by Berg’s program [33] based on Ref. [34]. [18] H.X. Kondoand M. Taiji, J. Chem. Phys.138, 244113 (2013). [19] E. P. Bernard, W. Krauth, and D. B. Wilson, Phys.Rev.E 80, 056704 (2009). [20] K.S. Turitsyn,M. Chertkov, and M. Vucelja, Physica D 240, 410 (2011). [21] H.Suzuki,J. Imura, Y. Horio, and K. Aihara, Sci. Rep. 3(2013). [22] H.Suzuki,Phys. Rev.E 88, 052144 (2013). [23] R.Urano and Y. Okamoto, e-printarXiv:1412.6959 (2014). [24] Y.G. Spill, G. Bouvier, and M. Nilges, J. Comput. Chem. 34, 132 (2013). [25] A.M. Ferrenberg and R. H.Swendsen, Phys.Rev. Lett.63, 1195 (1989). [26] S.Kumar, J. M. Rosenberg, D. Bouzida, R. H.Swendsen, and P. A.Kollman, J. Comput. Chem. 13, 1011 (1992). [27] A.Mitsutake, Y. Sugita, and Y. Okamoto, J. Chem. Phys. 118, 6664 (2003). [28] M. R.Shirts and J. D.Chodera, J. Chem. Phys.129, 124105 (2008). [29] Y.Sugita, A. Kitao, and Y.Okamoto, J. Chem. Phys.113, 6042 (2000). [30] Y.M. Rhee and V.S. Pande, Biophys. J. 84, 775 (2003). [31] T. Vogel, Y. W. Li, T. Wu¨st, and D. P.Landau, Phys. Rev.Lett 110, 210603 (2013). [32] T. Vogel, Y. W. Li, T. Wu¨st, and D. P.Landau, Phys.l Rev.E 90, 023302 (2014). [33] B. A. Berg, Markov Chain Monte Carlo Simulations and Their Statistical Analysis (World Scientific, Singapore, 2004). [34] A.E. Ferdinand and M. E. Fisher, Phys.Rev.185, 832 (1969).