View source for Rep:Mod:DMS3053-5

    import numpy as np
    from math import exp, sqrt

    class IsingLattice:

        E = 0.0
        E2 = 0.0
        M = 0.0
        M2 = 0.0
        E_list = []
        M_list = []

        n_cycles = 0
        skip_region = 0 # Defines the number of initial cycles not to be averaged for above quantities

        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 later modifications
            summed_spins = np.sum(np.multiply(self.lattice, np.roll(self.lattice, 1, axis = 0)))
            # Multiplies a spin with its vertical neighbour, via rolling the array, then sums these results
            summed_spins += np.sum(np.multiply(self.lattice, np.roll(self.lattice, 1, axis = 1)))
            # Multiplies a spin with its horizontal neighbour, then sums these results
            energy = - 1.0 * J * summed_spins
            return energy 

        def magnetisation(self):
            # Returns the total magnetisation of the current lattice configuration.
            magnetisation = np.sum(self.lattice) # Sums all of the elements of the generated lattice
            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]
            if self.n_cycles >= self.skip_region: # Begins adding values to running sums after self.skip_region cycles
                self.E += energy
                self.E2 += energy**2
                self.M += magnetisation
                self.M2 += magnetisation**2 # self.X updates the running total for property X
                self.E_list.append(energy)
                self.M_list.append(magnetisation) # Appends values to list for standard error calculation
            else:
                pass
            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
            averaging_cycles = self.n_cycles - self.skip_region
            avgE = self.E / averaging_cycles
            avgE2 = self.E2 / averaging_cycles
            avgM = self.M / averaging_cycles
            avgM2 = self.M2 / averaging_cycles
            std_errorE = np.std(self.E_list) / sqrt(averaging_cycles)
            std_errorM = np.std(self.M_list) / sqrt(averaging_cycles)
            return avgE, avgE2, avgM, avgM2, self.n_cycles, std_errorE, std_errorM