Rep:Mod:DMS3053-2

CMP Programming Experiment - Daniel Spencer, 00736964 Python IsingLattice.py script - Introduction to Monte Carlo simulation

   import numpy as np
   from math import exp

   class IsingLattice:

       n_cycles = 0

       def __init__(self, n_rows, n_cols):
           self.n_rows = n_rows
           self.n_cols = n_cols
           self.lattice = np.random.choice([-1,1], size=(n_rows, n_cols))

       def energy(self):
           # Returns the total energy of the current lattice configuration.
           J = 1.0 # Including this factor allows for potential modifications
           summed_spins = 0
           for i in range(len(self.lattice)): # Cycles through the rows of the spin array
               for j in range(len(self.lattice[i])): # Cycles through the elements in the ith row of the spin array
                   if i + 1 < self.n_rows and j + 1 < self.n_cols: # Multiplies i,jth element with i+1,j and i,j+1 and adds to sum (applies to elements with two neighbours in the array)
                       summed_spins += self.lattice[i,j] * self.lattice[i+1,j]
                       summed_spins += self.lattice[i,j] * self.lattice[i,j+1]
                   elif i + 1 == self.n_rows and j + 1 < self.n_cols: # Multiplies i,jth element with 0,j (due to periodic nature of the lattice) and i,j+1 (applies to elements in final row)
                       summed_spins += self.lattice[i,j] * self.lattice[0,j]
                       summed_spins += self.lattice[i,j] * self.lattice[i,j+1]
                   elif i + 1 < self.n_rows and j + 1 == self.n_cols: # Multiplies i,jth element with i+1,j and i,0 (applies elements in final column)
                       summed_spins += self.lattice[i,j] * self.lattice[i+1,j]
                       summed_spins += self.lattice[i,j] * self.lattice[i,0]
                   else: # Multiplied i,jth element with 0,j and i,0 (final array element - bottom right corner)
                       summed_spins += self.lattice[i,j] * self.lattice[0,j]
                       summed_spins += self.lattice[i,j] * self.lattice[i,0]
           energy = - J * summed_spins
           return energy

       def magnetisation(self):
           # Returns the total magnetisation of the current lattice configuration.
           magnetisation = 0
           for i in range(len(self.lattice)): # Cycles through the rows of the spin array
               for j in range(len(self.lattice[i])): # Cycles through the elements in the ith row of the spin array
                   magnetisation += self.lattice[i,j] # Updates the running total of the magnetisation by adding i,jth value
           return magnetisation
       
       def montecarlostep(self, T):
           # This function performs a single Monte Carlo step
           energy = self.energy() # Energies calculated using this function are in units of the Boltzmann constant
           magnetisation = self.magnetisation()
           # This section randomly selects the coordinate of the spin to be flipped, then flips it 
           #- generates new configuration
           random_i = np.random.choice(range(0, self.n_rows))
           random_j = np.random.choice(range(0, self.n_cols))
           self.lattice[random_i, random_j] = - self.lattice[random_i, random_j]
           energy1 = self.energy() # Calculates the energy of the new configuration
           magnetisation1 = self.magnetisation() # Calculates the magnetisation of the new configuration
           # The following line selects a random number in the range [0,1) for step 4 below
           random_number = np.random.random()
           # This section performs step 4 of the Monte Carlo cycle
           deltaE = energy1 - energy # Calculates energy difference between 
           #new and original spin arrays
           if deltaE < 0:# Accepts the new configuration and sets energy and magnetisation equal to the properties for the new arrangement
               energy = energy1 
               magnetisation = magnetisation1 
           elif random_number <= exp(- deltaE / T): # Allows for statistical inclusion of arrangements with deltaE > 0, with a 'large' Boltzmann factor
               energy = energy1
               magnetisation = magnetisation1
           else: # Rejects the new arrangement - flips the spin back to that in the original self.lattice
               self.lattice[random_i, random_j] = - self.lattice[random_i, random_j]
           self.E += energy
           self.E2 += energy**2
           self.M += magnetisation
           self.M2 += magnetisation**2 # self.X updates the running total for property X
           self.n_cycles += 1 # Updates the running total for number of Monte Carlo cycles
           return energy, magnetisation
       
       
       def statistics(self):
           # This function calculates the values for the averages of E, E*E (E2), M, M*M (M2), and returns them
           avgE = self.E / self.n_cycles
           avgE2 = self.E2 / self.n_cycles
           avgM = self.M / self.n_cycles
           avgM2 = self.M2 / self.n_cycles # These are the average values for each property (divided by number of MC cycles)
           return avgE, avgE2, avgM, avgM2