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):". 509
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. 490