Classical Monte-Carlo simulations of the Ising model on square lattice¶

The Ising model considers the following Hamiltonian

$\hat{\mathcal{H}} = -J\displaystyle\sum_{<i,j>}\sigma_i \sigma_j$,

in which spins $\sigma_i$ can take only two discrete values $\sigma_i \in \{-1,+1\}$. Below, we will consider spins on two-dimensional square lattice taking into account couplings $J$ only for the nearest neighbor sites. For $J > 0$ the ferromagnetic transition temperature is known from the analytic solution

$T_C = \frac{2J}{k_B \ln(1 + \sqrt{2})} \approx 2.27 J/k_B$ .

Ensemble averages of physical quantities such as magnetization $M$ can be calculated using the Monte-Carlo simulation technique. Below, we will use the Metropolis method to perform Monte-Carlo simulations on our model. The procedure can be summarized as follow:

  1. Start with a random configuration of spins in a large supercell (e.g. $50 \times 50$).
  2. Flip a spin at random site $i$ into a new configuration.
  3. Calculate the energy difference $\Delta E$ between the old and the new configuration.
  4. Generate a random number $0 < r < 1$.
  5. Accept the new configuration only if $r < e^{-\frac{\Delta E}{k_B T}}$.
  6. Go back to step 2.

At the end of each series of such spin flips (in fact, done in sequential order in this demonstration) we will be calculating average magnetization $M$ over the system. This quantity is plotted at the end of the simulation, as well as its average $<M>$ over the entire simulation.

In [1]:
import numpy as np
import matplotlib as mat
import matplotlib.pyplot as plt
import ipywidgets as ipw

from ipywidgets import interact
from random import choice
%matplotlib inline

First define two functions that, when being called, map a single dimension array to the lattice and vice versa.

In [2]:
# Define a function to map a 1D array to a 2D/3D lattice, only works on lattice of equal-lateral size in Cartesian system
def map(index, size, dimension): 
    if dimension == 1:
        return index;
    if dimension == 2:
        x = index%size;
        y = np.floor(index/size);
       # return vp.vector(x,y,0);
        return [x, y];
In [3]:
# Define a function to map the 3D vector back into the 1D array index
def antimap(array, size, dimension):
    if dimension == 1:
        return array[0];
    if dimension == 2:
        return array[0] + array[1]*size;

Then define another function that, when being called, returns an array that includes the indices of the current site's nearest neighbours

In [4]:
# Define a function to find the index of the nearest neighbours
def nNbr(index, size, dimension):
    if dimension == 1:
        return [(index+size-1)%size, (index+size+1)%size]; # apply periodic condition
    
    if dimension == 2:
        pos = map(index,size,dimension);
        north = np.mod(np.add(pos,[0,size+1]), size);
        south = np.mod(np.add(pos,[0,size-1]), size);
        east = np.mod(np.add(pos,[size+1,0]), size);
        west = np.mod(np.add(pos,[size-1,0]),size);
        return [antimap(north, size, 2), antimap(south, size, 2), antimap(east, size, 2), antimap(west, size, 2)];

Finally, the data is output frame by frame using interactive figures.

In [5]:
# Define the function that updates the plot
def update(step):
    fig = plt.figure(figsize=(15,9)); # Initialize a figure
    figa = fig.add_subplot(121);
    plt.subplots_adjust(left=0, bottom=0.25)
    Ising = figa.scatter(lattice[:,0],lattice[:,1], s=90, c = spinC[step,:], marker = 's', edgecolor = None);
    figb = fig.add_subplot(122);
    figb.plot(tt, mag,'o-', c='k');
    figb.plot(tt[step], mag[step],'o',c='r');
    plt.subplots_adjust(left=0,bottom=0.25);
    plt.xlabel('Step');
    plt.ylabel('Magnetization M');
    plt.show()

Define the parameters defining your system and conditions such as temperature, dimensionality, size of the lattice, as well as Ising coupling strength measured in the units of Boltzmann constant.

In [6]:
# Set up the system parameters
temp = 1.0; # temperature in [K] 
L = 50; # Define the lattice size
dims = 2; # Define the dimensionality of the lattice (up to 2D)
N = L**dims; # Compute the total number of the sites
kB = 8.617E-5 # Boltzmann constant in [eV/K]
J = -1.0*kB; # Define the coupling strength between spins
beta = 1/(temp*kB);

Initialize the lattice, use spinV to store spin configurations, and spinC to store the corresponding coloring. Iterate through the entire lattice.

In [7]:
# Set up array containers for data manipulation and initial conditions
#lattice = np.empty((N*N,3),int); # Create a container for the lattice sites
lattice = np.empty((N,dims),float); # Create a container for site positions in the lattice
spinV = np.empty(N,int); # Create a container for the spins for all the lattice sites
spinC = np.chararray(N); # Create a container to store color values (spin visualization)
nIdx = np.empty((N,dims*2), np.intp); # Create a container to store the indices of the nearest neighbours of each site
probs = np.empty(N,float); # Create a container to store in each iteration spin-flip probability
for ii in range(N):
    lattice[ii] = np.subtract(map(ii,L,dims),L/2);
    spinV[ii] = 1; # Create a container to store the values of the spins
    spinC[ii] = 'r';
    nIdx[ii] = nNbr(ii, L, dims);
mag = [sum(spinV)/N]; # Initializat a variable for the total magnitization of the lattice
In [8]:
# Optional: Randomization of the initial condition
for ii in range(int(N/2)): # Randomly choose half of the lattice to flip the spin 
    kk = choice(range(N));
    spinV[kk] = -1;
    spinC[kk] = 'k'; # color all the spin down as black
mag = [sum(spinV)/N]; # Calculate the current total magnetization
spinC = np.expand_dims(spinC, axis = 0);

Set the maximum iteration limit, and start the simulation using the Metropolis algorithm. In the current case, we use arrays of dimension one to store all the spins (spinV), the corresponding coloring (spinC) as well as the random number generated (ruler) for the Metropolis algorithm. At the very end after all iterations having been carried out, it calculates the average magnetization over entire simulation excluding the first 50 steps (defined as variable "exlu_lim").

In [9]:
# !!This step will take a while depending on the number of iterations and the speed of the computer!!
# Pre-define the total iteration number and pre-calculate each step to save time when plotting
MaxIteration = 200; # Start the iteration count
tt = range(MaxIteration+1); # Initiate a array for plotting the magnetization
ruler = np.random.rand(MaxIteration,len(spinV)); # Initiate a matrix to store pre-calculated random numbers for all iterations
tempC = np.empty(N,str); # create an empty array to store colors (spins) at each iteration to update the lattice
for step in range(MaxIteration):
    tempSpin = spinV; # Create a temporary container to store updated spins
    for ii in range(len(spinV)):
        nIdx[ii] = nNbr(ii, L, dims);
        probs[ii] = np.exp(beta*J*sum(spinV[nIdx[ii]])*spinV[ii])/np.exp(-beta*J*sum(spinV[nIdx[ii]])*spinV[ii]);
        tempSpin[ii] = -spinV[ii] if probs[ii] > ruler[step,ii] else spinV[ii];
        tempC[ii] = 'r' if tempSpin[ii] > 0 else 'k';
    spinC = np.append(spinC, np.expand_dims(tempC,axis = 0), axis = 0); 
    spinV = tempSpin;
    mag.append(sum(tempSpin)/N);
#Output the interactive plot and the final magnetization
exlu_lim = 50;
print('Average magnetization <M>: ' + f"{np.mean(mag[exlu_lim:]):.3f}"); 

Average magnetization <M>: 0.999

Finally, output the spin map

In [10]:
# Start the interactive plot to visualizate the (pseudo-)magnetization process
interact(update,step=ipw.IntSlider(min=0, max=MaxIteration, step=1, value=0));

interactive(children=(IntSlider(value=0, description='step', max=200), Output()), _dom_classes=('widget-intera…
In [ ]: