View source for Talk:MIKE123

'''<span style="color:red">JC: General comments: All tasks answered and data looks good.  Try to make your explanations of your results clearer and more focused and bring in some of the background theory where appropriate to show your understanding.</span>'''

Transition States and Reactivity

=Introduction to the Icing Model=

==The Model==

===Task 1: Lowest Energy of the Icing model===

The interaction energy is defined as <math>- \frac{1}{2} J \sum_i^N \sum_{j\  \in\  \mathrm{neighbours}\left(i\right)} s_i s_j</math>, where the first sum loops over the number of spins in the system and the second sum loops over the number of neighbours. When keeping all parameters constant apart from the number of dimensions, the only thing changing is the second sum. 

For <math>D =1</math> the number of neighbours is <math>2</math>, for <math>D=2</math> it's <math>4</math> and for <math>D=3</math> it's <math>6</math>. In general one can say that it the number of neighbours is <math>2D</math>. 

In the lowest energy state all spins will be aligned leading to maximum multiplicity and thus a scenario where <math>s_i = s_j</math> for all <math>i</math> and<math>j</math>. As the value of <math>s</math> can only have the values <math>1</math> or <math>-1</math>, the product under the sum must always equal to <math>1</math>.

This interaction energy equation therefore simplifies to <math>- \frac{1}{2} J \sum_1^N \sum_1^{2D} 1</math>. The <math>\frac{1}{2}</math> cancels out with the <math>2</math> from <math>2D</math> and once the sums are evaluated they simplify to <math>E\ =\ -DNJ</math>. q.e.d.

The entropy of this system is given by <math>S = k_B \ln W</math>, where <math>W</math> is the multiplicity. In this system the multiplicity is one as there is only one configuration for a particular state when all spins are aligned. Thus the entropy equals to zero.

'''<span style="color:red">JC: Correct. Did you consider that the ground state can be all spin up or all spin down?</span>'''

==Phase Transitions==

===Task 2: Energy change caused by a change in state and the entropy involved===

Using the equation for the interaction energy for the lowest state from above, the lowest energy configuration has a value of <math>E=-3000 J</math>.
In order to obtain an energy for the system when one spin is antiparallel to the others one can examine the general equation. The antiparallel spin will occur as <math>s_i</math> once and as <math>s_j</math> 6 times (for 6 different <math>s_i</math> neighbours once). This will lead to terms adding to the energy and terms cancelling out each other eventually leading to a value of <math>E = -2988 J</math>. So the energy change is <math>dE = 12 J</math>.

'''<span style="color:red">JC: Can you show more of the maths working for this step.</span>'''

The flipping of the spin can happen on any point of the 1000 sites in the lattice. This demonstrates that there are 1000 configurations where one spin is antiparallel to the others. This leads to an entropy value gain of <math>dS=k_B \ln 1000=9.53 \times 10^{-23}  J  K^{-1}</math>

'''<span style="color:red">JC: Correct.</span>'''

===Task 3: Magnetisation in a 3-D lattice===

As there is always one more positive spin in both the 1D and 2D lattice, both systems have a magnetisation of <math>+1</math>.

Absolute zero indicates that all spins in a system are aligned, so they either all have an value of <math>+1</math> or <math>-1</math>, leading to a total magnetisation of either <math>+1000</math> or <math>-1000</math> respectively.

=Calculating the energy and magnetisation=

===Task 4: completing the energy() and magnetisation() function ===
Please refer to the comments in the code:

<pre>
    def energy(self):        
        energy = 0.0
#copying the boundaries of the lattice and extending/appending them on the opposite ends in all four directions
        new_lattice = np.append(self.lattice, [self.lattice[0]], axis=0)
        new_lattice2 = np.append([self.lattice[-1]], new_lattice, axis = 0)
        new_lattice3 = np.concatenate((np.transpose([new_lattice2[:,-1]]),new_lattice2), axis = 1)
        new_lattice4 = np.concatenate((new_lattice3, np.transpose([new_lattice3[:,1]])), axis = 1)       
#looping through the new lattice 
        for ycor in range(len(new_lattice4)):
            for xcor in range(len(new_lattice4)):
#this makes sure that only the points from the original lattice are counted
                if ycor >0 and ycor < len(new_lattice4)-1 and xcor >0 and xcor < len(new_lattice4[0])-1:
#checks the energy in every direction and adds the contribution
                    if new_lattice4[ycor][xcor] + new_lattice4[ycor-1][xcor] == 0:
                        energy += +0.5
                    else:
                        energy += -0.5
                    if new_lattice4[ycor][xcor] + new_lattice4[ycor+1][xcor] == 0:
                        energy += +0.5
                    else:
                        energy += -0.5
                    if new_lattice4[ycor][xcor] + new_lattice4[ycor][xcor-1] == 0:
                        energy += +0.5
                    else:
                        energy += -0.5
                    if new_lattice4[ycor][xcor] + new_lattice4[ycor][xcor+1] == 0:
                        energy += +0.5
                    else:
                        energy += -0.5                
        return energy

    def magnetisation(self):
#sum() was assumed common knowledge
        magnetisation = np.sum(self.lattice)
        return magnetisation
</pre>

'''<span style="color:red">JC: Code looks good and is well commented.</span>'''

===Task 5: Checking the code using the ILcheck.py script===
The values seemed to match!

[[File:Task4EnergyandMagnetisationMIKE123.png]]

=Introduction to the Monte Carlo simulation=

==Average Energies and Magnetisation==

===Task 6: Configuration of 100 spins and the estimated time it takes to evaluate the average magnetisation at one temperature ===

In a system of 100 spins there are <math>2^{100}</math> configurations. In base 10 this approximates to about <math>1.27 * 10^{30}</math> configurations. So dividing this number by the computational rate of <math>10^9</math> configurations per seconds gives us a time of <math>1.28 * 10^{21}</math> seconds or about <math>4 * 10^{13}</math> years.

'''<span style="color:red">JC: So evaluating every configuration is not possible.</span>'''

==Modifying IsingLattice.py==

===Task 7: Completing the montecarlostep function and the statistics function to obtain values for <math><E></math> and <math><M></math>===

Please refer to the comments in the code below:
<pre>
    def montecarlostep(self, T):
#saving the energy before the change
        energy0 = self.energy()
        magnet0 = self.magnetisation()
#the following two lines will select the coordinates of the random spin for you
        random_i = np.random.choice(range(0, self.n_rows))
        random_j = np.random.choice(range(0, self.n_cols))
#flipping the spins in the original lattice
        if self.lattice[random_i][random_j] == -1:
            self.lattice[random_i][random_j] += 2
        else:
            self.lattice[random_i][random_j] -= 2
#saving the energy after the change
        energy1 = self.energy()
        magnet1 = self.magnetisation()
        energy_d = energy1 - energy0
#the following line will choose a random number in the range [0,1) for you
        random_number = np.random.random()
#checking whether the original lattice should have been kept (rather than checking if we should change)
        if energy_d > 0 and random_number > np.e**(-energy_d/T):
#if yes then the spin is flipped back and the original energy is stored in energy1
            if self.lattice[random_i][random_j] == -1:
                self.lattice[random_i][random_j] += 2
            else:
                self.lattice[random_i][random_j] -= 2
            energy1 = energy0
            magnet1 = magnet0
#always append the energy1 from the new system/corrected-back system
        self.E += energy1
        self.E2 += energy1**2
        self.M += magnet1
        self.M2 += magnet1**2
        self.n_cycles += 1
        return energy1, magnet1

    def statistics(self):
#divide by cycles, have an if condition so that it doesn't crash when n_cycles = 0
        if self.n_cycles > 0:
            avg_e = self.E/self.n_cycles
            avg_e2 = self.E2/self.n_cycles
            avg_m = self.M/self.n_cycles
            avg_m2 = self.M2/self.n_cycles
            return avg_e, avg_e2, avg_m, avg_m2, self.n_cycles
        else:
            return self.E, self.E2, self.M, self.M2, self.n_cycles </pre>

'''<span style="color:red">JC: Looks correct, to flip the spin you could just multiply the spin by -1.</span>'''

===Task 8: Discussing the the behaviour of a system when <math>T < T_C</math>===

As given in the instructions, the system under investigation is set to be below the <math>T_C</math>, so that <math>T < T_C</math> is fulfilled. In this case thermal energy is too low to influence any entropic forces, so that it is very unlikely that higher spin energy states are populated spontaneously. Thus what will happen is that the system tries to minimise it's energy by aligning all spins. This will lead to all spins being parallel and thus <math><M> \neq 0</math>. This was also observed when running ILanim.py on the system:

[[File:CMPWIKIPICTUREOFTHEENERGYSTEPS.png | 360px]]

{| class="wikitable" style="text-align: left;"
|-
! style="width:20em;" |Quantity after 445 steps
! style="width:20em;" |Value after 445 steps
|-
|<math><E></math>
| -1.329485
|-
|<math><E^2></math>
| 2.120480
|-
|<math><M></math>
| 0.664477
|-
|<math><M^2></math>
| 0.536613
|-
|}

=Accelerating the code=

===Task 9: Duration of ILtimetrial.py for 2000  Monte Carlo steps===
Duration for 2000 Monte Carlo Steps:
{| class="wikitable" style="text-align: left;"
|-
! style="width:3em;" |Trial
! style="width:10em;" |Time (s)
|-
|1
|17.46
|-
|2
|17.55
|-
|3
|17.52
|-
|}
So one can safely assume the code executes within about 17.5 +/- 0.2 seconds.

===Task 10: Improving the energy() and magnetisation() function===
<pre>
    def energy(self):       
        energy = 0.0
        s_ilattice = self.lattice
#creates four lattices which are shifted with respect to s_ilattice in the four directions by one item (they wrap around edges)
        s_jlattice1 = np.roll(s_ilattice, 1, axis=0) #down
        s_jlattice2 = np.roll(s_ilattice, -1, axis=0) #up
        s_jlattice3 = np.roll(s_ilattice, 1, axis=1) #right
        s_jlattice4 = np.roll(s_ilattice, -1, axis=1) #left
#they are added essentially creating a latticesum of neighbours for every item in s_ilattice
        s_jlatticesum = s_jlattice1 + s_jlattice2 + s_jlattice3 + s_jlattice4
#the coordinates of an item in s_ilattice correspond to the coordinates of the sum of its neighbours in s_jlatticesum
#multiplying them with each other gives s_i *s_j for all neighbours
#taking the sum in order to get a total energy throughout the lattice
        energy += np.sum(np.multiply(s_ilattice, s_jlatticesum))*(-0.5)
        return energy

#as before
    def magnetisation(self):
        magnetisation = np.sum(self.lattice)
        return magnetisation
</pre>

'''<span style="color:red">JC: You actually only need to use roll twice (if you remove the factor of 0.5 from the energy) as this prevents double counting of interactions.</span>'''

===Task 11: Duration of ILtimetrial.py for 2000  Monte Carlo steps after function improvement ===
{| class="wikitable" style="text-align: left;"
|-
! style="width:3em;" |Trial
! style="width:10em;" |Time before improvement (s)
! style="width:10em;" |Time after improvement (s)
|-
|1
|17.46
|0.42
|-
|2
|17.55
|0.43
|-
|3
|17.52
|0.43
|-
|}

Thus an average of about 0.43 +/- 0.1 seconds for the 2000 Monte Carlo Steps can be assumed. This is a major improvement of about 17 seconds. A major contributor to this improvement in run time is certainly the removal of for loops and the implementation of built in functions.

=The effect of temperature=

==Correcting the Averaging Code==
=== Task 12: Adjusting the code to account only for fluctuations around the equilibrium state ===
<pre>

# global variable has been implemented which can be altered at wish
    ignored_cyc = 225000

    def montecarlostep(self, T):
        energy0 = self.energy()
        magnet0 = self.magnetisation()
        random_i = np.random.choice(range(0, self.n_rows))
        random_j = np.random.choice(range(0, self.n_cols))
        if self.lattice[random_i][random_j] == -1:
            self.lattice[random_i][random_j] += 2
        else:
            self.lattice[random_i][random_j] -= 2
        energy1 = self.energy()
        magnet1 = self.magnetisation()
        energy_d = energy1 - energy0
        random_number = np.random.random()
        if energy_d > 0 and random_number > np.e**(-energy_d/T):
            if self.lattice[random_i][random_j] == -1:
                self.lattice[random_i][random_j] += 2
            else:
                self.lattice[random_i][random_j] -= 2
            energy1 = energy0
            magnet1 = magnet0
#energies are only added if the number of cycles is larger than the number of ignored cycles
        if self.n_cycles > self.ignored_cyc:
            self.E += energy1
            self.E2 += energy1**2
            self.M += magnet1
            self.M2 += magnet1**2
        self.n_cycles += 1
        return energy1, magnet1
        
    def statistics(self):
        if self.n_cycles > 0:
#statistical values are calculated with the corrected number of cycles
            avg_cycles = self.n_cycles - self.ignored_cyc
            avg_e = self.E/avg_cycles
            avg_e2 = self.E2/avg_cycles
            avg_m = self.M/avg_cycles
            avg_m2 = self.M2/avg_cycles
            return avg_e, avg_e2, avg_m, avg_m2
        else:
            return self.E, self.E2, self.M, self.M2
</pre>

The system was analysed for a grid of <math>16 \times 16</math> and  <math>32 \times 32</math> at temperatures of  <math>T = 2.0</math> and  <math>T = 4.0</math>.
It was discovered that as the critical temperature lies somewhere between those two temperatures, as below  <math>T_C</math> the system took time (up to 225000 steps) to reach equilibrium while when  <math>T > T_C</math> the system had enough thermal energy to easily overcome the energy barrier and thus reach equilibrium really fast (less than 1000 steps) . The latter observation is due to the fact that a above the Curie Temperature random magnetisation is very unlikely to happen. 
It was further observed that at bigger lattice sizes the time it took for the system to reach equilibrium was larger.
It was thus concluded that in order to be prepared for the worst case scenario (<math>32 \times 32</math> with <math>T</math> just under <math>T_C</math>), a total of 225000 cycles will be ignored. This number can be obtained from observing the last graph.

<gallery widths=280px heights=200px>
File:MonteCarlo0ignored2temp150000.png|16x16 at T = 2.0, 150000 steps
File:MonteCarlo0ignored4temp150000.png|16x16 at T = 4.0, 150000 steps
File:MonteCarlo32320ignored2temp500000.png|32x32 at T = 2.0, 500000 steps
File:MonteCarlo32320ignored4temp500000.png|32x32 at T = 4.0, 500000 steps
</gallery>

'''<span style="color:red">JC: Good analysis and correctly implemented in code.</span>'''

==Running over a range of temperatures==

===Task 13: Average Energy and Magnetisation for each Temperature===

It was decided that the simulations will run in the temperature range of<math> T = 0.5 </math>and <math>T=5.0</math> in steps of <math>dT = 0.05</math> for a total number of 350000 steps with 225000 steps being disregarded in the energy and magnetisation averages. The results are presented in the next section.

=The effect of system size=
==Scaling the System Size==
===Task 14: Energy per spin versus temperature diagrams for different lattice sizes ===

In order to plot the functions from the obtained data files the following function was defined to be called for every lattice size and the output is shown below.

<pre>
import numpy as np
from matplotlib import pylab as plt


def plotfunct(filename, lattice_side_size):
    data = np.loadtxt(filename)
    temps = data[:,0]
    energies=data[:,1]
    energies2=data[:,2]
    magnet=data[:,3]
    magnet2=data[:,4]
    spins = lattice_side_size**2
#errorbars as square root of variance, energies are given per spin
    energies_error = np.sqrt(energies2-energies**2)/spins
    magnet_error = np.sqrt(magnet2-magnet**2)/spins

    plt.ylabel("Energy per spin")
    plt.xlabel("Temperature")
    plt.title('Energy per spin vs Temperature (' + str(lattice_side_size)+'x'+str(lattice_side_size)+')')
    plt.errorbar(temps, energies/spins,yerr = energies_error, ecolor = "green")
    plt.show()
    plt.ylabel("Magnetisation per spin")
    plt.xlabel("Temperature")
    plt.title('Magnetisation per spin vs Temperature (' + str(lattice_side_size)+'x'+str(lattice_side_size)+')')
    plt.errorbar(temps, magnet/spins, yerr = magnet_error, ecolor = "green")
    pl.show()
</pre>

As lattice size increases, fluctuations and error bars decrease massively. It is further observed that the point at which equilibrium is reached in the  magnetisation graph shifts from left to right (to increasing temperature) as lattice size is increasing. The first observation can be explained through the way data is processed. In small systems like 2x2, flipping of one spin has a huge impact on the absolute change in energy (1 in 4 spins), whereas the same flip in a 32x32 system has a minor contribution (1 in 1024). Latter will thus have a much smaller impact on the statistical means that are used to obtain the graphs. The second observation can also be explained through the quantity of spins in each lattice. As the Monte Carlo algorithm assigns multiple random numbers when it is called, it is much more likely to totally reverse a small system than a large one. I.e. completely reversing a system is much more likely to happen in a 2x2 lattice (4 different flips needed) versus a 32x32 lattice (1024 different flips needed!!!).

The quality of the data (reduced noise) does not seem to linearly depend on lattice size (judged by eye). However it does definitely increase with increasing lattice size. While 32x32 seems adequate to perform calculations with, 64x64 would very likely improve the data even more.

<gallery widths=320px heights=300px>
File:Xxx2x2E.png
File:Xxx4x4E.png
File:Xxx8x8E.png
File:Xxx16x16E.png
File:Xxx32x32E.png
</gallery>
<gallery widths=320 heights=300px>
File:Xxx2x2M.png
File:Xxx4x4M.png
File:Xxx8x8M.png
File:Xxx16x16M.png
File:Xxx32x32M.png
</gallery>

'''<span style="color:red">JC: Graphs look good, did you also plot the energy per spin for each system size on the same axes to compare them?</span>'''

=Determining the heat capacity=
==Calculating the heat Capacity==
===Task 15: Deriving an expression for heat capacity from the partial derivative expression===

1. The heat capacity is given as <math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math> , with 
<math>\left\langle E\right\rangle = \frac{\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_\mathrm BT}\right\}}{\sum_\alpha  \exp \left\{-\frac{E\left(\alpha\right)}{k_\mathrm BT}\right\}}</math> 

2. Defining <math>Z=\sum_\alpha  \exp \left\{-\frac{E\left(\alpha\right)}{k_\mathrm BT}\right\}</math> and <math>y= \frac{1}{k_\mathrm BT}</math> for simplification purposes and noting that <math>\frac{\partial y}{\partial T} = -\frac{1}{k_BT^2}</math>

3. The expression which should be derived is <math>C = \frac{\mathrm{Var}[E]}{k_B T^2}</math>, where <math>\mathrm{Var}[E] = \left\langle E^2\right\rangle - \left\langle E\right\rangle ^2 </math>. 

4. So from that an expression for <math>\left\langle E\right\rangle </math> and <math>\left\langle E^2\right\rangle </math> in terms of <math>Z</math> and <math>y</math> shall be found.

5. Differentiating <math> Z </math> with respect to <math>y</math> gives:
<math>\frac{\partial Z}{\partial y}= - \sum_\alpha E\left(\alpha\right) \exp \left\{-yE\left(\alpha\right)\right\} </math>.

6. Using this information <math>\left\langle E\right\rangle </math> can be rewritten as  <math>\left\langle E\right\rangle = -\frac{1}{Z}\frac{\partial Z}{\partial y}</math>

7. Differentiating <math> Z </math> with respect to <math>y</math> twice gives:
<math>\frac{\partial^2 Z}{\partial y^2} = \sum_\alpha E\left(\alpha\right)^2 \exp \left\{-\frac{E\left(\alpha\right)}{k_\mathrm BT}\right\} </math>

8. Thus <math>\left\langle E^2\right\rangle </math> can be rewritten as <math>\left\langle E^2\right\rangle = \frac{1}{Z} \frac{\partial^2 Z}{\partial y^2}  </math>

9. Evaluating the original expression using the product rule gives:
<math>C  = \frac{\partial}{\partial T} \left(- \frac {1}{Z} \frac{\partial Z}{\partial y}\right) =  \frac{\partial \left(-\frac{1}{Z}\right)}{\partial T} \frac{\partial Z}{\partial y} - \frac{1}{Z} \frac{\partial}{\partial T} \left( \frac{\partial Z}{\partial y}\right) </math>

10. Applying the chain rule twice in the first term and once in the second term gives:
<math> C =  \frac{\partial \left(-\frac{1}{Z}\right)}{\partial Z} \frac{\partial Z}{\partial y} \frac{\partial y}{\partial T}\frac{\partial Z}{\partial y} - \frac{1}{Z} \frac{\partial \left( \frac{\partial Z}{\partial y}\right)}{\partial y}\frac{\partial y}{\partial T} =\frac{\partial y}{\partial T} \left( \left(\frac{1}{Z}\frac{\partial Z}{\partial y}\right)^2-\frac{1}{Z}\frac{\partial ^2Z}{\partial y^2}\right)</math>

11. Substituting the expressions for the partial derivatives from 2., 6. and 8., hence gives:
<math>C = \frac{\left\langle E^2\right\rangle - \left\langle E\right\rangle ^2} {k_B T^2}</math>

'''<span style="color:red">JC: Derivation is correct, but it might be easier to follow if you can just start from one side of the equation and show the steps needed to reach the other side.</span>'''

===Task 16: Plotting Heat Capacity vs Temperature for different lattice sizes===
<pre>
import numpy as np
from matplotlib import pylab as plt


def plotheatcap(filename, lattice_side_size):
    data = np.loadtxt(filename)
    temps = data[:,0]
    spins = lattice_side_size**2
    energies=data[:,1]
    energies2=data[:,2]
    varian = (energies2-energies**2)/spins

    plt.ylabel("Heat Capacity")
    plt.xlabel("Temperature")
    plt.plot(temps, varian/temps**2, label = str(lattice_side_size) + "x"+ str(lattice_side_size))
    

def plotmultiple(filenamearray, lattice_side_array):
    for i in range(len(filenamearray)):
        plotheatcap(filenamearray[i], lattice_side_array[i])
    plt.title("Heat capactiy per spin vs temperature for various lattice sizes")
    plt.legend(loc="best",title="Lattice size",fontsize="8")
    plt.show()
    
plotmultiple(['2x2.dat','4x4.dat','8x8.dat','16x16.dat','32x32.dat'],[2,4,8,16,32])

</pre>

This code leads to the graph below. The heat capacity per spin seem to slightly increase with the size of the lattice. The peak does also seem to become more narrow. This directly relates to the previous discussion about lattice size affecting statistical measures, as heat capacity is directly related to the variance of the obtained energies. I.e. larger spins system lead to smaller error / variance and thus a more narrow peak of the heat capacity. Ideally at very large lattice size this peak will become as narrow as possible.

'''<span style="color:red">JC: For an infinite lattice the peak diverges.</span>'''
 
<gallery widths = 500px heights=380px>
File:HeatCapxxx123.png
</gallery>

=Locating the Curie temperature=
==C++==
===Task 17: Comparing the quality of data obtained using python versus C++===
<pre>
import numpy as np
from matplotlib import pylab as plt


def plottwocap(filenamearray, lattice_side_array):
    plotheatcap(filenamearray[0],lattice_side_array[0])

    data2 = np.loadtxt(filenamearray[1])
    temps2 = data2[:,0]
    spins2 = lattice_side_array[1]**2
    energies2=data2[:,1]
    energies22=data2[:,2]
    varian2 = (energies22-energies2**2)
    capa = data2[:,5]
    
    plt.ylabel("Heat Capacity")
    plt.xlabel("Temperature")
    plt.plot(temps2, capa, label = "C++")
    plt.title("Heat capactiy per spin vs temperature (32x32)")
    plt.legend(loc="best",title="Lattice size",fontsize="8")
    plt.show()

def plottwo2(filenamearray, lattice_side_array):
    data1 = np.loadtxt(filenamearray[0])
    temps1 = data1[:,0]
    spins1 = lattice_side_array[0]**2
    energies = data1[:,1]
    magnet = data1[:,3]

    data2 = np.loadtxt(filenamearray[1])
    temps2 = data2[:,0]
    energies2=data2[:,1]
    magnet2 = data2[:,3]
    
    plt.ylabel("Energy per spin")
    plt.xlabel("Temperature")
    plt.title('Energy per spin vs Temperature (' + str(lattice_side_array[0])+'x'+str(lattice_side_array[0])+')')
    plt.plot(temps1, energies/spins1,label="python")
    plt.plot(temps2, energies2,label="C++")
    plt.legend(loc="best",title="Language",fontsize="8")
    plt.show()
    
    plt.ylabel("Magnetisation per spin")
    plt.xlabel("Temperature")
    plt.title('Magnetisation per spin vs Temperature (' + str(lattice_side_array[0])+'x'+str(lattice_side_array[0])+')')
    plt.plot(temps1, magnet/spins1,label="python")
    plt.plot(temps2, magnet2,label="C++")
    plt.legend(loc="best",title="Language",fontsize="8")
    plt.show()
    
plottwo2(['32x32.dat','C++_data/32x32C.dat'],[32,32])

plottwocap(['32x32.dat','C++_data/32x32C.dat'],[32,32])
</pre>

The data from the C++ files was collected in steps of <math>dT = 0.01T </math> rather than <math>dT = 0.05T </math> as was done in the previous section. This lead to a higher data rate and higher resolution in all graphs. It was especially noticeable in the region around the Curie Temperature for the magnetisation per spin graph as fluctuations were made visible which were originally assumed completely absent.

<gallery widths=300px heights=250px>
File:ENEEE.png
File:MAGMM.png
File:HEATTT.png
</gallery>

==Polynomial fitting==
===Task 18: Polyfitting Heat Capacity to a polynomial===
<pre>
import numpy as np
from matplotlib import pylab as plt

def plotpolyfit(filename):
    dataA = np.loadtxt(filename)
    temps = dataA[:,0]
    capa = dataA[:,5]
    
    fit = np.polyfit(temps,capa,5)
    fit2 = np.polyfit(temps,capa,100)
    T_min=np.min(temps)
    T_max=np.max(temps)
    T_range=np.linspace(T_min,T_max,1000)
    fitted_C_values = np.polyval(fit,T_range)
    fitted_C_values2 = np.polyval(fit2,T_range)

    plt.plot(temps,capa,label='Collected Data')
    plt.plot(T_range,fitted_C_values,label='Fitted Data (Degree 5)')
    plt.plot(T_range,fitted_C_values2,label='Fitted Data (Degree 100)')
    plt.legend(loc ="best", fontsize= 8)
    plt.title('Heat Capacity per Spin Curve Fitting (32x32)')
    plt.xlabel('Temperature')
    plt.ylabel("Heat Capacity")
    plt.show()

plotpolyfit('C++_data/32x32C.dat')
</pre>
Some initial trials to polyfit a polynomial to the whole data range only showed little success. Even when increasing the polynomial degree to 100, it was not possible to approximate the peak accurately.
<gallery widths = 500px heights=450px>
File:Adadasdasdasdc.png
</gallery>

===Task 19: Polyfit the heat capacity to a polynomial in a specific temperature range===
<pre>
import numpy as np
from matplotlib import pylab as plt

def plotpolyfitinrange(filename, dimension, minT, maxT, degree):
    dataA = np.loadtxt(filename)
    temps = dataA[:,0]
    capa = dataA[:,5]

    selection = np.logical_and(temps > minT, temps < maxT) #choose only those rows where both conditions are true
    peak_T_values = temps[selection]
    peak_C_values = capa[selection]
    
    fit = np.polyfit(peak_T_values,peak_C_values,degree)
    T_min=np.min(temps)
    T_max=np.max(temps)
    T_range=np.linspace(T_min,T_max,1000)
    fitted_C_values = np.polyval(fit,T_range)
    fitted_C_values2 = np.polyval(fit,peak_T_values)
    Cmax = np.max(fitted_C_values2)
    Tmax = peak_T_values[fitted_C_values2 == Cmax]
    print(Cmax, Tmax)
    plt.plot(temps,capa,label='Collected Data')
    plt.plot(T_range,fitted_C_values,label='Fitted Data between ['+ str(minT) +','+str(maxT)+'] (Degree '+str(degree)+')')
    plt.legend(loc ="best", fontsize= 8)
    plt.title('Heat Capacity per Spin Curve Fitting (' + str(dimension)+'x'+str(dimension)+')')
    plt.xlabel('Temperature')
    plt.ylabel("Heat Capacity")
    plt.ylim([-0.5,2])
    plt.show()

plotpolyfitinrange('C++_data/2x2C.dat',2,1.5,3.4,4)
plotpolyfitinrange('C++_data/4x4C.dat',4,1.5,3.4,4)
plotpolyfitinrange('C++_data/8x8C.dat',8,1.7,3,6)
plotpolyfitinrange('C++_data/16x16C.dat',16,2,2.7,6)
plotpolyfitinrange('C++_data/32x32C.dat',32,2.2,2.4,5)
</pre>

A big improvement was made when small regions were selected to run through the polyfit function. The last parameter in the function "plotpolyfztinrange" above indicates the polynomial degree used for further calculations. By reducing the domain of the function, even polynomial of degree 4 were enough to accurately fit the heat capacity data (for 2z2 and 4x4). 

<gallery widths = 500px heights = 400px>
File:POLYFIT2njsdnfksndk.png
</gallery>

'''<span style="color:red">JC: Fit looks good.</span>'''

===Task 20: Plotting the scaling relation between the Curie Temperature and the lattice side length===

The data for the table below was obtained when executing the previous function for every lattice. It showed that the Curie Temperature decreased with increasing lattice size as was expected in the provided function of <math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math>. 

{| class="wikitable"
! Size!! Temperature
|-
| 2x2|| 2.52
|-
| 4x4|| 2.48
|-
| 8x8|| 2.36
|-
| 16x16|| 2.32
|-
| 32x32|| 2.30
|}

In order to obtain the Currier temperature at infinite lattice the temperatures (table above) at which the heat capacities were maximised were plotted against the inverse of the respective lattice side length, <math>\frac{1}{L}</math>.

<pre>
from matplotlib import pylab as plt
import numpy as np

sidelength = np.array([2,4,8,16,32])
curieT = np.array([2.52,2.48,2.36,2.32,2.30])

fit = np.polyfit(1/sidelength, curieT, 1, cov=True)
params = fit[0]
params_cov = fit[1] 

x_range = np.linspace(0, 0.5, 1000) 
fitted_curie = params[0]*x_range + params[1]
plt.xlabel("1/L")
plt.ylabel("Curie temperature")
plt.plot(1/sidelength, curieT, linestyle='',marker='x',label="Collected Data")
plt.plot(x_range, fitted_curie, label="Linear fit")
plt.title("Curie temperature versus 1/L for LxL lattices")
plt.legend(loc="best",fontsize=10)
plt.show()
</pre>

<gallery widths =500px heights= 400px>
File:Sclaingrelation123123.png
</gallery>


A value of <math> 2.301 \pm 0.044</math> for the Currier Temperature at infinite lattice was obtained (y-intercept). This result lies within the boundaries of the literature value of 2.269.<ref name="Gangsta"/>. However, when analysing the graph below it seems quite obvious that the point at 2x2 is an outlier. The reasons for that were discussed thoroughly in previous sections above as the uncertainty of the 2x2 lattice's data continuously showed up throughout this report. 
Additional minor factors that contribute to this deviation from the literature value will come from the method used, i.e. the Monte Carlo algorithm rejects complete randomness and is thus an approximation. However, it is actually rather a personal surprise that such a good close result was obtained using this algorithm, that cut down run time from about <math>10^{13}</math> years to just around 30 minutes. So essentially having a fairly big lattice size of about <math>128 \times 128</math> should be ideal to increase accuracy while sacrificing little additional computational time. 

'''<span style="color:red">JC: Considering the problems with the 2x2 system, you could remove this data point from the fit and see if it improves the agreement with the theoretical value from the literature.</span>'''

==References==
<references>
<ref name="Gangsta">L. Onsager, Crystal Statistics. I. A Two-Dimensional Model with an Order-Disorder Transition, <i>Physical Review</i>, 1944, <b>65</b>, 117-149.</ref>
</references>