Ising model with metropolis algorithm around critical temperature

I'm trying to simulate Ising model using metropolis algorithm. Boundary conditions are periodic. I know how the algorithm works and I have written the code myself. Everything works as it should except around critical temperature.

I started my calculation at temperature T=8 (Boltzmann constant is 1). And then decreased the temperature in steps by 0.05. At T=2.35 I shortened the step to 0.005 till T=2.2 and after 2.2 the step was again 0.05. At each temperature I waited 15 million flips before I started sampling magnetization. Then I flipped the spins 25 million times and every 50th flip I sampled magnetization. This is necessary in order to have statistically uncorrelated samples. Btw, I am using Mersenne twister from gsl to generate random numbers.

Still around critical temperature the algorithm does not work good enough. Magnetization is jumping up and down. I don't really understand why that is, because I think I took enough samples to average out extremes. If someone can help me, I'll be really happy. Below are pictures of magnetization, susceptibility and specific heat.   • You're simply not averaging over enough runs. Apr 17 '15 at 19:19
• The algorithm is taking quiet a long time to finish, so I can't really increase the number of samples by a lot. Now I'm sampling magnetization after 25th flip of spin, and I also increased 25 million to 35 million flips when sampling and it doesn't seem to help. Another question I have is, whether do I have to average absolute value of magnetization or no. Now I just sum magnetization without changing the sign to +, around critical temperature it seem to be jumping from - to +. The graph above shows absolute of average magnetization.
– Jur
Apr 17 '15 at 20:26
• From figure 1 and especially figure 2 I can see you are incorrectly sampling magnetization. Are you using some sort of cluster algorithm i.e. Wolf algorithm ? You should sample the absolute value of magnetization, since it can change sign. P.S. From figure 3 it can be seen you are correctly sampling energy.
– user78758
Apr 25 '15 at 21:06
• Those graphs are actually from one single run of algorithm. I sampled <E>, <E^2>, <M> and <M^2> and then calculated different things. After I started sampling absolute value of magnetization things got better, however, the lines are still not smooth - still some ups and down occur at critical temperature. I tried wolff algorithm and (for me) it works a lot better than metropolis.
– Jur
Apr 27 '15 at 7:34
• This is called critical slow down. One way to get better statistics and speed up the simulation is the Swendsen-Wang algorithm for flipping clusters of spins rather than individual spins. May 11 '15 at 14:01