Simulating the Ising Model in Python

Question

I taught myself the Metropolis Algorithm and decided to try code it in Python. I chose to simulate the Ising model. I have an amateur understanding of Python and with that here is what I came up with -

import numpy as np, matplotlib.pyplot as plt, matplotlib.animation as animation

def Ising_H(x,y):

    s = L[x,y] * (L[(x+1) % l,y] + L[x, (y+1) % l] + L[(x-1) % l, y] + L[x,(y-1) % l])
    H = -J * s
    return H

def mcstep(*args): #One Monte-Carlo Step - Metropolis Algorithm

    x = np.random.randint(l)
    y = np.random.randint(l)
    i = Ising_H(x,y)
    L[x,y] *= -1
    f = Ising_H(x,y)
    deltaH = f - i
    if(np.random.uniform(0,1) > np.exp(-deltaH/T)):
        L[x,y] *= -1

    mesh.set_array(L.ravel())
    return mesh,

def init_spin_config(opt):

    if opt == 'h':
        #Hot Start
        L = np.random.randint(2, size=(l, l)) #lxl Lattice with random spin configuration
        L[L==0] = -1
        return L

    elif opt =='c':
        #Cold Start
        L = np.full((l, l), 1, dtype=int) #lxl Lattice with all +1
        return L

if __name__=="__main__":

     l = 15 #Lattice dimension
     J = 0.3 #Interaction strength
     T = 2.0 #Temperature
     N = 1000 #Number of iterations of MC step
     opt = 'h' 

     L = init_spin_config(opt) #Initial spin configuration

     #Simulation Vizualization
     fig = plt.figure(figsize=(10, 10), dpi=80)
     fig.suptitle("T = %0.1f" % T, fontsize=50)
     X, Y = np.meshgrid(range(l), range(l))
     mesh = plt.pcolormesh(X, Y, L, cmap = plt.cm.RdBu)
     a = animation.FuncAnimation(fig, mcstep, frames = N, interval = 5, blit = True)
     plt.show()

Apart from a 'KeyError' from a Tkinter exception and white bands when I try a 16x16 or anything above that, it looks and works fine. Now what I want to know is if this is right because -

I am uncomfortable with how I have used FuncAnimation to do the Monte Carlo simulation AND animate my mesh plot - does that even make sense?

And How about that cold start? All I am getting is a red screen.

Also, please tell me about the KeyError and the white banding.

The 'KeyError' came up as -

Exception in Tkinter callback
Traceback (most recent call last):
   File "/usr/lib/python2.7/lib-tk/Tkinter.py", line 1540, in __call__
      return self.func(*args)
   File "/usr/lib/python2.7/lib-tk/Tkinter.py", line 590, in callit
   func(*args)
   File "/usr/local/lib/python2.7/dist-packages/matplotlib/backends/backend_tkagg.py", line 147, in _on_timer
      TimerBase._on_timer(self)
   File "/usr/local/lib/python2.7/dist-packages/matplotlib/backend_bases.py", line 1305, in _on_timer
      ret = func(*args, **kwargs)
   File "/usr/local/lib/python2.7/dist-packages/matplotlib/animation.py", line 1049, in _step
      still_going = Animation._step(self, *args)
   File "/usr/local/lib/python2.7/dist-packages/matplotlib/animation.py", line 855, in _step
      self._draw_next_frame(framedata, self._blit)
   File "/usr/local/lib/python2.7/dist-packages/matplotlib/animation.py", line 873, in _draw_next_frame
      self._pre_draw(framedata, blit)
   File "/usr/local/lib/python2.7/dist-packages/matplotlib/animation.py", line 886, in _pre_draw
      self._blit_clear(self._drawn_artists, self._blit_cache)
   File "/usr/local/lib/python2.7/dist-packages/matplotlib/animation.py", line 926, in _blit_clear
      a.figure.canvas.restore_region(bg_cache[a])
KeyError: <matplotlib.axes._subplots.AxesSubplot object at 0x7fd468b2f2d0>

Answer 1

You are asking a lot of questions at a time.

• KeyError: cannot be reproduced. It's strange that it should only occur for some array sizes and not others. Possibly something is wrong with the backend, you may try to use a different one by placing those lines at the top of the script
  import matplotlib
matplotlib.use("Qt4Agg")
• white bands: cannot be reproduced either, but possibly they come from an automated axes scaling. To avoid that, you can set the axes limits manually plt.xlim(0,l-1) plt.ylim(0,l-1)

• Using FuncAnimation to do the Monte Carlo simulation is perfectly fine. f course it's not the fastest method, but if you want to follow your simulation on the screen, there is nothing wrong with it. One may however ask the question why there would be only one spin flipping per time unit. But that is more a question on the physics than about programming.

• Red screen for cold start: In the case of the cold start, you initialize your grid with only 1s. That means the minimum and maximum value in the grid is 1. Therefore the colormap of the pcolormesh is normalized to the range [1,1] and is all red. In general you want the colormap to span [-1,1], which can be done using vmin and vmax arguments.

  mesh = plt.pcolormesh(X, Y, L, cmap = plt.cm.RdBu, vmin=-1, vmax=1)
  This should give you the expected behaviour also for the "cold start".