View source for Rep:Mh3413CMP

=Monte Carlo Simulations of the Ising Model=

==Abstract==

By applying a Monte Carlo algorithm to a square lattice of magnetic atoms in the Ising model, the properties of a two-dimensional ferromagnetic system were investigated. It was found that below the Curie temperature the whole lattice would spontaneously align to produce a net magnetization in the presence of no external magnetic field. At a reduced temperature of T≈2.0 for the majority of lattice sizes, a sharp decrease in the magnetisation, along with an increase in the volatility of magnetization and energy was observed, corresponding to the <math>T_c</math> Curie temperature. As the temperature of the system increased for all lattice sizes, the energy per spin plateaued to a maximum value and the magnetization equilibrated to a value of zero, resembling a <math>y = -\sqrt[3]{|x|}</math> function. Using curve polynomial fitting techniques for the experimental data, the <math>T_c</math> for all the lattice sizes were able to be established and used to establish the value of <math>T_{C,\infty}</math>. The <math>T_c</math> Curie temperature was observed to decrease with increasing lattice size.

==Overview of Ising Model==
The Ising model is a simple physical model describing the ferromagnetic properties of a regular square lattice system, in the presence of no external magnetic field. In this model, the vectors of the atomic magnetic dipoles are restricted to states of pointing up(+1) or down(-1) in the two dimensional structure. The equation in Figure 1 describes the interaction energy in a system upon the assumption of:
* Only short range, nearest neighbor interactions between the adjacent lattice cells contributing to the total energy 
* Periodic boundary conditions applied to the lattice structures, causing them to wrap on them themselves
* No external magnetic field applied 

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
! <math> E=-\frac{1}{2} J \sum_i^N \sum_{j\  \in\  \mathrm{neighbours}\left(i\right)} s_i s_j</math>

|-
|'''Figure 1''': The energy of the interactions with the adjacent spins of lattice cell, when no magnetic field is applied

|}
The J is a coupling constant and is a quantitative description of the level of interaction between the adjacent spin cells. It was assumed that J=1.0 in the follow sections, resulting in the energy being in reduced units of <math>k_B</math>.
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
! [[image:Nearestneighbourmodel.png|center|500px|molecule1]]

|-
|'''Figure 2''': Number of nearest neighbours for each lattice dimension

|}

'''<big>TASK:</big>''' For the one dimensional case, there are two nearest neighbours with respect to the central atom. For 2- and 3-dimensional lattice systems, there are 4 and 6 respective neighbors. Overall, in a lattice with  periodic boundary conditions, each lattice cell would have <math> 2d </math>, nearest neighbours (where d=the number of dimensions). In the following model, it was assumed that the lattice was a 2-dimensional system. 
<center>
<math>{Number\ of\ neighbouring \ atoms \ = 2\times D}</math>
</center>

<math> {where \ D = \ the \ number \ of \ dimensions}  </math>

Using the equation in figure 1, '''the lowest energy state would correspond to a system where all the spins of the atoms are in the same state (either up or down)''': The factor of 1/2 is to account for the double counting of the interactions between the lattice atoms.

<center>
<math> Energy \ = \left(-\frac{1}{2}J \ for \ each \ interaction \right)\times\left(the \ number \ of \ neighbours \right)\times\left(number \ of \ atoms \ in \ lattice\right) </math>

<math> Energy \ = \left(-\frac{1}{2}J\right)\times\left( 2D\right)\times\left( N\right) </math>

<math>Energy_{ all\ spins \ aligned}\ =\ -DNJ</math>

</center>
'''<big>TASK:</big>''' '''As all the atoms have the same spin state, the multiplicity of the system is 1''', there is only one way to arrange the system. This system would be present at a temperature of 0K. 
<center>
<math>\Omega\ =\ 1</math>

<math>Entropy \ = \ k_B \ln 1  </math>

<math>Entropy \ =\ 0</math>
</center>
 
This result is consistent with the 3rd law of thermodynamics which states that '''"the entropy of a perfectly crystalline system is 0 at a temperature of 0K"'''. 

'''<big>TASK:</big>'''  The change in spin from -1 to +1 for one of the cells can be thought of as two processes: Change from -1 to 0 and 0 to +1 spin
For a (<math>D=3,\ N=1000,\ ie\ L=10 </math>) system in the lowest energy configuration, when one of the spins spontaneously change direction the energy change would correspond to:

[[File:ThirdYearCMPExpt-IsingSketch.png|300px|thumb|left|'''Figure 3''': Illustration of an Ising lattice in one (N=5), two (N=5x5), and three (N=5x5x5) dimensions. Red cells indicate the "up" spin", and blue cells indicate the "down" spin.]]
<center>
<math>\Delta Energy\ to\ change\ to\ 0\ spin =\left(\frac{1}{2}J\times 4D\right) </math>

<math>\Delta Energy\ to\ change\ to\ 0\ spin =\left(\frac{1}{2}J\times 4 \times 3\right)  </math>

<math>\Delta Energy\ to\ change\ to\ 0\ spin =\ +6\ J\ </math>

<math>\Delta Energy\ to\ change\ to\ -1\ from\ 1 =\ +6\times 2\ J\ </math>

<math>\Delta Energy\ =\ +12 J\ </math>
</center>

The entropy of this system would increase, as there would be an increase in the number of ways to produce the state (increase in multiplicity).

<center>
<math>\Omega\ =\ 1.00\times10^{3}</math>

<math>\Delta Entropy \ = \ k_B \ln 1.00\times10^{3}\ =\ 9.54\times10^{-23}\ J\ K^{-1} </math>
</center>


The magnetisation of the system is calculated by the sum of the parallel spins in the system. In the lowest energy configuration, where all the atoms in the lattice have the same direction of spin, the system would expect to have a maximum value of magnetisation. This corresponds to high ferromagnetic state<ref>Peter Atkins and Julio De Paula, Physical Chemistry, Oxford University Press, 5th edition, 2009, page 415<br>
<br>
</ref>. 

<center>
<math>Total Magnestisation = \sum_i s_i</math>
</center>

'''<big>TASK:</big>''' For figure 3, the magnetisation for the 1D and 2D case would be expected to be <math>\ Total\ Magnestisation\ =\ +1\ up\ spin</math>, as the number of up spin states are greater than the number of down spin states. As there is no external magnetic field that is present in the system, the sign of magnetization is not significant. 


'''<big>TASK:</big>''' For a (<math>D=3,\ N=1000,\ ie\ L=10 </math>) system where the atomic magnetic dipoles do not have sufficient energy to repel against their neighbours, it would be expected that the cells would be all aligned in the same direction. This would be most favourable in energetic terms and lead to a maximum magnetisation value of <math>Total\ Magnestisation = 1000</math>

==Calculating the energy and magnetisation==
'''<big>TASK:</big>''' It was assumed that J=1.0 in the follow sections. When J>0 there is a energetic drive for the spins to align as this minimizes the energy of the system.<ref>Maginn EJ, Shah JK. Ising Lattice [Internet]. Department of Chemical and Biomolecular Engineering, University of Notre Dame; [cited 2016]. Available from: <nowiki>http://www.peq.coppe.ufrj.br/thermo/lecture_notes/ising.pdf</nowiki></ref>
The figure 4 below, presents the code that '''initally''' was used to calculate the energy and the magnetisation of the system. Please refer to the '''Accelerating Code''' section for the modified version of the script.
===Energy===
The figure 5 gives a visual illustration of the actions undertaken on the randomly generated original matrix. These actions were undertaken to ensure the final energy calculations of the matrix took into consideration the assumption of the periodic boundary conditions, and for the simplification of the energy calculations.

The 1st matrix that is printed in figure 5 is the original matrix random generated from [-1,1] states. The energy function takes copies of the 1st and last rows of the original matrix and creates a new matrix, where the copied last row of the original matrix, becomes the 1st row of the new matrix. Likewise, the 1st row of the original matrix becomes the last row of the new matrix. The original x matrix within this additional 'wall' of numbers was left untouched. This same operation was applied to the columns, such that a 6x6 matrix was formed. The calculation of the interaction energies was only applied to the original 4x4 matrix within the 'wall' of the copied values, using the equation in figure 1.

===Magnetisation===
As the matrix values represent the magnetic spin of atoms in a lattice, the values are restricted to [-1,1]. The sum of the spins values in the lattice is equal to the total magnetisation of the system; a simple sum(sum(x)) was applied to the matrix to return the total summation value. The model assumes that in the presence of no external magnetic field, the lattice can still generate a net magnetisation through the long range correlations.

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
! [[image:EnergyMagnetisationmh3413CodingCalculations.PNG|center|600px|molecule1]]
! [[image:EnergyMagnetisationmh3413CodingManipluation.PNG|center|500px|molecule1]]
! [[image:EnergyMagnetisationmh3413.PNG|center|500px|molecule1]]

|-
|'''Figure 4''': Method used to calculate the energy of the system
|'''Figure 5''': Initial 4x 4 matrix used to create a 6x6 matrix, which was used to calculate the total energy
|'''Figure 6''': Running the ILcheck.py script on the Ising Lattice model

|}
'''<big>TASK:</big>''' The main engine of the model that initally was used to calculate the energy and magnetisation. Part of the IsingLatticeModel script.
<center>
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|   <pre>
    def energy(self):
        "Return the total energy of the current lattice configuration."
        
        "Addition the lattice rows to take into 4 neighbours per atom"
        a=(np.insert(self.lattice, self.n_rows, self.lattice[0], axis=0))
        b=(np.insert(a, 0, self.lattice[self.n_rows-1], axis=0))
        
        "Addition the lattice columns to take into 4 neighbours per atom"
        c=(np.insert(b,self.n_cols,[row[0] for row in b], axis=1))
        EditedLattice=(np.insert(c,0,[row[self.n_cols-1] for row in b], axis=1))
        
        #print(self.lattice)
        #print(EditedLattice)
        #print(range(1,self.n_rows))
        #print(range(1,self.n_cols))
        
        J=1.0
        Energy=0
        
        for i in range(1,self.n_rows+1):
            for j in range(1,self.n_cols+1):
                Energy=Energy-0.5*J*EditedLattice[i][j]*EditedLattice[i][j+1]+-0.5*J*EditedLattice[i][j]*EditedLattice[i][j-1]+-0.5*J*EditedLattice[i][j]*EditedLattice[i+1][j]+-0.5*J*EditedLattice[i][j]*EditedLattice[i-1][j]
        
        return Energy

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

</pre>
|}
</center>

==Introduction to Monte Carlo simulation==
'''<big>TASK:</big>''' '''Please refer to the Accelerated Code section to see the modified IsingLattice code'''

'''<big>TASK:</big>''' Due to the spin up(+1) and spin down(-1) possibilities for each atom, there are a increasing number of configurations possible with the increasing number of atoms.
<center>
<math> Configurations\ Possible = 2^N </math>
</center>

For a 100 spin system and computing power of <math>1\times 10^9</math> calculations per second: 
<center>
<math> Processing\ Time\ = \frac{2^{100}}{1\times 10^9}</math>

<math> Processing\ time\ for\ one\ value\ of\left\langle M\right\rangle_T\ ={1.27\ (2\ dp)\times 10^{21}\ seconds}</math>
</center>
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
! [[image:energyoptimizationmontecarlocode.png|center|350px|molecule1]]
! [[image:mh3413energyoptimizationforonestep.PNG|center|300px|molecule1]]
! [[image:Energyrunningmh3413v2.png|center|330px|molecule1]]
! [[image:energyrunningmh3413.png|center|330px|molecule1]]

|-
|'''Figure 7''': Populated Monte Carlo and Statistics() functions 
|'''Figure 8''': A single cycle of the Monte Carlo algorithm was implemented along with the statistics() function. The following quantities <math><E>, <E^2>, <M>, <M^2></math>, number of Monte Carlo steps that have elapsed, were returned by the function.
|'''Figure 9''': Plot of reduced energy per spin and reduced magnetisation per spin, against number Monte Carlo of steps
|'''Figure 10''': Plot of reduced energy per spin and reduced magnetisation per spin, against number Monte Carlo of steps -  equilibrium system

|}

'''<big>TASK:</big>''' Figure 7 to 10 captures the development and testing of the montecarlocycle() and statistics() function. Figures 9 and 10 ultilised the ILanim.py script to animate a real time Monte Carlo simulation. Figure 9 presents the minimum energy and maximum magnetisation value of the system in equilibrium state (ie graphs of the energy and magnetisation per spin). It is expected that when <math>T < T_C</math>, there would be spontaneous magnetisation (i.e. <math>\left\langle M\right\rangle \neq 0</math>) as the exchange and long range correlation effects would outweigh the thermal fluctuations, causing the neighboring atoms to favour aligning in the same direction. Above the <math>T_C</math>, the thermal energy in the system is sufficient to overcome any exchange and long range correlation effects.

Figure 10 is a system in an equilibrium state.

<math>Averaged\ quantities:</math>
<center>
<math><E> =  -1.70271773681 </math>

<math><E^2> =  3.02165944851 </math>

<math><M> =  -0.800301970757 </math>

<math><M^2> =  0.70215154164 </math>

<math>Number\ of\ Steps\ =1274 </math>
</center>

==Accelerating the code==
'''<big>TASK:</big>''' Using the script ILtimetrial.py, the time required to perform 2000 Monte Carlo steps in the IsingLattice.py script was calculated.
On average the original code took:
<center>
<math> Time\ to\ calculate\ 2000\ steps =\ 40.31605199s \plusmn\ 0.271653419s</math>
</center>

The code was optimized using the roll() and multiply() functions. The original matrix was multiplied with the original matrix shifted by one unit in the x-axis and the original matrix shifted by one unit in the y-axis. This operation was able to capture the same neighbouring interactions as the previous code, but was more time efficient. 
On average the optimized code took:
<center>
<math> Time\ to\ calculate\ 2000\ steps =\ 0.635318786s \plusmn\ 0.009924127s</math>
</center>

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
! [[image:mh3413Runtimebeforeoptimization.PNG|center|500px|molecule1]]
! [[image:Mh3413Runtimeafteroptimization.PNG|center|500px|molecule1]]

|-
|'''Figure 11''': Before Optimization of energy and magnetization code 
|'''Figure 12''': After Optimization of energy and magnetization code 

|}

===Optimized Ising Lattice Code===
'''<big>TASK:</big>''' The main engine of the model that applies the Monte Carlo steps and calculates the energy and magnetisation of the system
<center>
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|   <pre>
    import numpy as np
    import math as math

    class IsingLattice:

        E = 0.0
        E2 = 0.0
        M = 0.0
        M2 = 0.0
        NoIgnore = 160000
        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."

            J=1.0
            Energy=0

            NewL=np.multiply(np.roll(self.lattice, 1, axis=1),self.lattice)
            NewL2=np.multiply(np.roll(self.lattice, 1, axis=0),self.lattice)

            Energy=(sum(sum(NewL))+sum(sum(NewL2)))*-J

            return Energy

        def magnetisation(self):
            "Return the total magnetisation of the current lattice configuration."

            magnetisation=sum(sum(self.lattice))

            return magnetisation

        def montecarlostep(self, T):
            # complete this function so that it performs a single Monte Carlo step

            latticecopy=np.copy(self.lattice)
            kb=1

            energy0 = self.energy()

            #the following two lines will select the coordinates of the random spin for you
            #Picking the random spin system in the lattice
            random_i = np.random.choice(range(0, self.n_rows))
            random_j = np.random.choice(range(0, self.n_cols))

            #Flipping operation of the spin of this random system
            if self.lattice[random_i][random_j]<0:
                self.lattice[random_i][random_j]=1
            else:
                self.lattice[random_i][random_j]=-1

            energy1 = self.energy()

            #Calculating the energy difference in the E1 and E0 states
            #the following line will choose a random number in the rang e[0,1) for you
            if (energy1-energy0)>0:
                random_number=np.random.random()
                if random_number<=math.exp(-(energy1-energy0)/(kb*T)):
                    pass
                else:
                    self.lattice=latticecopy
            else:
                pass

            self.n_cycles = self.n_cycles+1
            if self.n_cycles>=self.NoIgnore:
                self.E = self.E+self.energy()
                self.E2 = self.E2+(self.energy())**2
                self.M = self.M+self.magnetisation()
                self.M2 = self.M2+(self.magnetisation())**2
            else:
                self.E = 0
                self.E2 = 0
                self.M = 0
                self.M2 = 0

            return(self.energy(),self.magnetisation())

        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
            if self.n_cycles>=self.NoIgnore:
                AverageE=self.E/(self.n_cycles-self.NoIgnore)
                AverageE2=self.E2/(self.n_cycles-self.NoIgnore)
                AverageM=self.M/(self.n_cycles-self.NoIgnore)
                AverageM2=self.M2/(self.n_cycles-self.NoIgnore)
            else: #if not above the number of cycles ignored, the system will ignore the value
                AverageE=0
                AverageE2=0
                AverageM=0
                AverageM2=0

            # take into consideration the ignored cycles in the printing of the number of steps undertaken
            return (AverageE,AverageE2,AverageM,AverageM2,(self.n_cycles-self.NoIgnore))
</pre>
|}
</center>

==The effect of temperature==
As the model starts from a random set of lattice spins, there is a certain number of steps required in order for the system to reach an energetic and magnetic thermal equilibrium. There are two main variables that effect this value:
*Lattice Size
*Temperature
'''Note: In each of the captions of the figures below, there is an approximated number of steps required to reach the equilibrium state'''

For figures 13 to 16 the temperature was kept constant at T=1.85, while the lattice size was increased. As seen below, with the increasing lattice sizes, the number of steps required to reach the equilibrium position increased. For example, the magnetisation and energy in the 10x10 lattice structure equilibrised <10000 steps, while in the 64x64 system the energy reached equilibrium at approx 160,000 steps.

'''Changing Lattice Size'''

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
! [[image:Figure_1.10by10t1.85K.png|center|400px|molecule1]]
! [[image:Figure_1.16by16t1.85K.png|center|400px|molecule1]]
! [[image:Figure_1.32by32t1.85K.png|center|400px|molecule1]]
! [[image:Figure 1.64by64t1.85K.png|center|400px|molecule1]]
|-
|'''Figure 13''': Final state of the 10x10 lattice at T=1.85, runtime = 7500, equilibrium state<10000
|'''Figure 14''': Final state of the 16x16 lattice at T=1.85, runtime = 7500, equilibrium state≈10000
|'''Figure 15''': Final state of the 32x32 lattice at T=1.85, runtime = 30,000, equilibrium state≈120000
|'''Figure 16''': Final state of the 64x64 lattice at T=1.85, runtime = 50,000, equilibrium state≈140000
|}

For figures 17 to 20, the lattice size was kept constant at 32x32, while the reduced temperatures were varied between T=0.25-5.0. With increasing temperatures of the system, the number of steps required to reach the equilibrium position in the system increased. For figure 20, T=5 is significantly greater than the Curie temperature for the system and thus the system reaches magnetic and energy equilibrium instantaneously as the thermal fluctuations eliminate all the correlation, enthaplic and entropic factors. 

'''Changing Temperature'''
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"

! [[image:32x32T0.25run175000mh3413CMP.png|center|400px|molecule1]]
! [[image:32x32T0.8run175000mh3413CMP.png|center|400px|molecule1]]
! [[image:32x32T1.85run175000mh3413CMP.png|center|400px|molecule1]]
! [[image:32x32T5run175000mh3413CMP.png|center|400px|molecule1]]
|-
|'''Figure 17''': Final state of the 32x32 lattice at T=0.25, runtime = 175,000, equilibrium state≈90,000
|'''Figure 18''': Final state of the 32x32 lattice at T=0.80, runtime = 175,000, equilibrium state≈155,000
|'''Figure 19''': Final state of the 32x32 lattice at T=1.85, runtime = 175,000, equilibrium state≈160,000
|'''Figure 20''': Final state of the 32x32 lattice at T=5, runtime = 175,000, equilibrium state≈Immediately goes to equilibrium
|}

'''<big>TASK:</big>'''  It was found that for each lattice size and temperature, the system has a different number of cycles needed to obtain an equilibrium state. Based on the information above, for the worst case scenarios (T= approx 1.85 and very large lattice sizes eg 64x64), '''all the lattices reach an equilibrium by 160,000 steps'''. Therefore, the '''first 160,000 cycles of the Monte Carlo steps were chosen to be ignored in the calculation of the averages'''. This result led to the average values of magnetisation and energy having '''greater magnitudes and more accurate values''' for approximating the equilibrium system. A greater number of steps could have been ignored, and this may have lead to increased accuracy in the calculation of the variables, however in terms of the computational time required to calculate the these variables (eg for runtime=240,000) and the benefit obtained, it was decided that an optimal value of steps to ignore was 160,000.

The Ising Lattice script model no longer considers the energy and magnetisation values for the steps up to 160,000 ie until the system reaches an equilibrium point.

<math>Averaged\ quantities:</math>
<center>
<math><E> =   -1.99620390625 </math>

<math><E^2> =  3.98527241211 </math>

<math><M> =  0.999049609375 </math>

<math><M^2> =  0.998119226074 </math>

<math>Number\ of\ steps\ set\ =240,000 </math>
</center>

==Running over a range of temperatures/Scaling the System Size==

'''<big>TASK:</big>''' Using the ILtemperaturerange.py script, the average energy and magnetisation was plotted against the reduced temperature. Taking into account the 160,000 ignored simulations, the '''runtime was set to 240,000''' and the '''dt steps adjusted to a lower value of 0.05''' over a '''temperature range from T=0.25-5.0'''. This allowed the fine details of the fluctuations in energy and magnetisation for the system to be captured. The error bars at each temperature point, are calculated using the fluctuations in values from the runs conducted. The variability originates from the variations that occur in the <math>k_BT</math>.

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"

! [[File:4by4figure_1-1mh3413.png|center|400px|molecule1]]
! [[File:8by8mh3413.PNG|center|400px|molecule1]]
! [[File:16by16mh3413.PNG|center|400px|molecule1]]
! [[File:32x32dt0.05mh3413.png|center|400px|molecule1]]

|-
|'''Figure 21''': Average energy and magnetisation simulations against temperature for the 4x4 lattice
|'''Figure 22''': '''<big>TASK:</big>''' Average energy and magnetisation simulations against temperature for the 8x8 lattice
|'''Figure 23''': Average energy and magnetisation simulations against temperature for the 16x16 lattice
|'''Figure 24''': Average energy and magnetisation simulations against temperature for the 32x32 lattice
|}

The <math>T_c</math> Curie Temperature is the '''temperature at which the system undergoes a phase transition from an ordered to a disordered system''' and loses the net magnetization. The magnetisation changes rapidly at this point. Above the <math>T_c</math> Curie Temperature, the changing of the spins from an ordered system to a completely disordered system – there is no energy penalty involved in this. In addition, there are long range correlation interactions present in the system which result a spin change, increasing the probability of triggering a spin change in another atom. Figure 25 captures the energy per spin against the reduced temperature. The energy rapdily increases upon reaching the Curie temperature; this graph resembles a <math>y = \sqrt[3]{|x|}</math> function. With increasing lattice size, the energy per spin was observed to increase, with the 8x8,16x16 and 32x32 having very similar shapes.

The Curie temperature was observed to increase with increasing lattice sizes in figure 26, but contradicted the results of figure 25 in which the point of inflexion (steepest gradient for dE/dT) corresponds to the Curie temperature; the Curie temperature decrease with increasing lattice size. It is believed that for figure 26, the small size of the lattice caused inaccuracies in the simulations, due to the each lattice cell being so close to one another. The change in a spin of one cell in this lattice would dramatically increase the probability of the neighboring lattice cells changing spins. Furthermore, if the lattice cell spends half of its time in a total spin up state and the other half of its time in a total spin down state, the average of the system would end up being around zero.

'''<big>TASK:</big>''' Furthermore, the increased volatility around the <math>T_c</math> Curie Temperature point decreases in significance with increasing lattice sizes (Figure 26). The temperature range at which the phase transition takes place, decreases with increasing lattice size. For a 32x32 system, the phase transition takes place almost instantly with insignificant volatility, while in the 4x4 system the transition takes place over approximately T=3 due to the fluctuation within the system. '''It is believed that for a lattice size of 16x16, the system is large enough to capture the long range fluctuations due to the sharp transition and the fluctuations observed'''.  

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"

! [[File:Figure_1energyperspinmh3413.png|center|500px|molecule1]]
! [[File:Figure_2magnetisationmh3413.png|center|500px|molecule1]]


|-
|'''Figure 25''': Reduced energy per spin versus reduced temperature
|'''Figure 26''': Reduced magnetisation per spin versus reduced temperature

|}

===Running over a range of temperatures===
'''<big>TASK:</big>'''  Plotting of the average energy and magnetisation for each temperature, with error bars, for an 8x8 lattice
<center>
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|   <pre>
from IsingLattice import *
from matplotlib import pylab as pl
import numpy as np

n_rows = 8
n_cols = 8
il = IsingLattice(n_rows, n_cols)
il.lattice = np.ones((n_rows, n_cols))
spins = n_rows*n_cols
runtime = 240000
times = range(runtime)
temps = np.arange(0.25, 5.0, 0.05)

#Empty arrays to be filled
energies = []
magnetisations = []
energysq = []
magnetisationsq = []
stdenergy = []
stdmagnetisation=[]
for t in temps:
    for i in times:
        if i % 100 == 0:
            print(t, i)
        energy, magnetisation = il.montecarlostep(t)
    aveE, aveE2, aveM, aveM2, n_cycles = il.statistics()
    energies.append(aveE)
    energysq.append(aveE2)
    magnetisations.append(aveM)
    magnetisationsq.append(aveM2)
    
    #reset the IL object for the next cycle
    il.E = 0.0
    il.E2 = 0.0
    il.M = 0.0
    il.M2 = 0.0
    il.n_cycles = 0
#Calculation in the energy and magnetisation errors of the system by calculating the Var(x) and then the standard deviations
energyerror=np.array(energysq)-(np.array(energies)**2)
magnetisationserror=np.array(magnetisationsq)-(np.array(magnetisations)**2)
stdmagnetisation=np.sqrt(magnetisationserror)
stdenergy=np.sqrt(energyerror)

#Ploting figure and formatting issues
fig = pl.figure()
enerax = fig.add_subplot(2,1,1)
enerax.set_ylabel("Energy per spin")
enerax.set_xlabel("Temperature")
enerax.set_ylim([-2.1, 2.1])
magax = fig.add_subplot(2,1,2)
magax.set_ylabel("Magnetisation per spin")
magax.set_xlabel("Temperature")
magax.set_ylim([-1.1, 1.1])

#Addition of error bars for each of the temperature points run at
enerax.errorbar(temps, np.array(energies)/spins,yerr=np.array(stdenergy)/spins)
magax.errorbar(temps, np.array(magnetisations)/spins,yerr=np.array(stdmagnetisation)/spins)
pl.show()
final_data = np.column_stack((temps, energies, energysq, magnetisations, magnetisationsq,stdmagnetisation,stdenergy))
np.savetxt("32x32-larger dt.dat", final_data)
</pre>
|}
</center>

===Scaling the System Size===
'''<big>TASK:</big>''' The code below was used to plot energy and magnetisation per spin against temperature. Currently, it plots the energy but can be easily adapted for magnestism.
<center>
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|   <pre>
import matplotlib.pyplot as plt
import numpy as np

#Extracting the data from each data file and catergorising it into variables
data2by2 = np.loadtxt("2x2.dat")
spins2by2=4
temps2by2=data2by2[:,0]
energies2by2=data2by2[:,1]
energysq2by2=data2by2[:,2]
magnetisations2by2=data2by2[:,3]
magnetisationsq2by2=data2by2[:,4]

data4by4 = np.loadtxt("4x4.dat")
spins4by4=16
temps4by4=data4by4[:,0]
energies4by4=data4by4[:,1]
energysq4by4=data4by4[:,2]
magnetisations4by4=data4by4[:,3]
magnetisationsq4by4=data4by4[:,4]

data8by8 = np.loadtxt("8x8.dat")
spins8by8=64
temps8by8=data8by8[:,0]
energies8by8=data8by8[:,1]
energysq8by8=data8by8[:,2]
magnetisations8by8=data8by8[:,3]
magnetisationsq8by8=data8by8[:,4]

data16by16 = np.loadtxt("16x16.dat")
spins16by16=256
temps16by16=data16by16[:,0]
energies16by16=data16by16[:,1]
energysq16by16=data16by16[:,2]
magnetisations16by16=data16by16[:,3]
magnetisationsq16by16=data16by16[:,4]

data32by32 = np.loadtxt("32x32.dat")
temps32by32=data32by32[:,0]
spins32by32=1024
energies32by32=data32by32[:,1]
energysq32by32=data32by32[:,2]
magnetisations32by32=data32by32[:,3]
magnetisationsq32by32=data32by32[:,4]

#graph curves colour selection
plt.gca().set_color_cycle(['red', 'green', 'blue', 'yellow', 'purple'])

#Calculation of the properties per spin and plotted
plt.plot(temps2by2, energies2by2/spins2by2)
plt.plot(temps4by4, energies4by4/spins4by4)
plt.plot(temps8by8, energies8by8/spins8by8)
plt.plot(temps16by16, energies16by16/spins16by16)
plt.plot(temps32by32, energies32by32/spins32by32)

plt.xlabel('Reduced Temperature')
plt.ylabel('Reduced Energy per spin')

plt.legend(['2x2', '4x4', '8x8', '16x16','32x32'], loc='upper right')

plt.show()
plt.savefig('test.jpg')
</pre>
|}
</center>

==Determining the heat capacity==
'''<big>TASK:</big>''' Proof of: 

<math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math> = <math>\frac{\mathrm{Var}[E]}{k_B T^2}</math>

LDS:

<math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math>

<math><E> = \frac{\sum_i^N E_i}{N}</math> Or <math><E> = \sum_iE_iP(E_i)</math>

<math><E> = \frac{1}{\sum_ie^{-\frac{E_i}{k_BT}}}\sum_i(E_ie^{-\frac{E_i}{k_BT}}) </math>

where <math>P(E_i) = \frac{1}{\sum_ie^{-\frac{E_i}{k_BT}}}\sum_i(e^{-\frac{E_i}{k_BT}}) </math>

<math>C = \frac{\partial}{\partial T}\frac{\sum_i^N E_i}{N} = \frac{\partial}{\partial T}(\frac{1}{e^{-\frac{E_i}{k_BT}}}\sum_i(E_ie^{-\frac{E_i}{k_BT}}))</math>

Using the quotient rule:

<math> \frac{\partial \left\langle E\right\rangle}{\partial T} =\frac{(\sum_ie^{-\frac{E_i}{k_BT}})\times \sum_i(\frac{(E_i)^2}{k_BT^2}e^{-\frac{E_i}{k_BT}})-\sum_i(E_ie^{-\frac{E_i}{k_BT}})\times \sum_i(\frac{E_i}{k_BT^2}e^{-\frac{E_i}{k_BT}})}{\sum_i(e^{-\frac{E_i}{k_BT}})^2}</math>

<math> \frac{\partial \left\langle E\right\rangle}{\partial T} =\frac{(\sum_ie^{-\frac{E_i}{k_BT}})\sum_i((E_i)^2e^{-\frac{E_i}{k_BT}})-\sum_i(E_ie^{-\frac{E_i}{k_BT}})^2}{\sum_i(e^{-\frac{E_i}{k_BT}})^2(k_BT^2)}</math>

<math>\frac{\partial \left\langle E\right\rangle}{\partial T} = \frac{1}{(k_BT^2)}(\frac{(\sum_i((E_i)^2e^{-\frac{E_i}{k_BT}})}{\sum_i(e^{-\frac{E_i}{k_BT}})}-\frac{\sum_i(E_ie^{-\frac{E_i}{k_BT}})^2}{\sum_i(e^{-\frac{E_i}{k_BT}})^2})
</math>

<math>\frac{\partial \left\langle E\right\rangle}{\partial T} = \frac{\sum_iE_i^2P(E_i)-(\sum_iE_iP(E_i))^2}{(k_BT^2)}</math>

<math>\frac{\partial \left\langle E\right\rangle}{\partial T} = \frac{<E^2>-<E>^2}{k_BT^2}</math>

RDS:

<math>C = \frac{Var(E)}{k_BT^2}</math>

==Calculating the heat capacity==
[[File:Figure_1-1v2mh3413.png|550px|thumb|right|'''Figure 26''': '''<big>TASK:</big>''' Heat capacity versus temperature for the different lattice sizes]]
<center>
<math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math> = <math>\frac{\mathrm{Var}[E]}{k_B T^2}</math>
</center>
'''<big>TASK:</big>''' The increased detail and volatility around the curve at the <math>T_c</math> Curie Temperature, is due to the lower value of dt=0.05 that was run in the 'Running over a range of temperatures' section. This section utilized the .dat files that were produced from the previous tasks. The heat capacity vs reduced temperature graph presents a strong peak in vicinity of the critical temperature, caused by the <math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math>, the gradient of the energy per spin vs reduced temperature rapidly increasing near the Curie temperature. The point of inflexion in this graph (energy per spin vs reduced temperature) corresponds to the maxima peak observed in the heat capacity versus temperature graph. It is shown that with increasing array size, an increase in the peak amplitude is observed.

===Calculating the heat capacity===
'''<big>TASK:</big>''' Plot showing the heat capacity versus temperature for each of the lattice sizes from the previous section
<center>
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|   <pre>
import matplotlib.pyplot as plt
import numpy as np

#Load the data from the temperature range files and calculate the variance of each of the energies using the quoted formula
data2by2 = np.loadtxt("2x2.dat")
spins2by2=4
temps2by2=data2by2[:,0]
Varenergies2by2=data2by2[:,2]-data2by2[:,1]**2

data4by4 = np.loadtxt("4x4.dat")
spins4by4=16
temps4by4=data4by4[:,0]
Varenergies4by4=data4by4[:,2]-data4by4[:,1]**2

data8by8 = np.loadtxt("8x8.dat")
spins8by8=64
temps8by8=data8by8[:,0]
Varenergies8by8=data8by8[:,2]-data8by8[:,1]**2

data16by16 = np.loadtxt("16x16.dat")
spins16by16=256
temps16by16=data16by16[:,0]
Varenergies16by16=data16by16[:,2]-data16by16[:,1]**2

data32by32 = np.loadtxt("32x32.dat")
temps32by32=data32by32[:,0]
spins32by32=1024
Varenergies32by32=data32by32[:,2]-data32by32[:,1]**2

plt.gca().set_color_cycle(['red', 'green', 'blue', 'yellow', 'purple'])

#plotting of the Variance/kT2 against the reduced temperature
plt.plot(temps2by2,Varenergies2by2/np.multiply(1.38064852e-23,temps2by2**2))
plt.plot(temps4by4, Varenergies4by4/np.multiply(1.38064852e-23,temps2by2**2))
plt.plot(temps8by8,Varenergies8by8/np.multiply(1.38064852e-23,temps2by2**2))
plt.plot(temps16by16, Varenergies16by16/np.multiply(1.38064852e-23,temps2by2**2))
plt.plot(temps32by32, Varenergies32by32/np.multiply(1.38064852e-23,temps2by2**2))
plt.xlim((0,5))
plt.ylim((-1E+25,0.8E+26))
plt.xlabel('Reduced Temperature')
plt.ylabel('Reduced Heat Capacity')

plt.legend(['2x2', '4x4', '8x8', '16x16','32x32'], loc='upper left')

plt.show()
</pre>
|}
</center>

==Locating the Curie temperature==

'''Figures 17 to 21 are the plots of the 32x32 lattice C++ data''' with the computational experimental data obtained for the <math>E, E^2, M, M^2, C</math> per spin variables. The 32x32 lattice was used for the data comparison, as the volatility and error associated with the energy and magnetisation calculations for this lattice size was significantly lower than for the smaller arrays. A legend was added to the graphs to indicate the source of data. Overall, there is a good agreement between the C++ and the experimental data, however the C++ calculation method does show advantages in figure 31. The C++ data produces a pronounced peak, indicating a clear <math>T_C</math>, while the experimental data has large variability around this temperature point. This accuracy can be attributed to the large number of steps that the C++ data undertakes. 

The figures 32-34 present the 'trial and error' method was applied to a 64x64 lattice, to produce the polynominal fit to the C++ data. It was found that a maximum polynominal value of 400 produced the best fit for the whole temperature range. However, due to the nature of the peak being very abrupt in the function, the polynominal fit function was unable to produce a satisfactory fit. When the temperature range of the polynominal fit function was reduced to the regions of the peak, the polynominal was able to produce a good fit (figure 34).

'''<big>TASK:</big>'''
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"

! [[image:32x32_energy_comparisonmh3413.png|center|330px|molecule1]]
! [[image:32x32_energy_square_comparisonmh3413.png|center|330px|molecule1]]
! [[image:32x32_magnetisation_comparisonmh3413.png|center|330px|molecule1]]
! [[image:32x32_magnetisation_squared_comparisonmh3413.png|center|330px|molecule1]]
! [[image:32x32_Heat_Capacity_comparisonmh3413.png|center|330px|molecule1]]
|-
|'''Figure 27''': '''<big>TASK:</big>''' Plotting reduced energy per spin against reduced temperature
|'''Figure 28''': '''<big>TASK:</big>''' Plotting reduced energy squared per spin against reduced temperature
|'''Figure 29''': '''<big>TASK:</big>''' Plotting reduced magnetisation per spin against reduced temperature
|'''Figure 30''': '''<big>TASK:</big>''' Plotting reduced magnetisation squared per spin against reduced temperature
|'''Figure 31''': '''<big>TASK:</big>''' Plotting reduced heat capacity per spin against reduced temperature
|}
===Comparing the C++ data and the experimental data===
'''<big>TASK:</big>''' For each lattice size, plotting the C++ data against experimental data. This current code is optimized for comparing the heat capacity values across the temperature range, but can be made to plot for <math> E, E^2, M, M^2, C</math>
<center>
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|   <pre>
import matplotlib.pyplot as plt
import numpy as np

#My experimental data,  calculation of the variance of the energy to use in the equation for the heat capacity
data32by32 = np.loadtxt("32x32.dat")
spins32by32=32*32
temps32by32=data32by32[:,0]
Varenergies32by32=data32by32[:,2]-data32by32[:,1]**2

#C++ data provided
Cdata32by32 = np.loadtxt("32x32c++.dat")
Cspins32by32=32*32
Ctemps32by32=Cdata32by32[:,0]
Cenergies32by32=Cdata32by32[:,5]
#print(Cenergies32by32)

plt.gca().set_color_cycle(['red', 'green', 'blue', 'yellow', 'pink'])

#calculation of the actual heat capacity ie mulitplied by kb + plotting the data on the same graph
plt.plot(temps32by32, (1.38064852e-23*Varenergies32by32/np.multiply(1.38064852e-23,temps32by32**2))/spins32by32)
plt.plot(Ctemps32by32,Cenergies32by32)

plt.xlabel('Reduced Temperature')
plt.ylabel('Heat Capacity per spin')

plt.legend(['32x32 Experimental', '32x32 C++ model'], loc='upper right')

plt.show()
</pre>
|}
</center>

===Comparing C++ data and experimental data + polynominal fitting===
'''<big>TASK:</big>''' A script to read the data from a particular file, and plot C vs T, as well as a fitted polynomial. The 400 degree polynomial was found to be the best fit without restricting the temperature ranges. Figures 32-33

'''<big>TASK:</big>''' Modified version of the previous script to fit the polynomial to the peak of the heat capacity ie restricting the temperature range figures 34


{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"

! [[image:64x64Polynominal_8_fitmh3413.png|center|450px|molecule1]]
! [[image:64x64Polynominal_400_fitmh3413.png|center|450px|molecule1]]
! [[image:64x64Polynominal_400_fit_Restricted_T2K_to_T3.5Kmh3413.png|center|450px|molecule1]]

|-
|'''Figure 32''': '''<big>TASK:</big>''' Polynomial 8 for a 64x64 lattice
|'''Figure 33''': '''<big>TASK:</big>''' Polynomial 400 for a 64x64 lattice
|'''Figure 34''': '''<big>TASK:</big>''' Polynomial 400 for a 64x64 lattice restricted to fit between T2 to T3.5
|
|}

<center>
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|   <pre>
import matplotlib.pyplot as plt
import numpy as np

#List for the values of the maximum heat capacity and temperature to go into
Temperature=[]
HeatCapacity=[]

length=[2,4,8,16,32,64]
name=["2x2c++","4x4c++","8x8c++","16x16c++","32x32c++","64x64c++"]

#Loop to export data intp python from various c++ data files
for i in name:
    naming= i+".dat"
    data = np.loadtxt(str(naming)) #assume data is now a 2D array containing two columns, T and C
    T = data[:,0] #get the first column
    C = data[:,5] # get the 6th column
    
    #The region which the polynominal will fit to
    Tmin = 1.5
    Tmax = 4.0
    
    #Restricting the peak picking to the temeprature peak region due to peaks at the edge of boundaries
    #Restricting Tmin and Tmax to other values if produces a better fit to the data
    if i=="2x2c++":
        Tminforpeak = 1.5
        Tmaxforpeak = 3.0
    elif i=="4x4c++":
        Tminforpeak = 1.7
        Tmaxforpeak = 3.3
    elif i=="8x8c++":
        Tminforpeak = 1.8
        Tmaxforpeak = 3.3
    elif i=="16x16c++":
        Tmin = 1.5
        Tmax = 2.9
        Tminforpeak = 1.5
        Tmaxforpeak = 2.8
    elif i=="32x32c++":
        Tmin = 1.5
        Tmax = 2.8
        Tminforpeak = 1.5
        Tmaxforpeak = 2.7
    else:
        Tmin = 1.5
        Tmax = 2.75
        Tminforpeak = 1.5
        Tmaxforpeak = 2.7

    selection = np.logical_and(T >=Tmin, T <= Tmax) #choose only those rows where both conditions are true
    peak_T_values = T[selection]
    peak_C_values = C[selection]

    #first we fit the polynomial to the data
    fit = np.polyfit(peak_T_values, peak_C_values, 400) # fit a third order polynomial

    #now we generate interpolated values of the fitted polynomial over the range of our function
    T_range = np.linspace(Tminforpeak , Tmaxforpeak, 1000) #generate 1000 evenly spaced points between T_min and T_max
    fitted_C_values = np.polyval(fit, T_range) # use the fit object to generate the corresponding values of C
    #print(T)

    Cmax = np.max(fitted_C_values)
    Tmax = T_range[fitted_C_values == Cmax]
    
    #Exporting the Heat Capacity maximum, temperature value and the corresponding lattice size to a .dat file
    Temperature.append(Tmax)
    HeatCapacity.append(Cmax)

final_data = np.column_stack((np.array(length),np.array(Temperature), np.array(HeatCapacity)))
np.savetxt("Heat Capacity Ising Model.dat", final_data)
</pre>
|}
</center>

===Finding the peak in C===
'''<big>TASK:</big>''' 
<center>
<math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math>
</center>

Using the equation above, the 1/(lattice length) was plotted against the reduced temperature. The resulting linear relationship was then extrapolated and calculated at T=0; corresponds to the value of <math>T_{C,\infty}</math>. This was compared with Onsager(1944<ref name=":0">Onsager, L. Crystal Statistics I A Two-Dimensional
Model with an Order-Disorder Transition. Phys Rev. [Online] 1944;65(3-4): 117-149. Available from:
<nowiki>http://link.aps.org/doi/10.1103/PhysRev.65.117, doi :10.1103/PhysRev.65.117, [Accessed
15 November 2016].</nowiki></ref>) literature values and had very good agreement. 

The major sources of error are due to the assumptions made in the model:
* The noise at the peak value creates uncertainty in the plot fitting
* Assumption that the relationship in figure 35 is linear
* Required to extrapolation data for T=0 - uncertainty in what the actual data would be in this region
* There are more than 4 nearest neighbours around a lattice cell; the longer range of interactions are not considered

<center>
<math>T_{C,\infty}(C^{++} experimental) = 2.27837987241</math>

<math>T_{C,\infty}(literature reports) = 2.27</math><ref name=":0" />
</center>

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"

! [[image:Figure_1-1CalculationTCinfinite.png|center|550px|molecule1]]
! [[image:CaptureTCinfintevalue&code.PNG|center|600px|molecule1]]


|-
|'''Figure 35''': Plot of the reduced temperature against 1/(lattice length)
|'''Figure 36''': Code used to produce Figure 35, in addition to the extrapolation of the <math>T_{C,\infty}</math> 

|
|}

'''<big>TASK:</big>''' A plot that uses the scaling relation to determine <math>T_{C,\infty}</math>
<center>
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|   <pre>
import matplotlib.pyplot as plt
import numpy as np

data = np.loadtxt("Heat Capacity Ising Model.dat")
L = data[:,0] #get the first column, Size of the Lattice
T = data[:,1] #get the second column, Temperature of Maximum C
C = data[:,2] #get the third column, Maximum value of C 

plt.plot(1/L,T)

plt.xlabel('1/(Lattice Size)')
plt.ylabel('Reduced Heat Capacity')

#extrapolation for 1/lattice size values
Latticemax = np.ndarray.max(1/L)
Latticemin = 0

#extrapolation using an approximate linear equation
fit = np.polyfit(1/L, T, 1) # fit a first order polynomial ie linear
T_range = np.linspace(Latticemin , Latticemax, 1000) #generate 1000 evenly spaced points between T_min and T_max
fitted_Tinfinite_values = np.polyval(fit, T_range) # use the fit object to generate the corresponding values of C

#Plotting the fitted equation
plt.plot(T_range, fitted_Tinfinite_values)

plt.legend(['Experimental C++ Data','Fitted Extrapolated Equation'], loc='upper left')

plt.show()

#Printing out the Tinfinite value
print(np.polyval(fit, 0))
</pre>
|}
</center>