Developments in Wang-Landau simulations of a simple continuous homopolymer

Developments in Wang-Landau simulations of a simple continuous homopolymer

D. T. Seaton; S. J. Mitchell; D. P. Landau

Center for Simulational Physics, University of Georgia, Athens, GA 30602, USA

ABSTRACT

The Wang-Landau method is used to study thermodynamic properties of a three-dimensional flexible homopolymer chain with continuous monomer positions. Results describing the coil-globule and solid-liquid transitions are presented for chain lengths up to N = 100. In order to elucidate the thermodynamic behavior, finite chain length effects and the influence of the energy range over which the density of states is determined are carefully analyzed. Simulation efficiency is also studied and it is shown that setting the natural logarithm of the final modification factor equal to 10^-6 is an appropriate choice for this model.

Keywords: Homopolymer chain; Wang-Landau method; Coil-globule transition; Solid-liquid transition

I. INTRODUCTION

Polymer systems have provoked much interest in an array of different disciplines, and in recent decades computer simulations have become powerful and mature. In particular, Monte Carlo methods [1] have shown great promise in the investigation of thermodynamic properties of these systems [2]. However, traditional Monte Carlo methods have had trouble sampling the rough energy landscapes associated with some polymer models. The Wang-Landau method [3-5] has proven quite useful in overcoming these inefficiencies. Since the target of our current study is understanding the behavior of a flexible homopolymer when the temperature is lowered and multiple states begin to compete, Wang-Landau sampling is very well suited for this problem.

In two recent studies, the Wang-Landau method was used to investigate the phase behavior of two related models for a single homopolymer. Rampf et al. [6-8] used the bond-fluctuation model [9] to study single chain properties, while Parsons and Williams [10, 11] used a continuous homopolymer model. Both models had simple interaction potentials and both studies found coil-globule and solid-liquid transitions. The coil-globule transition [12] represents the change from an unwound, random coil, to a more compact, globular state. The solid-liquid transition is the collapse of the polymer from a compact globular state into a nearly frozen, "crystalline" state. Representative configurations for these states are shown in Fig. 1. Additionally, these studies [6-8, 10, 11] provide some evidence of intermediate states for large chain lengths. Parsons and Williams also made comparisons between their results and the results from Rampf et al., highlighting differences between the two models. These differences were originally quite dramatic, but recent work from Rampf et al. [8] shows that by increasing the interaction distance of their non-bonded potential, the coil-globule and solid-liquid transitions behave more like those in the Parsons and Williams study. Nonetheless, the two studies still show significant differences in their thermodynamic behavior.

The aim of the work in this article is to understand the behavior of a simple continuous homopolymer, while also making comparisons to the results found in the previously discussed studies [6-8, 10, 11]. The results reported here show conflicting conclusions compared with the continuum model results from Parsons and Williams. The nature of the solid-liquid and the coil-globule transitions appear to be different, even though the models are nearly identical. In order to understand these differences, the work presented here focuses on making a systematic study of a single continuous homopolymer using Wang-Landau sampling.

II. MODEL

The polymer model we consider consists of a chain of N identical monomers, where each monomer position may vary continuously in three-dimensional space. There are two types of interactions: a bonded interaction between nearest-neighbor monomers and a non-bonded interaction between non-bonded monomers. The Hamiltonian for this model is

where the first term represents the bonded energy, summed over all N-1 bonds using the bond potential U^Bond(l_i), which is a combination of a finitely extensible nonlinear elastic (FENE) potential [2, 13] and a Lennard-Jones potential. The second term in represents the non-bonded energy, summed over all interactions between non-bonded monomer pairs (i,j), where U^LJ(r_ij) is a Lennard-Jones potential, and r_ij is the distance between two non-bonded monomers. The Lennard-Jones potential is given by

where the subscript ij represents an i-j pair, r_ij is the distance between these two monomers, e is the well depth, and s is the distance at which the potential function is a minimum. Dimensionless units are used, where e = 1 and s = 2^-1/6, and the minimum is then located at r_ij= 1.

The FENE potential is given by

where R_o is the finite extensibility and is set to R_o = 1.2, k is the stiffness constant and is set to k = 2, l_iis the bond length and has an equilibrium distance of l_o = 1. All parameters were chosen such that the minima of both U^Bond(l_i)and U^LJ(r_ij) are at a distance of 1. For comparison, the FENE and Lennard-Jones potentials are plotted in Fig. 2.

III. METHOD

Wang-Landau sampling is a Monte Carlo technique that was developed as a temperature independent, iterative method with the ability to sample rough energy landscapes [5, 14]. The technique was originally implemented for discrete models [3], and in order to perform simulations of a continuous model a few issues must be addressed. (Recent work has focused on optimizing the Wang-Landau method for continuous systems [15, 16]; the implementation for this study was directed at understanding the effects of finite length and restricted energy range for the random walkers, so an implementation close to the original method was adopted.)

The Wang-Landau sampling approach is based on calculating the density of states g(E), where E is the energy of the system. An initial estimate of the density of states g(E) = 1.0 is made. The algorithm proceeds by generating trial states, using appropriate "moves", and accepting them with transition probability

where E_i is the initial energy and E_f is the energy of the trial state. If the trial move is accepted, the density of states is multiplied by a modification factor f, i.e. g(E_f) ® g(E_f)×f, and a histogram H(E) is also updated, H(E_f)® H(E_f) + 1. If the trial state is rejected, the same procedure is followed for the initial energy E_i. This process is repeated until the histogram H(E) is sufficiently flat, at which point the histogram is reset to zero, the modification factor is reduced by f ® , and the random walk continues. The initial value of the modification factor is set to f = e¹ and is eventually reduced to a value which is f_final » 1. The flatness of H(E) is defined as its minimum divided by its average, and in this study, the flatness criterion is set to 0.8. Because g(E) becomes extremely large in the actual computation, the logarithm of the density of states, ln[g(E)], is used to prevent overflow during the actual simulation. More details about Wang-Landau sampling can be found elsewhere [5, 14].

The energy range over which the random walk is allowed to proceed is an important issue. Some studies have split the range into intervals [7, 17] in order to make the random walk over each range converge faster. The "pieces" of the individual g(E)'s are then strung together to produce a comprehensive distribution for the overall g(E). The simulations in this report will be the foundation for future work, so the simplest approach was implemented by first establishing a lower and upper bound for the Wang-Landau method to sample. (The choice of a lower energy bound is important because with too low a value the histogram never becomes sufficiently flat.) In some continuous models the ground state energy is not known, and the value of the upper energy bound, beyond which contributions to the behavior near any polymer collapse are small, may also not be well known. In this study, simulations have been performed over a number of different energy ranges in order to determine how they affect the thermodynamic properties.

Once an energy range has been defined, the energy is divided into a discrete number of bins, which defines the binning resolution. Testing showed that a bin resolution of dE = 0.1 was generally a good choice. The only exception being for the chain length of N = 100, where dE = 0.2 was used to prevent the number of bins from becoming too large. (For example, for N = 100, a resolution of dE = 0.1 and an energy range of E/N = [-4.0,3.0] requires 7000 bins.) Since the effects of bin resolution have still not been well studied, attempts were made to keep the number of bins reasonably low for computational efficiency.

In each simulation the chain starts from a straight configuration. Changes are made to the polymer configurations using the following types of trial moves: single monomer displacement ("diffusion"), reptation, crank-shaft, and random pivot. One Monte Carlo sweep (MCS) consists of N attempted monomer displacement moves and a single attempt of each of the other types of moves. The moves are identical to those described in reference [18]. It is also important to note that an attempt was made to locate the approximate ground state energy for each chain length, insuring that the energy ranges used were reasonably appropriate in each simulation. This was done by simply allowing the Wang-Landau method to sample a large energy space with a small bin-resolution. Wang-Landau sampling has the ability to be optimized to find ground states in various systems[15], but in this work, only testing was performed and only approximate ground state energies are reported.

IV. SIMULATION RESULTS

In all figures that follow, data points represent averages over 10 or more independent runs. Fig. 3 shows the density of states for four different chain lengths: N = 10, 25, 50, and 100, where the maximum of each function has been normalized to ln(g(E)) = 0. The range of g(E) increases rapidly with chain length, and the dependence upon energy becomes particularly important at low energies. These are the results from early simulations where it was still unclear what the energy range should be for different chain lengths. Each data set is shown over a distinct energy range, where the minimum energy boundary, E_min/N, varies with chain length, but the maximum energy boundary is fixed at E_max/N = 0.0 for all chain lengths. The maximum energy was initially set to zero based on the work from Rampf et al. [6-8] and Parsons and Williams [10, 11]. The minimum energy was "pushed" as low as it was possible to simulate efficiently, i.e. with a maximum run time of 10 days.

A. Effect of Energy Boundaries

The choice of the minimum energy boundary, E_min/N, can dramatically affect thermodynamic quantities, e.g. the specific heat, C_V. In Fig.4, the specific heat is plotted versus temperature, T, using a chain length of N = 10, for three different values of the minimum energy boundary, E_min/N = -1.5, -1.7, and -1.9. The curve with the highest E_min/N shows a single rounded peak at approximately T » 0.5. As the minimum energy is lowered, noticeable differences in the specific heat appear. The curve with E_min/N = -1.7 begins to show two interesting features, a taller peak and a shoulder. As the energy boundary is lowered to E_min/N = -1.9, a clear peak appears at T » 0.2, while the shoulder remains at T » 0.5. These results using N = 10 show how dramatically the minimum energy boundary can affect thermodynamic quantities. This chain length allowed the minimum energy to be pushed very close to the ground state with, E_gs » -1.93.

The maximum energy boundary can also affect thermodynamic quantities more substantially than might be expected. In Fig. 5, the maximum energy boundary is varied for a chain length of N = 10. Each of the four curves in this plot has a minimum energy of E_min/N = -1.9, which is the lowest energy boundary from the previous graph, Fig. 4. Each set of data points represents a different maximum energy boundary, E_max/N = 0.0, 1.0, 2.0, and 3.0. Changing E_max/N does not affect the specific heat at low temperatures, but does affect the quantity at higher temperatures. The shoulder seen in Fig. 4 has now been stretched out slightly, compared to the curve with E_max/N = 0.0. The three curves with larger values for this boundary have smaller, but still noticeable, differences. The inset in Fig. 5 gives a better view of these differences at higher T. Metropolis sampling was also performed for the temperature range in the inset, and the results of these simulations agreed with the data from the E_max/N = 3.0 curve to within the respective error bars.

B. Chain Length and Transitions

The results for the small chain length shown in the previous figures illustrate some of the issues one must deal with when performing Wang-Landau simulations on this type of model. Once these features were better understood, larger chain lengths were simulated, where similar energy range studies were performed for each chain length. In Fig. 6, the specific heat is plotted versus temperature for four different chain lengths, N = 10, 25, 50, and 100. The two features which became apparent while varying the minimum energy boundary are now very different. The small peak noticed for N = 10 has become quite sharp for larger lengths, suggesting that a first-order transition might occur in the thermodynamic limit of infinite chain length, and the small shoulder has now been stretched out to even higher temperatures, indicating an even broader 'intermediate' region. The N = 100 curve displays two prominent features. The first is the low-temperature peak at T » 0.4, which is believed to represent a solid-liquid transition. The second feature is the coil-globule transition which is represented by the shoulder at T » 2.0. Recent work with larger chain lengths [8, 11] has also shown that the area in between these transitions can contain other interesting meta-stable states; however, for the chain lengths presented above, such features were not observed. Simulations of larger chain lengths, accompanied by rigorous structural analyses, are still needed to better understand all of these interesting regions.

Note that in Fig. 6, the energy ranges for the respective ln[g(E)] have been adjusted to follow the results seen in the previous figures. e.g., for N = 100, the energy range is E/N = [-4.0,3.0]. E_max/N = 3.0 for all of the chain lengths, while E_min/N are the same for the chain lengths given in Fig. 4.

C. Simulation Efficiency and the Modification Factor

In order to improve efficiency, much of the recent work on Wang-Landau simulations of polymer systems has involved determining how to best implement the method[16]. An important issue is the choice of the final value of the modification factor f, i.e. its value at the end of each simulation. The original study on discrete spin models [3] used a final value of ln[f] = 10^-8; but in most of the recent work on polymer models, a final modification factor of ln[f] = 10^-6 was used. The aim of this section is to better understand why 10^-6 is an appropriate choice for this model.

In order to resolve whether the selection of ln[f_min] = 10^-6 made in some studies was warranted, 50 independent runs were made for each of the following chain lengths: N = 10, 25, and 50. In each run, g(E) was recorded for each value of the modification factor, i.e. at the completion of each level of iteration. The final value of ln[f] was reduced to at least ln[f_min] = 10^-8, and the final values of ln[g(E)] were averaged together. In Fig. 7, the average values of ln[g(E)] are plotted as a function of energy for each of the chain lengths. The energy ranges follow the same criterion described in the previous section regarding energy boundaries, where the maximum and minimum cutoffs have been extended far enough to insure that the important thermodynamic features are reproduced correctly. Using these average values of ln[g(E)], a reference value of the specific heat, , is calculated for each chain length. These values will now serve as a reference which will be compared to specific heat curves calculated from the ln[g(E)] values at each iteration level of ln[f]. When comparing the reference with specific heat curves calculated at a particular value of ln[f], the absolute differences, , are then averaged and stored for each value of ln[f]. These difference values, for each ln[f], are then averaged for all 50 independent runs and error bars are calculated.

The top of Fig. 8 shows the absolute difference between and the calculated C_V for each value of ln[f], and the bottom of Fig. 8 shows the average number of MC sweeps versus ln[f]. The same 50 independent runs for each chain length were used to produce the curves in both plots. The top plot shows how the error decreases during the course of a simulation. For the N = 10 system length, the curve decreases as ln[f] decreases, but then levels off at ln[f] » 10^-6. The simulations for this system size were continued to a final value of the modification factor of ln[f] = 10^-10, and it is clear that any extra time spent running a single simulation past ln[f] » 10^-6 does not necessarily improve the final result (at least for our choice of flatness criterion). The other two chain lengths plotted were simulated to a value of ln[f] = 10^-8. Again, the average deviation levels off at ln[f] » 10^-6. However, the bottom plot reveals that a huge price is payed in terms of the number of MC sweeps needed to iterate ln[f] from 10^-6 to 10^-8. For N = 25, the number of MC sweeps executed after ln[f] = 10^-6 is twice the number executed before ln[f] = 10^-6. For N = 50 we observe a similar behavior, only the number of sweeps needed to reach ln[f] = 10^-8 is nearly 3 times larger than the number needed to reach ln[f] = 10^-6. Thus, when considering both efficiency and the error associated with Wang-Landau simulations, we find that ln[f_min] = 10^-6 is a very reasonable choice.

V. CONCLUSIONS

Results have been presented for Wang-Landau simulations of a simple flexible homopolymer for chain lengths up to N = 100. It was found that the energy range used by the Wang-Landau random walk can dramatically affect the shape of the specific heat. For instance, the solid-liquid transition can be easily hidden if the minimum energy is too high. Along those same lines, if the maximum energy boundary is too low, the specific heat at higher temperatures indicates too low a coil-globule transition temperature. (Although not shown in this article, similar effects are seen in other thermodynamic quantities.) We also find that a final value of the modification factor equal to ln[f] = 10^-6 is appropriate, given our choice of flatness criterion (0.8), and iterating to a still lower value only increases the CPU time with very little improvement in the final results.

The results presented here differ considerably from those reported by Parsons and Williams [10, 11], even though the models used are nearly identical. Interestingly, our results tend to have many more similarities with the work from Rampf et al. [6-8] using the bond fluctuation model, especially with regards to the their latest work [8]. This is particularly true for the solid-liquid transition, where a sharp peak in the specific heat is seen at low temperature. The exact values of the transition temperatures still differ, but the general behavior is congruent. The results from Parsons and Williams do not show a sharp peak associated with the solid-liquid transition. The effects due to limited energy ranges may be part of the reason why Parsons and Williams do not see a "first order" peak. The coil-globule transition is likely also affected. Moreover, Parsons and Williams used a modified implementation of the Wang-Landau method in order to improve efficiency. However, it is still unclear whether this method has any errors associated with the calculation of the density of states. The results presented in this article were generated using the standard approach of the Wang-Landau method, and should provide a solid foundation for future work on larger chain lengths.

Acknowledgments

We thank T. Wüst, C. Gervais, S.-H. Tsai, and W. Paul for helpful discussions. This work is supported by DMR-0307082 and DMR-0341874.

Received on 17 October, 2007

• [1] D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, second edition (Cambridge University Press, New York, NY, 2005).
• [2] L. Yelash, M. Müller, W. Paul, and K. Binder, J. Chem. Theory Comput. 2, 588 (2006).
• [3] F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2000).
• [4] F. Wang and D. P. Landau, Phys. Rev. E 64, 056101 (2001).
• [5] D. P. Landau and F. Wang, Braz. J. Phys. 34, 354 (2004).
• [6] F. Rampf, W. Paul, and K. Binder, Europhys. Lett. 70, 628 (2005).
• [7] F. Rampf, K. Binder, and W. Paul, J. Polym. Sci. Pol. Phys. 44, 2542 (2006).
• [8] W. Paul, T. Strauch, F. Rampf, and K. Binder, Phys. Rev. E 75, 060801 (2007).
• [9] I. Carmesin and K. Kremer, Macromolecules 21, 2819 (1988).
• [10] D. F. Parsons and D. R. M. Williams, J. Chem. Phys. 124, 221103 (2006).
• [11] D. F. Parsons and D. R. M. Williams, Phys. Rev. E 74, 041804 (2006).
• [12] B. M. Baysal and F. E. Karasz, Macromol. Theory Simul. 12, 627 (2003).
• [13] A. Milchev, A. Bhattacharya, and K. Binder, Macromolecules 34, 1881 (2000).
• [14] D. P. Landau, S.-H. Tsai, and M. Exler, Am. J. Phys. 72, 1294 (2004).
• [15] C. Zhou, T. C. Schulthess, S. Torbr¨ugge, and D. P. Landau, Phys. Rev. Lett. 96, 120201 (2006).
• [16] P. Poulain, F. Calvo, R. Antoine, M. Broyer, and P. Dugourd, Phys. Rev. E 73, 056704 (2006).
• [17] A. Malakis, A. Peratzakis, and N. G. Fytas, Phys. Rev. E 70, 066128 (2004).
• [18] D. T. Seaton, S. J. Mitchell, and D. P. Landau, Braz. J. Phys. 36, 623 (2006).

Publication Dates

• Publication in this collection
  27 Mar 2008
• Date of issue
  Mar 2008

History

• Received
  17 Oct 2007