import numpy as np
from math import exp
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):
# 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]
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
