# Why the magnetisation shows abrupt behaviour for this 3D ising spin system

I am trying to simulate a 3D Ising spin system (+1 & -1) using Monte Carlo Metropolis Algorithm. I want to get different physical quantities from this simulation like magnetization, Average Energy, Susceptibility, Specific Heat. I have used exactly similar code for 2 dimension that gives perfect result but for 3D it fails. The magnetization curve is quite noisy. In spite of increasing trials it always shows sudden valleys. I have made all the required modifications for 3D but can't find out what is going wrong. The specific heat also does not show the discontinuity in the right manner.

The following code simulates for a 3d Ising spin system 5x5x5 with 1250 Monte Carlo steps (excluding 2000 steps to let the system reach equilibrium). Currently, it only calculates magnetization.Below is the plot I am currently getting from the code.



# cost            change in energy
# J               spin coupling constant, will be normalised
# spin_states     3D array, the lattice
# M               total magnetisation of the lattice
# m               magnetisation per spin { m = M / (N * N * N) }
# N               size of the lattice

from numpy.random import rand

import numpy as np
import matplotlib.pyplot as plt
import time     # for timing

def initialise(N):
''' generates a random spin spin_statesuration for initial condition'''
spin_states = np.random.choice([1, -1], size=(N, N,N))
#spin_states=np.ones([N,N,N])
return spin_states

def calcMag(spin_states):
'''Magnetization of a given spin_statesuration'''
mag = abs(np.sum(spin_states))
return mag

def calcEnergy(spin_states):
'''Energy of a given spin_statesuration'''
energy = 0
for i in range(len(spin_states)):
for j in range(len(spin_states)):
for k in range(len(spin_states)) :
s = spin_states[i,j,k]
energy += -s*find_neighbours(spin_states,N,i,j,k)
return energy/8.

def find_neighbours(spin_states,N,x,y,z):
left   =spin_states[x,(y-1)%N,z]
right  =spin_states[x,(y+1)%N,z]
top    =spin_states[(x-1)%N,y,z]
bottom =spin_states[(x+1)%N,y,z]
front  =spin_states[x,y,(z+1)%N]
back   =spin_states[x,y,(z-1)%N]

tot_spin=left+right+top+bottom+front+back

return (tot_spin)

def mcmove(spin_states, beta):
'''Monte Carlo move using Metropolis algorithm '''
x = np.random.randint(0, N)
y = np.random.randint(0, N)
z = np.random.randint(0, N)
s = spin_states[x, y, z]

cost = 2*s*find_neighbours(spin_states,N,x,y,z)
if cost < 0:
s *= -1
elif rand() < np.exp(-cost*beta):
s *= -1
spin_states[x, y,z] = s
return spin_states

nt      = 80         #  number of temperature points
N       = 5        #  size of the lattice, size x size
eqSteps = 2000       #  number of MC sweeps for equilibration
mcSteps = 1250       #  number of MC sweeps for calculation

spin_states = initialise(N)

m = calcMag(spin_states)
print ("m =", m)

#k = mcSteps * N * N * N

dataM = []

T       = np.linspace(1, 7, nt)
E,M,C,X = np.zeros(nt), np.zeros(nt), np.zeros(nt), np.zeros(nt)
n1, n2  = 1.0/(mcSteps*N*N*N), 1.0/(mcSteps*mcSteps*N*N*N)
for tt in range(nt):
E1 = M1 = E2 = M2 = 0
spin_states = initialise(N)
iT=1.0/T[tt]; iT2=iT*iT;   # inverse temperature for beta

for i in range(eqSteps):         # equilibrate
mcmove(spin_states, iT)           # Monte Carlo moves

for i in range(mcSteps):
mcmove(spin_states, iT)
Ene = calcEnergy(spin_states)     # calculate the energy
Mag = calcMag(spin_states)        # calculate the magnetisation

E1 = E1 + Ene
M1 = M1 + Mag
M2 = M2 + Mag*Mag
E2 = E2 + Ene*Ene

E[tt] = n1*E1
M[tt] = n1*M1
C[tt] = (n1*E2 - n2*E1*E1)*iT2
X[tt] = (n1*M2 - n2*M1*M1)*iT

print(M)
plt.plot(T,M)
plt.xlabel('temperature')
plt.ylabel('magnetisation')
plt.title('3D Ising Model')
plt.grid(True)
plt.show()

plt.figure()
plt.plot(T,C,'+')



Please can anybody help me to find out what is going wrong? As far as I can think, the metropolis algorithm is really doing not much change. I have done the simulation for diff. lattice sizes with diff no. of steps. It never changes much from the initial value. I am new to simulation and this is my second simulation project. Any help will be much appreciated. Thank you in advance.

• Hi and welcome SciComp SE. Can you add some plots of the magnetisation as a function of temperature and add some more details about what you are simulating (size of the lattice, number of steps etc.) Mar 21 '20 at 8:22
• @Ir1985 Thanks for your help. I have added the plot and other details as you suggested. Please take a look when you are free. Mar 21 '20 at 13:00

Your lattice consists of 5 x 5 x 5 = 125 spins, so your number of Montecarlo steps to reach equilibrium should be >> 125, because you randomly picking a site and flipping it, so random numbers should uniformly generated so that it will cover whole lattice. For much finer measurement of thermodynamic quantities, you should take more number of points between the temperature range and more number of Montecarlo steps >> 25.

So I modified your code a bit, instead of picking a site randomly and flipping, I swept through each site starting row and column wise.

def mcmove(spin_states, beta):
'''Monte Carlo move using Metropolis algorithm '''
for x in range(len(spin_states)):
for y in range(len(spin_states)):
for z in range(len(spin_states)):
s = spin_states[x,y,z]
cost = 2*s*find_neighbours(spin_states,N,x,y,z)
if cost < 0:
s *= -1
elif rand() < np.exp(-cost*beta):
s *= -1
spin_states[x, y,z] = s
return spin_states


I ran the simulations on a 16 x 16 x 16 lattice, temperature ranging from 0.1 - 8.0 with 500 points in-between, with 5000 steps to reach equilibrium and 5000 steps to take measurement (you can get similar results with less number of Montecarlo steps, say 100 + 100, since it will sweep through whole lattice). These are my results