Papers in Physics, vol. 7, art. 070001 (2015)

Received: 19 January 2015, Accepted: 25 February 2015
Edited by: C. S. O’Hern
Reviewed by: M. Pica Ciamarra, Nanyang Technological University, Singapore.
Licence: Creative Commons Attribution 3.0
DOI: http://dx.doi.org/10.4279/PIP.070001

www.papersinphysics.org

ISSN 1852-4249

Wang–Landau algorithm for entropic sampling of arch-based microstates
in the volume ensemble of static granular packings

D. Slobinsky,1, 2∗ Luis A. Pugnaloni1, 2†

We implement the Wang–Landau algorithm to sample with equal probabilities the static
configurations of a model granular system. The “non-interacting rigid arch model” used
is based on the description of static configurations by means of splitting the assembly
of grains into sets of stable arches. This technique allows us to build the entropy as a
function of the volume of the packing for large systems. We make a special note of the
details that have to be considered when defining the microstates and proposing the moves
for the correct sampling in these unusual models. We compare our results with previous
exact calculations of the model made at moderate system sizes. The technique opens a
new opportunity to calculate the entropy of more complex granular models.

I. Introduction

In the study of static packings there exists still a
lack of predictive capabilities of the available the-
ories. Assemblies of objects that pack (such as
grains) can generally sample such packed configura-
tions only by the external excitation of the system.
These packings can be built by repeating a given
packing protocol (e.g., homogeneous compression
or deposition under an external field against a con-
fining boundary) on an initial random configura-
tion. Also, a Markovian or non-Markovian series
can be constructed by exciting the system from
the previous packing configuration. To what extent
the series of packings obtained (using either type of
protocol) can be modeled without information on

∗E-mail: dslobinsky@frlp.utn.edu.ar
†E-mail: luis.pugnaloni@frlp.utn.edu.ar

1 Departamento de Ingenieŕıa Mecánica, Facultad Re-
gional La Plata, Universidad Tecnológica Nacional, Av.
60 Esq. 124, 1900 La Plata, Argentina.

2 Consejo Nacional de Investigaciones Cient́ıficas y
Técnicas (CONICET), Argentina.

the dynamics that drives the system to the packed
configuration is still uncertain. The main reason
for this is that the few statistical approaches that
attempt to do this are strongly hinder by the poor
current ability to generate such packed structures
without using a dynamics to build the packings.

One might expect that the packing fraction and
its fluctuations, among other properties, could be
obtained from basic statistics without resourcing
to a full molecular dynamic type of simulations
(also known as “discrete element method”, DEM).
Although these types of simulations are powerful
enough to predict the behaviour of most systems
that pack, it is desirable to find a description that
could neglect the detailed dynamics between con-
secutive packed configurations.

The use of the tools provided by the ensemble
theory of statistical mechanics in problems of gran-
ular matter is still limited. Although many studies
perform statistical analysis of configurations of a
granular sample obtained during careful prepara-
tions in the laboratory or in molecular dynamic-
type simulations, very rarely sampling in a partic-
ular statistical ensemble is carried out either via an

070001-1



Papers in Physics, vol. 7, art. 070001 (2015) / D. Slobinsky et al.

analytic or a computational calculation of a model
system. This has prevented a direct assessment as
to whether ensemble theory is appropriate to de-
scribe the behaviour of these peculiar systems.

One pioneering contribution to the topic is the
idea that granular systems at mechanical equi-
librium could be treated as ensemble members,
putting forward the conjecture that the mean val-
ues of measurable quantities could be calculated us-
ing statistical mechanics for these ensembles [1]. In
this scheme, each of the thermodynamic variables
finds a counterpart, the volume taking the place of
the energy, and the compactivity that of the tem-
perature, amongst other transformations. However
beautiful this description may seem, the computa-
tional challenges to generate ensemble samples in
this context are extraordinaries [2–5].

One complication that prevents to a large extent
the use of the machinery of statistical mechanics in
this case is the fact that configurations, unlike in
traditional liquid theories, have to be checked for
the constraint of mechanical equilibrium. In a pre-
vious work [6], we have made a proposal on how to
deal with this, at least in a first approximation, by
describing the excitations of static granular systems
under gravity in terms of its arches. Since arches
are sets of grains that stabilize each other, these
are the basic units of mechanically stable struc-
tures in the packing. Any static configuration can
be described in terms of the arches formed by its
grains, their arch shape, position, orientations, etc.
We have considered, as an example, a model of a
two-dimensional (2D) granular system composed of
disks where arches are assumed to take a single pos-
sible structure and the arch–arch interactions (due
to the interlocking of arches) is neglected. We cal-
culated the exact entropy of this model (the non-
interacting rigid arch model, NIRA) by construct-
ing all possible configurations for moderate system
sizes. Of course, generating each state is a cumber-
some task if the system size is increased or if the
number of degrees of freedoms (DoF) is increased
by using more realistic models. Therefore, an alter-
native approach based on sampling the phase space
with the desired probability is necessary.

In the present work, we will calculate the entropy
of the NIRA model in the microcanonical ensem-
ble [7] using entropic sampling through the Wang–
Landau (WL) algorithm [8–11]. This approach al-
lows us to obtain the entire entropy function for all

possible volumes of the system in one single simu-
lation for larger systems and potentially for more
complex models. All derived properties, such as
compactivity and volume fluctuations, can then be
calculated through numerical differentiation. We
pay particular attention to the different descrip-
tions that can be realized for the NIRA model.
Some of these representations do not provide a di-
rect way of sampling the configuration space uni-
formly.

This work is organized as follows: in section II.,
we will review the WL algorithm. In section III.,
we will review the NIRA model and discuss differ-
ent ways of representing it, along with the issues
related to uniform sampling of the configurations.
We then present a representation that allows very
fast calculations of the entropy and we compare
the results with the exact counting of all config-
urations for systems of moderate size. Finally, we
discuss future directions to refine the arch-based en-
semble volume function towards capturing detailed
features of more realistic systems.

II. Wang-Landau algorithm

The WL method has revolutionized computational
statistical mechanics [8–11]. WL is a pure statisti-
cal method that can retrieve the density of states
(DoS) (hence the entropy) over a bounded region of
the energy spectrum from the sole knowledge of the
energy function. The spectacular computational
performance achieved by this method stems from
the fact that it presents no limitations for the sys-
tem to tunnel between potential barriers, in stark
contrast with classical Monte Carlo methods that
underperform when they encounter deep valleys in
the energy landscape.

WL finds the entropy of the system by means of
a Markov chain in the energy landscape which is
conveniently biased towards the less probable ener-
gies in a strongly history dependent manner. This
is achieved by using the multicanonical approach
[12] in which each possible configuration is sampled
with a probability given by the inverse of the den-
sity of states for its given energy. Specifically, WL
aims to obtain a flat histogram of visited energies
E by forcing the system to go through all config-
urations with a probability which is inverse to the
previous occurrence of that energy in the Markov

070001-2



Papers in Physics, vol. 7, art. 070001 (2015) / D. Slobinsky et al.

chain. The method is ergodic and asymptotically
fulfils the detailed balance condition [13]. There ex-
ists extensive literature on the WL algorithm but
here, we only summarize its most relevant steps. In
the following sections, we will refer to a configura-
tion of the system as a fixed set of values of all its
DoF.

In WL, one defines two histograms that are con-
tinuously updated as the Markov chain proceeds.
These histograms are the entropy S(E), which is
the output of the algorithm, and a control his-
togram H(E). After initializing H(E) = 0, S(E) =
1 and a starting configuration with energy E0, the
rules to update these histograms are:

I Propose a new configuration and calculate its
energy E1. The new configuration is generally
derived from the previous configuration by a
change in the value of one of its DoF.

II Accept the new configuration ac-
cording to a probability given by:
min [1, exp(S(E0) −S(E1))].

III Update the two histograms in the correct en-
ergy bin E (E1 if the new configuration is
accepted and E0 if otherwise), accordingly:
H(E) = H(E) + 1 for the control histogram,
and S(E) = S(E) + f for the entropy. Here, f
is a correction that controls the precision of the
algorithm, which will be decreased (see next
step), usually starting at f = 1.

IV If the control histogram H(E) is flat enough
according to some arbitrary criterion, decrease
f (for instance by making f = f/2) and reset
all entries of H(E) to zero.

V If f > � (with � a prescribed tolerance), return
to step I, otherwise stop.

After each reduction of the correction term f, the
entropy histogram is built with a finer grained pre-
cision. However, to speed up the initial estimates of
S(E), f is set to a high initial value. Different ap-
proaches are followed to accelerate the final stages
of refinement by decreasing f with alternative cri-
teria [14].

In Monte Carlo approaches, the detailed balance
condition ensures that the Markov chain has a lim-
iting distribution [15]. Detailed balance can be

stated as follows:

pµP(µ → ν) = pνP(ν → µ) (1)

where pµ is the probability distribution of configu-
ration µ and P(µ → ν) the transition probability
from configuration µ to configuration ν which can
be written as:

P(µ → ν) = SP (µ → ν)AP (µ → ν) (2)

with SP (µ → ν) the selection probability, which
is the probability that the algorithm generates a
trial configuration ν starting from configuration µ;
and AP (µ → ν) the acceptance probability, i.e., the
probability that the algorithm will accept the trial
configuration ν. Hence, since the target distribu-
tion in the entropic sampling is the inverse of the
DoS, i.e., pµ ∝ g(Eµ)−1 = exp[−S(Eµ)], Eqs. (1)
and (2) implies

AP (µ → ν)
AP (ν → µ)

= exp[S(Eµ)−S(Eν)]
SP (ν → µ)
SP (µ → ν)

(3)

In WL, detailed balance is not fulfilled in general.
During the construction of entropy, the acceptance
probability given in step II above (i.e., AP (µ →
ν) = min [1, exp(S(Eµ) −S(Eν))]) evolves and it is
only when the entropy gets sufficiently refined that
detailed balance is met.

Notice that the form chosen for AP requires
that the selection probability be symmetric (i.e.,
SP (µ → ν) = SP (ν → µ)). Although this condi-
tion is simple to comply with in models like Ising,
for arch-based descriptions of static granular packs
this is non-trivial. Therefore, one must be care-
ful to represent a system in such a way that trial
moves between different configurations have the
same backward and forward selection probability.

After reviewing the main characteristics of the
arch-based ensemble in the next section, we will
show that some natural representations of the mi-
crostates lead to non-symmetric selection probabil-
ity schemes. We present, however, a way of repre-
senting the configurations that does allow for the
direct use of WL. In all the WL simulations, we
have used 300 bins for the histograms. The tol-
erance for the f correction was set to � = 2−15.
The histogram H(E) is considered flat (see step IV)
whenever (Hmax −Hmin)/Hmax < 0.2, with Hmax
and Hmin the maximum and minimum height of
the histogram.

070001-3



Papers in Physics, vol. 7, art. 070001 (2015) / D. Slobinsky et al.

III. Arch-based microstates

In Ref. [6], we have introduced a way of describing
the microstate of a static granular system under
gravity by considering the arches that the grains
form instead of the more traditional approach of
using the particle positions. Arches are defined as
sets of mutually stable grains. All other particles
being fixed, the removal of any of the grains in the
arch would induce the destabilization of the rest
of the set. Any assembly of grains, static under
gravity, can be split into a number of arches which
are mutually exclusive [16, 17].

The major difficulty in sampling static granular
configurations is the fact that these are sparse (with
zero measure) in the overwhelming number of pos-
sible particle positions. Moreover, there are not
recipes to generate a static configuration from an-
other static configuration by simply moving a grain
from its position.

Since each arch is stable on its own right, the
arch-based description warrants that any configu-
ration proposed fulfils a basic stability constraint;
i.e., that each set of grains identified as an arch has
internal contacts that keep it stable. The problem
is now moved to generating all possible combina-
tion of arches, including how many of them are of
a given size, shape, orientation and position and
also generating all arrangements of these that can
be stable resting on each other. Of course, all of
these DoF can be represented with different levels
of approximation.

In Ref. [6], we have described the five general
steps necessary to carry out an Edwards entropy
calculation (i.e., the number of states associated to
each given volume of the packing) within an arch-
based scheme. These are:

1. Define the microstate of the system in terms
of arches.

2. Define the external constraints imposed to
arches.

3. Define a volume function that yields the to-
tal volume of the microstate in terms of the
arches.

4. Define an algorithm to generate all microstates
defined in step 1 that comply with the external
constraints of step 2. Or sample microstates
with equal probabilities.

5. Calculate the volume of each microstate gen-
erated in step 4 using the function in step 3
and build the DoS.

Step 4 may constitute a significant limitation to
the real possibility of calculating the density of
states for systems of reasonable size. Generating
all configurations is certainly impractical for most
models (especially if they have continuous DoF).
This paper demonstrates how to sample configura-
tions by using WL (rather than generating all of
them) to accomplish step 4.

We will focus on a model we have already solved
exactly by counting every single possible configu-
rations so as to have a reference system to com-
pare with the WL algorithm results. This is the
“non-interacting rigid arch model” (NIRA). In this
model, only the number of arches ni of each size i
(in number of grains) that form part of the pack-
ing is used to describe the microstate. All arches
consisting of the same number of grains are consid-
ered to occupy the same volume (hence one single
possible shape is assumed for each arch size; this
is implied in the word “rigid” used to name the
model). The total volume of the system is assumed
to be the sum of the volume of the individual arches
and the arch–arch interlocking DoF is not taken
into account (hence the word “non-interacting”).
To represent a two-dimensional pack of equal-sized
disks, we have taken the volume vi of an arch of
i grains to be the area under the regular polygon
that inscribes all disks in a “regular” arch. The
special cases of “arches” of one particle and two-
particle arches are considered separately [6]. The
total volume V [{ni}] of the system is

V [{ni}] =
N∑
i=1

vini where v1 =

√
3d2

2
, v2 = 2.1v1,

(4)

vi =
Nd2

4
tan

π

2N

[
1 +

(
tan

π

2N

)−1]2
for i > 2.

Where N is the total number of disks of diameter
d in the system.

It is important to mention that, typically, the
maximum size that an arch can take is physically
bounded (e.g., due to the size of the container that
holds the granular sample). Hence, we will put a
cutoff C to the largest arch allowed in the system.

070001-4



Papers in Physics, vol. 7, art. 070001 (2015) / D. Slobinsky et al.

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.9 1 1.1 1.2 1.3 1.4 1.5

E
n

tr
o

p
y
/P

a
rt

ic
le

Volume/Particle

cut 3
cut 4
cut 5
cut 6

(a)

-15

-10

-5

0

5

10

15

0.9 1 1.1 1.2 1.3 1.4 1.5

C
o

m
p

a
c
ti
v
it
y

Volume/Particle

cut 3
cut 4
cut 5
cut 6

(b)

Figure 1: (a) Entropy as a function of volume for
the NIRA model in 2D calculated by counting all
possible configurations for 500 disks [6] using differ-
ent arch size cutoffs. (b) The corresponding com-
pactivity calculated by numerical differentiation.

The cutoff C is an external constraint (that is im-
posed in step 2, above). Importantly, this cutoff
imposes a limit to the correlations in the system
which leads to an extensive entropy (see Ref. [6])

In counting microstates, one has to bear in mind
that arches of the same size are indistinguishable in
the NIRA model, whereas arches of a different size
can be distinguished. Hence, if there are ni arches
of i grains in a configuration (with i = 1, ...C, be-
ing C the size cutoff), then the number of permuta-
tion of arches that yield distinguishable microstates
with these arches is

NA!/(n1!...nC!), (5)

where NA =
∑
i ni is the total number of arches of

the configuration (including those “arches” of size
1).

Despite all the simplifications, the model applied
to a 2D system of equal-sized disks yields quali-
tative agreement with DEM simulations of tapped
disks [6]. The NIRA model is in many respects sim-
ilar to an ideal gas of excitations or quasi-particles

(the arches) with a single DoF (their size).
Figure 1(a) shows the entropy S calculated by

counting all possible microstates for 500 disks us-
ing different C for the largest arch allowed [6]. As it
is expected, if C increases, looser configurations are
possible and hence states with higher volumes be-
come more significant. However, the entropy for
low volumes rapidly converges to a well defined
curve. Part (b) of Fig. 1 shows the compactivity
χ defined as χ−1 = ∂S/∂V . This is the analogue
to the temperature in thermal systems [1]. The
entropy presents a maximum as observed by others
[2,3]. States for volumes beyond this maximum cor-
respond to negative compactivities [see Fig. 1(b)].
This is caused by the inversion population of these
volume bounded systems. Some authors suggest
these negative χ macrostates may be inaccessible,
however, this does not need to be the case; some
preparation protocols may indeed lead to very low
packing fractions [3, 18]. An interesting prediction
of the NIRA model is that systems constrained by
different C achieve the same χ at different specific
volume V/N (i.e., at different packing fractions).
As a consequence, two samples of grains “equili-
brated” to the same compactivity will show distinct
packing fractions if the maximum arch size possi-
ble in each sample is different. Different values of C
in practice may be achieved by using narrow con-
tainers or by changing the static friction coefficient
of the grains. There have been some progress in
the study of the equilibration of vibrated granular
samples in “contact” [19]. However, there are still
no attempts to couple static granular packs under
gravity. Further developments in this direction may
help validating this prediction of the NIRA model.

The configurations of the NIRA model are com-
patible with different representations. In the fol-
lowing subsections we will discuss some of these
representations and their suitability for the imple-
mentation of the WL algorithm.

i. The arch size distribution representation

In our previous paper [6], we have used a vector
{ni} that represents the number of arches consist-
ing of i grains in the configuration, i.e., {ni} =
(n1,n2, ...,nC), with C the largest arch allowed in
the system and n1 the number of grains not forming
arches. As an example, a possible configuration in
a system of N = 10 grains and a cutoff arch length

070001-5



Papers in Physics, vol. 7, art. 070001 (2015) / D. Slobinsky et al.

of C = 6 represented in this way could be:

(3, 1, 0, 0, 1, 0). (6)

In this example, there are three grains not forming
arches, two grains forming an arch of size two, and
five grains forming another arch of size 5.

In Ref. [6], we have swept all possible configura-
tions and multiplied each by its analytical degen-
eration due to the different permutations of arches
with repetitions given by Eq. (5).

It is difficult to propose an algorithm to move
between configurations represented in this way and
yet comply with the symmetry of the selection
probability required by the WL algorithm. For ex-
ample, consider a move that consists in removing
a grain from one arch of size k and adding it to
another arch of size k′. Such move would require
subtracting 1 from the coordinate k of {ni} and
adding 1 to the coordinate k − 1, since this arch
will now be smaller by one grain. Additionally, the
coordinate k′ of {ni} needs to be reduced in 1 and
the coordinate k′ + 1 increased in 1, since this arch
will now be part of the set of arches larger by one
grain. There are different ways of selecting a grain
to be moved and to select its arch of destination.

Let us consider, for instance, that all grains and
destination arches are chosen with same probabil-
ity. Now consider a move in the Markov chain that
takes configuration (6) [i.e., µ = (3, 1, 0, 0, 1, 0)]
into configuration ν = (2, 1, 0, 0, 0, 1). This cor-
responds to taking one grain that was not forming
an arch and inserting it in the five-particle arch to
make it a six-particle arch. The probability of se-
lecting a particle from an “arch of size one” in this
case is 3/10. The probability of choosing the arch
of size five as the destination is 1/5 (there are four
other arches in configuration µ plus the possibil-
ity of leaving the grain on its own without forming
an arch with others). Hence the selection proba-
bility is SP (µ → ν) = 3/50. A similar analysis
shows that to return to the original configuration
SP (ν → µ) = 6/10×1/4 = 3/20 (there are 6 grains
out of ten that can be taken from the six-grain
arch and there are three possible other destination
arches plus the case with the grain not forming an
arch in the new configuration). Clearly, this selec-
tion probabilities are non-symmetric as they should
be to apply the algorithm of section II.

A possible workaround to the previous represen-
tation is to multiply each coordinate of the vec-

tor {ni} by the corresponding arch size. Therefore,
each coordinate now indicates how many grains are
involved in all arches of the given size. In this new
representation, the configuration of Eq. (6) (10 par-
ticles with a cutoff C = 6) is written as

{n′i} = (3, 2, 0, 0, 5, 0) (7)

In this case, a trial move may consist in ran-
domly transferring a grain from one coordinate k
to another coordinate k′. This is done simply by
subtracting 1 to {n′k} and adding 1 to {n

′
k′}. In this

representation, the selection probability of moving
a particle from one arch of size k to another of size
k′ is SP (µ → ν) = 1/N×1/(C−1) irrespective of k
and k′. Hence, the selection probability is symmet-
ric and suitable to implement the entropic sampling
using WL.

Unfortunately, the representation of Eq. (7) does
not tell apart two microstates that differ only in the
permutation of two arches of different sizes. There-
fore, the corresponding degeneracy given by Eq. (5)
cannot be accounted for 1. More importantly, the
change of representations from Eq. (6) to Eq. (7)
is not an exact mapping because these newly pro-
posed moves allow for “fractions of arches” to exist
since we do not request that the coordinates of {n′i}
be a multiple of i after each move. For instance,
the configuration (0, 0, 0, 0, 0, 10) is allowed in rep-
resentation (7) but does not exist in representation
(6). In the new representation such configuration
implies that there is one arch of size 6 plus a frac-
tion of an arch of size six. We can still treat these
“fracrtions of arches” by assigning to them a frac-
tion of the arch volume. However, there is a large
number of new unrealistic configurations (there are
many ways of chosing a numbers that are not con-
mensurate with the arch size) in this representation
that will bias the result for the entropy.

Figure 2 shows the entropy for the NIRA model
using representation (7) for different maximum
arch sizes C compared with the exact result ob-
tained by counting all configurations and all per-
mutations [6]. As we can see, not including all dis-
tinguishable permutations and including the new
“fractional arches” gives a wrong entropy function.

1One is tempted to add to the degeneracy factor (5) to
correct the entropy in step III of the algorithm. However,
the algorithm compensates this factor in order to obtain a
flat histogram. Therefore this is not a viable solution.

070001-6



Papers in Physics, vol. 7, art. 070001 (2015) / D. Slobinsky et al.

0

0.5

1

1.5

2

0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8

E
n

tr
o

p
y
/P

a
rt

ic
le

Volume/Particle

cut 5
cut 6
cut 7
cut 8

cut 4

Figure 2: Entropy as a function of volume for the
NIRA model calculated using the WL algorithm
(symbols) for representation (7) for 4 < C < 8.
The full lines represent the exact results for 200
grains.

ii. The arch listing representation

As previously discussed, other ways of represent-
ing the system should be carefully chosen in order
to ensure that there exist moves that sample the
different configurations uniformly. An alternative
natural representation consists of using a vector,
{mi}, with N coordinates, where each coordinate
mi can take any value from 0 to C, the cutoff for the
arch size, provided that

∑
i mi = N. The content

of each coordinate indicates that the configuration
has an arch of that size. The configuration of Eq.
(6) in these representations can be expressed, for
example, as

(5, 0, 0, 1, 0, 1, 2, 0, 1, 0). (8)

There are, of course, multiple permutations in Eq.
(8) that lead to the same distribution of arch sizes
compatible with Eq. (6).

Trial moves in the system represented in this way
can be done by subtracting one from a non-zero mi
and adding one to any other coordinate that has a
value smaller than C. This is equivalent to reducing
the size of one arch in one particle and either creat-
ing a new arch of size one (if the new coordinate had
a zero value) or increasing the size of another arch.
Unfortunately, this algorithm has a non-symmetric
selection probability. SP = 1/NA × 1/(N − NC),
where NA is the number of arches (i.e., the number
of non-zero coordinates in the configuration) and

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9

E
n
tr

o
p
y
/P

a
rt

ic
le

Volume/Particle

cut4
cut5
cut6
cut7
cut8

Figure 3: Entropy as a function of volume for the
NIRA model with cutoff 4 ≤ C ≤ 8 (symbols as
in Fig. 2) using the binary arch representation to
carry out the entropic sampling through the WL
algorithm for 200 grains. The solid lines correspond
to the exact counting of microstates from Ref. [6].

NC is the number of arches with the maximum al-
lowed size C. Therefore, SP will depend on the
total number of arches and the number of arches of
size C.

Besides the non-symmetric SP , the system rep-
resented in this way and sampled with these trial
moves clearly overestimates the number of states
of a given volume. This is because, apart from the
NA!/(n1!...nC!) permutations of the distinguish-
able arches [see Eq. (5)], there are N!/(N −NA)!
additional permutations due to the zeros in a given
vector in Eq. (8) (there are N −NA zeros).

One can, in principle, resolve these issues. The
selection probability may be turned into a sym-
metric one by adapting the acceptance probability
in Eq. (3). The degeneracy due to the presence
of zeros in the representation (8) can be handled
by switching to a representation without including
these, and allowing for a vector of variable length.
However, although less intuitive, there is a much
simpler, suitable representation that we discuss in
the next section.

iii. The binary arch representation

Finally, we present a representation that comply
with the symmetric selection probability and si-
multaneously account for the permutation of dis-
tinguishable arches.

070001-7



Papers in Physics, vol. 7, art. 070001 (2015) / D. Slobinsky et al.

In this case, the system is chosen to be repre-
sented by a vector of N coordinates with binary
values (zeros and ones). These coordinates are not
associated to specific particles 1, 2, 3, etc. Rather,
as we move along the vector from left to right, we
can think the ones as representing a first grain in
an arch (whatever its identity) and the following
zeroes as the remaining particles. For instance, the
configuration in equation 6 in this new representa-
tion can be given by

(1, 0, 0, 0, 0︸ ︷︷ ︸
5

, 1︸︷︷︸
1

, 1︸︷︷︸
1

, 1, 0︸︷︷︸
2

, 1︸︷︷︸
1

). (9)

In this representation, all permutations of arches
of different sizes are accounted for naturally. The
sections of the vector representing an arch [under-
braced in Eq. (9)] can be permuted to yield all
distinguishable configurations [see Eq. (5)], which
correspond to different vectors in this binary rep-
resentation. Indistinguishable configurations corre-
sponding to permutations of arches of same size are
also indistinguishable for the binary vector.

The number NA of arches in a configuration is
simply the sum of all the elements of the vector.
Note that the underbraced numbers coincide in
value and order with the non-zero figures of the
vector described by Eq. (8).

The trial moves consist in picking a coordinate
and changing the state of that coordinate (if 1
change to 0 and vice versa). This results in a
symmetric selection probability of SP (µ → ν) =
SP (ν → µ) = 1/N. In each move, the constraint
imposed by the cutoff C must be checked and the
trial configuration must be rejected whenever the
constraint is not complied with.

In Fig. 3, the result for this representation and
sampling strategy is plotted along with the exact
result showing a remarkable agreement.

IV. Conclusions

We have been able to compute the entropy of a
system of non-interacting rigid arches using a WL
algorithm in the volume ensemble in different rep-
resentations.

We have exposed the difficulties in dealing with
different representations of the configurations of
arches and the mechanisms used to propose trial

moves for the WL algorithm. These difficulties ap-
pear during the choice of a simple sampling scheme
that ensure a symmetric selection probability of
configurations. Additionally, the degeneracy due
to distinguishable permutations of arches pose a
further complication in the use of WL. The most
suitable representation that we found for a non-
interacting system of rigid arches resulted in a bi-
nary vector.

We believe that entropic sampling of arches
through the WL algorithm has a great potential
for testing the granular statistical mechanics hy-
pothesis (such as equiprobability and ergodicity).
Having a sampling algorithm like WL adapted for
these types of models is crucial to continue the
road map towards the refinement of an arch-based
framework for static granular packs. In particular,
the non-interacting condition is clearly a crude ap-
proximation and should be lifted, along with the
introduction of a more accurate volume function.

[1] S Edwards, R Oakeshott, Theory of powders,
Physica A 157, 1080 (1989).

[2] S McNamara, P Richard, S K De Richter, G Le
Caër, R Delannay, Measurement of granular
entropy, Phys. Rev. E 80, 031301 (2009).

[3] M P Ciamarra, A Coniglio, Random very loose
packings, Phys. Rev. Lett. 101, 128001 (2008).

[4] D Asenjo, F Paillusson, D Frenkel, Numeri-
cal calculation of granular entropy, Phys. Rev.
Lett. 112, 098002 (2014).

[5] C S O’Hern, L E Silbert, A J Liu, S R Nagel,
Jamming at zero temperature and zero applied
stress: The epitome of disorder, Phys. Rev. E
68, 011306 (2003).

[6] D Slobinsky, L A Pugnaloni, Arch-based con-
figurations in the volume ensemble of static
granular systems, J. Stat. Mech. P02005
(2015).

[7] J Lee, New Monte Carlo algorithm: Entropic
sampling, Phys. Rev. Lett. 71, 211 (1993).

[8] F Wang, D P Landau, Efficient, multiple-
range random walk algorithm to calculate the
density of states, Phys. Rev. Lett. 86, 2050
(2001).

070001-8



Papers in Physics, vol. 7, art. 070001 (2015) / D. Slobinsky et al.

[9] F Wang, D P Landau, Determining the den-
sity of states for classical statistical models: A
random walk algorithm to produce a flat his-
togram, Phys. Rev. E 64, 056101 (2001).

[10] C Zhou, R N Bhatt, Understanding and im-
proving the Wang-Landau algorithm, Phys.
Rev. E 72, 025701 (2005).

[11] S Trebst, D A Huse, M Troyer, Optimizing the
ensemble for equilibration in broad-histogram
Monte Carlo simulations, Phys. Rev. E 70,
046701 (2004).

[12] B A Berg, T Neuhaus, Multicanonical en-
semble: A new approach to simulate first-
order phase transitions, Phys. Rev. Lett. 68,
9 (1992).

[13] D P Landau, K Binder, A Guide to Monte
Carlo Simulations in Statistical Physics, 2nd
ed., Cambridge University Press (2005).

[14] R E Belardinelli, V D Pereyra, Fast algorithm
to calculate density of states, Phys. Rev. E 75,
046701 (2007).

[15] M E Newman, G T Barkema, Monte Carlo
methods in statistical physics, Vol. 13, pp. 36-
42, Clarendon Press Oxford (1999).

[16] L A Pugnaloni, G Barker, A Mehta, Multi-
particle structures in non-sequentially reor-
ganized hard sphere deposits, Adv. Complex
Syst. 4, 289 (2001).

[17] R Arévalo, D Maza, L A Pugnaloni, Identifi-
cation of arches in two-dimensional granular
packings, Phys. Rev. E 74, 021303 (2006).

[18] L A Pugnaloni, M Mizrahi, C M Carlevaro,
F Vericat, Nonmonotonic reversible branch in
four model granular beds subjected to vertical
vibration, Phys. Rev. E 78, 051305 (2008).

[19] J G Puckett, K E Daniels, Equilibrating tem-
peraturelike variables in jammed granular sub-
systems, Phys. Rev. Lett. 110, 058001 (2013).

070001-9