View source for Rep:Wm1415CMP

'''<big>This is William Micou's report on the compulsory CMP experiment: Ising Model.</big>'''

==Introduction to the Ising Model==

'''<big>TASK</big>: Show that the lowest possible energy for the Ising model is <math>E\ =\ -DNJ</math>, where <math>D</math> is the number of dimensions and <math>N</math> is the total number of spins. What is the multiplicity of this state? Calculate its entropy.'''

In the absence of an external magnetic field, the energy of the lattice depends on the interaction energy <math>E</math>between adjacent sites<sup>[1]</sup>:

<center><math>E = - \frac{1}{2} J \sum_i^N \sum_{j\  \in\  \mathrm{neighbours}\left(i\right)} s_i s_j</math></center>

Thus, the energy is minimised when neighbouring spins are aligned.  <math>J</math> is the exchange energy and comes about from the Pauli Exclusion Principle<sup>[2]</sup>: electrons cannot occupy the same quantum state.  Electrons on neighbouring atoms with spins aligned cannot occupy the same space - electron repulsion is reduced and the total energy is lowered<sup>[2]</sup>.

The state with the lowest possible energy is the state with all spins aligned.  To calculate its energy, the total number of interactions must be counted.

For an infinite lattice, there are 2 nearest neighbours per dimension (degree of freedom): number of pairs of interactions <math> =\ 2ND</math>.

Hence, total number of interactions <math>=\frac{1}{2}\times2ND = ND</math>.

Since all spins are aligned, every interaction will contribute <math>-J</math> to the energy:

<math>E=-J\times ND=-DNJ</math> as required.

This lowest energy state is comprised of 2 microstates: all spins up, and all spins down.  Thus, the entropy <math>S</math> is:

<math>S=k_{B}\ln\Omega=k_{B}\ln2</math>

'''<big>TASK</big>: Imagine that the system is in the lowest energy configuration. To move to a different state, one of the spins must spontaneously change direction ("flip"). What is the change in energy if this happens (<math>D=3,\ N=1000</math>)? How much entropy does the system gain by doing so?'''

Lowest energy state is -3000J

The spin being flipped has 6 nearest neighbours.  Its 6 interactions with the lattice contribute -6J to the energy.  After flipping, the spin now contributes +6J to the energy.

Energy of one spin flipped state: -2988J

Difference of +12J

There are 2000 microstates for this system: 1000 all but one spin up states, or 1000 all but one spin down states.

<math>S=k_{B}\ln\Omega=k_{B}\ln2000-k_{B}\ln2=k_{B}\ln1000</math>

'''<big>TASK</big>: Calculate the magnetisation of the 1D and 2D lattices in figure 1. What magnetisation would you expect to observe for an Ising lattice with <math>D = 3,\ N=1000</math> at absolute zero?'''

1D lattice: magnetisation = +1

2D lattice: magnetisation = +1

At absolute zero, entropic effects vanish and the system will remain in the lowest energy state, with all spins aligned.  Thus the magnetisation would  be either +N or -N (+1000 or -1000).

==Calculating the energy and magnetisation==

'''<big>TASK</big>: complete the functions energy() and magnetisation(), which should return the energy of the lattice and the total magnetisation, respectively. In the energy() function you may assume that <math>J=1.0</math> at all times (in fact, we are working in ''reduced units'' in which <math>J=k_B</math>, but there will be more information about this in later sections). Do not worry about the efficiency of the code at the moment &mdash; we will address the speed in a later part of the experiment.'''

Below is the initial, inefficient code for energy() and magnetisation() using indented for loops.

<pre>
    def energy(self):
        "Return the total energy of the current lattice configuration."
        energy = 0.0
        for i in range(0,self.n_rows):
            for j in range(0,self.n_cols):
        		start = int(self.lattice[i][j])
        		top = int(self.lattice[i-1][j])
        		bot = int(self.lattice[(i+1)%self.n_rows][j])
        		left = int(self.lattice[i][j-1])
        		right = int(self.lattice[i][(j+1)%self.n_cols])
        		total_interaction = (start*top)+(start*bot)+(start*left)+(start*right)
         		energy = energy + (-0.5*total_interaction)
        return energy

    def magnetisation(self):
        "Return the total magnetisation of the current lattice configuration."
        magnetisation = 0
        for i in range(0,self.n_rows):
        	for j in range(0,self.n_cols):
        		magnetisation = magnetisation + int(self.lattice[i][j])
        return magnetisation
</pre>

The energy() method loops over every spin in the lattice and adds up every interaction.  The modulus allows for boundary conditions for spins at the edges of our lattice.  The magnetisation() method simply adds up every spin value in the lattice.

'''<big>TASK</big>: Run the ILcheck.py script from the IPython Qt console using the command <pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

[[File:wm1415cmp_first_check.png|500px|thumb|center|ILcheck.py generates 3 lattices and compares the expected energy and magnetisation with the values calculated from the above methods.  The initial energy() and magnetisation() methods are in agreement with expected values: the code works.]]

==Introduction to the Monte Carlo simulation==

'''<big>TASK</big>: How many configurations are available to a system with 100 spins? To evaluate these expressions, we have to calculate the energy and magnetisation for each of these configurations, then perform the sum. Let's be very, very, generous, and say that we can analyse <math>1\times 10^9</math> configurations per second with our computer. How long will it take to evaluate a single value of <math>\left\langle M\right\rangle_T</math>?'''

For a system of 100 spins, there are <math>2^{100}</math> configurations.  Evaluating <math>10^9</math> calculations per second, it would take around <math>1.3\times10^{21}</math> seconds to finish the calculation - around 4 orders of magnitude larger than the age of the Universe.  Clearly, this is impractical: importance sampling is a method of circumventing this.

'''<big>TASK</big>: Implement a single cycle of the above algorithm in the montecarlocycle(T) function. This function should return the energy of your lattice and the magnetisation at the end of the cycle. You may assume that the energy returned by your energy() function is in units of <math>k_B</math>! Complete the statistics() function. This should return the following quantities whenever it is called: <math><E>, <E^2>, <M>, <M^2></math>, and the number of Monte Carlo steps that have elapsed.'''

<pre>
    def montecarlostep(self, T):
        "Performs a single cycle of the Monte Carlo algorithm presented in the lab script."
        # first: increase the number of steps elapsed
        self.n_cycles = self.n_cycles + 1
        # introduce a variable that will indicate whether the new configuration should be rejected or kept
        accept_configuration = True
        # energy of the original lattice (E_0)
        energy_0 = self.energy()
        # select coordinates for a random spin in the lattice
        random_i = np.random.choice(range(0, self.n_rows))
        random_j = np.random.choice(range(0, self.n_cols))
        # flip the spin in the lattice
        self.lattice[random_i, random_j] = -1 * self.lattice[random_i, random_j]
        # calculate new lattice energy (E_1)
        energy_1 = self.energy()
        # calculate the difference in energy
        delta = energy_1 - energy_0
        # if the new configuration is lower in energy, keep the new configuration
        if delta < 0:
            accept_configuration = True

        # choose a random number in the range [0,1)
        random_number = np.random.random()
        boltzmann_factor = np.exp(-delta / T)
        # if the new configuration is higher in energy, compare random number to Boltzmann factor
        if delta >= 0:
            if random_number <= boltzmann_factor:
                accept_configuration = True
            else:
                accept_configuration = False

        if accept_configuration == False:
            # return to original lattice
            self.lattice[random_i, random_j] = -1 * self.lattice[random_i, random_j]

        # keep running totals of data for n_cycles > 250000 steps
        discount = 250000
        if self.n_cycles > discount:
            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
        # return the post-MMC step energy and magnetisation as a list
        return [self.energy(), self.magnetisation()]

    def statistics(self):
        "Calculates the average values of E, E*E (E2), M, M*M(M2) and returns them."
        # Discount the first 250000 steps.
        discount = 250000
        avg_E = self.E / (self.n_cycles-discount)
        avg_E2 = self.E2 / (self.n_cycles-discount)
        avg_M = self.M / (self.n_cycles-discount)
        avg_M2 = self.M2 / (self.n_cycles-discount)
        return [avg_E, avg_E2, avg_M, avg_M2, self.n_cycles]
</pre>

This code comes from later in the assignment - when the methods are modified to discard the first 250000 steps from the running averages of energy and magnetisation.

'''<big>TASK</big>: If <math>T < T_C</math>, do you expect a spontaneous magnetisation (i.e. do you expect <math>\left\langle M\right\rangle \neq 0</math>)? When the state of the simulation appears to stop changing (when you have reached an equilibrium state), use the controls to export the output to PNG and attach this to your report. You should also include the output from your statistics() function.'''

The original ILanim.py expected only 4 entries from the statistics() method:

<pre>
    E, E2, M, M2 = il.statistics()
    print("Averaged quantities:")
    print("E = ", E/spins)
    print("E*E = ", E2/spins/spins)
    print("M = ", M/spins)
    print("M*M = ", M2/spins/spins)
</pre>

I modified the code to expect 5 entries from the statistics() method, which also outputs the number of cycles elapsed:

<pre>
    E, E2, M, M2, n = il.statistics()
    print("Averaged quantities:")
    print("E = ", E/spins)
    print("E*E = ", E2/spins/spins)
    print("M = ", M/spins)
    print("M*M = ", M2/spins/spins)
    print("No. cycles = ", n)
</pre>

The figure below shows the outcome of running ILanim.py for an 8x8 lattice, at temperature = 1.  This is below the critical temperature: the figure shows the system converged to the lowest energy state, with all spins aligned.  The atoms do not have sufficient thermal energy to flip their spins: the exchange effect dominates over the entropic effect at low temperatures - spontaneous magnetisation is expected.

<center>[[File:wm1415_Second_animpy.png|600px]]</center>

The output of statistics():

<pre>
    Averaged quantities:
    E =  -1.44239045383
    E*E =  2.44085607394
    M =  0.732834507042
    M*M =  0.641262287265
    No. cycles =  639
</pre>

These are averaged quantities over the course of the rearrangement of the lattice: no cycles were discarded in the calculation of the averages (discount variable set to zero).

==Accelerating the code==

'''<big>TASK</big>: Use the script ILtimetrial.py to record how long your ''current'' version of IsingLattice.py takes to perform 2000 Monte Carlo steps. This will vary, depending on what else the computer happens to be doing, so perform repeats and report the error in your average!'''
 
9.74 <math>\pm</math> 0.26 seconds before optimisation (10 trials)

'''<big>TASK</big>: Look at the documentation for the [http://docs.scipy.org/doc/numpy/reference/generated/numpy.sum.html NumPy sum function]. You should be able to modify your magnetisation() function so that it uses this to evaluate M. The energy is a little trickier. Familiarise yourself with the NumPy [http://docs.scipy.org/doc/numpy/reference/generated/numpy.roll.html roll] and [http://docs.scipy.org/doc/numpy/reference/generated/numpy.multiply.html multiply] functions, and use these to replace your energy double loop (you will need to call roll and multiply twice!).'''

<pre>
    def energy(self):
        "Return the total energy of the current lattice configuration."
        # Roll the lattice vertically and multiply with original to find total vertical interactions
        vertical_interactions = np.multiply(self.lattice, np.roll(self.lattice, 1, axis=0))
        # Roll the lattice horizontally and multiply with original to find total horizontal interactions
        horizontal_interactions = np.multiply(self.lattice, np.roll(self.lattice, 1, axis=1))
        return (-1)*(np.sum(vertical_interactions + horizontal_interactions))

    def magnetisation(self):
        "Return the total magnetisation of the current lattice configuration."
        # Sum every spin in the lattice to obtain the total magnetisation.
        return np.sum(self.lattice)
</pre>


'''<big>TASK</big>: Use the script ILtimetrial.py to record how long your ''new'' version of IsingLattice.py takes to perform 2000 Monte Carlo steps. This will vary, depending on what else the computer happens to be doing, so perform repeats and report the error in your average!'''

0.21 <math>\pm</math> 0.01 seconds after second optimisation (10 trials) (new code)

==The effect of temperature==

'''<big>TASK</big>: The script ILfinalframe.py runs for a given number of cycles at a given temperature, then plots a depiction of the ''final'' lattice state as well as graphs of the energy and magnetisation as a function of cycle number. This is much quicker than animating every frame! Experiment with different temperature and lattice sizes. How many cycles are typically needed for the system to go from its random starting position to the equilibrium state? Modify your statistics() and montecarlostep() functions so that the first N cycles of the simulation are ignored when calculating the averages. You should state in your report what period you chose to ignore, and include graphs from ILfinalframe.py to illustrate your motivation in choosing this figure.'''
{| class="wikitable"
!Lattice
!Temperature
!ILfinalframe.py
!Comments
|-
|6x6
|1.0
|[[File:wm1415_Final_frame5_1.png|200px|thumb|left]][[File:wm1415_Final_frame5_2.png|200px|thumb|right]]
|System stabilises after ~170 steps
|-
|8x8
|1.0
|[[File:wm1415_Final_frame1_1.png|200px|thumb|left]][[File:wm1415_Final_frame1_2.png|200px|thumb|right]]
|System stabilises after 600-700 steps
|-
|8x8
|2.0
|[[File:wm1415_Final_frame_9_1.png|200px|thumb|left]][[File:wm1415_Final_frame_9_2.png|200px|thumb|right]]
|System stabilises after ~2600 steps
|-
|8x8
|3.0
|[[File:wm1415_Final_frame8_1.png|200px|thumb|center]]
|Above the critical temperature: spins flip rapidly and the system starts fluctuating from the first step.
|-
|10x10
|0.5
|[[File:wm1415_Final_frame4_1.png|200px|thumb|left]][[File:wm1415_Final_frame4_2.png|200px|thumb|right]]
|System stabilises after ~1800 steps
|-
|10x10
|1.5
|[[File:wm1415_Final_frame3_1.png|200px|thumb|left]][[File:wm1415_Final_frame3_2.png|200px|thumb|right]]
|System stabilises after ~2200 steps
|-
|10x10
|3.0
|[[File:wm1415_Final_frame2_1.png|200px|thumb|center]]
|Above the critical temperature: spins flip rapidly and the system starts fluctuating from the first step.
|-
|15x15
|1.0
|[[File:wm1415_Final_frame7_1.png|200px|thumb|center]]
|Magnetisation appears to converge after 1-2000 steps, but then much later on the final state is found.
|-
|32x32
|0.5
|[[File:wm1415_32x32_t05_clump.png|200px|thumb|left|'Clumping']][[File:wm1415_32x32_t05_noclump2.png|200px|thumb|right|No 'clumping']]
|The larger 32x32 lattice can either find the lowest possible energy state,

or get stuck in a metastable clumped state - even with the same initial conditions.
|-
|32x32
|1.0
|[[File:wm1415_32x32_t1_clump2.png|200px|thumb|left|'Clumping']][[File:Wm1415_32x32_t1_noclump.png|200px|thumb|right|No 'clumping']]
|The larger 32x32 lattice can either find the lowest possible energy state,

or get stuck in a metastable clumped state - even with the same initial conditions.
|-
|32x32
|3.0
|[[File:wm1415_32x32_T3.png|200px|thumb|center]]
|Above the critical temperature, the larger 32x32 lattice reaches a

'final' state after around 5000 steps: many more steps required cf smaller lattices.
|}

The number of steps required to converge to the 'final' state increases with lattice size.  Consider the case below the critical temperature - we expect to observe spontaneous magnetisation and for all spins to align.  For larger lattices, it may take many steps to flip the remaining sites since they are chosen entirely at random - and the probability of choosing the 'correct' site to flip decreases with lattice size.  For the smaller lattice sizes (<16x16), the system behaves as expected within at most a few thousand steps.  

Larger lattices are prone to 'clumping' - the entire system is not in the same spin state, but large regions of spin-aligned atoms form.  These are metastable states: spins have aligned within the clump to lower the energy, but there is not enough thermal energy in the system to break up neighbouring clumps of opposing spins.  A local minimum of the PES has been found, but the global minimum requires crossing an activation barrier.  A good example of this is the 32x32 lattice above: using exactly the same conditions, the system can either reach the lowest possible energy state or show clumping.

For the effect of temperature, consider the 8x8 lattice shown above.  Well below the critical temperature, there is not sufficient thermal energy to cause large fluctuations in the spin states of the system and spontaneous magnetisation occurs quickly.  Approaching the critical temperature, more fluctuations are possible, and more steps are required for the system to finally reach the lowest possible energy state.  Above the critical temperature, spontaneous magnetisation no longer occurs, and entropic effects start to dominate - from the first step the entire system is fluctuating rapidly.

'''<big>TASK</big>: Use ILtemperaturerange.py to plot the average energy and magnetisation for each temperature, ''with error bars'', for an <math>8\times 8</math> lattice. Use your initution and results from the script ILfinalframe.py to estimate how many cycles each simulation should be. The temperature range 0.25 to 5.0 is sufficient. Use as many temperature points as you feel necessary to illustrate the trend, but do not use a temperature spacing larger than 0.5. T NumPy function savetxt() stores your array of output data on disk &mdash; you will need it later. Save the file as ''8x8.dat'' so that you know which lattice size it came from.'''


The experiment runs for 500000 steps, with the first 250000 steps discarded from the calculation of the average energies.  This is a compromise between allowing the larger lattice systems a better chance of converging to the final state, and running the experiment on a reasonable timescale.  The results from the 8x8 lattice, with error bars:

{| align=center
|-
|[[File:wm1415_New8x8energyplot.png|400px]]
|[[File:wm1415_New8x8magplot.png|400px]]
|}

==The effect of system size==

'''<big>TASK</big>: Repeat the final task of the previous section for the following lattice sizes: 2x2, 4x4, 8x8, 16x16, 32x32. Make sure that you name each datafile that your produce after the corresponding lattice size! Write a Python script to make a plot showing the energy ''per spin'' versus temperature for each of your lattice sizes. Hint: the NumPy loadtxt function is the reverse of the savetxt function, and can be used to read your previously saved files into the script. Repeat this for the magnetisation. As before, use the plot controls to save your a PNG image of your plot and attach this to the report. How big a lattice do you think is big enough to capture the long range fluctuations?'''

{| class="wikitable"
|-
!Lattice Size
!Energy plot
!Magnetisation plot
|-
|2x2
|[[File:wm1415_New2x2energyplot.png|400px|center]]
|[[File:wm1415_New2x2magplot.png|400px|center]]
|-
|4x4
|[[File:wm1415_New4x4energyplot.png|400px|center]]
|[[File:wm1415_New4x4magplot.png|400px|center]]
|-
|8x8
|[[File:wm1415_New8x8energyplot.png|400px|center]]
|[[File:wm1415_New8x8magplot.png|400px|center]]
|-
|16x16
|[[File:wm1415_New16x16energyplot.png|400px|center]]
|[[File:wm1415_New16x16magplot.png|400px|center]]
|-
|32x32
|[[File:wm1415_New32x32energyplot.png|400px|center]]
|[[File:wm1415_New32x32magplot.png|400px|center]]
|-
|64x64
|[[File:wm1415_New64x64energyplot.png|400px|center]]
|[[File:wm1415_New64x64magplot.png|400px|center]]
|-
|All
|[[File:wm1415_Newener_comp.png|400px|center]]
|[[File:wm1415_Newmag_comp.png|400px|center]]
|}

For each finite lattice size plot, there is a noticeable trend in the size of the error bars.  These represent the number of energy and magnetisation states available at each temperature.  Far below <math>T_C</math>, the error bars are almost negligible: there is spontaneous magnetisation and there is not sufficient thermal energy in the system to overcome the activation barrier in flipping the spins.  From around T = 1.5, more states become accessible to the system and the error bars steadily increase in size until the critical temperature is reached, where the system undergoes a secondary phase transition and the error bars are at their largest: the states available at this temperature span the largest range of energies and magnetisations.  At <math> T > T_C</math>, the system has enough thermal energy to undergo fluctuations and the error bars remain considerably large.

As the size of the lattice increases, the size of the error bars decrease.  For a 2x2 lattice, a single spin flip represents a huge change in the energy and magnetisation of the whole system.  For a 64x64 lattice, the difference in energy and magnetisation for the whole system when a single spin flips is almost negligible.  The fluctuations result in smaller fractional changes in the system as the lattice size increases, hence the error bars decrease.

When plotted on the same graph, the energy and magnetisation plots appear to converge for the higher lattice sizes.  The 16x16 and higher lattices are big enough to capture the long-range fluctuations: a single spin flip brings about a finer change in the whole system's energy.

==Determining the heat capacity==
'''<big>TASK</big>: By definition,'''

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

'''From this, show that'''

<math>C = \frac{\mathrm{Var}[E]}{k_B T^2}</math>

(Where <math>\mathrm{Var}[E]</math> is the variance in <math>E</math>.)

The average energy can be related to the partition function<sup>[3]</sup>, <math>Z=\sum\exp\left(-\beta\varepsilon_{i}\right)</math>, with <math>\beta=\frac{1}{k_{B}T}</math> and <math>\varepsilon_{i}</math> being the energy of each state <math>i</math>:


<center><math>\langle E\rangle=\frac{E}{N}=\frac{1}{N}\sum N_{i}\varepsilon_{i}=\frac{1}{Z}\sum\varepsilon_{i}\exp\left(-\beta\varepsilon_{i}\right)=-\frac{1}{Z}\frac{\partial Z}{\partial\beta}=-\frac{\partial\ln Z}{\partial\beta}</math></center>


Similarly, the expectation value for <math>\langle E^{2}\rangle</math> can be calculated from the partition function:


<center><math>\langle E^{2}\rangle=\sum\varepsilon_{i}^{2}p_{i}=\frac{\sum\varepsilon_{i}^{2}\exp\left(-\beta\varepsilon_{i}\right)}{\sum\exp\left(-\beta\varepsilon_{i}\right)}=\frac{1}{Z}\frac{\partial^{2}}{\partial\beta^{2}}\left[\sum\exp\left(-\beta\varepsilon_{i}\right)\right]=\frac{1}{Z}\frac{\partial^{2}Z}{\partial\beta^{2}}</math></center>


Variance, <math>Var\left[E\right]</math>, can be expressed as:


<center><math>Var\left[E\right]=\langle E^{2}\rangle-\langle E\rangle^{2}</math></center>


<center><math>Var\left[E\right]=\frac{1}{Z}\frac{\partial^{2}Z}{\partial\beta^{2}}-\left(-\frac{1}{Z}\frac{\partial Z}{\partial\beta}\right)^{2}</math></center>


Let us show that this expression for variance is equivalent to <math>\frac{\partial}{\partial\beta}\left(\frac{1}{Z}\frac{\partial Z}{\partial\beta}\right)</math>:


<center><math>\frac{\partial}{\partial\beta}\left(\frac{1}{Z}\frac{\partial Z}{\partial\beta}\right)=\frac{1}{Z}\frac{\partial^{2}Z}{\partial\beta^{2}}+\frac{\partial\left(\frac{1}{Z}\right)}{\partial\beta}\frac{\partial Z}{\partial\beta}=\frac{1}{Z}\frac{\partial^{2}Z}{\partial\beta^{2}}+\left(\frac{\partial\left(\frac{1}{Z}\right)}{\partial Z}\frac{\partial Z}{\partial\beta}\right)\frac{\partial Z}{\partial\beta}=\frac{1}{Z}\frac{\partial^{2}Z}{\partial\beta^{2}}-\frac{1}{Z^{2}}\left(\frac{\partial Z}{\partial\beta}\right)^{2}=Var\left[E\right]</math></center>


Hence:

<center><math>Var\left[E\right]=\frac{\partial}{\partial\beta}\left(\frac{1}{Z}\frac{\partial Z}{\partial\beta}\right)=\frac{\partial^{2}\ln Z}{\partial\beta^{2}}</math></center>


But, since <math>\langle E\rangle=-\frac{\partial\ln Z}{\partial\beta}</math>, and <math>\frac{\partial}{\partial\beta}=\frac{\partial T}{\partial\beta}\frac{\partial}{\partial T}=\left(\frac{1}{\frac{\partial\beta}{\partial T}}\right)\frac{\partial}{\partial T}=\left(-k_{B}T^{2}\right)\frac{\partial}{\partial T}</math>, the above can be rewritten as<sup>[3]</sup>:


<center><math>Var\left[E\right]=\frac{\partial^{2}\ln Z}{\partial\beta^{2}}=-\frac{\partial}{\partial\beta}\left(-\frac{\partial\ln Z}{\partial\beta}\right)=-\frac{\partial}{\partial\beta}\langle E\rangle=-\left(-k_{B}T^{2}\frac{\partial}{\partial T}\right)\langle E\rangle=k_{B}T^{2}C_{V}</math></center>


Therefore:

<center><math>C_{V}=\frac{Var\left[E\right]}{k_{B}T^{2}}</math></center>

'''<big>TASK</big>: Write a Python script to make a plot showing the heat capacity versus temperature for each of your lattice sizes from the previous section. You may need to do some research to recall the connection between the variance of a variable, <math>\mathrm{Var}[X]</math>, the mean of its square <math>\left\langle X^2\right\rangle</math>, and its squared mean <math>\left\langle X\right\rangle^2</math>. You may find that the data around the peak is very noisy &mdash; this is normal, and is a result of being in the critical region. As before, use the plot controls to save your a PNG image of your plot and attach this to the report.'''

{| class="wikitable"
|-
!2x2 Heat Capacity
!4x4 Heat Capacity
!8x8 Heat Capacity
!16x16 Heat Capacity
!32x32 Heat Capacity
!64x64 Heat Capacity
!Combined
|-
|[[File:wm1415_Newheatcap_2x2.png|200px|center]]
|[[File:wm1415_Newheatcap_4x4.png|200px|center]]
|[[File:wm1415_Newheatcap_8x8.png|200px|center]]
|[[File:wm1415_Newheatcap_16x16.png|200px|center]]
|[[File:wm1415_Newheatcap_32x32.png|200px|center]]
|[[File:wm1415_Newheatcap_64x64.png|200px|center]]
|[[File:wm1415_Newheatcap_comp.png|200px|center]]
|}

These plots show the heat capacities for each lattice size, using the same data from the previous section.  In all cases, the regions <math>T << T_C</math> and <math>T >> T_C</math> are well-resolved, and in the next section it will be shown that they map to the reference C++ data perfectly.  However, in the region of the critical temperature, larger lattices give a very noisy plot.  Further repeats of larger lattices at the critical temperature will be explored later in the report.

The last plot, showing the heat capacity curves for all lattice sizes, somewhat demonstrates an equation stated in the lab script, relating the temperature at which heat capacity is maximised with lattice size:

<center><math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math></center>

From this, one would expect the peak in the heat capacity to shift to lower temperatures as the lattice size increases. This can be seen up to the 8x8 lattice: above this, only one repeat of the experiment is not sufficient to accurately map the critical region and the trend is not clear for this example.

==Locating the Curie temperature==

'''<big>TASK</big>: A C++ program has been used to run some much longer simulations than would be possible on the college computers in Python. You can view its source code [https://github.com/niallj/ducking-avenger/tree/master/Ising here] if you are interested. Each file contains six columns: <math>T, E, E^2, M, M^2, C</math> (the final five quantities are per spin), and you can read them with the NumPy loadtxt function as before. For each lattice size, plot the C++ data against your data. For ''one'' lattice size, save a PNG of this comparison and add it to your report &mdash; add a legend to the graph to label which is which. To do this, you will need to pass the label="..." keyword to the plot function, then call the legend() function of the axis object (documentation [http://matplotlib.org/api/axes_api.html#matplotlib.axes.Axes.legend here]).'''


{| class="wikitable"
|-
!Lattice Size
!Energy: Python vs C++
!Magnetisation: Python vs C++
!Heat Capacity: Python vs C++
|-
|2x2
|[[File:wm1415_New2_energy_ref_comp.png|400px|center]]
|[[File:wm1415_New2_magnetisation_ref_comp.png|400px|center]]
|[[File:wm1415_New2_capacity_ref_comp.png|400px|center]]
|-
|4x4
|[[File:wm1415_New4_energy_ref_comp.png|400px|center]]
|[[File:wm1415_New4_magnetisation_ref_comp.png|400px|center]]
|[[File:wm1415_New4_capacity_ref_comp.png|400px|center]]
|-
|8x8
|[[File:wm1415_New8_energy_ref_comp.png|400px|center]]
|[[File:wm1415_New8_magnetisation_ref_comp.png|400px|center]]
|[[File:wm1415_New8_capacity_ref_comp.png|400px|center]]
|-
|16x16
|[[File:wm1415_New16_energy_ref_comp.png|400px|center]]
|[[File:wm1415_New16_magnetisation_ref_comp.png|400px|center]]
|[[File:wm1415_New16_capacity_ref_comp.png|400px|center]]
|-
|32x32
|[[File:wm1415_New32_energy_ref_comp.png|400px|center]]
|[[File:wm1415_New32_magnetisation_ref_comp.png|400px|center]]
|[[File:wm1415_New32_capacity_ref_comp.png|400px|center]]
|-
|64x64
|[[File:wm1415_New64_energy_ref_comp.png|400px]]
|[[File:wm1415_New64_magnetisation_ref_comp.png|400px|center]]
|[[File:wm1415_New64_capacity_ref_comp.png|400px]]
|}

For all lattice sizes, the energy and magnetisation are in agreement with the reference C++ data.  The heat capacity plots also match the reference for smaller lattice sizes, producing a smooth curve with almost no noise: the system is small enough to probe all the possible configurations.  However, there is substantial noise in the heat capacity for larger lattices (16x16 and above) in the region of the critical temperature: one repeat of the experiment is not sufficient to accurately model the many possible mechanisms of phase transition.

The relaxation process for a random lattice configuration can be described by a free energy surface.  The path whereby the system relaxes depends on the initial configuration - which, in our model, is random.  Thus, for each run of the experiment, a different section of the free energy surface is explored: the local minimum where the system stabilises will have a different energy in the different paths taken.  The activation barriers between local minima are mediated by 'clumping': regions of spin-aligned lattice sites form, which are unfavourable to disrupt. In the critical temperature region, there is not sufficient thermal energy to eliminate clumping: as a result, there are many metastable states the system can relax to, each with different energies.  This results in a large number of possible configurations, and hence a noisy plot of energy, magnetisation and heat capacity.

To reduce the noise, more repeats must be performed and the results averaged to sample more relaxation paths on the free energy surface.  The figures below show the results of further simulations run in the critical temperature range for larger lattices, with a comparison between a single run of the experiment and an averaged plot from 10 repeats.  Error bars in the heat capacity plot were obtained by calculating the standard deviation across repeats.

{| class="wikitable"
|-
!Lattice Size
!Heat Capacity: Single run
!Heat Capacity: 10 repeats
|-
|16x16
|[[File:wm1415_16x16_heatcap_1run.png|400px]]
|[[File:wm1415_16x16_heatcap_rpts.png|400px]]
|-
|32x32
|[[File:wm1415_32x32_heatcap_1run.png|400px]]
|[[File:wm1415_32x32_heatcap_rpts.png|400px]]
|-
|64x64
|[[File:wm1415_64x64_heatcap_1run.png|400px]]
|[[File:wm1415_64x64_heatcap_rpts.png|400px]]
|}

Even with these repeat experiments, the values for the heat capacity for the larger 32x32 and 64x64 lattices seem to be consistently smaller than the reference data at the critical temperature.  Perhaps this is because the system is not being given enough steps to fully explore the potential energy surface.  All these plots have been using the same method for this report - discounting the first 250,000 steps and recording the next 250,000.  To test this, an additional run of the 32x32 lattice was performed: discounting the first 250,000 steps and recording the next 1,250,000 steps:

<center>[[File:wm1415_Large_32x32_heatcap_1run.png|400px]]</center>

As this was computationally very expensive, only one repeat was manageable for the timescale of this experiment.  While it may not be scientifically rigorous to draw any conclusions from one repeat of an experiment which depends on chance, this figure is certainly promising - the extra steps seem to have alleviated the underestimation of the heat capacity - and a closer match to the C++ data may be obtained with further repeats.  Unfortunately, this was not investigated further due to time constraints of the experiment.

To extract the Curie temperature, the reference data - with clearly resolved peaks in the heat capacity - will be examined first, and then compared to the value obtained from the data generated in Python.

===Locating the Curie temperature - Reference C++ data===

'''<big>TASK</big>: write a script to read the data from a particular file, and plot C vs T, as well as a fitted polynomial. Try changing the degree of the polynomial to improve the fit &mdash; in general, it might be difficult to get a good fit! Attach a PNG of an example fit to your report.'''

<center>[[File:wm1415_Initial_polyfits.png|400px|thumb|In general, a higher order polyfit, which includes higher order terms, gives the best fit.  However, without restricting the range of the fit, the peak cannot be located with this approach.]]</center>

'''<big>TASK</big>: Modify your script from the previous section. You should still plot the whole temperature range, but fit the polynomial only to the peak of the heat capacity! You should find it easier to get a good fit when restricted to this region.'''

{| align=center
|-
|[[File:wm1415_Polyfit_max.png|400px]]
|[[File:wm1415_Polyfit_max_zoom.png|400px]]
|}

For the remainder of the report, a polyfit order of 8 has been chosen: high enough to obtain a good fit, but sufficiently low to avoid poorly optimised high-order term fits.

'''<big>TASK</big>: find the temperature at which the maximum in C occurs for each datafile that you were given. Make a text file containing two colums: the lattice side length (2,4,8, etc.), and the temperature at which C is a maximum. This is your estimate of <math>T_C</math> for that side length. Make a plot that uses the scaling relation given above to determine <math>T_{C,\infty}</math>. By doing a little research online, you should be able to find the theoretical exact Curie temperature for the infinite 2D Ising lattice. How does your value compare to this? Are you surprised by how good/bad the agreement is? Attach a PNG of this final graph to your report, and discuss briefly what you think the major sources of error are in your estimate.'''
{| class="wikitable"
|-
!2x2 Polyfit
!4x4 Polyfit
!8x8 Polyfit
!16x16 Polyfit
!32x32 Polyfit
!64x64 Polyfit
|-
|[[File:wm1415_Polyfit_2x2.png|200px|center]]
|[[File:wm1415_Polyfit_4x4.png|200px|center]]
|[[File:wm1415_Polyfit_8x8.png|200px|center]]
|[[File:wm1415_Polyfit_16x16.png|200px|center]]
|[[File:wm1415_Polyfit_32x32.png|200px|center]]
|[[File:wm1415_Polyfit_64x64.png|200px]]
|}

<center>[[File:wm1415_Curie_t_inf.png|500px]]
<math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math></center>

A plot of <math>T_{C, L}</math> against <math>\frac{1}{L}</math> yields a straight line, with intercept <math>T_{C, \infty}</math> and gradient <math>A</math> (a constant).  Fitting the C++ data resulted in a value for <math>T_{C,\infty} = 2.277 \pm 0.017 \frac{J}{k_B}</math>, with the uncertainty calculated from the error in the linear fit.  This is in good agreement with the reference value<sup>[1]</sup> of <math>T_{C,\infty} \thickapprox 2.269 \frac{J}{k_B}</math>.  The main source of error is the fitting of the straight line, which is quantifiable and represents a 1% uncertainty.  There is also an uncertainty in the extraction of <math>T_{C,L}</math> values from the polynomial fitting of the peaks in each heat capacity plot, as the fits are not exact.  The simulations themselves also introduce an uncertainty: for example, if the lattice did not equilibrate within the limit of the number of simulation steps.

However, this experiment has reproduced the theoretical value within the error of the straight line fit, demonstrating the power of the Monte Carlo methods used.

===Locating the Curie temperature - Python data===

Below are the fitted functions for the Python heat capacity data.  For the 16x16, 32x32 and 64x64 lattices, the averaged data points over 10 repeats were used (see above).

{| class="wikitable"
|-
!2x2 Polyfit
!4x4 Polyfit
!8x8 Polyfit
!16x16 Polyfit
!32x32 Polyfit
!64x64 Polyfit
|-
|[[File:wm1415_Pypolyfit_2x2.png|200px|center]]
|[[File:wm1415_Pypolyfit_4x4.png|200px|center]]
|[[File:wm1415_Pypolyfit_8x8.png|200px|center]]
|[[File:wm1415_Pypolyfit_16x16.png|200px|center]]
|[[File:wm1415_Pypolyfit_32x32.png|200px|center]]
|[[File:wm1415_Pypolyfit_64x64.png|200px|center]]
|}

Compiling this data to plot the straight line, as outlined previously, yielded very similar results to the C++ data:

<center>[[File:wm1415_PyCurie_t_inf.png|400px]]</center>

There is a larger uncertainty in the fitting of the straight line but the resulting curie temperature, <math>T_{C,\infty} = 2.28 \pm 0.03 \frac{J}{k_B}</math>, is nearly identical to the result from the C++ data and is in good agreement with the reference value<sup>[1]</sup>.  Due to the noise in the heat capacity plots in the critical temperature region, the polynomial fits of the peaks were not as convincing as the fits for the C++ data.  The <math>T_{C,L}</math> values that were subsequently extracted from these polynomial fits were shifted from the C++ data, resulting in the greater spread of data points in the straight line plot above.  Furthermore, there are contributing errors from the simulations themselves - for example, the lattice not reaching equilibrium or getting stuck in a local minimum (the 'clumping' effect discussed earlier) affecting the heat capacity.

==References==

1 B. Liu, M. Gitterman, ''American Journal of Physics'', '''71''', 806 (2003), pp. 1-4

2 R. Fitzpatrick, 2006, The Ising Model, ''Computational Physics'', retrieved from http://farside.ph.utexas.edu/teaching/329/lectures/node110.html in November 2017.

3 P. Atkins and J. de Paula, ''Atkins' Physical Chemistry'', Oxford University Press, UK, 8th edn, 2006, pp. 564-573