用蒙特卡洛法求解二维Ising模型的平均能量出现了一些问题

首先这个图像就完全不对
附上代码
import numpy as np
import random
def Ising(T,N,J=1,B=0):
S = np.ones((N,N),dtype=int)
for i in range(100*N**2):
pos_x = np.random.randint(0,N)
pos_y = np.random.randint(0,N)
dE = get_dE(N,S,pos_x,pos_y,J,B)
E_tot = Energy(N,S,pos_x,pos_y,J,B)
if dE<0:
S[pos_x,pos_y] = -S[pos_x,pos_y] E_tot = Energy(N,S,pos_x,pos_y,J,B)+dE
elif np.random.rand()<np.exp(-dE/T): S[pos_x,pos_y] = -S[pos_x,pos_y]
E_tot = Energy(N,S,pos_x,pos_y,J,B)-dE
m = np.sum(S)/(N**2)
e=E_tot/(N**2)
return m,S,e
def get_neighbor(N,pos_x,pos_y):
return [((pos_x+1)%N,pos_y),(pos_x-1,pos_y),(pos_x,(pos_y+1)%N),(pos_x,pos_y-1)]
def Energy(N,S,pos_x,pos_y,J,B):
E_tot=0.0 for pos_x in range(N):
for pos_y in range(N): E_tot=-J*S[pos_x,pos_y]*(np.sum(get_neighbor(N,pos_x,pos_y)))-B*S[pos_x,pos_y]+E_tot
return E_tot
def get_dE(N,S,pos_x,pos_y,J,B):
dE = 2*B*S[pos_x,pos_y]
for i,j in get_neighbor(N,pos_x,pos_y):
dE = dE +2*J*S[pos_x,pos_y]*S[i,j]
return dE
N=4
J=1
B=0
for T in np.arange(0.1, 10.0, 0.2):
m,S,e =Ising(T,N,J,B)
m_avg=(abs(m))
print("%f %f\n" %(T,e))