Computing the Ising Model for NiO

Question

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.

Answer 1

Your Metropolis algorithm's loop is wrong as this is not a Metropolis procedure.

You're looping through all the spins, flipping them is the flip reduces the energy when you actually have to pick $x,y \in [0,N]$ randomly at each step and check if flipping is reducing the energy. If it does, you keep the flip, if it doesn't you keep it with the probability distribution you're using (i.e. keep if $\chi > 1/k_BT$ with $\chi$ uniformly distributed on $[0,1]$).

I would write something like:

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

Note that this loop only return one matrix. If you want to compute any statistically interesting value, you have to compute regularly values like energy or mean orientation.

Edit : I tried to make your code work and beside this algorithm problem, a few other comments came to my mind. Here is my code with a few comments:

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

Simple note, I didn't need 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)

Global variable are rarely a good idea (in short code like that it's no big deal). Use conditions to use/define them like "if (main):". More importantly, when using physical units (here Ev and K), convert everything in the same set of units. $k_B$ has a dimension of energy per kelvin, use it. If you want to use dimensionless units, you have to compute $J1/k_BT$ and $J2/k_BT$ first.

I've made a function to compute the energy of a single configuration:

#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

It is important to note the evolution of your instantaneous observables for 2 reasons: 1. If your outputs are wrong, you need to check how the different variables evolve during the simulation 2. If you want to compute macroscopic observables, you need to make some statistics on instantaneous observables. Therefore your initial value of the energy is important to check the convergence.

This is my version of your Metropolis algorithm:

#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

I modified your algorithm in several ways. Your version was looping in the atoms in a deterministic way as I commented first. This one pick the atoms in a pseudo-random way which is the heart of the Metropolis algorithm. Also, you can use intermediate variables for condition in order to make the code easier to read. Writing lines like:

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

Quickly make things confusing. Python is not really designed for direct functional evaluation like this. (Also note that you got the signs wrong in your initial version.) You also used random.random() library which doesn't explicitly give you uniform distribution. random.uniform() is here for you.

In the end:

#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

Don't hesitate to stock intermediate variable in this kind of situation. If you pass a function like metropolis(...) as a variable of another function (say print) it will cause a new evaluation of the function. Again, python is not well suited for this kind of functional programming style (unless you're using generators with lazy evaluation).

A final comment on your expected behavior: MC-metropolis algorithm is random. You expected to observe the configuration at equilibrium but after 100 random flips (not all of them accepted), it is not said you'll converge to equilibrium. the shortest convergence I observed with my algorithm was around 250 flips but don't hesitate to set them to 10000. Such a simple code must be quite fast. Be aware though that the convergence will require more and more cycles as the system grows.