Ly2412: /* Locating the Curie temperature */

==Introduction to the Ising model==
'''<big>TASK1</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.'''

For a spin in a system of dimension <math>D</math>, there are <math>2D</math> interactions associated with it due to its <math>2D</math> number of neighbouring spins. To minimise the energy of these interactions, this spin and its <math>2D</math> neighbours must have the same sign, i.e. all up or all down. Therefore the minimum energy of these interactions associated with this one spin is <math>E\ =-0.5\times2DJ = -DJ</math>. If the system has total number of spins of <math>N</math>, then the overall lowest energy must be <math>E\ =\ -DNJ</math>. The multiplicity of such a system should be 2 since the whole system can either have all of its spins up or down. The entropy of such a system is therefore <math>S\ =\ K_b\ln 2=9.57\times 10^{-24} JK^{-1} </math>.

'''<big>TASK2</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?'''

For a system of <math>D=3,\ N=1000</math> in its lowest energy configuration, the energy of this configuration is <math>E\ =-3000J</math> and the entropy of this system is <math>S\ =9.57\times 10^{-24} JK^{-1} </math>. If one of the spins change its direction, there are 6 unique interaction energies changes from -1 to +1 and the overall energy changes to <math>E\ =-2988J</math>. In this excited state, the entropy is <math>S\ =\ K_b\ln 2000=1.05\times 10^{-22} JK^{-1} </math> because there are 2000 ways of arranging this system. Therefore the change in energy of this process is <math>\Delta E\ =12J</math> and the change in entropy is <math>\Delta S\ =9.54\times 10^{-23} JK^{-1}</math>.

'''<big>TASK3</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?'''

The magnetisations of the 1D and 2D lattices in figure 1 are both +1. For a Ising lattice with <math>D = 3,\ N=1000</math> at absolute zero, its magnetisation can be either +1000 or -1000.

==Calculating the energy and magnetisation==
In this section, two functions from the IsingLattice.py file, energy() and magnetisation(), were completed to calculate the energy and magnetisation of any 2D Ising lattice and a screenshot of these codes is shown below in '''Figure 1'''.

[[File:LY IL code.png|thumb|center|500px|'''Figure 1. The original code for energy() and magnetisation() functions''']]

The validity of these two functions were checked using a given ILcheck.py file and the output is given below in '''Figure 2'''.

[[File:LY_ILcheck.png|thumb|center|400px|'''Figure 2. The output from the ILcheck.py file''']]

==Introduction to Monte Carlo simulation==

'''<big>TASK1</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, since each spin can be either up or down, the total number of configurations is therefore <math>W\ =2^{100}=1.27\times 10^{30}</math> and it will take about <math>4.02\times 10^{13}</math> years to calculate a single value of <math>\left\langle M\right\rangle_T</math>.

'''<big>TASK2</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.'''

A montecarlostep(T) function and a statistics() function were written in this section to perform a single step Monte Carlo simulation and calculate average properties of the system, respectively. The output from ILanim.py file is shown below in '''Figure 3'''.
{{multiple image
| image1 = LY_ILanim_out.png

| image2 = LY_ILanim_out_2.png

| direction = vertical
| align = center
| width = 350
| footer = '''Figure 3. The output from ILanim.py and statistics()'''
}}

==Accelerating the code==
'''<big>TASK1</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!'''

'''<big>TASK2</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!'''

In this section, the energy() and statistics() functions were modified using numpy.sum(), numpy.roll() and numpy.multiply() functions. After the modification, the calculation time has been shortened significantly in performing 2000 Monte Carlo steps. The comparison is shown below in '''Table 1'''.

{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 1. The effect of accelerating codes on computational time
! !! before acceleration !! after acceleration
|-
| screen shot || [[File:LY IL MC 2000 original.png|400px]] || [[File:LY IL MC 2000 accelerated final.png|400px]]
|-
| average (s)|| 21.6 || 0.283
|-
| standard error (s)|| 0.162 || 0.00140
|}

==The effect of temperature==
'''<big>TASK1</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.'''

From the top picture in '''Figure 3''', we can see that the first few hundreds step of the Monte Carlo calculation are relaxation steps. Therefore, these steps should be excluded during the calculation of average properties of the system. Some outputs from the ILfinalframe.py file for the 8x8 lattice at different temperatures are shown below in '''Table 2'''. From this table, it is clear that the first 2000 steps should be excluded and a total run-time of 50000 should be appropriate in the calculation of average proprieties of this lattice size.

{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 2. The effect of temperature on the total number of relaxation steps of the 8x8 lattice
! Temperature(K) !! 0.5 !! 1.5 !! 2.0 !! 2.5
|-
| Output Screenshot || [[File:LY IL MC 2000 original.png|300px]] || [[File:LY IL MC 2000 accelerated final.png|300px]] || ||
|}

'''<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 — you will need it later. Save the file as ''8x8.dat'' so that you know which lattice size it came from.'''

From '''Table 2''', it is clear that the first 2000 steps should be excluded and a total run-time of 50000 should be appropriate in the calculation of average proprieties of the 8x8 lattice. The output from the IL.temperaturerange.py file of a 8x8 lattice is shown below in '''Figure 4'''. The temperature range is from 0.2 to 5.0 and step-size is 0.1. Errors bars of average energies and magnetisations at each temperature point are shown as well. In fact these error bars are too small to be seen clearly. However, this does not mean that the calculations are accurate enough......
[[File:LY_88.png|thumb|center|500px|'''Figure 4. The output from the ILtemperature.py file for 8x8 lattice''']]

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

The number of relaxation and equilibration steps needed to calculate reliable average properties depend strongly on the system size. Therefore, several ILfinalframe.py calculations were performed for different lattice sizes to determine the number of relaxation steps and total run-time required. The results are summarised in '''Table 3''' below and these values are used when ILtemperaturerange.py calculations were performed for lattice sizes: 2x2, 4x4, 8x8, 16x16, 32x32.
{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 3. The relaxation and equilibration steps for different lattice sizes
! lattice size !! 2x2 !! 4x4!!8x8!!16x16!!32x32
|-
| number of relaxation steps || 500 || 1000||2000 ||20000 || 75000
|-
| total run time||20000 ||30000 || 50000||150000 ||300000
|}

The ILtemperaturerange.py script was then used to calculate average energies and magnetisations of the following lattice sizes: 2x2, 4x4, 16x16, 32x32 and results are shown in '''Table 4'''.
{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 4. The output of ILtemperaturerange.py for different lattice sizes
! lattice size !! 2x2 !! 4x4
|-
| output || [[File:LY_22.png|400px]] || [[File:LY_44.png|400px]]
|-
| <center>'''lattice size'''</center> ||<center> '''16x16'''</center> ||<center> '''32x32'''</center>
|-
| output || [[File:LY_1616.png|400px]] || [[File:LY_3232.png|400px]]
|}

''Note:'' The temperature ranges in all the above calculations are all from 0.2 to 5.0 with the same step size 0.1.

The energies and magnetisations per spin versus temperature graphs for each lattice size were then plotted together for comparisons, as shown in '''Figure 5'''.
{{multiple image
| image1 = LY_energyallin1.png

| image2 = LY_magallin1.png

| direction = vertical
| align = center
| width = 450
| footer = '''Figure 5. The energies and magnetisations vs. temperature graphs for different lattice sizes shown in one figure'''
}}

From '''Figure 5''', it is clear to see that by increasing the lattice size from 2x2 to 4x4 and from 4x4 to 8x8, both energies and magnetisations per temperature curves change in shapes significantly. In other words, by increasing the lattice size from 2x2 to 8x8, the quality of the calculations has been greatly improved and more details of the behaviors of both energies and magnetisations per temperature curves have been revealed. However, by changing the lattice size from 8x8 to 16x16 and 32x32, both plots in '''Figure 5''' do not change in general features, where as the computational time increases significantly. Therefore, the 8x8 lattice can be seen as the minimum lattice size required to capture the long range fluctuations.

==Determining the heat capacity==
'''<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.'''

Since the heat capacity is defined as <math>C = \frac{\partial E}{\partial T} = \frac{\mathrm{Var}[E]}{k_B T^2}</math>, by extracting the <math>\left\langle E^2\right\rangle</math> and <math>\left\langle E\right\rangle</math> data from the ''.dat'' files generated previously, heat capacity versus temperature curves for each lattice size were computed and plotted in '''Figure 6'''.

[[File:LY_heatcapacityallin1.png|thumb|center|500px|'''Figure 6. Heat capacity versus temperature curves for all the lattice sizes''']]

In '''Figure 6''', as lattice size increases, the temperature as which maximum heat capacity occurs shifts to the left gradually. This is expected due to the finite size effect, which is very common in the simulations of systems bearing long ranged interactions. The trend in '''Figure 6''' is actually in good agreement with the relationship <math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math>, where <math>T_{C, L}</math> is the Curie temperature of an <math>L\times L</math>lattice, <math>T_{C,\infty}</math> is the Curie temperature of an infinite lattice, and <math>A</math> is a constant.

==Locating the Curie temperature==
'''<big>TASK1</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]).'''

The heat capacity versus temperature curves calculated using data generated by Phthon and C++ for the 32x32 lattice size were plotted in '''Figure 7'''.

[[File:LY_3232compare.png|thumb|center|500px|'''Figure 7. Comparisons between heat capacity calculated from C++ and Python data for the lattice size of 32x32''']]

'''<big>TASK2</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.'''

'''<big>TASK3</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.'''

'''<big>TASK4</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:LY 3232compare.png

2015-02-11T15:20:47Z

Ly2412:

Rep:Mod:PYTHONLY2412

2015-02-11T15:18:13Z

Ly2412: /* The effect of system size */

==Introduction to the Ising model==
'''<big>TASK1</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.'''

For a spin in a system of dimension <math>D</math>, there are <math>2D</math> interactions associated with it due to its <math>2D</math> number of neighbouring spins. To minimise the energy of these interactions, this spin and its <math>2D</math> neighbours must have the same sign, i.e. all up or all down. Therefore the minimum energy of these interactions associated with this one spin is <math>E\ =-0.5\times2DJ = -DJ</math>. If the system has total number of spins of <math>N</math>, then the overall lowest energy must be <math>E\ =\ -DNJ</math>. The multiplicity of such a system should be 2 since the whole system can either have all of its spins up or down. The entropy of such a system is therefore <math>S\ =\ K_b\ln 2=9.57\times 10^{-24} JK^{-1} </math>.

'''<big>TASK2</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?'''

For a system of <math>D=3,\ N=1000</math> in its lowest energy configuration, the energy of this configuration is <math>E\ =-3000J</math> and the entropy of this system is <math>S\ =9.57\times 10^{-24} JK^{-1} </math>. If one of the spins change its direction, there are 6 unique interaction energies changes from -1 to +1 and the overall energy changes to <math>E\ =-2988J</math>. In this excited state, the entropy is <math>S\ =\ K_b\ln 2000=1.05\times 10^{-22} JK^{-1} </math> because there are 2000 ways of arranging this system. Therefore the change in energy of this process is <math>\Delta E\ =12J</math> and the change in entropy is <math>\Delta S\ =9.54\times 10^{-23} JK^{-1}</math>.

'''<big>TASK3</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?'''

The magnetisations of the 1D and 2D lattices in figure 1 are both +1. For a Ising lattice with <math>D = 3,\ N=1000</math> at absolute zero, its magnetisation can be either +1000 or -1000.

==Calculating the energy and magnetisation==
In this section, two functions from the IsingLattice.py file, energy() and magnetisation(), were completed to calculate the energy and magnetisation of any 2D Ising lattice and a screenshot of these codes is shown below in '''Figure 1'''.

[[File:LY IL code.png|thumb|center|500px|'''Figure 1. The original code for energy() and magnetisation() functions''']]

The validity of these two functions were checked using a given ILcheck.py file and the output is given below in '''Figure 2'''.

[[File:LY_ILcheck.png|thumb|center|400px|'''Figure 2. The output from the ILcheck.py file''']]

==Introduction to Monte Carlo simulation==

'''<big>TASK1</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, since each spin can be either up or down, the total number of configurations is therefore <math>W\ =2^{100}=1.27\times 10^{30}</math> and it will take about <math>4.02\times 10^{13}</math> years to calculate a single value of <math>\left\langle M\right\rangle_T</math>.

'''<big>TASK2</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.'''

A montecarlostep(T) function and a statistics() function were written in this section to perform a single step Monte Carlo simulation and calculate average properties of the system, respectively. The output from ILanim.py file is shown below in '''Figure 3'''.
{{multiple image
| image1 = LY_ILanim_out.png

| image2 = LY_ILanim_out_2.png

| direction = vertical
| align = center
| width = 350
| footer = '''Figure 3. The output from ILanim.py and statistics()'''
}}

==Accelerating the code==
'''<big>TASK1</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!'''

'''<big>TASK2</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!'''

In this section, the energy() and statistics() functions were modified using numpy.sum(), numpy.roll() and numpy.multiply() functions. After the modification, the calculation time has been shortened significantly in performing 2000 Monte Carlo steps. The comparison is shown below in '''Table 1'''.

{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 1. The effect of accelerating codes on computational time
! !! before acceleration !! after acceleration
|-
| screen shot || [[File:LY IL MC 2000 original.png|400px]] || [[File:LY IL MC 2000 accelerated final.png|400px]]
|-
| average (s)|| 21.6 || 0.283
|-
| standard error (s)|| 0.162 || 0.00140
|}

==The effect of temperature==
'''<big>TASK1</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.'''

From the top picture in '''Figure 3''', we can see that the first few hundreds step of the Monte Carlo calculation are relaxation steps. Therefore, these steps should be excluded during the calculation of average properties of the system. Some outputs from the ILfinalframe.py file for the 8x8 lattice at different temperatures are shown below in '''Table 2'''. From this table, it is clear that the first 2000 steps should be excluded and a total run-time of 50000 should be appropriate in the calculation of average proprieties of this lattice size.

{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 2. The effect of temperature on the total number of relaxation steps of the 8x8 lattice
! Temperature(K) !! 0.5 !! 1.5 !! 2.0 !! 2.5
|-
| Output Screenshot || [[File:LY IL MC 2000 original.png|300px]] || [[File:LY IL MC 2000 accelerated final.png|300px]] || ||
|}

'''<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 — you will need it later. Save the file as ''8x8.dat'' so that you know which lattice size it came from.'''

From '''Table 2''', it is clear that the first 2000 steps should be excluded and a total run-time of 50000 should be appropriate in the calculation of average proprieties of the 8x8 lattice. The output from the IL.temperaturerange.py file of a 8x8 lattice is shown below in '''Figure 4'''. The temperature range is from 0.2 to 5.0 and step-size is 0.1. Errors bars of average energies and magnetisations at each temperature point are shown as well. In fact these error bars are too small to be seen clearly. However, this does not mean that the calculations are accurate enough......
[[File:LY_88.png|thumb|center|500px|'''Figure 4. The output from the ILtemperature.py file for 8x8 lattice''']]

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

The number of relaxation and equilibration steps needed to calculate reliable average properties depend strongly on the system size. Therefore, several ILfinalframe.py calculations were performed for different lattice sizes to determine the number of relaxation steps and total run-time required. The results are summarised in '''Table 3''' below and these values are used when ILtemperaturerange.py calculations were performed for lattice sizes: 2x2, 4x4, 8x8, 16x16, 32x32.
{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 3. The relaxation and equilibration steps for different lattice sizes
! lattice size !! 2x2 !! 4x4!!8x8!!16x16!!32x32
|-
| number of relaxation steps || 500 || 1000||2000 ||20000 || 75000
|-
| total run time||20000 ||30000 || 50000||150000 ||300000
|}

The ILtemperaturerange.py script was then used to calculate average energies and magnetisations of the following lattice sizes: 2x2, 4x4, 16x16, 32x32 and results are shown in '''Table 4'''.
{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 4. The output of ILtemperaturerange.py for different lattice sizes
! lattice size !! 2x2 !! 4x4
|-
| output || [[File:LY_22.png|400px]] || [[File:LY_44.png|400px]]
|-
| <center>'''lattice size'''</center> ||<center> '''16x16'''</center> ||<center> '''32x32'''</center>
|-
| output || [[File:LY_1616.png|400px]] || [[File:LY_3232.png|400px]]
|}

''Note:'' The temperature ranges in all the above calculations are all from 0.2 to 5.0 with the same step size 0.1.

The energies and magnetisations per spin versus temperature graphs for each lattice size were then plotted together for comparisons, as shown in '''Figure 5'''.
{{multiple image
| image1 = LY_energyallin1.png

| image2 = LY_magallin1.png

| direction = vertical
| align = center
| width = 450
| footer = '''Figure 5. The energies and magnetisations vs. temperature graphs for different lattice sizes shown in one figure'''
}}

From '''Figure 5''', it is clear to see that by increasing the lattice size from 2x2 to 4x4 and from 4x4 to 8x8, both energies and magnetisations per temperature curves change in shapes significantly. In other words, by increasing the lattice size from 2x2 to 8x8, the quality of the calculations has been greatly improved and more details of the behaviors of both energies and magnetisations per temperature curves have been revealed. However, by changing the lattice size from 8x8 to 16x16 and 32x32, both plots in '''Figure 5''' do not change in general features, where as the computational time increases significantly. Therefore, the 8x8 lattice can be seen as the minimum lattice size required to capture the long range fluctuations.

==Determining the heat capacity==
'''<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.'''

Since the heat capacity is defined as <math>C = \frac{\partial E}{\partial T} = \frac{\mathrm{Var}[E]}{k_B T^2}</math>, by extracting the <math>\left\langle E^2\right\rangle</math> and <math>\left\langle E\right\rangle</math> data from the ''.dat'' files generated previously, heat capacity versus temperature curves for each lattice size were computed and plotted in '''Figure 6'''.

[[File:LY_heatcapacityallin1.png|thumb|center|500px|'''Figure 6. Heat capacity versus temperature curves for all the lattice sizes''']]

In '''Figure 6''', as lattice size increases, the temperature as which maximum heat capacity occurs shifts to the left gradually. This is expected due to the finite size effect, which is very common in the simulations of systems bearing long ranged interactions. The trend in '''Figure 6''' is actually in good agreement with the relationship <math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math>, where <math>T_{C, L}</math> is the Curie temperature of an <math>L\times L</math>lattice, <math>T_{C,\infty}</math> is the Curie temperature of an infinite lattice, and <math>A</math> is a constant.

==Locating the Curie temperature==
'''<big>TASK1</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]).'''

The heat capacity versus temperature curves calculated using data generated by Phthon and C++ for the 32x32 lattice size were plotted in '''Figure 7'''.

'''<big>TASK2</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.'''

'''<big>TASK3</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.'''

'''<big>TASK4</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.'''

Rep:Mod:PYTHONLY2412

2015-02-11T15:16:17Z

Ly2412: /* Locating the Curie temperature */

==Introduction to the Ising model==
'''<big>TASK1</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.'''

For a spin in a system of dimension <math>D</math>, there are <math>2D</math> interactions associated with it due to its <math>2D</math> number of neighbouring spins. To minimise the energy of these interactions, this spin and its <math>2D</math> neighbours must have the same sign, i.e. all up or all down. Therefore the minimum energy of these interactions associated with this one spin is <math>E\ =-0.5\times2DJ = -DJ</math>. If the system has total number of spins of <math>N</math>, then the overall lowest energy must be <math>E\ =\ -DNJ</math>. The multiplicity of such a system should be 2 since the whole system can either have all of its spins up or down. The entropy of such a system is therefore <math>S\ =\ K_b\ln 2=9.57\times 10^{-24} JK^{-1} </math>.

'''<big>TASK2</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?'''

For a system of <math>D=3,\ N=1000</math> in its lowest energy configuration, the energy of this configuration is <math>E\ =-3000J</math> and the entropy of this system is <math>S\ =9.57\times 10^{-24} JK^{-1} </math>. If one of the spins change its direction, there are 6 unique interaction energies changes from -1 to +1 and the overall energy changes to <math>E\ =-2988J</math>. In this excited state, the entropy is <math>S\ =\ K_b\ln 2000=1.05\times 10^{-22} JK^{-1} </math> because there are 2000 ways of arranging this system. Therefore the change in energy of this process is <math>\Delta E\ =12J</math> and the change in entropy is <math>\Delta S\ =9.54\times 10^{-23} JK^{-1}</math>.

'''<big>TASK3</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?'''

The magnetisations of the 1D and 2D lattices in figure 1 are both +1. For a Ising lattice with <math>D = 3,\ N=1000</math> at absolute zero, its magnetisation can be either +1000 or -1000.

==Calculating the energy and magnetisation==
In this section, two functions from the IsingLattice.py file, energy() and magnetisation(), were completed to calculate the energy and magnetisation of any 2D Ising lattice and a screenshot of these codes is shown below in '''Figure 1'''.

[[File:LY IL code.png|thumb|center|500px|'''Figure 1. The original code for energy() and magnetisation() functions''']]

The validity of these two functions were checked using a given ILcheck.py file and the output is given below in '''Figure 2'''.

[[File:LY_ILcheck.png|thumb|center|400px|'''Figure 2. The output from the ILcheck.py file''']]

==Introduction to Monte Carlo simulation==

'''<big>TASK1</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, since each spin can be either up or down, the total number of configurations is therefore <math>W\ =2^{100}=1.27\times 10^{30}</math> and it will take about <math>4.02\times 10^{13}</math> years to calculate a single value of <math>\left\langle M\right\rangle_T</math>.

'''<big>TASK2</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.'''

A montecarlostep(T) function and a statistics() function were written in this section to perform a single step Monte Carlo simulation and calculate average properties of the system, respectively. The output from ILanim.py file is shown below in '''Figure 3'''.
{{multiple image
| image1 = LY_ILanim_out.png

| image2 = LY_ILanim_out_2.png

| direction = vertical
| align = center
| width = 350
| footer = '''Figure 3. The output from ILanim.py and statistics()'''
}}

==Accelerating the code==
'''<big>TASK1</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!'''

'''<big>TASK2</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!'''

In this section, the energy() and statistics() functions were modified using numpy.sum(), numpy.roll() and numpy.multiply() functions. After the modification, the calculation time has been shortened significantly in performing 2000 Monte Carlo steps. The comparison is shown below in '''Table 1'''.

{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 1. The effect of accelerating codes on computational time
! !! before acceleration !! after acceleration
|-
| screen shot || [[File:LY IL MC 2000 original.png|400px]] || [[File:LY IL MC 2000 accelerated final.png|400px]]
|-
| average (s)|| 21.6 || 0.283
|-
| standard error (s)|| 0.162 || 0.00140
|}

==The effect of temperature==
'''<big>TASK1</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.'''

From the top picture in '''Figure 3''', we can see that the first few hundreds step of the Monte Carlo calculation are relaxation steps. Therefore, these steps should be excluded during the calculation of average properties of the system. Some outputs from the ILfinalframe.py file for the 8x8 lattice at different temperatures are shown below in '''Table 2'''. From this table, it is clear that the first 2000 steps should be excluded and a total run-time of 50000 should be appropriate in the calculation of average proprieties of this lattice size.

{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 2. The effect of temperature on the total number of relaxation steps of the 8x8 lattice
! Temperature(K) !! 0.5 !! 1.5 !! 2.0 !! 2.5
|-
| Output Screenshot || [[File:LY IL MC 2000 original.png|300px]] || [[File:LY IL MC 2000 accelerated final.png|300px]] || ||
|}

'''<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 — you will need it later. Save the file as ''8x8.dat'' so that you know which lattice size it came from.'''

From '''Table 2''', it is clear that the first 2000 steps should be excluded and a total run-time of 50000 should be appropriate in the calculation of average proprieties of the 8x8 lattice. The output from the IL.temperaturerange.py file of a 8x8 lattice is shown below in '''Figure 4'''. The temperature range is from 0.2 to 5.0 and step-size is 0.1. Errors bars of average energies and magnetisations at each temperature point are shown as well. In fact these error bars are too small to be seen clearly. However, this does not mean that the calculations are accurate enough......
[[File:LY_88.png|thumb|center|500px|'''Figure 4. The output from the ILtemperature.py file for 8x8 lattice''']]

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

The number of relaxation and equilibration steps needed to calculate reliable average properties depend strongly on the system size. Therefore, several ILfinalframe.py calculations were performed for different lattice sizes to determine the number of relaxation steps and total run-time required. The results are summarised in '''Table 3''' below and these values are used when ILtemperaturerange.py calculations were performed for lattice sizes: 2x2, 4x4, 8x8, 16x16, 32x32.
{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 3. The relaxation and equilibration steps for different lattice sizes
! lattice size !! 2x2 !! 4x4!!8x8!!16x16!!32x32
|-
| number of relaxation steps || 500 || 1000||2000 ||20000 || 75000
|-
| total run time||20000 ||30000 || 50000||150000 ||300000
|}

The ILtemperaturerange.py script was then used to calculate average energies and magnetisations of the following lattice sizes: 2x2, 4x4, 16x16, 32x32 and results are shown in '''Table 4'''.
{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 4. The output of ILtemperaturerange.py for different lattice sizes
! lattice size !! 2x2 !! 4x4
|-
| output || [[File:LY_22.png|400px]] || [[File:LY_44.png|400px]]
|-
| <center>'''lattice size'''</center> ||<center> '''16x16'''</center> ||<center> '''32x32'''</center>
|-
| output || [[File:LY_1616.png|400px]] || [[File:LY_3232.png|400px]]
|}

The energies and magnetisations per spin versus temperature graphs for each lattice size were then plotted together for comparisons, as shown in '''Figure 5'''.
{{multiple image
| image1 = LY_energyallin1.png

| image2 = LY_magallin1.png

| direction = vertical
| align = center
| width = 450
| footer = '''Figure 5. The energies and magnetisations vs. temperature graphs for different lattice sizes shown in one figure'''
}}

From '''Figure 5''', it is clear to see that by increasing the lattice size from 2x2 to 4x4 and from 4x4 to 8x8, both energies and magnetisations per temperature curves change in shapes significantly. In other words, by increasing the lattice size from 2x2 to 8x8, the quality of the calculations has been greatly improved and more details of the behaviors of both energies and magnetisations per temperature curves have been revealed. However, by changing the lattice size from 8x8 to 16x16 and 32x32, both plots in '''Figure 5''' do not change in general features, where as the computational time increases significantly. Therefore, the 8x8 lattice can be seen as the minimum lattice size required to capture the long range fluctuations.

==Determining the heat capacity==
'''<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.'''

Since the heat capacity is defined as <math>C = \frac{\partial E}{\partial T} = \frac{\mathrm{Var}[E]}{k_B T^2}</math>, by extracting the <math>\left\langle E^2\right\rangle</math> and <math>\left\langle E\right\rangle</math> data from the ''.dat'' files generated previously, heat capacity versus temperature curves for each lattice size were computed and plotted in '''Figure 6'''.

[[File:LY_heatcapacityallin1.png|thumb|center|500px|'''Figure 6. Heat capacity versus temperature curves for all the lattice sizes''']]

In '''Figure 6''', as lattice size increases, the temperature as which maximum heat capacity occurs shifts to the left gradually. This is expected due to the finite size effect, which is very common in the simulations of systems bearing long ranged interactions. The trend in '''Figure 6''' is actually in good agreement with the relationship <math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math>, where <math>T_{C, L}</math> is the Curie temperature of an <math>L\times L</math>lattice, <math>T_{C,\infty}</math> is the Curie temperature of an infinite lattice, and <math>A</math> is a constant.

==Locating the Curie temperature==
'''<big>TASK1</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]).'''

The heat capacity versus temperature curves calculated using data generated by Phthon and C++ for the 32x32 lattice size were plotted in '''Figure 7'''.

'''<big>TASK2</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.'''

'''<big>TASK3</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.'''

'''<big>TASK4</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.'''

Rep:Mod:PYTHONLY2412

2015-02-11T15:09:22Z

Ly2412: /* Determining the heat capacity */

==Introduction to the Ising model==
'''<big>TASK1</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.'''

For a spin in a system of dimension <math>D</math>, there are <math>2D</math> interactions associated with it due to its <math>2D</math> number of neighbouring spins. To minimise the energy of these interactions, this spin and its <math>2D</math> neighbours must have the same sign, i.e. all up or all down. Therefore the minimum energy of these interactions associated with this one spin is <math>E\ =-0.5\times2DJ = -DJ</math>. If the system has total number of spins of <math>N</math>, then the overall lowest energy must be <math>E\ =\ -DNJ</math>. The multiplicity of such a system should be 2 since the whole system can either have all of its spins up or down. The entropy of such a system is therefore <math>S\ =\ K_b\ln 2=9.57\times 10^{-24} JK^{-1} </math>.

'''<big>TASK2</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?'''

For a system of <math>D=3,\ N=1000</math> in its lowest energy configuration, the energy of this configuration is <math>E\ =-3000J</math> and the entropy of this system is <math>S\ =9.57\times 10^{-24} JK^{-1} </math>. If one of the spins change its direction, there are 6 unique interaction energies changes from -1 to +1 and the overall energy changes to <math>E\ =-2988J</math>. In this excited state, the entropy is <math>S\ =\ K_b\ln 2000=1.05\times 10^{-22} JK^{-1} </math> because there are 2000 ways of arranging this system. Therefore the change in energy of this process is <math>\Delta E\ =12J</math> and the change in entropy is <math>\Delta S\ =9.54\times 10^{-23} JK^{-1}</math>.

'''<big>TASK3</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?'''

The magnetisations of the 1D and 2D lattices in figure 1 are both +1. For a Ising lattice with <math>D = 3,\ N=1000</math> at absolute zero, its magnetisation can be either +1000 or -1000.

==Calculating the energy and magnetisation==
In this section, two functions from the IsingLattice.py file, energy() and magnetisation(), were completed to calculate the energy and magnetisation of any 2D Ising lattice and a screenshot of these codes is shown below in '''Figure 1'''.

[[File:LY IL code.png|thumb|center|500px|'''Figure 1. The original code for energy() and magnetisation() functions''']]

The validity of these two functions were checked using a given ILcheck.py file and the output is given below in '''Figure 2'''.

[[File:LY_ILcheck.png|thumb|center|400px|'''Figure 2. The output from the ILcheck.py file''']]

==Introduction to Monte Carlo simulation==

'''<big>TASK1</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, since each spin can be either up or down, the total number of configurations is therefore <math>W\ =2^{100}=1.27\times 10^{30}</math> and it will take about <math>4.02\times 10^{13}</math> years to calculate a single value of <math>\left\langle M\right\rangle_T</math>.

'''<big>TASK2</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.'''

A montecarlostep(T) function and a statistics() function were written in this section to perform a single step Monte Carlo simulation and calculate average properties of the system, respectively. The output from ILanim.py file is shown below in '''Figure 3'''.
{{multiple image
| image1 = LY_ILanim_out.png

| image2 = LY_ILanim_out_2.png

| direction = vertical
| align = center
| width = 350
| footer = '''Figure 3. The output from ILanim.py and statistics()'''
}}

==Accelerating the code==
'''<big>TASK1</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!'''

'''<big>TASK2</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!'''

In this section, the energy() and statistics() functions were modified using numpy.sum(), numpy.roll() and numpy.multiply() functions. After the modification, the calculation time has been shortened significantly in performing 2000 Monte Carlo steps. The comparison is shown below in '''Table 1'''.

{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 1. The effect of accelerating codes on computational time
! !! before acceleration !! after acceleration
|-
| screen shot || [[File:LY IL MC 2000 original.png|400px]] || [[File:LY IL MC 2000 accelerated final.png|400px]]
|-
| average (s)|| 21.6 || 0.283
|-
| standard error (s)|| 0.162 || 0.00140
|}

==The effect of temperature==
'''<big>TASK1</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.'''

From the top picture in '''Figure 3''', we can see that the first few hundreds step of the Monte Carlo calculation are relaxation steps. Therefore, these steps should be excluded during the calculation of average properties of the system. Some outputs from the ILfinalframe.py file for the 8x8 lattice at different temperatures are shown below in '''Table 2'''. From this table, it is clear that the first 2000 steps should be excluded and a total run-time of 50000 should be appropriate in the calculation of average proprieties of this lattice size.

{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 2. The effect of temperature on the total number of relaxation steps of the 8x8 lattice
! Temperature(K) !! 0.5 !! 1.5 !! 2.0 !! 2.5
|-
| Output Screenshot || [[File:LY IL MC 2000 original.png|300px]] || [[File:LY IL MC 2000 accelerated final.png|300px]] || ||
|}

'''<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 — you will need it later. Save the file as ''8x8.dat'' so that you know which lattice size it came from.'''

From '''Table 2''', it is clear that the first 2000 steps should be excluded and a total run-time of 50000 should be appropriate in the calculation of average proprieties of the 8x8 lattice. The output from the IL.temperaturerange.py file of a 8x8 lattice is shown below in '''Figure 4'''. The temperature range is from 0.2 to 5.0 and step-size is 0.1. Errors bars of average energies and magnetisations at each temperature point are shown as well. In fact these error bars are too small to be seen clearly. However, this does not mean that the calculations are accurate enough......
[[File:LY_88.png|thumb|center|500px|'''Figure 4. The output from the ILtemperature.py file for 8x8 lattice''']]

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

The number of relaxation and equilibration steps needed to calculate reliable average properties depend strongly on the system size. Therefore, several ILfinalframe.py calculations were performed for different lattice sizes to determine the number of relaxation steps and total run-time required. The results are summarised in '''Table 3''' below and these values are used when ILtemperaturerange.py calculations were performed for lattice sizes: 2x2, 4x4, 8x8, 16x16, 32x32.
{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 3. The relaxation and equilibration steps for different lattice sizes
! lattice size !! 2x2 !! 4x4!!8x8!!16x16!!32x32
|-
| number of relaxation steps || 500 || 1000||2000 ||20000 || 75000
|-
| total run time||20000 ||30000 || 50000||150000 ||300000
|}

The ILtemperaturerange.py script was then used to calculate average energies and magnetisations of the following lattice sizes: 2x2, 4x4, 16x16, 32x32 and results are shown in '''Table 4'''.
{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 4. The output of ILtemperaturerange.py for different lattice sizes
! lattice size !! 2x2 !! 4x4
|-
| output || [[File:LY_22.png|400px]] || [[File:LY_44.png|400px]]
|-
| <center>'''lattice size'''</center> ||<center> '''16x16'''</center> ||<center> '''32x32'''</center>
|-
| output || [[File:LY_1616.png|400px]] || [[File:LY_3232.png|400px]]
|}

The energies and magnetisations per spin versus temperature graphs for each lattice size were then plotted together for comparisons, as shown in '''Figure 5'''.
{{multiple image
| image1 = LY_energyallin1.png

| image2 = LY_magallin1.png

| direction = vertical
| align = center
| width = 450
| footer = '''Figure 5. The energies and magnetisations vs. temperature graphs for different lattice sizes shown in one figure'''
}}

From '''Figure 5''', it is clear to see that by increasing the lattice size from 2x2 to 4x4 and from 4x4 to 8x8, both energies and magnetisations per temperature curves change in shapes significantly. In other words, by increasing the lattice size from 2x2 to 8x8, the quality of the calculations has been greatly improved and more details of the behaviors of both energies and magnetisations per temperature curves have been revealed. However, by changing the lattice size from 8x8 to 16x16 and 32x32, both plots in '''Figure 5''' do not change in general features, where as the computational time increases significantly. Therefore, the 8x8 lattice can be seen as the minimum lattice size required to capture the long range fluctuations.

==Determining the heat capacity==
'''<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.'''

Since the heat capacity is defined as <math>C = \frac{\partial E}{\partial T} = \frac{\mathrm{Var}[E]}{k_B T^2}</math>, by extracting the <math>\left\langle E^2\right\rangle</math> and <math>\left\langle E\right\rangle</math> data from the ''.dat'' files generated previously, heat capacity versus temperature curves for each lattice size were computed and plotted in '''Figure 6'''.

[[File:LY_heatcapacityallin1.png|thumb|center|500px|'''Figure 6. Heat capacity versus temperature curves for all the lattice sizes''']]

In '''Figure 6''', as lattice size increases, the temperature as which maximum heat capacity occurs shifts to the left gradually. This is expected due to the finite size effect, which is very common in the simulations of systems bearing long ranged interactions. The trend in '''Figure 6''' is actually in good agreement with the relationship <math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math>, where <math>T_{C, L}</math> is the Curie temperature of an <math>L\times L</math>lattice, <math>T_{C,\infty}</math> is the Curie temperature of an infinite lattice, and <math>A</math> is a constant.

==Locating the Curie temperature==
'''<big>TASK1</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]).'''

'''<big>TASK2</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.'''

'''<big>TASK3</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.'''

'''<big>TASK4</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.'''

Rep:Mod:PYTHONLY2412

2015-02-11T14:55:40Z

Ly2412: /* Determining the heat capacity */

==Introduction to the Ising model==
'''<big>TASK1</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.'''

For a spin in a system of dimension <math>D</math>, there are <math>2D</math> interactions associated with it due to its <math>2D</math> number of neighbouring spins. To minimise the energy of these interactions, this spin and its <math>2D</math> neighbours must have the same sign, i.e. all up or all down. Therefore the minimum energy of these interactions associated with this one spin is <math>E\ =-0.5\times2DJ = -DJ</math>. If the system has total number of spins of <math>N</math>, then the overall lowest energy must be <math>E\ =\ -DNJ</math>. The multiplicity of such a system should be 2 since the whole system can either have all of its spins up or down. The entropy of such a system is therefore <math>S\ =\ K_b\ln 2=9.57\times 10^{-24} JK^{-1} </math>.

'''<big>TASK2</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?'''

For a system of <math>D=3,\ N=1000</math> in its lowest energy configuration, the energy of this configuration is <math>E\ =-3000J</math> and the entropy of this system is <math>S\ =9.57\times 10^{-24} JK^{-1} </math>. If one of the spins change its direction, there are 6 unique interaction energies changes from -1 to +1 and the overall energy changes to <math>E\ =-2988J</math>. In this excited state, the entropy is <math>S\ =\ K_b\ln 2000=1.05\times 10^{-22} JK^{-1} </math> because there are 2000 ways of arranging this system. Therefore the change in energy of this process is <math>\Delta E\ =12J</math> and the change in entropy is <math>\Delta S\ =9.54\times 10^{-23} JK^{-1}</math>.

'''<big>TASK3</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?'''

The magnetisations of the 1D and 2D lattices in figure 1 are both +1. For a Ising lattice with <math>D = 3,\ N=1000</math> at absolute zero, its magnetisation can be either +1000 or -1000.

==Calculating the energy and magnetisation==
In this section, two functions from the IsingLattice.py file, energy() and magnetisation(), were completed to calculate the energy and magnetisation of any 2D Ising lattice and a screenshot of these codes is shown below in '''Figure 1'''.

[[File:LY IL code.png|thumb|center|500px|'''Figure 1. The original code for energy() and magnetisation() functions''']]

The validity of these two functions were checked using a given ILcheck.py file and the output is given below in '''Figure 2'''.

[[File:LY_ILcheck.png|thumb|center|400px|'''Figure 2. The output from the ILcheck.py file''']]

==Introduction to Monte Carlo simulation==

'''<big>TASK1</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, since each spin can be either up or down, the total number of configurations is therefore <math>W\ =2^{100}=1.27\times 10^{30}</math> and it will take about <math>4.02\times 10^{13}</math> years to calculate a single value of <math>\left\langle M\right\rangle_T</math>.

'''<big>TASK2</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.'''

A montecarlostep(T) function and a statistics() function were written in this section to perform a single step Monte Carlo simulation and calculate average properties of the system, respectively. The output from ILanim.py file is shown below in '''Figure 3'''.
{{multiple image
| image1 = LY_ILanim_out.png

| image2 = LY_ILanim_out_2.png

| direction = vertical
| align = center
| width = 350
| footer = '''Figure 3. The output from ILanim.py and statistics()'''
}}

==Accelerating the code==
'''<big>TASK1</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!'''

'''<big>TASK2</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!'''

In this section, the energy() and statistics() functions were modified using numpy.sum(), numpy.roll() and numpy.multiply() functions. After the modification, the calculation time has been shortened significantly in performing 2000 Monte Carlo steps. The comparison is shown below in '''Table 1'''.

{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 1. The effect of accelerating codes on computational time
! !! before acceleration !! after acceleration
|-
| screen shot || [[File:LY IL MC 2000 original.png|400px]] || [[File:LY IL MC 2000 accelerated final.png|400px]]
|-
| average (s)|| 21.6 || 0.283
|-
| standard error (s)|| 0.162 || 0.00140
|}

==The effect of temperature==
'''<big>TASK1</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.'''

From the top picture in '''Figure 3''', we can see that the first few hundreds step of the Monte Carlo calculation are relaxation steps. Therefore, these steps should be excluded during the calculation of average properties of the system. Some outputs from the ILfinalframe.py file for the 8x8 lattice at different temperatures are shown below in '''Table 2'''. From this table, it is clear that the first 2000 steps should be excluded and a total run-time of 50000 should be appropriate in the calculation of average proprieties of this lattice size.

{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 2. The effect of temperature on the total number of relaxation steps of the 8x8 lattice
! Temperature(K) !! 0.5 !! 1.5 !! 2.0 !! 2.5
|-
| Output Screenshot || [[File:LY IL MC 2000 original.png|300px]] || [[File:LY IL MC 2000 accelerated final.png|300px]] || ||
|}

'''<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 — you will need it later. Save the file as ''8x8.dat'' so that you know which lattice size it came from.'''

From '''Table 2''', it is clear that the first 2000 steps should be excluded and a total run-time of 50000 should be appropriate in the calculation of average proprieties of the 8x8 lattice. The output from the IL.temperaturerange.py file of a 8x8 lattice is shown below in '''Figure 4'''. The temperature range is from 0.2 to 5.0 and step-size is 0.1. Errors bars of average energies and magnetisations at each temperature point are shown as well. In fact these error bars are too small to be seen clearly. However, this does not mean that the calculations are accurate enough......
[[File:LY_88.png|thumb|center|500px|'''Figure 4. The output from the ILtemperature.py file for 8x8 lattice''']]

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

The number of relaxation and equilibration steps needed to calculate reliable average properties depend strongly on the system size. Therefore, several ILfinalframe.py calculations were performed for different lattice sizes to determine the number of relaxation steps and total run-time required. The results are summarised in '''Table 3''' below and these values are used when ILtemperaturerange.py calculations were performed for lattice sizes: 2x2, 4x4, 8x8, 16x16, 32x32.
{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 3. The relaxation and equilibration steps for different lattice sizes
! lattice size !! 2x2 !! 4x4!!8x8!!16x16!!32x32
|-
| number of relaxation steps || 500 || 1000||2000 ||20000 || 75000
|-
| total run time||20000 ||30000 || 50000||150000 ||300000
|}

The ILtemperaturerange.py script was then used to calculate average energies and magnetisations of the following lattice sizes: 2x2, 4x4, 16x16, 32x32 and results are shown in '''Table 4'''.
{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 4. The output of ILtemperaturerange.py for different lattice sizes
! lattice size !! 2x2 !! 4x4
|-
| output || [[File:LY_22.png|400px]] || [[File:LY_44.png|400px]]
|-
| <center>'''lattice size'''</center> ||<center> '''16x16'''</center> ||<center> '''32x32'''</center>
|-
| output || [[File:LY_1616.png|400px]] || [[File:LY_3232.png|400px]]
|}

The energies and magnetisations per spin versus temperature graphs for each lattice size were then plotted together for comparisons, as shown in '''Figure 5'''.
{{multiple image
| image1 = LY_energyallin1.png

| image2 = LY_magallin1.png

| direction = vertical
| align = center
| width = 450
| footer = '''Figure 5. The energies and magnetisations vs. temperature graphs for different lattice sizes shown in one figure'''
}}

From '''Figure 5''', it is clear to see that by increasing the lattice size from 2x2 to 4x4 and from 4x4 to 8x8, both energies and magnetisations per temperature curves change in shapes significantly. In other words, by increasing the lattice size from 2x2 to 8x8, the quality of the calculations has been greatly improved and more details of the behaviors of both energies and magnetisations per temperature curves have been revealed. However, by changing the lattice size from 8x8 to 16x16 and 32x32, both plots in '''Figure 5''' do not change in general features, where as the computational time increases significantly. Therefore, the 8x8 lattice can be seen as the minimum lattice size required to capture the long range fluctuations.

==Determining the heat capacity==
'''<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.'''

Since the heat capacity is defined as <math>C = \frac{\partial E}{\partial T} = \frac{\mathrm{Var}[E]}{k_B T^2}</math>, by extracting the <math>\left\langle E^2\right\rangle</math> and <math>\left\langle E\right\rangle</math> data from the ''.dat'' files generated previously, heat capacity versus temperature curves for each lattice size were computed and plotted in '''Figure 6'''.

[[File:LY_heatcapacityallin1.png|thumb|center|500px|'''Figure 6. Heat capacity versus temperature curves for all the lattice sizes''']]

==Locating the Curie temperature==
'''<big>TASK1</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]).'''

'''<big>TASK2</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.'''

'''<big>TASK3</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.'''

'''<big>TASK4</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:LY heatcapacityallin1.png

2015-02-11T14:55:23Z

Ly2412:

Rep:Mod:PYTHONLY2412

2015-02-11T14:52:04Z

Ly2412: /* Determining the heat capacity */

==Introduction to the Ising model==
'''<big>TASK1</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.'''

For a spin in a system of dimension <math>D</math>, there are <math>2D</math> interactions associated with it due to its <math>2D</math> number of neighbouring spins. To minimise the energy of these interactions, this spin and its <math>2D</math> neighbours must have the same sign, i.e. all up or all down. Therefore the minimum energy of these interactions associated with this one spin is <math>E\ =-0.5\times2DJ = -DJ</math>. If the system has total number of spins of <math>N</math>, then the overall lowest energy must be <math>E\ =\ -DNJ</math>. The multiplicity of such a system should be 2 since the whole system can either have all of its spins up or down. The entropy of such a system is therefore <math>S\ =\ K_b\ln 2=9.57\times 10^{-24} JK^{-1} </math>.

'''<big>TASK2</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?'''

For a system of <math>D=3,\ N=1000</math> in its lowest energy configuration, the energy of this configuration is <math>E\ =-3000J</math> and the entropy of this system is <math>S\ =9.57\times 10^{-24} JK^{-1} </math>. If one of the spins change its direction, there are 6 unique interaction energies changes from -1 to +1 and the overall energy changes to <math>E\ =-2988J</math>. In this excited state, the entropy is <math>S\ =\ K_b\ln 2000=1.05\times 10^{-22} JK^{-1} </math> because there are 2000 ways of arranging this system. Therefore the change in energy of this process is <math>\Delta E\ =12J</math> and the change in entropy is <math>\Delta S\ =9.54\times 10^{-23} JK^{-1}</math>.

'''<big>TASK3</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?'''

The magnetisations of the 1D and 2D lattices in figure 1 are both +1. For a Ising lattice with <math>D = 3,\ N=1000</math> at absolute zero, its magnetisation can be either +1000 or -1000.

==Calculating the energy and magnetisation==
In this section, two functions from the IsingLattice.py file, energy() and magnetisation(), were completed to calculate the energy and magnetisation of any 2D Ising lattice and a screenshot of these codes is shown below in '''Figure 1'''.

[[File:LY IL code.png|thumb|center|500px|'''Figure 1. The original code for energy() and magnetisation() functions''']]

The validity of these two functions were checked using a given ILcheck.py file and the output is given below in '''Figure 2'''.

[[File:LY_ILcheck.png|thumb|center|400px|'''Figure 2. The output from the ILcheck.py file''']]

==Introduction to Monte Carlo simulation==

'''<big>TASK1</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, since each spin can be either up or down, the total number of configurations is therefore <math>W\ =2^{100}=1.27\times 10^{30}</math> and it will take about <math>4.02\times 10^{13}</math> years to calculate a single value of <math>\left\langle M\right\rangle_T</math>.

'''<big>TASK2</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.'''

A montecarlostep(T) function and a statistics() function were written in this section to perform a single step Monte Carlo simulation and calculate average properties of the system, respectively. The output from ILanim.py file is shown below in '''Figure 3'''.
{{multiple image
| image1 = LY_ILanim_out.png

| image2 = LY_ILanim_out_2.png

| direction = vertical
| align = center
| width = 350
| footer = '''Figure 3. The output from ILanim.py and statistics()'''
}}

==Accelerating the code==
'''<big>TASK1</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!'''

'''<big>TASK2</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!'''

In this section, the energy() and statistics() functions were modified using numpy.sum(), numpy.roll() and numpy.multiply() functions. After the modification, the calculation time has been shortened significantly in performing 2000 Monte Carlo steps. The comparison is shown below in '''Table 1'''.

{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 1. The effect of accelerating codes on computational time
! !! before acceleration !! after acceleration
|-
| screen shot || [[File:LY IL MC 2000 original.png|400px]] || [[File:LY IL MC 2000 accelerated final.png|400px]]
|-
| average (s)|| 21.6 || 0.283
|-
| standard error (s)|| 0.162 || 0.00140
|}

==The effect of temperature==
'''<big>TASK1</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.'''

From the top picture in '''Figure 3''', we can see that the first few hundreds step of the Monte Carlo calculation are relaxation steps. Therefore, these steps should be excluded during the calculation of average properties of the system. Some outputs from the ILfinalframe.py file for the 8x8 lattice at different temperatures are shown below in '''Table 2'''. From this table, it is clear that the first 2000 steps should be excluded and a total run-time of 50000 should be appropriate in the calculation of average proprieties of this lattice size.

{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 2. The effect of temperature on the total number of relaxation steps of the 8x8 lattice
! Temperature(K) !! 0.5 !! 1.5 !! 2.0 !! 2.5
|-
| Output Screenshot || [[File:LY IL MC 2000 original.png|300px]] || [[File:LY IL MC 2000 accelerated final.png|300px]] || ||
|}

'''<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 — you will need it later. Save the file as ''8x8.dat'' so that you know which lattice size it came from.'''

From '''Table 2''', it is clear that the first 2000 steps should be excluded and a total run-time of 50000 should be appropriate in the calculation of average proprieties of the 8x8 lattice. The output from the IL.temperaturerange.py file of a 8x8 lattice is shown below in '''Figure 4'''. The temperature range is from 0.2 to 5.0 and step-size is 0.1. Errors bars of average energies and magnetisations at each temperature point are shown as well. In fact these error bars are too small to be seen clearly. However, this does not mean that the calculations are accurate enough......
[[File:LY_88.png|thumb|center|500px|'''Figure 4. The output from the ILtemperature.py file for 8x8 lattice''']]

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

The number of relaxation and equilibration steps needed to calculate reliable average properties depend strongly on the system size. Therefore, several ILfinalframe.py calculations were performed for different lattice sizes to determine the number of relaxation steps and total run-time required. The results are summarised in '''Table 3''' below and these values are used when ILtemperaturerange.py calculations were performed for lattice sizes: 2x2, 4x4, 8x8, 16x16, 32x32.
{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 3. The relaxation and equilibration steps for different lattice sizes
! lattice size !! 2x2 !! 4x4!!8x8!!16x16!!32x32
|-
| number of relaxation steps || 500 || 1000||2000 ||20000 || 75000
|-
| total run time||20000 ||30000 || 50000||150000 ||300000
|}

The ILtemperaturerange.py script was then used to calculate average energies and magnetisations of the following lattice sizes: 2x2, 4x4, 16x16, 32x32 and results are shown in '''Table 4'''.
{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 4. The output of ILtemperaturerange.py for different lattice sizes
! lattice size !! 2x2 !! 4x4
|-
| output || [[File:LY_22.png|400px]] || [[File:LY_44.png|400px]]
|-
| <center>'''lattice size'''</center> ||<center> '''16x16'''</center> ||<center> '''32x32'''</center>
|-
| output || [[File:LY_1616.png|400px]] || [[File:LY_3232.png|400px]]
|}

The energies and magnetisations per spin versus temperature graphs for each lattice size were then plotted together for comparisons, as shown in '''Figure 5'''.
{{multiple image
| image1 = LY_energyallin1.png

| image2 = LY_magallin1.png

| direction = vertical
| align = center
| width = 450
| footer = '''Figure 5. The energies and magnetisations vs. temperature graphs for different lattice sizes shown in one figure'''
}}

From '''Figure 5''', it is clear to see that by increasing the lattice size from 2x2 to 4x4 and from 4x4 to 8x8, both energies and magnetisations per temperature curves change in shapes significantly. In other words, by increasing the lattice size from 2x2 to 8x8, the quality of the calculations has been greatly improved and more details of the behaviors of both energies and magnetisations per temperature curves have been revealed. However, by changing the lattice size from 8x8 to 16x16 and 32x32, both plots in '''Figure 5''' do not change in general features, where as the computational time increases significantly. Therefore, the 8x8 lattice can be seen as the minimum lattice size required to capture the long range fluctuations.

==Determining the heat capacity==
'''<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.'''

Since the heat capacity is defined as <math>C = \frac{\partial E}{\partial T} = \frac{\mathrm{Var}[E]}{k_B T^2}</math>, by extracting the <math>\left\langle E^2\right\rangle</math> and <math>\left\langle E\right\rangle</math> data from the ''.dat'' files generated previously, heat capacity versus temperature curves for each lattice size were computed and plotted in '''Figure 6'''.

==Locating the Curie temperature==
'''<big>TASK1</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]).'''

'''<big>TASK2</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.'''

'''<big>TASK3</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.'''

'''<big>TASK4</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.'''

Rep:Mod:PYTHONLY2412

2015-02-11T14:51:45Z

Ly2412: /* Determining the heat capacity */

==Introduction to the Ising model==
'''<big>TASK1</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.'''

For a spin in a system of dimension <math>D</math>, there are <math>2D</math> interactions associated with it due to its <math>2D</math> number of neighbouring spins. To minimise the energy of these interactions, this spin and its <math>2D</math> neighbours must have the same sign, i.e. all up or all down. Therefore the minimum energy of these interactions associated with this one spin is <math>E\ =-0.5\times2DJ = -DJ</math>. If the system has total number of spins of <math>N</math>, then the overall lowest energy must be <math>E\ =\ -DNJ</math>. The multiplicity of such a system should be 2 since the whole system can either have all of its spins up or down. The entropy of such a system is therefore <math>S\ =\ K_b\ln 2=9.57\times 10^{-24} JK^{-1} </math>.

'''<big>TASK2</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?'''

For a system of <math>D=3,\ N=1000</math> in its lowest energy configuration, the energy of this configuration is <math>E\ =-3000J</math> and the entropy of this system is <math>S\ =9.57\times 10^{-24} JK^{-1} </math>. If one of the spins change its direction, there are 6 unique interaction energies changes from -1 to +1 and the overall energy changes to <math>E\ =-2988J</math>. In this excited state, the entropy is <math>S\ =\ K_b\ln 2000=1.05\times 10^{-22} JK^{-1} </math> because there are 2000 ways of arranging this system. Therefore the change in energy of this process is <math>\Delta E\ =12J</math> and the change in entropy is <math>\Delta S\ =9.54\times 10^{-23} JK^{-1}</math>.

'''<big>TASK3</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?'''

The magnetisations of the 1D and 2D lattices in figure 1 are both +1. For a Ising lattice with <math>D = 3,\ N=1000</math> at absolute zero, its magnetisation can be either +1000 or -1000.

==Calculating the energy and magnetisation==
In this section, two functions from the IsingLattice.py file, energy() and magnetisation(), were completed to calculate the energy and magnetisation of any 2D Ising lattice and a screenshot of these codes is shown below in '''Figure 1'''.

[[File:LY IL code.png|thumb|center|500px|'''Figure 1. The original code for energy() and magnetisation() functions''']]

The validity of these two functions were checked using a given ILcheck.py file and the output is given below in '''Figure 2'''.

[[File:LY_ILcheck.png|thumb|center|400px|'''Figure 2. The output from the ILcheck.py file''']]

==Introduction to Monte Carlo simulation==

'''<big>TASK1</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, since each spin can be either up or down, the total number of configurations is therefore <math>W\ =2^{100}=1.27\times 10^{30}</math> and it will take about <math>4.02\times 10^{13}</math> years to calculate a single value of <math>\left\langle M\right\rangle_T</math>.

'''<big>TASK2</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.'''

A montecarlostep(T) function and a statistics() function were written in this section to perform a single step Monte Carlo simulation and calculate average properties of the system, respectively. The output from ILanim.py file is shown below in '''Figure 3'''.
{{multiple image
| image1 = LY_ILanim_out.png

| image2 = LY_ILanim_out_2.png

| direction = vertical
| align = center
| width = 350
| footer = '''Figure 3. The output from ILanim.py and statistics()'''
}}

==Accelerating the code==
'''<big>TASK1</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!'''

'''<big>TASK2</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!'''

In this section, the energy() and statistics() functions were modified using numpy.sum(), numpy.roll() and numpy.multiply() functions. After the modification, the calculation time has been shortened significantly in performing 2000 Monte Carlo steps. The comparison is shown below in '''Table 1'''.

{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 1. The effect of accelerating codes on computational time
! !! before acceleration !! after acceleration
|-
| screen shot || [[File:LY IL MC 2000 original.png|400px]] || [[File:LY IL MC 2000 accelerated final.png|400px]]
|-
| average (s)|| 21.6 || 0.283
|-
| standard error (s)|| 0.162 || 0.00140
|}

==The effect of temperature==
'''<big>TASK1</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.'''

From the top picture in '''Figure 3''', we can see that the first few hundreds step of the Monte Carlo calculation are relaxation steps. Therefore, these steps should be excluded during the calculation of average properties of the system. Some outputs from the ILfinalframe.py file for the 8x8 lattice at different temperatures are shown below in '''Table 2'''. From this table, it is clear that the first 2000 steps should be excluded and a total run-time of 50000 should be appropriate in the calculation of average proprieties of this lattice size.

{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 2. The effect of temperature on the total number of relaxation steps of the 8x8 lattice
! Temperature(K) !! 0.5 !! 1.5 !! 2.0 !! 2.5
|-
| Output Screenshot || [[File:LY IL MC 2000 original.png|300px]] || [[File:LY IL MC 2000 accelerated final.png|300px]] || ||
|}

'''<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 — you will need it later. Save the file as ''8x8.dat'' so that you know which lattice size it came from.'''

From '''Table 2''', it is clear that the first 2000 steps should be excluded and a total run-time of 50000 should be appropriate in the calculation of average proprieties of the 8x8 lattice. The output from the IL.temperaturerange.py file of a 8x8 lattice is shown below in '''Figure 4'''. The temperature range is from 0.2 to 5.0 and step-size is 0.1. Errors bars of average energies and magnetisations at each temperature point are shown as well. In fact these error bars are too small to be seen clearly. However, this does not mean that the calculations are accurate enough......
[[File:LY_88.png|thumb|center|500px|'''Figure 4. The output from the ILtemperature.py file for 8x8 lattice''']]

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

The number of relaxation and equilibration steps needed to calculate reliable average properties depend strongly on the system size. Therefore, several ILfinalframe.py calculations were performed for different lattice sizes to determine the number of relaxation steps and total run-time required. The results are summarised in '''Table 3''' below and these values are used when ILtemperaturerange.py calculations were performed for lattice sizes: 2x2, 4x4, 8x8, 16x16, 32x32.
{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 3. The relaxation and equilibration steps for different lattice sizes
! lattice size !! 2x2 !! 4x4!!8x8!!16x16!!32x32
|-
| number of relaxation steps || 500 || 1000||2000 ||20000 || 75000
|-
| total run time||20000 ||30000 || 50000||150000 ||300000
|}

The ILtemperaturerange.py script was then used to calculate average energies and magnetisations of the following lattice sizes: 2x2, 4x4, 16x16, 32x32 and results are shown in '''Table 4'''.
{| class="wikitable" style="margin: 1em auto 1em auto;"
|+ Table 4. The output of ILtemperaturerange.py for different lattice sizes
! lattice size !! 2x2 !! 4x4
|-
| output || [[File:LY_22.png|400px]] || [[File:LY_44.png|400px]]
|-
| <center>'''lattice size'''</center> ||<center> '''16x16'''</center> ||<center> '''32x32'''</center>
|-
| output || [[File:LY_1616.png|400px]] || [[File:LY_3232.png|400px]]
|}

The energies and magnetisations per spin versus temperature graphs for each lattice size were then plotted together for comparisons, as shown in '''Figure 5'''.
{{multiple image
| image1 = LY_energyallin1.png

| image2 = LY_magallin1.png

| direction = vertical
| align = center
| width = 450
| footer = '''Figure 5. The energies and magnetisations vs. temperature graphs for different lattice sizes shown in one figure'''
}}

From '''Figure 5''', it is clear to see that by increasing the lattice size from 2x2 to 4x4 and from 4x4 to 8x8, both energies and magnetisations per temperature curves change in shapes significantly. In other words, by increasing the lattice size from 2x2 to 8x8, the quality of the calculations has been greatly improved and more details of the behaviors of both energies and magnetisations per temperature curves have been revealed. However, by changing the lattice size from 8x8 to 16x16 and 32x32, both plots in '''Figure 5''' do not change in general features, where as the computational time increases significantly. Therefore, the 8x8 lattice can be seen as the minimum lattice size required to capture the long range fluctuations.

==Determining the heat capacity==
'''<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.'''

Since the heat capacity is defined as <math>C = \frac{\partial E}{\partial T} = \frac{\mathrm{Var}[E]}{k_B T^2}</math>, by extracting the <math>\left\langle E^2\right\rangle</math> and <math>\left\langle E\right\rangle</math> data from the ''.dat''' files generated previously, heat capacity versus temperature curves for each lattice size were computed and plotted in '''Figure 6'''.

==Locating the Curie temperature==
'''<big>TASK1</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]).'''

'''<big>TASK2</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.'''

'''<big>TASK3</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.'''

'''<big>TASK4</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.'''

Rep:Mod:PYTHONLY2412

2015-02-11T14:45:19Z

Ly2412: /* The effect of system size */