Ising Model - 2D : Monte Carlo Simulation

First of all we import the required libraries:

In [28]:
import matplotlib.pyplot as plt
import math as math
import random as random
import pylab as pylab
import numpy as np
import seaborn as sns
sns.set()

Global parameters

Then we set some global variables:

In [29]:
Kb = 1.0
JJ = 1.0
mu = 1.0

Cold Start and Hot Start

Before begining simulation, if we select all spin up or all spin down then it is called cold configuration.

In [30]:
def cold_start(L):
    U = [[1.0 for k in range(L)]for l in range(L)]
    return U   

We can make a plot of this cold configuration by the help of Uplotter :

In [31]:
CU = cold_start(10)
sns.heatmap(CU,annot =True)
Out[31]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fda178deeb8>

if we select spin randomly up or down then it is called cold configuration

In [32]:
def hot_start(L): 
    U = [[0.0 for i in range(L)]for j in range(L)]
    for i in range(L):
          for j in range(L):
            t = random.sample([-1,1],1)
            U[i][j] = t[0]  
    return U 

Similarly we can plot thishot configuration as well.

In [33]:
HU = hot_start(10)
sns.heatmap(HU,annot =True)
Out[33]:
<matplotlib.axes._subplots.AxesSubplot at 0x7fda17899ac8>

Hamiltonian

Hamiltonian of L by L lattice is given by

$ H = - J \sum_{i\neq j} S_{i}S_{j} - \mu B \sum_{i} S_{i} $

In [34]:
def Hamiltonian(U,B):
    H = 0.0
    L = len(U)
    for i in range(L):
        for j in range(L):
                
            ni = 0;nj =0;
            if i == 0: ni = L
            if j == 0: nj = L
           
            H = H -  0.5*JJ*U[i][j]*(U[i][(j-1)+nj] \
                        + U[(i+1)%L][j] \
                        + U[i][(j+1)%L] \
                        + U[(i-1)+ni][j])  - mu*B*U[i][j]
    return H

One can check what are hamiltonian for hot and cold start:

In [35]:
Hamiltonian(CU, B=1.0)
Out[35]:
-300.0
In [36]:
Hamiltonian(HU, B=0)
Out[36]:
-12.0

Mangnetization

One can calculate magnetization by simply taking average over all spins:

In [37]:
def magnetization(U):
    return np.array(U).sum()/float(len(U*len(U)))          
In [38]:
magnetization(HU)
Out[38]:
-0.04
In [39]:
magnetization(CU)
Out[39]:
1.0

Spin Flipper

In [40]:
def ld(k,L):
    if k == 0:
        return L
    else: return 0
In [41]:
def spin_flipper(U,B, printkey):
        L = len(U)
        
        i = random.randint(0, L-1) 
        j = random.randint(0, L-1)
       
        if printkey ==1:print("flipped at", i,j)
            
        U[i][j] = -U[i][j]
       
           
        dH =  -2.0*JJ*U[i][j]*(U[i][(j-1)+ld(j,L)] \
                            + U[(i+1)%L][j] \
                            + U[i][(j+1)%L] \
                            + U[(i-1)+ld(i,L)][j])  - mu*B*U[i][j]   
        
        return U,dH,i,j 

Thermalization

In [44]:
def Thermalization(U,T,nrun,printkey) :
        M = [0.0 for k in range(nrun)]
        irun = 0       
        B = 0.0
        HH = [0.0 for k in range(nrun)]
       
        while irun < nrun:
                    V = U
                    w = magnetization(U)
                    h = Hamiltonian(U,B)
                    
                    U,dH,p,q = spin_flipper(U,B,printkey)
                   
                   
                    if dH < 0:
                        
                        if printkey ==1: print(irun, "E decreased! You are accepted !",dH)
                            
                        M[irun] = magnetization(U)
                        HH[irun] = Hamiltonian(U,B)
                       
                    else:
                        
                        if printkey ==1:print(irun, "E increased!",dH)
                            
                        frac = math.exp(-dH/(Kb*T))
                        b = random.uniform(0.0,1.0)
                            
                        if printkey ==1:print("frac =",frac,"b=",b,"dH = ",dH)
                                
                        if  b < frac:
                                    
                            if printkey ==1:print(irun, " You Lucky!")
                                
                            M[irun] = magnetization(U)
                            HH[irun] = Hamiltonian(U,B)
                             
                        else:
                            if printkey ==1: print(irun, "Loser!")
                            if printkey ==1: print("spin restablished at",p,q)
                            U[p][q] = -U[p][q]
                            M[irun] = w 
                            HH[irun] = h
                           
                            
                    for i in range(L):
                        for j in range(L):
                            if U[i][j] != V[i][i]:
                                 if printkey ==1: print("Warning!spin is changed!", i,j) 
                            
                    
                    
                    
                    if printkey ==2 : print(irun, M[irun])   
                    irun = irun +1
        return M,U,HH

Lets print out some measurements of m

In [45]:
nrun = 20
T = 1.0
L = 10
U = cold_start(L)  
M,U,HH = Thermalization(U,T,nrun, 1)
flipped at 4 2
0 E increased! 8.0
frac = 0.00033546262790251185 b= 0.5497168174143898 dH =  8.0
0 Loser!
spin restablished at 4 2
flipped at 9 8
1 E increased! 8.0
frac = 0.00033546262790251185 b= 0.8708650724871093 dH =  8.0
1 Loser!
spin restablished at 9 8
flipped at 5 5
2 E increased! 8.0
frac = 0.00033546262790251185 b= 0.6435985276056473 dH =  8.0
2 Loser!
spin restablished at 5 5
flipped at 6 5
3 E increased! 8.0
frac = 0.00033546262790251185 b= 0.24779278590316134 dH =  8.0
3 Loser!
spin restablished at 6 5
flipped at 1 7
4 E increased! 8.0
frac = 0.00033546262790251185 b= 0.13621648458611013 dH =  8.0
4 Loser!
spin restablished at 1 7
flipped at 3 8
5 E increased! 8.0
frac = 0.00033546262790251185 b= 0.12510920623692545 dH =  8.0
5 Loser!
spin restablished at 3 8
flipped at 7 0
6 E increased! 8.0
frac = 0.00033546262790251185 b= 0.39719764293621906 dH =  8.0
6 Loser!
spin restablished at 7 0
flipped at 0 2
7 E increased! 8.0
frac = 0.00033546262790251185 b= 0.40943034843557125 dH =  8.0
7 Loser!
spin restablished at 0 2
flipped at 0 5
8 E increased! 8.0
frac = 0.00033546262790251185 b= 0.17964706089284155 dH =  8.0
8 Loser!
spin restablished at 0 5
flipped at 4 6
9 E increased! 8.0
frac = 0.00033546262790251185 b= 0.6938205246000504 dH =  8.0
9 Loser!
spin restablished at 4 6
flipped at 9 2
10 E increased! 8.0
frac = 0.00033546262790251185 b= 0.01488254643068998 dH =  8.0
10 Loser!
spin restablished at 9 2
flipped at 9 6
11 E increased! 8.0
frac = 0.00033546262790251185 b= 0.24448004198505613 dH =  8.0
11 Loser!
spin restablished at 9 6
flipped at 2 4
12 E increased! 8.0
frac = 0.00033546262790251185 b= 0.5531985218511174 dH =  8.0
12 Loser!
spin restablished at 2 4
flipped at 6 9
13 E increased! 8.0
frac = 0.00033546262790251185 b= 0.5381882188668542 dH =  8.0
13 Loser!
spin restablished at 6 9
flipped at 4 5
14 E increased! 8.0
frac = 0.00033546262790251185 b= 0.9961179807143019 dH =  8.0
14 Loser!
spin restablished at 4 5
flipped at 5 0
15 E increased! 8.0
frac = 0.00033546262790251185 b= 0.3269242304046608 dH =  8.0
15 Loser!
spin restablished at 5 0
flipped at 9 6
16 E increased! 8.0
frac = 0.00033546262790251185 b= 0.6448877624644563 dH =  8.0
16 Loser!
spin restablished at 9 6
flipped at 3 0
17 E increased! 8.0
frac = 0.00033546262790251185 b= 0.8958953479801371 dH =  8.0
17 Loser!
spin restablished at 3 0
flipped at 4 0
18 E increased! 8.0
frac = 0.00033546262790251185 b= 0.8491684424879403 dH =  8.0
18 Loser!
spin restablished at 4 0
flipped at 1 1
19 E increased! 8.0
frac = 0.00033546262790251185 b= 0.2378031735888625 dH =  8.0
19 Loser!
spin restablished at 1 1
In [46]:
L = 20
nrun = 10000
T = 2.4
U = cold_start(L) 
M,U,HH = Thermalization(U,T,nrun, 0)
X = np.arange(0,len(M),1)
plt.figure(figsize = [15,6])
plt.plot(X,M,"-")
plt.show() 

We can plot both run with hot and cold start together:

In [47]:
L =20
nrun = 10000
T = 5.0
U1 = cold_start(L) 
U2 = hot_start(L)
M1,U1,HH = Thermalization(U1,T,nrun,0)
M2,U2,HH = Thermalization(U2,T,nrun,0)
X = np.arange(0,len(M1),1)

plt.figure(figsize = [15,6])
plt.plot(X,M1,"-")
plt.plot(X,M2,"-")
plt.show()

Phase Transition

In [49]:
L = 32
nrun = 10000
Tn = 100
avm = []
stdh = []
KT = []

for t in range(1,Tn+1):
        T = 0.1*t
        KT.append(T)
        U = cold_start(L)
        M,U,HH = Thermalization(U,T,nrun,0)
        nM = M[1000:nrun-1]
        nH = HH[1000:nrun-1]
        stdh.append(np.std(nH))
        avm.append(np.mean(nM))
In [50]:
plt.figure(figsize = [10,8])
plt.scatter(KT,avm)
plt.xlabel("Temperature")
plt.ylabel("Average magnetization")
plt.show() 

Specific Heat Capacity

In [52]:
plt.figure(figsize = [10,8])
plt.scatter(KT,stdh)
plt.xlabel("Temperature")
plt.ylabel("Specific Heat Capacity")
plt.show() 

Effect of Lattice Size

In [54]:
S = []
for L in [4,8,12,16,32]:
    nrun = 1000
    Tn = 100
    avm = []
    stdh = []
    KT = []

    for t in range(1,Tn+1):
        T = 0.1*t
        KT.append(T)
        U = cold_start(L)
        M,U,HH = Thermalization(U,T,nrun,0)
        nH = HH[100:nrun-1]
        stdh.append(np.std(nH))
    S.append(max(stdh))
In [56]:
Tx = [4,8,12,16,32]
plt.figure(figsize = [10,8])
plt.scatter(Tx,S)
plt.plot(Tx,S)
plt.xlabel("Lattice size")
plt.ylabel("C_max")
plt.show() 
In [ ]: