View source for Rep:ModMonteCarlo:K98631

'''<u>Peter Wright 3<sup>rd</sup> year CMP python- Monte Carlo simulation report</u>'''

<u>'''Introduction to the Ising model'''</u>

a) In the Ising model we consider a lattice, of variable dimension D, with N spins arranged on individual lattice domain sites, with spins s<sub>i,j</sub>=+-1.For a configuration ‘v’, the total interaction energy due to spin interaction, J can be written using the following equation;

[[File:PWX4.png |thumb|centre| Total energy of Ising Lattice when no external magnetic field is present.]]

J is a spin-coupling constant, and interactions occur between neighboring spins. Periodic boundary conditions are set such that the
lattice site at the end of a row/column interacts with the lattice site at the opposite end. This happens in both directions. Parralel spins have a –J
interaction and opposing spins have no energetic interaction.

Because J>0, spin alignment is energetically
favorable, and hence spin alignment is favoured. Below the Curie temperature;
where the energetic component of the systems free energy is dominant, we expect ''spontaneous magnetization, ''leading to
long-range ordering of spins. The lowest energy ground state of the system has
E=-DNJ.

This is shown by considering a generous square 2D
lattice; the total number of edges= 4N/2, or 2N if there are N sites. As the
ground equilibrium state of the lattice correspond to either all spins being
+1(up) or -1(down), i.e in the same state, the ground state energy will be –J
for each edge. Every spin, in an nxn lattice of dimensionality D, is connected
to its neighbor by 2D (four) edges, so the total contribution, taking into account
periodic boundary conditions is&nbsp; -4NJ/2.
The factor of 0.5 is to account for double counting that occurs when cycling
through each lattice site as each edge appears twice in the sum. Therefore for
dimensionality D, the ground state energy= -DNJ.

b) For an nxn lattice consisting of N=1000 spins, and
D=3. If we consider the lattice to be in its ground state, whereby E=-DNJ, and
we now randomly flip one spin such that it is now aligned oppositely to all
other lattice spins, the system configuration will increase in energy.

As D=3, each lattice site has 6 edges, therefore when
a single spin is flipped the total number of energetically favorable
interactions lost is +6J. This does not propagate as the system wants to
minimize the energy loss, and it would be rare, due to Boltzmann probability
that this new configuration would be adopted. Nevertheless, for the lattice in
question; &nbsp;= -2994J - -3000J = +6J. This spin flipping
leads to greater system disorder, which is entropically favorable. The ground
state configuration, which is favoured at T<T<sub>C</sub>, has only one ‘microstate
(w)’ or permutation, and therefore the third law of thermodynamics; S=K<sub>B</sub>In(0) &nbsp;= 0. The new, higher energy configuration has
N lattice permutations, and therefore S=K<sub>B</sub>In(1000) ≈ 9.54x10<sup>-23</sup>JK<sup>-1</sup>mol<sup>-1</sup>.

This was also shown using calculations made in python-
source code found below for 1, 2 and 3D lattices. '''(EXTRA CODE, JUST FOR DEMONSTRATION)'''

<nowiki><code></nowiki>

import numpy as np

f = lambda x, y : 1 if np.sign(x) == np.sign(y) else 0

J = -1;

D = 3; N = 3; dim = np.array([N] * D); e_min=0; a_min = np.zeros(N ** D); a_min = np.reshape(a_min, dim)


for k in range(1):

&nbsp; &nbsp; a = np.random.choice([-1,1],N ** D); a = np.reshape(a, dim)

&nbsp; &nbsp; a = np.ones(N ** D); a[0] = -1

&nbsp; &nbsp; a = np.reshape(a, dim)

&nbsp; &nbsp; if D == 1: b = a; d = [[f(j, b[i-1]) + f(j, b[(i+1) % N]) for i,j in enumerate(b)]]

&nbsp; &nbsp; if D == 2:

&nbsp; &nbsp; &nbsp; &nbsp; d1 = [[f(j, b[i-1]) + f(j, b[(i+1) % N]) for i,j in enumerate(b)] for b in a]&nbsp;

&nbsp; &nbsp; &nbsp; &nbsp; d2 = [[f(j, b[i-1]) + f(j, b[(i+1) % N]) for i,j in enumerate(b)] for b in a.T]

&nbsp; &nbsp; &nbsp; &nbsp; d = np.sum(d1) + np.sum(d2)

&nbsp; &nbsp; if D == 3:

&nbsp; &nbsp; &nbsp; &nbsp; c = a

&nbsp; &nbsp; &nbsp; &nbsp; d1 = [[[f(j, b[i-1]) + f(j, b[(i+1) % N]) for i,j in enumerate(b)] for b in a] for a in c]

&nbsp; &nbsp; &nbsp; &nbsp; d2 = [[[f(j, b[i-1]) + f(j, b[(i+1) % N]) for i,j in enumerate(b)] for b in a.T] for a in c]

&nbsp; &nbsp; &nbsp; &nbsp; d3 = [[[f(j, b[i-1]) + f(j, b[(i+1) % N]) for i,j in enumerate(b)] for b in a] for a in c.T]

&nbsp; &nbsp; &nbsp; &nbsp; d = np.sum(d1) + np.sum(d2) + np.sum(d3)

&nbsp; &nbsp; e = 0.5 * J * np.sum(d)

&nbsp; &nbsp; if e < e_min: e_min = e; a_min = a

print(e_min, a_min)

<nowiki></code></nowiki>

It is shown that a single spin flip, leads to an
energy change of '''+2DJ.'''

c) A lattices total magnetization is calculated via a simple sum of all lattice site spins;

Figure 1 PNG (see CMP Intro to Ising Lattice page1).

1D- 3up spins and 2 down spins = +1

2D- 13up spins and 12 down spins= +1

At absolute zero, T-0K, thermal energy available to a
lattice, K<sub>B</sub>T=0J. Considering this, the lattice cannot flip a spin in
order to reach a higher energy configuration other than its ground state. Furthermore
there will be no entropic contribution to the free energy, making the process
impossible. Therefore a lattice, D=3 and N=1000; all spins will be in the same
state.

=+1000 or -1000.

<u>'''Calculating the energy and
magnetisation'''</u>

Working in reduced units,
where J=K<sub>B</sub> python IsingLattice functions are completed. Source code found below.

<nowiki><code></nowiki>

'''def energy(self):'''

&nbsp; &nbsp; &nbsp; &nbsp; "Return the total energy of the current lattice configuration."

&nbsp; &nbsp; &nbsp; &nbsp; f = lambda x, y : 1 if np.sign(x) == np.sign(y) else 0 #lamdba function to compute energy of interaction between neighbours x and y

&nbsp; &nbsp; &nbsp; &nbsp;J = 1

&nbsp; &nbsp; &nbsp; &nbsp; d1 = [[f(j, b[i-1]) + f(j, b[(i+1) % self.n_rows]) for i,j in enumerate(b)] for b in self.lattice]&nbsp;#enumerate rows and calculate energy with modulus calc on index to satisfy the periodic boundary condition

&nbsp; &nbsp; &nbsp; &nbsp; d2 = [[f(j, b[i-1]) + f(j, b[(i+1) % self.n_cols]) for i,j in enumerate(b)] for b in self.lattice.T] #enumerate columns (transpose) and calculate energy

&nbsp; &nbsp; &nbsp; &nbsp; d = np.sum(d1) + np.sum(d2) #add row and column energies

&nbsp; &nbsp; &nbsp; &nbsp; energy = 0.5 * J * np.sum(d) #0.5 factor adjusts for double counting in neighbour interactions

&nbsp; &nbsp; &nbsp; &nbsp; return energy

'''def magnetisation(self):'''

&nbsp; &nbsp; &nbsp; &nbsp; "Return the total magnetisation of the current lattice configuration."

&nbsp; &nbsp; &nbsp; &nbsp; j = 0

&nbsp; &nbsp; &nbsp; &nbsp; for i in np.hstack(self.lattice): #use hstack to vectorise matrix and enumerate through each element and accumulate total

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; j += i

&nbsp; &nbsp; &nbsp; &nbsp; return j

&nbsp; &nbsp; &nbsp; &nbsp; <nowiki></code></nowiki>

Functions are
tested to ensure they work using the ILCheck.py script.

[[File:PW1.png|thumb|center|IlCheck.py energy min, random, and energy max configurations]]

[[File:PW2.png|thumb|center|IlCheck.py example lattice plot; configurations vs. energies]]

<u>'''Introduction to Monte Carlo simulation'''</u>

a)&nbsp;&nbsp;&nbsp;&nbsp; A lattice system of 100 spins, where each can either
be in the +1 or -1 state has 2<sup>100</sup> possible configurations or approximately&nbsp;1.268x10<sup>30</sup>. This is a lot of
configurations to analyse. For example if we consider that our computer can
analyse 1x10<sup>9</sup> configurations per second, to evaluate an average
value such as that of <M><sub>T</sub> over all the configurations would
take roughly 1.269x10<sup>21</sup>
seconds to analyse. This is obviously not feasible and we need to consider a
way to prioritise configurations that can actually be formed. Most of the
configurations will be of a high energy, if we consider a ground state lattice
the probability of obtaining these configurations is given by the Boltzmann
statistical thermodynamic probability. This factor will become very small for high-energy configurations. We can
consider ‘importance sampling’ and the Monte-Carlo simulation, which filters
out configurations where the Boltzmann factor is very small, according to its
probability weighing. 

b)&nbsp;&nbsp;&nbsp;&nbsp;
MonteCarloStep(T)
and statistics() functions are completed see source code below.

<nowiki><code></nowiki>

'''def montecarlostep(self, temperature):'''

&nbsp; &nbsp; &nbsp; &nbsp; # complete this function so that it performs a single Monte Carlo step

&nbsp; &nbsp; &nbsp; &nbsp; #the following two lines will select the coordinates of the random spin for you

&nbsp; &nbsp; &nbsp; &nbsp; random_i = np.random.choice(range(0, self.n_rows))

&nbsp; &nbsp; &nbsp; &nbsp; random_j = np.random.choice(range(0, self.n_cols))

&nbsp; &nbsp; &nbsp; &nbsp; E0_energy = self.energy()

&nbsp; &nbsp; &nbsp; &nbsp; E0_magnetisation = self.magnetisation()

&nbsp; &nbsp; &nbsp; &nbsp; self.lattice[random_i, random_j]*=-1 #flip the randomly selected lattice site

&nbsp; &nbsp; &nbsp; &nbsp; E1_energy = self.energy()

&nbsp; &nbsp; &nbsp; &nbsp; E1_magnetisation = self.magnetisation()

&nbsp; &nbsp; &nbsp; &nbsp; energy_diff = (E1_energy - E0_energy)

&nbsp; &nbsp; &nbsp; &nbsp; if energy_diff <= 0: #accept new configuration if lower energy

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; E0_energy = E1_energy

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; E0_magnetisation = E1_magnetisation

&nbsp; &nbsp; &nbsp; &nbsp; elif energy_diff > 0:

&nbsp; &nbsp; &nbsp; &nbsp; #the following line will choose a random number in the rang e[0,1) for you

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; random_number = np.random.random()

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; if random_number <= np.exp(-1 * energy_diff / (temperature)): # if new config is higher then do Boltzmann test to determine whether to accept a higher energy state with probability exponentially decaying with respect to the energy difference

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; E0_energy = E1_energy

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; E0_magnetisation = E1_magnetisation

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; else:

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; self.lattice[random_i, random_j] *= -1 #else reject new higher energy configuration

&nbsp; &nbsp; &nbsp; &nbsp;

&nbsp; &nbsp; &nbsp; &nbsp; self.n_cycles += 1

&nbsp; &nbsp; &nbsp; &nbsp; if self.n_cycles > 2000: #ignore first 2000 cycles which corresponds to the equilibriating period- starting from a purely random state is not a realistic physical state

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; self.equicycles += 1

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; self.E += E0_energy

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; self.E2 += (E0_energy ** 2)

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; self.M += E0_magnetisation

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; self.M2 += (E0_magnetisation ** 2)

&nbsp; &nbsp; &nbsp; &nbsp; return E0_energy, E0_magnetisation

&nbsp; &nbsp; &nbsp; &nbsp;
'''def statistics(self):'''

&nbsp; &nbsp; &nbsp; &nbsp; self.AE = self.E/self.equicycles

&nbsp; &nbsp; &nbsp; &nbsp; self.AE2 = self.E2/self.equicycles

&nbsp; &nbsp; &nbsp; &nbsp; self.AM = self.M/self.equicycles

&nbsp; &nbsp; &nbsp; &nbsp; self.AM2 = self.M2/self.equicycles

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; return self.AE, self.AE2, self.AM, self.AM2, self.equicycles

&nbsp; &nbsp; &nbsp; &nbsp; # 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

<nowiki></code></nowiki>

Functions are tested using ILanim.py scipt and three outputs were validated with outputs just after the system T=1<Tc , L=8x8, reached its expected ground state equilibrium.

[[File:PW3.png|thumb|center|<span lang="EN-GB">Statistics()
output; Averaged quantities: <E> =&nbsp;
-1.34013513514 <E*E> =&nbsp;
2.01581925676 <M> =&nbsp;
0.428986486486 <M*M> =&nbsp;
0.35966849662</span>]]

[[File:PW4.png|thumb|center|Statistics()
output; Averaged quantities: <E> =&nbsp;
-1.34004160637 <E*E> =&nbsp; 1.8864445776
<M> =&nbsp; 0.351325072359 <M*M> =&nbsp; 0.192099453916]]

[[File:PW5.png|thumb|center|Statistics()
output; Averaged quantities: <E> =&nbsp;
-1.48104434907 <E*E> =&nbsp;
2.44324481402 <M> =&nbsp; -0.704175608011 <M*M> =&nbsp; 0.613123379381]]

<span lang="EN-GB">Observations: It
can be seen that the system reached its expected equilibrium state where all
spins are aligned. The number of MonteCarlosteps() function runs needed in
order to reach this state varied from 563-1374. This is explained via the
randomness factor incorporated in the montecarlostep() function, whereby a np.random number between 0-1 is selected and then tested against <= the
energy change of configurations incorporated in its Boltzmann probability; <math>\exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\}</math>Furthermore the systems initial state is a randomly configured numpy array and
has no physical meaning. The ‘random walk’ toward the ground state is therefore
bound to vary due to these two factors. However reaching this state validates
the functions.</span>

This was expected as when T<Tc spontaneous
magnetization occurs. This physical phenomenon is explained by thermodynamics
as the Curie temperature represents the highest temperature at which non-zero
magnetization occurs. It represents a phase transition point for the Ising
Lattice at which the energetic and entropic driving forces of free energy are
balanced. ; Below Tc the energetic factor dominates and it
is therefore favourable for spins to be aligned to produce to maximum number of
spin interactions and hence minimise system disorder and have spontaneous
magnetization. However above which the entropic factor dominates and disorder
is favoured leading to domain misalignment and a tendency to zero lattice
magnetisation, as checkerboard spins lead to maximum disorder and a new
equilibrium for the system is observed and an energy minimum=0, is observed.&nbsp; 

<span lang="EN-GB"><u>'''Accelerating the code'''</u></span>

&nbsp;

<span lang="EN-GB">Both energy() and
magnetisation functions are re-coded.</span>

<span lang="EN-GB">Energy()- Using a
for loop and enumeration inside the python interface to using a numpy.roll
function. This is significantly quicker as it utilises an internal function within
python, compiled in C++, a much lower level language. Instead of looping
through each lattice site and counting the number of interactions with
neighbours individually, the roll function allows the user to quickly measure
all lattice site interactions simultaneously. Furthermore I have demonstrated
the possibility of avoiding lattice-site double counting. Each roll either
along a row/column (row.T), already takes into account two neighbouring edges
interactions. The comparion can be seen in the source code evolution……..</span>

<span lang="EN-GB">Magnetisation()-
Again a python for loop is replaces, this time with a simple np.sum of the
array output from the lattice.&nbsp; The comparison can be seen in the source code evolution; Note the when utilising the roll function, it was found that two rolls were adequate as these account for neighbouring site interactions, without double-counting (old code commented out)</span>

<nowiki><code></nowiki>

'''def energy(self):'''

&nbsp; &nbsp; &nbsp; &nbsp; "Return the total energy of the current lattice configuration."

&nbsp; &nbsp; &nbsp; &nbsp; #roll1 = np.roll(self.lattice, 1, axis = 1) #use bumpy roll in place of inefficient enumeration in base python

&nbsp; &nbsp; &nbsp; &nbsp; #roll2 = np.roll(self.lattice, -1, axis = 1)

&nbsp; &nbsp; &nbsp; &nbsp; #roll3 = np.roll(self.lattice, 1, axis = 0)

&nbsp; &nbsp; &nbsp; &nbsp; #roll4 = np.roll(self.lattice, -1, axis = 0)

&nbsp; &nbsp; &nbsp; &nbsp; #mult1 = np.multiply(self.lattice, roll1) #efficiently compute energy of neighbour interaction by matching up neighbours and multiplying energy states (neighbours with same sign give a 1 and neighbours with different signs -1) - we can then sum positive results to compute the energy 

&nbsp; &nbsp; &nbsp; &nbsp; #mult2 = np.multiply(self.lattice, roll2)

&nbsp; &nbsp; &nbsp; &nbsp; #mult3 = np.multiply(self.lattice, roll3)

&nbsp; &nbsp; &nbsp; &nbsp; #mult4 = np.multiply(self.lattice, roll4)

&nbsp; &nbsp; &nbsp; &nbsp; #energy = -0.5 * J * (np.sum(mult1[mult1>0]) + np.sum(mult2[mult2>0]) + np.sum(mult3[mult3>0]) + np.sum(mult4[mult4>0]))

&nbsp; &nbsp; &nbsp; &nbsp; target_energy = -1.0*np.sum(np.multiply(self.lattice, np.roll(self.lattice, 1, 0))) #more efficient approach that eliminates 2 rolls and 2 multiplies and avoids the double counting of neighbour interactions

&nbsp; &nbsp; &nbsp; &nbsp; target_energy -= 1.0*np.sum(np.multiply(self.lattice, np.roll(self.lattice, 1, 1)))

&nbsp; &nbsp; &nbsp; &nbsp; energy = target_energy

&nbsp; &nbsp; &nbsp; &nbsp; #d = np.sum(d1) + np.sum(d2)

&nbsp; &nbsp; &nbsp; &nbsp; #energy2 = 0.5 * J * np.sum(d)

&nbsp; &nbsp; &nbsp; &nbsp; return energy&nbsp; &nbsp; &nbsp; &nbsp; 

&nbsp; &nbsp; &nbsp; &nbsp;<nowiki> </nowiki>'''def magnetisation(self):'''

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; "Return the total magnetisation of the current lattice configuration."

&nbsp; &nbsp; &nbsp; &nbsp; return np.sum(self.lattice)

<nowiki></code></nowiki>

<span lang="EN-GB">The ILtimetrial.py
scipt is utilised to quantitatively determine the length of time the computer
took to run 2000 MonteCarlo cycles. The results were determined within a 30min time
frame in order to minimise error, due to varying CPU usage throughout the day.
The results for the individual, and combined function accelerations, the
average, and the error involved in these measurements are presented below; each
trial included 10 repeat runs [See Appendix].</span>
{| class="MsoTableGrid"
 
  |  
  
  |
<nowiki>  </nowiki><span lang="EN-GB">Average time (s)</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">Standard  Deviation</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">Standard Error</span>  
 
 |-
  |
<nowiki>  </nowiki><span lang="EN-GB">Slow E(), Slow
<nowiki> </nowiki> M()</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">64.2217106</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">0.190598169129</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">0.0602724332306</span>  
 
 |-
  |
<nowiki>  </nowiki><span lang="EN-GB">Slow E(), Fast
<nowiki> </nowiki> M()</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">63.6513610</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">0.248100941936</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">0.078456406615</span>  
 
 |-
  |
<nowiki>  </nowiki><span lang="EN-GB">Fast E(), Slow
<nowiki> </nowiki> M()</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">2.12632640</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">0.00857236824791</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">0.00271082086051</span>  
 
 |-
  |
<nowiki>  </nowiki><span lang="EN-GB">FAST E(), Fast
<nowiki> </nowiki> M()</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">1.11812940</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">0.00750646746627</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">0.00237375343754</span>  
 
|}
<span lang="EN-GB">As can be seen the
major acceleration was due to the energy() function. This increase is
approximately x30.2, compared to magnetisation()’s x1.009. This difference is
significant and can be attributed to the fact that the original energy function
required two extensive for loops and calculation of interactions, whereas the original
magnetisation function had only one for loop and no enumerate operator.</span>

<span lang="EN-GB"><u>'''The effect of temperature'''</u></span>

<span lang="EN-GB">a)&nbsp;&nbsp;&nbsp;&nbsp;
The
IsingLattice model used so far, originally starts from a random lattice nxn
configuration. This configuration has no physical meaning. When running the
ILanim.py script is can easily be seen that below Tc, which decreases with
increased lattice size (2D), after a certain number of runs tends to an
equilibrium state where all spins are aligned, either up, or down with total
energy being equal to the minimum ground state, E=-DNJ whereby energy per
spin=-2, and magnetisation per spin =. Furthermore it was seen when T>Tc,
whereby system disorder and domain disaligment is favoured that a new
equilibrium state is reached; total energy(/spin) and magnetisation per spin
both equal to zero. Both of these states are equally valid physical states of
the Ising Lattice. It can also be seen from the IL.Finalframe that in the
energy and magnetisation per spin plots vs. number of MC cycles that this
equilibrium stage takes a certain number of runs to reach. Beforehand the
system was in a physically meaningless random state and is in an equilibrating
period. To study the systems physical behaviour we opt to ignore data
corresponding to this region using self.equicycles counter seen. This is seen in the MonteCarloStep() and statistics() function's in the '''Monte-Carlo intro section'''.</span>

[[File:PW6.png|thumb|center|T=0.5 8x8]]

[[File:PW7.png|thumb|center|T=0.5 16x16]]

[[File:PW8.png|thumb|center|T=0.5 32x32]]

Cutoff; <span lang="EN-GB">8x8; <500/480/750, 16x16;
<7100/10000/9700, 32x32; <15000</span>

[[File:PW9.png|thumb|center|T=1.0 8x8]]

[[File:PW10.png|thumb|center|T=1.0 16x16]]

[[File:PW11.png|thumb|center|T=1.0 32x32]]

Cutoff; 8x8; <370/900/650, 16x16; <6300/10200/9500, 32x32; <20000

[[File:PW12.png|thumb|center|T=1.5 8x8]]

[[File:PW13.png|thumb|center|T=1.5 16x16]]

[[File:PW14.png|thumb|center|T=1.5 32x32]]

Cutoff; <span lang="EN-GB">8x8; <1800/1400/1400, 16x16;
<12000/13500/15000, 32x32; <25000</span>

[[File:PW15.png|thumb|center|T=2.0 8x8]]

[[File:PW16.png|thumb|center|T=1.5 16x16]]

[[File:PW17.png|thumb|center|T=1.5 32x32]]

Cutoff; <span lang="EN-GB">8x8; <1250/2050/1000, 16x16;
<25000/14000/15000, 32x32; <75000</span>

[[File:PW18.png|thumb|center|T=2.5 8x8]]

[[File:PW19.png|thumb|center|T=2.5 16x16]]

[[File:PW20.png|thumb|center|T=2.5 32x32]]

Cutoff; <span lang="EN-GB">8x8; <2100/2000/1400, 16x16; <2000, 32x32;
<10000</span>

<span lang="EN-GB">Three repeats were undertaken for each lattice size at each temperature in order to validate any trends using the ILfinalframe.py script. Two trends are immediately seen; as lattice size increases the number of cycles taken to reach a physical
state equilibrium increases. This is due to the random initial configuration
and, how via the MonteCarlo simulation a random walk is assigned in order to
reach these states. For more lattice sites this will take longer as there are
more spins to flip and test. This generally leads to an increase in the order
of magnitude of steps needed in order to reach equilibrium, this is only seen
for&nbsp; T<Tc. This is because above the
Curie temperature the randomly assigned initial configuration is very near the
actual physical disordered equilibrium state, and therefore easy to reach via a
random walk. However below this, the number of cycles taken in the
equilibrating period behaves extensively in relation to the lattice size.
Secondly once the lattice systems have reaches equilibrium the fluctuations
above the equilibrium values increase with temperature. As a result of this
analysis a cut-off point of 4500 cycles was assigned and implemented in the
MonteCarloStep() function using an if statement relation the no. of cycles
elapsed to ones >4500, and only these would be taken into account in the
average data in the statistics() function. This number would eliminate the main
source of error because for the maximum lattice size analysed, 32x32, where the
variation is a maximum, this cut-off value is feasible and therefore is for smaller lattices by virtue of a lesser number of spins needing to be flipped.
Moreover 150000 cycles are taken for averages such that any error would be
diluted in output averages.</span>

<span lang="EN-GB">Furthermore it should be noted that metastable states
do form and can last for upwards of 10000 cycles and increase in likelihood
with lattice size. Because once one starts to form starts it can be very hard
to flip back, especially at large lattice sizes where it’s less probable as there are more
possibilities geometrically for regions to form and be closed in at the edges
ie. DeltaE=8J to change=unlikely.. These states seem to be temperature independent
as they occurred at low and well as T>Tc. Tc appears to occur around T=2.0K
as entropic and energetic driving forces are most balanced here leading to
larger numbers of cycles being needed to reach equilibrium than at other
temperatures. This observation is validated by the later result that Tc is
around 2.24K.</span>

Furthermore 2x2 and 4x4 lattices were analysed even though they are too finite for useful results. In general the Curie temperature appears earlier and the same overall trends occur as for the other lattices, for plots [Appendix].

<span lang="EN-GB">These states do have physical meaning, but will
interfere with averaging data and are a source of error. Some examples of these
include; T=0.5K (8x8, 32x32) and T=2.5K (32x32), these are local minima, not
equilibrium states of the system.</span>

[[File:PW21.png|thumb|center|T=0.5 8x8]]

[[File:PW22.png|thumb|center|T=2.5 32x32]]

[[File:PW23.png|thumb|center|T=2.5 32x32]]

b)     <span lang="EN-GB">a)&nbsp;&nbsp;&nbsp;&nbsp;ILtemperaturerange.py scipt utilised to run 150000 MonteCarlo cycles over a temperature range (0.5,
5.0, 0.1), with error bars to determine size of error and validate the cutoff point of 4500 for an 8x8 lattice. The data outputs are then saved using
np.savetxt(“”, final_data) where final data contains 5 volumns; energy per
spin, energy per spin squared, magnetisation per spin, magnetisation per spin
squared and variance per spin. The output is shown below. This demonstrates the
phase transition at around T=2.2K where energy/s tends from -2 to 0, and
magnetisation per spin tends from +1 to 0. This demonstrates graphically the
two physical equilibrium states of the Ising Lattice, one for ferromagnetic energetic
spin alignment, and T>Tc where entropic disorder of a paramagnet is
favoured. The 8x8 lattice size plot is shown below. Furthermore error bars are plotted (updated ILtemperaturerange source code next section). I determined that each simulation should be 150000 cycles as this would lead to a lot of post-cutoff data that when averaged would dilute any errors.</span>

[[File:PW24.png|thumb|center|ILtemperaturerange.py 8x8]]

<span lang="EN-GB"><u>'''Scaling the system size'''</u></span>

<span lang="EN-GB">This process was
repeated for L=2, 4, 16, 32, in order to obtain the output data in a .data file
for all lattice sizes under analysis. </span>A Python script to make a plot showing the energy ''per
spin'' versus temperature for each of your lattice sizes, whereby the
numpy.loadtxt function was utilized to read and index the column data in the
files and plot the energy per spin vs. temp and magnetization per spin vs. temp
for all lattice sizes.

[[File:PW25.png|thumb|center|2x2]]

[[File:PW26.png|thumb|center|4x4]]

[[File:PW27.png|thumb|center|16x16]]

[[File:PW28.png|thumb|center|32x32]]

<nowiki><code></nowiki>

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 = 150000

times = range(runtime)

temps = np.arange(0.5, 5.0, 0.2)

energies = []

magnetisations = []

energysq = []

magnetisationsq = []

for t in temps:

&nbsp; &nbsp; for i in times:

&nbsp; &nbsp; &nbsp; &nbsp; if i % 100 == 0:

&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; print(t, i)

&nbsp; &nbsp; &nbsp; &nbsp; energy, magnetisation= il.montecarlostep(t)

&nbsp; &nbsp; aveE, aveE2, aveM, aveM2, n_cycles = il.statistics()

&nbsp; &nbsp; energies.append(aveE)

&nbsp; &nbsp; energysq.append(aveE2)

&nbsp; &nbsp; magnetisations.append(aveM)

&nbsp; &nbsp; magnetisationsq.append(aveM2)

&nbsp; &nbsp; #reset the IL object for the next cycle

&nbsp; &nbsp; il.E = 0.0

&nbsp; &nbsp; il.E2 = 0.0

&nbsp; &nbsp; il.M = 0.0

&nbsp; &nbsp; il.M2 = 0.0

&nbsp; &nbsp; il.n_cycles = 0

&nbsp; &nbsp; il.equicycles = 0

yerr = ((((np.array(energysq) - np.array(energies)**2)**0.5)/spins)/(n_cycles**0.5)) #error determination for error bars

yerr2 = ((((np.array(magnetisationsq) - np.array(magnetisations)**2)**0.5)/spins)/(n_cycles**0.5))

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])

enerax.errorbar(temps, np.array(energies)/spins, yerr, fmt='--o')

magax = fig.add_subplot(2,1,2)

magax.set_ylabel("Magnetisation per spin")

magax.set_xlabel("Temperature")

magax.set_ylim([-1.1, 1.1])

magax.errorbar(temps, np.array(magnetisations)/spins, yerr2, fmt='--o')

enerax.plot(temps, np.array(energies)/spins)

magax.plot(temps, np.array(magnetisations)/spins)

pl.show()

Var_ET = ((np.array(energysq) - np.array(energies)**2)/(temps**2))

print(spins)

print(energies)

print(magnetisations)

print(np.array(energysq)/(spins**2))

print(magnetisationsq)

final_data = np.column_stack((temps, (np.array(energies)/(spins)), (np.array(energysq)/(spins**2)), (np.array(magnetisations)/(spins)), (np.array(magnetisationsq)/(spins**2)), (Var_ET/(spins))))

np.savetxt("8x8NEW2.dat", final_data)

<nowiki></code></nowiki>

<span lang="EN-GB">From these graphs it can be seen that the error bars are extremely small, this shows the
variation in my data points due to initial equilibrating-state data values are
insignificant and the cut-off point is justified. Furthermore these lattice
sizes follow the same general phase transition trend as for the 8x8 lattice.
However it can be seen that as the lattice size increases the point of
transition tends to higher temperature values.</span>

<span lang="EN-GB"><u>These lattice size data plots are easily compared utilising the following plot;</u></span>

[[File:PW29.png |thumb|center| Superposition of Lattice energy energy and magnetisation per spin vs. temperature]]

<span lang="EN-GB">From this plot it
is obvious that lattice sizes 2x2 and 4x4 are too small to give reasonable
results, and this is also seen in the linear regression to determine the Tc for
an infinite lattice. Phase transition theory states that an infinite system is
required, in a system like this the heat capacity diverges to infinity and
becomes discontinuous. The plot above shows the region where the finite
lattices converge to resemble infinite systems behaviour. This happens as the
length of the 2D lattice becomes 8 and continues for 16x16 and 32x32 systems.
Therefore the minimum lattice size needed to capture long-range fluctuations is
8x8.&nbsp; Data obtained for intermediate
lattice sizes would help to pinpoint the exact length at which this behaviour
starts to occur.</span>

This required using the np.savetxt() function to import my data and superimpose the plots ontop of each other, source code as follows;

<nowiki><code></nowiki>

from IsingLattice import *

from matplotlib import pylab as pl

import numpy as np

a = np.loadtxt("2x2NEW.dat")

b = np.loadtxt("4x4NEW.dat")

c = np.loadtxt("8x8NEW.dat")

d = np.loadtxt("16x16NEW.dat")

e = np.loadtxt("32x32NEW.dat")

j = [a, b, c, d, e] # incorporates imported data into a list

fig = pl.figure()

for i in j:

&nbsp; &nbsp; enerax = fig.add_subplot(2, 1, 1)

&nbsp; &nbsp; enerax.set_ylabel("Energy per spin (J)")

&nbsp; &nbsp; enerax.set_xlabel("Temperature (K)")

&nbsp; &nbsp; magax = fig.add_subplot(2, 1, 2)

&nbsp; &nbsp; magax.set_ylabel("Magnetisation per spin")

&nbsp; &nbsp; magax.set_xlabel("Temperature (K)")

&nbsp; &nbsp; enerax.plot(i[:, 0], i[:, 1]) # indexing relevant colums from data, 1st col temp, 2nd energy per spin

&nbsp; &nbsp; magax.plot(i[:, 0], i[:, 3]) # as above, but 3rd column magnetisation per spin

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

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

pl.show()

<nowiki></code></nowiki>

Furthermore we can now conclusively talk about the qualitative behaviour of our Ising Lattice. 

T<Tc; Low temperature ferromagnetic phase. There is a dominance of one spin throughout the Monte-Carlo cycles e.g. s<sub>i </sub>= +1. Also flucuations of the other spin occurring only happens in isolated clusters. These increase in frequency with lattice size, but the system always tends to it ground state after sufficient cycles. Furthermore the number and size of these clusters also increases as the temperature approaches Tc. In addition it is seen on the magnetisation per spin vs. temp plots that fluctuations from an average <M> value of +/- 1 are extremely rare. Finally it is seen on the energy per spin vs. time plot that the energy always tends to -2, this is because below the critical temperature the energetic factor of the systems free energy dominates leading to long range spin ordering and spontaneous magnetisation as a result.

T≈Tc; Near the critical Curie temperature. This is the point at which the temperature of the system is at the highest it can be where spontaneous magnetisation can occur. Significant fluctuations are seen at all lattice sizes (T≈2.5K) . Furthermore there is no clear dominance of a particular spin state. This behaviour is due to the balance of energetic and entropic component's of the systems free energy, with no clear factor dominating. And system takes a long time to reach a stable equilibrium state.

T>Tc; High temperature paramagnetic phase. No dominance of either spin state and their are no obvious patterns to domain alignment. This is because the entropic component of the systems free energy is now dominant and the system favors a disordered state resulting in no long-range spin alignment and no spontaneous magnetisation. <M> and <E> ≈ 0. This results because the system has passed through a order-disorder phase transition point at its Curie temperature, this is lattice size dependant- smaller lattice, lower Tc.

'''Physical reason for this behaviour:''' For a general 2D Ising Lattice, an ordered state has a net energy of -NJ and a net magnetisation of +1. A disordered state resulting from a randomly induced spin flip as a result of the Monte Carlo function leads to a configuration of energy (-N + 4)J, which is N<sup>1/2</sup> parts out of N higher than previously. For large lattice sizes and higher dimensions this energy because of the boltzmann factors inverse exponential relation with energy change is sufficient to stabilise the more stable state. However this stabilisation at higher temperatures(>Tc) is overcome due to the entropic contribution to the free energy.

<span lang="EN-GB"><u>'''Calculating the heat capacity'''</u></span>

<span lang="EN-GB">This is required as it follows from Lars Onsager’s observation that the heat capacity is very
strongly peaked at the Curie temperature of the 2D Ising Lattice, we therefore can use the variance of our average energy data about its equilibrium value to
calculate how our systems heat capacity varies with temperature and therefore locate the T<sub>C,L</sub> for each lattice size and then derive Tc, </span><sub>∞. </sub>However first it is needed to show the relationship between variance and the mean of the average energy data squared and its squared mean. This is derived below.

[[File:PWX.png |thumb|center| Derivation of the variance of a closed system in the canonical ensemble]]

<span lang="EN-GB">Using this relation a plot of heat capacity vs.
temperature can be made for each of the lattice sizes using python; source code and plot output as follows.</span>

<nowiki><code></nowiki>

from IsingLattice import *

from matplotlib import pylab as pl

import numpy as np



a = np.loadtxt("2x2vee.dat")

b = np.loadtxt("4x4vee.dat")

c = np.loadtxt("8x8vee.dat")

d = np.loadtxt("16x16vee.dat")

e = np.loadtxt("32x32vee.dat")


j = [a, b, c, d, e]

for i in j

&nbsp; &nbsp; pl.plot(i[:, 0], i[:, 5]) # same as above, column 5 is heat capacity

pl.legend()

pl.title("Heat Capacity vs. Temperature for various lattice sizes")

pl.xlabel("Temperature (K)")

pl.ylabel("Heat Capacity(J/K)")

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

pl.show()

<nowiki></code></nowiki>

[[File:PW30.png |thumb|center|Heat Capacity vs. temperature, spiking occurs at the Curie Temperate for the various finite lattice sizes]]

Observations: As predicted by Onsager,
the heat capacity spikes for each lattice size around a particular region where
the Curie temperature for each lattice is. Overall this is in the region 2.0-2.8K.
Furthermore it can be seen that for larger lattice sizes this spike is more
pronounced and is accounted for by the fact that the heat capacity for a larger
lattice will be higher at it is an extensive property by definition. More
importantly this spike occurs at slightly lesser temperatures as lattice size
increases. [Note] Heat capacity values plotted are not ‘per spin’ in order to
give a realistic value for the entire lattice systems. It is seen that as the
lattice size increases the peak in the heat capacity during the phase
transition at Tc for the finite lattices is diverging.

<span lang="EN-GB">'''<u>Locating the Curie temperature</u>'''</span>

a) <span lang="EN-GB">As a validity
check for the data obtained through my IsingLattice () model, data is plotted
against C++ data obtained from much more extensive cycles. This was done for
all lattice sizes, the resulting plot for the 8x8 lattice data set is shown
below and the in the appendix for the other lattices.</span>

<span lang="EN-GB">This utilised a
tuple and a figure label list indented in a for loop and plotted 5 subplots per
figure, each data file was imported and given and assigned a specific variable to be used in a list. This was in top to bottom order; energy per spin vs. t, energy squared per spin
vs. t, magnetisation per spin vs. t, magnetisation squared per spin vs. t and
variance per spin vs. t, where t=temperature.</span>

[[File:PW31.png |thumb|center| 8x8 Python data vs. C++ data: Graph order top to bottom; E, E<sup>2</sup>, M, M<sup>2</sup>and heat capacity. All of these are 'per spin' and all are plotted in the temperature range with contains the Curie temperature.]]

Observations; &nbsp;It can be seen that the fit of the python data
compared to the C++ data is good, but deviates as the Curie temperature for
larger lattice sizes. Furthermore the energy squared per spin data from python
does not pick up on long-range fluctuations as the system approaches the curie
temperature and this is seen for all lattice sizes. Other lattice plots are
seen in the appendix. Overall the data fit between the two files is good,
leading to confidence in the data obtained via the IsingLattice python function. The code utilised for this is shown below. Furthermore a spike in the variance/fluctuations is shown very clearly in the C++ data especially in the magnetisation vs. temp graphs, with more runs of the MonteCarlo() cycles the Python script would be able to demonstrate this also.

<nowiki><code></nowiki>

from IsingLattice import *

from matplotlib import pylab as pl

import numpy as np


a = np.loadtxt("2x2NEW2.dat")

b = np.loadtxt("4x4NEW2.dat")

c = np.loadtxt("8x8NEW2.dat")

d = np.loadtxt("16x16NEW2.dat")

e = np.loadtxt("32x32NEW2.dat")


path_cpp = "../C++_data/"

f = np.loadtxt(path_cpp + "2x2.dat")

g = np.loadtxt(path_cpp + "4x4.dat")

h = np.loadtxt(path_cpp + "8x8.dat")

i = np.loadtxt(path_cpp + "16x16.dat")

j = np.loadtxt(path_cpp + "32x32.dat")


x = [(a, f), (b, g), (c, h), (d, i), (e, j)] #create list of tuples for each combination of python and C++ lattice size datasets

latt_size = ['2x2', '4x4', '8x8', '16x16', '32x32']

for m, (n, o) in enumerate(x):

&nbsp; &nbsp; fig = pl.figure(m + 1)

&nbsp; &nbsp; fig.suptitle(latt_size[m], fontsize=20) #index's the tuple with the latt size list in the for loop

&nbsp; &nbsp; for p in range(5): # p is used enumerate through each of the 5 measures (starting at p+1 since p starts at zero)

&nbsp; &nbsp; &nbsp; &nbsp; pl.subplot(5, 1, p)

&nbsp; &nbsp; &nbsp; &nbsp; pl.plot(n[:, 0], n[:, p + 1])

&nbsp; &nbsp; &nbsp; &nbsp; pl.plot(o[:, 0], o[:, p + 1])

&nbsp; &nbsp; &nbsp; &nbsp; pl.legend(['Python', 'C++'], loc='upper left')

&nbsp; &nbsp; pl.show()

&nbsp; &nbsp; &nbsp; &nbsp; <nowiki></code></nowiki>

b) Further analysis
is undertaken by fitting the heat capacity vs. t data to a polynomial in order
to obtain a functional representation of the graphs via polynomial
coefficients. This is firstly done for the whole of the temperature range, via
importing the data and indexing the columns appropriately and using np.min/max
functions and the np.polyfit(x, y, degree) operator and then unpacking the
coefficients using np.polyval to evaluate the coefficients across the entire
temperature range; 0.5-5K. The fit was plotted using a np.linspace function to
generate 1000 data points in the temperature interval specified. &nbsp;The script is written is a way that uses an‘i’
variable to pick out specific lattice size data. Source code found below.

<nowiki><code></nowiki>

from matplotlib import pylab as pl

import numpy as np


a = np.loadtxt("2x2NEW2.dat")

b = np.loadtxt("4x4NEW2.dat")

c = np.loadtxt("8x8NEW2.dat")

d = np.loadtxt("16x16NEW2.dat")

e = np.loadtxt("32x32NEW2.dat")

fig = pl.figure()

i = c #can pick out a specific lattice data set for analysis

T = i[:, 0]

HC = i[:, 5]

fit = np.polyfit(T, HC, 4)#polyfit fits a polynomial a + bx + cx^2 + dx^3 to the data points, saving the coefficients

T_min = np.min(T)

T_max = np.max(T)

T_range = np.linspace(T_min, T_max, num = 1000)

fitted_C_values = np.polyval(fit, T_range)#polyval then evaluates the polynomial values in a data set range T_range- giving a new set of data points&nbsp;

enerax = fig.add_subplot(2, 1, 1)

enerax.set_ylabel("Heat capacity no-fit (J/K)")

enerax.set_xlabel("Temperature (K)")

magax = fig.add_subplot(2, 1, 2)

magax.set_ylabel("Heat capacity fitted (J/K)")

magax.set_xlabel("Temperature (K)")

enerax.plot(T, HC)

magax.plot(T_range ,fitted_C_values)

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

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

label = enerax.legend

pl.show()

<nowiki></code></nowiki>

<span lang="EN-GB">This was done and the optimal polynomial fit was of order 4, as the original data resembles a parabola, however the fit needs increases polynomial order to account for the sharper peak. The result for the 2x2 lattice is shown below.</span>

[[File:PW32.png |thumb|centre| Full temperature range polynomial fit of 2x2 lattice]]

Using a polyfit
function for an entire curve is useful, but we are only interested in
evaluating the peak region, and therefore limiting it to inside this region,
and obtaining the output can optimise the fit. Modifying the code to a selected
temperature range, 2.0-2.8k as this is a temperature range where the
heat capacity for all systems peaks, &nbsp;for
heat capacity data fitting and outputting the fitted heat capacity values only
in this region optimises the fit. &nbsp;This
is because the polyfit function is not trying to fit polynomial coefficients to
the whole temperature range, but only to the one of interest. &nbsp;This is achieved by generating the variables
peak_T_range, and fitted_C_values. The source code for this is found below.

<nowiki><code></nowiki>

from matplotlib import pylab as pl

import numpy as np


a = np.loadtxt("2x2NEW2.dat")

b = np.loadtxt("4x4NEW2.dat")

c = np.loadtxt("8x8NEW2.dat")

d = np.loadtxt("16x16NEW2.dat")

e = np.loadtxt("32x32NEW2.dat")


fig = pl.figure()

i = a


T = i[:, 0]

HC = i[:, 5]

T_Cmin = 2.0

T_Cmax = 2.8

T_min = np.min(T)

T_max = np.max(T)

selection = np.logical_and(T > T_Cmin, T < T_Cmax) # boolean indexing to select T values in the range Cmin to Cmax

peak_T_values = T[selection]

peak_C_values = HC[selection]

fit = np.polyfit(peak_T_values, peak_C_values, 3)

peak_T_range = np.linspace(T_Cmin, T_Cmax, num = 1000) # linspace creates 1000 data points in temp range chosen to fit polynomial to.

fitted_C_values = np.polyval(fit, peak_T_range)


enerax = fig.add_subplot(2, 1, 1)

enerax.set_ylabel("Heat capacity no-fit (J/K)")

enerax.set_xlabel("Temperature (K)")

magax = fig.add_subplot(2, 1, 2)

magax.set_ylabel("Heat capacity fitted (J/K)")

magax.set_xlabel("Temperature (K)")

enerax.plot(T, HC)

magax.plot(peak_T_range, fitted_C_values)

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

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

label = enerax.legend


C_max = np.max(fitted_C_values)

T_max = peak_T_range[fitted_C_values == C_max]


pl.show()

print(C_max, T_max)&nbsp;

<nowiki></code></nowiki>&nbsp;

<span lang="EN-GB">This was done and the optimal polynomial fit was of order 4 for the same reasons as before. The output for the 2x2 lattice is shown below.</span>

[[File:PW33.png |thumb|centre| Restricted phase transition peak temperature range polynomial fit of 2x2 lattice]]

<span lang="EN-GB">Using the new optimised polynomial peak fits;
the maximum heat capacity and its corresponding temperature data points can be
abstracted. This is achieved incorporating Cmax = np.max(fitted_C _values)
followed by </span>Tmax = peak_T_range[fitted_C_values == Cmax]. This data was outputted and saved in a text file, "python curie temp approx. txt", containing two columns; lattice length l and Tmax. Source code utilised is as used previously in this section, but with the following added functionality to obtain the key data. Text file data shown in table below.
{| class="MsoTableGrid"
 
  |
<nowiki>  </nowiki><span lang="EN-GB">2D Ising lattice
<nowiki> </nowiki> side length</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">Temperature at
<nowiki> </nowiki> which C is a max (Tc)</span>  
 
 |-
  |
<nowiki>  </nowiki><span lang="EN-GB">2</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">2.46526527</span>  
 
 |-
  |
<nowiki>  </nowiki><span lang="EN-GB">4</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">2.38598599</span>  
 
 |-
  |
<nowiki>  </nowiki><span lang="EN-GB">8</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">2.34434434</span>  
 
 |-
  |
<nowiki>  </nowiki><span lang="EN-GB">16</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">2.30190190</span>  
 
 |-
  |
<nowiki>  </nowiki><span lang="EN-GB">32</span>
  
  |
<nowiki>  </nowiki><span lang="EN-GB">2.19539540</span>  
 
|}

<nowiki><code></nowiki>

C_max = np.max(fitted_C_values)

T_max = peak_T_range[fitted_C_values == C_max]

print(C_max, T_max)&nbsp;

<nowiki></code></nowiki>&nbsp;

Utilising the following scaling relation relating Tc for a finite
lattice to Tc for an infinite lattice;&nbsp;

[[File:PWX2.png |thumb|centre| Scaling relation of Tc for a finite 2D Ising Lattice to theoretical value for an infinite lattice as determined by Onsager]]

<span lang="EN-GB">A linear regression script was written in python via importing the txt file data and
using code </span>m,b = np.polyfit(x, y, 1), and then outputting b the y-intercept
which is equal to . Source code and output line plot below.&nbsp;

[[File:PW34.png |thumb|centre| Linear regression plot using Tc scaling relation to determine the theoretical infinite 2DIL Curie temperature.]]

<nowiki><code>&nbsp;</nowiki>

from matplotlib import pylab as pl&nbsp;

import numpy as np

data = np.loadtxt("python curie temp approx.txt")

fig = pl.figure()

Tc = data[:, 1]

Lattlen = 1/(np.array(data[:, 0])) # determination of data sets used in the y = mx + c plot

pl.plot(Lattlen, Tc)

m,b = np.polyfit(Lattlen, Tc, 1)

print(b)

pl.title("Linear regression to determine Tc for an infinite lattice")#dont need A constant as we are only interested in the y-intercept for Tc infinite 2D latt

pl.ylabel("Tc,L (K)")

pl.xlabel("1/L")

pl.show()

<nowiki></code></nowiki>

<span lang="EN-GB">Y-intercept = 2.2463797125K. Lars Onsager’s determined the theoretical Curie temperature for an infinite 2D Ising lattice
to be 2.269185314K, the error in this result is  +/- 1%. Onsager utilised the following relation [1].</span>

[[File:PWX3.png |thumb|centre| Onsager relation.]]

<span lang="EN-GB">'''<u>Error analysis</u>'''</span>

<span lang="EN-GB">The error in the
result for the theoretical infinite 2D Ising lattice value can be attributed to
how the T<sub>C</sub> data for the finite lattices l= 2, 4, 8, 16, 32 was obtained.
Firstly the cut-off point of 4500 Monte Carlo runs may be not be exactly equal
to the equilibrating period for all lattice sizes in particular the 32x32
lattice where this period is longest. However the error bars for the data were
small so this can only be a small source of error. Secondly intermediate local
minima meta-stable states may have formed throughout the runs where the
final_data was obtained. This would lead to some of the incorporated data not
being fluctuations about the systems equilibrium point. In addition more runs
could have been taken when calculating the average data values.</span>

<span lang="EN-GB">More significant
sources of error include using a polynomial fit to the heat capacity data at
the phase transition point. For the Ising lattice transitioning between its
ferromagnetic and paramagnetic state causes the heat capacity to diverge to
infinity and become discontinuous as it peaks at Tc. This is dampened for
finite lattices but fitting a polynomial to this as the lattice size increases
and converges to infinite behaviour where long-range fluctuations are visible,
after 8x8, leads to a bad fit. Finally as it is known that the 2x2 and 4x4
lattice sizes weren’t large enough to lead to reasonable physical results, data
values for these sets should be eliminated from the linear regression, as it is
seen that they deviate from linearity in the scaling relation. Omitting these
sets and running trials for intermediate lattice sizes between 8 and 32 would
lead to a better linear fit and therefore &nbsp;data
value that was closer to the theoretical value determined via Onsager’s
relation. However upon removing the 2x2 and 4x4 data sets, the data deviates
from the theoretical &nbsp;value downwards. Upon inspection of the plots
of the python vs. C++ data&nbsp; in the
appendix, it is seen that the heat capacity per spin values for the 32x32
lattice deviate a lot&nbsp; and this is the
source of error. This can be accounted for by the fact that the cut-off despite
low error bars, actually is too short. It is known 32x32 lattices at T=2.0K
take a very long time to equilibrate on the 100000’s of runs and this was not
accounted for in the cut-off. This resulted in a value of 2.012K, a worse approximation for the
theoretical value. However it must be noted this is purely due to error in the
16x16 and 32x32 data sets as my plots show these in fact converge to show
long-range fluctuations as they should.</span>

<span lang="EN-GB">'''<u>References</u>'''</span>

[1] - &nbsp;R. Albert, and A-L. Barabási, Rev. Mod. Phys.
74 (2002) 4797

<span lang="EN-GB">'''<u>Appendix</u>'''</span>

<span lang="EN-GB">'''Accelerating
code run speed data (raw)'''</span>

<span lang="EN-GB">'''Slow mag, slow
energy'''</span>

<span lang="EN-GB">Took 64.217159s</span>

<span lang="EN-GB">Took
64.27629900000001s</span>

<span lang="EN-GB">Took
64.19219800000002s</span>

<span lang="EN-GB">Took 64.670774s</span>

<span lang="EN-GB">Took
63.92942700000003s</span>

<span lang="EN-GB">Took 64.19099799999998s</span>

<span lang="EN-GB">Took 64.314731s</span>

<span lang="EN-GB">Took
64.12510699999996s</span>

<span lang="EN-GB">Took
64.09174399999995s</span>

<span lang="EN-GB">Took
64.20866899999999s</span>

<span lang="EN-GB">AVERAGE=
64.2217106</span>

<span lang="EN-GB">'''Fast mag slow
energy(new)'''</span>

<span lang="EN-GB">Took
64.06209899999999s</span>

<span lang="EN-GB">Took
63.507130000000004s</span>

<span lang="EN-GB">Took
63.52340700000002s</span>

<span lang="EN-GB">Took
63.70164199999999s</span>

<span lang="EN-GB">Took 63.567950999999994s</span>

<span lang="EN-GB">Took
64.03124500000001s</span>

<span lang="EN-GB">Took
63.29038899999995s</span>

<span lang="EN-GB">Took
63.80433400000004s</span>

<span lang="EN-GB">Took
63.47986399999991s</span>

<span lang="EN-GB">Took
63.54554900000005s</span>

<span lang="EN-GB">AVERAGE= 63.651361</span>

<span lang="EN-GB">'''Slow mag fast
energy'''</span>

<span lang="EN-GB">Took 2.122013s</span>

<span lang="EN-GB">Took 2.12318s</span>

<span lang="EN-GB">Took
2.1139779999999995s</span>

<span lang="EN-GB">Took 2.128309s</span>

<span lang="EN-GB">Took 2.1437489999999997s</span>

<span lang="EN-GB">Took
2.125437999999999s</span>

<span lang="EN-GB">Took
2.1378489999999992s</span>

<span lang="EN-GB">Took
2.1242959999999975s</span>

<span lang="EN-GB">Took 2.121333s</span>

<span lang="EN-GB">Took
2.123118999999999s</span>

<span lang="EN-GB">AVERAGE= 2.1263264</span>

<span lang="EN-GB">'''Fast mag fast
energy'''</span>

<span lang="EN-GB">Took
1.1121390000000002s</span>

<span lang="EN-GB">Took
1.1222689999999997s</span>

<span lang="EN-GB">Took
1.1211180000000005s</span>

<span lang="EN-GB">Took 1.1128610000000005s</span>

<span lang="EN-GB">Took
1.1268589999999996s</span>

<span lang="EN-GB">Took
1.1132969999999993s</span>

<span lang="EN-GB">Took
1.1050430000000002s</span>

<span lang="EN-GB">Took
1.1151069999999998s</span>

<span lang="EN-GB">Took 1.124898s</span>

<span lang="EN-GB">Took 1.127703s</span>

<span lang="EN-GB">AVERAGE= 1.1181294</span>

'''<u>2x2 &4x4 Cutoff analysis plots</u>'''

[[File:PW39.png | thumb | left | T=3.0 2x2 ]]

[[File:PW40.png | thumb | left | T=3.0 4x4]]

<span lang="EN-GB">'''<u>Other C++ vs.↵python data plots</u>'''</span>

[[File:PW35.png | thumb | center| 2x2 Python vs. C++]]

[[File:PW36.png | thumb | center | 4x4 Python vs. C++]]

[[File:PW37.png | thumb | center | 16x16 Python vs. C++]]

[[File:PW38.png | thumb | center | 32x32 Python vs. C++]]