- Email: [email protected]
Enhanced sampling in simulations of dense systems
Computer Physics Communications 146 (2002) 113–117 www.elsevier.com/locate/cpc
Enhanced sampling in simulations of dense systems W. Paul ∗ , M. Müller Institut für Physik, Johannes-Gutenberg-Universität, Staudingerweg 7, D-55099 Mainz, Germany
Abstract In the simulations of a variety of systems we encounter the problem of large relaxation times due to the dense packing of the systems constituents. We propose an algorithm to overcome this slowing down by temporarily allowing the constituents of a 3d systems to escape into a 4th space coordinate. The idea will be exemplified for the problem of a homopolymer collapse. 2002 Published by Elsevier Science B.V. PACS: 83.10.Nn; 02.70.Lq; 66.20.+d
1. Introduction In the computer simulation of condensed matter systems we are often faced with the problem of prohibitively large relaxation times of the systems due to the dense packing of their constituents. Because of their mutual excluded volume the motion of one particle requires the simultaneous cooperative rearrangement of many of its neighbors. This jamming or caging plays a central role in as diverse situations as rafting logs down a river or the structural glass transition, studied experimentally, for example, for hard-sphere colloids [1]. Another example where this problem occurs and which therefore has inspired the invention of new and efficient simulation methods is the coil-globule transition of polymers [2,3] or protein folding [4–6] in the biophysical area. In this situation and also for the equilibration of a dense polymer melt, the packing constraints get further amplified by the connectivity * Corresponding author.
E-mail address: [email protected] (W. Paul).
of the chains and the resulting topological interactions between them. All these problems have been with the simulation community for several decades already, but still pose a big challenge. For instance, the simulations of the structural glass transition upon cooling of model systems (Lennard–Jones mixtures, bead-spring polymers) were generally restricted to the slightly super-cooled temperature region due to equilibration problems [7, 8]. Although these are partly density related and partly temperature related, a method to reduce the slowing down due to dense packing would be very helpful. Also the equilibrium phase behavior of collapsed polymer chains below the Θ-point [9] is poorly understood and we will use this problem to exemplify our idea. In the next section we will explain the basic physical idea and how to cast it into an expanded ensemble Monte Carlo algorithm and introduce the bond-fluctuation model of polymers. Section 3 will then detail some results and Section 4 will present some conclusions and outlook.
0010-4655/02/$ – see front matter 2002 Published by Elsevier Science B.V. PII: S 0 0 1 0 - 4 6 5 5 ( 0 2 ) 0 0 4 4 2 - 3
114
W. Paul, M. Müller / Computer Physics Communications 146 (2002) 113–117
2. Physical idea and expanded ensemble realization The basic idea of our approach [10] can be stated very simply: if the high density of the system is the reason for the strong slowing down of the relaxation, one has to—at least temporarily—reduce the density in order to speed up equilibration. This is also the basic idea underlying Monte Carlo algorithms proposed in the literature to measure the chemical potential in dense systems [11,12], to study the protein folding problem [4] or simply to equilibrate dense polymer melts [13]. We propose to avoid the density slowing down by allowing all particles of a 3-dimensional system to move in 4 space dimensions. Instead of simulating a cubic box of linear size L with periodic boundary conditions, we can simulate a system S = L3 × [0, L4 ]. We will employ L4 = 1 in the following. Every state of our extended system with all particles having the same coordinate in the fourth direction is then a valid configuration of the physical system we want to study. Since due to entropic reasons the system will not on its own assemble in one three-dimensional cut of our extended ensemble, we have to drive the simulated system into this state using an external field in the fourth direction. The expanded system Hamiltonian is H = H0 +
N
hx4 (i),
(1)
i=1
where H0 is the Hamiltonian of the original model and x4 (i) is the coordinate of particle i in the fourth dimension. The strength of the external field, h, determines the distribution of the particles in the fourth dimension. We will completely ignore the excluded volume interaction of particles having different fourth coordinates and thus the effective density in each threedimensional cut of the four-dimensional simulation volume is strongly reduced and the density slowing down (glass transition) is avoided. This also reduces the topological interactions of the chains, which in turn helps to speed up the conformational equilibration.
We define an expanded ensemble [14–16] which we will treat in a Monte Carlo simulation through the following partition function 1 H0 + N i=1 hx4 (i) . (2) exp − Z= W (h) kB T h
{c}
Here {c} stands for the set of 4N coordinates of the particles and we assume a discrete set of values 0 = h1 < h2 < · · · , hm for the strength of the external field. The Monte Carlo procedure will contain random jumps hi → hi±1 . In order for these to be accepted with sufficient probability, the weight factors W (h) are determined self-consistently such that the h-value (ideally) performs a symmetric random walk between 0 and hm . If we would allow for swaps of complete configurations between different h-values, we could also realize this idea in the form of a parallel tempering algorithm [16]. In applications to other problems like, for instance, dense liquids or polymer melts, one would need to explore the freedom in the choice of Monte Carlo moves to find the most efficient realization of our idea. For example, in a simulation of a continuum model with Lennard–Jones interactions 12 6 σ σ (3) − U (r) = 4 r r r could mean the norm of the 4d inter-particle distance vector. The systems can then be simulated either by a combination of MD and MC propagation or by MD propagation alone, employing an external potential driving the particles into a 3-dimensional cut and using the parallel tempering approach. We will use the bond-fluctuation model for this work, where each repeat unit of a polymer chain is defined as a unit cube on a simple cubic lattice [17]. Consecutive monomers are connected by bonds from the set {[2, 0, 0], [2, 1, 0], [2, 1, 1], [2, 2, 1], [3, 0, 0], [3, 1, 0]}, where the square brackets stand for all vectors obtainable from the given ones through lattice symmetry operations. The choice of a maximum √ length of b = 10 together with the selection of allowed bond vectors ensures non-crossability of the chains in the course of random hopping movements of the repeat units. To induce a collapse transition into the model we introduce an attractive square well potential √ acting between all monomers which are less than 6 lattice
W. Paul, M. Müller / Computer Physics Communications 146 (2002) 113–117
units apart in space, i.e. on all those sitting on one of the positions {[2, 0, 0], [2, 1, 0], [2, 1, 1]} around a central one. The Hamiltonian of the model can be written as βH0 = −n
(4)
where n is the number of monomer pairs according to the above definition and β = 1/kB T . We employ neither excluded volume nor attractive interactions between monomers with x4 (i) = 0 and those with x4 (i) = 1. The expanded system Hamiltonian can therefore be written as βH = −n + hn4d ,
(5)
where n4d is the number of monomers with x4 (i) = 1 and where we have absorbed temperature into the definition of the interaction energy and the external potential h. The expanded ensemble partition function is 1 exp{n − hn4d }. (6) Z= W (h) h
{c}
Chain conformations are changed through random hops, where a randomly chosen monomer hops in a randomly chosen lattice direction—excluded volume and bond constraints allowing—and slithering snake moves, where one tries to attach a randomly chosen bond at one randomly chosen chain end and cuts the bond at the opposite end of the chain. Furthermore monomers are moved up and down the fourth dimension randomly, and also h is changed into neighboring values randomly. All the moves are, of course, subjected to a Metropolis acceptance criterion. Four Monte Carlo step (MCS) consist of three attempted slithering snake moves per chain, one attempted random hop and up/down move per monomer and one attempted h-change. The weight factors W (h) are determined iteratively from the requirement exp{n − hn4d } W (h) = Zh =
115
For the model we are studying here, the coilglobule transition has been detemined to be at Θ = 0.476 ± 0.005 (Θ = 2.1 ± 0.01) [2,18].
3. Results As a first step we analyze the efficiency of the proposed new algorithm in exploring the configuration space of a chain of length N = 256 at = 0.644, i.e. deep in its collapsed state. In Fig. 1 we show the autocorrelation function of the number of pair contacts in a comparison between a standard 3d simulation and the new proposed 4d algorithm. In the 3d case, the decorrelation on the time scale displayed comes only about by fluctuations in a single phase. The time scale for phase changes is about 107 MCS [10]. In Fig. 2 the probability distribution of the energy (or number of pair contacts) at coexistence is shown for chain lengths from N = 64, where the double peak structure is barely visible, to N = 512, where we have a deep minimum around the dividing nmin = 3.125 between the two phases. Up to N = 256 a simpler form of the algorithm, where the interaction strength is kept constant (h) = 0 and which therefore is not adapted to sample both peaks in the energy histogram efficiently, is still capable to provide an equilibrium sampling of the two energetically separated phases. Performing its random
{c}
= exp −F (h) − hn4d .
(7)
To efficiently sample the coexisting phases for collapsed long chains, we choose the weights W (h) to generate equal weights between the phases, so that we will work at c (N) [10].
Fig. 1. Autcorrelation function of the number of pair contacts per monomer using a standard 3d simulation algorithm (full line) and in an expanded ensemble simulation described in the text (dashed line). Both simulations are for a chain of length N = 256 and a temperature 35% below the Θ-temperature.
116
W. Paul, M. Müller / Computer Physics Communications 146 (2002) 113–117
Fig. 2. Probability distributions of the coordination number (internal energy) at the coexistence value of the interaction energy for the indicated chain lengths.
Fig. 3. Average radial density profile of liquid and solid globules for N = 512 at c .
walk up and down in h, the system is able to sample both phases of the three-dimensional configuration. We can describe this as a way of “moving around the barrier” between the liquid and the solid state using the extended degree of freedom, whereas the standard multicanonical reweighting would flatten the barrier between the two phases. The two coexisting phases are a liquid and a solid globule, respectively. This can be seen in Fig. 3 where we show the density as a function of distance from the center of gravity of the globule for the two phases. In the liquid phase the center of the globule assumes a density typical for melts in the bond-fluctuation model and there is a broad interfacial region to the bad solvent around the globule. In the
Fig. 4. Extrapolation of the finite chain length transition energies to the infinite chain length limit. Opaque circles are the simulation result and the filled square is the linear regression extrapolation to infinite chain length.
solid state the internal density is very close to the maximum possible value of one and there is a sharp interface to the solvent. The division between the two phases also holds for the large scale chain structure and the distributions for the radius of gyration in the two phases. The solid globules display a very sharp distribution corresponding to the sharp interface observed in the radial density profile in Fig. 3. The liquid globules are on average much larger with a much broader distribution. Both distributions are well separated, especially for the longer chains. From the phase transition points c (N) we can finally determine the phase transition coupling in the infinite chain length limit. Phenomenologically we find an effective scaling behavior (see Fig. 4) c (N) = c (∞) − aN −2/3 with c = 0.510 ± 0.006 and a = 5.4 ± 0.2. This corresponds to a transition temperature of Tc = 1.96 ± 0.01 which is roughly 7% below the Θ-temperature.
4. Summary and outlook We have shown that an expanded ensemble simulation algorithm in which we artificially allow the particles to temporarily explore a fourth space dimension is efficient in overcoming the density slowing down in the globular state of a collapsed homopolymer chain. In this way we could study the equilibrium phase be-
W. Paul, M. Müller / Computer Physics Communications 146 (2002) 113–117
havior of the globule and analyze a first order liquidsolid transition for collapsed polymer chains in the bond-fluctuation model. The solid phase of the globule has a density close to the maximum dense packing (ρ = 1) on the lattice. Slowing down due to dense packing also leads to the structural glass transition of liquids. In these liquids one envisions a particle as being trapped in a cage formed by it’s neighbors. The mode-coupling critical temperature in the slightly supercooled liquid marks the density/temperature below which activated processes are required for a particle to leave that cage [19]. The ability to equilibrate these systems below the mode-coupling critical temperature would mark a major advance in the computer simulation studies of the glass transition. We believe that the proposed algorithm will prove helpful also in this case, at least in the context of simulations of dense polymer melts [13].
Acknowledgements
We gratefully acknowledge support the INTAS program under grant number 97-0678, the the German Ministery for Education and Research under grant number 03N6015 and a generous grant of computer time at the Computer Center of the University of Mainz. We are grateful to K. Binder for helpful comments.
117
References [1] W. van Megen, T.C. Mortensen, S.R. Williams, J. Müller, Phys. Rev. E 58 (1998) 6073. [2] P. Grassberger, Phys. Rev. E 56 (1996) 3682. [3] Y. Zhou, C.K. Hall, M. Karplus, Phys. Rev. Lett. 77 (1996) 2822. [4] G. Chikenji, M. Kikuchi, Y. Iba, Phys. Rev. Lett. 83 (1999) 1886. [5] H. Frauenkorn, U. Bastolla, E. Gerstner, P. Grassberger, W. Nadler, Phys. Rev. Lett. 80 (1998) 3149. [6] Y. Okamoto, Recent Res. Devel. in Pure & Applied Chem. 2 (1998) 1; A. Mitsutake, Y. Sugita, Y. Okamoto, preprint, condmat/0012021. [7] W. Kob, J. Phys.: Condens. Matter 11 (1999) R85–R115. [8] K. Binder, J. Baschnagel, W. Paul, Prog. Polym. Sci., submitted. [9] Y. Zhou, M. Karplus, J.M. Wichert, C.K. Hall, J. Chem. Phys. 107 (1997) 10 691. [10] W. Paul, M. Müller, J. Chem. Phys. 115 (2001) 630. [11] M. Müller, W. Paul, J. Chem. Phys. 100 (1994) 719. [12] N.B. Wilding, M. Müller, J. Chem. Phys. 101 (1994) 4324. [13] A. Bunker, B. Dünweg, Phys. Rev. E 63 (2001) 16 701. [14] E. Marinari, G. Parisi, Europhys. Lett. 19 (1992) 451. [15] A.P. Lyubartsev, A.A. Martinovski, S.V. Shevkunov, P.N. Vorontsov-Velyaminov, J. Chem. Phys. 96 (1992) 1776. [16] J.J. de Pablo, F.A. Escobedo, Adv. Chem. Phys. 105 (1999) 337. [17] H.P. Deutsch, K. Binder, J. Chem. Phys. 94 (1991) 2294; H.P. Wittmann, K. Kremer, Comput. Phys. Commun. 61 (1990) 309; W. Paul, K. Binder, D.W. Heermann, K. Kremer, J. Phys. II (France) 1 (1991) 37. [18] N.B. Wilding, M. Müller, K. Binder, J. Chem. Phys. 105 (1996) 802. [19] W. Götze, J. Phys.: Condens. Matter 11 (1999) A1–A45.
Enhanced sampling in simulations of dense systems W. Paul ∗ , M. Müller Institut für Physik, Johannes-Gutenberg-Universität, Staudingerweg 7, D-55099 Mainz, Germany
Abstract In the simulations of a variety of systems we encounter the problem of large relaxation times due to the dense packing of the systems constituents. We propose an algorithm to overcome this slowing down by temporarily allowing the constituents of a 3d systems to escape into a 4th space coordinate. The idea will be exemplified for the problem of a homopolymer collapse. 2002 Published by Elsevier Science B.V. PACS: 83.10.Nn; 02.70.Lq; 66.20.+d
1. Introduction In the computer simulation of condensed matter systems we are often faced with the problem of prohibitively large relaxation times of the systems due to the dense packing of their constituents. Because of their mutual excluded volume the motion of one particle requires the simultaneous cooperative rearrangement of many of its neighbors. This jamming or caging plays a central role in as diverse situations as rafting logs down a river or the structural glass transition, studied experimentally, for example, for hard-sphere colloids [1]. Another example where this problem occurs and which therefore has inspired the invention of new and efficient simulation methods is the coil-globule transition of polymers [2,3] or protein folding [4–6] in the biophysical area. In this situation and also for the equilibration of a dense polymer melt, the packing constraints get further amplified by the connectivity * Corresponding author.
E-mail address: [email protected] (W. Paul).
of the chains and the resulting topological interactions between them. All these problems have been with the simulation community for several decades already, but still pose a big challenge. For instance, the simulations of the structural glass transition upon cooling of model systems (Lennard–Jones mixtures, bead-spring polymers) were generally restricted to the slightly super-cooled temperature region due to equilibration problems [7, 8]. Although these are partly density related and partly temperature related, a method to reduce the slowing down due to dense packing would be very helpful. Also the equilibrium phase behavior of collapsed polymer chains below the Θ-point [9] is poorly understood and we will use this problem to exemplify our idea. In the next section we will explain the basic physical idea and how to cast it into an expanded ensemble Monte Carlo algorithm and introduce the bond-fluctuation model of polymers. Section 3 will then detail some results and Section 4 will present some conclusions and outlook.
0010-4655/02/$ – see front matter 2002 Published by Elsevier Science B.V. PII: S 0 0 1 0 - 4 6 5 5 ( 0 2 ) 0 0 4 4 2 - 3
114
W. Paul, M. Müller / Computer Physics Communications 146 (2002) 113–117
2. Physical idea and expanded ensemble realization The basic idea of our approach [10] can be stated very simply: if the high density of the system is the reason for the strong slowing down of the relaxation, one has to—at least temporarily—reduce the density in order to speed up equilibration. This is also the basic idea underlying Monte Carlo algorithms proposed in the literature to measure the chemical potential in dense systems [11,12], to study the protein folding problem [4] or simply to equilibrate dense polymer melts [13]. We propose to avoid the density slowing down by allowing all particles of a 3-dimensional system to move in 4 space dimensions. Instead of simulating a cubic box of linear size L with periodic boundary conditions, we can simulate a system S = L3 × [0, L4 ]. We will employ L4 = 1 in the following. Every state of our extended system with all particles having the same coordinate in the fourth direction is then a valid configuration of the physical system we want to study. Since due to entropic reasons the system will not on its own assemble in one three-dimensional cut of our extended ensemble, we have to drive the simulated system into this state using an external field in the fourth direction. The expanded system Hamiltonian is H = H0 +
N
hx4 (i),
(1)
i=1
where H0 is the Hamiltonian of the original model and x4 (i) is the coordinate of particle i in the fourth dimension. The strength of the external field, h, determines the distribution of the particles in the fourth dimension. We will completely ignore the excluded volume interaction of particles having different fourth coordinates and thus the effective density in each threedimensional cut of the four-dimensional simulation volume is strongly reduced and the density slowing down (glass transition) is avoided. This also reduces the topological interactions of the chains, which in turn helps to speed up the conformational equilibration.
We define an expanded ensemble [14–16] which we will treat in a Monte Carlo simulation through the following partition function 1 H0 + N i=1 hx4 (i) . (2) exp − Z= W (h) kB T h
{c}
Here {c} stands for the set of 4N coordinates of the particles and we assume a discrete set of values 0 = h1 < h2 < · · · , hm for the strength of the external field. The Monte Carlo procedure will contain random jumps hi → hi±1 . In order for these to be accepted with sufficient probability, the weight factors W (h) are determined self-consistently such that the h-value (ideally) performs a symmetric random walk between 0 and hm . If we would allow for swaps of complete configurations between different h-values, we could also realize this idea in the form of a parallel tempering algorithm [16]. In applications to other problems like, for instance, dense liquids or polymer melts, one would need to explore the freedom in the choice of Monte Carlo moves to find the most efficient realization of our idea. For example, in a simulation of a continuum model with Lennard–Jones interactions 12 6 σ σ (3) − U (r) = 4 r r r could mean the norm of the 4d inter-particle distance vector. The systems can then be simulated either by a combination of MD and MC propagation or by MD propagation alone, employing an external potential driving the particles into a 3-dimensional cut and using the parallel tempering approach. We will use the bond-fluctuation model for this work, where each repeat unit of a polymer chain is defined as a unit cube on a simple cubic lattice [17]. Consecutive monomers are connected by bonds from the set {[2, 0, 0], [2, 1, 0], [2, 1, 1], [2, 2, 1], [3, 0, 0], [3, 1, 0]}, where the square brackets stand for all vectors obtainable from the given ones through lattice symmetry operations. The choice of a maximum √ length of b = 10 together with the selection of allowed bond vectors ensures non-crossability of the chains in the course of random hopping movements of the repeat units. To induce a collapse transition into the model we introduce an attractive square well potential √ acting between all monomers which are less than 6 lattice
W. Paul, M. Müller / Computer Physics Communications 146 (2002) 113–117
units apart in space, i.e. on all those sitting on one of the positions {[2, 0, 0], [2, 1, 0], [2, 1, 1]} around a central one. The Hamiltonian of the model can be written as βH0 = −n
(4)
where n is the number of monomer pairs according to the above definition and β = 1/kB T . We employ neither excluded volume nor attractive interactions between monomers with x4 (i) = 0 and those with x4 (i) = 1. The expanded system Hamiltonian can therefore be written as βH = −n + hn4d ,
(5)
where n4d is the number of monomers with x4 (i) = 1 and where we have absorbed temperature into the definition of the interaction energy and the external potential h. The expanded ensemble partition function is 1 exp{n − hn4d }. (6) Z= W (h) h
{c}
Chain conformations are changed through random hops, where a randomly chosen monomer hops in a randomly chosen lattice direction—excluded volume and bond constraints allowing—and slithering snake moves, where one tries to attach a randomly chosen bond at one randomly chosen chain end and cuts the bond at the opposite end of the chain. Furthermore monomers are moved up and down the fourth dimension randomly, and also h is changed into neighboring values randomly. All the moves are, of course, subjected to a Metropolis acceptance criterion. Four Monte Carlo step (MCS) consist of three attempted slithering snake moves per chain, one attempted random hop and up/down move per monomer and one attempted h-change. The weight factors W (h) are determined iteratively from the requirement exp{n − hn4d } W (h) = Zh =
115
For the model we are studying here, the coilglobule transition has been detemined to be at Θ = 0.476 ± 0.005 (Θ = 2.1 ± 0.01) [2,18].
3. Results As a first step we analyze the efficiency of the proposed new algorithm in exploring the configuration space of a chain of length N = 256 at = 0.644, i.e. deep in its collapsed state. In Fig. 1 we show the autocorrelation function of the number of pair contacts in a comparison between a standard 3d simulation and the new proposed 4d algorithm. In the 3d case, the decorrelation on the time scale displayed comes only about by fluctuations in a single phase. The time scale for phase changes is about 107 MCS [10]. In Fig. 2 the probability distribution of the energy (or number of pair contacts) at coexistence is shown for chain lengths from N = 64, where the double peak structure is barely visible, to N = 512, where we have a deep minimum around the dividing nmin = 3.125 between the two phases. Up to N = 256 a simpler form of the algorithm, where the interaction strength is kept constant (h) = 0 and which therefore is not adapted to sample both peaks in the energy histogram efficiently, is still capable to provide an equilibrium sampling of the two energetically separated phases. Performing its random
{c}
= exp −F (h) − hn4d .
(7)
To efficiently sample the coexisting phases for collapsed long chains, we choose the weights W (h) to generate equal weights between the phases, so that we will work at c (N) [10].
Fig. 1. Autcorrelation function of the number of pair contacts per monomer using a standard 3d simulation algorithm (full line) and in an expanded ensemble simulation described in the text (dashed line). Both simulations are for a chain of length N = 256 and a temperature 35% below the Θ-temperature.
116
W. Paul, M. Müller / Computer Physics Communications 146 (2002) 113–117
Fig. 2. Probability distributions of the coordination number (internal energy) at the coexistence value of the interaction energy for the indicated chain lengths.
Fig. 3. Average radial density profile of liquid and solid globules for N = 512 at c .
walk up and down in h, the system is able to sample both phases of the three-dimensional configuration. We can describe this as a way of “moving around the barrier” between the liquid and the solid state using the extended degree of freedom, whereas the standard multicanonical reweighting would flatten the barrier between the two phases. The two coexisting phases are a liquid and a solid globule, respectively. This can be seen in Fig. 3 where we show the density as a function of distance from the center of gravity of the globule for the two phases. In the liquid phase the center of the globule assumes a density typical for melts in the bond-fluctuation model and there is a broad interfacial region to the bad solvent around the globule. In the
Fig. 4. Extrapolation of the finite chain length transition energies to the infinite chain length limit. Opaque circles are the simulation result and the filled square is the linear regression extrapolation to infinite chain length.
solid state the internal density is very close to the maximum possible value of one and there is a sharp interface to the solvent. The division between the two phases also holds for the large scale chain structure and the distributions for the radius of gyration in the two phases. The solid globules display a very sharp distribution corresponding to the sharp interface observed in the radial density profile in Fig. 3. The liquid globules are on average much larger with a much broader distribution. Both distributions are well separated, especially for the longer chains. From the phase transition points c (N) we can finally determine the phase transition coupling in the infinite chain length limit. Phenomenologically we find an effective scaling behavior (see Fig. 4) c (N) = c (∞) − aN −2/3 with c = 0.510 ± 0.006 and a = 5.4 ± 0.2. This corresponds to a transition temperature of Tc = 1.96 ± 0.01 which is roughly 7% below the Θ-temperature.
4. Summary and outlook We have shown that an expanded ensemble simulation algorithm in which we artificially allow the particles to temporarily explore a fourth space dimension is efficient in overcoming the density slowing down in the globular state of a collapsed homopolymer chain. In this way we could study the equilibrium phase be-
W. Paul, M. Müller / Computer Physics Communications 146 (2002) 113–117
havior of the globule and analyze a first order liquidsolid transition for collapsed polymer chains in the bond-fluctuation model. The solid phase of the globule has a density close to the maximum dense packing (ρ = 1) on the lattice. Slowing down due to dense packing also leads to the structural glass transition of liquids. In these liquids one envisions a particle as being trapped in a cage formed by it’s neighbors. The mode-coupling critical temperature in the slightly supercooled liquid marks the density/temperature below which activated processes are required for a particle to leave that cage [19]. The ability to equilibrate these systems below the mode-coupling critical temperature would mark a major advance in the computer simulation studies of the glass transition. We believe that the proposed algorithm will prove helpful also in this case, at least in the context of simulations of dense polymer melts [13].
Acknowledgements
We gratefully acknowledge support the INTAS program under grant number 97-0678, the the German Ministery for Education and Research under grant number 03N6015 and a generous grant of computer time at the Computer Center of the University of Mainz. We are grateful to K. Binder for helpful comments.
117
References [1] W. van Megen, T.C. Mortensen, S.R. Williams, J. Müller, Phys. Rev. E 58 (1998) 6073. [2] P. Grassberger, Phys. Rev. E 56 (1996) 3682. [3] Y. Zhou, C.K. Hall, M. Karplus, Phys. Rev. Lett. 77 (1996) 2822. [4] G. Chikenji, M. Kikuchi, Y. Iba, Phys. Rev. Lett. 83 (1999) 1886. [5] H. Frauenkorn, U. Bastolla, E. Gerstner, P. Grassberger, W. Nadler, Phys. Rev. Lett. 80 (1998) 3149. [6] Y. Okamoto, Recent Res. Devel. in Pure & Applied Chem. 2 (1998) 1; A. Mitsutake, Y. Sugita, Y. Okamoto, preprint, condmat/0012021. [7] W. Kob, J. Phys.: Condens. Matter 11 (1999) R85–R115. [8] K. Binder, J. Baschnagel, W. Paul, Prog. Polym. Sci., submitted. [9] Y. Zhou, M. Karplus, J.M. Wichert, C.K. Hall, J. Chem. Phys. 107 (1997) 10 691. [10] W. Paul, M. Müller, J. Chem. Phys. 115 (2001) 630. [11] M. Müller, W. Paul, J. Chem. Phys. 100 (1994) 719. [12] N.B. Wilding, M. Müller, J. Chem. Phys. 101 (1994) 4324. [13] A. Bunker, B. Dünweg, Phys. Rev. E 63 (2001) 16 701. [14] E. Marinari, G. Parisi, Europhys. Lett. 19 (1992) 451. [15] A.P. Lyubartsev, A.A. Martinovski, S.V. Shevkunov, P.N. Vorontsov-Velyaminov, J. Chem. Phys. 96 (1992) 1776. [16] J.J. de Pablo, F.A. Escobedo, Adv. Chem. Phys. 105 (1999) 337. [17] H.P. Deutsch, K. Binder, J. Chem. Phys. 94 (1991) 2294; H.P. Wittmann, K. Kremer, Comput. Phys. Commun. 61 (1990) 309; W. Paul, K. Binder, D.W. Heermann, K. Kremer, J. Phys. II (France) 1 (1991) 37. [18] N.B. Wilding, M. Müller, K. Binder, J. Chem. Phys. 105 (1996) 802. [19] W. Götze, J. Phys.: Condens. Matter 11 (1999) A1–A45.