Heuristic Monte Carlo Method Applied to Cooperative Motion
Algorithm for Binary Lattice Fluid
P. Knychala1, ∗ and M. Banaszak2

arXiv:1410.3015v1 [cond-mat.stat-mech] 11 Oct 2014

1

Chair of Computer Science, The President Stanislaw

Wojciechowski Higher Vocational State School in Kalisz
ul. Nowy Swiat 4, 62-800 Kalisz, Poland
2

Faculty of Physics, A. Mickiewicz University

ul. Umultowska 85, 61-614 Poznan, Poland
(Dated: March 1, 2018)

Abstract
The Cooperative Motion Algorithm is an efficient lattice method to simulate dense polymer
systems and is often used with two different criteria to generate a Markov chain in the configuration
space. While the first method is the well-established Metropolis algorithm, the other one is an
heuristic algorithm which needs justification. As an introductory step towards justification for the
3D lattice polymers, we study a simple system which is the binary equimolar fluid on a 2D triangular
lattice. Since all lattice sites are occupied only selected type of motions are considered, such
the vacancy movements, swapping neighboring lattice sites (Kawasaki dynamics) and cooperative
loops. We compare both methods, calculating the energy as well as heat capacity as a function
of temperature. The critical temperature, which was determined using the Binder cumulant,
was the same for all methods with the simulation accuracy and in agreement with the exact
critical temperature for the Ising model on the 2D triangular lattice. In order to achieve reliable
results at low temperatures we employ the parallel tempering algorithm which enables simultaneous
simulations of replicas of the system in a wide range of temperatures.
Keywords: Cooperative Motion Algorithm, Heuristic Monte Carlo

∗

p.knychala@pwsz.kalisz.pl

1

I.

INTRODUCTION

Computer simulations are one of the most important research tools in modern science [1].
In particular, Monte Carlo method, which generates a Markov chain in the space of states,
is used in physics and chemistry problems since 1953 [2]. In this algorithm, various states of
system are generated in such a way and the next state depends exclusively on the previous
state. A random walk in the state space can be generated in order to achieve a given density
of states in such a manner that some new states are accepted and the other ones rejected.
One of the most celebrated examples (of using a lattice Monte Carlo method) is the Ising
model which was developed for magnetic systems [3, 4]. In this model each spin can be
directed up or down and the number of the up and down spins is not preserved, so the order
parameter (total magnetization) is not conserved. The Ising model with conserved order
parameter can also be employed to study the separations in liquid mixtures [5, 6].
Monte Carlo method is also used to investigate properties of polymers. Those materials
are increasingly used in many areas of life and have the fascinating ability to self assemble
into periodic nanostructures. Polymer systems consist of many macromolecules and each
one is composed of numerous repeat units, that are called monomers or segments. Computer
simulations of such large system are not trivial so that different optimizations and simplifications are used. For example, in coarse-grained methods details of the chemical structure
of polymer are neglected and a few chemical monomers are replaced by single segment.
Another simplification used in computer simulations of polymer is the use of lattice, where
most floating-point calculations are eliminated.
A simulation of dense lattice polymer models is a considerable challenge because the lattice sites are completely filled with chain segments and two segments cannot occupy the same
space, and therefore it is not a straightforward task to generate a Markov chain in the configuration space. For lattice fluid, the Kawasaki dynamics [7] can be used, in which random
pairs of neighboring spins are swapped. This method, however, cannot be easily extended to
polymer systems, because the linear connectivity is usually violated by random swapping.
An alternative approach is the vacancy dynamics in which the Monte Carlo movements are
driven by diffusion of a single vacancy [5]. There are also other approaches, using more
than one vacancy, proposed by Yaldram and Binder [8, 9]. Finally, there are methods which
use cooperative loops, either in small incremental steps in a single cooperative loop. For
2

lattice polymer, Pakula developed two methods (employing cooperative loops, one way or
the other) to generate the Monte Carlo movements: Cooperative Motion Algorithm [10–12]
(CMA) and Dynamic Lattice liquid (DLL) Algorithm [13].
The CMA seems to be the most efficient method to simulate dense polymer system. Representative states of system are generated by cooperative movements of segments along the
random walk loops on the face-centered cubic lattice (FCC). The FCC lattice is very suitable for simulation of dense polymer system since the coordination number, z, of this lattice
is 12, which is similar to liquid where average number of nearest neighbors is about 12.
Acceptance of the new state of the system can be achieved by using well known Metropolis
rules [2] or using a heuristic method proposed by Pakula [14]. In the first case, the difference
between energy of new and previous states is considered, and in case of heuristic method only
the energy of new state is taken into account. The CMA method was used to investigated
various polymer system [15–27] but the comparison of the different criteria of state acceptance (the Metropolis and heuristic approach) were not systematically studied and described
in literature. For example, Vilgis and Weyersberg [15] applied CMA with the Metropolis
test to analyze phase behavior of symmetric diblock copolymer and their results were favorably compared with random-phase approximation. The same version of CMA was also
used to investigate static properties of a dendrimer with positively charged terminal groups
by Majtyka and Klos [24]. The CMA with heuristic method was employed to investigate
properties of copolymer system with various distributions of monomers along the chain [17].
Static and dynamic properties of diblock copolymer melts near the order - disorder transition were investigated using heuristic method in ref. [18]. Both above methods were used
to determine microphase separation in triblock copolymers [19] and in ion-containing block
copolymers [27].
In this paper, we compare two methods:
• standard Metropolis (with Kawasaki dynamics and cooperative loops)
• heuristic method (with cooperative loops)
The heuristic method was often used with CMA, but to best of our knowledge, has not
been systematically compared with the Metropolis algorithm. We use simple 2D triangular
lattice model of liquid consisting of two type of elements which is similar to the Ising model
3

(however the order parameter is conserved) instead of complex polymer system and seek to
answer the following questions:
1. How different criteria of new state acceptance and Monte Carlo movements (full loop
vs. loops divided into incremental steps) influence the critical temperature below
which the system separates into two phases?
2. How the above criteria influence the energy or heat capacity; are they quantitatively
consistent?

II.

SIMULATION METHODS

Since the main aim of this paper is to investigate algorithms, we used simple 2D triangular
lattice liquid instead of 3D lattice. The lattice of size L × L consists of L2 /2 sites and each
lattice site contains one of two types of segments: A or B. The volume fraction f of segment
A is equal to 0.5 and for segment B it is 1 − f . The number of each type of segments is
constant during simulations so that it is a model with conserved order parameter. An
attempt to move a single segment defines a single Monte Carlo step. There are no vacancies
on the lattice so that segments can move via cooperative motion in which a few elements
are moved at the same time, as shown in fig. 1. This single cooperative motion is called a
cooperative loop. Figure 1a presents a single loop before segments are moved and the final
state is presented in fig. 1b. In the athermal state (without interaction between segments)
this cooperative loop is created in the following way:
1. a segment from random lattice site is removed (temporary vacancy is created) and
this site is labeled both INITIAL and CURRENT,
2. one of six neighbors of CURRENT lattice site is selected at random,
3. segment from this selected site is moved to CURRENT site and selected site becomes
the CURRENT one,
4. points 2 and 3 are iterated until the random neighbor is not the same as the INITIAL
site,
5. the segment which was initially removed is put into the vacancy.
4

After athermal simulation, consisting of such cooperative loops, the segments are randomly
distributed on the lattice.
Next, interactions between segments are introduced and described through ij parameters.
The interaction between different type of segment AB is equal to  and between the same
type of segment is 0 (AA = BB = 0). The interaction are limited to nearest neighbors and
P
energy of single segment i is Ei∗ = 6j=1 ij ( parameter is also used as energy unit). The
reduced temperature T ∗ is kT /. It means that segments prefer segments of the same type
and below the critical temperature, Tc , the system separates in two phases, first is rich in
segments A and the second in segments B.
In case of this thermal simulation various states of system are generated using the same
cooperative motions (and also Kawasaki dynamics motions) and probability transition between current state X and new state X 0 is described by following formula:



∆E
0
.
P (X → X ) = min 1, exp −
kT

(1)

In this paper, three approaches employing the Metropolis criterion and also the heuristic
method are used.
In the first case, the loop is created using single vacancy which moves as a random walk.
The difference of energy is calculated for full loop and ∆E = EX 0 − EX , where EX is
energy of the initial state (state before cooperative loop) and EX 0 is energy of the state
after cooperative loop. We use standard Metropolis criteria [2], and refer to this method as
GM , the motions are non-local since more then one segment is moved in a single step of
simulation. In this case some mismatch with theory can occur due to possible large energy
changes involved with acceptance of long loops, as discussed in next chapter.
In the second case (LM ), the loop is divided into incremental steps in which the movement
of a single segment is analyzed. We refer to it as Metropolis with single vacancy. At the
beginning of the loop, the segment from random lattice site is removed without any condition
since it decrease energy of system (vacancy which does not interact with segments is created).
At the end of the loop, the removed segment is inserted with probability proportional to
Boltzmann factor since this motion always increases energy of system. In case of failure,
the whole loop is rejected. Inside the loop, EX describes energy of sliding segment at initial
position as shown in fig 2a and EX 0 describes energy of the same segment at new position
(fig. 2b). In a single step, we try to move only one segment.
5

Moreover, the Metropolis method was employed for the Kawasaki dynamics [28] (KA),
swapping neighboring lattice sites, which generated results that we use as a background for
other methods.
The last approach is similar to LM method with one exception, the acceptance condition
is always based on energy of single segment in new position (energy of previous state is not
included in Boltzmann factor). It means that ∆E = EX 0 . We call this method as heuristic
Monte Carlo method (HE) since the correctness of this approach is not obvious and compare
it to commonly used Metropolis criterion can be interesting and provide additional and
intriguing insights into CMA algorithm.
We use five different lattice size in presented simulations: 40×40, 50×50, 60×60, 80×80,
and 100 × 100. First of all, the athermal simulations is performed for at least 2 × 106 MCS.
Next, we perform thermal run from 3 × 106 to 6 × 106 MCS depending on lattice size and
the first half of run is needed to equilibrate the system, and the second half is used to collect
the data. All simulation experiments are repeated at least 3 times starting with different
initial athermal states.
In order to improve statistics (particularly at low temperatures) the parallel tempering [29, 30] (PT) algorithm was also used in simulations. In the PT method, M replicas of
system are simulated in parallel (in this case it was 48), each in different temperature Ti∗ ,
with i ranging from 1 to M . After a number of Monte Carlo steps (MCS) (in our case it
was 500 MCS), we try to exchange replicas with neighboring Ti∗ in random order with the
following probability:
P (Ri → Ri+1 ) = min [1, exp (−(βi − βi+1 )(Ei+1 − Ei ))]

(2)

where βi = 1/kT .

III.

RESULTS AND DISCUSSION

First, we carry out simulations for lattice size equal to 50×50. In fig. 3, energy per lattice
site, E ∗ , and heat capacity, Cv in terms of temperature for different methods are presented.
The heat capacity is calculated from the equation:
Cv =

(E ∗ − hE ∗ i)2
kT ∗2
6

(3)

One can observe single peak in Cv at 1/T ∗ = 0.55 (T ∗ = 1.82) for all Metropolis approaches
as well as heuristic method which is related to the separation of the system into two phases.
Above this critical temperatures, Tc , the segments of type A and B are uniformly distributed
on the lattice (like in fig 5b), while below Tc system start to separate into two phases (fig 5a).
The first one is rich in segments of type A and the second one mostly consists of segments B.
A significant decrease in energy is observed also at Tc , as shown in fig. 3b. This temperature
Tc is consistent with the theoretical value of critical temperature for the Ising model at 2D
triangular lattice (with coordination number z = 6) calculated by Houtappel in 1950 [31]
√
and equals to 1/ ln( 3) = 1.821.
In fig. 4a, we present probability of state acceptance, P (X), in term of loop length, l, for
the GM method. It is obvious that this probability decreases with loop length increasing
and for loop length greater than 200 is close to 0. We introduces additional parameter that
describes contribution of different loop lengths. It determines what length of the loop takes
most time of simulation and is described as:
C(l) = lP (l),

(4)

where P (l) is the probability of the loop of length l. This parameter is shown in fig. 4b.
The contribution decreases with loop length increasing and reaches a minimum at about
200, and then begins to rise sharply. The contribution of short loop length (shorter than 50)
is significant but the contribution of long loops (greater than 200) is more significant. As
shown in fig. 4a the probability of acceptance of long loops is close to 0 and this affects the
efficiency of GM method. Since the GM method is the least efficient, requiring the longest
time of simulation, we do not take this approach into account in further analysis.
Moreover, it is possible to determine of the critical temperature with greater precision by
calculating the Binder cumulant [32] (UL ):
h(m − 0.5)4 i
UL = 1 −
,
3 h(m − 0.5)2 i2
where m is the order parameter described by equation [33]:
Z 1
m=2
xA P (xA ) dx.

(5)

(6)

0.5

The P (xA ) is calculated from the histogram of the probability of finding a particle of a
given type in a given area. It can be calculated for different temperatures in terms of local
7

concentration xA of segments type A defined as NA /NR , where NA is the number of segment
type A in specified area and NR is the total number of segments in the same area. In contrast
to the Ising model where the number of spins is not preserved, in presented simulations the
volume fraction of segments of different types is constant so P (xA ) can not be calculated for
the whole lattice (in that case UL will be 0). In this paper, the radius of this area depends on
size of lattice and is equal to 0.3L as shown in fig. 5. Such approach was previously used and
described in ref. [34]. The value 0.3 was chosen arbitrary from a few selected tested values
because it provided the best statistics. Moreover, we observed that Tc does not depend on
this parameter.
The sample P (xA ) for different temperatures are presented in fig. 6a. For T > Tc the
P (xA ) has Gaussian distribution with maximum at xA = 0.5 and the peak is getting narrower
as the temperature is increased. This means that, statistically, in selected area the same
number of segment A and B is observed (fig. 6d). At Tc system starts to separate into
two phases as shown in fig. 6c. Below critical temperature P (xA ) has bimodal distribution
with two maxima at xA = 0 and 1 corresponding to two coexisting phases (fig. 6b). Those
maxima become narrower with decreasing temperature since the clusters grow and become
more homogeneous.
The value of the Binder cumulant depends on size of lattice for different values of temperatures except of Tc where UL tends towards an universal value and not depending on
L. For T > Tc segments are uniformly mixed (system is disordered) and UL → 0 while for
T < Tc system is separated into two phases (ordered) and UL → 2/3 as shown in fig. 7.
The simulations for different size of lattice were carried out to determine the Tc using
the Binder cumulant. We use four values of lattice size: 40 × 40, 60 × 60, 80 × 80, and
100 × 100 which are presented in fig. 8. Results of simulations for two Metropolis approaches
and heuristic method are presented in fig. 9-11. The summary of the results is presented
in table I. The critical temperatures presented in table I are determined on the basis of
intersection of the Binder cumulant calculated for different lattice size. We use 48 replicas,
different temperatures, in parallel tempering simulations in range form T ∗ = 1.62 to 2.02.
The distance between neighboring temperatures is 0.0085 so the simulation error (accuracy
of the simulation) is about 0.5% (0.0085/1.821). In the case of Kawasaki dynamics which
is used as the background for other methods, Tc = 1.816 and the error is 0.24%. For LM
step acceptance method the average Tc is 1.826 and the error is 0.31%. The Tc = 1.823 with
8

error equals to 0.16% is obtained for HE method. The critical temperature obtained for
heuristic method is consistent with both Metropolis approaches. In all cases the value of
error is below the accuracy of the simulation.

IV.

CONCLUDING REMARKS

We carried out Monte Carlo simulations of simple binary liquid using two methods:
• standard Metropolis (with different approaches),
• heuristic method.
The parameters determined during simulation with different approaches are qualitatively
and quantitatively consistent. The critical temperature obtained for all methods is the same
with the simulation accuracy. Moreover, the Metropolis with vacancy as well as heuristic
method (where loops are divided into incremental steps) are more efficient than standard
Metropolis method with full loop acceptance what is important for the simulation of dense
polymer system.

ACKNOWLEDGMENTS

P.K. and M.B. gratefully acknowledge the research grant from the Polish NCN No. DEC2012/07/N/ST4/00293 and the computational grant from the Supercomputing and Networking Center (PSNC) in Poznan, Poland.

9

[1] D. Frenkel and B. Smit, Understanding Molecular Simulation, Second Edition: From Algorithms to Applications (Computational Science) (Academic Press, 2001).
[2] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller, Journal of Chemical
Physics 21, 1087 (1953).
[3] K. Binder, Reports on Progress in Physics 60, 487 (1997).
[4] M. E. J. Newman and G. T. Barkema, Monte Carlo Methods in Statistical Physics (Oxford
University Press, 1999).
[5] P. Fratzl and O. Penrose, Physical Review B 50, 3477 (1994).
[6] K. Halagan and P. Polanowski, Journal of Non-Crystalline Solids 355, 1318 (2009).
[7] K. Kawasaki, T. Nagai, and K. Nakashima, Philosophical Magazine B: Physics of Condensed
Matter; Electronic, Optical and Magnetic Properties 60, 399 (1989).
[8] K. Yaldram and K. Binder, Acta Metallurgica Et Materialia 39, 707 (1991).
[9] K. Yaldram and K. Binder, Journal of Statistical Physics 62, 161 (1991).
[10] T. Pakula, Macromolecules 20, 679 (1987).
[11] T. Pakula and S. Geyler, Macromolecules 20, 2909 (1987).
[12] T. Pakula and S. Geyler, Macromolecules 21, 1665 (1988).
[13] T. Pakula and J. Teichmann, Materials Research Society Symposium - Proceedings 455, 211
(1997).
[14] A. Gauger, A. Weyersberg, and T. Pakula, Macromolecular Theory and Simulations 2, 531
(1993).
[15] A. Weyersberg and T. Vilgis, Physical Review E 48, 377 (1993).
[16] A. Gauger, A. Weyersberg, and T. Pakula, Die Makromolekulare Chemie, Theory and Simulations 2, 531 (1993).
[17] T. Pakula and K. Matyjaszewski, Macromolecular Theory and Simulations 5, 987 (1996).
[18] T. Pakula, K. Karatasos, S. H. Anastasiadis, and G. Fytas, Macromolecules 30, 8463 (1997).
[19] M. Banaszak, S. Woloszczuk, T. Pakula, and S. Jurga, Physical Review E 66, 031804 (2002).
[20] Y. Qin, H.-L. Liu, and Y. Hu, Molecular Simulation 29, 649 (2003).
[21] M. Majtyka and J. Klos, Journal of Physics Condensed Matter 18, 3581 (2006).
[22] J. Jeszka and T. Pakula, Polymer 47, 7289 (2006).

10

[23] P. Polanowski and J. Jeszka, Langmuir 23, 8678 (2007).
[24] M. Majtyka and J. Klos, Physical Chemistry Chemical Physics 9, 2284 (2007).
[25] A. Sikorski, P. Polanowski, P. Adamczyk, and S. Zerko, Journal of Molecular Modeling 17,
2209 (2011).
[26] P. Polanowski, E. Wawrzyska, and A. Sikorski, Macromolecular Theory and Simulations 22,
238 (2013).
[27] P. Knychala, M. Banaszak, and N. Balsara, Macromolecules 47, 2529 (2014).
[28] K. Kawasaki, Phase transitions and critical phenomena, vol. 2 (Academic Press, 1972).
[29] D. J. Earl and M. W. Deem, Physical Chemistry Chemical Physics 7, 3910 (2005).
[30] R. H. Swendsen and J. S. Wang, Physical Review Letters 57, 2607 (1986).
[31] R. Houtappel, Physica 16, 425 (1950).
[32] K. Binder and D. W. Heermann, Monte Carlo Simulation in Statistical Physics (Springer
Verlag, 1992).
[33] S. Das, J. Horbach, K. Binder, M. Fisher, and J. Sengers, Journal of Chemical Physics 125
(2006).
[34] W. Kwak and D. Landau, International Journal of Modern Physics C 17, 15 (2006).

11

intersection

KA

LM

HE

40 × 40 − 60 × 60

1.824 1.815 1.809

40 × 40 − 80 × 80

1.819 1.822 1.826

40 × 40 − 100 × 100 1.816 1.825 1.823
60 × 60 − 80 × 80

1.810 1.828 1.837

60 × 60 − 100 × 100 1.814 1.831 1.826
80 × 80 − 100 × 100 1.815 1.836 1.820
average

1.816 1.826 1.823

standard deviation 0.005 0.007 0.009
error

0.24% 0.31% 0.16%

TABLE I. Critical temperatures for different methods obtained as intersection points of the Binder
cumulant from different lattice sizes.

12

Figure captions:
Fig. 1: Example of single cooperative motion for triangular 2D lattice. Grey arrow show
direction of segment motion.
Fig. 2: Nearest neighbors using to calculate energy of sliding segment (green) before (a)
and after (b) move of a single segment.
Fig. 3: Results of simulations (averaged over 5 runs) for 50 × 50 lattice size and four
different method: (a) heat capacity, Cv , with error bars (b) energy per lattice site, E ∗ .
Fig. 4: Probability of loop acceptance, P (X), (a) and contribution of different length of
loop into simulation, C(l), (b) in term of loop length, l, for the GM method at T ∗ = 1.82.
Fig. 5: Size of area (equal to 0.3L) which is used in simulation to calculation the density
profiles.
Fig. 6: Density profiles for 9 different temperatures (a) and snapshots from simulation
for 50 × 50 lattice size at T ∗ = 1.0 (b), 1.82 (c), 3.5 (d).
Fig. 7: The sample Binder cumulant in terms of 1/T ∗ for lattice size 100 × 100.
Fig. 8: Four sizes of box used in simulations to calculate the Binder cumulant.
Fig. 9: The Binder cumulant (average over 3 runs) for the Kawasaki method and four
different size of box. The inset shows intersection of the Binder cumulant.
Fig. 10: The Binder cumulant (average over 3 runs) for the LM method and four different
size of box. The inset shows intersection of the Binder cumulant.
Fig. 11: The Binder cumulant (average over 3 runs) for the heuristic method and four
different size of box. The inset shows intersection of the Binder cumulant.

13

3
2

1

2

9

3

8

5
8

1

9

4

4
7

6
7

5
6

FIG. 1.

14

b)

a)
FIG. 2.

15

0.0012
KA
LM
GM
HE

0.0010

Cv

0.0008
0.0006

a)

0.0004
0.0002
0.0000
0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.4000
1.2000
1.0000

E*

0.8000
b)
0.6000
0.4000
0.2000
0.0000
0.2

0.3

0.4

0.5

0.6
1/T

FIG. 3.

16

*

0.7

0.8

0.9

1

P(X)

1.0
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.0

a)

1

10

100

1000

0.22
0.20

C(l)

0.18
0.16

b)

0.14
0.12
0.10
0.08
1

10

100
l

FIG. 4.

17

1000

b)

a)
FIG. 5.

18

0.05
1.00
1.28
1.60
1.82
1.93
2.26
2.58
2.91
3.50

P(xA)

0.04
0.03
0.02

a)

0.01
0
0

0.2

0.4

0.6

0.8

1

xA

b)

c)

FIG. 6.

19

d)

0.7
0.65
0.6

UL

0.55

disorder

order

0.5
0.45
0.4
0.35
0.3
0.25
0.45

0.5

0.55

0.6
1/T*

FIG. 7.

20

0.65

0.7

40x40 60x60

80x80
FIG. 8.

21

100x100

0.66
0.64
0.62

L=40
L=60
L=80
L=100

0.6
UL

0.65
0.58
0.56
0.6
0.54
0.52

0.55
0.54

0.5
0.48
0.52

0.525

0.53

0.535

0.54

0.545
1/T*

FIG. 9.

22

0.55

0.555

0.55
0.56

0.56
0.565

0.57

UL

0.66
L=40
0.64
L=60
0.62
L=80
L=100
0.6
0.58
0.56
0.54
0.52
0.5
0.48
0.46
0.52 0.525

0.65
0.6
0.55
0.54
0.53

0.535

0.54

0.545
1/T*

FIG. 10.

23

0.55

0.555

0.55
0.56

0.56
0.565

0.57

0.66
0.64
0.62

L=40
L=60
L=80
L=100

0.6
UL

0.65
0.58
0.56
0.6
0.54
0.52

0.55
0.54

0.5
0.48
0.52

0.525

0.53

0.535

0.54

0.545
1/T*

FIG. 11.

24

0.55

0.555

0.55
0.56

0.56
0.565

0.57

