File:Dr8sb1016.png

2018-11-21T11:16:14Z

Sb1016:

Rep:MonteCarlosim sb1016

2018-11-21T11:14:32Z

File:Dr3sb1016.png

2018-11-21T11:07:05Z

Sb1016:

File:Dr2sb1016.png

2018-11-21T11:01:18Z

Sb1016:

Rep:MonteCarlosim sb1016

2018-11-21T10:59:26Z

Sb1016:

== Introduction ==

== Introduction to 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 Ising model, the interaction energy is given by <math>- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>.
In one dimension, each atom/cell has two neighbours which repeats due to the periodicity of the lattice.
{| border="1" style="border-collapse:collapse"
|3
| 1
| 2
| 3
|}

From the table above, cell 1 interacts with cell 3 and 2, whereas 2 interacts with 1 and 3. For the lowest energy configuration, the neighbouring cells will have the same spin. Therefore, the interaction energy is given by <math>- \frac{1}{2} J N( 2 s_i s_j ) </math>, which is equal to <math>E\ =\ -NJ</math> where D=1. In 2D there are 4 interactions per cell and the lowest energy is given by <math>E\ =\ -2NJ</math> where D=2.

In the lowest energy configuration, as all the cells have the same spin, it can be either +1 or -2, giving two possible states. The multiplicity is 2. Entropy is therefore <math>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?'''

In the lowest energy configuration, the energy of the system is -3000J. When a spin is flipped, firstly the system loses 6 favourable interactions, gaining 6J. Moreover, it then has 6 unfavourable interactions which also adds 6J. In total the energy of the system increases by 12J.

In a system where N=1000, any of the 1000 spins can flip giving 1000 configurations. Therefore entropy is <math>K_b ln(1000)</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: +1

2D Lattice: +1

For the 3D Ising Lattice with 1000 cells: At 0K, lowest possible energy configuration with total magnetisation. Therefore the lattice can either have all spin up or spin down. This will give a total magnetisation of ±1000.

== Calculating energy and magnetisation ==

[[File:ILchecksb1016.PNG|thumb|500px|Fig 1.The result from ILcheck.py]]

'''<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 — we will address the speed in a later part of the experiment.'''

'''<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.'''

The result from ILcheck.py in Fig 1. shows a 4x4 lattices from a low energy state with all parallel spins to a high energy state with anti-parallel spins.

==Introduction to 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 with 100 spins, there are <math> 2^100</math> configurations. To evaluate a single <math>\left\langle M\right\rangle_T</math> value, it would take <math> 1.467\times10^16</math> days.

'''<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.'''

Results from montecarlostep() for a 8x8 lattice at T=1.0: (8.0, 2.0)

Results from the statistics() function for a 8x8 lattice at T=1.0 : (8.0, 64.0, 2.0, 4.0, 1)

[[File:ILanimsb1016.png|thumb|400px|Fig2. The result from ILanim.py]]

'''<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.'''

<math>F = U - TS</math>, where <math>F</math> is the Helmholtz free energy, <math>U</math> is the internal energy of the system, <math>T</math> is temperature and <math>S</math> is entropy.

At a given temperature, the system is controls the balance between the entropically favoured configuration and the energetically favoured configuration. When <math>T < T_C</math>, the entropic driving force is not big enough to cause spins to flip to anti-parallel spin configuration. Therefore, a large number of parallel spins configuration and spontaneous magnetisation can be observed below the Curie temperature.

From Fig2. it can be seen that after around 600 steps, the system seems to reach equilibrium.

==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!'''

{| class="wikitable"
|Time Taken(s)
|9.646776585
|9.098013143
|9.025683231
|9.186795471
|9.219164678
|-
|Average time taken(s)
|colspan="5"|9.235286621
|-
|Standard Deviation
|colspan="5"|0.242235881
|}

'''<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."
energy = -(np.sum(np.multiply(self.lattice(np.roll(self.lattice(1,0))))))-
(np.sum(np.multiply(self.lattice(np.roll(self.lattice(1,1))))))
return energy

def magnetisation(self):
"Return the total magnetisation of the current lattice configuration."
magnetisation=np.sum(self.lattice)
return magnetisation
</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!'''

With the new code, the average time for Monte Carlo steps to run has decreased by 8.73s.

{| class="wikitable"
|Time Taken(s)
|0.577154774
|0.500357554
|0.493485749
|0.486324092
|0.485946439
|-
|Average time taken(s)
|colspan="5"|0.508653721
|-
|Standard Deviation
|colspan="5"|0.038747648
|}

==The effect of temperature==

[[File:20200.0004sb1016.png|thumb|250px|Fig3. 20x20 Lattice at T=0.0004]] [[File:50501sb1016.png|thumb|250px|Fig4. 50x50 Lattice at T=1]]

'''<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.'''

As lattice size and temperature of the system increases, the Monte Carlo steps needed for equilibrium increases. For a 2x2 lattice, at T=5.0, it reached equilibrium before 20,000 steps. However, a 32x32 lattice requires nearly 50,000 steps.

It can be seen from Fig 3. and Fig 4. that 20x20 lattice at a low temperature requires far lesser steps than 50x50 at a higher temperature. In order to not waste computational energy, the minimum number of steps was set to 20,000 without a upper limit. The montecarlostep function was modified so that the monetcarlostep result is only recorded of if the standard deviation is less than 5% of the mean. The code is further explained below. No change was made to the statistics function.

<pre>
def montecarlostep(self, T):
energy_list=[]
self.n_cycles = self.n_cycles + 1
self.steps=self.steps+1
energy = self.energy()

#First montecarlostep is run and the result is added to the energy_list
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
self.lattice[random_i,random_j] = -self.lattice[random_i,random_j]
energy_diff = self.energy() - energy
random_number = np.random.random()
if energy_diff > 0:
if random_number > np.exp(-energy_diff/T):
self.lattice[random_i,random_j] = -self.lattice[random_i,random_j]
energy_list=energy_list+[energy]

#More values are added to the energy list by running the monetcarlostep in order to calculate standard deviation
while (self.steps<20000):
energy = self.energy()
#
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
self.lattice[random_i,random_j] = -self.lattice[random_i,random_j]
energy_diff = self.energy() - energy
random_number = np.random.random()
#
if energy_diff > 0:
if random_number > np.exp(-energy_diff/T):
self.lattice[random_i,random_j] = -self.lattice[random_i,random_j]
energy_list=energy_list+[energy]
self.steps=self.steps+1

#To make sure the following step only runs the first time the IF condition is specified
i=0
if self.steps==20000:
#Montecarlostep is run till the standard deviation of the energy list is less than 5% of the mean.
while i==0:
if np.std(energy_list)>abs(0.05*np.mean(energy_list)):
#
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
self.lattice[random_i,random_j] = -self.lattice[random_i,random_j]
energy_diff = self.energy() - energy
random_number = np.random.random()
#this will choose a random number in the range[0,1)
if energy_diff > 0:
if random_number > np.exp(-energy_diff/T):
self.lattice[random_i,random_j] = -self.lattice[random_i,random_j]
energy_list=energy_list+[energy]
self.steps=self.steps+1
if np.std(energy_list)<abs(0.05*np.mean(energy_list)):
i=1
#The result from the last monte carlo step is recorded
self.E = self.E + self.energy()
self.M = self.M + self.magnetisation()
self.E2 = self.E2 + self.energy()**2
self.M2 = self.M2 + self.magnetisation()**2
return(self.energy(), self.magnetisation())
</pre>

'''<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 intuition 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. The NumPy function savetxt() stores your array of output data on disk — you will need it later. Save the file as ''8x8.dat'' so that you know which lattice size it came from.'''

[[File:88tempsb1016.png|250px|alt=Alt text|Fig5. 8x8 Lattice results with error bars]]
[[File:Errorsb1016.png|250px|alt=Alt text|Fig6. Error bars after zooming in on Fig 5.]]

Temperature spacing was set to 0.01 to effectively capture long range fluctuations. After the montecarlostep is run for the first time for the lattice at a specific temperature, it is already reached the equilibrium value. The runtime was therefore set 100,000 cycles to make sure the average calculated has the smallest standard deviation and standard error possible. Due to this, the error bars being plotted (using standard error) are extremely small and cannot be seen without zooming in.

As temperature increases, energy/spin remains relatively constant at -2.0 till the critical region. Near <math>T_c</math> (between T=2.0 and T=3.0), energy/spin increases to 0.0 and remains relatively constant after.

The initial magnetisation is +1 with all parallel spins. As temperature increases close to the <math>T_c</math>, the entropic driving force takes over and cause spin flipping to maximise entropy. This creates an anti-parallel spin configuration with equal and alternating spin ups and downs, giving 0 magnetisation.

==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?'''

As the lattice size increases, long range fluctuations decrease. 16x16 is large enough to capture these fluctuations. In the critical region around <math>T_c</math>, a single change in spin can have a large impact on magnetisation. In a smaller lattice, the impact of this single spin is greater due to the small number of configurations available for spin flipping. However in a larger lattice, this impact is smaller and therefore fluctuations are also much smaller.

{| class="wikitable"
|Lattice size
|Energy/spin vs Temperature
|Magnetisation/spin vs Temperature
|-
|2x2
|[[File:22energysb1016.png|250px|Fig6.]]
|[[File:22magsb1016.png|250px|Fig11.]]
|-
|4x4
|[[File:44energysb1016.png|250px|Fig7.]]
|[[File:44magsb1016.png|250px|Fig12.]]
|-
|8x8
|[[File:88energysb1016.png|250px|Fig8.]]
|[[File:88magsb1016.png|250px|Fig13.]]
|-
|16x16
|[[File:1616energysb1016.png|250px|Fig9.]]
|[[File:1616magsb1016.png|250px|Fig14.]]
|-
|32x32
|[[File:3232energysb1016.png|250px|Fig10.]]
|[[File:3232magsb1016.png|250px|Fig15.]]
|}

==Determining 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>.)

Expectation of energy can be calculated by:
[[File:Dr1sb1016.png|250px|Fig16.]]

'''<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 — 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"
|Lattice size
|Heat capacity/spin vs Temperature
|-
|2x2
|[[File:22hcsb1016.png|250px|Fig16.]]
|-
|4x4
|[[File:44hcsb1016.png|250px|Fig17.]]
|-
|8x8
|[[File:88hcsb1016.png|250px|Fig18.]]
|-
|16x16
|[[File:1616hcsb1016.png|250px|Fig19.]]
|-
|32x32
|[[File:3232hcsb1016.png|250px|Fig20.]]
|}

==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 — 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
|32x32
|-
|Energy/spin vs Temperature
|[[File:3232Ecomparesb1016.png|250px|Fig21.]]
|-
|Magnetisation/spin vs Temperature
|[[File:3232Mcomparesb1016.png|250px|Fig22.]]
|-
|Heat Capacity/spin vs Temperature
|[[File:3232HCcomparesb1016.png|250px|Fig23.]]
|}

'''<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 — in general, it might be difficult to get a good fit! Attach a PNG of an example fit to your report.'''

[[File:441fittingsb1016.png|250px|Fig23. Initial polynomial fitting for 4x4 lattice.]]

'''<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.'''

[[File:442fittingsb1016.png|250px|Fig23.Polynomial fitting of the peak for 4x4 lattice.]]

'''<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.'''

[[File:curietempsb1016.png|250px|Fig23.Polynomial fitting of the peak for 4x4 lattice.]]

==Conclusion==

2018-11-21T08:21:19Z

Sb1016:

Rep:MonteCarlosim sb1016

2018-11-21T02:06:05Z

Sb1016:

Introduction

'''<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 Ising model, the interaction energy is given by <math>- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>.
In one dimension, each atom/cell has two neighbours which repeats due to the periodicity of the lattice.
{| border="1" style="border-collapse:collapse"
|3
| 1
| 2
| 3
|}

From the table above, cell 1 interacts with cell 3 and 2, whereas 2 interacts with 1 and 3. For the lowest energy configuration, the neighbouring cells will have the same spin. Therefore, the interaction energy is given by <math>- \frac{1}{2} J N( 2 s_i s_j ) </math>, which is equal to <math>E\ =\ -NJ</math> where D=1.
In 2D there are 4 interactions per cell and the lowest energy is given by <math>E\ =\ -2NJ</math> where D=2.

In the lowest energy configuration, as all the cells have the same spin, it can be either +1 or -2, giving two possible states. The multiplicity is 2. Entropy is therefore <math>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?'''

'''<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: +1
2D Lattice: +1
For the 3D Ising Lattice with 1000 cells: At 0K, lowest possible energy configuration with total magnetisation. Therefore the lattice can either have all spin up or spin down. This will give a total magnetisation of ±1000.

'''Calculating energy and magnetisation'''

TASK: 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 J=1.0 at all times (in fact, we are working in reduced units in which J=k_B, but there will be more information about this in later sections). Do not worry about the efficiency of the code at the moment — we will address the speed in a later part of the experiment.

TASK: Run the ILcheck.py script from the IPython Qt console using the command

[[File:ILchecksb1016.PNG|250px|thumb|Fig 1. The result from ILcheck.py]]

'''Introduction to 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 with 100 spins, there are <math> 2^100</math> configurations. To evaluate a single <math>\left\langle M\right\rangle_T</math> value, it would take <math> 1.467\times10^16</math> days.

TASK: 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 k_B! Complete the statistics() function. This should return the following quantities whenever it is called: <E>, <E^2>, <M>, <M^2>, and the number of Monte Carlo steps that have elapsed.

TASK: If T < T_C, do you expect a spontaneous magnetisation (i.e. do you expect \left\langle M\right\rangle \neq 0)? 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.

From Fig2. it can be seen that after around 600 steps, the system seems to reach equilibrium.

[[File:ILanimsb1016.png|500px|thumb|Fig2. The result from ILanim.py]]

'''Accelerating the code'''

TASK: 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!

{| class="wikitable"
|Time Taken(s)
|9.646776585
|9.098013143
|9.025683231
|9.186795471
|9.219164678
|-
|Average time taken(s)
|colspan="5"|9.235286621
|-
|Standard Deviation
|colspan="5"|0.242235881
|}

TASK: Look at the documentation for the 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 roll and 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."
energy = -(np.sum(np.multiply(self.lattice(np.roll(self.lattice(1,0))))))-
(np.sum(np.multiply(self.lattice(np.roll(self.lattice(1,1))))))
return energy

def magnetisation(self):
"Return the total magnetisation of the current lattice configuration."
magnetisation=np.sum(self.lattice)
return magnetisation
</pre>

TASK: 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!

{| class="wikitable"
|Time Taken(s)
|0.577154774
|0.500357554
|0.493485749
|0.486324092
|0.485946439
|-
|Average time taken(s)
|colspan="5"|0.508653721
|-
|Standard Deviation
|colspan="5"|0.038747648
|}

'''The effect of temperature'''

TASK: 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.

As lattice size and temperature of the system increases, the Monte Carlo steps needed for equilibrium increases. For a 2x2 lattice, at T=5.0, it reached equilibrium before 20,000 steps. However, a 32x32 lattice requires nearly 50,000 steps.

[[File:20200.0004sb1016.png|500px|thumb|Fig3. 20x20 Lattice at T=0.0004]]
[[File:50501sb1016.png|500px|thumb|Fig4. 50x50 Lattice at T=1]]

It can be seen from Fig 3. and Fig 4. that 20x20 lattice at a low temperature requires far lesser steps than 50x50 at a higher temperature. In order to not waste computational energy, the minimum number of steps was set to 20,000 without a upper limit. The montecarlostep function was modified so that the monetcarlostep result is only recorded of if the standard deviation is less than 5% of the mean. The code is further explained below. No change was made to the statistics function.

<pre>
def montecarlostep(self, T):
energy_list=[]
self.n_cycles = self.n_cycles + 1
self.steps=self.steps+1
energy = self.energy()

#First montecarlostep is run and the result is added to the energy_list
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
self.lattice[random_i,random_j] = -self.lattice[random_i,random_j]
energy_diff = self.energy() - energy
random_number = np.random.random()
if energy_diff > 0:
if random_number > np.exp(-energy_diff/T):
self.lattice[random_i,random_j] = -self.lattice[random_i,random_j]
energy_list=energy_list+[energy]

#More values are added to the energy list by running the monetcarlostep in order to calculate standard deviation
while (self.steps<20000):
energy = self.energy()
#
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
self.lattice[random_i,random_j] = -self.lattice[random_i,random_j]
energy_diff = self.energy() - energy
random_number = np.random.random()
#
if energy_diff > 0:
if random_number > np.exp(-energy_diff/T):
self.lattice[random_i,random_j] = -self.lattice[random_i,random_j]
energy_list=energy_list+[energy]
self.steps=self.steps+1

#To make sure the following step only runs the first time the IF condition is specified
i=0
if self.steps==20000:
#Montecarlostep is run till the standard deviation of the energy list is less than 5% of the mean.
while i==0:
if np.std(energy_list)>abs(0.05*np.mean(energy_list)):
#
random_i = np.random.choice(range(0, self.n_rows))
random_j = np.random.choice(range(0, self.n_cols))
self.lattice[random_i,random_j] = -self.lattice[random_i,random_j]
energy_diff = self.energy() - energy
random_number = np.random.random()
#this will choose a random number in the range[0,1)
if energy_diff > 0:
if random_number > np.exp(-energy_diff/T):
self.lattice[random_i,random_j] = -self.lattice[random_i,random_j]
energy_list=energy_list+[energy]
self.steps=self.steps+1
if np.std(energy_list)<abs(0.05*np.mean(energy_list)):
i=1
#The result from the last monte carlo step is recorded
self.E = self.E + self.energy()
self.M = self.M + self.magnetisation()
self.E2 = self.E2 + self.energy()**2
self.M2 = self.M2 + self.magnetisation()**2
return(self.energy(), self.magnetisation())
</pre>

TASK: Use ILtemperaturerange.py to plot the average energy and magnetisation for each temperature, with error bars, for an 8\times 8 lattice. Use your intuition 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. The NumPy function savetxt() stores your array of output data on disk — you will need it later. Save the file as 8x8.dat so that you know which lattice size it came from.

Temperature spacing was set to 0.01 to effectively capture long range fluctuations. After the montecarlostep is run for the first time for the lattice at a specific temperature, it is already reached the equilibrium value. The runtime was therefore set 100,000 cycles to make sure the average calculated has the smallest standard deviation and standard error possible.

[[File:88tempsb1016.png|250px|Fig5.]]

'''The effect of system size'''

TASK: 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
|2x2
|4x4
|8x8
|16x16
|32x32
|-
|Energy/spin vs Temperature
|[[File:22energysb1016.png|250px|Fig6.]]
|[[File:44energysb1016.png|250px|Fig7.]]
|[[File:88energysb1016.png|250px|Fig8.]]
|[[File:1616energysb1016.png|250px|Fig9.]]
|[[File:3232energysb1016.png|250px|Fig10.]]
|-
|Magnetisation/spin vs Temperature
|[[File:22magsb1016.png|250px|Fig11.]]
|[[File:44magsb1016.png|250px|Fig12.]]
|[[File:88magsb1016.png|250px|Fig13.]]
|[[File:1616magsb1016.png|250px|Fig14.]]
|[[File:3232magsb1016.png|250px|Fig15.]]
|}

'''determining heat capacity'''

TASK: 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, \mathrm{Var}[X], the mean of its square \left\langle X^2\right\rangle, and its squared mean \left\langle X\right\rangle^2. You may find that the data around the peak is very noisy — 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"
|Lattice size
|2x2
|4x4
|8x8
|16x16
|32x32
|-
|Heat capacity/spin vs Temperature
|[[File:22hcsb1016.png|250px|Fig16.]]
|[[File:44hcsb1016.png|250px|Fig17.]]
|[[File:88hcsb1016.png|250px|Fig18.]]
|[[File:1616hcsb1016.png|250px|Fig19.]]
|[[File:3232hcsb1016.png|250px|Fig20.]]
|}

'''Locating the Curie temperature'''

TASK: 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 here if you are interested. Each file contains six columns: T, E, E^2, M, M^2, C (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 — 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 here).

{| class="wikitable"
|Lattice size
|32x32
|-
|Energy/spin vs Temperature
|[[File:3232Ecomparesb1016.png|250px|Fig21.]]
|-
|Magnetisation/spin vs Temperature
|[[File:3232Mcomparesb1016.png|250px|Fig22.]]
|-
|Heat Capacity/spin vs Temperature
|[[File:3232HCcomparesb1016.png|250px|Fig23.]]
|}

TASK: 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 — in general, it might be difficult to get a good fit! Attach a PNG of an example fit to your report.

[[File:441fittingsb1016.png|250px|Fig23. Initial polynomial fitting for 4x4 lattice.]]

TASK: 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.

[[File:442fittingsb1016.png|250px|Fig23.Polynomial fitting of the peak for 4x4 lattice.]]

TASK: 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 T_C for that side length. Make a plot that uses the scaling relation given above to determine T_{C,\infty}. 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.

[[File:curietempsb1016.png|250px|Fig23.Polynomial fitting of the peak for 4x4 lattice.]]

Conclusion