How Can I Fix the IndexError in My 2D Ising Model Simulation?

[fluidistic]

Hi people,
I'd like to simulate a 2d array of atoms/ions under a magnetic field H via the Ising's model. I've read http://fraden.brandeis.edu/courses/phys39/simulations/AsherIsingModelReport.pdf to get an idea on how to proceed. I'm not sure I got the correct $\Delta E$ values though.
I'm having a problem, namely :

runfile('/home/user/Python-Stuff/ising_model/ising_model.py', wdir='/home/user/Python-Stuff/ising_model')
Traceback (most recent call last):

  File "<ipython-input-6-0c1af05397ef>", line 1, in <module>
    runfile('/home/user/Python-Stuff/ising_model/ising_model.py', wdir='/home/user/Python-Stuff/ising_model')

  File "/usr/lib/python3.5/site-packages/spyderlib/widgets/externalshell/sitecustomize.py", line 699, in runfile
    execfile(filename, namespace)

  File "/usr/lib/python3.5/site-packages/spyderlib/widgets/externalshell/sitecustomize.py", line 88, in execfile
    exec(compile(open(filename, 'rb').read(), filename, 'exec'), namespace)

  File "/home/user/Python-Stuff/ising_model/ising_model.py", line 78, in <module>
    if delta_energy(i,j) < 0:

  File "/home/user/Python-Stuff/ising_model/ising_model.py", line 41, in delta_energy
    delta_energy = -2 * ( J * (lattice[i+1][j] + lattice[i+1][j-1] + lattice[i][j+1] + lattice[i][j-1] ) + H)

IndexError: index 10 is out of bounds for axis 0 with size 10

Which I tried to fix by creating 2 functions (at lines 54 and 62 of my code), to no avail.

Here's the full code:

# 2 dimensional Ising's model. The goal is to obtain a hysteresis curve/loop for both ferromagnetic (J>0) and ferrimagnetic (J<0 and mu_i != mu_j) materials.
from random import choice, random, randint
import numpy as np
import math

J = 1
H = 1
beta = 1
# # of rows
n = 10
# # of columns
m = 10
'''# Generate the 2d array of atoms/ions as an nxm "matrix"
lattice = [[randint(0,1) for i in range(n)] for j in range(m)]

# Periodic boundary conditions: the first and last rows and columns must be equal
lattice.append(lattice[0])
for row in lattice: '''

# Try with numpy
lattice = np.empty(shape=(n,m))
np.shape(lattice)

for x in range(n):
  for y in range(m):
  lattice[x][y] = choice((-1,1))

# Set the last row equal to the 1st one (boundary condition)
for y in range(m): 
  lattice[n-1][y] = lattice[0][y]
# Set the last column equal to the 1st
for x in range(n):
  lattice[x][m-1] = lattice[x][0]
 
 
# Create a function that computes the delta energy resulting of flipping the spin state of lattice[i][j] 
def delta_energy(i, j):
  if lattice[i][j] ==1:
  delta_energy = 2 * ( J * (lattice[i+1][j] + lattice[i+1][j-1] + lattice[i][j+1] + lattice[i][j-1] ) + H) 
  else:
  delta_energy = -2 * ( J * (lattice[i+1][j] + lattice[i+1][j-1] + lattice[i][j+1] + lattice[i][j-1] ) + H) 
  return delta_energy

# Create a function that does the spin flipping 
def spin_flipper():
  if lattice[i][j] == 1:
  lattice[i][j] = -1
  else:
  lattice[i][j] = 1
  return lattice[i][j]
 
 
# Fix the problem of the boundary condition when we pick for instance i=0 and we have to calculate lattice[i-1][something], since we go out of range.
def boundary_conditions_fixer_row(i):
  if i-1 < 0:
  return n-1
  if i+1 > n:
  return 0
  else:
  return i

def boundary_conditions_fixer_column(j):
  if j-1 < 0:
  return m-1
  if j+1 > m:
  return 0
  else:
  return j
 
# Now the Monte Carlo method. We run many successive times the same algorithm in hopes to convering to the equilibrium state. 
for trial in range(15):
  # Pick a random atom/ion
  i = randint(0, n-1)
  j = randint(0, m-1) 
 

  # if delta_energy <0, flip the spin. Else, flip it with probability exp(-beta * delta_energy)
  if delta_energy(i,j) < 0:
  spin_flipper()
  else:
  if random() < exp(-beta * delta_energy):
  spin_flipper()

# Now let's see how the spins are
print(lattice)

So I don't know how to deal with the problem of picking a random element in the nxm matrix if that element is either on the 1st or last row or first or last column. Because in order to calculate $\Delta E$, I need to use lattice[i+1][j] and if i=n, I get out of range. I'd like this out of range to be fixed and if such thing happen, make the element lattice[i+1][j] be worth lattice[0][j].
Maybe I could add 2 new rows and 2 new columns to artificially fix this problem but that looks an ugly way to do it. Any comment is appreciated.

 

[Mark44]

If n is set to 10, range(n) results in 0, 1, 2, 3, ..., 9

If your array has 10 elements, their indexes run from 0 through 9, not 1 through 10. I believe this is what's causing your problem.

LaTeX tip: use two $ characters at start and end (standalone stuff) or two # characters at start and end (inline stuff).

For $\Delta E$, type ##\Delta E##.

 

[fluidistic]

Thanks for the latex tips, I edited my post.
I don't think that range(n) is the problem. The matrix generated is indeed "correct", filled randomly with -1's and 1's and the last row is the same as the first row, and the last column is the same as the first one.
I must face the problem of picking a random element and if that element is on an edge of the matrix, I'm not sure how to tell the program not to look outside the matrix to get to the nearest neighbors, but to get to the correct edge.

Edit: I fixed the errors by using the modulus operator %. The code now runs fine, although I don't know whether the output is correct!

 

[fluidistic]

I fixed everything (tons of errors and debugging).
Result:

. That's a part of a hysteresis loop, so I am happy. This means I have chances not to have goofed.Edit: Unfortunately my code is broken somehow. The output depends whether or not I include positive and negative values for H, which should not be the case.
Here's my latest version:

# 2 dimensional Ising's model. The goal is to obtain a hysteresis curve/loop for both ferromagnetic (J>0) and ferrimagnetic (J<0 and mu_i != mu_j) materials.
from random import choice, random, randint
import numpy as np
from math import exp
import pylab

J = 1
#H = -10
beta = 1
# # of rows
n = 15
# # of columns
m = 15
'''# Generate the 2d array of atoms/ions as an nxm "matrix"
lattice = [[randint(0,1) for i in range(n)] for j in range(m)]

# Periodic boundary conditions: the first and last rows and columns must be equal
lattice.append(lattice[0])
for row in lattice: '''

# Try with numpy
def set_lattice():
  lattice = np.empty(shape=(n,m))
  np.shape(lattice)

  for x in range(n):
  for y in range(m):
  lattice[x][y] = choice((-1,1))

  # Set the last row equal to the 1st one (boundary condition)
  for y in range(m):   
  lattice[n-1][y] = lattice[0][y]
  # Set the last column equal to the 1st
  for x in range(n):
  lattice[x][m-1] = lattice[x][0]
  return lattice

# Store a base lattice in memory to re-use at every trial of Monte Carlo method
base_lattice = set_lattice()
   
# Create a function that computes the delta energy resulting of flipping the spin state of lattice[i][j]   
def delta_energy(i, j , H):
  if lattice[i][j] ==1:
  delta_energy = 2 * ( J * (lattice[(i+1) % n][j % m] + lattice[(i+1) % n ][(j-1) % m] + lattice[i % n][(j+1) % m] + lattice[i % n][(j-1) % m] ) + H)   
  else:
  delta_energy = -2 * ( J * (lattice[(i+1) % n][j % m] + lattice[(i+1) % n][(j-1) % m] + lattice[i % n][(j+1) % m] + lattice[i % n][(j-1) % m] ) + H)   
  return delta_energy

# Create a function that does the spin flipping   
def spin_flipper(i,j):
  if lattice[i][j] == 1:
  lattice[i][j] = -1
  else:
  lattice[i][j] = 1
  return lattice[i][j]
   
   
'''# Fix the problem of the boundary condition when we pick for instance i=0 and we have to calculate lattice[i-1][something], since we go out of range.
def boundary_conditions_fixer_row(i):
  if i-1 < 0:
  return n-1
  if i+1 > n:
  return 0
  else:
  return i

def boundary_conditions_fixer_column(j):
  if j-1 < 0:
  return m-1
  if j+1 > m:
  return 0
  else:
  return j
'''

def run_simulation(num_trials, H):   
  # Now the Monte Carlo method. We run many successive times the same algorithm in hopes to convering to the equilibrium state.   
  for trial in range(num_trials):
  # reset the lattice after each simulation!
  base_lattice
  # Pick a random atom/ion
  i = randint(0, n-1)
  j = randint(0, m-1)   
   

  # if delta_energy <0, flip the spin. Else, flip it with probability exp(-beta * delta_energy)
  if delta_energy(i,j, H) < 0:
  spin_flipper(i,j)
  else:
  if random() < exp(-beta * delta_energy(i, j, H)):
  spin_flipper(i,j)
  magnetization = np.sum(lattice)
  return magnetization
# Now let's see how the spins are
#print(lattice)   def showPlot(num_trials, Hfields, title = 'M vs H', x_label = 'H field', y_label = 'Magnetization'):
#  Hfields = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
  y_values = []
  for H in Hfields:
  print ("Plotting", H, "H fields...")
  y_values.append(run_simulation(num_trials, H))
  print(y_values)
  pylab.plot(Hfields, y_values)
  pylab.title(title)
  pylab.legend(('Magnetization'))
  pylab.xlabel(x_label)
  pylab.ylabel(y_label)
  pylab.show()  
  # call like so: showPlot(25000, [3,-0.5, -0.56, -0.57, -0.6,-0.63, -0.67, -0.7])

 

[JorisL]

As an aside, you can use helical boundary conditions which are easier to program and generally faster.
The first time I took a computational physics class I programmed very naive $\to$ runs up to 24 hours (on a cluster nonetheless) for a reasonably small lattice.
The second time we got some tips $\to$ I think the runs took no longer than 30 minutes.

If you like I can take a look if I can find the code and pass it on, the results were correct. The method of obtaining them not so much (best coding practices as a metric).

 

[Fightfish]

Edit: Unfortunately my code is broken somehow. The output depends whether or not I include positive and negative values for H, which should not be the case.

Mind elaborating?
As another aside, I would highly recommend the book Monte Carlo Methods in Statistical Physics by Newman and Barkema. It gives a really detailed discussion of the various subtypes of Monte Carlo methods (beyond just the Metropolis algorithm) as applied to Ising-type models and also introduces many tips and tricks.

 

[fluidistic]

Mind elaborating?
As another aside, I would highly recommend the book Monte Carlo Methods in Statistical Physics by Newman and Barkema. It gives a really detailed discussion of the various subtypes of Monte Carlo methods (beyond just the Metropolis algorithm) as applied to Ising-type models and also introduces many tips and tricks.

Nevermind, the output looked always different and seemd to be linked on what values for H I picked. For example if I picked a lot of different negative values for H and only 1 positive, I'd get an answer that's totally different than in the case of picking all negative values for H. (As if the single value of H positive could influence the magnetization for when H was negative). Apparently I get this behavior because I had picked too few steps (15-25k). With a higher number of steps this behavior is less pronounced (albeit it still is for 150k steps, that is, when I pick a 15x15 lattice).
And ok, thank you for the suggestion of the book.