I am working on Metropolis-Montecarlo algorithm for 2D Ising model in python partly based on this document. I ran the simulation for 100 times on a 25 x 25 lattice with external magnetic field B = 0. I couldn't get the phase transition when I plot the Temperature vs Magnetization, rather it's more fluctuating.


Initial Lattice

Final Lattice

Temperature vs Magnetization


This code will generate the same output as above.

# required packages
import numpy as np
from numpy.random import rand, random, randint
import matplotlib.pyplot as plt

# definitions
def spin_field(rows, cols):
    ''' generates a configuration with spins -1 and +1'''
    return np.random.choice([-1, 1], size = (rows, cols))

def neighbours(x, y, lattice, dim):
    ''' finds the neigbours of a particular lattice point with periodic boundery conditons '''
    left   = (x - 1, y)
    right  = ((x + 1) % dim, y)
    top    = (x, y - 1)
    bottom = (x, (y + 1) % dim)

    return [lattice[left[0], left[1]],
            lattice[right[0], right[1]],
            lattice[top[0], top[1]],
            lattice[bottom[0], bottom[1]]

def energy_calc(lattice, dim, B):
    ''' calulates the energy of the whole configuration '''
    # coupling constant
    j = -0.44
    en = 0
    for x in range(0, dim):
        for y in range(0, dim):
             en += -j * lattice[x, y] * np.sum(neighbours(x, y, lattice, dim)) - B * lattice[x, y]
    return en

# number of monte carlo simulations
mcs = 100

# square lattice dimensions
dim = 25

# external magnetic field
B = 0

mag = []

# initialisation of lattice with random spins
lattice_0 = spin_field(dim, dim)

for temp in np.linspace(0.1,4,50):
    for i in range(0, mcs):
        # picks a random lattice point
        x = randint(dim)
        y = randint(dim)
        lattice_1 = lattice_0

        # flip the spin
        lattice_1[x, y] *= -1
        H0 = energy_calc(lattice_0, dim, B)
        H1 = energy_calc(lattice_1, dim, B)

        # difference in energy of the configuations
        dE = H1 - H0
        if (dE <= 0):
            lattice_0 = lattice_1
        elif (np.exp(-1 * dE / temp) >= random()):
            lattice_0 = lattice_1

    mag.append(abs(sum(sum(lattice_0))) / (dim * dim))

T = np.linspace(0.1,4,50)
plt.plot(T, mag)

2 Answers 2


I think that you have a bug in the code where you put $j=-0.44$ and also calculates the energy as $E=-J\sum S_i S_j$ which means that you actually calculates the Anti-ferromagnetic Ising model. If I am correct then this should go a long way to explaing why the magnetization per site keeps close to zero.

Another point is that I believe that $100$ MC steps for a $25\times 25$ lattice is a very small number, which does not allow the algorithm to scan properly the available configurations. You have 625 lattice sites, which means that at least 525 of them will remain at their initial values. You should either increase the number of steps (to the order of $10^3-10^4$) or decrease the size of your system, which is quite large I think, or both of course. At any rate I would recommend that the number of MC steps would be at an order of magnitude larger than the number of lattice sites.


You main issue is that you are doing far too few Monte Carlo steps. Typically, one Monte Carlo sweep is defined as $N=L^2$ attempted spin flips, so instead of for i in range(0, mcs): you should do for i in range(0, mcsdimdim):, that way the number of attempted updates scales with the system size.mcs` should be at least $10^3$, probably $10^4$ to be safe (this will still be very fast).

You have a few other issues that might come up once you fix the main problem:

  • As yu-v pointed out, you have set j = -0.44 and you have a minus sign in front of -j *...., so you are doing the antiferromagnetic Ising model.
  • You do no averaging for your magnetization measurement. This isn't necessarily problem for a demonstration code, but it will mean your results will have large statistical error.
  • You calculate the energy change of the entire system for your acceptance probability. This is not wrong, but it is extremely inefficient. Since you have only nearest-neighbor interactions, you only need to account for the neighbors of the spin you are attempting to flip, which will save a ton of CPU time.
  • Finally, you can get substantially better results if you were to go from high to low temperature and not reinitialize the spin state for each new temperature. Then you will effectively be doing simulated annealing, which should give you much better results, especially at low temperature.

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.