# 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 [ ]: