Monte Carlo Simulation of 2D Ising Model

Monte Carlo Simulation of 2D Ising Model using Metropolis algorithm Python and Numba.


During my first year PhD at UC, I get the chance to take statistical mechanics class from a condensed matter physicist. It was a new experience when he focused too much on different side that I did not explore before. The final project in this class was about using Monte-Carlo method to simulate the Ising model. Ising model is a model of a system that undergoes a phase transition. There are many reasons why this model is interesting when it is not representing any actual physical systems.

  1. It can be solved exactly, We can calculate the thermodynamical quantities and interpret them. In another can calculate the partition function and put it in an exact form.
  2. Ising model exhibits symmetry breaking in low-temperature phase.
  3. It is a simple model that can be used to understand the physics of phase transition.
  4. The Ising model is simple, yet it can be applied to many physical systems. This gives us a feature about universality, a feature that tells us that in such case the same theory applies to all sorts of different phase transitions.

The Ising Model

We can think of Ising model as a mathematical model where everything consists of small squares inside a bigger one. It is usually called lattices where we then have a big lattice of sites (small boxes) where each site can be in one of two states. Those states as we will see will refer to the spin state of the electron inside this site. Each site will be labeled by index ii. The Hamiltonian for the Ising model can be written as

=JijSiSj\mathcal{H}=-J \sum_{\langle i j\rangle} S_i S_j

Lattice square representation of Ising model,
  • Spin SiS_i has value of +1 or -1.
  • ij i j implying that we are interested in nearest-neighbor interaction only.
  • JJ is the interaction strength where ( J>0 ) assuming that strength is positive.
  • There are many systems that can be represented using this model like “lattice gas” with each of the sites is the possible location of a particle; means that site is empty and means that site is occupied by a particle.

    This system will undergo a 2nd order phase transition at some temperature T=TcT=T_c which is called critical temperature. For cases when T<TcT < T_c, the system magnetizes, and the state is called the ferromagnetic state which is an ordered state. For cases when T>Tc T > T_c ​, the system is disordered or in a paramagnetic state. We can define the order parameter for this system to be the average magnetization mm.

    m=SNm=\frac{\langle S\rangle}{N}

    This parameter is used to distinguish between the two phases that this system can have. It is zero when we are in disordered (paramagnetic) state and takes a non-zero value in ordered (ferromagnetic) state.

    Monte-Carlo Simulation

    In this section, I’m going through the details of implementing the Monte-Carlo (MC) method for our Ising model. I apply the MC methods using Metropolis Algorithm to Ising model and extract physical parameters (Energy, Specific heat and Magnetization).

    Physical Model

    Lattice is a periodical structure of points that align one by one. 2D lattice can be plotted as:

    * * * * * * * *   
    * * * * * * * * 
    * * * * * * * *
    * * * * * * * *
    * * * * * * * *

    The points in lattice are called lattice points, nearest lattice points of point ^ are those lattice points denoted by (*) shown in the graph below:

    * * *(*)* * * *
    * *(*)^(*)* * *
    * * *(*)* * * *
    * * * * * * * *

    Each lattice point is denoted by a number i in the Hamiltonian.

    The expression for the Energy of the total system is

    H=Ji=0N1j=0N1(si,jsi,j+1+si,jsi+1,j) H=-J \sum_{i=0}^{N-1} \sum_{j=0}^{N-1}\left(s_{i, j} s_{i, j+1}+s_{i, j} s_{i+1, j}\right)

    * * * * * * * * 
    * * * * * * * *
    * * * * * * * * <-the i-th lattice point
    * * * * * * * *
    * * * * * * * *

    Periodical structure means that lattice point at(1,1) is the same as that at(1,9) if the lattice is 5 x 8. For example (1,1)<=>(6,1), (2,3)<=>(2,11). A 2D lattice can be any NxN_x by Ny N_y . The location (x,y)(x,y) here is another denotion of lattice point that is fundamentally same as i-th lattice point denotation above.

    * * * * * * * * 4
    * * * * * * * * 3
    * * * * * * * * 2
    * * * * * * * * 1
    1 2 3 4 5 6 7 8 

    Python Implementation


    1. Prepare some initial configurations of N spins.
    2. Flip spin of a lattice site chosen randomly
    3. Calculate the change in energy due to that
    4. If this change is negative, accept such move. If change is positive, accept it with probability exp^{-dE/kT}
    5. Repeat 2-4.
    6. Calculate Other parameters and plot them

    Python code for Ising model

    ##  Import needed python libraries 
    import matplotlib.pyplot as plt
    from numba import jit # for faster computation
    import numpy as np  # for array and matrix operations
    import random # for random number generation
    import time # for timing the code 
    from tqdm import tqdm # for progress bar
    #from tqdm.notebook import tqdm # for progress bar in jupyter notebook

    Firrst, we need to define the parameters of the system. We will use the following parameters:

    # Define parameters 
    B = 0; # Magnetic field strength
    L = 50; # Lattice size (width)
    s = np.random.choice([1,-1],size=(L,L)) # Begin with random spin sites with values (+1 or -1) for up or down spins. 
    n= 1000 * L**2 # number of MC sweeps 
    Temperature = np.arange(1.6,3.25,0.01) # Initlaize temperature range, takes form np.arange(start,stop,step)

    And after that we need to define the main function to calculate the energy of the system. We will use the following function:

    # Define functions to compute energy and magnetization
    Energy of the lattice calculations. 
    The energy here is simply the sum of interactions between spins divided by the total number of spins
    @jit(nopython=True, cache=True) # wonderful jit optimization compiler in its high performance mode
    def calcE(s):
        E = 0
        for i in range(L):
            for j in range(L):
                E += -dE(s,i,j)/2
        return E/L**2
    Calculate the Magnetization of a given configuration
    Magnetization is the sum of all spins divided by the total number of spins
    @jit(nopython=True, cache=True) 
    def calcM(s):
        m = np.abs(s.sum())
        return m/L**2

    The next step is to define the function to calculate the change in energy due to a spin flip. We will use the following function:

    # Calculate interaction energy between spins. Assume periodic boundaries.
    # Interaction energy will be the difference in energy due to flipping spin i,j 
    # (Example: 2*spin_value*neighboring_spins)
    @jit(nopython=True, cache=True) 
    def dE(s,i,j): # change in energy function
        if i == 0:
            t = s[L-1,j]  # periodic boundary (top)
            t = s[i-1,j]
        if i == L-1:
            b = s[0,j]  # periodic boundary (bottom)
            b = s[i+1,j]
        if j == 0:
            l = s[i,L-1]  # periodic boundary (left)
            l = s[i,j-1]
        if j == L-1:
            r = s[i,0]  # periodic boundary  (right)
            r = s[i,j+1]
        return 2*s[i,j]*(t+b+r+l)  # difference in energy is i,j is flipped

    Now we implement the Metropolis Algorithm to flip the spin of a lattice site chosen randomly. We will use the following function:

    # Monte-carlo sweep implementation
    @jit(nopython=True, cache=True) 
    def mc(s,Temp,n):   
        for m in range(n):
            i = random.randrange(L)  # choose random row
            j = random.randrange(L)  # choose random column
            ediff = dE(s,i,j)
            if ediff <= 0: # if the change in energy is negative
                s[i,j] = -s[i,j]  # accept move and flip spin
            elif random.random() < np.exp(-ediff/Temp): # if not accept it with probability exp^{-dU/kT}
                s[i,j] = -s[i,j]
        return s

    And finally we implement the main function to calculate the energy and magnetization of the system. We will use the following function:

    # Compute physical quantities
    @jit(nopython=True, cache=True) 
    def physics(s,T,n):
        En = 0
        En_sq = 0
        Mg = 0
        Mg_sq = 0
        for p in range(n):
            s = mc(s,T,1)
            E = calcE(s)
            M = calcM(s)
            En += E
            Mg += M
            En_sq += E*E
            Mg_sq += M*M
        En_avg = En/n
        mag = Mg/n
        CV = (En_sq/n-(En/n)**2)/(T**2)
        return En_avg, mag, CV

    And we initialize the arrays to store the values of the energy, magnetization and specific heat capacity. And we use time() function to calculate the time taken by the code to run. We will use the following code:

    # Inititalize magnetization, average energy and heat capacity
    mag = np.zeros(len(Temperature))
    En_avg = np.zeros(len(Temperature))
    CV = np.zeros(len(Temperature))
    start = time.time()

    Now we are ready to write for loop to sweep over the temperature range and calculate the energy, magnetization and specific heat capacity. We will use the following code:

    # Simulate at particular temperatures (T) and compute physical quantities
    for ind, T in enumerate(track(Temperature)):
        # Sweeps spins
        s = mc(s,T,n)
        # Compute physical quanitites with 1000 sweeps per spin at temperature T
        En_avg[ind], mag[ind], CV[ind] = physics(s,T,n)
    end = time.time()
    print("Time it took in seconds is = %s" % (end - start))
    time = (end - start)/60
    print('It took ' + str(time) + ' minutes to execute the code')

    Now let’s plot the results. We will use the following code that plots the energy, magnetization and specific heat capacity as a function of temperature:

    f = plt.figure(figsize=(18, 10)); # plot the calculated values
    sp =  f.add_subplot(2, 2, 1 );
    plt.plot(Temperature, En_avg, marker='.', color='IndianRed')
    plt.xlabel("Temperature (T)", fontsize=20);
    plt.ylabel("Energy ", fontsize=20);         plt.axis('tight');
    sp =  f.add_subplot(2, 2, 2 );
    plt.plot(Temperature, abs(mag), marker='.', color='RoyalBlue')
    plt.xlabel("Temperature (T)", fontsize=20); 
    plt.ylabel("Magnetization ", fontsize=20);   plt.axis('tight');
    sp =  f.add_subplot(2, 2, 3 );
    plt.plot(Temperature, CV, marker='.', color='IndianRed')
    plt.xlabel("Temperature (T)", fontsize=20);  
    plt.ylabel("Specific Heat ", fontsize=20);   plt.axis('tight');   
    plt.subplots_adjust(0.12, 0.11, 0.90, 0.81, 0.26, 0.56)
    plt.suptitle("Simulation of 2D Ising Model by Metropolis Algorithm\n" + "Lattice Dimension:" + str(L) + "X" + str(
        L) + "\n" + "External Magnetic Field(B)=" + str(B) + "\n" + "Metropolis Step=" + str(n)) # function to show the plots


    These are plots of the physical quantities for different MC steps.

    Lattice square representation of Ising model,


    You can run the Jupyter Notebook provided on Colab directly, or you can download from my project GitHub repository and run locally. Also, there is python script that you can run. I’m using Numba cache so that it produces cache folder (specific to machine CPU and configuration) so that you will need to produce yours by running it for one time and subsequent runs will be about twice faster.

    To run the script

    1. clone the repository
    git clone
    1. run the script
    python IsingModel/

    Appendix I

    I spent many nights working on this work, most of time I needed to optimize my code, I even tried to move to Matlab (last time I used it was like 5 years ago). But I learned a nice thing from my desire to optimize code speed. It is the usage of Numba’s JIT compiler. Read more about that here. I also instead of using multiple nested loops I dragged all these into just one. Imagine running 50x50 lattice simulation in my older codes for hours (one took 6 hours) vs 15 minutes for the current script. (On my Mac m1 Machine). Also, I made the code available on colab and can be accessed here (without many comments).

    Appendix II

    I have used this codebase as a benchmark test for python releases. I plan on maintain this test (available as repo on github). Initialy I used numpy-1.23.4 for all different version of python used in the benchmark. This is the latest version of numpy. I got the following results.


    The benchmarks are executed on Macbook pro Apple Silicon M1 version. The python version are installed using homebrew

    The arguments of the script are L(length of the lattice), n(number of Monte Carlo cycles).

    python 3.11

    time python3.11 10 10000
    Time it took in seconds is = 109.2954261302948
    python3.11 10 10000  110.82s user 0.79s system 100% cpu 1:51.38 total

    python 3.10

    time python3.10 10 10000
    Time it took in seconds is = 123.83638191223145
    python3.10 10 10000  124.80s user 0.69s system 101% cpu 2:04.06 total

    python 3.9

    time python3.9 10 10000
    Time it took in seconds is = 123.56476402282715
    python3.9 10 10000  123.94s user 1.08s system 101% cpu 2:03.67 total

    python 3.8

    time python3.8 10 10000
    Time it took in seconds is = 135.1741383075714
    python3.8 10 10000  137.05s user 1.00s system 100% cpu 2:17.24 total

    Other versions

    Python started supporting apple silicon starting from python3.8 and I try to make things as consistent as possible. Not to mention that numpy optimization for different version might be unaccounted difference.

    #Monte Carlo   #Ising Model   #Python   #Numba   #Metropolis Algorithm