I am trying to compute the Ising model for NiO. As O carries no magnetic moment, I only need to consider the case of Ni which requires a second nearest neighbour Ising model. As can be seen in the figure below, the Ni atoms interact with the their nearest neighbours with a coupling constant J1 = 2.3 meV and their second nearest neighbours with a coupling constant J2 = -21 meV.
I have created some code that generates a matrix that alternates 1 and -1 (spin up/spin down) in every second entry and 0 for every other entry (representing oxygen). I have also defined functions that will flip the spin for every nearest neighbour and second nearest neighbour. As the dominating coupling constant J2 < 0, the system should be antiferromagnetic so the spins should align diagonally repeating the pattern (1,0,-1,0) e.g:
[ 0 -1 0 1 0 -1 0 1 0 -1]
[ 1 0 -1 0 1 0 -1 0 1 0]
[ 0 1 0 -1 0 1 0 -1 0 1]
[-1 0 1 0 -1 0 1 0 -1 0]
[ 0 -1 0 1 0 -1 0 1 0 -1]
[ 1 0 -1 0 1 0 -1 0 1 0]
[ 0 1 0 -1 0 1 0 -1 0 1]
[-1 0 1 0 -1 0 1 0 -1 0]
[ 0 -1 0 1 0 -1 0 1 0 -1]
[ 1 0 -1 0 1 0 -1 0 1 0]
However when I run the code, I am not able to achieve that. I can reach a certain amount of order at low temperatures (T~2) but not total ferromagnetism as can be seen below. Going lower (e.g. T~0.01) will yield disorder:
Here is my code:
#!/usr/bin/env python
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import random
#constants
N = 10 #dimensions of matrix
J1 = 2.3 #coupling constant
J2 = -21
h = 0 #magnetic field, must be set to 0 to compute observables
counts = 100
T = 2 #temperature
k=1 #boltzmann constant
class initial_lattice:
def __init__(self,N): #create initial matrix of size NxN
self.N=N
self.matrix_lattice()
def matrix_lattice(self):
self.lattice = np.random.choice([-1, 1], (N, N))
self.lattice[::2, ::2] = 0
self.lattice[1::2, 1::2] = 0
lattice1=initial_lattice(N)
#function that sums up all neighbouring sites of the inital position. %N imposes a boundary condition so the function knows when to stop
def diagonal_neighbours(matrix,x,y,N):
d1 = matrix[(x+1)%N, (y+1) %N]
d2 = matrix[(x+1)%N, (y-1) %N]
d3 = matrix[(x-1)%N, (y+1)%N]
d4 = matrix[(x-1) %N, (y-1)%N]
return d1 + d2 + d3 + d4
def lateral_neighbours(matrix,x,y,N):
l1 = matrix[x, (y+2) %N]
l2 = matrix[x, (y-2) %N]
l3 = matrix[(x+2) %N, y]
l4 = matrix[(x-2) %N, y]
return l1 + l2 + l3 + l4
#function for change in energy
def deltaE(matrix, x, y, N, J1, J2, h):
return -(2*J1*matrix[x,y]*(diagonal_neighbours(matrix,x,y,N)))-(2*J2*matrix[x,y]*(lateral_neighbours(matrix,x,y,N)))+2*h*matrix[x,y]
#metropolis algorithim
def metropolis(matrix, counts,N, T, J1,J2, h, k):
for n in range (counts):
for y in range(0, N):
for x in range(0,N):
if deltaE(matrix, x, y, N, J1, J2, h)>=0:
matrix[x,y] *= -1 #if energy change is greater than/equal to 0, flips spin
else:
r = random.random() #generates random number
if r<np.exp(deltaE(matrix, x, y, N, J1, J2, h)/(k*T)):
matrix[x,y] *= -1 #if random number generated between 0 and 1 is less than exp^dE/k*T flips spin
return matrix
print metropolis(lattice1.lattice, counts,N, T, J1, J2, h, k)
plt.imshow(metropolis(lattice1.lattice, counts,N, T, J1, J2, h, k),cmap='bwr',interpolation="none")
plt.show() #plots Ising model in equilibrium
Any help would be greatly appreciated.