Rep:JJ1516Ising

2018-11-21T09:53:39Z

Jj1516: /* Finding the Peak C Values for Each Lattice Size */

Ex Script Link:
https://wiki.ch.ic.ac.uk/wiki/index.php?title=Third_year_CMP_compulsory_experiment

= Introduction to the Ising Model =

== The Model ==

=== Minimum Energy State ===
====Minimum Interaction Energy====
The interaction energy of an Ising lattice is given by:

<math>E = -\frac{1}{2}J\displaystyle\sum_{i}^{N} \displaystyle\sum_{j\;\in\;{neighbours}\;(i)} \!\!\!\!\!\!\!\!\!\!\!\!s_{i}s_{j}</math>

In an Ising lattice with D dimensions and periodic boundary conditions, a given spin <math> s_i</math> has 2D neighbours.

The minimum interaction energy is achieved when all spins in the lattice are aligned (all ↑ or all ↓). As such, any spin <math> s_i</math> must be aligned with all 2D of its neighbours. That is, <math>s_{i}s_{j} = 1</math> for all <math>j</math>.

Therefore, for any <math> s_i</math>, the sum of its interactions is <math>\displaystyle\sum_{j\;\in\;{neighbours}\;(i)} \!\!\!\!\!\!\!\!\!\!\!\!s_{i}s_{j} = 2D </math>

As this is true for all N spins in the lattice, it follows that <math>\displaystyle\sum_{i}^{N} 2D = 2DN</math>.

Thus the minimum interaction energy is given by:

:<math>E_{min} = -\frac{1}{2}J\cdot2DN = - DNJ </math>

==== Entropy of the Minimum Energy State====

The entropy of a state is given by:

:<math>S = k\; ln\; \Omega </math>

Where <math> \Omega </math> is the multiplicity of the state, the number of degenerate configurations of said state. The minimum energy state is achievable with all spins ↑ or with all spins ↓. In the absence of an external magnetic field, these configurations are degenerate and indistinguishable giving <math> \Omega = 2 </math> for this state. Therefore,

<math> S_{E_{min}} = k_{B}\; ln\; 2 = 0.693k_{B} = 9.57 \times 10^{-24} J\!\cdot\!\! K^{-1}</math>

== Phase Transitions ==

=== Changes to Minimum Energy Configuration ===

====Energy Change ====
The energy change <math> \Delta E</math> associated with a single spin 'flip' from the lowest energy onfiguration may be calculated by considering the change in the interaction energy with the spin's neighbours.

Let <math> D = 3, N = 1000 </math>

<math>\Delta{E} = E_1 - E_0</math>

<math>E_1</math> Interaction energy after one spin has been 'flipped'

<math>E_0</math> Lowest energy configuration

[[Image:Task 2 table.PNG]]

<math>\displaystyle\sum_{interactions} = -2D + 2D(2D-2) + 2D[N-(2D+1)] </math>

<math> = 2DN - 8D </math>

<math>E_1 = \displaystyle-\frac{1}{2}J\;\sum_{interactions} = 4DJ - DNJ</math>

<math>\Delta{E} = E_1 - E_0 = 4DJ - DNJ + DNJ = +4DJ = +12J</math>

This energy change is independent of <math> N</math> and depends only on the dimension size and the coupling constant <math> J </math>.

In this calculation, the following assumption has been made.

* The lattice dimensions are infinite relative to the size of each cell

This allows '''any''' cell to be treated as the centre of the lattice. As every cell interacts with 2D neighbours, flipping any cell in the minimum energy configuration results in the same energy change. Alternatively, the periodic boundary conditions of the Ising model remove the necessity for such an assumption. As a consequence, the absence of <math> N </math> from the resulting energy change is justified.

====Entropy Change====

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

<math>\Omega = \displaystyle\frac{N!}{{N_{\uparrow}!}{N_{\downarrow}!}}</math>

Where <math>N_{x}!</math> denotes the factorial of the number of spins in spin state <math>x</math>.

<math>\Delta{S} = k_{B}\;ln\;\left(\displaystyle\frac{\Omega_{2}}{\Omega_{1}}\right)</math>

<math>\Omega_{2} = \frac{1000!}{999!1!} = 1000</math>

<math> \Omega_{1} = \frac{1000!}{1000!} = 1</math>

<math>\Delta{S} = k_{B}\;ln\;1000 = 9.54 \times 10^{-23} JK^{-1}</math>

<math>\Delta{S_{mol}} = N_{A}k_{B}\;ln\;1000 = 57.4 \,JK^{-1}mol^{-1}</math>

===Magnetisation Values of 1D and 2D Lattices===

<math>M = \sum_i s_i</math>

====1D Lattice Magnetisation====

<math>M_{1D} =\;\;+1</math>

====2D Lattice Magnetisation====

<math>M_{2D} =\;\;+1</math>

====Magnetisation at 0K====

At 0K, no thermal energy is available to overcome the energetic barriers to spin flipping. Considering the Hemholtz free energy:

<math> A = E - TS </math>

A temperature of 0K nullifies the entropic contribution to the free energy. Thus the energetic driving force (which maximises parallel spin interactions) dominates the entropic driving force (which maximises multiplicity of lattice states). For an Ising Lattice with D = 3, N = 1000 at 0K, the alignment of all spins results in the minimum energy configuration and an overall magnetisation of <math> M = \pm 1000 </math>.

= Calculating the Energy and Magnetisation =
==Modifying the Files==

===Energy and Magnetisation Functions===

====Energy===
[[File:First_Ene_Method.PNG]]

The above energy method segments a given lattice into lattice cells on the edges, corners or in the bulk of the lattice and calculates the energy interactions with its neighbours as appropriate.

====Magnetisation====
[[File:First_Mag_Method.PNG]]

The magnetisation method loops over all elements in the array representing the Ising lattice and adds the value of the spin to the magnetisation variable.

==Testing the Files==

===ILCheck.py===

[[File:ILChecks.PNG]]

Above, the results of the ILcheck.py script are presented.

= Introduction to the Monte Carlo Simulation =
==Average Energy and Magnetisation==
===Configurations of a 100 Spin System===
For a system of two spins, 4 configurations are possible (↑↑, ↓↓, ↑↓, ↓↑) i.e. <math>2^{2}</math>. For a system of 100 spins, there are <math>2^{100}</math> configurations. In the first example, the ↑↑ and ↓↓ are degenerate, as are the ↑↓ and ↓↑ configurations. Therefore, dividing the number of possible configurations by the multiplicity of each configuration (<math>2</math>) gives the number of non-degenerate configurations. As previously stated, degenerate configurations are indistinguishable in the absence of an external magnetic field.

For 100 spins, this gives <math>\displaystyle\frac{2^{100}}{2} = 6.338 \times 10^{29} </math> distinguishable configurations.

If processing speed permits <math> 1 \times 10^{9}</math> calculations each second, the evaluation of a single value of <math>\left\langle M\right\rangle_T</math> would last <math> 6.338 \times 10^{20} </math> seconds or <math> 1458 </math> times the age of the Universe!

==Modifying IsingLattice.py==

===Monte Carlo and Statistics Methods===
[[File:Monte_Carlo.PNG]]

[[File:stats.PNG]]

The Metropolis Monte Carlo (MMC) algorithm provides a computationally inexpensive path from random initial states to the minimum free energy states of the Ising lattices generated. First, a random 2D Ising lattice is generated. A single random spin within this lattice is "flipped" and the resulting (internal) energy change between these states is measured. Recalling the relationship between the Hemholtz free energy, internal energy and entropy, it follows that a decrease in internal energy also lowers the free energy of the system: thus the configuration change is accepted.

In the case that the spin "flip" does not lower the internal energy, the MMC algorithm tests this configuration change against the Boltzman factor. In the canonical (NVT) ensemble, the free energy of the system is given by:

:<math>A = -k_{B}T\,lnZ_{NVT} </math>

where Z is the canonical partition function which depends upon <math> e^{-\frac{\Delta E}{k_BT}} </math>, the Boltzman factor. By choosing a random value in the interval [0,1)from a continuous uniformm probability distribution and rejecting any value of the Boltzman factor less than this value, ###

===Spontaneous Magnetisation===

Below the Curie Temperature, <math> T_{C} </math>, the energetic driving force dominates the entropic driving force, resulting in a highly ordered system and an overall magnetisation. The spins of the system align in order to maximise favourable coupling interactions rather than adopting a less ordered configuration of higher multiplicity. Therefore, when <math>T < T_C</math>, the most probable configurations have an overall magnetisation resulting in <math>\left\langle M\right\rangle \neq 0</math>.

The plots of magnetisation per spin and energy per spin obtained from ILanim.py support this prediction. Below, the plots show that upon equilibration, the system adopts the minimum energy per spin: <math>E_{min_{spin}} = \frac{-NDJ}{N} = -DJ = -2</math> and exhibits the maximum magnitude of magnetisation per spin: <math>|M_{spin}| = 1</math>

[[File:2T=1.0.png]]
[[File:T=1.0.png]]
[[File:3T=1.0.png]]

An unknown system error prevented the statistics() function from printing the average values.

=Accelerating the Code=

===Initial Method Speeds===
Taking 10 repeat measurements for the time taken to perform 2000 Monte Carlo steps resulted in

* An average time of <math>11.7s</math>
* A standard error of <math>0.020s</math>

=== Updating Measurement===

NOTE: ILtimetrial.py has been modified to:
* Perform a requested number of '2000 Monte Carlo Step' measurements
* Record these measurements in a CSV file
* Return the mean measuremd time and the Standard Error i.e. <math> \frac{\sigma}{n} </math>

[[File:TimeTrialMethod.PNG]]

===Updating Energy and Magnetisation Methods===
[[File:New_Mag_Method.PNG]]
[[File:New_Ene_Method.PNG]]

===Updated Method Speeds===
[[File:TimeTrial.PNG]]

Taking 30 measurements of the time taken to perform 2000 Monte Carlo Steps resulted in:

* An average time of <math>\displaystyle{0.87s}</math>
* A standard error of <math>2.5 \times 10^{-3}s</math>

This represents a 93% decrease in the average computation time from the initial methods.

=The Effect of Temperature=

==Correcting the Averaging Code==

===Finding a Delay Value===
To quickly find a value for <math> N_{d}</math>, the minimum number of cycles typically needed for the equilibration of a system of within the lattice and temperature ranges of interest, it was assumed that:

* Larger lattices would require more Monte Carlo steps to reach equilibrium from their initial, random configurations

* Lattices at lower temperatures would require more Monte Carlo steps to reach equilibrium from their initial, random configurations.

Therefore, an initial value for <math> N_{d}</math> was found for a <math> 128\times128 </math> lattice at a temperature of <math>0.5</math> reduced temperature units.
*<math> N_{d} = 185000 </math>

By this number of steps, all smaller lattices at all higher temperatures were found to have equilibrated also. Thus, this value of<math> N_{d}</math> was accepted as an appropriate delay.

[[File:AFinFrame 128 3e5 T= 0.5line.png|400px]]

For the 128x128 lattices at T=0.5, equilibration was ensured from 185000 steps onwards.

Smaller lattices at the same temperature did indeed equilibrate more slowly, but all equilibrated within <math> N_{d}</math> MC steps.

[[File:AFinFrame 64 3e5 T= 0.5line.png|400px]]
[[File:AFinFrame 32 3e5 T= 0.5line2.png|400px]]
[[File:AFinFrame 16 3e5 T= 0.5line2.png|400px]]
[[File:AFinFrame 8 3e5 T= 0.5line.png|400px]]

The figure below shows the modifications made to IsingLattice.py to incorporate the <math> N_{d} </math> step delay.

==Running over a Range of Temperatures==

To examine the effect of temperature on the average energy <math> \left\langle E\right\rangle </math> and average magnetisation <math> \left\langle M\right\rangle </math> of Ising lattices at equilibrium, the following parameters were used:

* Temperature Range: <math> 0.5 - 5 </math> reduced units
* Temperature Interval: <math>0.01</math> reduced units (<math>450</math> temperature values)
* Number of MMC Steps per Temperature value: <math> 2 \times 10^5 </math> with the aforementioned value of <math> N_{d} = 1.85 \times 10^5 </math> Monte Carlo steps
* Lattice Sizes: <math> 2 \times 2 </math>, <math> 4 \times 4 </math> , <math> 8 \times 8 </math> , <math> 16 \times 16 </math> , <math> 32 \times 32 </math> , <math> 64 \times 64 </math>

Plots of <math> \left\langle E\right\rangle </math> and <math> \left\langle M\right\rangle </math> per spin as versus temperature are presented below for each lattice size. Errorbars have been implemented to demonstrate how ''well-behaved'' lattice systems of each size are at each temperature, that is, how temperature and lattice size affect the size of the variations in the expected values. For the error calculations, the standard errors, <math> \sigma_\bar{X} </math>, in <math> \left\langle E\right\rangle </math> and <math> \left\langle M\right\rangle </math> per spin have been used:

====8x8====

[[File:8x8e+error.png]] [[File:8x8m+error.png]]

A Python script has been written to generate a plot of <math>\left\langle E\right\rangle</math> or <math>\left\langle M\right\rangle</math> against temperature based upon user input. A screenshot of this code has been inserted below.

[[File:ExpectedvsTemp.PNG|200px]]

=The Effect of System Size=

==Scaling the System Size==

===More Lattice Sizes===

====2x2====
[[File:2x2e+error.png]] [[File:2x2m+error.png]]

====4x4====
[[File:4x4e+error.png]] [[File:4x4m+error.png]]

====16x16====
[[File:16x16e+error.png]] [[File:16x16m+error.png]]

====32x32====
[[File:32x32e+error.png]] [[File:32x32m+error.png]]

====64x64====
[[File:64x64e+error.png]] [[File:64x64m+error.png]]

Scrutiny of the plots of average energy per spin for each lattice size reveals an increase in the errorbar width with increasing lattice size. This is most notable in the low temperature region of the plots. Between the <math> 8\times8 </math> and <math> 16\times16 </math> lattices, the long range fluctuations become suddenly prominent. The <math>16\times16</math> lattice is the first (in the range of sizes tested)to clearly capture these fluctuations.

=Determining the Heat Capacity=

==Calculating the Heat Capacity==

===Heat Capacity Derivation===
Working in the canonical ensemble (constant volume), the heat capacity <math> C_V </math> is given by:

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

where <math>\left\langle E\right\rangle</math> is the average total energy or expected value of the internal energy of the system.

: <math> \left\langle E\right\rangle = - \frac{\partial ln Q}{\partial \beta} = -\frac{1}{Q}\frac{\partial Q}{\partial \beta} </math>
where <math> Q</math> is the partition function for the system of N interacting spins.

By evaluating the second form of the partial derivative with respect to <math> \beta</math>:

<math> \left\langle E\right\rangle = -\frac{1}{Q}\frac{\partial Q}{\partial \beta} = \frac{\Sigma \epsilon_{i} e^{-\beta \epsilon_{i}}}{\Sigma e^{-\beta \epsilon_{i}}} </math>

Applying the product rule enables the differentiation of this expression for <math>\left\langle E\right\rangle</math> with respect to temperature, giving:

<math>\frac{\partial \left\langle E\right\rangle}{\partial T} = \frac{\Sigma \epsilon_{i} \frac{\epsilon_i}{k_BT^2}e^{-\beta \epsilon_{i}}}{\Sigma e^{-\beta \epsilon_{i}}} - \frac{\Sigma \epsilon_i e^{-\beta \epsilon_i}\cdot \frac{\epsilon_i}{k_BT^2}e^{-\beta \epsilon_i}}{(\Sigma e^{-\beta \epsilon_{i}})^2}</math>

<math>\frac{\partial \left\langle E\right\rangle}{\partial T} = \frac{1}{k_BT^2}\left[\frac{\Sigma \epsilon_i^2 e^{-\beta \epsilon_i}}{\Sigma e^{-\beta \epsilon_i}} - \left(\frac{\Sigma \epsilon_{i} e^{-\beta \epsilon_{i}}}{\Sigma e^{-\beta \epsilon_{i}}}\right)^2\right] </math>

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

Thus,

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

===Heat Capacity versus Temperature===
Below, a screenshot of the Python script written to produce plots of <math> C_V </math> against <math> T </math>, together with the output images of this script.

[[File:HeatCapScript.PNG]]

====2x2====
[[File:HC2.png]]

====4x4====
[[File:HC4.png]]

====8x8====
[[File:HC8.png]]

====16x16====
[[File:HC16.png]]

====32x32====
[[File:HC32.png]]

= Locating the Curie Temperature =

===C++ Data Comparison===
Below a comparison of the <math> C_V </math> against <math> T </math> plot generated for the <math> 2 /times 2</math> to the corresponding C++ data is offered.

[[File:CHC2.png]]

The data generated by the ILtemperaturerange.py script is a faithful match to the C++ data. However, this is an exceptional case; the plots for all other lattice sizes were found to be very dissimilar, necessitating polynomial fitting.

==Polyfitting of C++ Data==

===Fitting a Polynomial to the Entire Range===
The following Python script was written to perform a polynomial fit to the entire range of the graph of <math> C_V </math> against <math> T </math>.

[[File:HeatCapPoly.PNG]]

Below, the result of the polyfit to the <math> 2/times2</math> lattice data is offered (for comparison to the previous, non-fitted example).

[[File:15Polyc2x2.png]]

Polynomials with 5 terms and above were found to match the data in the vicinity of the critical region considerably well. Below 7 terms, the ''tails'' of the plot were matched poorly. In the above example, a polynomial of 15 terms was fit to the data.

===Fitting a Polynomial to the Critical Region===

The following Python script was written to perform a polynomial fit to the region of interest (the ''critical'' region) of the plots of <math> C_V </math> against <math> T </math> to facilitate the discovery of <math> T_C </math>.

[[File:CapPolyCrit.PNG]]

An example of the result returned by the script is offered below for the C++ data for the <math> 16\times16</math> lattice.

[[File:HHCPC16,2.6-2.1,5.png]]

===Finding the Peak C Values for Each Lattice Size===
Below, a screenshot of the code used to obtain the temperatures at which the maximum values of C occured is presented for each lattice size. The temperatures themselves are also featured on the image.

[[File:Tmaxes.PNG]]

The range of lattice sizes and their associated <math>T_{C} </math> values were exported to a [[Media:Tmax.txt|text file]], enabling a plot to be created.

[[File:CritTempPlot.png]]

Subsequently, the scipy.optimize.curve_fit function was used to fit <math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math> function to the data obtained. This enabled the estimation the constant <math> A </math> and the value of <math>T_{C} </math> at <math> \frac{A}{L} = 0</math>, otherwise <math> T_{C,\infty}</math>.

[[File:FindingTc.PNG]]

A value of <math> T_{C,\infty} = 2.2799 </math> was obtained. Comparing this to the [https://en.wikipedia.org/wiki/Square-lattice_Ising_model#Critical_temperature|literature] value of <math>2.269</math>, reveals an error in the obtained value of only <math> 0.47% </math>. Surprisingly, the investigation of lattices in the given size range was enough to capture the theoretical behaviour of an infinite lattice. As explained by Frenkel and Smit, the implementation of periodic boundary conditions enables the simulation of an infinite bulk (or infinite adjacent lattices) surrounding the single lattice under scrutiny. The agreement of the experimental and theoretical values for <math> T_{C,\infty} </math> prove the effectiveness of this idea.
A likely source of error arises from the polynomial fitting of the data. The peaks in <math> C_V </math> against <math> T </math> were far more sharp than the broad peaks of the polynomial functions used to fit them. In several cases, the polynomial peak was a slight, yet discernible overestimate of the true peak. This has resulted in an overestimation of the final <math> T_{C,\infty} </math> value. Fitting a function with a sharper peak to the critical region could overcome this overestimation and provide a more faithful final value.

== References ==
* Frenkel, Daan, et al. Understanding Molecular Simulation : From Algorithms to Applications, Elsevier Science & Technology, 2001

* Krieger, M. Constitutions of matter : Mathematically modeling the most everyday of physical phenomena. Chicago ; London: University of Chicago Press, 1996

* https://en.wikipedia.org/wiki/Square-lattice_Ising_model Accessed 13-11-18

Rep:JJ1516Ising

2018-11-21T09:52:22Z

2018-11-20T12:05:25Z

Jj1516:

Ex Script Link:
https://wiki.ch.ic.ac.uk/wiki/index.php?title=Third_year_CMP_compulsory_experiment

= Introduction to the Ising Model =

== Calculating the minimum energy and entropy ==
Show that the lowest possible lattice energy under the Ising Model is <math> E = -DNJ </math>

The energy of an Ising lattice is given by: <math> E = -DNJ</math>

In D dimensions, each lattice cell within a lattice with periodic boundary conditions has 2D neighbours.
<math>E = -\frac{1}{2}J\displaystyle\sum_{i}^{N} \displaystyle\sum_{j\;\in\;{neighbours}\;(i)} \!\!\!\!\!\!\!\!\!\!\!\!s_{i}s_{j}</math>

In D dimensions with periodic boudnary conditions, each lattice cell has 2D neighbours.

For the minimum energy interactions, all <math>s_{i}s_{j}</math> pairs are spin aligned such that <math>s_{i}s_{j} = 1</math> for all <math>i</math> and <math>j</math>.

Therefore, <math>\displaystyle\sum_{j\;\in\;{neighbours}\;(i)} \!\!\!\!\!\!\!\!\!\!\!\!s_{i}s_{j} = 2D </math>

It follows that <math>\displaystyle\sum_{i}^{N} 2D = 2DN</math>

<math>E_{min} = -\frac{1}{2}J2DN = - DNJ </math>

The entropy of a given lattice configuration is given by:

<math>S_{min} = k\; ln\; \Omega </math>

For a given lattice energy, <math> \Omega </math> is the total number of configurations which will give this energy. configuration may be achieved. In the case of the minimum energy configuration, which may be achieved with all spins or all spins down, <math> \Omega = 2 </math>. Therefore,

<math> S_{E_{min}} = k_{B}\; ln\; 2 = 0.693k_{B} </math>

== Calculating energy and entropy changes ==
'''What is the energy change from the lowest energy configuration when a single spin "flips"? What is the entropy change? (<math> D = 3, N = 1000 </math>)'''

<math>\Delta{E} = E_1 - E_0</math>

<math>E_1</math> One lattice cell spin flipped

<math>E_0</math> Lowest energy configuration

[[Image:Task 2 table.PNG]]

<math>\displaystyle\sum_{interactions} = -2D + 2D(2D-2) + 2D[N-(2D+1)] </math>

<math> = 2DN - 8D </math>

<math>E_1 = \displaystyle-\frac{1}{2}J\;\sum_{int} = 4DJ - DNJ</math>

<math>\Delta{E} = +4DJ = +12J</math>

Energy difference is independent of the number of lattice cells

In this calculation, the following assumption has been made.

* The lattice dimensions are infinite relative to the size of each cell

This allows '''any''' cell to be treated as the centre of the lattice. In other words, it is possible to shift the lattice such that a cell initially at the lattice boundary is moved to the centre. This is justified by the periodic boundary conditions of the Ising model and in turn justifies this calculation's independence of N.

Entropy Change

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

<math>\Omega = \displaystyle\frac{N!}{{N_{\uparrow}!}{N_{\downarrow}!}}</math>

Where <math>N_{i}!</math> denotes the factorial of the number of spins in spin state <math>i</math>.

<math>\Delta{S} = k_{B}\;ln\;\left(\displaystyle\frac{\Omega_{2}}{\Omega_{1}}\right)</math>

<math>\Omega_{2} = \frac{1000!}{999!1!} = 1000</math>

<math> \Omega_{1} = \frac{1000!}{1000!} = 1</math>

<math>\Delta{S} = k_{B}\;ln\;1000 = 9.54 \times 10^{-23} JK^{-1}</math>

<math>\Delta{S_{mol}} = N_{A}k_{B}\;ln\;1000 = 57.4 \,JK^{-1}mol^{-1}</math>

===Magnetisation Values of 1D and 2D Lattices===

<math>M_{1D} =\;\;+1</math>

<math>M_{2D} =\;\;+1</math>

At 0K, no thermal energy is available to overcome the energetic barriers to spin flipping. Thus the energetic driving forces dominates the entropic driving force; for an Ising Lattice with D = 3, N = 1000 at 0K, the lowest energy configuration is achieved when all spins are parallel. This resutls in the maximum number of coupling interactions between spins.

<math>M_{expected} = 1000</math>

= Calculating the Energy and Magnetisation =

[[File:First_Ene_Method.PNG]]

[[File:First_Mag_Method.PNG]]

[[File:ILChecks.PNG]]

= Introduction to the Monte Carlo Simulation =
===Configurations Available to a 100 Spin System===
For a system of two spins, 4 configurations are possible (↑↑, ↓↓, ↑↓, ↓↑) i.e. <math>2^{2}</math>. For a system of 100 spins, <math>2^{100}</math> configurations are expected. In the first example, the ↑↑ and ↓↓ are degenerate, as are the ↑↓ and ↓↑ configurations. Therefore, dividing the number of possible configurations by the multiplicity of each configuration (<math>2</math>) gives the number of distinguishable configurations.
For 100 spins, this gives <math>\displaystyle\frac{2^{100}}{2} = 6.338 \times 10^{29} </math> distinguishable configurations or states.

If our computational speed allowed for <math> 1 \times 10^{9}</math> calculations each second, the evaluation of a single value of <math>\left\langle M\right\rangle_T</math> would last <math> 6.338 \times 10^{20} </math> seconds or <math> 1458 </math> times the age of the Universe!

===Modifying IsingLattice.py===

[[File:Monte_Carlo.PNG]]

The Metropolis Monte Carlo (MMC) algorithm provides a computationally inexpensive path from random initial states to the minimum free energy states of the Ising lattices generated. First, a random 2D Ising lattice is generated. A single random spin within this lattice is "flipped" and the resulting energy change between these states is measured. If the energy change is negative, this represents a decrease in the internal energy of the system which lowers the free energy: this configuration change is accepted.

In the case that the spin "flip" does not lower the free energy, the MMC algorithm tests this configuration change against the Boltzman factor, enabling a determination of the resulting entropy change. In the canonical (NVT) ensemble, the entropy of the system is given by

###

[[File:stats.PNG]]

===Spontaneous Magnetisation===

Below the Curie Temperature, <math> T_{C} </math>, energetic driving forces dominate entropic driving forces. Therefore, the spins of the system align in order to maximise favourable interactions rather than adopting a less ordered configuration, which would maximise multiplicity (and thus entropy).
Below <math>T < T_C</math>, the most probable configurations have an overall magnetisation resulting in <math>\left\langle M\right\rangle \neq 0</math>.

The plots of magnetisation per spin and energy per spin obtained from ILanim.py upon system equilibration show that at minimum energy per spin <math>E_{min_{spin}} = \frac{-NDJ}{N} = -DJ = -2</math>, the magnitude of the overall magnetisation per spin is at a maximum <math>|M_{spin}| = 1</math>

[[File:2T=1.0.png]]
[[File:T=1.0.png]]
[[File:3T=1.0.png]]

==Accelerating the Code==

===Initial Method Speeds===
Taking 10 repeat measurements for the time taken to perform 2000 Monte Carlo steps resulted in

* An average time of <math>11.7s</math>
* A standard error of <math>0.020s</math>

=== Updating Measurement===

NOTE: ILtimetrial.py has been modified to:
* Perform a requested number of '2000 Monte Carlo Step' measurements
* Record these measurements in a CSV file
* Return the mean measuremd time and the Standard Error i.e. <math> \frac{\sigma}{n} </math>

[[File:TimeTrialMethod.PNG]]

===Updating Energy and Magnetisation Methods===
[[File:New_Mag_Method.PNG]]
[[File:New_Ene_Method.PNG]]

===Updated Method Speeds===
[[File:TimeTrial.PNG]]

Taking 30 measurements of the time taken to perform 2000 Monte Carlo Steps resulted in:

* An average time of <math>\displaystyle{0.87s}</math>
* A standard error of <math>2.5 \times 10^{-3}s</math>

This represents a '''INSERT'''% decrease in computation time from the initial methods.

==The Effect of Temperature==

===Correcting the Averaging Code===

To quickly find a value for <math> N_{d}</math>, the minimum number of cycles typically needed for the equilibration of a system of within the lattice and temperature ranges of interest, it was assumed that:

* Larger lattices would require more Monte Carlo steps to reach equilibrium from their initial, random configurations

* Lattices at lower temperatures would require more Monte Carlo steps to reach equilibrium from their initial, random configurations.

Therefore, an initial value for <math> N_{d}</math> was found for a <math> 128\times128 </math> lattice at a temperature of <math>0.5</math> reduced temperature units.
*<math> N_{d} = 185000 </math>

By this number of steps, all smaller lattices at all higher temperatures were found to have equilibrated also. Thus, this value of<math> N_{d}</math> was accepted as an appropriate delay.

[[File:AFinFrame 128 3e5 T= 0.5line.png|400px]]

For the 128x128 lattices at T=0.5, equilibration was ensured from 185000 steps onwards.

Smaller lattices at the same temperature did indeed equilibrate more slowly, but all equilibrated within <math> N_{d}</math> MC steps.

[[File:AFinFrame 64 3e5 T= 0.5line.png|400px]]
[[File:AFinFrame 32 3e5 T= 0.5line2.png|400px]]
[[File:AFinFrame 16 3e5 T= 0.5line2.png|400px]]
[[File:AFinFrame 8 3e5 T= 0.5line.png|400px]]

The figure below shows the modifications made to IsingLattice.py to incorporate the <math> N_{d} </math> step delay.

===Running over a Range of Temperatures===

To examine the effect of temperature on the average energy <math> \left\langle E\right\rangle </math> and average magnetisation <math> \left\langle M\right\rangle </math> of Ising lattices at equilibrium, the following parameters were used:

* Temperature Range: <math> 0.5 - 5 </math> reduced units
* Temperature Interval: <math>0.01</math> reduced units (<math>450</math> temperature values)
* Number of MMC Steps per Temperature value: <math> 2 \times 10^5 </math> with the aforementioned value of <math> N_{d} = 1.85 \times 10^5 </math> Monte Carlo steps
* Lattice Sizes: <math> 2 \times 2 </math>, <math> 4 \times 4 </math> , <math> 8 \times 8 </math> , <math> 16 \times 16 </math> , <math> 32 \times 32 </math> , <math> 64 \times 64 </math>

Plots of <math> \left\langle E\right\rangle </math> and <math> \left\langle M\right\rangle </math> per spin as versus temperature are presented below for each lattice size. Errorbars have been implemented to demonstrate how ''well-behaved'' lattice systems of each size are at each temperature, that is, how temperature and lattice size affect the size of the variations in the expected values. For the error calculations, the standard errors, <math> \sigma_\bar{X} </math>, in <math> \left\langle E\right\rangle </math> and <math> \left\langle M\right\rangle </math> per spin have been used:

====8x8====

[[File:8x8e+error.png]] [[File:8x8m+error.png]]

A Python script has been written to generate a plot of <math>\left\langle E\right\rangle</math> or <math>\left\langle M\right\rangle</math> against temperature based upon user input. A screenshot of this code has been inserted below.

[[File:ExpectedvsTemp.PNG|200px]]

===Scaling the System Size===

====2x2====
[[File:2x2e+error.png]] [[File:2x2m+error.png]]

====4x4====
[[File:4x4e+error.png]] [[File:4x4m+error.png]]

====16x16====
[[File:16x16e+error.png]] [[File:16x16m+error.png]]

====32x32====
[[File:32x32e+error.png]] [[File:32x32m+error.png]]

====64x64====
[[File:64x64e+error.png]] [[File:64x64m+error.png]]

==Determining the Heat Capacity==

===Calculating the Heat Capacity===
Working in the canonical ensemble (constant volume), the heat capacity <math> C_V </math> is given by:

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

where <math>\left\langle E\right\rangle</math> is the average total energy or expected value of the internal energy of the system.

: <math> \left\langle E\right\rangle = - \frac{\partial ln Q}{\partial \beta} = -\frac{1}{Q}\frac{\partial Q}{\partial \beta} </math>
where <math> Q</math> is the partition function for the system of N interacting spins.

By evaluating the second form of the partial derivative with respect to <math> \beta</math>:

<math> \left\langle E\right\rangle = -\frac{1}{Q}\frac{\partial Q}{\partial \beta} = \frac{\Sigma \epsilon_{i} e^{-\beta \epsilon_{i}}}{\Sigma e^{-\beta \epsilon_{i}}} </math>

Applying the product rule enables the differentiation of this expression for <math>\left\langle E\right\rangle</math> with respect to temperature, giving:

<math>\frac{\partial \left\langle E\right\rangle}{\partial T} = \frac{\Sigma \epsilon_{i} \frac{\epsilon_i}{k_BT^2}e^{-\beta \epsilon_{i}}}{\Sigma e^{-\beta \epsilon_{i}}} - \frac{\Sigma \epsilon_i e^{-\beta \epsilon_i}\cdot \frac{\epsilon_i}{k_BT^2}e^{-\beta \epsilon_i}}{(\Sigma e^{-\beta \epsilon_{i}})^2}</math>

<math>\frac{\partial \left\langle E\right\rangle}{\partial T} = \frac{1}{k_BT^2}\left[\frac{\Sigma \epsilon_i^2 e^{-\beta \epsilon_i}}{\Sigma e^{-\beta \epsilon_i}} - \left(\frac{\Sigma \epsilon_{i} e^{-\beta \epsilon_{i}}}{\Sigma e^{-\beta \epsilon_{i}}}\right)^2\right] </math>

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

Thus,

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

===Heat Capacity versus Temperature===

====2x2====
[[File:HC2.png]]

====4x4====
[[File:H4.png]]

====8x8====
[[File:HC8.png]]

====16x16====
[[File:HC16.png]]

====32x32====
[[File:HC32.png]]

=== Python Scripts ===