There are two (possible) problems to this code I could use help with.
I have written the following code to calculate the average energy, magnetisation and specific heat of a material with lattice sizes 10x10, 50x50 and 100x100.
The graphs shown have significant fluctuations, which I feel is for two reasons.
- Firstly, I can only seem to run the program with around 100 'relax' steps before performing the actual calculation. Answers to similar questions on here suggest doing around 1000 steps, but my program takes extremely long to do this.
Secondly, there was a hint given to us by our professor which is the following:
"If you start the loop over the temperature with the lowest temperature, use as initial state the configuration with all spins up. When you change the temperature, use the last spin configuration as the first configuration for the next Nwait MMC steps. – This will reduce the fluctuations on the data"
I've been trying to implement this by for example, putting
if m == 1:
s = np.ones(10,10)
else:
config = initialstate(N)
in the for loop with N = 10x10.
But that only partially solves it, as I'm not using the last spin configuration as the first configuration for the next relax steps.
The code:
import numpy as np
from numpy.random import rand
import matplotlib.pyplot as plt
#----------------------------------------------------------------------
## BLOCK OF FUNCTIONS USED IN THE MAIN CODE
#----------------------------------------------------------------------
def initialstate(N):
''' generates a random spin configuration for initial condition'''
state = 2*np.random.randint(2, size=(N,N),dtype='int64')-1
return state
def mcmove(config, beta):
'''Monte Carlo move using Metropolis algorithm '''
for i in range(N):
for j in range(N):
a = np.random.randint(0, N,dtype='int64')
b = np.random.randint(0, N,dtype='int64')
s = config[a, b]
nb = config[(a+1)%N,b] + config[a,(b+1)%N] + config[(a-1)%N,b] + config[a,(b-1)%N]
energy = 2*s*nb
if energy < 0:
s *= -1
elif rand() < np.exp(-energy*beta):
s *= -1
config[a, b] = s
return config
def calcEnergy(config):
'''Energy of a given configuration'''
energy = 0
for i in range(len(config)):
for j in range(len(config)):
S = config[i,j]
nb = config[(i+1)%N, j] + config[i,(j+1)%N] + config[(i-1)%N, j] + config[i,(j-1)%N]
energy += -nb*S
return energy/4.
def calcMag(config):
'''Magnetization of a given configuration'''
mag = np.sum(config)
return mag
#----------------------------------------------------------------------
## MAIN PART OF THE CODE
#----------------------------------------------------------------------
nt = 20 # number of temperature points
eqSteps = 100 # number of MC sweeps for equilibration
mcSteps = 100 # number of MC sweeps for calculation
T = np.linspace(1, 4, nt) #temperature
Energy_10 = np.zeros(nt,dtype='int64')
Magnetization_10 = np.zeros(nt,dtype='int64')
SpecificHeat_10 = np.zeros(nt,dtype='int64')
Energy_20 = np.zeros(nt,dtype='int64')
Magnetization_20 = np.zeros(nt,dtype='int64')
SpecificHeat_20 = np.zeros(nt,dtype='int64')
Energy_30 = np.zeros(nt,dtype='int64')
Magnetization_30 = np.zeros(nt,dtype='int64')
SpecificHeat_30 = np.zeros(nt,dtype='int64')
for m in range(len(T)): #for loop for 10x10 lattice
E1 = M1 = 0
N = 10 # size of the lattice, N x N
config = initialstate(N)
for i in range(eqSteps):
mcmove(config, 1.0/T[m])
for i in range(mcSteps):
mcmove(config, 1.0/T[m]) # monte carlo moves
Ene = calcEnergy(config) # calculate the energy
Mag = calcMag(config) # calculate the magnetisation
E1 = E1 + Ene
M1 = M1 + Mag
Energy_10[m] = E1/(N**2)
Magnetization_10[m] = M1/(N**2)
SpecificHeat_10[m] = (((T[m]**2)*(E1/(N**2))*(E1)-(E1/(N**2))**2))/(N**2)
for m in range(len(T)): #for loop for 50x50 lattice
E1 = M1 = 0
N = 50
config = initialstate(N)
for i in range(eqSteps):
mcmove(config, 1.0/T[m])
for i in range(mcSteps):
mcmove(config, 1.0/T[m]) # monte carlo moves
Ene = calcEnergy(config) # calculate the energy
Mag = calcMag(config) # calculate the magnetisation
E1 = E1 + Ene
M1 = M1 + Mag
Energy_20[m] = E1/(N**2)
Magnetization_20[m] = M1/(N**2)
SpecificHeat_20[m] = ((T[m]**2)*(E1/(N**2))*(E1)-(E1/(N**2))**2)/(N**2)
for m in range(len(T)): #for loop for 100x100 lattice
E1 = M1 = 0
N = 100
config = initialstate(N)
for i in range(eqSteps):
mcmove(config, 1.0/T[m])
for i in range(mcSteps):
mcmove(config, 1.0/T[m]) # monte carlo moves
Ene = calcEnergy(config) # calculate the energy
Mag = calcMag(config) # calculate the magnetisation
E1 = E1 + Ene
M1 = M1 + Mag
Energy_30[m] = E1/(N**2)
Magnetization_30[m] = M1/(N**2)
SpecificHeat_30[m] = ((T[m]**2)*(E1/(N**2))*(E1)-(E1/(N**2))**2)/(N**2)
#----------------------------------------------------------------------
# Plot the Energy, Magnetization and Specific Heat
#----------------------------------------------------------------------
f = plt.figure(figsize=(18, 10), dpi=80, facecolor='w', edgecolor='k');
sp = f.add_subplot(2, 2, 1 );
plt.plot(T, Energy_10, 'g',label = 'N=10')
plt.plot(T, Energy_20, 'r',label = 'N=50')
plt.plot(T, Energy_30, 'b',label = 'N=100')
plt.legend()
plt.xlabel("Temperature (T)", fontsize=20);
plt.ylabel("Energy ", fontsize=20);
plt.show()
sp = f.add_subplot(2, 2, 2 );
plt.plot(T, abs(Magnetization_10), 'g',label = 'N=10')
plt.plot(T, abs(Magnetization_20), 'r',label = 'N=50')
plt.plot(T, abs(Magnetization_30), 'b',label = 'N=100')
plt.legend()
plt.xlabel("Temperature (T)", fontsize=20);
plt.ylabel("Magnetization ", fontsize=20);
plt.show()
sp = f.add_subplot(2, 2, 3 );
plt.plot(T, SpecificHeat_10, 'g',label = 'N=10')
plt.plot(T, SpecificHeat_20, 'r',label = 'N=50')
plt.plot(T, SpecificHeat_30, 'b',label = 'N=100')
plt.legend()
plt.xlabel("Temperature (T)", fontsize=20);
plt.ylabel("Specific Heat ", fontsize=20);
plt.show()
#----------------------------------------------------------------------
