05_tempering Obias and Banzon 12 Tempering and Annealing in a Verdier-Stockmayer Polymer E. R. Obias and R. S. Banzon Structure and Dynamics Group, National Institute of Physics, University of the Philippines, Diliman, Quezon City 1101 E-mail: erobias@up.edu.ph; rbanzon@nip.upd.edu.ph ABSTRACT Science Diliman (July-December 2004) 16:2, 12–16 Two Monte Carlo methods, simulated annealing and parallel tempering, were applied to a Verdier- Stockmayer polymer. The efficiency of the two algorithms in exploring the lowest energy state possible for the model polymers was measured by the number of energy-degenerate configurations (configurations that have the same energy but are structurally different). Parallel tempering consistently explored more energy-degenerate configurations as compared with simulated annealing. INTRODUCTION Monte Carlo (MC) methods are widely used in a number of simulations in statistical physics (Liang & Wong, 1997; Marinari & Parisi, 1992; Hansmann, 1997). This paper considers the applications and limitations of two MC methods, simulated annealing (SA) (Kirkpatrick et al., 1983) and parallel tempering (PT) (Frenkel & Smit, 2002), applied to a two- dimensional (2D) polymer that follows the Verdier- Stockmayer (Verdier & Stockmayer, 1962; Gould & Tobochnik, 1996) model. The main task of numerical simulations on this matter is to study the polymer’s low-energy conformations, which depend on how the phase space is comprehensively explored. However, systems of this type have energy landscapes characterized by many deep local minima separated by high-energy barriers. At low temperatures, traditional Monte Carlo and molecular-dynamics simulations tend to get “trapped” in one of these minima, thereby hampering simulations to explore the polymer’s ground-state configurations. The two algorithms discussed in this paper, which are modifications of the original Metropolis algorithm (Metropolis et al., 1953), hope to alleviate the problem by minimizing the probability of settling to these local minima and continue to search for the global minimum. In SA, the temperature is initially high so that a relatively large percentage of the random steps that result in an increase in the energy will be accepted, thus making the polymer move freely. The temperature is then gradually lowered, hence the term annealing, until the target temperature is reached and the polymer is able to relax at its most stable conformation. In PT, simulation is done over n systems at different temperatures. Multiple copies of the polymer run in parallel until a “swap” is imposed, hence the term parallel tempering. Temperature differences are relatively small enough to allow the occasional swapping of two neighboring systems with different temperatures. METHODOLOGY The polymer used in the simulation follows a 2D Verdier-Stockmayer model. Each monomer occupies one of the vertices in a 2D lattice. Initially, the polymer is completely unfolded. Local movement is one monomer at a time and self-avoiding. Figure 1 Tempering and Annealing 13 illustrates the two types of movements: end move (monomers 1 and 8) and corner move (monomers 4 and 7). Total energy will be computed as the sum of all nonbonded nearest-neighbor interactions. In both SA and PT, a trial move is done by randomly selecting a monomer and moving it to a valid site subject to the types of movements depicted in Fig. 1. The move is accepted if a random number (0,1) is less than the Boltzmann probability exp(–DE/T ) {Boltzmann factor set to unity}, where DE is the change in total energy and T is the temperature. If the change in energy is small enough or the temperature is high enough, then there is a high probability of acceptance. In PT, an additional trial move attempts to “swap” neighboring temperature polymers. This trial move is done by randomly selecting two neighboring temperature polymers and attempting to swap their configurations. The swap move is accepted if a random number (0,1) is less than the probability exp[(Ei–Ej)(bi– bj)], where b = 1/T. If there is enough overlap between these two systems, then there is a high probability of accepting the swap. The annealing schedule for simulated annealing is linear, every 5x105 MC steps, while the swapping ratio adapted for PT is 90% MC moves: 10% swapping moves (Frenkel & Smit, 2002). For a fair comparison, both algorithms have the same steps per temperature and the same total steps. The experiment consists of two parts. The first part had the interaction energies randomly generated [–4, – 2] and assigned to each monomer type (i.e., a hypothetical polymer is generated by choosing a sequence of N integers at random from the range 1– 20, corresponding to the 20 different possible amino acids) and temperatures ranged [0.5,10] in increments of 0.5 (Giordano, 1997) to determine the transition region using both methods. The second part had the interaction energies fixed to –2 and temperatures were limited in the range [1,5] to investigate further the ability of the two algorithms to explore the lowest energy state possible. The lengths of the polymers used were 15 and 30 monomer units. RESULTS AND DISCUSSION Determination of the transition region Figures 2 and 3 show average energy and end-to-end length of the 15-monomer chain using SA and PT. A transition region is observed in the region around T~2, both for the energy and end-to-end length graphs, as seen in the sharp curves in the plots. Each point is an average of 5x105 Monte Carlo steps. The polymer has its lowest energy state at the lowest temperature. Average end-to-end length decreases with temperature as seen in Fig. 3. This means that the polymer has a relatively greater degree of folding at low temperatures. Note that in PT, each point in the graph, which is also an average of 5x105 MC and swapping steps, is in a different system. The fluctuations can be eliminated by averaging more MC steps and performing more experiments. 1’ 1 2 1’ 4 3 4 5 6 7 8’ 8 7’ Fig. 1. Possible movements of a polymer with eight monomers. Initial position shown with solid lines and unprimed monomer numbers; possible moves are indicated by dotted lines and primed monomer numbers. 0 1 2 3 4 5 6 7 8 9 10 Temperature A ve ra g e en er g y -22 -20 -18 -16 -14 -12 -10 -8 -6 -4 Fig. 2. Average energy as a function of temperature. Obias and Banzon 14 Investigation in the transition region To further investigate the algorithm’s ability to explore the lowest energy state configurations, we must be able to easily identify different configurations in a region where, expectedly, a large number of configurations may occur. One way to do this is to fix the interaction energies to some number, in this case –2 so there is a high probability of acceptance, and limit the thermodynamical states in the vicinity of the transition region in the temperature range [1,5]. Figure 6 shows a 1x103 point representative of energies using SA with interaction energies fixed to –2 and temperatures annealed from T = 5 to T = 1 in increments of 1 every 5x105 MC steps for a total of 2.5x106 MC Fig. 3. Average end-to-end length as a function of temperature. 0 1 2 3 4 5 6 7 8 9 10 Temperature A ve ra ge e nd -t o- en d le ng th 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 Fig. 4. Energy at T = 1 using SA. 0 1000 2000 3000 4000 5000 MC steps (x100) E n er g y -26 -24 -22 -20 -18 -16 -14 -12 -10 -8 -6 Fig. 5. Energy at T = 1 using PT. 0 1000 2000 3000 4000 5000 MC steps (x100) E n er g y -25 -20 -15 -10 -5 Although the two methods produced very similar results, they are different in terms of their ability to explore the lowest energy states. Figures 4 and 5 show a 5x103 points representative of the 5x105 energies sampled at a temperature equal to 1 for the two algorithms. Each point represents one energy state. In Fig. 4, it is noticeable that some states are visited repeatedly as shown by consistent values lined horizontally in the graph even though we had initially set the interaction energies randomly. The polymer is somehow “trapped” at some states that it tends to repeatedly fold at this configuration. This repeated visits in some states are probable since the probability that they are degenerate is small considering the short length of the polymer. This is not seen in Fig. 5, as there are more energy states represented at the low-energy states. Fig. 6. 1x103 point representative of energies when interaction is fixed to –2 using SA. Temperature 1 2 3 4 5 -16 -14 -12 -10 -8 -6 -4 -2 0 E n er g y Tempering and Annealing 15 steps. One can clearly identify the nine possible energy states from the plot with –16 as the lowest possible energy. Figure 7 shows a similar pattern for parallel tempering with –16 also as the lowest possible energy. In Fig. 7, five systems with temperatures in the range [1,5] are represented with 5x105 MC steps (including the swaps) for each system. With no distinct criteria for comparison on their ability to explore the lowest energy states, we resort to counting structural differences called energy-degenerate configurations (EDC) explored in the lowest energy state. Figure 8 shows six examples of EDCs in 15 monomers with eight nearest-neighbor interactions and energy = –16. Using the same parameters, Fig. 9 presents a comparison of 20 independent experiments on the number of EDCs explored by the two algorithms at T = 1 in the simulation, where the energy of interaction is fixed to –2 and temperature ranged [1,5]. Rotation and mirror images are considered EDCs so the number may be overestimated by a factor of 8 assuming all types of configurations are sampled. As seen from the graph, PT is consistently able to explore more distinct configurations as compared with SA. To verify the consistency of the performance in exploring EDCs at the lowest energy state, we simulated a longer polymer consisting of 30 monomer units. Trial experiments showed, however, that the lowest energy state is difficult to achieve using the two algorithms as the polymer is considerably long, thus giving it a much larger energy landscape and more local minima to overcome. Figures 10 and 11 show a comparison of two lowest energy configurations obtained from trial simulations, which are –38 and –40. In Fig. 10, experiment 4 using SA did not reach the energy state equal to –38, while all experiments in PT did. In Fig. 11, 9 out of 20 experiments in PT managed to find the lowest energy state possible with the energy Fig. 7. 1x103 point representative of energies when interaction is fixed to –2 using PT. Temperature 1 2 3 4 5 -16 -14 -12 -10 -8 -6 -4 -2 0 E ne rg y Fig. 8. Examples of energy-degenerate configurations (EDC) in 15 monomers with energy = –16. Experiment 2 0 4 300 6 450 8 10 12 14 150 16 18 20 N um be r of E D C s Fig. 9. Comparison on the number of EDCs explored at the lowest possible energy state (T = 1) in 15 monomers with energy = –16. Experiment 2 0 4 200 6 250 8 10 12 14 100 16 18 20 N um be r of E D C s Fig. 10. Comparison on the number of EDCs explored at the “second” lowest possible energy state (T = 1) in 30 monomers with energy < –38. 150 50 Obias and Banzon 16 equal to –40 as compared with simulated annealing which only had three. This further supplements the advantage of PT as a better energy landscape explorer. CONCLUSION Information on the number of EDCs demonstrate to be a useful criterion in comparing algorithms whose objective is to study details of the energy landscape of a system. Parallel tempering is consistently able to search more EDCs at the lowest energy state as compared with SA. However, both algorithms find it difficult to find the lowest energy state for longer polymers due to its larger energy landscape and more local minima to overcome. Further investigations in the annealing schedule and swapping ratio will be made. REFERENCES Frenkel, D. & B. Smit, Understanding molecular simulation. 2nd ed. Academic Press, San Diego: Chap. 14. Giordano, N., 1997. Computational Physics. Prentice-Hall Inc, New Jersey: Chap. 11. Gould, H. & J. Tobochnik, 1996. An introduction to computer simulation methods. 2nd ed. Addison-Wesley Publishing Company Inc., New York: Chap. 12. Hansmann, U., 1997. Parallel tempering algorithm for conformational studies of biological molecules. Chem. Phys. Lett. 281: 140–150. Kirkpatrick, S., C.D. Gelatt, Jr., & M.P. Vecchi, 1983. Optimization by simulated annealing. Science. 220: 671–680. Liang, L. & W.H. Wong, 1997. Dynamic weighting in Monte Carlo and optimization. Proceedings of the National Academy of Sciences. 94: 14220–14224. Marinari, E. & G. Parisi, 1992. Simulated tempering: A new Monte Carlo scheme. Europhys. Lett. 19: 451–455. Metropolis, N., A. Rosenbluth, M. Rosenbluth, A. Teller, & E. Teller, 1953. Equation of state calculations by fast computing machines. J. Chem. Phys. 21: 1087–1092. Verdier, P.H., & W.H. Stockmayer, 1962. Monte Carlo calculations on the dynamics of polymers in dilute solution. J. Chem. Phys. 36: 227. Fig. 11. Comparison on the number of EDCs explored at the lowest possible energy state (T = 1) in 30 monomers with energy = –40. Experiment 2 0 4 5 6 6 8 10 12 14 3 16 18 20 N um be r of E D C s 4 2 1