Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Python 2d Ising's model simulation

  1. Apr 29, 2016 #1

    fluidistic

    User Avatar
    Gold Member

    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 :
    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
     
    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:
    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)

     
    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.
     
    Last edited: Apr 29, 2016
  2. jcsd
  3. Apr 29, 2016 #2

    Mark44

    Staff: Mentor

    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##.
     
  4. Apr 29, 2016 #3

    fluidistic

    User Avatar
    Gold Member

    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!
     
    Last edited: Apr 29, 2016
  5. Apr 29, 2016 #4

    fluidistic

    User Avatar
    Gold Member

    I fixed everything (tons of errors and debugging).
    Result: kqKFma0.png . 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:
    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
    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])
     
     
    Last edited: Apr 29, 2016
  6. Apr 30, 2016 #5
    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).
     
  7. May 1, 2016 #6
    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.
     
  8. May 1, 2016 #7

    fluidistic

    User Avatar
    Gold Member

    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: 2d Ising's model simulation
  1. Quantum Simulation (Replies: 5)

  2. Physics simulation (Replies: 5)

  3. Simulating the world (Replies: 14)

  4. Modeling and simulaton (Replies: 5)

Loading...