计算科学 - Ising NiO模型能量 - 吾爱随笔录

我正在为 NiO 模拟 Ising 模型。我已经模拟了 2d、3d、三角形晶格，并尝试对 NiO 模型做同样的事情。

有论文说基态能量约为每个粒子 -36 mev。每个粒子我得到-21 mev（我应该不考虑氧原子吗？，那么它将是-42mev）。我做错了什么吗。以下是代码。

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:
                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/=2
    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:
            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

此处显示 J1 和 J2

此外，如果您可以评论磁化强度、容量热、磁化率，那也会很有帮助