AcceleratedEnergyAndMagnetisation.py

import numpy as np

class IsingLattice:
    E = 0.0
    E2 = 0.0
    M = 0.0
    M2 = 0.0
    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):
        "Return the total energy of the current lattice configuration."
        energy = 0.0
  
        #don't double count, and don't x1/2
        i_forward = np.roll (self.lattice,1,1) 
  
        j_forward = np.roll (self.lattice,1,0) 
                
        sum_sisj = np.sum(np.multiply(i_forward,self.lattice))
        sum_sisj+=np.sum(np.multiply(j_forward,self.lattice))
        
        E = -1.0 * sum_sisj
                
        return E

    def magnetisation(self):
        "Return the total magnetisation of the current lattice configuration."
        magnetisation = np.sum (self.lattice)           
        return magnetisation

     def montecarlostep(self, T):
        #complete this function so that it performs a single Monte Carlo step
        #STEP 1
        E0 = self.energy()

        #STEP 2
        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]*=-1
    
        #STEP 3
        E1 = self.energy()
    
        #STEP 4
        delta_E = E1 - E0
    
        r = np.random.random()
        if delta_E > 0:
            if r > np.exp(-delta_E/T):
                self.lattice[random_i, random_j]*=-1
            
        #for statistics function

        M = self.magnetisation()
        E = self.energy()
        E2 = E*E
        M2 = M*M

        self.E+=E
        self.M+=M
        self.E2+=E2
        self.M2+=M2
        self.n_cycles+=1

        return E, M

            
    def statistics(self):
            # complete this function so that it calculates the correct values for the averages of E, E*E (E2), M, M*M (M2), and returns them    
        av_M = self.M/self.n_cycles
        av_E = self.E/self.n_cycles
        av_M2 = self.M2/self.n_cycles
        av_E2 = self.E2/self.n_cycles       
        return av_E, av_E2, av_M, av_M2, self.n_cycles