Ising 2D ferromagnetic material

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 = 2
H*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./(N
2) - 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 = E
E
Mmean=sum(M)/numiter
Emean=sum(E)/numiter
Mmean2=sum(MM)/numiter
Emean2=sum(E
E)/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’);

In Python, indentation is significant.

The body of the function must be indented. (We normally use four spaces per indent.) That is how the Python interpreter tells where the body begins and ends.

After a def, the body starts with an indent, and ends when the indentation level pops back.

# Indent level = 0, so not inside a function body.
# def tells the interpreter to expect an indented block
def function():
    # indent level = 1
    body of the function

# and back to indent level = 0 again
# so no longer part of the function body

The return statement must be inside the function, so it must be indented.

You must write:

def ising(numiter,N,T,H):
    "the ising function"
    J = 1.
    k = 1.
    M = np.array([])
    E = np.array([])
    # ... lots more code ...
    Mmean2 = sum(M*M)/numiter
    Emean2 = sum(E*E)/numiter
    # return must be indented
    return Mmean, Emean, Mmean2, Emean2

# Here we are outside the function again.

This is sometimes called “the offside rule”:

These are the commands which expect the indentation level to increase by one:

  • class and def

  • for and while

  • if, elif and else

  • try, except and finally

  • with

  • match, case in Python 3.10.

At the end of each block, you decrease the indentation by one level again.

By the way, you don’t need semicolons at the end of each line in Python. Using them is legal, but is a strong sign that you don’t know the language :slight_smile: