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
