Ising NiO model energy

Question

I am simulating Ising model for NiO. I have simulated for 2d,3d,triangular lattices, and have tried to do the same with NiO model.

There are papers which say that the ground state energy is around -36 mev per particle. Where as I am getting -21 mev per particle (should I not consider Oxygen atoms?, then it will be -42mev). Is there something wrong I am doing. Following is the code.

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 and J2 are shown here

Also, if you can comment on magnetisation, capacity heat, magnetic susceptibility, that will be helpful too

Answer 1

As, $J2 = -21 meV$ is more dominant than $J1 = 2.3meV$, the system is antiferromagnetic in nature. The expected ground state energy per moelcule(NiO) is $-42meV$. When the state reaches equilibirum, two of the nearest neighbors are alike and two are opposite, cancelling the energy contirbutions of each other. The energy contribution is from the second nearest negihbours giving $-4*21 = -84 meV$, but as each bond(energy) is counted twice, it should be halved. The above calculations are for Ni atom which has magnetic moment. As oxygen does not couple with other spins, the average energy per site(atom) is $\frac{-42}{2} = -21 meV$.

I have run through the code. It is correct.