Dear all,
I am doing project on modelling 2D ferromagnetic material using Ising energy model, I have some problems with the code. I introduced ising function (this is a requirement) but it says that return operator is outside the function. I am sharing the code. Any help will be appreciated. Thank you!
import numpy as np
from math import *
import matplotlib.pyplot as plt
from numpy import random
N = 25
numiter = 10
H = 1
#============================================FUNCTION ISING========================================
def ising(numiter,N,T,H):
“the ising function”
J = 1.;
k = 1.;
M=np.array();
E=np.array();
Generate a random initial configuration
grid = (random.rand(N,N) > 0.5)*2 - 1;
#Original configuration will all spins up
#grid = -np.ones(s=(N,N));
#Evolve the system for a fixed number of steps
plt.figure(1)
plt.imshow((grid+1)*128,cmap=‘bone’,interpolation=“nearest”);
plt.show()
for i in range(numiter+1):
#---------------------------------------------------------------------
T=3.+0.25*(i-1);
x1=np.roll(grid, +1, axis=1)
x2=np.roll(grid, -1, axis=1)
x3=np.roll(grid, +1, axis=0)
x4=np.roll(grid, -1, axis=0)
Sum = x1 + x2 + x3 + x4
DeltaE1 = 2JgridSum #array that gives an energy variation for flipping
DeltaE2 = 2H*grid
E = 0
for j in range(len(grid)):
for k in range(len(grid)):
if DeltaE1[j][k] < 0:
grid[i][j] *= -1
elif np.exp(-1./T*DeltaE1[j][k]) > np.random.rand():
grid[i][j] *= -1
E += -Sum*grid
M = np.sum(grid)
E = E/2.
#-------------------------------------------------------------------
#Store values of moment per site and energy per site
M=np.append(M,sum(sum(grid))/N2);
E=np.append(E,-sum(sum(DeltaE1))/4./(N2) - sum(sum(DeltaE2))/2./N**2);
Display the current state of the system (optional)
#figure(1)
#if (i%10000 == 0):
plt.subplot(221)
plt.xlabel(“i= %7d, T = %0.2f, M = %0.2f, E = %0.2f” %(i, T, M[i], E[i]))
plt.imshow((grid+1)*128,cmap=‘bone’,interpolation=“nearest”);
plt.subplot(222)
tgrid = range(1,i+1);
plt.plot(tgrid, M[:-1],‘bo’)
plt.xlabel(“iteration”)
plt.ylabel(“Instantaneaous Momentum”)
plt.show()
#plt.figure(3)
plt.subplot(223)
tgrid = range(1,i+1);
plt.plot(tgrid, E[:-1],‘go’)
plt.xlabel(“iteration”)
plt.ylabel(“Instantaneaous Energy”)
plt.show()
plt.pause(2)
plt.pause(0.000000001)
M2 = MM
E2 = EE
Mmean=sum(M)/numiter
Emean=sum(E)/numiter
Mmean2=sum(MM)/numiter
Emean2=sum(EE)/numiter
return Mmean, Emean, Mmean2, Emean2
#===============================================================================
#Initial Configuration
n_grid = 50; # Size of grid
numiter = 80000;
tgrid = range(1,numiter+1);
Ms = np.array();
Ts = np.array();
Es = np.array();
chi=np.array();
cv=np.array();
Hs=np.array();
#Mmean=np.array();
#Emean=np.array();
#Mmean2=np.array();
#Emean2=np.array();
for i in range(1,5):
H=0.;
T=3.+0.25*(i-1);
[Mmean, Emean,Mmean2,Emean2] = ising(numiter, n_grid, T,H);
plt.close(1);
mychi = 1/T*(Mmean2 - MmeanMmean);
mycv = 1/T**2(Emean2 - Emean*Emean);
Ms=np.append(Ms,abs(Mmean));
Es=np.append(Es,Emean);
chi=np.append(chi,mychi);
cv=np.append(cv,mycv);
Ts=np.append(Ts,T)
plt.figure(2)
plt.subplot(221)
plt.plot(Ts, Es, ‘ro’);
plt.ylabel(‘energy per site’);
plt.xlabel(‘Temperature’);
plt.subplot(222)
plt.plot(Ts, Ms, ‘bo’);
plt.ylabel(‘magnetization per site’);
plt.xlabel(‘Temperature’);
plt.subplot(223)
plt.plot(Ts, chi, ‘go’);
plt.ylabel(‘magnetization susceptibility per site’);
plt.xlabel(‘Temperature’);
plt.subplot(224)
plt.plot(Ts, cv, ‘mo’);
plt.ylabel(‘specific heat per site’);
plt.xlabel(‘Temperature’);
plt.savefig(“results/ex1.pdf”) # saving the result as a pdf
#fout1 = open(“output.dat”,“w”)
myrestable = np.vstack((Ts, Es, chi,cv)).T
np.savetxt(“results/output.dat”, myrestable, delimiter="\t")
#pbaspect([2 1 1]);
#print(gcf, ‘-depsc2’, ‘ising-energy’);