I am trying to use the metropolis algorithm to simulate the Ising model. The problem that I am having is that the code will not settle all the time. For high beta values, the energy preference should be that the spins are all in the same direction, however in my code, even though it works most of the time, occasionally there will be horizontal or vertical bands of spin up and then spin down in the grid (grid2). This then makes the average spin count settle at a value which is not 1 or -1 which it should be for high beta. I have some idea where the problem lies - my delta energy at the boundary when this occurs is 4 and so at high beta this makes it really unlikely that the spin will flip and so it ends up being stuck. I was thinking that maybe I am not running it long enough but looking at the average_spin_count on the variable explorer it seems to settle after maximum 200 turns. I was also thinking maybe using diagonals but other codes and simulations of this don't seem to use this idea and so I was hoping someone could point out where I am going wrong. (You have to run the code multiple times until you get a case where there are not all spin up/down). Thanks
import numpy as np
import numpy.random as rnd
N = 20
#grid size
random_number = np.random.randint(0,1000)
prob = 0.5
#probability of starting with spin down.
np.random.seed(random_number)
grid =np.random.choice(list(range(-1, 0)) + list(range(1, 2)), size=(N, N), p=[prob,1-prob])
def find_energy(grid, r, c):
if r < N-1 and c < N-1:
energy = -((grid[r,c]*grid[r+1,c])+(grid[r,c]*grid[r-1,c])+
(grid[r,c]*grid[r,c+1])+(grid[r,c]*grid[r,c-1]))
if r == (N-1) and c < N-1:
energy = -((grid[r,c]*grid[0,c])+(grid[r,c]*grid[r-1,c])+
(grid[r,c]*grid[r,c+1])+(grid[r,c]*grid[r,c-1]))
if c == (N-1) and r < N-1:
energy = -((grid[r,c]*grid[r+1,c])+(grid[r,c]*grid[r-1,c])+
(grid[r,c]*grid[r,0])+(grid[r,c]*grid[r,c-1]))
if c== (N-1) and r==(N-1):
energy = -((grid[r,c]*grid[0,c])+(grid[r,c]*grid[r-1,c])+
(grid[r,c]*grid[r,0])+(grid[r,c]*grid[r,c-1]))
return energy
def get_coordinate():
"""Returns a random coordinate as a tuple"""
return tuple(np.random.randint(0, high=N, size=(1, 2))[0])
def probability( delta_energy, beta):
probability = np.exp(-1*beta*delta_energy)
return probability
def compare_energies(n, beta):
total_number_spin_up=np.array([])
total_number_spin_down=np.array([])
for i in range(n):
original_point = get_coordinate()
original_energy = find_energy(grid, original_point[0], original_point[1])
new_energy = original_energy * -1
delta_energy = new_energy-original_energy
x = (original_point[0])
y = (original_point[1])
if delta_energy <= 0:
grid[x,y] = grid[x, y]*-1
#flips the spin
elif delta_energy > 0:
if probability(delta_energy, beta) > rnd.uniform(0,1):
grid[x,y] = (grid[x, y])*-1
number_of_spin_up = np.count_nonzero(grid == 1)
total_number_spin_up = np.append(total_number_spin_up,number_of_spin_up)
total_number_spin_down = np.append(total_number_spin_down, pow(N,2)-number_of_spin_up)
spin_up = total_number_spin_up[-1]
spin_down = total_number_spin_down[-1]
average_spin = (spin_up-spin_down)/N**2
return grid, average_spin
beta = 10
#runs the code to settle it.
average_spin_count = np.array([])
for i in range(300):
grid_2, average_spin = compare_energies(N**2, beta)
average_spin_count = np.append(average_spin_count, average_spin)
y = average_spin_count[-1]
If you plot the final grid of such a case, you will see what happens, use matplotlib for example, see below:
import matplotlib.pyplot as plt
plt.imshow(grid)
plt.show()
You will find that the space is separated in two regions. Alternatively, you can initialize your grid as follows instead of random numbers.
grid = np.zeros((N, N))
grid[:N//2] = 1
You will see a similar behavior. The solution is very simple. If you calculate the energy and the corresponding probability for a spin flip at the boundary, you will find an energy value of -3 leading to a probability of 9.76e-27 which is basically 0 for your given beta. Therefore, this state is very stable and you will likely not see a change within 300 iterations.