A New Monte Carlo Method and Its Implications for Generalized Cluster
  Algorithms
A New Monte Carlo Method and Its Implications for Generalized Cluster Algorithms
C. H. Mak and Arun K. Sharma
Department of Chemistry, University of Southern California,
Los Angeles, California 90089-0482, USA
(Dated: November 13, 2018)
We describe a novel switching algorithm based on a “reverse” Monte Carlo method, in which the
potential is stochastically modified before the system configuration is moved. This new algorithm
facilitates a generalized formulation of cluster-type Monte Carlo methods, and the generalization
makes it possible to derive cluster algorithms for systems with both discrete and continuous de-
grees of freedom. The roughening transition in the sine-Gordon model has been studied with this
method, and high-accuracy simulations for system sizes up to 10242 were carried out to examine the
logarithmic divergence of the surface roughness above the transition temperature, revealing clear
evidence for universal scaling of the Kosterlitz-Thouless type.
PACS numbers: 05.10.Ln, 05.50.+q, 64.60.Ht, 75.40.Mg
Large-scale Monte Carlo (MC) simulations are often
plagued by slow sampling problems. These problems are
especially severe in systems near the critical point or in
those with strong correlations. Slow sampling problems
manifest themselves as poor scaling of the dynamic relax-
ation time with the system size, making large-size sim-
ulations extremely slow to converge. The cause of these
problems is that most MC simulations are based on lo-
cal moves, and when the correlation length of the system
grows or as relaxation modes of the system become heav-
ily entangled, local moves become increasingly inefficient.
But if nonlocal MC moves are used [1], their acceptance
ratios are often found to be exceedingly low when system
correlations are strong.
One way to circumvent these problems was suggested
by Swendsen and Wang [2], who devised a clever scheme
where large-scale nonlocal MC moves may be constructed
to achieve high sampling efficiencies by exploiting certain
geometric symmetries in the system. This algorithm led
to a marked reduction in the dynamical scaling exponent
in the 2-dimensional Ising model near criticality. Since
the nonlocal moves in this algorithm update a large num-
ber of degrees of freedom at the same time, the Swendsen-
Wang method and others inspired by it are also often
referred to as “cluster Monte Carlo” methods.
Since Swendsen and Wang’s paper in 1987, many
cluster-type MC algorithms have appeared [3, 4, 5, 6,
7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. But
the success of cluster MC has not been universal because
the proper cluster moves needed seem to be highly de-
pendent on the system, and efficient cluster MC meth-
ods have been found for only a small number of models
[3, 5, 6, 7, 11, 14, 15, 18, 19] so far. The difficulty of
formulating a generalized MC method that could work
for any system seems to be associated with the appar-
ent geometric nature of the cluster-type MC methods –
all existing cluster MC methods have been derived in
one way or another by using certain geometric features
of the system. For example, in the original Swendsen-
Wang formulation a mapping between the Ising model
and the percolation model originally described by For-
tuin and Kasteleyn [20] was exploited to effect cluster
spin flips. In other models, the requisite mapping is not
always obvious, making cluster MC methods difficult to
implement for general systems.
In this letter, we will show that the derivations of clus-
ter MC methods do not have to be based on geometric
features of the systems. Instead, they may be more con-
veniently formulated based on algebraic features of the
system potential V (C). We will exploit this algebraic for-
mulation and suggest a way to generalize cluster Monte
Carlo methods to systems with any potential.
We focus on classical systems with partition function
Z = Tr e−V (C), where V (C) is the potential in units of
the temperature T . Acceptable Monte Carlo methods
to sample the system configurations can be constructed
using any transition probabilities W (C → C′) as long as
the detailed balance condition
W (C → C′)e−V (C) = W (C′ → C)e−V (C
′) (1)
is satisfied. Conventional MC methods such as Metropo-
lis [21] accomplishes this in two steps: a trial move is
made from C to C′ with some transition probability, and
the move is then accepted or rejected according to an
acceptance probability based on V (C′), V (C) or both, so
that the composite process satisfies Eqn.(1). This way
of constructing the Markov chain – trial moves followed
by acceptance/rejection – has been the accepted “stan-
dard” method for doing MC since the inception of the MC
method [22]. Other MC methods do exist, such as the
heat bath algorithm [23], which follow alternative strate-
gies, but by far the standard method is conceptually the
simplest and most convenient in practice.
In the Monte Carlo method we are proposing, we will
reverse the order of the two steps in the standard method.
That is, we will first determine an acceptable way to mod-
ify the potential V and then find a transition C → C′ that
is consistent with the new potential. To our knowledge,
http://arxiv.org/abs/0704.1539v1
the basic elements of this “reverse MC” idea were first
suggested by Kandel et al. [4], who used it to stochasti-
cally remove interaction terms from the system’s poten-
tial in an Ising model to arrive at an alternative deriva-
tion of the Swendsen-Wang method. The formulation of
Kandel et al. was limited to discrete systems like the
Ising model. In the following, we will show how the re-
verse MC idea may be formulated more generally for any
system, discrete or continuous, and how it may then be
used as a framework to construct generalized cluster al-
gorithms.
Consider a system with potential V =
i vi+Vr, con-
sisting of a number of “interaction terms” vi plus a “resid-
ual” Vr . These interactions may be bonds between parti-
cles, interactions of the particles with a field, or any other
additive terms in V . We consider replacing each interac-
tion term vi by some pre-selected ṽi with a “switching”
probability
Si(C) = cie
∆vi(C), (2)
where ∆vi = vi− ṽi, ci = e
i and ∆v∗i = maxC vi(C).
The outcome of the switches defines two complemen-
tary sets of interactions – the switched ones σ̃ and the
unswitched ones σ̄. Using the outcome of the switches,
we define a stochastically modified potential Ṽ as follows:
ṽi +
v̄j + Vr, (3)
with v̄i = vi − ln(1 − Si). An MC pass starts with an
attempt to switch every interaction vi to the new ṽi us-
ing the Si defined in Eqn.(2). If the switch is successful,
the interaction is replaced by ṽi. If not, the interaction
is replaced by another interaction v̄i. This is followed
by an update in the configuration of the entire system
from C to C′ using a transition probability W̃ (C → C′)
that satisfies detailed balance on the modified potential
Ṽ . This constitutes one pass. The move from C to C′ can
of course be carried out using any conventional MC move
that satisfies detailed balance on the modified potential.
But the reverse formulation of the MC method now offers
possibilities that were previously unavailable to conven-
tional MC methods — if a simple scheme can be devised
to update the configuration of the entire system on the
stochastically modified potential, one can envision de-
signing global moves for the system to accelerate its sam-
pling, and our freedom in choosing the ṽi can be actively
exploited to facilitate this. Within this context, the origi-
nal formulation of Kandel et al. corresponds to switching
vi to ṽi = 0, i.e. simplifying the potential by deleting in-
teractions from it. Kandel et al. showed that for the Ising
model they could easily construct global moves on this
stochastically simplified potential and their formulation
regenerates the Swendsen-Wang method. But compared
to the deletion formulation of Kandel et al., the switch-
ing implementation of the reverse MC method now offers
a much wider set of possibilities because the form of the
“switch to” interactions is completely arbitrary. Whereas
previously there may not be an obvious way to globally
update the configuration of the system on the original po-
tential, with the proper choices for ṽi large-scale moves
may now become possible on the stochastically modified
potential. Indeed, we have shown that the switching idea
may be used to formulate a cluster MC algorithm for a
Lennard-Jones fluid [24].
Equations (2), (3) and the transition probability W̃
define the switching algorithm. To prove detailed bal-
ance Eqn.(1) for the switching algorithm, it is suffi-
cient to treat a case where there are only two interac-
tion terms. Extension to any number of interactions
is straightforward. Starting with C, with two interac-
tion terms v1 and v2, there are four possible outcomes
from the switch: I. both 1 and 2 are switched, which oc-
curs with probability PI = S1(C)S2(C), II. 1 is switched
and 2 is unswitched, with PII = S1(C)[1 − S2(C)], III.
1 is unswitched and 2 is switched, with PIII = [1 −
S1(C)]S2(C), and IV. both 1 and 2 are unswitched, with
PIV = [1 − S1(C)][1 − S2(C)]. After the switch, an up-
date C → C′ is made with a transition probability W̃
that satisfies detailed balance on the modified potential
Ṽ defined in Eqn.(3). Each of the four channels will
have a different W̃ : W̃I, W̃II, etc., and W (C → C
′) in
Eqn.(1) is the sum PIW̃I + PIIW̃II + PIIIW̃III + PIVW̃IV
over all four channels. For the reverse transition, we
start with C′ and consider switching v1(C
′) → ṽ1(C
′) and
′) → ṽ2(C
′). Again there are four possible outcomes
and we call these scenarios I′, II′, III′ and IV′ as for
the forward transition. W (C′ → C) in Eqn.(1) is again
the sum PI′W̃I′ +PII′W̃II′ +PIII′W̃III′ +PIV′W̃IV′ . Using
the choice of S and Ṽ in Eqs.(2) and (3), it is easy to
show that detailed balance is obeyed along each chan-
nel, i.e. PIW̃I = PI′W̃I′ , PIIW̃II = PII′W̃II′ , etc. Of
course, detailed balance only requires the total W to sat-
isfy Eqn.(1), and it is possible to choose alternate forms
of S and Ṽ to do that, which may provide further flexi-
bilities.
In the rest of this letter, we will illustrate the effec-
tiveness of the switching implementation of the reverse
MC method, and show how it can be used to easily de-
rive a cluster MC method in a system with continuous
degrees of freedom. Previously, it has been extremely dif-
ficult to design cluster MC algorithms for systems with
continuous degrees of freedom. The few that have been
reported to date [3, 5, 6, 7, 11, 14, 19] were mainly based
on embedding discrete degrees of freedom into continu-
ous ones. The only exception is the recent discovery of a
geometric MC algorithm by Liu and Luitjen [19] where
they formulated a rejection-free MC method to sample
the Lennard-Jones fluid at its critical point.
The switching algorithm we have proposed makes the
process of deriving cluster-type MC methods much more
straightforward compared to those based on geometric
features of the system. We will illustrate this using
the sine-Gordon model, which can be used to study the
roughening transition on 2-dimensional surfaces. The
sine-Gordon (SG) model has the potential
VSG = T
〈i,j〉
|φi − φj |
cos(φi)
 , (4)
where φi are continuous variables on a 2-dimensional
square lattice, the second sum is over all sites and the
first sum is over all nearest-neighbor pairs. The SG
model is often considered to be a coarse-grained ver-
sion of the discrete Gaussian (DG) model with poten-
tial VDG = T
〈i,j〉 |hi − hj |
2, where hi are integers.
The DG model can in turn be mapped directly onto
the Coulomb gas model [25], and as a result, the SG
model should belong in the same universality class as the
Kosterlitz-Thouless (KT) transition [26, 27].
Roughening is expected to be a weak transition. The
only easily discernible divergence is exhibited in a loga-
rithmic dependence of the surface roughness σ2 = 〈|φi −
〈φ〉|2〉 on the system size L at the roughening tempera-
ture TR. Below TR, σ
2 is expected to approach a finite
value as L → ∞. In addition to this, since the divergence
is slow, large lattice sizes are needed to reach the scaling
limit. All of these features of the SG model make it hard
to accurately study the roughening transition using MC
simulations. Previous simulations have been limited to
small systems [14, 28, 29, 30, 31, 32].
In order to locate TR and study the scaling behavior
at the roughening transition, we make use of the switch-
ing algorithm of the reverse MC method proposed above.
The essential difficulty in treating the SG model is due
to the nonlinear cosine terms in the potential in Eqn.(4).
If these nonlinearities could be removed, the residual po-
tential becomes a simple Gaussian and we could move
the system configuration effectively using uncoupled sur-
face modes. With this in mind, we separate the potential
into two parts and treat the cosine terms as interactions
vi = −T
−1 cosφi and the harmonic part as the resid-
ual Vr. Each of the interactions is switched to a uni-
form potential ṽi = −T
−1 with Si = e
[1−cosφi]/T . After
the switches, a number of φi would have effectively lost
their couplings to the cosine potential, while the rest have
their interactions with the cosine potential replaced by
v̄i = − ln[e
cosφi/T − e−1/T ]. In the ensuring MC move,
we can update the unswitched φi which are now cou-
pled to the replacement interactions v̄i using conventional
methods, but try to formulate an update scheme where
the rest of the φi, now forming a constrained Gaussian
field, may be updated globally. A Gaussian field sub-
ject to linear constraints is still Gaussian, and in prin-
ciple we can diagonalize the potential to obtain all the
normal modes and then move each one independently.
This problem is the subject of fracton dynamics and has
been studied previously [33]. However, the cost of ob-
64 128 256 512 1024
64 128 256 512 1024
(b)(a)
FIG. 1: (a) Surface thickness σ2 as a function of the log
of lattice size L for different temperatures T . (b) Expanded
view of (a) for several temperatures near TR shifted vertically
to coincide at L = 64. Dashed line is the expected KT slope
at TR, showing that TR is slightly above T = 25 but below
taining all the normal modes of the constrained surface
and their frequencies will grow rapidly with the size of
the lattice and will only be feasible for small-size simula-
tions. Since the scaling limit in the SG model can only
be reached with large system sizes, we will need an al-
ternative method. The method we have used to update
the constrained Gaussian fields is based on the method
of Hoffman and Ribak [34]. Since the statistics of the
fluctuations of a Gaussian field from its mean is indepen-
dent of the value of the mean field, the fluctuations from
a free Gaussian field can be transferred to a constrained
field with a different mean. Near the roughening tem-
perature, the switching procedure produces roughly 5%
unswitched field points, and the corresponding mean field
with these constraints can easily be determined using a
steepest descent molecular dynamics method. To ensure
ergodicity, a conventional Metropolis move is also carried
out with every reverse MC move.
Figure 1(a) shows simulation results for the scaling of
the surface roughness σ2 with the length L of the lattice
in simulations with different lattice sizes L2 up to 10242
and at several temperatures T from 16 to 30. KT theory
[26, 27] predicts a logarithmic divergence for σ2 with a
universal slope at TR
2(L) = σ20(TR) +
lnL, (5)
where a is the lattice constant of the surface, and in the
units of Eqn.(4), a = 2π. Therefore, at TR the slope of
Fig. 1(a) should be equal to 4. Above TR, the logarith-
mic behavior of σ2 continues to hold except the constant
σ20 as well as the slope both increase with T . The data
in Fig. 1(a) show that for T = 21 and below, σ2 appears
to approach a finite value as L → ∞. Therefore, it is
64 128 256 512 1024
10000
Metropolis
cluster MC
ξ = 1.4
ξ = 2.5
FIG. 2: Dynamic scaling for the relaxation time τ (in units of
MC passes) of σ2 as a function of lattice size L in Metropolis
versus cluster MC, with their corresponding exponent ξ.
clear that TR > 21. The most recent simulation of the
SG model by Sanchez et al. [32] (referred to as the “or-
dered SG model”, OSGM, in this paper) suggested that
TR ≈ 16. Our data show that this is incorrect, and their
error is likely due to slow sampling problems. Locating
the precise value of TR is more involved, since the data
for T > 21 show no obvious tendency toward a finite σ2.
There are two possibilities: either these temperatures are
above TR or the system size may not be large enough to
have reached the scaling limit for these temperatures.
To determine which one is the case, we must resort to
a comparison between the simulation data with KT the-
ory. Figure 1(b) shows an expanded view of Fig. 1(a)
for a few temperatures 23 < T < 30, but for each T
the curve has been shifted vertically to remove the off-
set σ20 so that they all coincide at L = 64. The heavy
dashed line indicates the KT slope at TR according to
Eqn.(5). The data therefore suggest that TR is slightly
larger than 25 but less than 26, which is consistent with
the RG prediction for TR = 8π in the continuum model
[35, 36]. The apparent lack of an asymptotic σ2 in the
data for 21 < T < TR implies that even for L = 1024,
these lattice sizes are not yet large enough to be in the
scaling limit for those temperatures. Finally, to compare
the dynamic scaling behavior of the switching algorithm
with Metropolis, Fig. 2 shows the relaxation time in the
measurement of σ2 with the lattice size L slightly above
TR. Compared with the dynamic exponent ξ ≈ 2.5 in
Metropolis, the switching algorithm shows a markedly
improved ξ ≈ 1.4.
This work was supported by the National Science
Foundation under grant CHE-9970766.
[1] See, for example: K.E. Schmidt, Phys. Rev. Lett. 51,
2175 (1983); G.G. Batrouni, G.R. Katz, A.S. Kronfield,
G.P. Lepage, B. Svetitsky and K.G. Wilson, Phys. Rev. D
32, 2736 (1985); E. Dagotto and J.B. Kogut, Phys. Rev.
Lett. 58, 299 (1987); and J. Goodman and A.D. Sokal,
Phys. Rev. Lett. 56, 1015 (1986).
[2] R.H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 58, 86
(1987).
[3] R.G. Edwards and A.D. Sokal, Phys. Rev. D 38, 2009
(1988).
[4] D. Kandel, E. Domany, D. Ron, A. Brandt and E. Loh,
Phys. Rev. Lett. 60, 1591 (1988).
[5] U. Wolff, Phys. Rev. Lett. 60, 1461 (1988).
[6] F. Niedermayer, Phys. Rev. Lett. 61, 2026 (1988).
[7] U. Wolff, Phys. Rev. Lett. 62, 361 (1989).
[8] D. Kandel, R. Ben-Av and E. Domany, Phys. Rev. Lett.
65, 941 (1990).
[9] J.S. Wang and R.H. Swendsen, Physica A 167, 565
(1990).
[10] D. Kandel and E. Domany, Phys. Rev. B 43, 8539 (1991).
[11] H.G. Evertz, M. Hasenbusch, M. Marcu, K. Pinn and S.
Solomon, Phys. Lett. B 254, 185 (1991).
[12] S. Liang, Phys. Rev. Lett. 69, 2145 (1992).
[13] M. Hasenbusch, L. Gideon, M. Marcu and K. Pinn, Phys.
Rev. B 46, 10472 (1992).
[14] M. Hasenbusch, M. Marcu and K. Pinn, Physica A 211,
255 (1994).
[15] C. Dress and W. Krauth, J. Phys. A 28, L597 (1995).
[16] J. Houdayer, European Phys. J. B 22, 479 (2001).
[17] H.W.J. Blote, J.R. Heringa and E. Luijten, Comp. Phys.
Comm. 147, 58 (2002).
[18] W. Krauth and R. Moessner, Phys. Rev. B 67, 064503
(2003).
[19] J.W. Liu and E. Luijten, Phys. Rev. Lett. 92, 035504
(2004).
[20] C.M. Fortuin and P.W. Kasteleyn, Physica 57, 536
(1972).
[21] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H.
Teller and E. Teller, J. Chem. Phys. 21, 1087 (1953).
[22] N. Metropolis and S. Ulam, J. Am. Stat. Assoc. 44, 335
(1949).
[23] K. Binder, Monte Carlo Methods in Statistical Physics,
Topics of Current Physics, Vol. 7, (Springer, Berlin,
1986).
[24] C.H. Mak, J. Chem. Phys. 122, 214110 (2005).
[25] S.T. Chui and J.D. Weeks, Phys. Rev. B 14, 4978 (1976).
[26] J.M. Kosterlitz and D.J. Thouless, J. Phys. C 6, 1181
(1973).
[27] J.M. Kosterlitz, J. Phys. C 7, 1046 (1974).
[28] R.H. Swendsen, Phys. Rev. B 15, 5421 (1977).
[29] W.J. Shugard, J.D. Weeks and G.H. Gilmer, Phys. Rev.
Lett. 41, 1399 (1978).
[30] F. Falo, A.R. Bishop, P.S. Lomdahl and B. Horovitz,
Phys. Rev. B 43, 8081 (1991).
[31] A. Sanchez, D. Cai, N. Gronbech-Jensen, A.R. Bishop
and Z.J. Wang, Phys. Rev. B 51, 14664 (1995).
[32] A. Sanchez, A.R. Bishop and E. Moro, Phys. Rev. E 62,
3219 (2000).
[33] T. Nakayama, K. Yakubo and R.L. Orbach, Rev. Mod.
Phys. 66, 381 (1994).
[34] Y Hoffman and E. Ribak, Astrophys. J. 380, L5 (1991).
[35] D.J. Amit, Y.Y. Goldschmidt and S. Grinstein, J. Phys.
A 13, 585 (1980).
[36] H.J.F. Knops and L.W.J. den Ouden, Physica A 103,
597 (1980).
