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.