Summary: In Wang Landau method, you have to calculate the probability $\exp (- \Delta S)$ but $\Delta S$ is generally large when the size of the system is not small. How can I make $\Delta S$ small?

I'm writing a monte carlo simulation with Wang and Landau algorithm. When the size of the system is small, the density of states is correctly gotten. However, when the size is large (*1), much more iterations are needed to flatten the energy histogram and thus the entropy of each energy bin becomes large. For example, when the entropy refinement parameter $f = 1.0$ (note this is the initial value), the energy histogram $H$ is given as below.

H[0] = 618479
H[1] = 619151
H[2] = 620313
H[3] = 620378
H[4] = 620409
H[5] = ...

And entropy $S$ is also given as below. (Since $f = 1.0$ now, $H$ and $S$ have the same values.)

S[0] = 618479
S[1] = 619151
S[2] = 620313
S[3] = 620378
S[4] = 620409
S[5] = ...

These histograms seem flat enough, but the difference of elements is of order $O(1000)$. We have to decide whether or not to accept the suggested new states with probability $\exp (- \Delta S)$, but this is zero when $\Delta S$ is large. So the only states whose $S$ is small are accepted and the algorithm doesn't proceed at all.

How can I solve this problem? I think it is not practical to wait until the histogram becomes completely flat. I wonder if there is a more feasible solution.


(*1): I assume $L \times L$ grid graph. I can compute $L \leq 9$ but not $L \geq 10$. You may think $10 \times 10$ grid graph is "not large", though.


  • Actually I'm trying to estimate the number of rook paths of a grid graph, using Wang Landau method rather than strict enumeration. The system is $(L + 1)^2$ nodes placed on a grid (Fig.1). Nearest neighbors are randomly linked. We propose a new state by choosing one (virtual) edge, and flip it (Fig.2) (i.e. if there is a edge there, cut it, and if there isn't, create a new edge). In the ground state, the edges constitute a path which links first node (upper left) and last node (lower right). This is just my challenge, and not an existing model.

  • The histograms above were gotten in a few minutes on my personal computer. I don't know this is long or short, but at least in $L \leq 9$ case (in which histograms are gotten in dozens of seconds or so), the estimated number of solutions are very similar to the true one. So I think my model isn't bad...

enter image description here

enter image description here

  • $\begingroup$ What model are you simulating? How long a first-stage run have you carried out, to reach the histogram that you give in your question? How does this run length compare with similar examples in the literature? Ideally, please edit your question to include these details. $\endgroup$
    – user197851
    Jan 30, 2019 at 7:48
  • $\begingroup$ @LonelyProf Thank you for your response. I just added the details. $\endgroup$
    – ynn
    Jan 30, 2019 at 8:16
  • $\begingroup$ I'm not familiar with the model. In particular, it would be interesting to know how you define the energy function, and how many distinct energies there are (as a function of $L$). But my current thinking is that you may just be running the first stage too long. What is the flatness parameter for your histogram: that is, the ratio of the minimum value in $H(E)$ to the average value of $H(E)$ across all the nonzero bins? $\endgroup$
    – user197851
    Jan 30, 2019 at 10:37
  • $\begingroup$ @LonelyProf Explaining how I define the energy seems a little bit difficult (it's complex), but there are $2^{2L(L+1)}$ states and about $2 L^2$ distinct energy states. Now I've found the solution by myself, but as you pointed out and since I'm a beginner of WL method, there must be other problems left in my program. Forgive me for asking elementary questions. (1)What is the usual length of the first step? $\endgroup$
    – ynn
    Jan 30, 2019 at 11:44
  • $\begingroup$ @LonelyProf (2)I judge $H$ is flat if $(H_{max} - H_{min}) / (H_{max} + H_{min}) < F$, where I take $F = 0.1$. This definition follows my textbook, but I sometimes see your definition. Is your one more usual? If so, what value should I set as the flatness parameter? I'll be glad if you could answer (a part of) my additional questions. Thank you in advance. $\endgroup$
    – ynn
    Jan 30, 2019 at 11:46

2 Answers 2


The effective entropy function $S(E)$ used in the Wang-Landau algorithm is, like the true thermodynamic entropy, an extensive quantity. In general terms, there's nothing you can do about that. You want the sampling to be inversely proportional to the density of states $\Omega(E)=\exp[S(E)]$, in $k_B=1$ units. Ultimately, the method will become impractical for large enough systems.

However, I don't think that's the problem here. I suspect that you are simply running the first stage of Wang-Landau for too long. A typical flatness criterion is when the minimum histogram value comes within, say, 95% of the average value: $\min[H(E)]>0.95 \,\overline{H(E)}$. For a typical test system (2D Ising model on a $10\times 10$ lattice using single spin updating) I would achieve this flatness after visiting each energy around $2000$ times (there are about $100$ distinct energies). The range of histogram values between max and min is around $350$, but the differences between nearby histogram entries are in the range $1\leq \Delta S\leq10$, and this is the quantity of interest to the MC algorithm. So, I've only done $\sim 200000$ spin flips at that stage.

At this point, it is necessary to reduce the entropy-updating factor, reset the histograms, and continue to stage 2, and then similarly onwards until you've reached the desired convergence criterion for $f$. At each stage, typically, the computational effort to achieve the 95% flatness criterion goes up, as the entropy histogram is refined. For my Ising model, after $40$ or so stages, I've invested of order $10^8$--$10^9$ spin flips in total, but the overwhelming majority of the effort comes in the latter stages, each of which might take $100$ times longer than the first stage.

You seem to be visiting each energy more than $600000$ times in the first stage! So, I think the statistical fluctuations on numbers of that size will be working against you. When the system becomes stuck, it needs to accumulate $\sim 1000$ "hits" before it can escape. But this may mean that it is $\sim 1000$ units adrift from other neighbouring states.

I recommend trying again, but implementing the flatness criterion in the usual way, to reduce the entropy updating factor and proceed to the next stage, so your system progressively refines the entropy weights.

  • $\begingroup$ Thank you for your detailed advise. Such information especially for concrete numerical values is really precious for me. I implemented WL method for $N \times N$ 2D Ising model before the grid graph model, so now I just checked if the result from that is consistent with yours. I set $N = 10$ and change the way the flatness is checked. And there is no contradiction, i.e. the first stage ended in a moment ($\simeq O(10^{-2})$ (sec)) and each energy states were visited $\simeq O(10^3)$ times. $\endgroup$
    – ynn
    Jan 30, 2019 at 13:55
  • $\begingroup$ Actually, with my grid graph model, I'm executing so-called rare event sampling - the state of $E = 0$ is very very rare and be seldom proposed (even when $L = 5$, the ground states are $10^{-9}$ (%) of the all states). So maybe the behaviors are quite different between my model and Ising model. Anyway, I'll rewrite the part of flatness check and see if there's any not-small change. Thank you again for your many kind responses. $\endgroup$
    – ynn
    Jan 30, 2019 at 13:57
  • 1
    $\begingroup$ Yes, I have also used WL to simulate systems with low-probability ground states, and you are right, the behaviour can be different. And it is also true that if the program proceeds to the later stages, and only then "discovers" the low-lying state, this can be even worse, since many more entropy increments are then needed to catch up. $\endgroup$
    – user197851
    Jan 30, 2019 at 15:18

Now I've just found the solution. What's wrong was the interval of the flatness check. As I noted in my question, when the size is large, much more iterations are needed to flatten the energy histogram ($=$ there are larger number of distinct energy states), so I had to dynamically change the interval according to the value of $L$ rather than fix it to some constant value. In my case, if $L < 10$, the interval was $10^4$ monte carlo steps. But this was not appropriate for $L \geq 10$. So I set $10^5$ or even $10^6$ monte carlo steps as the interval, and then the problem solved. (Note this modification doesn't change the calculation time since much more iteration are needed to flatten the histogram.) I'm calculating the case of $L = 11, 12, \cdots$ and I'm getting plausible results.

  • $\begingroup$ I am a bit surprised that extending the interval between flatness checks (for larger systems) is the solution. And also surprised that your first stage does not satisfy the flatness criterion much sooner in the first stage. But, as I said, this type of system is outside my usual experience (for instance, the Ising model, as described in my own answer). So, perhaps my answer is wrong; only you can decide. $\endgroup$
    – user197851
    Jan 30, 2019 at 12:27
  • $\begingroup$ You've been kind enough to accept my answer, but really, if your own answer is the correct one, you should accept it instead! $\endgroup$
    – user197851
    Jan 30, 2019 at 18:04
  • $\begingroup$ @LonelyProf I think you are right, so I've changed the accepted answer. Personally your answer was really really helpful though :) $\endgroup$
    – ynn
    Feb 1, 2019 at 7:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.