Hi people,(adsbygoogle = window.adsbygoogle || []).push({});

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.Code (Text):

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:

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].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.

**Physics Forums - The Fusion of Science and Community**

Dismiss Notice

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Python 2d Ising's model simulation

Have something to add?

Draft saved
Draft deleted

Loading...

Similar Threads - Ising's model simulation | Date |
---|---|

3-D simulation software with mathematical parameters | Feb 22, 2018 |

Python Topic Modeling takes too long | Nov 24, 2017 |

Solution of Bose Hubbard model using ARPACK | Nov 6, 2017 |

Lost in Hidden Markov Model | Jul 9, 2017 |

Unified Modeling Language | Nov 23, 2016 |

**Physics Forums - The Fusion of Science and Community**