- #1

fluidistic

Gold Member

- 3,655

- 101

## Main Question or Discussion Point

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 :

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:

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.

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 :

Code:

```
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
```

Here's the full code:

Python:

```
# 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)
```

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.

Last edited: