Using canvas to simulate the Ising model

Question

I had originally posted this in ComputationalScience but I think perhaps it is best here as I think the challenge I have is with using classes and creating figures rather than the actual mathematical / physics behind the simulation.

I am attempting to simulate a 2D ising model. The simulation consists of an nxn matrix, whereby each element of the matrix has a value of 1 or -1, which flips based on condidtions outlined in the code. I would like to produce a moving figure which essentially produces a different picture for each iteration. but my figure will not appear. I am using an identical method (fig.canvas.draw()) to a code I ran previously with no problems.

Here is my code:

import numpy as np
import matplotlib.pyplot as plt

#initialised values 

T = 1.0
N = 5
sweeps = 50

nx = N
ny = N #size of the NxN matrix

#Ising code

class Ising:
    """creating the ising model in this class"""
#creating a function which creates an initial matrix with randomly orientated spins of 1 and -1    
    def __init__(self, nx, ny, random=False):
        if random:
            self._lattice = np.random.choice(np.array([-1,1]),size=(N,N))
        else:
            self._lattice = np.ones(N,N)


#creating a function which flips the spin of the individual lattice points        
    def flipspin(self, x, y):
        if self._lattice[x,y] == 1.0:
            self._lattice[x,y] == -1.0
        else:
            self._lattice[x,y] == 1.0


    def nearestneighbours(self, N, x, y):
         neighbourList = [[x+1,y],[x-1,y],[x,y+1],[x,y-1]]
        #check if neighbors have walked off the edge of the box, apply periodic boundary conditions if so
         for i,element in enumerate(neighbourList):
            for j,entry in enumerate(element):
                if (entry < 0):
                    neighbourList[i][j] = self.shape()[j]-1
                elif (entry >= self.shape()[j]):
                    neighbourList[i][j] = 0
         return neighbourList



def Hamiltonianenergy(lattice, x ,y):
    return 2 * lattice[x, y] * sum(lattice.nearestneighbours(lattice, N, x, y))


def metropolis(lattice,temperature,eps=1):
    J = 1    
    t = 0
    switch = 0
    energies = 0
    energycount = []
    res = []
    fig = plt.figure()    
    #evaluate inverse temperature
    beta = 1/temperature
    #construct ising lattice
    lattice = Ising(nx,ny,random=True)
    #sample x,y point on lattice from uniform distribution over the integers
    #flip spin and calculate dH

    #accept or reject according to metropolis function
    for t in range(sweeps):    
        i = np.random.randint(N) #randomly selecting an individual lattice point
        j = np.random.randint(N) 
        dH = Hamiltonianenergy(lattice,(i, j),eps,flip=True)
        if (dH > 0):
            boltzmann = np.exp(-1*beta*dH)
            randNum = np.random.uniform()
            result = True
            if boltzmann <= randNum:
                lattice.flip(i, j)
                result = False
        else:
            boltzmann = 1
            randNum = 1
            result = True
            energycount.append(energies)

        plt.imshow(lattice)
        fig.canvas.draw()
        plt.pause(0.01)    
        plt.show()     
        t +=1
        print t
        print energycount
        res.append(switch)

I am receiving the warning message on my figure:

local variable image is assigned to but never used

Answer 1

I've managed to fix all the programming issue of your code, but as I stated in a comment, I believe there is a problem with your metropolis part which is nod doing any flipping. Particularly, I think you're not doing your dH properly.

Placing comments (hopefully all of them) where I did change things to be able to run the code.

import numpy as np
import matplotlib.pyplot as plt

#initialised values 

T = 1.0
N = 5
sweeps = 5

nx = N
ny = N #size of the NxN matrix

#Ising code

class Ising:
    """creating the ising model in this class"""
    #creating a function which creates an initial matrix with randomly orientated spins of 1 and -1    
    def __init__(self, nx, ny, random=False):
        self._lattice = np.random.choice(np.array([-1,1]),size=(N,N))

#creating a function which flips the spin of the individual lattice points        
    def flipspin(self, x, y):
        self._lattice[x,y] *= -1  # simplified the flipping

    def __getitem__(self, key):  # added a special method to be able to do the lattice[(i,j)] lookup
        x, y = key
        return self._lattice[x, y]

    def __str__(self):  # just for debugging purposes
        return str(self._lattice)

    def nearestneighbours(self, N, x, y):  # this was fine, you're doign the sum here already
       return self._lattice[(x - 1) % N, y] + self._lattice[(x + 1) % N, y] + \
           self._lattice[x, (y - 1) % N] + self._lattice[x, (y + 1) % N]

    def show(self):  # using a method to show graphically your 2D matrix
        fig = plt.figure()    
        image = plt.imshow(self._lattice)
        fig.canvas.draw()
        plt.show()    

def Hamiltonianenergy(lattice, x ,y):  # cleaned up to be able to call it without problems
    return 2 * lattice[(x, y)] *  lattice.nearestneighbours(N, x, y)

def metropolis(temperature,eps=1):
    J = 1    
    t = 0
    switch = 0
    energies = 0
    energycount = []
    res = []

    #evaluate inverse temperature
    beta = 1/temperature
    #construct ising lattice
    lattice = Ising(nx,ny,random=True)
    print(lattice)

    #accept or reject according to metropolis function
    for t in range(sweeps):    
        i = np.random.randint(N) #randomly selecting an individual lattice point
        j = np.random.randint(N) 
        dH = Hamiltonianenergy(lattice, i,  j)  # here should be the problem, you have to make the difference of energies before and after the flip!

        if (dH > 0):
            boltzmann = np.exp(-1*beta*dH)
            randNum = np.random.uniform()
            result = True
            if boltzmann <= randNum:
                lattice.flipspin(i, j)  # you where not calling the method properly
                result = False
        else:
            boltzmann = 1
            randNum = 1
            result = True
            energycount.append(energies)

        # some console output
        print(t, i, j, result)
        print("dh: ", dH)
        print(lattice)

        t +=1

        res.append(switch)
        lattice.show()

metropolis(1000)  # running the metropolis function when executing the file