Rep:MonteCarloAndStatistics.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
        
        sum_sisj  = []
        
        for i in range(self.n_rows): 
            for j in range(self.n_cols):
            
                if i < self.n_rows -1:
                    sum_sisj.append(self.lattice[i,j] * self.lattice[i+1,j])
                else:
                    sum_sisj.append(self.lattice[i,j] * self.lattice[0,j])

                if i > 0:
                    sum_sisj.append(self.lattice[i,j] * self.lattice[i-1,j])
                else:
                    sum_sisj.append(self.lattice[i,j] * self.lattice[self.n_rows -1,j])
                
                if j < self.n_cols -1:
                    sum_sisj.append(self.lattice[i,j] * self.lattice[i,j+1])
                else:
                    sum_sisj.append(self.lattice[i,j] * self.lattice[i,0])  
                
                if j > 0:
                    sum_sisj.append(self.lattice[i,j] * self.lattice[i,j-1])
                else:          
                    sum_sisj.append(self.lattice[i,j] * self.lattice[i,self.n_cols -1])
        
        sum_sum_sisj = sum(sum_sisj)
        
        E = -0.5 * 1.0 * sum_sum_sisj
                
        return E

    def magnetisation(self):
        "Return the total magnetisation of the current lattice configuration."
        
        magnetisation_list = []
        
        for i in range(self.n_rows): 
            for j in range(self.n_cols):
                magnetisation_list.append(self.lattice[i,j])

        magnetisation = sum (magnetisation_list)
            
        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