您的 Metropolis 算法的循环是错误的，因为这不是 Metropolis 程序。

您正在循环所有旋转，翻转它们是翻转会减少您实际必须选择的能量$x,y \in [0,N]$在每一步随机检查翻转是否会减少能量。如果是，你保持翻转，如果不是，你保持它与你正在使用的概率分布（即保持如果$\chi > 1/k_BT$和$\chi$均匀分布在$[0,1]$）。

我会写这样的东西：

import math
[...]
for n in range(counts):
    x = random.randint(0,N)
    y = random.randint(0,N)
    if deltaE(matrix, x, y, N, J1, J2, h)>=0:
    [...]
return matrix

请注意，此循环仅返回一个矩阵。如果要计算任何统计上有趣的值，则必须定期计算能量或平均方向等值。

编辑： 我试图让你的代码工作，除了这个算法问题，我想到了一些其他的评论。这是我的代码和一些评论：

#!/usr/bin/env python
import numpy as np
#import scipy as sp
import matplotlib.pyplot as plt
import random

简单的说明，我不需要 scipy。

#constants
N    = 10 #dimensions of matrix
J1   = 2.3*10**(-3) #coupling constant in Ev
J2   = -21*10**(-3)
h    = 0 #magnetic field, must be set to 0 to compute observables 
counts = 1000
T    = 1 #temperature in K
k    = 8.6173303*10**(-5) #boltzmann constant in Ev/K
beta = 1/(k*T)

全局变量很少是一个好主意（在这样的短代码中没什么大不了的）。使用条件来使用/定义它们，例如“if (main):”。更重要的是，当使用物理单位（此处为 Ev 和 K）时，将所有内容转换为同一组单位。$k_B$每开尔文有一个能量维度，使用它。如果要使用无量纲单位，则必须计算$J1/k_BT$和$J2/k_BT$第一的。

我制作了一个函数来计算单个配置的能量：

#Compute actual energy don't forget periodic boundary conditions
def Energy(matrix, J1, J2, h):
    E = 0
    for x in range(N):
        for y in range(N):
            N1 = (x+1)%N
            N2 = (y+1)%N
            N3 = (x-1)%N
            N4 = (y-1)%N
            for i in [N1,N3]:
                for j in [N2,N4]:
                    E += matrix[x,y]*J1*matrix[i,j]/2.
            N1 = (x+2)%N
            N2 = (y+2)%N
            N3 = (x-2)%N
            N4 = (y-2)%N
            for i in [N1,N3]:
                for j in [N2,N4]:
                    E += matrix[x,y]*J2*matrix[i,j]/2.
            E -= h*matrix[x,y]
    return E

重要的是要注意瞬时可观测量的演变，原因有两个： 1. 如果输出错误，则需要检查不同变量在模拟过程中如何演变 2. 如果要计算宏观可观测量，则需要做一些瞬时观测值的统计。因此，您的能量初始值对于检查收敛性很重要。

这是我的 Metropolis 算法版本：

#metropolis algorithm 
def metropolis(matrix, counts, N, beta, J1, J2, h):
# This array E will stock energy values, E_conf will be the instantaneous value.
    E = np.zeros(counts)
    E_conf = Energy(matrix, J1, J2, h)
    for n in range (counts):
# Coordinates of the flip have to be reset for each move and randomly chosen!
        x = 0
        y = 0
# I know that (0,0) gives an Oxygen coordinates, I randomly chose values until I get a Ni atom
        while (matrix[x,y] == 0):
            x = random.randint(0,N-1)
            y = random.randint(0,N-1)
        dE = deltaE(matrix, x, y, N, J1, J2, h)
        if (dE <= 0):
            matrix[x, y] *= -1
            E_conf += dE
        else:
            r = random.uniform(0,1)
            tau = np.exp(-dE*beta)
            if (r < tau) :
                matrix[x, y] *= -1
                E_conf += dE
        E[n] = E_conf
    plt.plot(E)
    plt.show()
    return matrix

我以多种方式修改了您的算法。正如我首先评论的那样，您的版本以一种确定的方式在原子中循环。这个以伪随机方式选择原子，这是 Metropolis 算法的核心。此外，您可以使用中间变量作为条件，以使代码更易于阅读。像这样写行：

if (function(x,y,z,a,b,c) > 0)

快速让事情变得混乱。Python 并不是真正为像这样的直接功能评估而设计的。（另请注意，您在初始版本中弄错了符号。）您还使用了 random.random() 库，它没有明确地为您提供统一分布。random.uniform() 在这里为您服务。

到底：

#Initialising
lattice1=initial_lattice(N)
mat = lattice1.lattice

print ("Starting from initial configuration:")
print (mat)

plt.imshow(mat,cmap='bwr',interpolation="none") 
plt.show() #plots initial configuration

mat = metropolis(mat, counts, N, beta, J1, J2, h)

print (mat)

plt.imshow(mat,cmap='bwr',interpolation="none") 
plt.show() #plots last configuration

在这种情况下，不要犹豫存储中间变量。如果您将诸如 metropolis(...) 之类的函数作为另一个函数的变量（例如 print）传递，它将导致对该函数进行新的评估。同样，python 不太适合这种函数式编程风格（除非您使用具有惰性求值的生成器）。

对您预期行为的最后评论：MC-metropolis 算法是随机的。您希望在平衡时观察配置，但在 100 次随机翻转（并非所有翻转都被接受）之后，并不是说您会收敛到平衡。我用我的算法观察到的最短收敛时间大约是 250 次翻转，但毫不犹豫地将它们设置为 10000。这样一个简单的代码必须相当快。请注意，随着系统的增长，收敛将需要越来越多的周期。

def initial_state_Nio(N): 
    state = np.random.choice([-1, 1], (N, N))
    state[::2, ::2] = 0
    state[1::2, 1::2] = 0
    return state
def diag_nbrs(i,j,N): 
    return [((i+1)%N,(j+1)%N),((i+1)%N,(j-1)%N),((i-1)%N,(j+1)%N),((i-1)%N,(j-1)%N)]

def lat_nbrs(i,j,N): 
    return [(i,(j+2)%N),(i,(j-2)%N),((i+2)%N, j),((i-2)%N, j)]
def Energy_Nio(state, J1, J2, H):
    J1   = 2.3*10**(-3) #diagonal
    J2   = -21*10**(-3) #lateral coupling
    E = 0
    N = state.shape[0]
    for x in range(N):
        for y in range(N):
            if (x-y)%2: #avoids oxygen atoms to save time
                nbrs = diag_nbrs(x,y,N)
                for nbr in nbrs:
                    E += -state[x,y]*J1*state[nbr[0],nbr[1]]
                nbrs = lat_nbrs(x,y,N)
                for nbr in nbrs:
                    E += -state[x,y]*J2*state[nbr[0],nbr[1]]
    E -= H*state.sum()
    return E
def calcMag(state):
    return np.sum(state)
def step_update_Nio(state, beta, J1, J2,H,energy,mag,N):
    J1   = 2.3*10**(-3) #diagonal
    J2   = -21*10**(-3) #lateral coupling
    for i in range(N**2): #1 step per state on average
        dE = 0
        x = random.randint(0,N-1)
        y = random.randint(0,N-1)
        if (x-y)%2: #avoids oxygen atoms to save time
            nbrs = diag_nbrs(x,y,N)
            for nbr in nbrs:
                dE += 2*state[x,y]*J1*state[nbr[0],nbr[1]]
            nbrs = lat_nbrs(x,y,N)
            for nbr in nbrs:
                dE += 2*state[x,y]*J2*state[nbr[0],nbr[1]]
            dE += 2*H*state[x,y]
            if (dE <= 0):
                if state[x,y] == 1:
                    mag-=2
                else:
                    mag+=2
                energy += dE
                state[x, y] *= -1
            else:
                r = random.uniform(0,1)
                tau = np.exp(-dE*beta)
                if (r < tau) :
                    if state[x,y] == 1:
                        mag-=2
                    else:
                        mag+=2
                    energy += dE
                    state[x, y] *= -1
    return state,energy,mag
def run_Nio(state, steps, N, beta, J1, J2,H):
    J1   = 2.3*10**(-3) 
    J2   = -21*10**(-3)
    E = np.zeros(steps)
    M = np.zeros(steps)
    energy = Energy_Nio(state, J1, J2,H)
    mag = calcMag(state)
    for i in range(steps):
        state,energy,mag = step_update_Nio(state, beta, J1, J2,H,energy,mag,N)
        E[i] = energy
        M[i]= mag
    plt.plot(E)
    plt.show()
    plt.plot(M)
    plt.show()
    return state,E,M

我得到的结果是反铁磁顺序。