I'm doing some numerical Monte Carlo analysis on the 2 dimensional Ising model at the critical point. I was using the Metropolis 'single flip' evolution at first with success, though it suffers from critical slow down and makes studying large lattices unlikely possible. I'm now looking at cluster flip algorithms, specifically the Wolff algorithm.

I managed to implement it, and it looks to be working as it should (flipping a unique spin at $T = +\infty$, the whole lattice at $T=0$, matches the right energy density in the thermodynamic limit...) but I don't get the right behaviour for the two point $<\sigma_i\sigma_j>$ correlation function.

According to CFT it should behave like:

$$<\sigma_i\sigma_j> \;\propto \;\frac{1}{|i-j|^{\frac{1}{4}}}$$

I'm more and more convinced that it has to do with boundary conditions, I use non periodic free boundaries. The literature on the subject doesn't say much on this point.

Am I missing a subtlety (or an evidence) in this procedure, or in the use of this algorithm?

  • 6
    $\begingroup$ No, this is appropriate for a physics site. Only a physicist can answer this. $\endgroup$
    – suresh
    Commented Mar 5, 2014 at 14:54
  • 1
    $\begingroup$ For any kind of subleties in Ising model simulation I can recommend "Monte Carlo Simulation in Statistical Physics" from Binder and Heermann. Other than that have you tried with periodic conditions instead? $\endgroup$
    – Alexander
    Commented Mar 5, 2014 at 19:52
  • 4
    $\begingroup$ @DavidZ Yes, it is about physics. $\endgroup$
    – suresh
    Commented Mar 7, 2014 at 2:24
  • 1
    $\begingroup$ @DavidZ it's hard to imagine how a question about correlations in the Ising model could be said to be not about physics. $\endgroup$
    – N. Virgo
    Commented Mar 7, 2014 at 4:54
  • 1
    $\begingroup$ @Learning how big is your grid? Since the spatial range of the correlations also diverges at the critical point, the boundary conditions will have a strong effect on the whole system. I guess this can be mitigated somewhat in practice by making it very large, which is why I'm curious about your system's size. $\endgroup$
    – N. Virgo
    Commented Mar 7, 2014 at 5:00

2 Answers 2


I have tried to reproduce your problem using my own code, but I couldn't: I got the correct value for the exponent. I can't tell you what went wrong in your calculation, but I can tell you exactly what I did!

You can inspect my code on GitHub. It was the first thing I ever wrote in C++ and today I would do many things differently, so please don't judge the code too harshly. There is a readme file in the repository that explains pretty much everything about it.

I uploaded my runscript, the relevant part of the corresponding results and the small Pyxplot script that does the fit into this GitHub Gist. My value for the exponent is: $$\beta = 0.243 \pm 0.001$$ You can see the fit plotted together with the Monte Carlo data in the sscorr_fit.pdf file.

Here is a list of things that I think might be relevant.

  • I'm sure you know that, but just for the record: In the equation in the question one is not supposed to take the absolute value of the difference of the indices $i$ and $j$ of the two sites. The indices you assign to the sites are completely arbitrary. What matters is the distance between the two sites, so I would prefer to write it like this: $$\langle s_i s_j \rangle \propto \frac{1}{|\vec r_j - \vec r_i|^\frac{1}{4}}$$

  • I only measured the correlation along the direction of the lattice vectors and not along any diagonal. That should not matter though since there is a proof that the spin-spin correlation is rotationally symmetric at the critical point.

  • I only fitted your formular in the range $1 \leq |\vec r_j - \vec r_i| \leq 50$ for a 512x512 model. I think if you go beyond that you get too much influence across the periodic boundary conditions. You can definitely not fit the entire $1 \leq |\vec r_j - \vec r_i| \leq 256$ range, as the correlation must have a minimum at $|\vec r_j - \vec r_i| = 256$.

  • The correlation function changes quite quickly with temperature around the critical point, so make sure you have the right temperature and everthing is equilibrated. In my experience it is best to start the Wolff algorithm with a model where all spins point in the same direction. If you take random spins and then cool down from $T=\infty$ to $T_c$ the Wolff algorithm is extremely inefficient in the beginning as it can not build large clusters in this noise of random spins.

Hope that helps!

  • $\begingroup$ Thank you very much for the input, I will look deeper at your code. Can I ask you how many systems did you simulate (ie. are you averaging on) ? I managed to derive a good coefficient in my case too (very close to 0.249 ) but my data was noisier, probably the sum used to fit (linearly in log-log scale) cancelled most of the noise. Also how many Wolff evolution steps do you perform to thermalize starting from ordered T = 0 configurations? Many thanks again ;) $\endgroup$ Commented Jul 7, 2014 at 16:41
  • 1
    $\begingroup$ @Learningisamess: Let me reiterate what I said as a comment to the other answer. You can completely avoid the issue of determining the number of steps required to thermalize: just use a perfect simulation algorithm; when the latter stops, you are guaranteed that the system has thermalized! One such (cluster) algorithm is described in this paper: M. Huber, a bounding chain for Swendsen-Wang, Random Structures & Algorithms, 22(1) :43–59, 2003. It is trickier to implement than the standard algorithm, but it is very efficient and you get all the benefits of perfect simulation. $\endgroup$ Commented Jul 7, 2014 at 17:44
  • $\begingroup$ Everything I posted was obtained from a single simulation, the runscript of which is part of the git repo I linked. So I did not do any averaging beyond averaging within the Markov chain itself. I did run a few simulations at slightly different temperatures, but I did not include their results in what I posted above. $\endgroup$ Commented Jul 7, 2014 at 23:14
  • $\begingroup$ For equilibration I used 100 Monte Carlo steps (MCS). One MCS is supposed to flip roughly a number of spins equal to the number of lattice sites. How many Wolff clusters are needed for that is determined at runtime by measuring the average size of the last clusters. The whole point of all this business is to make the MCS independent from the number of lattice sites, so that the autocorrelation time measured in MCS is just a constant number that does not depend on the system size. Check the settings in the runscript against the SSMC readme file to see what I did exactly. $\endgroup$ Commented Jul 7, 2014 at 23:24
  • $\begingroup$ Interesting, I have a different way of reaching equilibrium I would first monitor the value of the mean energy per spin (for each configuration) and read from there how many steps I needs for each systems size (I usually stick to 3,4 sizes also). I did check if starting from a totally ordered configuration is quicker than a totally disordered, and well amusingly it depends on the boundary conditions: with fixed boundaries, disordered is very very advantageous (the reason it that the Wolff algorithm needs to be modified with fixed boundaries with an accepting factor 'à la Metropolis' $\endgroup$ Commented Jul 8, 2014 at 14:24
  1. It is only exactly at the critical temperature that this CFT result works. You haven't mentioned if you have used the critical temperature when you did the monte-carlo.

  2. At/near critical point, autocorrelation time becomes huge. (If I am not mistaken, autocorrelation time must blow up exactly at critical temperature, however it is cut-off due to finiteness of system). So, its best to record measurements once every 3-5 autocorrelation times. Notice that, for a given random seed, autocorrelation effects can give a systematic error, while statistical error due to autocorrelation effect can be estimated by using different random seeds. Look into this book: http://www.amazon.com/Monte-Carlo-Methods-Statistical-Physics/dp/0198517971

  3. My guess is that, as long as u are 10-20 spins away from the boundary when u do the measurement of spin-spin corr., you wont see boundary effects. This is only a guess. Why dont you put in periodic boundary conditions and see? The neighbour of i,j is (i+1)%L, (j+1)%L and 3 others like this.

  • $\begingroup$ Is it easy to tune the temperature to be at $T_c$? $\endgroup$
    – suresh
    Commented Mar 7, 2014 at 2:28
  • 1
    $\begingroup$ For ising, we know the exact value (2.27...). In general, for 2nd order transitions, the critical temperature is best estimated by plotting the magnetic 4th order cumulant and noticing the crossing point. $\endgroup$ Commented Mar 7, 2014 at 4:36
  • $\begingroup$ First, thank you for the input Srivatsan. To answer your questions: 1. Yes, the temperature is set to the critical value: I start with a totally random lattice (so $T = \infty$) and then I 'cool it down' using Wolff algorithm evolution with probability to add the neighboring state of $1 - \text{e}^{-2*K_c}$ where Kc is value close enough to $\frac{1}{T_c}$ for the correlation length at criticality being larger than my system size (I'm using lattices between 512x512 to 4096x4096 so 6-7 digits precision in the coupling). $\endgroup$ Commented Mar 7, 2014 at 11:00
  • $\begingroup$ 2. Thank for the reference, actually I had already studied this book. Indeed they show how the autocorrelation time increases drastically at the critical point specially with single spin-flip (ie. Metropolis) "dynamics". That's why I'd like to make it work with the Wolff algorithm, I succeeded in calculating the right weight with the Metropolis algorithm (and free boundary conditions) but it took a lot of computation time and a more effective method will turn really useful in the next steps of my project. $\endgroup$ Commented Mar 7, 2014 at 11:08
  • $\begingroup$ 3. Since my first post, I put the system on periodic boundary conditions and I found a new value for the conformal weight of the spin operator: 0.05 . It still isn't the correct one (0.125) but it's on the "other side" as fixed boundary conditions had given me 0.4 . At least it tells me that boundary condition influence this derivation so we may on the right path, but I'm also not so sure of this calculation as it was only done over 5000 samples. I'll check that and maybe come back to you later. Again, I appreciate your input ;) $\endgroup$ Commented Mar 7, 2014 at 11:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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