ChemWiki - User contributions [en]

Programming a 2D Ising Model

2026-05-03T14:25:29Z

Gmallia:

'''<big>This is the Ising Model experiment (aka programming for simple simulations).</big>'''

'''<big>This is a core experiment for year 3 students on the Chemistry with Molecular Physics degree.</big>'''

= Programming a 2D Ising Model/ Main page =

==Introduction==

In this exercise, you are going to use the Python that you learned in the first year to write a code to perform Monte-Carlo simulations of the 2D [http://en.wikipedia.org/wiki/Ising_model Ising model], a set of spins on a lattice which is used to model ferromagnetic behaviour, and also to analyse the results of the simulation to find the heat capacity of the system and the Curie temperature — the temperature below which the system is able to maintain a spontaneous magnetisation.

A more extensive introduction to the Ising Model will be given in the next page. This page gives some practical information about how the experiment.

==Assessment==
At the end of this experiment you must submit a report, writing up your findings. Each section of the experiment has a number of tasks that you should complete, labelled '''<big>TASK</big>'''. If this is a mathematical exercise, your report should contain a short summary of the solution. If it is a graphical exercise, you should include the relevant image. Your report should explain briefly what you did in each stage of the experiment, and what your findings were.

The report should be submitted as a PDF file to Turnitin, created with your favourite word processor or LaTeX, and should contain the annotated source code of the Python scripts that your write during the experiment. For clarity, code included in the report should be written in a monospaced font (e.g. Courier), which in LaTeX can be easily achieved by placing your code inside a <code>\texttt{}</code> command. '''The code, except where relevant to discussion, can be placed in an appendix.'''

As well as including relevant source code in your report, you should also submit a bundle with all your scripts and analysis code as a separate Blackboard submission.

===On the use of AI tools===
The code you write for this experiment is assessed, and as per the [https://bb.imperial.ac.uk/bbcswebdav/pid-3344481-dt-content-rid-21115105_1/xid-21115105_1 department's guidelines on the use of AI], it is not acceptable to use generative AI tools to write computer code for this assignment.

Writing computer programs yourself enhances your logical thinking and problem solving skills. Using generative AI tools to write your code makes you dumb. Don't choose to be dumb, don't use AI to write code for you.

==Getting Help==

The Graduate Teaching Assistants demonstrating this experiment are Fionn Carman and Lei Li. Demonstrators will be available in room 232A at the following times:

{| class="wikitable"
|+ Demonstrator timetable
! !! Monday !! Tuesday !! Wednesday !! Thusrsday !! Friday
|-
| Morning || 10:00 - Experiment intro 10:00-13:00 - Li || 10:00-13:00 - Fionn|| || 10:00-12:00 - João || 10:00-13:00 - Fionn
|-
| Afternoon || 14:00-17:00 - Li|| 14:00-17:00 - Fionn|| || || 14:00-17:00 - Fionn
|}

If you have questions outside these times, please use the Discussion Board on Blackboard.

João and Giuseppe are the members of academic staff responsible for this experiment. You are welcome to address your questions and comments to them.

==Structure of this Experiment==

This experimental manual has been broken up into a number of subsections. Direct links to each of them may be found below. You should attempt them in order, and you should complete all of them to finish the experiment.

# [[Programming_a_2D_Ising_Model/Introduction_to_the_Ising_model|Introduction to the Ising model]]
# [[Programming_a_2D_Ising_Model/Calculating the energy and magnetisation|Calculating the energy and magnetisation]]
# [[Programming_a_2D_Ising_Model/Introduction to Monte Carlo simulation|Introduction to the Monte Carlo simulation]]
# [[Programming_a_2D_Ising_Model/Accelerating the code|Accelerating the code]]
# [[Programming_a_2D_Ising_Model/The effect of temperature|The effect of temperature]]
# [[Programming_a_2D_Ising_Model/The effect of system size|The effect of system size]]
# [[Programming_a_2D_Ising_Model/Determining the heat capacity|Determining the heat capacity]]
# [[Programming_a_2D_Ising_Model/Locating the Curie temperature|Locating the Curie temperature]]

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]]</big>'''

Programming a 2D Ising Model/Locating the Curie temperature

2026-05-03T14:08:30Z

Gmallia:

'''<big>This is the eighth (and final) section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Determining the heat capacity|7. Determining the heat capacity]]</big>'''

= 8. Locating the Curie temperature =
You should have seen in the previous section that the heat capacity becomes strongly peaked in the vicinity of the critical temperature and that the peak became increasingly sharply peaked as the system size was increased — in fact, Onsager proved that in an infinite system the heat capacity should diverge at ''T'' = ''TC''. In our finite systems, however, not only does the heat capacity not diverge, the Curie temperature changes with system size! This is known as a ''finite size effect''.

It can be shown, however, that the temperature at which the heat capacity has its maximum must scale according to <math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math>, where ''TC,L'' is the Curie temperature of an ''L''×''L'' lattice, ''TC,∞'' is the Curie temperature of an infinite lattice, and ''A'' is a constant in which we are not especially interested. This scaling equation comes from expanding in the limit of large ''L''.

==== '''<big>TASK 8a</big>:''' ====
'''Reference data has been generated for longer simulations than would be possible on the college computers in Python. Each file contains six columns: ''T'', ''E'', ''E''2, ''M'', ''M''2, ''C'' (NOTE: these final five quantities are normalised per spin, unlike your Python data), and you can read them with the NumPy <code>loadtxt</code> function as before. For each lattice size, plot the reference data against your data. For ''one'' lattice size, save a SVG or PNG image of this comparison and add it to your report — add a legend to the graph to label which is which (check if you need a [https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.legend.html reminder about use of legend]).'''

==Polynomial fitting==

To find the temperature at which the heat capacity and susceptibilities have their maxima, we are going to fit a polynomial to the data in the critical region. NumPy provides the useful [http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html <code>polyfit</code>] and [http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyval.html#numpy.polyval <code>polyval</code>] functions for this purpose. The following code is written for the heat capacities - try to also repeat this for the susceptibility! The usage is as follows:

<pre>
data = np.loadtxt("...") #assume data is now a 2D array containing two columns, T and C
T = data[:,0] #get the first column
C = data[:,1] # get the second column

#first we fit the polynomial to the data
fit = np.polyfit(T, C, 3) # fit a third order polynomial

#now we generate interpolated values of the fitted polynomial over the range of our function
T_min = np.min(T)
T_max = np.max(T)
T_range = np.linspace(T_min, T_max, 1000) #generate 1000 evenly spaced points between T_min and T_max
fitted_C_values = np.polyval(fit, T_range) # use the fit object to generate the corresponding values of C
</pre>

==== '''<big>TASK 8b</big>:''' ====
'''Write a script to read the data from a particular file, and plot ''C'' and ''χ'' 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 SVG or PNG image of an example fit to your report.'''

===Fitting in a particular temperature range===

Rather than fit to all of the data, we might want to fit in only a particular range. NumPy provides a very powerful way to index arrays based on certain conditions. For example, if we want to extract only those data points which are between a particular ''T''min and ''T''max, we can use the following syntax:

<pre>
data = np.loadtxt("...") #assume data is now a 2D array containing two columns, T and C
T = data[:,0] #get the first column
C = data[:,1] # get the second column

Tmin = 0.5 #for example
Tmax = 2.0 #for example

selection = np.logical_and(T > Tmin, T < Tmax) #choose only those rows where both conditions are true
peak_T_values = T[selection]
peak_C_values = C[selection]
</pre>

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

===Finding the peak in C===

Your fitting script should now generate two variables: <code>peak_T_range</code>, containing 1000 equally spaced temperature values between ''T''min and ''T''max (whatever you chose those values to be), and fitted_''C''_values, containing the fitted heat capacity at each of those points. Use the NumPy <code>max</code> function to find the maximum in ''C''. If you store the maximum value of ''C'' in the variable <code>Cmax</code>, you can use the following notation to find the corresponding temperature:

<pre>
Cmax = np.max(...)
Tmax = peak_T_range[fitted_C_values == Cmax]
</pre>

==== '''<big>TASK 8d</big>:''' ====
'''Find the temperature at which the maximum in ''C'' occurs for each datafile that you were given or generated with Python. Make a text file containing two columns: the lattice side length (2,4,8, etc.), and the temperature at which ''C'' is a maximum. This is your estimate of ''TC'' for that side length. Make a plot that uses the scaling relation given above to determine ''TC,∞''. By doing a little research, 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 SVG or PNG image of this final graph to your report, and discuss briefly what you think the major sources of error are in your estimate.'''

'''<big>This is the eighth (and final) section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Determining the heat capacity|7. Determining the heat capacity]]</big>'''

Programming a 2D Ising Model/Locating the Curie temperature

2026-05-03T14:07:39Z

Gmallia:

'''<big>This is the eighth (and final) section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Determining the heat capacity|7. Determining the heat capacity]]</big>'''

= 8. Locating the Curie temperature =
You should have seen in the previous section that the heat capacity becomes strongly peaked in the vicinity of the critical temperature and that the peak became increasingly sharply peaked as the system size was increased — in fact, Onsager proved that in an infinite system the heat capacity should diverge at ''T'' = ''TC''. In our finite systems, however, not only does the heat capacity not diverge, the Curie temperature changes with system size! This is known as a ''finite size effect''.

It can be shown, however, that the temperature at which the heat capacity has its maximum must scale according to <math>T_{C, L} = \frac{A}{L} + T_{C,\infty}</math>, where ''TC,L'' is the Curie temperature of an ''L''×''L'' lattice, ''TC,∞'' is the Curie temperature of an infinite lattice, and ''A'' is a constant in which we are not especially interested. This scaling equation comes from expanding in the limit of large ''L''.

==== '''<big>TASK 8a</big>:''' ====
'''Reference data has been generated for longer simulations than would be possible on the college computers in Python. Each file contains six columns: ''T'', ''E'', ''E''2, ''M'', ''M''2, ''C'' (NOTE: these final five quantities are normalised per spin, unlike your Python data), and you can read them with the NumPy <code>loadtxt</code> function as before. For each lattice size, plot the reference data against your data. For ''one'' lattice size, save a SVG or PNG image of this comparison and add it to your report — add a legend to the graph to label which is which (check if you need a [https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.legend.html reminder about use of legend]).'''

==Polynomial fitting==

To find the temperature at which the heat capacity and susceptibilities have their maxima, we are going to fit a polynomial to the data in the critical region. NumPy provides the useful [http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html <code>polyfit</code>] and [http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyval.html#numpy.polyval <code>polyval</code>] functions for this purpose. The following code is written for the heat capacities - try to also repeat this for the susceptibility! The usage is as follows:

<pre>
data = np.loadtxt("...") #assume data is now a 2D array containing two columns, T and C
T = data[:,0] #get the first column
C = data[:,1] # get the second column

#first we fit the polynomial to the data
fit = np.polyfit(T, C, 3) # fit a third order polynomial

#now we generate interpolated values of the fitted polynomial over the range of our function
T_min = np.min(T)
T_max = np.max(T)
T_range = np.linspace(T_min, T_max, 1000) #generate 1000 evenly spaced points between T_min and T_max
fitted_C_values = np.polyval(fit, T_range) # use the fit object to generate the corresponding values of C
</pre>

==== '''<big>TASK 8b</big>:''' ====

==== '''Write a script to read the data from a particular file, and plot ''C'' and ''χ'' 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 SVG or PNG image of an example fit to your report.''' ====

===Fitting in a particular temperature range===

Rather than fit to all of the data, we might want to fit in only a particular range. NumPy provides a very powerful way to index arrays based on certain conditions. For example, if we want to extract only those data points which are between a particular ''T''min and ''T''max, we can use the following syntax:

<pre>
data = np.loadtxt("...") #assume data is now a 2D array containing two columns, T and C
T = data[:,0] #get the first column
C = data[:,1] # get the second column

Tmin = 0.5 #for example
Tmax = 2.0 #for example

selection = np.logical_and(T > Tmin, T < Tmax) #choose only those rows where both conditions are true
peak_T_values = T[selection]
peak_C_values = C[selection]
</pre>

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

===Finding the peak in C===

Your fitting script should now generate two variables: <code>peak_T_range</code>, containing 1000 equally spaced temperature values between ''T''min and ''T''max (whatever you chose those values to be), and fitted_''C''_values, containing the fitted heat capacity at each of those points. Use the NumPy <code>max</code> function to find the maximum in ''C''. If you store the maximum value of ''C'' in the variable <code>Cmax</code>, you can use the following notation to find the corresponding temperature:

<pre>
Cmax = np.max(...)
Tmax = peak_T_range[fitted_C_values == Cmax]
</pre>

==== '''<big>TASK 8d</big>:''' ====
'''Find the temperature at which the maximum in ''C'' occurs for each datafile that you were given or generated with Python. Make a text file containing two columns: the lattice side length (2,4,8, etc.), and the temperature at which ''C'' is a maximum. This is your estimate of ''TC'' for that side length. Make a plot that uses the scaling relation given above to determine ''TC,∞''. By doing a little research, 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 SVG or PNG image of this final graph to your report, and discuss briefly what you think the major sources of error are in your estimate.'''

'''<big>This is the eighth (and final) section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Determining the heat capacity|7. Determining the heat capacity]]</big>'''

Programming a 2D Ising Model

2026-05-03T14:04:35Z

Gmallia: /* Structure of this Experiment */

'''<big>This is the Ising Model experiment (aka programming for simple simulations).</big>'''

'''<big>This is a core experiment for year 3 students on the Chemistry with Molecular Physics degree.</big>'''

= Programming a 2D Ising Model/ Main page =

==Introduction==

In this exercise, you are going to use the Python that you learned in the first year to write a code to perform Monte-Carlo simulations of the 2D [http://en.wikipedia.org/wiki/Ising_model Ising model], a set of spins on a lattice which is used to model ferromagnetic behaviour, and also to analyse the results of the simulation to find the heat capacity of the system and the Curie temperature — the temperature below which the system is able to maintain a spontaneous magnetisation.

A more extensive introduction to the Ising Model will be given in the next page. This page gives some practical information about how the experiment.

==Assessment==
At the end of this experiment you must submit a report, writing up your findings. Each section of the experiment has a number of tasks that you should complete, labelled '''<big>TASK</big>'''. If this is a mathematical exercise, your report should contain a short summary of the solution. If it is a graphical exercise, you should include the relevant image. Your report should explain briefly what you did in each stage of the experiment, and what your findings were.

The report should be submitted as a PDF file to Turnitin, created with your favourite word processor or LaTeX, and should contain the annotated source code of the Python scripts that your write during the experiment. For clarity, code included in the report should be written in a monospaced font (e.g. Courier), which in LaTeX can be easily achieved by placing your code inside a <code>\texttt{}</code> command. '''The code, except where relevant to discussion, can be placed in an appendix.'''

As well as including relevant source code in your report, you should also submit a bundle with all your scripts and analysis code as a separate Blackboard submission.

===On the use of AI tools===
The code you write for this experiment is assessed, and as per the [https://bb.imperial.ac.uk/bbcswebdav/pid-3344481-dt-content-rid-21115105_1/xid-21115105_1 department's guidelines on the use of AI], it is not acceptable to use generative AI tools to write computer code for this assignment.

Writing computer programs yourself enhances your logical thinking and problem solving skills. Using generative AI tools to write your code makes you dumb. Don't choose to be dumb, don't use AI to write code for you.

==Getting Help==

The Graduate Teaching Assistants demonstrating this experiment are Fionn Carman and Lei Li. Demonstrators will be available in room 232A at the following times:

{| class="wikitable"
|+ Demonstrator timetable
! !! Monday !! Tuesday !! Wednesday !! Thusrsday !! Friday
|-
| Morning || 10:00 - Experiment intro 10:00-13:00 - Li || 10:00-13:00 - Fionn|| || 10:00-12:00 - João || 10:00-13:00 - Fionn
|-
| Afternoon || 14:00-17:00 - Li|| 14:00-17:00 - Fionn|| || || 14:00-17:00 - Fionn
|}

If you have questions outside these times, please use the Discussion Board on Blackboard.

João and Giuseppe are the members of academic staff responsible for this experiment. You are welcome to address your questions and comments to them.

==Structure of this Experiment==

This experimental manual has been broken up into a number of subsections. Direct links to each of them may be found below. You should attempt them in order, and you should complete all of them to finish the experiment.

# [[Programming_a_2D_Ising_Model/Introduction_to_the_Ising_model|Introduction to the Ising model]]
# [[Programming_a_2D_Ising_Model/Calculating the energy and magnetisation|Calculating the energy and magnetisation]]
# [[Programming_a_2D_Ising_Model/Introduction to Monte Carlo simulation|Introduction to the Monte Carlo simulation]]
# [[Programming_a_2D_Ising_Model/Accelerating the code|Accelerating the code]]
# [[Programming_a_2D_Ising_Model/The effect of temperature|The effect of temperature]]
# [[Programming_a_2D_Ising_Model/The effect of system size|The effect of system size]]
# [[Programming_a_2D_Ising_Model/Determining the heat capacity|Determining the heat capacity]]
# [[Programming_a_2D_Ising_Model/Locating the Curie temperature|Locating the Curie temperature]]

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]]</big>'''

Programming a 2D Ising Model

2026-05-03T14:02:53Z

Gmallia:

'''<big>This is the Ising Model experiment (aka programming for simple simulations).</big>'''

'''<big>This is a core experiment for year 3 students on the Chemistry with Molecular Physics degree.</big>'''

= Programming a 2D Ising Model/ Main page =

==Introduction==

In this exercise, you are going to use the Python that you learned in the first year to write a code to perform Monte-Carlo simulations of the 2D [http://en.wikipedia.org/wiki/Ising_model Ising model], a set of spins on a lattice which is used to model ferromagnetic behaviour, and also to analyse the results of the simulation to find the heat capacity of the system and the Curie temperature — the temperature below which the system is able to maintain a spontaneous magnetisation.

A more extensive introduction to the Ising Model will be given in the next page. This page gives some practical information about how the experiment.

==Assessment==
At the end of this experiment you must submit a report, writing up your findings. Each section of the experiment has a number of tasks that you should complete, labelled '''<big>TASK</big>'''. If this is a mathematical exercise, your report should contain a short summary of the solution. If it is a graphical exercise, you should include the relevant image. Your report should explain briefly what you did in each stage of the experiment, and what your findings were.

The report should be submitted as a PDF file to Turnitin, created with your favourite word processor or LaTeX, and should contain the annotated source code of the Python scripts that your write during the experiment. For clarity, code included in the report should be written in a monospaced font (e.g. Courier), which in LaTeX can be easily achieved by placing your code inside a <code>\texttt{}</code> command. '''The code, except where relevant to discussion, can be placed in an appendix.'''

As well as including relevant source code in your report, you should also submit a bundle with all your scripts and analysis code as a separate Blackboard submission.

===On the use of AI tools===
The code you write for this experiment is assessed, and as per the [https://bb.imperial.ac.uk/bbcswebdav/pid-3344481-dt-content-rid-21115105_1/xid-21115105_1 department's guidelines on the use of AI], it is not acceptable to use generative AI tools to write computer code for this assignment.

Writing computer programs yourself enhances your logical thinking and problem solving skills. Using generative AI tools to write your code makes you dumb. Don't choose to be dumb, don't use AI to write code for you.

==Getting Help==

The Graduate Teaching Assistants demonstrating this experiment are Fionn Carman and Lei Li. Demonstrators will be available in room 232A at the following times:

{| class="wikitable"
|+ Demonstrator timetable
! !! Monday !! Tuesday !! Wednesday !! Thusrsday !! Friday
|-
| Morning || 10:00 - Experiment intro 10:00-13:00 - Li || 10:00-13:00 - Fionn|| || 10:00-12:00 - João || 10:00-13:00 - Fionn
|-
| Afternoon || 14:00-17:00 - Li|| 14:00-17:00 - Fionn|| || || 14:00-17:00 - Fionn
|}

If you have questions outside these times, please use the Discussion Board on Blackboard.

João and Giuseppe are the members of academic staff responsible for this experiment. You are welcome to address your questions and comments to them.

==Structure of this Experiment==

This experimental manual has been broken up into a number of subsections. Direct links to each of them may be found below. You should attempt them in order, and you should complete all of them to finish the experiment.

# [[Programming_a_2D_Ising_Model/Introduction_to_the_Ising_model|Introduction to the Ising model]]
# [[Programming_a_2D_Ising_Model/Calculating the energy and magnetisation|Calculating the energy and magnetisation]]
# [[Programming_a_2D_Ising_Model/Introduction to Monte Carlo simulation|Introduction to the Monte Carlo simulation]]
# [[Programming_a_2D_Ising_Model/Accelerating the code|Accelerating the code]]
# [[Programming_a_2D_Ising_Model/The effect of temperature|The effect of temperature]]
# [[Programming_a_2D_Ising_Model/The effect of system size|The effect of system size]]
# [[Programming_a_2D_Ising_Model/Determining the heat capacity|Determining the heat capacity]]
# [[Programming_a_2D_Ising_Model/Locating the Curie temperature|Locating the Curie temperature]]

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]]</big>'''

Programming a 2D Ising Model

2026-05-03T14:01:40Z

Gmallia:

'''<big>This is the Ising Model experiment (aka programming for simple simulations).</big>'''

'''<big>This is a core experiment for year 3 students on the Chemistry with Molecular Physics degree.</big>'''

= Programming a 2D Ising Model/ Main page =

==Introduction==

In this exercise, you are going to use the Python that you learned in the first year to write a code to perform Monte-Carlo simulations of the 2D [http://en.wikipedia.org/wiki/Ising_model Ising model], a set of spins on a lattice which is used to model ferromagnetic behaviour, and also to analyse the results of the simulation to find the heat capacity of the system and the Curie temperature — the temperature below which the system is able to maintain a spontaneous magnetisation.

A more extensive introduction to the Ising Model will be given in the next page. This page gives some practical information about how the experiment.

==Assessment==
At the end of this experiment you must submit a report, writing up your findings. Each section of the experiment has a number of tasks that you should complete, labelled '''<big>TASK</big>'''. If this is a mathematical exercise, your report should contain a short summary of the solution. If it is a graphical exercise, you should include the relevant image. Your report should explain briefly what you did in each stage of the experiment, and what your findings were.

The report should be submitted as a PDF file to Turnitin, created with your favourite word processor or LaTeX, and should contain the annotated source code of the Python scripts that your write during the experiment. For clarity, code included in the report should be written in a monospaced font (e.g. Courier), which in LaTeX can be easily achieved by placing your code inside a <code>\texttt{}</code> command. '''The code, except where relevant to discussion, can be placed in an appendix.'''

As well as including relevant source code in your report, you should also submit a bundle with all your scripts and analysis code as a separate Blackboard submission.

===On the use of AI tools===
The code you write for this experiment is assessed, and as per the [https://bb.imperial.ac.uk/bbcswebdav/pid-3344481-dt-content-rid-21115105_1/xid-21115105_1 department's guidelines on the use of AI], it is not acceptable to use generative AI tools to write computer code for this assignment.

Writing computer programs yourself enhances your logical thinking and problem solving skills. Using generative AI tools to write your code makes you dumb. Don't choose to be dumb, don't use AI to write code for you.

==Getting Help==

The Graduate Teaching Assistants demonstrating this experiment are Fionn Carman and Lei Li. Demonstrators will be available in room 232A at the following times:

{| class="wikitable"
|+ Demonstrator timetable
! !! Monday !! Tuesday !! Wednesday !! Thusrsday !! Friday
|-
| Morning || 10:00 - Experiment intro 10:00-13:00 - Li || 10:00-13:00 - Fionn|| || 10:00-12:00 - João || 10:00-13:00 - Fionn
|-
| Afternoon || 14:00-17:00 - Li|| 14:00-17:00 - Fionn|| || || 14:00-17:00 - Fionn
|}

If you have questions outside these times, please use the Discussion Board on Blackboard.

João and Giuseppe are the members of academic staff responsible for this experiment. You are welcome to address your questions and comments to them.

==Structure of this Experiment==

This experimental manual has been broken up into a number of subsections. Direct links to each of them may be found below. You should attempt them in order, and you should complete all of them to finish the experiment.

# [[Programming_a_2D_Ising_Model/Introduction_to_the_Ising_model|Introduction to the Ising model]]
# [[Programming_a_2D_Ising_Model/Calculating the energy and magnetisation|Calculating the energy and magnetisation]]
# [[Programming_a_2D_Ising_Model/Introduction to Monte Carlo simulation|Introduction to the Monte Carlo simulation]]
# [[Programming_a_2D_Ising_Model/Accelerating the code|Accelerating the code]]
# [[Programming_a_2D_Ising_Model/The effect of temperature|The effect of temperature]]
# [[Programming_a_2D_Ising_Model/The effect of system size|The effect of system size]]
# [[Programming_a_2D_Ising_Model/Determining the heat capacity|Determining the heat capacity]]
# [[Programming_a_2D_Ising_Model/Locating the Curie temperature|Locating the Curie temperature]]

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]]</big>'''

Programming a 2D Ising Model/Determining the heat capacity

2026-05-03T14:00:31Z

Gmallia:

'''<big>This is the seventh section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/The effect of system size|6. The effect of system size]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Locating the Curie temperature|8. Locating the Curie temperature]]</big>'''

==Calculating the heat capacity==

As we have seen, increasing the temperature above ''TC'' induces a phase transition — the magnetisation of the system rapidly drops, but it can be hard to use this information to pinpoint the Curie temperature itself. As well as demonstrating the closed form solution to the partition function that we mentioned in the introduction, Lars Onsager also demonstrated that the heat capacity of the 2D Ising model should become very strongly peaked at the phase transition temperature (in fact, when ''T'' = ''TC'' exactly, ''C'' diverges).

==== '''<big>TASK 7a</big>:''' ====
'''By definition, <math>C = \frac{\partial \left\langle E\right\rangle}{\partial T}</math>. From this, show that <math>C = \frac{\mathrm{Var}[E]}{k_B T^2}</math>, where <math>\mathrm{Var}[E]</math> is the variance in <math>E</math>.'''

Using this relation, and the data that you generated in the previous sections, you can now determine the heat capacity of your lattice as a function of temperature and system size. in a similar way, we can also define another quantity called the magnetic susceptibility, ''χ'', which depends on the average properties of the order parameter of the system:

<math>\chi = \beta(\left\langle M^2\right\rangle - \left\langle M\right\rangle^2)</math>

which will behave in a similar way to the heat capacity. Here <math>\beta = \frac{1}{k_BT}</math> is the inverse of the thermodynamic temperature, often used in statistical thermodynamics.

==== '''<big>TASK 7b</big>:''' ====
'''Write a Python script to make a plot showing the heat capacity and susceptibility 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, Var[''X''], the mean of its square &lang;''X''2&rang;, and its squared mean &lang;''X''&rang;2. You may find that the data around the peak is very noisy — this is normal, and is a result of being in the critical region. As before, use the plot controls to save a SVG or PNG image of your plot and attach this to the report.'''

'''<big>This is the seventh section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/The effect of system size|6. The effect of system size]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Locating the Curie temperature|8. Locating the Curie temperature]]</big>'''

Programming a 2D Ising Model/The effect of system size

2026-05-03T13:53:53Z

Gmallia:

'''<big>This is the sixth section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/The effect of temperature|5. The effect of temperature]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Determining the heat capacity|7. Determining the heat capacity]]</big>'''

= 6. The effect of system size =

==Scaling the System Size==

Your plots from the previous section showed how the magnetisation and energy vary with temperature, and you should be able to see the onset of the phase transition. In this region, the energetic and entropic driving forces are of almost equal importance, and large fluctuations in the state of the system can take place. In fact, it is a characteristic of phase transitions that fluctuations within a system start to take place over very long ranges. This is a big problem in our simulation — so far, we only have 64 spins (or "atoms", if you prefer), so our system may not be big enough for these long ranged fluctuations to be correctly modelled.

==== '''<big>TASK 6</big>:''' ====
'''Repeat the final task of the previous section for the following lattice sizes: 2×2, 4×4, 8×8, 16×16, 32×32. Make sure that you name each datafile that you 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 <code>loadtxt</code> function is the reverse of the <code>savetxt</code> 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 an SVG or 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 correlations?'''

'''<big>This is the sixth section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/The effect of temperature|5. The effect of temperature]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Determining the heat capacity|7. Determining the heat capacity]]</big>'''

Programming a 2D Ising Model/Accelerating the code

2026-05-03T13:51:45Z

Gmallia: /* 4. Accelerating the code */

'''<big>This is the fourth section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of temperature|5. The effect of temperature]]</big>'''

= 4. Accelerating the code =
So far, your energy() and magnetisation() functions probably use Python loops to go over every spin in your system (and a double loop for the energy). One of the drawbacks of Python's flexibility is a performance penalty — though developing code in Python tends to be
faster than other compiled programming languages, the resulting program is much, much slower.

Fortunately, the NumPy module provides an interface to a number of functions written to give high performance.
These are written in the C programming language under the hood and take numpy arrays, which are contiguous arrays of a single data type,
unlike Python lists, which may contain mixes of data types.

In this section, we are going to use these functions to improve the performance of our simulation.

When trying to optimise the performance of something, you should always measure first, rather than making assumptions. Although there are some very obvious candidates in our program, you should always measure properly.

The <code>ILtimetrial.py</code> script measures how fast it takes for 2000 iterations of your code to run.

The <code>ILProfiling.py</code> script is the same, except uses the python standard library package [https://docs.python.org/3/library/profile.html cProfile] to sample what functions are being called and give you information
about the time spent in each. Analysing the output of the script, you will see that by far the majority of the time is spent in the energy method, followed by the magnetisation one. This makes sense, as each Monte-Carlo step can call the energy function multiple times.

(If you want to profile your code in more detail, please seek advice from a demonstrator.)

In order to improve the performance of this, we need to:
* Run the energy call fewer times per step.
* Improve the performance of the energy function.

Fortunately, the NumPy module provides an interface to a number of functions written to give high performance. In this section, we are going to use these functions to improve the performance of our simulation.

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

==== '''<big>TASK 4c</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! What makes this new version faster?'''

==== '''<big>TASK 4d:</big>''' ====
'''<big>Making the code even faster </big>: Consider trying to make the Monte Carlo step function substantially faster by calculating the difference in energy directly, in a way that avoids looping over the whole lattice. Think back how you answered some of the introduction questions. You can make it so the energy and magnetisation only needs to be calculated explicitly at the beginning of the simulation.'''

'''<big>This is the fourth section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of temperature|5. The effect of temperature]]</big>'''

Programming a 2D Ising Model/Introduction to Monte Carlo simulation

2026-05-03T13:51:06Z

Gmallia: /* 3. Introduction to Monte Carlo simulation */

'''<big>This is the third section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]]</big>'''

= 3. Introduction to Monte Carlo simulation =

==Average Energy and Magnetisation==

Consider again the expressions for the average energy and magnetisation that we gave in the introduction.

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

Imagine we want to evaluate these at a particular temperature, in a system of 100 spins.

==== '''<big>TASK 3a</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 1×109 configurations per second with our computer. How long will it take to evaluate a single value of &lang;''M''&rang;''T''?'''

Clearly, we need to try a cleverer approach.

For the overwhelming majority of the possible configurations, the Boltzmann factor, <math>\exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\}</math>, will be very small, and that state will contribute very little to the average value. We can save an enormous amount of time, and make this problem tractable, if we only consider the states whose Boltzmann factor is not so vanishingly small. This technique is called ''importance sampling'' — instead of sampling every point in phase space, we sample only those that the system is likely to occupy.

==Importance Sampling==

This is easily stated, of course, but how can we know which states contribute a lot to the average without actually calculating the energy? So far, we have been imagining generating random states of the system; that is to say, we have been choosing ''α'' from a uniform distribution in which no arrangement of spins is any more likely to be chosen than any other. If we instead had a method which generated states randomly from the probability distribution <math>\exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\}</math>, then our problem would be solved. The method which allows us to do this was developed in the 1950s, and is called the [https://doi.org/10.1063/1.1699114 Metropolis Monte Carlo (MMC)] method, or more often simply "Monte Carlo simulation". This was one of the major breakthroughs of the early days of computational science, and in the year 2000, the IEEE named among the [https://doi.org/10.1109/MCISE.2000.814652 "Top 10 Algorithms of the 20th Century"].

The algorithm is as follows:

# Start from a given configuration of spins, ''α0'', with energy ''E0''.
# Choose a single spin '''at random''', and "flip" it, to generate a new configuration ''α1''.
# Calculate the energy of this new configuration, ''E1''.
# Calculate the energy difference between the states, ''ΔE = E1 - E0''.
## If the ''ΔE ≤ 0'' (the spin flipping decreased the energy), then we '''accept''' the new configuration.
##* We set ''α0 = α1'', and ''E0 = E1'', and then '''go to step 5'''.
## If ''ΔE > 0'', the spin flipping increased the energy. By considering the probability of observing the starting and final states, ''α0'' and ''α1'', it can be shown that the probability for the transition between the two to occur is <math>\exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>. To ensure that we only accept this kind of spin flip with the correct probability, we use the following procedure:
### Choose a random number, ''R'', in the interval [0,1).
### If <math>R \leq \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''accept''' the new configuration.
###* We set ''α0 = α1'', and ''E0 = E1'', and then '''go to step 5'''.
### If <math>R > \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''reject''' the new configuration.
###* ''α0'' and ''E0'' are left unchanged. '''Go to step 5'''.
# Update the running averages of the energy and magnetisation.
# Monte Carlo "step" complete, '''return to step 2'''.

Step 4 is the key to the method. By accepting moves which increase the energy only with a certain probability, we ensure that the sequence of states that we generate is correctly distributed. A transition which carries a very large energy penalty, ''ΔE'' is extremely unlikely to be selected. The use of random numbers in this step is the reason that the method acquired the name "Monte Carlo", after the casinos located there!

Note, unlike molecular dynamics, each step only moves a single atom. In Monte Carlo simulation, there is the notion of a cycle which is ''N''spins-steps, so that each spin on average has an attempted flip per cycle.

If you are interested in the mathematical details of why this procedure generates a sequence of states distributed in the correct way, consult the Monte Carlo chapter of [https://library-search.imperial.ac.uk/permalink/f/o1297h/TN_cdi_askewsholts_vlebooks_9780080519982 "Understanding Molecular Simulation"], by Frenkel and Smit, but a discussion of this is not required for the experiment. Chapter 7 of [https://library-search.imperial.ac.uk/permalink/44IMP_INST/mek6kh/alma997227534401591 "Statistical Mechanics: Theory and Molecular Simulation"] by Tuckerman also gives a good description.

==Modifying IsingLattice.py==

Our IsingLattice object contains a stub function, '''<code>montecarlostep(T)</code>''', which takes a single argument — the temperature at which we want to do the cycle. There are also attributes used to calculate the average energy and magnetisation and respective variances: <code>self.E_tally</code>, <code>self.E2_tally</code>, <code>self.M_tally</code>, <code>self.M2_tally</code>, and <code>self.n_steps</code>. You should use these values to keep a running sum of all the energies, squared energies, magnetisations, and squared magnetisations that your system has sampled, as well as the number of steps that have been performed. Finally, there is also a function <code>statistics()</code>, that takes no arguments and should return the average energy, squared energy, magnetisation, squared magnetisation and current step, '''in that order'''.

==== '''<big>TASK 3b</big>:''' ====
'''Implement a single step of the above algorithm in the <code>montecarlostep(T)</code> function. This function should return the energy of your lattice and the magnetisation at the end of the cycle. Complete the <code>statistics()</code> function. This should return the following quantities whenever it is called: &lang;''E''&rang;, &lang;''E''2&rang;, &lang;''M''&rang;, &lang;''M''2&rang;, and ''N''steps, the number of Monte Carlo steps that have elapsed. Put your code for these functions in the report.'''

''Note on units:'' In our simulation we will be using an unitless scaled energy <math>\tilde{E}=\frac{E}{J}</math>, and a unitless scaled temperature <math>\tilde{T}=\frac{k_B T}{J}</math>, where ''J'' is the interaction energy between spins discussed previously. The use of these scaled quantities in practice means that in the code we can set ''kB'' to numerically equal to 1. Going forward we will drop the &tilde; from the notation of energy and temperature.

You have been provided with a script, '''ILanim.py''', which will use your <code>IsingLattice</code> object to run a Monte Carlo simulation and display the output on the screen in real time. By default, this simulation will run an 8×8 lattice at ''T''=1.0, which is below the critical temperature. You will see a representation of the lattice, and graphs of the energy and magnetisation '''per spin'''. When you close the window, the script will use your '''<code>statistics()</code>''' function to print the average energy and magnetisation from the simulation. Note, you should run this script from the command line.

==== '''<big>TASK 3c</big>:''' ====
'''If ''T'' < ''TC'', do you expect a spontaneous magnetisation (i.e. do you expect &lang;''M''&rang; ≠ 0)? When the state of the simulation appears to stop changing (when you have reached an equilibrium state), use the controls to export the output to an SVG or PNG file and attach this to your report. You should also include the output from your statistics() function.'''

'''<big>This is the third section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]]</big>'''

Programming a 2D Ising Model/Accelerating the code

2026-05-03T13:50:27Z

Gmallia:

'''<big>This is the fourth section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of temperature|5. The effect of temperature]]</big>'''

= 4. Accelerating the code =
So far, your energy() and magnetisation() functions probably use Python loops to go over every spin in your system (and a double loop for the energy). One of the drawbacks of Python's flexibility is a performance penalty — though developing code in Python tends to be
faster than other compiled programming languages, the resulting program is much, much slower.

Fortunately, the NumPy module provides an interface to a number of functions written to give high performance.
These are written in the C programming language under the hood and take numpy arrays, which are contiguous arrays of a single data type,
unlike Python lists, which may contain mixes of data types.

In this section, we are going to use these functions to improve the performance of our simulation.

When trying to optimise the performance of something, you should always measure first, rather than making assumptions. Although there are some very obvious candidates in our program, you should always measure properly.

The <code>ILtimetrial.py</code> script measures how fast it takes for 2000 iterations of your code to run.

The <code>ILProfiling.py</code> script is the same, except uses the python standard library package [https://docs.python.org/3/library/profile.html cProfile] to sample what functions are being called and give you information
about the time spent in each. Analysing the output of the script, you will see that by far the majority of the time is spent in the energy method, followed by the magnetisation one. This makes sense, as each Monte-Carlo step can call the energy function multiple times.

(If you want to profile your code in more detail, please seek advice from a demonstrator.)

In order to improve the performance of this, we need to:
* Run the energy call fewer times per step.
* Improve the performance of the energy function.

Fortunately, the NumPy module provides an interface to a number of functions written to give high performance. In this section, we are going to use these functions to improve the performance of our simulation.

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

==== '''<big>TASK 4c</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! What makes this new version faster?'''

==== '''<big>Task 4d:</big>''' ====
'''<big>Making the code even faster </big>: Consider trying to make the Monte Carlo step function substantially faster by calculating the difference in energy directly, in a way that avoids looping over the whole lattice. Think back how you answered some of the introduction questions. You can make it so the energy and magnetisation only needs to be calculated explicitly at the beginning of the simulation.'''

'''<big>This is the fourth section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of temperature|5. The effect of temperature]]</big>'''

Programming a 2D Ising Model/The effect of temperature

2026-05-03T13:49:53Z

Gmallia:

'''<big>This is the fifth section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of system size|6. The effect of system size]]</big>'''

= 5. The effect of temperature =
We now have enough code to start to investigate the phase behaviour of our model. In the first instance, we are going to make plots of the average energy and magnetisation of the system as a function of temperature, to get an idea of where the change between energetic and entropic domination takes place, but we first have to make a small change to the code that averages the energy and magnetisation.

==Correcting the Averaging Code==
As you saw in the previous part of the experiment, when we begin from a random configuration, it takes a few Monte Carlo cycles to reach the equilibrium state of the system at a particular temperature. This manifests itself as a relatively sharp drop in the energy of the system in the first few hundred or thousand steps, before the energy starts to fluctuate around a constant value. If we want our average properties to be correct, we should only average in the region where the ''average'' energy and magnetisation are constant.

==== '''<big>TASK 5a</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 a few different temperatures and lattice sizes. Approximately how many cycles (remember, a cycle is ''N''spins steps) are needed for the system to go from its random starting position to the equilibrium state? Modify your <code>__init__()</code>, <code>statistics()</code> and <code>montecarlostep()</code> methods so that the first ''N'' cycles of the simulation are ignored when calculating the averages. Be careful of off-by-one errors. 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 number.'''

==Running over a range of temperatures==

The script ILtemperaturerange.py performs simulations of a given length over a range of temperatures. The temperature range is determined by the parameters given to the <code>np.arange</code> function near the top of the file. Starting from a uniformly aligned lattice, a simulation of length <code>n_steps</code> will be performed at temperature <code>T_cold</code>, and the average energy and magnetisation will be recorded. A new simulation of length <code>n_steps</code> will then be performed at temperature <code>T_cold + dT</code>, and so on until all temperatures in the range have been simulated. This results in a NumPy array of five columns: the temperature, then the four statistical quantities: <code>ave_energies</code>, <code>ave_sqenergies</code> (squared energies), <code>ave_magnetisations</code>, <code>ave_sqmagnetisations</code> (squared magnetisations).

==== '''<big>TASK 5b</big>:''' ====
'''Use ILtemperaturerange.py to plot the average energy and magnetisation for each temperature, ''with error bars'', for an 8 × 8 lattice. Use multiple repeats to calculate these error bars. Use your intuition and results from the script ILfinalframe.py to estimate how many cycles each simulation should be. The temperature range 0.25 to 5.0 is sufficient. Use as many temperature points as you feel necessary to illustrate the behaviour, but do not use a temperature spacing larger than 0.5. The NumPy function <code>savetxt()</code> 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.'''

'''<big>This is the fifth section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of system size|6. The effect of system size]]</big>'''

Programming a 2D Ising Model/Calculating the energy and magnetisation

2026-05-03T13:48:51Z

Gmallia: /* 2. Calculating the energy and magnetisation */

'''<big>This is the second section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

= 2. Calculating the energy and magnetisation =

==Getting the files for the experiment==

'''Various scripts have been prepared to assist you with this experiment, and you can obtain them by downloading the files from the Blackboard.'''

We recommend you open the extracted folder with JupyterLab, which has an integrated terminal where you can run the simulation scripts.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. This class will be the basis for our model of the system, and is restricted to modelling the two dimensional case. It contains a number of stub methods that you will complete as you work through the experiment.

== A note about Python Classes ==
The file IsingLattice.py makes use of a Python feature called ''class''. You should be familiar with the idea that all Python ''objects'' have an associated ''type'' -- you have already encountered objects such as integers, string, lists, and NumPy arrays. Classes allow us to define new types of object. In this case, we define a type called <code>IsingLattice</code>:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. We can see, for example, that each <code>IsingLattice</code> object has the methods <code>energy()</code>, <code>magnetisation()</code>, <code>montecarlostep(T)</code>, and <code>statistics()</code> (the method <code>__init__</code> will be explained shortly). Each of these methods takes a special first argument called self, which is a variable pointing to the current copy of IsingLattice. To use the new IsingLattice object that we have defined, we use code like:

<pre>
import IsingLattice

il = IsingLattice(5, 5) #create an IsingLattice object with 5 row and 5 columns

energy = il.energy() #call the energy() method of our IsingLattice object
</pre>

When writing code "within" the class definition, all attributes and methods belonging to the class must be accessed using the "self." notation. You can see an example of this in the <code>montecarlostep()</code> method:

<pre>
energy = self.energy()
</pre>

This line creates a variable called energy which contains whatever value was returned by the method <code>energy()</code> which is part of the <code>IsingLattice</code> class. If you just want to define a local variable instead, like an iterator in a loop, the <code>self.</code> notation is not required.

<pre>
def magnetisation(self):
# we have to use self.lattice, because the lattice variable is part of the object, not local to this function
for i in self.lattice:
...
# we use i, not self.i
# we don't need to access i outside this function
</pre>

Have a look at the other files that you have been given to see how this IsingLattice object is used.

== The Constructor ==
The only method currently filled in is <code>__init__</code>. This is a special method called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this method is called automatically. It takes two arguments (excluding "self"), which are simply the number of rows and columns that our lattice will have. The attribute <code>self.lattice</code> is then created, which is a NumPy array with the requested number of rows and columns. The function <code>np.random.choice</code> is used to create a random collection of spins for our initial configuration. You can see the documentation for the function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html numpy.random.choice].

==== '''<big>TASK 2a</big>:''' ====
'''Complete the functions <code>energy()</code> and <code>magnetisation()</code>, which should return the energy of the lattice and the total magnetisation, respectively. In the <code>energy()</code> function you may assume that ''J''=1.0 at all times (in fact, we are working in ''reduced units'' in which all energies are a multiple of ''J'', but there will be more information about this in later sections). Do not worry about the efficiency of the code at the moment — we will address the speed in a later part of the experiment. Paste your code in your report and describe how it works.'''

==Testing the files==

When you have completed the energy() and magnetisation() methods, we need to test that they work correctly.
We have provided you with a script, '''ILcheck.py''', which will create three IsingLattice objects and check their energies and magnetisations, displaying these on a graph. One configuration corresponds to the energy minimum, one to the energy maximum, and one to an intermediate state.

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

There are also some unit tests available in the <code>test_energy.py</code> file. These tests are designed to be run by simply running the [https://docs.pytest.org pytest] command in your terminal,
<pre>pytest</pre>
You should make sure your code passes the tests for the <code>energy()</code> and <code>magnetisation()</code> methods at this time. The file contains tests for methods that will only be implemented later on, so some failures are expected for now.
These sorts of test are useful because they are easier to write and faster to run and can be made to run automatically.

You can also debug and test your program on a Jupyter notebook or IPython console by importing your IsingLattice class, creating one object, and inspecting its properties. If you do this, we recommend that you start your session running the commands:

<pre>
%load_ext autoreload
%autoreload 2
</pre>

This will ensure that when you try executing your code, you always run the most recent version (please note this isn't fool proof and if you find problems you may need to restart the kernel)!

'''<big>This is the second section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

Programming a 2D Ising Model/Calculating the energy and magnetisation

2026-05-03T13:46:06Z

Gmallia: /* 2. Calculating the energy and magnetisation */

'''<big>This is the second section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

= 2. Calculating the energy and magnetisation =

==Getting the files for the experiment==

'''Various scripts have been prepared to assist you with this experiment, and you can obtain them by downloading the files from the Blackboard.'''

We recommend you open the extracted folder with JupyterLab, which has an integrated terminal where you can run the simulation scripts.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. This class will be the basis for our model of the system, and is restricted to modelling the two dimensional case. It contains a number of stub methods that you will complete as you work through the experiment.

===A note about Python Classes===

The file IsingLattice.py makes use of a Python feature called ''class''. You should be familiar with the idea that all Python ''objects'' have an associated ''type'' -- you have already encountered objects such as integers, string, lists, and NumPy arrays. Classes allow us to define new types of object. In this case, we define a type called <code>IsingLattice</code>:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. We can see, for example, that each <code>IsingLattice</code> object has the methods <code>energy()</code>, <code>magnetisation()</code>, <code>montecarlostep(T)</code>, and <code>statistics()</code> (the method <code>__init__</code> will be explained shortly). Each of these methods takes a special first argument called self, which is a variable pointing to the current copy of IsingLattice. To use the new IsingLattice object that we have defined, we use code like:

<pre>
import IsingLattice

il = IsingLattice(5, 5) #create an IsingLattice object with 5 row and 5 columns

energy = il.energy() #call the energy() method of our IsingLattice object
</pre>

When writing code "within" the class definition, all attributes and methods belonging to the class must be accessed using the "self." notation. You can see an example of this in the <code>montecarlostep()</code> method:

<pre>
energy = self.energy()
</pre>

This line creates a variable called energy which contains whatever value was returned by the method <code>energy()</code> which is part of the <code>IsingLattice</code> class. If you just want to define a local variable instead, like an iterator in a loop, the <code>self.</code> notation is not required.

<pre>
def magnetisation(self):
# we have to use self.lattice, because the lattice variable is part of the object, not local to this function
for i in self.lattice:
...
# we use i, not self.i
# we don't need to access i outside this function
</pre>

Have a look at the other files that you have been given to see how this IsingLattice object is used.

== The Constructor ==
The only method currently filled in is <code>__init__</code>. This is a special method called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this method is called automatically. It takes two arguments (excluding "self"), which are simply the number of rows and columns that our lattice will have. The attribute <code>self.lattice</code> is then created, which is a NumPy array with the requested number of rows and columns. The function <code>np.random.choice</code> is used to create a random collection of spins for our initial configuration. You can see the documentation for the function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html numpy.random.choice].

==== '''<big>TASK 2a</big>:''' ====
'''Complete the functions <code>energy()</code> and <code>magnetisation()</code>, which should return the energy of the lattice and the total magnetisation, respectively. In the <code>energy()</code> function you may assume that ''J''=1.0 at all times (in fact, we are working in ''reduced units'' in which all energies are a multiple of ''J'', but there will be more information about this in later sections). Do not worry about the efficiency of the code at the moment — we will address the speed in a later part of the experiment. Paste your code in your report and describe how it works.'''

==Testing the files==

When you have completed the energy() and magnetisation() methods, we need to test that they work correctly.
We have provided you with a script, '''ILcheck.py''', which will create three IsingLattice objects and check their energies and magnetisations, displaying these on a graph. One configuration corresponds to the energy minimum, one to the energy maximum, and one to an intermediate state.

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

There are also some unit tests available in the <code>test_energy.py</code> file. These tests are designed to be run by simply running the [https://docs.pytest.org pytest] command in your terminal,
<pre>pytest</pre>
You should make sure your code passes the tests for the <code>energy()</code> and <code>magnetisation()</code> methods at this time. The file contains tests for methods that will only be implemented later on, so some failures are expected for now.
These sorts of test are useful because they are easier to write and faster to run and can be made to run automatically.

You can also debug and test your program on a Jupyter notebook or IPython console by importing your IsingLattice class, creating one object, and inspecting its properties. If you do this, we recommend that you start your session running the commands:

<pre>
%load_ext autoreload
%autoreload 2
</pre>

This will ensure that when you try executing your code, you always run the most recent version (please note this isn't fool proof and if you find problems you may need to restart the kernel)!

'''<big>This is the second section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

Programming a 2D Ising Model/Calculating the energy and magnetisation

2026-05-03T13:44:56Z

Gmallia: /* Modifying the files */

'''<big>This is the second section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

= 2. Calculating the energy and magnetisation =

==Getting the files for the experiment==

'''Various scripts have been prepared to assist you with this experiment, and you can obtain them by downloading the files from the Blackboard.'''

We recommend you open the extracted folder with JupyterLab, which has an integrated terminal where you can run the simulation scripts.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. This class will be the basis for our model of the system, and is restricted to modelling the two dimensional case. It contains a number of stub methods that you will complete as you work through the experiment.

===A note about Python Classes===

The file IsingLattice.py makes use of a Python feature called ''class''. You should be familiar with the idea that all Python ''objects'' have an associated ''type'' -- you have already encountered objects such as integers, string, lists, and NumPy arrays. Classes allow us to define new types of object. In this case, we define a type called <code>IsingLattice</code>:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. We can see, for example, that each <code>IsingLattice</code> object has the methods <code>energy()</code>, <code>magnetisation()</code>, <code>montecarlostep(T)</code>, and <code>statistics()</code> (the method <code>__init__</code> will be explained shortly). Each of these methods takes a special first argument called self, which is a variable pointing to the current copy of IsingLattice. To use the new IsingLattice object that we have defined, we use code like:

<pre>
import IsingLattice

il = IsingLattice(5, 5) #create an IsingLattice object with 5 row and 5 columns

energy = il.energy() #call the energy() method of our IsingLattice object
</pre>

When writing code "within" the class definition, all attributes and methods belonging to the class must be accessed using the "self." notation. You can see an example of this in the <code>montecarlostep()</code> method:

<pre>
energy = self.energy()
</pre>

This line creates a variable called energy which contains whatever value was returned by the method <code>energy()</code> which is part of the <code>IsingLattice</code> class. If you just want to define a local variable instead, like an iterator in a loop, the <code>self.</code> notation is not required.

<pre>
def magnetisation(self):
# we have to use self.lattice, because the lattice variable is part of the object, not local to this function
for i in self.lattice:
...
# we use i, not self.i
# we don't need to access i outside this function
</pre>

Have a look at the other files that you have been given to see how this IsingLattice object is used.

===The Constructor===

The only method currently filled in is <code>__init__</code>. This is a special method called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this method is called automatically. It takes two arguments (excluding "self"), which are simply the number of rows and columns that our lattice will have. The attribute <code>self.lattice</code> is then created, which is a NumPy array with the requested number of rows and columns. The function <code>np.random.choice</code> is used to create a random collection of spins for our initial configuration. You can see the documentation for the function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html numpy.random.choice].

==== '''<big>TASK 2a</big>:''' ====
'''Complete the functions <code>energy()</code> and <code>magnetisation()</code>, which should return the energy of the lattice and the total magnetisation, respectively. In the <code>energy()</code> function you may assume that ''J''=1.0 at all times (in fact, we are working in ''reduced units'' in which all energies are a multiple of ''J'', but there will be more information about this in later sections). Do not worry about the efficiency of the code at the moment — we will address the speed in a later part of the experiment. Paste your code in your report and describe how it works.'''

==Testing the files==

When you have completed the energy() and magnetisation() methods, we need to test that they work correctly.
We have provided you with a script, '''ILcheck.py''', which will create three IsingLattice objects and check their energies and magnetisations, displaying these on a graph. One configuration corresponds to the energy minimum, one to the energy maximum, and one to an intermediate state.

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

There are also some unit tests available in the <code>test_energy.py</code> file. These tests are designed to be run by simply running the [https://docs.pytest.org pytest] command in your terminal,
<pre>pytest</pre>
You should make sure your code passes the tests for the <code>energy()</code> and <code>magnetisation()</code> methods at this time. The file contains tests for methods that will only be implemented later on, so some failures are expected for now.
These sorts of test are useful because they are easier to write and faster to run and can be made to run automatically.

You can also debug and test your program on a Jupyter notebook or IPython console by importing your IsingLattice class, creating one object, and inspecting its properties. If you do this, we recommend that you start your session running the commands:

<pre>
%load_ext autoreload
%autoreload 2
</pre>

This will ensure that when you try executing your code, you always run the most recent version (please note this isn't fool proof and if you find problems you may need to restart the kernel)!

'''<big>This is the second section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

Programming a 2D Ising Model/Introduction to the Ising model

2026-05-03T13:43:54Z

Gmallia: /* The Ising model */

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]]</big>'''

= 1. Introduction to the Ising model =

== The Ising model ==
The Ising model is a simple physical model proposed in the 1920s to understand the behaviour of ferromagnets. In most materials, the magnetic dipoles of the atoms have no preferred orientation, and the bulk material has no net magnetic moment. In certain materials, however, there is preferential alignment of the magnetic dipoles along certain directions. These materials, such as iron, exhibit an overall magnetic moment, and are called ferromagnetic.

The phenomenon is based on competition between two physical effects — energy minimisation and entropy maximisation. On the one hand, parallel alignment of magnetic moments is energetically favourable. On the other hand, aligning the system in this way reduces the entropy.

In the Ising model, we consider a collection of spins, ''si'', each of which can either have the value "up" (+1) or "down" (-1), arranged on a lattice. In one dimension, this takes the form of a "chain" of adjacent spins, ''s1'' to ''sN''. In two dimensions we have a grid of cells, each of which contains one spin. In three dimensions, the system is a cuboid, and so on. The arrangements for 1, 2 and 3 dimensions are illustrated in '''figure 1'''.

[[File:ThirdYearCMPExpt-IsingSketch.png|300px|thumb|right|'''Figure 1''': Illustration of an Ising lattice in one (N=5), two (N=5×5), and three (N=5×5×5) dimensions. Red cells indicate the "up" spin", and blue cells indicate the "down" spin.]]

The rules of the model are simple:

* When no magnetic field is applied, the only energetic contribution comes from the interaction between spins in ''adjacent'' lattice cells
* The total interaction energy is defined as <math>- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>, where ''J'' is a constant which controls the strength of interaction and the notation ''j'' &in; neighbours(''i'') indicates that spin ''j'' must lie in a lattice cell adjacent to spin ''i''. The factor of one half is included to make sure that we do not count each spin-spin interaction twice.
* Periodic boundary conditions are applied. In other words, the "final" spin in a particular direction is considered to be adjacent to the "first" spin in that direction. In '''figure 2''', for example, cell 1 has neighbour interactions with cells 2, 3, 4, and 7, while cell 9 interacts with cells 3, 6, 7, and 8.

[[File:ThirdYearCMPExpt-IsingNumberedCell.png|300px|thumb|right|'''Figure 2''': Illustration of a 2D Ising lattice, with numbered cells to indicate the neighbour interactions.]]

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

==Phase Transitions==

The Ising model is interesting because it is one of the simplest physical models to display a phase transition (so long as we are in two dimensions or more). The temperature of the system controls the balance between the energetic driving force, which wants the system to adopt the lowest energy configuration that you imagined in the previous section, and the entropic driving force, which seeks to maximise the number of configurations available to the system.

==== '''<big>TASK 1b</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 (''D''=3, ''N''=1000)? How much entropy does the system gain by doing so?'''

At low temperatures, little thermal energy is available to overcome such energy barriers, and we expect the system to adopt the lowest energy (the energetic driving force will dominate). Under these conditions, we expect a large number of parallel spins, and we can characterise this by considering the total ''magnetisation'' of the system: <math>M = \sum_i s_i</math>. At higher temperatures, however, enough thermal energy will be available that spins can flip relatively freely. Under these conditions, the entropic driving force should dominate. The transition temperature between the two regimes is called the ''Curie temperature'', ''TC''. Below ''TC'', we observe the highly ordered magnetised phase, while above ''TC'' a disordered, demagnetised, phase emerges.

==== '''<big>TASK 1c</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 ''D''=3, ''N''=1000 at absolute zero?'''

We are interested, therefore, in determining the magnetisation ''M'', and energy ''E'', of the system as functions of temperature ''T''. We can write statistical mechanical equations for these quantities:

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

where ''α'' is used as a shorthand to represent all spins in the system, and ''Z'' is the partition function,
<math>Z = \sum_\alpha \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>.

It was discovered relatively quickly after the model was first proposed that the partition function in the one dimensional case can be solved analytically (though doing so is beyond the scope of this experiment):

<math>Z\left(T, N\right) = \left[2\cosh\left(\frac{J}{k_BT}\right)\right]^N</math>

Much later, [http://en.wikipedia.org/wiki/Lars_Onsager Lars Onsager] established that the two dimensional partition function can also be found analytically, though the mathematics is significantly more involved and the resulting expression is correspondingly much more complicated. In dimensions greater than two, no analytical solution is known.

Luckily for us, the problem can instead be tackled by numerical methods, and it is this that will be the focus of the experiment. We will discuss how the equations for ''M'' and ''E'' can be solved, but first we must set up a Python script to model the lattice.

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]]</big>'''

Programming a 2D Ising Model/Introduction to the Ising model

2026-05-03T13:42:43Z

Gmallia: /* The Ising model */

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]]</big>'''

= 1. Introduction to the Ising model =

== The Ising model ==
The Ising model is a simple physical model proposed in the 1920s to understand the behaviour of ferromagnets. In most materials, the magnetic dipoles of the atoms have no preferred orientation, and the bulk material has no net magnetic moment. In certain materials, however, there is preferential alignment of the magnetic dipoles along certain directions. These materials, such as iron, exhibit an overall magnetic moment, and are called ferromagnetic.

The phenomenon is based on competition between two physical effects — energy minimisation and entropy maximisation. On the one hand, parallel alignment of magnetic moments is energetically favourable. On the other hand, aligning the system in this way reduces the entropy.

In the Ising model, we consider a collection of spins, ''si'', each of which can either have the value "up" (+1) or "down" (-1), arranged on a lattice. In one dimension, this takes the form of a "chain" of adjacent spins, ''s1'' to ''sN''. In two dimensions we have a grid of cells, each of which contains one spin. In three dimensions, the system is a cuboid, and so on. The arrangements for 1, 2 and 3 dimensions are illustrated in '''figure 1'''.

[[File:ThirdYearCMPExpt-IsingSketch.png|300px|thumb|right|'''Figure 1''': Illustration of an Ising lattice in one (N=5), two (N=5×5), and three (N=5×5×5) dimensions. Red cells indicate the "up" spin", and blue cells indicate the "down" spin.]]

The rules of the model are simple:

* When no magnetic field is applied, the only energetic contribution comes from the interaction between spins in ''adjacent'' lattice cells
* The total interaction energy is defined as <math>- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>, where ''J'' is a constant which controls the strength of interaction and the notation ''j'' &in; neighbours(''i'') indicates that spin ''j'' must lie in a lattice cell adjacent to spin ''i''. The factor of one half is included to make sure that we do not count each spin-spin interaction twice.
* Periodic boundary conditions are applied. In other words, the "final" spin in a particular direction is considered to be adjacent to the "first" spin in that direction. In '''figure 2''', for example, cell 1 has neighbour interactions with cells 2, 3, 4, and 7, while cell 9 interacts with cells 3, 6, 7, and 8.

[[File:ThirdYearCMPExpt-IsingNumberedCell.png|300px|thumb|right|'''Figure 2''': Illustration of a 2D Ising lattice, with numbered cells to indicate the neighbour interactions.]]

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

==Phase Transitions==

The Ising model is interesting because it is one of the simplest physical models to display a phase transition (so long as we are in two dimensions or more). The temperature of the system controls the balance between the energetic driving force, which wants the system to adopt the lowest energy configuration that you imagined in the previous section, and the entropic driving force, which seeks to maximise the number of configurations available to the system.

'''<big>TASK 1b</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 (''D''=3, ''N''=1000)? How much entropy does the system gain by doing so?'''

At low temperatures, little thermal energy is available to overcome such energy barriers, and we expect the system to adopt the lowest energy (the energetic driving force will dominate). Under these conditions, we expect a large number of parallel spins, and we can characterise this by considering the total ''magnetisation'' of the system: <math>M = \sum_i s_i</math>. At higher temperatures, however, enough thermal energy will be available that spins can flip relatively freely. Under these conditions, the entropic driving force should dominate. The transition temperature between the two regimes is called the ''Curie temperature'', ''TC''. Below ''TC'', we observe the highly ordered magnetised phase, while above ''TC'' a disordered, demagnetised, phase emerges.

'''<big>TASK 1c</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 ''D''=3, ''N''=1000 at absolute zero?'''

We are interested, therefore, in determining the magnetisation ''M'', and energy ''E'', of the system as functions of temperature ''T''. We can write statistical mechanical equations for these quantities:

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

where ''α'' is used as a shorthand to represent all spins in the system, and ''Z'' is the partition function,
<math>Z = \sum_\alpha \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>.

It was discovered relatively quickly after the model was first proposed that the partition function in the one dimensional case can be solved analytically (though doing so is beyond the scope of this experiment):

<math>Z\left(T, N\right) = \left[2\cosh\left(\frac{J}{k_BT}\right)\right]^N</math>

Much later, [http://en.wikipedia.org/wiki/Lars_Onsager Lars Onsager] established that the two dimensional partition function can also be found analytically, though the mathematics is significantly more involved and the resulting expression is correspondingly much more complicated. In dimensions greater than two, no analytical solution is known.

Luckily for us, the problem can instead be tackled by numerical methods, and it is this that will be the focus of the experiment. We will discuss how the equations for ''M'' and ''E'' can be solved, but first we must set up a Python script to model the lattice.

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]]</big>'''

Programming a 2D Ising Model/The effect of temperature

2026-05-03T13:41:04Z

Gmallia:

'''<big>This is the fifth section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of system size|6. The effect of system size]]</big>'''

= 5. The effect of temperature =
We now have enough code to start to investigate the phase behaviour of our model. In the first instance, we are going to make plots of the average energy and magnetisation of the system as a function of temperature, to get an idea of where the change between energetic and entropic domination takes place, but we first have to make a small change to the code that averages the energy and magnetisation.

==Correcting the Averaging Code==
As you saw in the previous part of the experiment, when we begin from a random configuration, it takes a few Monte Carlo cycles to reach the equilibrium state of the system at a particular temperature. This manifests itself as a relatively sharp drop in the energy of the system in the first few hundred or thousand steps, before the energy starts to fluctuate around a constant value. If we want our average properties to be correct, we should only average in the region where the ''average'' energy and magnetisation are constant.

'''<big>TASK 5a</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 a few different temperatures and lattice sizes. Approximately how many cycles (remember, a cycle is ''N''spins steps) are needed for the system to go from its random starting position to the equilibrium state? Modify your <code>__init__()</code>, <code>statistics()</code> and <code>montecarlostep()</code> methods so that the first ''N'' cycles of the simulation are ignored when calculating the averages. Be careful of off-by-one errors. 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 number.'''
'''

==Running over a range of temperatures==

The script ILtemperaturerange.py performs simulations of a given length over a range of temperatures. The temperature range is determined by the parameters given to the <code>np.arange</code> function near the top of the file. Starting from a uniformly aligned lattice, a simulation of length <code>n_steps</code> will be performed at temperature <code>T_cold</code>, and the average energy and magnetisation will be recorded. A new simulation of length <code>n_steps</code> will then be performed at temperature <code>T_cold + dT</code>, and so on until all temperatures in the range have been simulated. This results in a NumPy array of five columns: the temperature, then the four statistical quantities: <code>ave_energies</code>, <code>ave_sqenergies</code> (squared energies), <code>ave_magnetisations</code>, <code>ave_sqmagnetisations</code> (squared magnetisations).

'''<big>TASK 5b</big>: Use ILtemperaturerange.py to plot the average energy and magnetisation for each temperature, ''with error bars'', for an 8 × 8 lattice. Use multiple repeats to calculate these error bars. Use your intuition and results from the script ILfinalframe.py to estimate how many cycles each simulation should be. The temperature range 0.25 to 5.0 is sufficient. Use as many temperature points as you feel necessary to illustrate the behaviour, but do not use a temperature spacing larger than 0.5. The NumPy function <code>savetxt()</code> 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.'''

'''<big>This is the fifth section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of system size|6. The effect of system size]]</big>'''

Programming a 2D Ising Model/The effect of temperature

2026-05-03T13:40:11Z

Gmallia: /* Running over a range of temperatures */

'''<big>This is the fifth section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of system size|6. The effect of system size]]</big>'''

We now have enough code to start to investigate the phase behaviour of our model. In the first instance, we are going to make plots of the average energy and magnetisation of the system as a function of temperature, to get an idea of where the change between energetic and entropic domination takes place, but we first have to make a small change to the code that averages the energy and magnetisation.

==Correcting the Averaging Code==
As you saw in the previous part of the experiment, when we begin from a random configuration, it takes a few Monte Carlo cycles to reach the equilibrium state of the system at a particular temperature. This manifests itself as a relatively sharp drop in the energy of the system in the first few hundred or thousand steps, before the energy starts to fluctuate around a constant value. If we want our average properties to be correct, we should only average in the region where the ''average'' energy and magnetisation are constant.

'''<big>TASK 5a</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 a few different temperatures and lattice sizes. Approximately how many cycles (remember, a cycle is ''N''spins steps) are needed for the system to go from its random starting position to the equilibrium state? Modify your <code>__init__()</code>, <code>statistics()</code> and <code>montecarlostep()</code> methods so that the first ''N'' cycles of the simulation are ignored when calculating the averages. Be careful of off-by-one errors. 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 number.'''
'''

==Running over a range of temperatures==

The script ILtemperaturerange.py performs simulations of a given length over a range of temperatures. The temperature range is determined by the parameters given to the <code>np.arange</code> function near the top of the file. Starting from a uniformly aligned lattice, a simulation of length <code>n_steps</code> will be performed at temperature <code>T_cold</code>, and the average energy and magnetisation will be recorded. A new simulation of length <code>n_steps</code> will then be performed at temperature <code>T_cold + dT</code>, and so on until all temperatures in the range have been simulated. This results in a NumPy array of five columns: the temperature, then the four statistical quantities: <code>ave_energies</code>, <code>ave_sqenergies</code> (squared energies), <code>ave_magnetisations</code>, <code>ave_sqmagnetisations</code> (squared magnetisations).

'''<big>TASK 5b</big>: Use ILtemperaturerange.py to plot the average energy and magnetisation for each temperature, ''with error bars'', for an 8 × 8 lattice. Use multiple repeats to calculate these error bars. Use your intuition and results from the script ILfinalframe.py to estimate how many cycles each simulation should be. The temperature range 0.25 to 5.0 is sufficient. Use as many temperature points as you feel necessary to illustrate the behaviour, but do not use a temperature spacing larger than 0.5. The NumPy function <code>savetxt()</code> 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.'''

'''<big>This is the fifth section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of system size|6. The effect of system size]]</big>'''

Programming a 2D Ising Model/Introduction to Monte Carlo simulation

2026-05-03T13:37:27Z

Gmallia:

'''<big>This is the third section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]]</big>'''

= 3. Introduction to Monte Carlo simulation =

==Average Energy and Magnetisation==

Consider again the expressions for the average energy and magnetisation that we gave in the introduction.

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

Imagine we want to evaluate these at a particular temperature, in a system of 100 spins.

'''<big>TASK 3a</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 1×109 configurations per second with our computer. How long will it take to evaluate a single value of &lang;''M''&rang;''T''?'''

Clearly, we need to try a cleverer approach.

For the overwhelming majority of the possible configurations, the Boltzmann factor, <math>\exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\}</math>, will be very small, and that state will contribute very little to the average value. We can save an enormous amount of time, and make this problem tractable, if we only consider the states whose Boltzmann factor is not so vanishingly small. This technique is called ''importance sampling'' — instead of sampling every point in phase space, we sample only those that the system is likely to occupy.

==Importance Sampling==

This is easily stated, of course, but how can we know which states contribute a lot to the average without actually calculating the energy? So far, we have been imagining generating random states of the system; that is to say, we have been choosing ''α'' from a uniform distribution in which no arrangement of spins is any more likely to be chosen than any other. If we instead had a method which generated states randomly from the probability distribution <math>\exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\}</math>, then our problem would be solved. The method which allows us to do this was developed in the 1950s, and is called the [https://doi.org/10.1063/1.1699114 Metropolis Monte Carlo (MMC)] method, or more often simply "Monte Carlo simulation". This was one of the major breakthroughs of the early days of computational science, and in the year 2000, the IEEE named among the [https://doi.org/10.1109/MCISE.2000.814652 "Top 10 Algorithms of the 20th Century"].

The algorithm is as follows:

# Start from a given configuration of spins, ''α0'', with energy ''E0''.
# Choose a single spin '''at random''', and "flip" it, to generate a new configuration ''α1''.
# Calculate the energy of this new configuration, ''E1''.
# Calculate the energy difference between the states, ''ΔE = E1 - E0''.
## If the ''ΔE ≤ 0'' (the spin flipping decreased the energy), then we '''accept''' the new configuration.
##* We set ''α0 = α1'', and ''E0 = E1'', and then '''go to step 5'''.
## If ''ΔE > 0'', the spin flipping increased the energy. By considering the probability of observing the starting and final states, ''α0'' and ''α1'', it can be shown that the probability for the transition between the two to occur is <math>\exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>. To ensure that we only accept this kind of spin flip with the correct probability, we use the following procedure:
### Choose a random number, ''R'', in the interval [0,1).
### If <math>R \leq \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''accept''' the new configuration.
###* We set ''α0 = α1'', and ''E0 = E1'', and then '''go to step 5'''.
### If <math>R > \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''reject''' the new configuration.
###* ''α0'' and ''E0'' are left unchanged. '''Go to step 5'''.
# Update the running averages of the energy and magnetisation.
# Monte Carlo "step" complete, '''return to step 2'''.

Step 4 is the key to the method. By accepting moves which increase the energy only with a certain probability, we ensure that the sequence of states that we generate is correctly distributed. A transition which carries a very large energy penalty, ''ΔE'' is extremely unlikely to be selected. The use of random numbers in this step is the reason that the method acquired the name "Monte Carlo", after the casinos located there!

Note, unlike molecular dynamics, each step only moves a single atom. In Monte Carlo simulation, there is the notion of a cycle which is ''N''spins-steps, so that each spin on average has an attempted flip per cycle.

If you are interested in the mathematical details of why this procedure generates a sequence of states distributed in the correct way, consult the Monte Carlo chapter of [https://library-search.imperial.ac.uk/permalink/f/o1297h/TN_cdi_askewsholts_vlebooks_9780080519982 "Understanding Molecular Simulation"], by Frenkel and Smit, but a discussion of this is not required for the experiment. Chapter 7 of [https://library-search.imperial.ac.uk/permalink/44IMP_INST/mek6kh/alma997227534401591 "Statistical Mechanics: Theory and Molecular Simulation"] by Tuckerman also gives a good description.

==Modifying IsingLattice.py==

Our IsingLattice object contains a stub function, '''<code>montecarlostep(T)</code>''', which takes a single argument — the temperature at which we want to do the cycle. There are also attributes used to calculate the average energy and magnetisation and respective variances: <code>self.E_tally</code>, <code>self.E2_tally</code>, <code>self.M_tally</code>, <code>self.M2_tally</code>, and <code>self.n_steps</code>. You should use these values to keep a running sum of all the energies, squared energies, magnetisations, and squared magnetisations that your system has sampled, as well as the number of steps that have been performed. Finally, there is also a function <code>statistics()</code>, that takes no arguments and should return the average energy, squared energy, magnetisation, squared magnetisation and current step, '''in that order'''.

'''<big>TASK 3b</big>: Implement a single step of the above algorithm in the <code>montecarlostep(T)</code> function. This function should return the energy of your lattice and the magnetisation at the end of the cycle. Complete the <code>statistics()</code> function. This should return the following quantities whenever it is called: &lang;''E''&rang;, &lang;''E''2&rang;, &lang;''M''&rang;, &lang;''M''2&rang;, and ''N''steps, the number of Monte Carlo steps that have elapsed. Put your code for these functions in the report.'''

''Note on units:'' In our simulation we will be using an unitless scaled energy <math>\tilde{E}=\frac{E}{J}</math>, and a unitless scaled temperature <math>\tilde{T}=\frac{k_B T}{J}</math>, where ''J'' is the interaction energy between spins discussed previously. The use of these scaled quantities in practice means that in the code we can set ''kB'' to numerically equal to 1. Going forward we will drop the &tilde; from the notation of energy and temperature.

You have been provided with a script, '''ILanim.py''', which will use your <code>IsingLattice</code> object to run a Monte Carlo simulation and display the output on the screen in real time. By default, this simulation will run an 8×8 lattice at ''T''=1.0, which is below the critical temperature. You will see a representation of the lattice, and graphs of the energy and magnetisation '''per spin'''. When you close the window, the script will use your '''<code>statistics()</code>''' function to print the average energy and magnetisation from the simulation. Note, you should run this script from the command line.

'''<big>TASK 3c</big>: If ''T'' < ''TC'', do you expect a spontaneous magnetisation (i.e. do you expect &lang;''M''&rang; ≠ 0)? When the state of the simulation appears to stop changing (when you have reached an equilibrium state), use the controls to export the output to an SVG or PNG file and attach this to your report. You should also include the output from your statistics() function.'''

'''<big>This is the third section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]]</big>'''

Programming a 2D Ising Model/Calculating the energy and magnetisation

2026-05-03T13:36:45Z

Gmallia:

'''<big>This is the second section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

= 2. Calculating the energy and magnetisation =

==Getting the files for the experiment==

'''Various scripts have been prepared to assist you with this experiment, and you can obtain them by downloading the files from the Blackboard.'''

We recommend you open the extracted folder with JupyterLab, which has an integrated terminal where you can run the simulation scripts.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. This class will be the basis for our model of the system, and is restricted to modelling the two dimensional case. It contains a number of stub methods that you will complete as you work through the experiment.

===A note about Python Classes===

The file IsingLattice.py makes use of a Python feature called ''class''. You should be familiar with the idea that all Python ''objects'' have an associated ''type'' -- you have already encountered objects such as integers, string, lists, and NumPy arrays. Classes allow us to define new types of object. In this case, we define a type called <code>IsingLattice</code>:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. We can see, for example, that each <code>IsingLattice</code> object has the methods <code>energy()</code>, <code>magnetisation()</code>, <code>montecarlostep(T)</code>, and <code>statistics()</code> (the method <code>__init__</code> will be explained shortly). Each of these methods takes a special first argument called self, which is a variable pointing to the current copy of IsingLattice. To use the new IsingLattice object that we have defined, we use code like:

<pre>
import IsingLattice

il = IsingLattice(5, 5) #create an IsingLattice object with 5 row and 5 columns

energy = il.energy() #call the energy() method of our IsingLattice object
</pre>

When writing code "within" the class definition, all attributes and methods belonging to the class must be accessed using the "self." notation. You can see an example of this in the <code>montecarlostep()</code> method:

<pre>
energy = self.energy()
</pre>

This line creates a variable called energy which contains whatever value was returned by the method <code>energy()</code> which is part of the <code>IsingLattice</code> class. If you just want to define a local variable instead, like an iterator in a loop, the <code>self.</code> notation is not required.

<pre>
def magnetisation(self):
# we have to use self.lattice, because the lattice variable is part of the object, not local to this function
for i in self.lattice:
...
# we use i, not self.i
# we don't need to access i outside this function
</pre>

Have a look at the other files that you have been given to see how this IsingLattice object is used.

===The Constructor===

The only method currently filled in is <code>__init__</code>. This is a special method called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this method is called automatically. It takes two arguments (excluding "self"), which are simply the number of rows and columns that our lattice will have. The attribute <code>self.lattice</code> is then created, which is a NumPy array with the requested number of rows and columns. The function <code>np.random.choice</code> is used to create a random collection of spins for our initial configuration. You can see the documentation for the function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html numpy.random.choice].

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

==Testing the files==

When you have completed the energy() and magnetisation() methods, we need to test that they work correctly.
We have provided you with a script, '''ILcheck.py''', which will create three IsingLattice objects and check their energies and magnetisations, displaying these on a graph. One configuration corresponds to the energy minimum, one to the energy maximum, and one to an intermediate state.

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

There are also some unit tests available in the <code>test_energy.py</code> file. These tests are designed to be run by simply running the [https://docs.pytest.org pytest] command in your terminal,
<pre>pytest</pre>
You should make sure your code passes the tests for the <code>energy()</code> and <code>magnetisation()</code> methods at this time. The file contains tests for methods that will only be implemented later on, so some failures are expected for now.
These sorts of test are useful because they are easier to write and faster to run and can be made to run automatically.

You can also debug and test your program on a Jupyter notebook or IPython console by importing your IsingLattice class, creating one object, and inspecting its properties. If you do this, we recommend that you start your session running the commands:

<pre>
%load_ext autoreload
%autoreload 2
</pre>

This will ensure that when you try executing your code, you always run the most recent version (please note this isn't fool proof and if you find problems you may need to restart the kernel)!

'''<big>This is the second section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

Programming a 2D Ising Model/Calculating the energy and magnetisation

2026-05-03T13:36:33Z

Gmallia:

'''<big>This is the second section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]].</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]].</big>'''

= 2. Calculating the energy and magnetisation =

==Getting the files for the experiment==

'''Various scripts have been prepared to assist you with this experiment, and you can obtain them by downloading the files from the Blackboard.'''

We recommend you open the extracted folder with JupyterLab, which has an integrated terminal where you can run the simulation scripts.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. This class will be the basis for our model of the system, and is restricted to modelling the two dimensional case. It contains a number of stub methods that you will complete as you work through the experiment.

===A note about Python Classes===

The file IsingLattice.py makes use of a Python feature called ''class''. You should be familiar with the idea that all Python ''objects'' have an associated ''type'' -- you have already encountered objects such as integers, string, lists, and NumPy arrays. Classes allow us to define new types of object. In this case, we define a type called <code>IsingLattice</code>:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. We can see, for example, that each <code>IsingLattice</code> object has the methods <code>energy()</code>, <code>magnetisation()</code>, <code>montecarlostep(T)</code>, and <code>statistics()</code> (the method <code>__init__</code> will be explained shortly). Each of these methods takes a special first argument called self, which is a variable pointing to the current copy of IsingLattice. To use the new IsingLattice object that we have defined, we use code like:

<pre>
import IsingLattice

il = IsingLattice(5, 5) #create an IsingLattice object with 5 row and 5 columns

energy = il.energy() #call the energy() method of our IsingLattice object
</pre>

When writing code "within" the class definition, all attributes and methods belonging to the class must be accessed using the "self." notation. You can see an example of this in the <code>montecarlostep()</code> method:

<pre>
energy = self.energy()
</pre>

This line creates a variable called energy which contains whatever value was returned by the method <code>energy()</code> which is part of the <code>IsingLattice</code> class. If you just want to define a local variable instead, like an iterator in a loop, the <code>self.</code> notation is not required.

<pre>
def magnetisation(self):
# we have to use self.lattice, because the lattice variable is part of the object, not local to this function
for i in self.lattice:
...
# we use i, not self.i
# we don't need to access i outside this function
</pre>

Have a look at the other files that you have been given to see how this IsingLattice object is used.

===The Constructor===

The only method currently filled in is <code>__init__</code>. This is a special method called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this method is called automatically. It takes two arguments (excluding "self"), which are simply the number of rows and columns that our lattice will have. The attribute <code>self.lattice</code> is then created, which is a NumPy array with the requested number of rows and columns. The function <code>np.random.choice</code> is used to create a random collection of spins for our initial configuration. You can see the documentation for the function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html numpy.random.choice].

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

==Testing the files==

When you have completed the energy() and magnetisation() methods, we need to test that they work correctly.
We have provided you with a script, '''ILcheck.py''', which will create three IsingLattice objects and check their energies and magnetisations, displaying these on a graph. One configuration corresponds to the energy minimum, one to the energy maximum, and one to an intermediate state.

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

There are also some unit tests available in the <code>test_energy.py</code> file. These tests are designed to be run by simply running the [https://docs.pytest.org pytest] command in your terminal,
<pre>pytest</pre>
You should make sure your code passes the tests for the <code>energy()</code> and <code>magnetisation()</code> methods at this time. The file contains tests for methods that will only be implemented later on, so some failures are expected for now.
These sorts of test are useful because they are easier to write and faster to run and can be made to run automatically.

You can also debug and test your program on a Jupyter notebook or IPython console by importing your IsingLattice class, creating one object, and inspecting its properties. If you do this, we recommend that you start your session running the commands:

<pre>
%load_ext autoreload
%autoreload 2
</pre>

This will ensure that when you try executing your code, you always run the most recent version (please note this isn't fool proof and if you find problems you may need to restart the kernel)!

'''<big>This is the second section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]].</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]].</big>'''

Programming a 2D Ising Model/Introduction to the Ising model

2026-05-03T13:35:17Z

Gmallia:

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]]</big>'''

= 1. Introduction to the Ising model =

== The Ising model ==
The Ising model is a simple physical model proposed in the 1920s to understand the behaviour of ferromagnets. In most materials, the magnetic dipoles of the atoms have no preferred orientation, and the bulk material has no net magnetic moment. In certain materials, however, there is preferential alignment of the magnetic dipoles along certain directions. These materials, such as iron, exhibit an overall magnetic moment, and are called ferromagnetic.

The phenomenon is based on competition between two physical effects — energy minimisation and entropy maximisation. On the one hand, parallel alignment of magnetic moments is energetically favourable. On the other hand, aligning the system in this way reduces the entropy.

In the Ising model, we consider a collection of spins, ''si'', each of which can either have the value "up" (+1) or "down" (-1), arranged on a lattice. In one dimension, this takes the form of a "chain" of adjacent spins, ''s1'' to ''sN''. In two dimensions we have a grid of cells, each of which contains one spin. In three dimensions, the system is a cuboid, and so on. The arrangements for 1, 2 and 3 dimensions are illustrated in '''figure 1'''.

[[File:ThirdYearCMPExpt-IsingSketch.png|300px|thumb|right|'''Figure 1''': Illustration of an Ising lattice in one (N=5), two (N=5×5), and three (N=5×5×5) dimensions. Red cells indicate the "up" spin", and blue cells indicate the "down" spin.]]

The rules of the model are simple:

* When no magnetic field is applied, the only energetic contribution comes from the interaction between spins in ''adjacent'' lattice cells
* The total interaction energy is defined as <math>- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>, where ''J'' is a constant which controls the strength of interaction and the notation ''j'' &in; neighbours(''i'') indicates that spin ''j'' must lie in a lattice cell adjacent to spin ''i''. The factor of one half is included to make sure that we do not count each spin-spin interaction twice.
* Periodic boundary conditions are applied. In other words, the "final" spin in a particular direction is considered to be adjacent to the "first" spin in that direction. In '''figure 2''', for example, cell 1 has neighbour interactions with cells 2, 3, 4, and 7, while cell 9 interacts with cells 3, 6, 7, and 8.

[[File:ThirdYearCMPExpt-IsingNumberedCell.png|300px|thumb|right|'''Figure 2''': Illustration of a 2D Ising lattice, with numbered cells to indicate the neighbour interactions.]]

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

==Phase Transitions==

The Ising model is interesting because it is one of the simplest physical models to display a phase transition (so long as we are in two dimensions or more). The temperature of the system controls the balance between the energetic driving force, which wants the system to adopt the lowest energy configuration that you imagined in the previous section, and the entropic driving force, which seeks to maximise the number of configurations available to the system.

'''<big>TASK 1b</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 (''D''=3, ''N''=1000)? How much entropy does the system gain by doing so?'''

At low temperatures, little thermal energy is available to overcome such energy barriers, and we expect the system to adopt the lowest energy (the energetic driving force will dominate). Under these conditions, we expect a large number of parallel spins, and we can characterise this by considering the total ''magnetisation'' of the system: <math>M = \sum_i s_i</math>. At higher temperatures, however, enough thermal energy will be available that spins can flip relatively freely. Under these conditions, the entropic driving force should dominate. The transition temperature between the two regimes is called the ''Curie temperature'', ''TC''. Below ''TC'', we observe the highly ordered magnetised phase, while above ''TC'' a disordered, demagnetised, phase emerges.

'''<big>TASK 1c</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 ''D''=3, ''N''=1000 at absolute zero?'''

We are interested, therefore, in determining the magnetisation ''M'', and energy ''E'', of the system as functions of temperature ''T''. We can write statistical mechanical equations for these quantities:

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

where ''α'' is used as a shorthand to represent all spins in the system, and ''Z'' is the partition function,
<math>Z = \sum_\alpha \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>.

It was discovered relatively quickly after the model was first proposed that the partition function in the one dimensional case can be solved analytically (though doing so is beyond the scope of this experiment):

<math>Z\left(T, N\right) = \left[2\cosh\left(\frac{J}{k_BT}\right)\right]^N</math>

Much later, [http://en.wikipedia.org/wiki/Lars_Onsager Lars Onsager] established that the two dimensional partition function can also be found analytically, though the mathematics is significantly more involved and the resulting expression is correspondingly much more complicated. In dimensions greater than two, no analytical solution is known.

Luckily for us, the problem can instead be tackled by numerical methods, and it is this that will be the focus of the experiment. We will discuss how the equations for ''M'' and ''E'' can be solved, but first we must set up a Python script to model the lattice.

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]]</big>'''

Programming a 2D Ising Model/Accelerating the code

2026-05-03T13:34:59Z

Gmallia: /* 4. Accelerating the code */

'''<big>This is the fourth section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of temperature|5. The effect of temperature]]</big>'''

= 4. Accelerating the code =
So far, your energy() and magnetisation() functions probably use Python loops to go over every spin in your system (and a double loop for the energy). One of the drawbacks of Python's flexibility is a performance penalty — though developing code in Python tends to be
faster than other compiled programming languages, the resulting program is much, much slower.

Fortunately, the NumPy module provides an interface to a number of functions written to give high performance.
These are written in the C programming language under the hood and take numpy arrays, which are contiguous arrays of a single data type,
unlike Python lists, which may contain mixes of data types.

In this section, we are going to use these functions to improve the performance of our simulation.

When trying to optimise the performance of something, you should always measure first, rather than making assumptions. Although there are some very obvious candidates in our program, you should always measure properly.

The <code>ILtimetrial.py</code> script measures how fast it takes for 2000 iterations of your code to run.

The <code>ILProfiling.py</code> script is the same, except uses the python standard library package [https://docs.python.org/3/library/profile.html cProfile] to sample what functions are being called and give you information
about the time spent in each. Analysing the output of the script, you will see that by far the majority of the time is spent in the energy method, followed by the magnetisation one. This makes sense, as each Monte-Carlo step can call the energy function multiple times.

(If you want to profile your code in more detail, please seek advice from a demonstrator.)

In order to improve the performance of this, we need to:
* Run the energy call fewer times per step.
* Improve the performance of the energy function.

Fortunately, the NumPy module provides an interface to a number of functions written to give high performance. In this section, we are going to use these functions to improve the performance of our simulation.

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

'''<big>TASK 4c</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! What makes this new version faster?'''

'''<big>Task 4d: Making the code even faster </big>: Consider trying to make the Monte Carlo step function substantially faster by calculating the difference in energy directly, in a way that avoids looping over the whole lattice. Think back how you answered some of the introduction questions. You can make it so the energy and magnetisation only needs to be calculated explicitly at the beginning of the simulation.'''

'''<big>This is the fourth section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of temperature|5. The effect of temperature]]</big>'''

Programming a 2D Ising Model/Accelerating the code

2026-05-03T13:34:34Z

Gmallia:

'''<big>This is the fourth section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of temperature|5. The effect of temperature]]</big>'''

= 4. Accelerating the code =
So far, your energy() and magnetisation() functions probably use Python loops to go over every spin in your system (and a double loop for the energy). One of the drawbacks of Python's flexibility is a performance penalty — though developing code in Python tends to be
faster than other compiled programming languages, the resulting program is much, much slower.

Fortunately, the NumPy module provides an interface to a number of functions written to give high performance.
These are written in the C programming language under the hood and take numpy arrays, which are contiguous arrays of a single data type,
unlike Python lists, which may contain mixes of data types.

In this section, we are going to use these functions to improve the performance of our simulation.

When trying to optimise the performance of something, you should always measure first, rather than making assumptions. Although there are some very obvious candidates in our program, you should always measure properly.

The <code>ILtimetrial.py</code> script measures how fast it takes for 2000 iterations of your code to run.

The <code>ILProfiling.py</code> script is the same, except uses the python standard library package [https://docs.python.org/3/library/profile.html cProfile] to sample what functions are being called and give you information
about the time spent in each. Analysing the output of the script, you will see that by far the majority of the time is spent in the energy method, followed by the magnetisation one. This makes sense, as each Monte-Carlo step can call the energy function multiple times.

(If you want to profile your code in more detail, please seek advice from a demonstrator.)

In order to improve the performance of this, we need to:
* Run the energy call fewer times per step.
* Improve the performance of the energy function.

Fortunately, the NumPy module provides an interface to a number of functions written to give high performance. In this section, we are going to use these functions to improve the performance of our simulation.

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

'''<big>TASK 4c</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! What makes this new version faster?'''

'''<big>Task 4d: Making the code even faster </big>: Consider trying to make the Monte Carlo step function substantially faster by calculating the difference in energy directly, in a way that avoids looping over the whole lattice. Think back how you answered some of the introduction questions. You can make it so the energy and magnetisation only needs to be calculated explicitly at the beginning of the simulation.'''

'''<big>This is the fourth section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of temperature|5. The effect of temperature]]</big>'''

Programming a 2D Ising Model/Accelerating the code

2026-05-03T13:33:53Z

Gmallia:

'''<big>This is the fourth section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of temperature|5. The effect of temperature]]</big>'''

= 4. Accelerating the code =
So far, your energy() and magnetisation() functions probably use Python loops to go over every spin in your system (and a double loop for the energy). One of the drawbacks of Python's flexibility is a performance penalty — though developing code in Python tends to be
faster than other compiled programming languages, the resulting program is much, much slower.

Fortunately, the NumPy module provides an interface to a number of functions written to give high performance.
These are written in the C programming language under the hood and take numpy arrays, which are contiguous arrays of a single data type,
unlike Python lists, which may contain mixes of data types.

In this section, we are going to use these functions to improve the performance of our simulation.

When trying to optimise the performance of something, you should always measure first, rather than making assumptions. Although there are some very obvious candidates in our program, you should always measure properly.

The <code>ILtimetrial.py</code> script measures how fast it takes for 2000 iterations of your code to run.

The <code>ILProfiling.py</code> script is the same, except uses the python standard library package [https://docs.python.org/3/library/profile.html cProfile] to sample what functions are being called and give you information
about the time spent in each. Analysing the output of the script, you will see that by far the majority of the time is spent in the energy method, followed by the magnetisation one. This makes sense, as each Monte-Carlo step can call the energy function multiple times.

(If you want to profile your code in more detail, please seek advice from a demonstrator.)

In order to improve the performance of this, we need to:
* Run the energy call fewer times per step.
* Improve the performance of the energy function.

Fortunately, the NumPy module provides an interface to a number of functions written to give high performance. In this section, we are going to use these functions to improve the performance of our simulation.

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

'''<big>TASK 4c</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! What makes this new version faster?'''

'''<big> Optional: Making the code even faster </big>: Consider trying to make the Monte Carlo step function substantially faster by calculating the difference in energy directly, in a way that avoids looping over the whole lattice. Think back how you answered some of the introduction questions. You can make it so the energy and magnetisation only needs to be calculated explicitly at the beginning of the simulation.'''

'''<big>This is the fourth section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]]</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/The effect of temperature|5. The effect of temperature]]</big>'''

Programming a 2D Ising Model/Calculating the energy and magnetisation

2026-05-03T13:29:52Z

Gmallia: /* Testing the files */

'''<big>This is the second section of the Ising model experiment.</big>'''

'''<big>You can return to the previous page, [[Programming a 2D Ising Model/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

= 2. Calculating the energy and magnetisation =

==Getting the files for the experiment==

'''Various scripts have been prepared to assist you with this experiment, and you can obtain them by downloading the files from the Blackboard.'''

We recommend you open the extracted folder with JupyterLab, which has an integrated terminal where you can run the simulation scripts.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. This class will be the basis for our model of the system, and is restricted to modelling the two dimensional case. It contains a number of stub methods that you will complete as you work through the experiment.

===A note about Python Classes===

The file IsingLattice.py makes use of a Python feature called ''class''. You should be familiar with the idea that all Python ''objects'' have an associated ''type'' -- you have already encountered objects such as integers, string, lists, and NumPy arrays. Classes allow us to define new types of object. In this case, we define a type called <code>IsingLattice</code>:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. We can see, for example, that each <code>IsingLattice</code> object has the methods <code>energy()</code>, <code>magnetisation()</code>, <code>montecarlostep(T)</code>, and <code>statistics()</code> (the method <code>__init__</code> will be explained shortly). Each of these methods takes a special first argument called self, which is a variable pointing to the current copy of IsingLattice. To use the new IsingLattice object that we have defined, we use code like:

<pre>
import IsingLattice

il = IsingLattice(5, 5) #create an IsingLattice object with 5 row and 5 columns

energy = il.energy() #call the energy() method of our IsingLattice object
</pre>

When writing code "within" the class definition, all attributes and methods belonging to the class must be accessed using the "self." notation. You can see an example of this in the <code>montecarlostep()</code> method:

<pre>
energy = self.energy()
</pre>

This line creates a variable called energy which contains whatever value was returned by the method <code>energy()</code> which is part of the <code>IsingLattice</code> class. If you just want to define a local variable instead, like an iterator in a loop, the <code>self.</code> notation is not required.

<pre>
def magnetisation(self):
# we have to use self.lattice, because the lattice variable is part of the object, not local to this function
for i in self.lattice:
...
# we use i, not self.i
# we don't need to access i outside this function
</pre>

Have a look at the other files that you have been given to see how this IsingLattice object is used.

===The Constructor===

The only method currently filled in is <code>__init__</code>. This is a special method called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this method is called automatically. It takes two arguments (excluding "self"), which are simply the number of rows and columns that our lattice will have. The attribute <code>self.lattice</code> is then created, which is a NumPy array with the requested number of rows and columns. The function <code>np.random.choice</code> is used to create a random collection of spins for our initial configuration. You can see the documentation for the function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html numpy.random.choice].

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

==Testing the files==

When you have completed the energy() and magnetisation() methods, we need to test that they work correctly.
We have provided you with a script, '''ILcheck.py''', which will create three IsingLattice objects and check their energies and magnetisations, displaying these on a graph. One configuration corresponds to the energy minimum, one to the energy maximum, and one to an intermediate state.

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

There are also some unit tests available in the <code>test_energy.py</code> file. These tests are designed to be run by simply running the [https://docs.pytest.org pytest] command in your terminal,
<pre>pytest</pre>
You should make sure your code passes the tests for the <code>energy()</code> and <code>magnetisation()</code> methods at this time. The file contains tests for methods that will only be implemented later on, so some failures are expected for now.
These sorts of test are useful because they are easier to write and faster to run and can be made to run automatically.

You can also debug and test your program on a Jupyter notebook or IPython console by importing your IsingLattice class, creating one object, and inspecting its properties. If you do this, we recommend that you start your session running the commands:

<pre>
%load_ext autoreload
%autoreload 2
</pre>

This will ensure that when you try executing your code, you always run the most recent version (please note this isn't fool proof and if you find problems you may need to restart the kernel)!

'''<big>This is the second section of the Ising model experiment.</big>'''

'''<big>Main Page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Programming a 2D Ising Model/Introduction to the Ising model|1. Introduction to the Ising model]].</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]].</big>'''

Programming a 2D Ising Model/Calculating the energy and magnetisation

2026-05-03T12:39:52Z

Gmallia:

'''<big>This is the second section of the Ising model experiment.</big>'''

'''<big>You can return to the previous page, [[Programming a 2D Ising Model/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

= 2. Calculating the energy and magnetisation =

==Getting the files for the experiment==

'''Various scripts have been prepared to assist you with this experiment, and you can obtain them by downloading the files from the Blackboard.'''

We recommend you open the extracted folder with JupyterLab, which has an integrated terminal where you can run the simulation scripts.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. This class will be the basis for our model of the system, and is restricted to modelling the two dimensional case. It contains a number of stub methods that you will complete as you work through the experiment.

===A note about Python Classes===

The file IsingLattice.py makes use of a Python feature called ''class''. You should be familiar with the idea that all Python ''objects'' have an associated ''type'' -- you have already encountered objects such as integers, string, lists, and NumPy arrays. Classes allow us to define new types of object. In this case, we define a type called <code>IsingLattice</code>:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. We can see, for example, that each <code>IsingLattice</code> object has the methods <code>energy()</code>, <code>magnetisation()</code>, <code>montecarlostep(T)</code>, and <code>statistics()</code> (the method <code>__init__</code> will be explained shortly). Each of these methods takes a special first argument called self, which is a variable pointing to the current copy of IsingLattice. To use the new IsingLattice object that we have defined, we use code like:

<pre>
import IsingLattice

il = IsingLattice(5, 5) #create an IsingLattice object with 5 row and 5 columns

energy = il.energy() #call the energy() method of our IsingLattice object
</pre>

When writing code "within" the class definition, all attributes and methods belonging to the class must be accessed using the "self." notation. You can see an example of this in the <code>montecarlostep()</code> method:

<pre>
energy = self.energy()
</pre>

This line creates a variable called energy which contains whatever value was returned by the method <code>energy()</code> which is part of the <code>IsingLattice</code> class. If you just want to define a local variable instead, like an iterator in a loop, the <code>self.</code> notation is not required.

<pre>
def magnetisation(self):
# we have to use self.lattice, because the lattice variable is part of the object, not local to this function
for i in self.lattice:
...
# we use i, not self.i
# we don't need to access i outside this function
</pre>

Have a look at the other files that you have been given to see how this IsingLattice object is used.

===The Constructor===

The only method currently filled in is <code>__init__</code>. This is a special method called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this method is called automatically. It takes two arguments (excluding "self"), which are simply the number of rows and columns that our lattice will have. The attribute <code>self.lattice</code> is then created, which is a NumPy array with the requested number of rows and columns. The function <code>np.random.choice</code> is used to create a random collection of spins for our initial configuration. You can see the documentation for the function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html numpy.random.choice].

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

==Testing the files==

When you have completed the energy() and magnetisation() methods, we need to test that they work correctly.
We have provided you with a script, '''ILcheck.py''', which will create three IsingLattice objects and check their energies and magnetisations, displaying these on a graph. One configuration corresponds to the energy minimum, one to the energy maximum, and one to an intermediate state.

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

There are also some unit tests available in the <code>test_energy.py</code> file. These tests are designed to be run by simply running the [https://docs.pytest.org pytest] command in your terminal,
<pre>pytest</pre>
You should make sure your code passes the tests for the <code>energy()</code> and <code>magnetisation()</code> methods at this time. The file contains tests for methods that will only be implemented later on, so some failures are expected for now.
These sorts of test are useful because they are easier to write and faster to run and can be made to run automatically.

You can also debug and test your program on a Jupyter notebook or IPython console by importing your IsingLattice class, creating one object, and inspecting its properties. If you do this, we recommend that you start your session running the commands:

<pre>
%load_ext autoreload
%autoreload 2
</pre>

This will ensure that when you try executing your code, you always run the most recent version (please note this isn't fool proof and if you find problems you may need to restart the kernel)!

'''<big>This is the second section of the Ising model experiment. You can return to the previous page, [[Programming a 2D Ising Model/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]].</big>'''

Programming a 2D Ising Model

2026-05-03T12:37:12Z

Gmallia:

Programming a 2D Ising Model

2026-05-03T12:36:49Z

Gmallia:

'''<big>This is the Ising Model experiment (aka programming for simple simulations).</big>'''

'''<big>This is a core experiment for year 3 students on the Chemistry with Molecular Physics degree.</big>'''

= Programming a 2D Ising Model =

==Introduction==

In this exercise, you are going to use the Python that you learned in the first year to write a code to perform Monte-Carlo simulations of the 2D [http://en.wikipedia.org/wiki/Ising_model Ising model], a set of spins on a lattice which is used to model ferromagnetic behaviour, and also to analyse the results of the simulation to find the heat capacity of the system and the Curie temperature — the temperature below which the system is able to maintain a spontaneous magnetisation.

A more extensive introduction to the Ising Model will be given in the next page. This page gives some practical information about how the experiment.

==Assessment==
At the end of this experiment you must submit a report, writing up your findings. Each section of the experiment has a number of tasks that you should complete, labelled '''<big>TASK</big>'''. If this is a mathematical exercise, your report should contain a short summary of the solution. If it is a graphical exercise, you should include the relevant image. Your report should explain briefly what you did in each stage of the experiment, and what your findings were.

The report should be submitted as a PDF file to Turnitin, created with your favourite word processor or LaTeX, and should contain the annotated source code of the Python scripts that your write during the experiment. For clarity, code included in the report should be written in a monospaced font (e.g. Courier), which in LaTeX can be easily achieved by placing your code inside a <code>\texttt{}</code> command. '''The code, except where relevant to discussion, can be placed in an appendix.'''

As well as including relevant source code in your report, you should also submit a bundle with all your scripts and analysis code as a separate Blackboard submission.

===On the use of AI tools===
The code you write for this experiment is assessed, and as per the [https://bb.imperial.ac.uk/bbcswebdav/pid-3344481-dt-content-rid-21115105_1/xid-21115105_1 department's guidelines on the use of AI], it is not acceptable to use generative AI tools to write computer code for this assignment.

Writing computer programs yourself enhances your logical thinking and problem solving skills. Using generative AI tools to write your code makes you dumb. Don't choose to be dumb, don't use AI to write code for you.

==Getting Help==

The Graduate Teaching Assistants demonstrating this experiment are Fionn Carman and Lei Li. Demonstrators will be available in room 232A at the following times:

{| class="wikitable"
|+ Demonstrator timetable
! !! Monday !! Tuesday !! Wednesday !! Thusrsday !! Friday
|-
| Morning || 10:00 - Experiment intro 10:00-13:00 - Li || 10:00-13:00 - Fionn|| || 10:00-12:00 - João || 10:00-13:00 - Fionn
|-
| Afternoon || 14:00-17:00 - Li|| 14:00-17:00 - Fionn|| || || 14:00-17:00 - Fionn
|}

If you have questions outside these times, please use the Discussion Board on Blackboard.

João and Giuseppe are the members of academic staff responsible for this experiment. You are welcome to address your questions and comments to them.

==Structure of this Experiment==

This experimental manual has been broken up into a number of subsections. Direct links to each of them may be found below. You should attempt them in order, and you should complete all of them to finish the experiment.

# [[Programming_a_2D_Ising_Model/Introduction_to_the_Ising_model|Introduction to the Ising model]]
# [[Programming_a_2D_Ising_Model/Calculating the energy and magnetisation|Calculating the energy and magnetisation]]
# [[Programming_a_2D_Ising_Model/Introduction to Monte Carlo simulation|Introduction to the Monte Carlo simulation]]
# [[Programming_a_2D_Ising_Model/Accelerating the code|Accelerating the code]]
# [[Programming_a_2D_Ising_Model/The effect of temperature|The effect of temperature]]
# [[Programming_a_2D_Ising_Model/The effect of system size|The effect of system size]]
# [[Programming_a_2D_Ising_Model/Determining the heat capacity|Determining the heat capacity]]
# [[Programming_a_2D_Ising_Model/Locating the Curie temperature|Locating the Curie temperature]]

Programming a 2D Ising Model

2026-05-03T12:36:07Z

Gmallia:

'''<big>This is the Ising Model experiment (aka programming for simple simulations).</big>'''

'''<big>This is a core experiment for year 3 students on the Chemistry with Molecular Physics degree.</big>'''

==Introduction==

In this exercise, you are going to use the Python that you learned in the first year to write a code to perform Monte-Carlo simulations of the 2D [http://en.wikipedia.org/wiki/Ising_model Ising model], a set of spins on a lattice which is used to model ferromagnetic behaviour, and also to analyse the results of the simulation to find the heat capacity of the system and the Curie temperature — the temperature below which the system is able to maintain a spontaneous magnetisation.

A more extensive introduction to the Ising Model will be given in the next page. This page gives some practical information about how the experiment.

==Assessment==
At the end of this experiment you must submit a report, writing up your findings. Each section of the experiment has a number of tasks that you should complete, labelled '''<big>TASK</big>'''. If this is a mathematical exercise, your report should contain a short summary of the solution. If it is a graphical exercise, you should include the relevant image. Your report should explain briefly what you did in each stage of the experiment, and what your findings were.

The report should be submitted as a PDF file to Turnitin, created with your favourite word processor or LaTeX, and should contain the annotated source code of the Python scripts that your write during the experiment. For clarity, code included in the report should be written in a monospaced font (e.g. Courier), which in LaTeX can be easily achieved by placing your code inside a <code>\texttt{}</code> command. '''The code, except where relevant to discussion, can be placed in an appendix.'''

As well as including relevant source code in your report, you should also submit a bundle with all your scripts and analysis code as a separate Blackboard submission.

===On the use of AI tools===
The code you write for this experiment is assessed, and as per the [https://bb.imperial.ac.uk/bbcswebdav/pid-3344481-dt-content-rid-21115105_1/xid-21115105_1 department's guidelines on the use of AI], it is not acceptable to use generative AI tools to write computer code for this assignment.

Writing computer programs yourself enhances your logical thinking and problem solving skills. Using generative AI tools to write your code makes you dumb. Don't choose to be dumb, don't use AI to write code for you.

==Getting Help==

The Graduate Teaching Assistants demonstrating this experiment are Fionn Carman and Lei Li. Demonstrators will be available in room 232A at the following times:

{| class="wikitable"
|+ Demonstrator timetable
! !! Monday !! Tuesday !! Wednesday !! Thusrsday !! Friday
|-
| Morning || 10:00 - Experiment intro 10:00-13:00 - Li || 10:00-13:00 - Fionn|| || 10:00-12:00 - João || 10:00-13:00 - Fionn
|-
| Afternoon || 14:00-17:00 - Li|| 14:00-17:00 - Fionn|| || || 14:00-17:00 - Fionn
|}

If you have questions outside these times, please use the Discussion Board on Blackboard.

João and Giuseppe are the members of academic staff responsible for this experiment. You are welcome to address your questions and comments to them.

==Structure of this Experiment==

This experimental manual has been broken up into a number of subsections. Direct links to each of them may be found below. You should attempt them in order, and you should complete all of them to finish the experiment.

# [[Programming_a_2D_Ising_Model/Introduction_to_the_Ising_model|Introduction to the Ising model]]
# [[Programming_a_2D_Ising_Model/Calculating the energy and magnetisation|Calculating the energy and magnetisation]]
# [[Programming_a_2D_Ising_Model/Introduction to Monte Carlo simulation|Introduction to the Monte Carlo simulation]]
# [[Programming_a_2D_Ising_Model/Accelerating the code|Accelerating the code]]
# [[Programming_a_2D_Ising_Model/The effect of temperature|The effect of temperature]]
# [[Programming_a_2D_Ising_Model/The effect of system size|The effect of system size]]
# [[Programming_a_2D_Ising_Model/Determining the heat capacity|Determining the heat capacity]]
# [[Programming_a_2D_Ising_Model/Locating the Curie temperature|Locating the Curie temperature]]

Programming a 2D Ising Model/Introduction to Monte Carlo simulation

2026-05-03T12:34:45Z

Gmallia:

'''<big>This is the third section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]].</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]].</big>'''

= 3. Introduction to Monte Carlo simulation =

==Average Energy and Magnetisation==

Consider again the expressions for the average energy and magnetisation that we gave in the introduction.

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

Imagine we want to evaluate these at a particular temperature, in a system of 100 spins.

'''<big>TASK 3a</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 1×109 configurations per second with our computer. How long will it take to evaluate a single value of &lang;''M''&rang;''T''?'''

Clearly, we need to try a cleverer approach.

For the overwhelming majority of the possible configurations, the Boltzmann factor, <math>\exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\}</math>, will be very small, and that state will contribute very little to the average value. We can save an enormous amount of time, and make this problem tractable, if we only consider the states whose Boltzmann factor is not so vanishingly small. This technique is called ''importance sampling'' — instead of sampling every point in phase space, we sample only those that the system is likely to occupy.

==Importance Sampling==

This is easily stated, of course, but how can we know which states contribute a lot to the average without actually calculating the energy? So far, we have been imagining generating random states of the system; that is to say, we have been choosing ''α'' from a uniform distribution in which no arrangement of spins is any more likely to be chosen than any other. If we instead had a method which generated states randomly from the probability distribution <math>\exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\}</math>, then our problem would be solved. The method which allows us to do this was developed in the 1950s, and is called the [https://doi.org/10.1063/1.1699114 Metropolis Monte Carlo (MMC)] method, or more often simply "Monte Carlo simulation". This was one of the major breakthroughs of the early days of computational science, and in the year 2000, the IEEE named among the [https://doi.org/10.1109/MCISE.2000.814652 "Top 10 Algorithms of the 20th Century"].

The algorithm is as follows:

# Start from a given configuration of spins, ''α0'', with energy ''E0''.
# Choose a single spin '''at random''', and "flip" it, to generate a new configuration ''α1''.
# Calculate the energy of this new configuration, ''E1''.
# Calculate the energy difference between the states, ''ΔE = E1 - E0''.
## If the ''ΔE ≤ 0'' (the spin flipping decreased the energy), then we '''accept''' the new configuration.
##* We set ''α0 = α1'', and ''E0 = E1'', and then '''go to step 5'''.
## If ''ΔE > 0'', the spin flipping increased the energy. By considering the probability of observing the starting and final states, ''α0'' and ''α1'', it can be shown that the probability for the transition between the two to occur is <math>\exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>. To ensure that we only accept this kind of spin flip with the correct probability, we use the following procedure:
### Choose a random number, ''R'', in the interval [0,1).
### If <math>R \leq \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''accept''' the new configuration.
###* We set ''α0 = α1'', and ''E0 = E1'', and then '''go to step 5'''.
### If <math>R > \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''reject''' the new configuration.
###* ''α0'' and ''E0'' are left unchanged. '''Go to step 5'''.
# Update the running averages of the energy and magnetisation.
# Monte Carlo "step" complete, '''return to step 2'''.

Step 4 is the key to the method. By accepting moves which increase the energy only with a certain probability, we ensure that the sequence of states that we generate is correctly distributed. A transition which carries a very large energy penalty, ''ΔE'' is extremely unlikely to be selected. The use of random numbers in this step is the reason that the method acquired the name "Monte Carlo", after the casinos located there!

Note, unlike molecular dynamics, each step only moves a single atom. In Monte Carlo simulation, there is the notion of a cycle which is ''N''spins-steps, so that each spin on average has an attempted flip per cycle.

If you are interested in the mathematical details of why this procedure generates a sequence of states distributed in the correct way, consult the Monte Carlo chapter of [https://library-search.imperial.ac.uk/permalink/f/o1297h/TN_cdi_askewsholts_vlebooks_9780080519982 "Understanding Molecular Simulation"], by Frenkel and Smit, but a discussion of this is not required for the experiment. Chapter 7 of [https://library-search.imperial.ac.uk/permalink/44IMP_INST/mek6kh/alma997227534401591 "Statistical Mechanics: Theory and Molecular Simulation"] by Tuckerman also gives a good description.

==Modifying IsingLattice.py==

Our IsingLattice object contains a stub function, '''<code>montecarlostep(T)</code>''', which takes a single argument — the temperature at which we want to do the cycle. There are also attributes used to calculate the average energy and magnetisation and respective variances: <code>self.E_tally</code>, <code>self.E2_tally</code>, <code>self.M_tally</code>, <code>self.M2_tally</code>, and <code>self.n_steps</code>. You should use these values to keep a running sum of all the energies, squared energies, magnetisations, and squared magnetisations that your system has sampled, as well as the number of steps that have been performed. Finally, there is also a function <code>statistics()</code>, that takes no arguments and should return the average energy, squared energy, magnetisation, squared magnetisation and current step, '''in that order'''.

'''<big>TASK 3b</big>: Implement a single step of the above algorithm in the <code>montecarlostep(T)</code> function. This function should return the energy of your lattice and the magnetisation at the end of the cycle. Complete the <code>statistics()</code> function. This should return the following quantities whenever it is called: &lang;''E''&rang;, &lang;''E''2&rang;, &lang;''M''&rang;, &lang;''M''2&rang;, and ''N''steps, the number of Monte Carlo steps that have elapsed. Put your code for these functions in the report.'''

''Note on units:'' In our simulation we will be using an unitless scaled energy <math>\tilde{E}=\frac{E}{J}</math>, and a unitless scaled temperature <math>\tilde{T}=\frac{k_B T}{J}</math>, where ''J'' is the interaction energy between spins discussed previously. The use of these scaled quantities in practice means that in the code we can set ''kB'' to numerically equal to 1. Going forward we will drop the &tilde; from the notation of energy and temperature.

You have been provided with a script, '''ILanim.py''', which will use your <code>IsingLattice</code> object to run a Monte Carlo simulation and display the output on the screen in real time. By default, this simulation will run an 8×8 lattice at ''T''=1.0, which is below the critical temperature. You will see a representation of the lattice, and graphs of the energy and magnetisation '''per spin'''. When you close the window, the script will use your '''<code>statistics()</code>''' function to print the average energy and magnetisation from the simulation. Note, you should run this script from the command line.

'''<big>TASK 3c</big>: If ''T'' < ''TC'', do you expect a spontaneous magnetisation (i.e. do you expect &lang;''M''&rang; ≠ 0)? When the state of the simulation appears to stop changing (when you have reached an equilibrium state), use the controls to export the output to an SVG or PNG file and attach this to your report. You should also include the output from your statistics() function.'''

'''<big>This is the third section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]].</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]].</big>'''

Programming a 2D Ising Model/Introduction to Monte Carlo simulation

2026-05-03T12:32:35Z

Gmallia: /* Modifying IsingLattice.py */

'''<big>This is the third section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]].</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]].</big>'''
==Average Energy and Magnetisation==

Consider again the expressions for the average energy and magnetisation that we gave in the introduction.

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

Imagine we want to evaluate these at a particular temperature, in a system of 100 spins.

'''<big>TASK 3a</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 1×109 configurations per second with our computer. How long will it take to evaluate a single value of &lang;''M''&rang;''T''?'''

Clearly, we need to try a cleverer approach.

For the overwhelming majority of the possible configurations, the Boltzmann factor, <math>\exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\}</math>, will be very small, and that state will contribute very little to the average value. We can save an enormous amount of time, and make this problem tractable, if we only consider the states whose Boltzmann factor is not so vanishingly small. This technique is called ''importance sampling'' — instead of sampling every point in phase space, we sample only those that the system is likely to occupy.

==Importance Sampling==

This is easily stated, of course, but how can we know which states contribute a lot to the average without actually calculating the energy? So far, we have been imagining generating random states of the system; that is to say, we have been choosing ''α'' from a uniform distribution in which no arrangement of spins is any more likely to be chosen than any other. If we instead had a method which generated states randomly from the probability distribution <math>\exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\}</math>, then our problem would be solved. The method which allows us to do this was developed in the 1950s, and is called the [https://doi.org/10.1063/1.1699114 Metropolis Monte Carlo (MMC)] method, or more often simply "Monte Carlo simulation". This was one of the major breakthroughs of the early days of computational science, and in the year 2000, the IEEE named among the [https://doi.org/10.1109/MCISE.2000.814652 "Top 10 Algorithms of the 20th Century"].

The algorithm is as follows:

# Start from a given configuration of spins, ''α0'', with energy ''E0''.
# Choose a single spin '''at random''', and "flip" it, to generate a new configuration ''α1''.
# Calculate the energy of this new configuration, ''E1''.
# Calculate the energy difference between the states, ''ΔE = E1 - E0''.
## If the ''ΔE ≤ 0'' (the spin flipping decreased the energy), then we '''accept''' the new configuration.
##* We set ''α0 = α1'', and ''E0 = E1'', and then '''go to step 5'''.
## If ''ΔE > 0'', the spin flipping increased the energy. By considering the probability of observing the starting and final states, ''α0'' and ''α1'', it can be shown that the probability for the transition between the two to occur is <math>\exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>. To ensure that we only accept this kind of spin flip with the correct probability, we use the following procedure:
### Choose a random number, ''R'', in the interval [0,1).
### If <math>R \leq \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''accept''' the new configuration.
###* We set ''α0 = α1'', and ''E0 = E1'', and then '''go to step 5'''.
### If <math>R > \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''reject''' the new configuration.
###* ''α0'' and ''E0'' are left unchanged. '''Go to step 5'''.
# Update the running averages of the energy and magnetisation.
# Monte Carlo "step" complete, '''return to step 2'''.

Step 4 is the key to the method. By accepting moves which increase the energy only with a certain probability, we ensure that the sequence of states that we generate is correctly distributed. A transition which carries a very large energy penalty, ''ΔE'' is extremely unlikely to be selected. The use of random numbers in this step is the reason that the method acquired the name "Monte Carlo", after the casinos located there!

Note, unlike molecular dynamics, each step only moves a single atom. In Monte Carlo simulation, there is the notion of a cycle which is ''N''spins-steps, so that each spin on average has an attempted flip per cycle.

If you are interested in the mathematical details of why this procedure generates a sequence of states distributed in the correct way, consult the Monte Carlo chapter of [https://library-search.imperial.ac.uk/permalink/f/o1297h/TN_cdi_askewsholts_vlebooks_9780080519982 "Understanding Molecular Simulation"], by Frenkel and Smit, but a discussion of this is not required for the experiment. Chapter 7 of [https://library-search.imperial.ac.uk/permalink/44IMP_INST/mek6kh/alma997227534401591 "Statistical Mechanics: Theory and Molecular Simulation"] by Tuckerman also gives a good description.

==Modifying IsingLattice.py==

Our IsingLattice object contains a stub function, '''<code>montecarlostep(T)</code>''', which takes a single argument — the temperature at which we want to do the cycle. There are also attributes used to calculate the average energy and magnetisation and respective variances: <code>self.E_tally</code>, <code>self.E2_tally</code>, <code>self.M_tally</code>, <code>self.M2_tally</code>, and <code>self.n_steps</code>. You should use these values to keep a running sum of all the energies, squared energies, magnetisations, and squared magnetisations that your system has sampled, as well as the number of steps that have been performed. Finally, there is also a function <code>statistics()</code>, that takes no arguments and should return the average energy, squared energy, magnetisation, squared magnetisation and current step, '''in that order'''.

'''<big>TASK 3b</big>: Implement a single step of the above algorithm in the <code>montecarlostep(T)</code> function. This function should return the energy of your lattice and the magnetisation at the end of the cycle. Complete the <code>statistics()</code> function. This should return the following quantities whenever it is called: &lang;''E''&rang;, &lang;''E''2&rang;, &lang;''M''&rang;, &lang;''M''2&rang;, and ''N''steps, the number of Monte Carlo steps that have elapsed. Put your code for these functions in the report.'''

''Note on units:'' In our simulation we will be using an unitless scaled energy <math>\tilde{E}=\frac{E}{J}</math>, and a unitless scaled temperature <math>\tilde{T}=\frac{k_B T}{J}</math>, where ''J'' is the interaction energy between spins discussed previously. The use of these scaled quantities in practice means that in the code we can set ''kB'' to numerically equal to 1. Going forward we will drop the &tilde; from the notation of energy and temperature.

You have been provided with a script, '''ILanim.py''', which will use your <code>IsingLattice</code> object to run a Monte Carlo simulation and display the output on the screen in real time. By default, this simulation will run an 8×8 lattice at ''T''=1.0, which is below the critical temperature. You will see a representation of the lattice, and graphs of the energy and magnetisation '''per spin'''. When you close the window, the script will use your '''<code>statistics()</code>''' function to print the average energy and magnetisation from the simulation. Note, you should run this script from the command line.

'''<big>TASK 3c</big>: If ''T'' < ''TC'', do you expect a spontaneous magnetisation (i.e. do you expect &lang;''M''&rang; ≠ 0)? When the state of the simulation appears to stop changing (when you have reached an equilibrium state), use the controls to export the output to an SVG or PNG file and attach this to your report. You should also include the output from your statistics() function.'''

'''<big>This is the third section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Previous section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]].</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]].</big>'''

Programming a 2D Ising Model/Introduction to the Ising model

2026-05-03T11:36:42Z

Gmallia:

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]].</big>'''

= 1. Introduction to the Ising model =

== The Ising model ==
The Ising model is a simple physical model proposed in the 1920s to understand the behaviour of ferromagnets. In most materials, the magnetic dipoles of the atoms have no preferred orientation, and the bulk material has no net magnetic moment. In certain materials, however, there is preferential alignment of the magnetic dipoles along certain directions. These materials, such as iron, exhibit an overall magnetic moment, and are called ferromagnetic.

The phenomenon is based on competition between two physical effects — energy minimisation and entropy maximisation. On the one hand, parallel alignment of magnetic moments is energetically favourable. On the other hand, aligning the system in this way reduces the entropy.

In the Ising model, we consider a collection of spins, ''si'', each of which can either have the value "up" (+1) or "down" (-1), arranged on a lattice. In one dimension, this takes the form of a "chain" of adjacent spins, ''s1'' to ''sN''. In two dimensions we have a grid of cells, each of which contains one spin. In three dimensions, the system is a cuboid, and so on. The arrangements for 1, 2 and 3 dimensions are illustrated in '''figure 1'''.

[[File:ThirdYearCMPExpt-IsingSketch.png|300px|thumb|right|'''Figure 1''': Illustration of an Ising lattice in one (N=5), two (N=5×5), and three (N=5×5×5) dimensions. Red cells indicate the "up" spin", and blue cells indicate the "down" spin.]]

The rules of the model are simple:

* When no magnetic field is applied, the only energetic contribution comes from the interaction between spins in ''adjacent'' lattice cells
* The total interaction energy is defined as <math>- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>, where ''J'' is a constant which controls the strength of interaction and the notation ''j'' &in; neighbours(''i'') indicates that spin ''j'' must lie in a lattice cell adjacent to spin ''i''. The factor of one half is included to make sure that we do not count each spin-spin interaction twice.
* Periodic boundary conditions are applied. In other words, the "final" spin in a particular direction is considered to be adjacent to the "first" spin in that direction. In '''figure 2''', for example, cell 1 has neighbour interactions with cells 2, 3, 4, and 7, while cell 9 interacts with cells 3, 6, 7, and 8.

[[File:ThirdYearCMPExpt-IsingNumberedCell.png|300px|thumb|right|'''Figure 2''': Illustration of a 2D Ising lattice, with numbered cells to indicate the neighbour interactions.]]

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

==Phase Transitions==

The Ising model is interesting because it is one of the simplest physical models to display a phase transition (so long as we are in two dimensions or more). The temperature of the system controls the balance between the energetic driving force, which wants the system to adopt the lowest energy configuration that you imagined in the previous section, and the entropic driving force, which seeks to maximise the number of configurations available to the system.

'''<big>TASK 1b</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 (''D''=3, ''N''=1000)? How much entropy does the system gain by doing so?'''

At low temperatures, little thermal energy is available to overcome such energy barriers, and we expect the system to adopt the lowest energy (the energetic driving force will dominate). Under these conditions, we expect a large number of parallel spins, and we can characterise this by considering the total ''magnetisation'' of the system: <math>M = \sum_i s_i</math>. At higher temperatures, however, enough thermal energy will be available that spins can flip relatively freely. Under these conditions, the entropic driving force should dominate. The transition temperature between the two regimes is called the ''Curie temperature'', ''TC''. Below ''TC'', we observe the highly ordered magnetised phase, while above ''TC'' a disordered, demagnetised, phase emerges.

'''<big>TASK 1c</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 ''D''=3, ''N''=1000 at absolute zero?'''

We are interested, therefore, in determining the magnetisation ''M'', and energy ''E'', of the system as functions of temperature ''T''. We can write statistical mechanical equations for these quantities:

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

where ''α'' is used as a shorthand to represent all spins in the system, and ''Z'' is the partition function,
<math>Z = \sum_\alpha \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>.

It was discovered relatively quickly after the model was first proposed that the partition function in the one dimensional case can be solved analytically (though doing so is beyond the scope of this experiment):

<math>Z\left(T, N\right) = \left[2\cosh\left(\frac{J}{k_BT}\right)\right]^N</math>

Much later, [http://en.wikipedia.org/wiki/Lars_Onsager Lars Onsager] established that the two dimensional partition function can also be found analytically, though the mathematics is significantly more involved and the resulting expression is correspondingly much more complicated. In dimensions greater than two, no analytical solution is known.

Luckily for us, the problem can instead be tackled by numerical methods, and it is this that will be the focus of the experiment. We will discuss how the equations for ''M'' and ''E'' can be solved, but first we must set up a Python script to model the lattice.

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]].</big>'''

Programming a 2D Ising Model/Introduction to Monte Carlo simulation

2026-05-03T11:34:56Z

Gmallia: /* Modifying IsingLattice.py */

'''<big>This is the third section of the Ising model experiment. You can return to the previous page, [[Programming a 2D Ising Model/Calculating the energy and magnetisation|Calculating the energy and magnetisation]], or jump ahead to the next section, [[Programming a 2D Ising Model/Accelerating the code|Accelerating the code]].</big>'''
==Average Energy and Magnetisation==

Consider again the expressions for the average energy and magnetisation that we gave in the introduction.

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

Imagine we want to evaluate these at a particular temperature, in a system of 100 spins.

'''<big>TASK 3a</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 1×109 configurations per second with our computer. How long will it take to evaluate a single value of &lang;''M''&rang;''T''?'''

Clearly, we need to try a cleverer approach.

For the overwhelming majority of the possible configurations, the Boltzmann factor, <math>\exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\}</math>, will be very small, and that state will contribute very little to the average value. We can save an enormous amount of time, and make this problem tractable, if we only consider the states whose Boltzmann factor is not so vanishingly small. This technique is called ''importance sampling'' — instead of sampling every point in phase space, we sample only those that the system is likely to occupy.

==Importance Sampling==

This is easily stated, of course, but how can we know which states contribute a lot to the average without actually calculating the energy? So far, we have been imagining generating random states of the system; that is to say, we have been choosing ''α'' from a uniform distribution in which no arrangement of spins is any more likely to be chosen than any other. If we instead had a method which generated states randomly from the probability distribution <math>\exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\}</math>, then our problem would be solved. The method which allows us to do this was developed in the 1950s, and is called the [https://doi.org/10.1063/1.1699114 Metropolis Monte Carlo (MMC)] method, or more often simply "Monte Carlo simulation". This was one of the major breakthroughs of the early days of computational science, and in the year 2000, the IEEE named among the [https://doi.org/10.1109/MCISE.2000.814652 "Top 10 Algorithms of the 20th Century"].

The algorithm is as follows:

# Start from a given configuration of spins, ''α0'', with energy ''E0''.
# Choose a single spin '''at random''', and "flip" it, to generate a new configuration ''α1''.
# Calculate the energy of this new configuration, ''E1''.
# Calculate the energy difference between the states, ''ΔE = E1 - E0''.
## If the ''ΔE ≤ 0'' (the spin flipping decreased the energy), then we '''accept''' the new configuration.
##* We set ''α0 = α1'', and ''E0 = E1'', and then '''go to step 5'''.
## If ''ΔE > 0'', the spin flipping increased the energy. By considering the probability of observing the starting and final states, ''α0'' and ''α1'', it can be shown that the probability for the transition between the two to occur is <math>\exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>. To ensure that we only accept this kind of spin flip with the correct probability, we use the following procedure:
### Choose a random number, ''R'', in the interval [0,1).
### If <math>R \leq \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''accept''' the new configuration.
###* We set ''α0 = α1'', and ''E0 = E1'', and then '''go to step 5'''.
### If <math>R > \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''reject''' the new configuration.
###* ''α0'' and ''E0'' are left unchanged. '''Go to step 5'''.
# Update the running averages of the energy and magnetisation.
# Monte Carlo "step" complete, '''return to step 2'''.

Step 4 is the key to the method. By accepting moves which increase the energy only with a certain probability, we ensure that the sequence of states that we generate is correctly distributed. A transition which carries a very large energy penalty, ''ΔE'' is extremely unlikely to be selected. The use of random numbers in this step is the reason that the method acquired the name "Monte Carlo", after the casinos located there!

Note, unlike molecular dynamics, each step only moves a single atom. In Monte Carlo simulation, there is the notion of a cycle which is ''N''spins-steps, so that each spin on average has an attempted flip per cycle.

If you are interested in the mathematical details of why this procedure generates a sequence of states distributed in the correct way, consult the Monte Carlo chapter of [https://library-search.imperial.ac.uk/permalink/f/o1297h/TN_cdi_askewsholts_vlebooks_9780080519982 "Understanding Molecular Simulation"], by Frenkel and Smit, but a discussion of this is not required for the experiment. Chapter 7 of [https://library-search.imperial.ac.uk/permalink/44IMP_INST/mek6kh/alma997227534401591 "Statistical Mechanics: Theory and Molecular Simulation"] by Tuckerman also gives a good description.

==Modifying IsingLattice.py==

Our IsingLattice object contains a stub function, '''<code>montecarlostep(T)</code>''', which takes a single argument — the temperature at which we want to do the cycle. There are also attributes used to calculate the average energy and magnetisation and respective variances: <code>self.E_tally</code>, <code>self.E2_tally</code>, <code>self.M_tally</code>, <code>self.M2_tally</code>, and <code>self.n_steps</code>. You should use these values to keep a running sum of all the energies, squared energies, magnetisations, and squared magnetisations that your system has sampled, as well as the number of steps that have been performed. Finally, there is also a function <code>statistics()</code>, that takes no arguments and should return the average energy, squared energy, magnetisation, squared magnetisation and current step, '''in that order'''.

'''<big>TASK 3b</big>: Implement a single step of the above algorithm in the <code>montecarlostep(T)</code> function. This function should return the energy of your lattice and the magnetisation at the end of the cycle. Complete the <code>statistics()</code> function. This should return the following quantities whenever it is called: &lang;''E''&rang;, &lang;''E''2&rang;, &lang;''M''&rang;, &lang;''M''2&rang;, and ''N''steps, the number of Monte Carlo steps that have elapsed. Put your code for these functions in the report.'''

''Note on units:'' In our simulation we will be using an unitless scaled energy <math>\tilde{E}=\frac{E}{J}</math>, and a unitless scaled temperature <math>\tilde{T}=\frac{k_B T}{J}</math>, where ''J'' is the interaction energy between spins discussed previously. The use of these scaled quantities in practice means that in the code we can set ''kB'' to numerically equal to 1. Going forward we will drop the &tilde; from the notation of energy and temperature.

You have been provided with a script, '''ILanim.py''', which will use your <code>IsingLattice</code> object to run a Monte Carlo simulation and display the output on the screen in real time. By default, this simulation will run an 8×8 lattice at ''T''=1.0, which is below the critical temperature. You will see a representation of the lattice, and graphs of the energy and magnetisation '''per spin'''. When you close the window, the script will use your '''<code>statistics()</code>''' function to print the average energy and magnetisation from the simulation. Note, you should run this script from the command line.

'''<big>TASK 3c</big>: If ''T'' < ''TC'', do you expect a spontaneous magnetisation (i.e. do you expect &lang;''M''&rang; ≠ 0)? When the state of the simulation appears to stop changing (when you have reached an equilibrium state), use the controls to export the output to an SVG or PNG file and attach this to your report. You should also include the output from your statistics() function.'''

'''<big>This is the third section of the Ising model experiment. You can return to the previous page, [[Programming a 2D Ising Model/Calculating the energy and magnetisation|Calculating the energy and magnetisation]], or jump ahead to the</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Accelerating the code|4. Accelerating the code]].</big>'''

Programming a 2D Ising Model/Introduction to the Ising model

2026-05-03T11:33:53Z

Gmallia: /* Phase Transitions */

'''<big>This is the first section of the Ising model experiment. You can return to the [[Programming_a_2D_Ising_Model|introduction page]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]].</big>'''

= 1. Introduction to the Ising model =

== The Ising model ==
The Ising model is a simple physical model proposed in the 1920s to understand the behaviour of ferromagnets. In most materials, the magnetic dipoles of the atoms have no preferred orientation, and the bulk material has no net magnetic moment. In certain materials, however, there is preferential alignment of the magnetic dipoles along certain directions. These materials, such as iron, exhibit an overall magnetic moment, and are called ferromagnetic.

The phenomenon is based on competition between two physical effects — energy minimisation and entropy maximisation. On the one hand, parallel alignment of magnetic moments is energetically favourable. On the other hand, aligning the system in this way reduces the entropy.

In the Ising model, we consider a collection of spins, ''si'', each of which can either have the value "up" (+1) or "down" (-1), arranged on a lattice. In one dimension, this takes the form of a "chain" of adjacent spins, ''s1'' to ''sN''. In two dimensions we have a grid of cells, each of which contains one spin. In three dimensions, the system is a cuboid, and so on. The arrangements for 1, 2 and 3 dimensions are illustrated in '''figure 1'''.

[[File:ThirdYearCMPExpt-IsingSketch.png|300px|thumb|right|'''Figure 1''': Illustration of an Ising lattice in one (N=5), two (N=5×5), and three (N=5×5×5) dimensions. Red cells indicate the "up" spin", and blue cells indicate the "down" spin.]]

The rules of the model are simple:

* When no magnetic field is applied, the only energetic contribution comes from the interaction between spins in ''adjacent'' lattice cells
* The total interaction energy is defined as <math>- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>, where ''J'' is a constant which controls the strength of interaction and the notation ''j'' &in; neighbours(''i'') indicates that spin ''j'' must lie in a lattice cell adjacent to spin ''i''. The factor of one half is included to make sure that we do not count each spin-spin interaction twice.
* Periodic boundary conditions are applied. In other words, the "final" spin in a particular direction is considered to be adjacent to the "first" spin in that direction. In '''figure 2''', for example, cell 1 has neighbour interactions with cells 2, 3, 4, and 7, while cell 9 interacts with cells 3, 6, 7, and 8.

[[File:ThirdYearCMPExpt-IsingNumberedCell.png|300px|thumb|right|'''Figure 2''': Illustration of a 2D Ising lattice, with numbered cells to indicate the neighbour interactions.]]

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

==Phase Transitions==

The Ising model is interesting because it is one of the simplest physical models to display a phase transition (so long as we are in two dimensions or more). The temperature of the system controls the balance between the energetic driving force, which wants the system to adopt the lowest energy configuration that you imagined in the previous section, and the entropic driving force, which seeks to maximise the number of configurations available to the system.

'''<big>TASK 1b</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 (''D''=3, ''N''=1000)? How much entropy does the system gain by doing so?'''

At low temperatures, little thermal energy is available to overcome such energy barriers, and we expect the system to adopt the lowest energy (the energetic driving force will dominate). Under these conditions, we expect a large number of parallel spins, and we can characterise this by considering the total ''magnetisation'' of the system: <math>M = \sum_i s_i</math>. At higher temperatures, however, enough thermal energy will be available that spins can flip relatively freely. Under these conditions, the entropic driving force should dominate. The transition temperature between the two regimes is called the ''Curie temperature'', ''TC''. Below ''TC'', we observe the highly ordered magnetised phase, while above ''TC'' a disordered, demagnetised, phase emerges.

'''<big>TASK 1c</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 ''D''=3, ''N''=1000 at absolute zero?'''

We are interested, therefore, in determining the magnetisation ''M'', and energy ''E'', of the system as functions of temperature ''T''. We can write statistical mechanical equations for these quantities:

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

where ''α'' is used as a shorthand to represent all spins in the system, and ''Z'' is the partition function,
<math>Z = \sum_\alpha \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>.

It was discovered relatively quickly after the model was first proposed that the partition function in the one dimensional case can be solved analytically (though doing so is beyond the scope of this experiment):

<math>Z\left(T, N\right) = \left[2\cosh\left(\frac{J}{k_BT}\right)\right]^N</math>

Much later, [http://en.wikipedia.org/wiki/Lars_Onsager Lars Onsager] established that the two dimensional partition function can also be found analytically, though the mathematics is significantly more involved and the resulting expression is correspondingly much more complicated. In dimensions greater than two, no analytical solution is known.

Luckily for us, the problem can instead be tackled by numerical methods, and it is this that will be the focus of the experiment. We will discuss how the equations for ''M'' and ''E'' can be solved, but first we must set up a Python script to model the lattice.

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Main page: [[Programming a 2D Ising Model]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]].</big>'''

Programming a 2D Ising Model/Calculating the energy and magnetisation

2026-05-03T11:32:00Z

Gmallia: /* Testing the files */

'''<big>This is the second section of the Ising model experiment. You can return to the previous page, [[Programming a 2D Ising Model/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

= 2. Calculating the energy and magnetisation =

==Getting the files for the experiment==

'''Various scripts have been prepared to assist you with this experiment, and you can obtain them by downloading the files from the Blackboard.'''

We recommend you open the extracted folder with JupyterLab, which has an integrated terminal where you can run the simulation scripts.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. This class will be the basis for our model of the system, and is restricted to modelling the two dimensional case. It contains a number of stub methods that you will complete as you work through the experiment.

===A note about Python Classes===

The file IsingLattice.py makes use of a Python feature called ''class''. You should be familiar with the idea that all Python ''objects'' have an associated ''type'' -- you have already encountered objects such as integers, string, lists, and NumPy arrays. Classes allow us to define new types of object. In this case, we define a type called <code>IsingLattice</code>:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. We can see, for example, that each <code>IsingLattice</code> object has the methods <code>energy()</code>, <code>magnetisation()</code>, <code>montecarlostep(T)</code>, and <code>statistics()</code> (the method <code>__init__</code> will be explained shortly). Each of these methods takes a special first argument called self, which is a variable pointing to the current copy of IsingLattice. To use the new IsingLattice object that we have defined, we use code like:

<pre>
import IsingLattice

il = IsingLattice(5, 5) #create an IsingLattice object with 5 row and 5 columns

energy = il.energy() #call the energy() method of our IsingLattice object
</pre>

When writing code "within" the class definition, all attributes and methods belonging to the class must be accessed using the "self." notation. You can see an example of this in the <code>montecarlostep()</code> method:

<pre>
energy = self.energy()
</pre>

This line creates a variable called energy which contains whatever value was returned by the method <code>energy()</code> which is part of the <code>IsingLattice</code> class. If you just want to define a local variable instead, like an iterator in a loop, the <code>self.</code> notation is not required.

<pre>
def magnetisation(self):
# we have to use self.lattice, because the lattice variable is part of the object, not local to this function
for i in self.lattice:
...
# we use i, not self.i
# we don't need to access i outside this function
</pre>

Have a look at the other files that you have been given to see how this IsingLattice object is used.

===The Constructor===

The only method currently filled in is <code>__init__</code>. This is a special method called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this method is called automatically. It takes two arguments (excluding "self"), which are simply the number of rows and columns that our lattice will have. The attribute <code>self.lattice</code> is then created, which is a NumPy array with the requested number of rows and columns. The function <code>np.random.choice</code> is used to create a random collection of spins for our initial configuration. You can see the documentation for the function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html numpy.random.choice].

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

==Testing the files==

When you have completed the energy() and magnetisation() methods, we need to test that they work correctly.
We have provided you with a script, '''ILcheck.py''', which will create three IsingLattice objects and check their energies and magnetisations, displaying these on a graph. One configuration corresponds to the energy minimum, one to the energy maximum, and one to an intermediate state.

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

There are also some unit tests available in the <code>test_energy.py</code> file. These tests are designed to be run by simply running the [https://docs.pytest.org pytest] command in your terminal,
<pre>pytest</pre>
You should make sure your code passes the tests for the <code>energy()</code> and <code>magnetisation()</code> methods at this time. The file contains tests for methods that will only be implemented later on, so some failures are expected for now.
These sorts of test are useful because they are easier to write and faster to run and can be made to run automatically.

You can also debug and test your program on a Jupyter notebook or IPython console by importing your IsingLattice class, creating one object, and inspecting its properties. If you do this, we recommend that you start your session running the commands:

<pre>
%load_ext autoreload
%autoreload 2
</pre>

This will ensure that when you try executing your code, you always run the most recent version (please note this isn't fool proof and if you find problems you may need to restart the kernel)!

'''<big>This is the second section of the Ising model experiment. You can return to the previous page, [[Programming a 2D Ising Model/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|3. Introduction to Monte Carlo simulation]].</big>'''

Programming a 2D Ising Model/Calculating the energy and magnetisation

2026-05-03T11:31:42Z

Gmallia: /* Testing the files */

'''<big>This is the second section of the Ising model experiment. You can return to the previous page, [[Programming a 2D Ising Model/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

= 2. Calculating the energy and magnetisation =

==Getting the files for the experiment==

'''Various scripts have been prepared to assist you with this experiment, and you can obtain them by downloading the files from the Blackboard.'''

We recommend you open the extracted folder with JupyterLab, which has an integrated terminal where you can run the simulation scripts.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. This class will be the basis for our model of the system, and is restricted to modelling the two dimensional case. It contains a number of stub methods that you will complete as you work through the experiment.

===A note about Python Classes===

The file IsingLattice.py makes use of a Python feature called ''class''. You should be familiar with the idea that all Python ''objects'' have an associated ''type'' -- you have already encountered objects such as integers, string, lists, and NumPy arrays. Classes allow us to define new types of object. In this case, we define a type called <code>IsingLattice</code>:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. We can see, for example, that each <code>IsingLattice</code> object has the methods <code>energy()</code>, <code>magnetisation()</code>, <code>montecarlostep(T)</code>, and <code>statistics()</code> (the method <code>__init__</code> will be explained shortly). Each of these methods takes a special first argument called self, which is a variable pointing to the current copy of IsingLattice. To use the new IsingLattice object that we have defined, we use code like:

<pre>
import IsingLattice

il = IsingLattice(5, 5) #create an IsingLattice object with 5 row and 5 columns

energy = il.energy() #call the energy() method of our IsingLattice object
</pre>

When writing code "within" the class definition, all attributes and methods belonging to the class must be accessed using the "self." notation. You can see an example of this in the <code>montecarlostep()</code> method:

<pre>
energy = self.energy()
</pre>

This line creates a variable called energy which contains whatever value was returned by the method <code>energy()</code> which is part of the <code>IsingLattice</code> class. If you just want to define a local variable instead, like an iterator in a loop, the <code>self.</code> notation is not required.

<pre>
def magnetisation(self):
# we have to use self.lattice, because the lattice variable is part of the object, not local to this function
for i in self.lattice:
...
# we use i, not self.i
# we don't need to access i outside this function
</pre>

Have a look at the other files that you have been given to see how this IsingLattice object is used.

===The Constructor===

The only method currently filled in is <code>__init__</code>. This is a special method called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this method is called automatically. It takes two arguments (excluding "self"), which are simply the number of rows and columns that our lattice will have. The attribute <code>self.lattice</code> is then created, which is a NumPy array with the requested number of rows and columns. The function <code>np.random.choice</code> is used to create a random collection of spins for our initial configuration. You can see the documentation for the function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html numpy.random.choice].

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

==Testing the files==

When you have completed the energy() and magnetisation() methods, we need to test that they work correctly.
We have provided you with a script, '''ILcheck.py''', which will create three IsingLattice objects and check their energies and magnetisations, displaying these on a graph. One configuration corresponds to the energy minimum, one to the energy maximum, and one to an intermediate state.

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

There are also some unit tests available in the <code>test_energy.py</code> file. These tests are designed to be run by simply running the [https://docs.pytest.org pytest] command in your terminal,
<pre>pytest</pre>
You should make sure your code passes the tests for the <code>energy()</code> and <code>magnetisation()</code> methods at this time. The file contains tests for methods that will only be implemented later on, so some failures are expected for now.
These sorts of test are useful because they are easier to write and faster to run and can be made to run automatically.

You can also debug and test your program on a Jupyter notebook or IPython console by importing your IsingLattice class, creating one object, and inspecting its properties. If you do this, we recommend that you start your session running the commands:

<pre>
%load_ext autoreload
%autoreload 2
</pre>

This will ensure that when you try executing your code, you always run the most recent version (please note this isn't fool proof and if you find problems you may need to restart the kernel)!

'''<big>This is the second section of the Ising model experiment. You can return to the previous page, [[Programming a 2D Ising Model/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the</big>'''

'''<big>Next section: [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|2. Introduction to Monte Carlo simulation]].</big>'''

Programming a 2D Ising Model/Introduction to the Ising model

2026-05-03T09:31:28Z

Gmallia: /* 1. Introduction to the Ising model */

'''<big>This is the first section of the Ising model experiment. You can return to the [[Programming_a_2D_Ising_Model|introduction page]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]].</big>'''

= 1. Introduction to the Ising model =

== The Ising model ==
The Ising model is a simple physical model proposed in the 1920s to understand the behaviour of ferromagnets. In most materials, the magnetic dipoles of the atoms have no preferred orientation, and the bulk material has no net magnetic moment. In certain materials, however, there is preferential alignment of the magnetic dipoles along certain directions. These materials, such as iron, exhibit an overall magnetic moment, and are called ferromagnetic.

The phenomenon is based on competition between two physical effects — energy minimisation and entropy maximisation. On the one hand, parallel alignment of magnetic moments is energetically favourable. On the other hand, aligning the system in this way reduces the entropy.

In the Ising model, we consider a collection of spins, ''si'', each of which can either have the value "up" (+1) or "down" (-1), arranged on a lattice. In one dimension, this takes the form of a "chain" of adjacent spins, ''s1'' to ''sN''. In two dimensions we have a grid of cells, each of which contains one spin. In three dimensions, the system is a cuboid, and so on. The arrangements for 1, 2 and 3 dimensions are illustrated in '''figure 1'''.

[[File:ThirdYearCMPExpt-IsingSketch.png|300px|thumb|right|'''Figure 1''': Illustration of an Ising lattice in one (N=5), two (N=5×5), and three (N=5×5×5) dimensions. Red cells indicate the "up" spin", and blue cells indicate the "down" spin.]]

The rules of the model are simple:

* When no magnetic field is applied, the only energetic contribution comes from the interaction between spins in ''adjacent'' lattice cells
* The total interaction energy is defined as <math>- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>, where ''J'' is a constant which controls the strength of interaction and the notation ''j'' &in; neighbours(''i'') indicates that spin ''j'' must lie in a lattice cell adjacent to spin ''i''. The factor of one half is included to make sure that we do not count each spin-spin interaction twice.
* Periodic boundary conditions are applied. In other words, the "final" spin in a particular direction is considered to be adjacent to the "first" spin in that direction. In '''figure 2''', for example, cell 1 has neighbour interactions with cells 2, 3, 4, and 7, while cell 9 interacts with cells 3, 6, 7, and 8.

[[File:ThirdYearCMPExpt-IsingNumberedCell.png|300px|thumb|right|'''Figure 2''': Illustration of a 2D Ising lattice, with numbered cells to indicate the neighbour interactions.]]

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

==Phase Transitions==

The Ising model is interesting because it is one of the simplest physical models to display a phase transition (so long as we are in two dimensions or more). The temperature of the system controls the balance between the energetic driving force, which wants the system to adopt the lowest energy configuration that you imagined in the previous section, and the entropic driving force, which seeks to maximise the number of configurations available to the system.

'''<big>TASK 1b</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 (''D''=3, ''N''=1000)? How much entropy does the system gain by doing so?'''

At low temperatures, little thermal energy is available to overcome such energy barriers, and we expect the system to adopt the lowest energy (the energetic driving force will dominate). Under these conditions, we expect a large number of parallel spins, and we can characterise this by considering the total ''magnetisation'' of the system: <math>M = \sum_i s_i</math>. At higher temperatures, however, enough thermal energy will be available that spins can flip relatively freely. Under these conditions, the entropic driving force should dominate. The transition temperature between the two regimes is called the ''Curie temperature'', ''TC''. Below ''TC'', we observe the highly ordered magnetised phase, while above ''TC'' a disordered, demagnetised, phase emerges.

'''<big>TASK 1c</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 ''D''=3, ''N''=1000 at absolute zero?'''

We are interested, therefore, in determining the magnetisation ''M'', and energy ''E'', of the system as functions of temperature ''T''. We can write statistical mechanical equations for these quantities:

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

where ''α'' is used as a shorthand to represent all spins in the system, and ''Z'' is the partition function,
<math>Z = \sum_\alpha \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>.

It was discovered relatively quickly after the model was first proposed that the partition function in the one dimensional case can be solved analytically (though doing so is beyond the scope of this experiment):

<math>Z\left(T, N\right) = \left[2\cosh\left(\frac{J}{k_BT}\right)\right]^N</math>

Much later, [http://en.wikipedia.org/wiki/Lars_Onsager Lars Onsager] established that the two dimensional partition function can also be found analytically, though the mathematics is significantly more involved and the resulting expression is correspondingly much more complicated. In dimensions greater than two, no analytical solution is known.

Luckily for us, the problem can instead be tackled by numerical methods, and it is this that will be the focus of the experiment. We will discuss how the equations for ''M'' and ''E'' can be solved, but first we must set up a Python script to model the lattice.

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Previous section: [[Programming_a_2D_Ising_Model|Introduction page]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]].</big>'''

Programming a 2D Ising Model/Calculating the energy and magnetisation

2026-05-03T09:30:18Z

Gmallia: /* Getting the files for the experiment */

'''<big>This is the second section of the Ising model experiment. You can return to the previous page, [[Programming a 2D Ising Model/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

= 2. Calculating the energy and magnetisation =

==Getting the files for the experiment==

'''Various scripts have been prepared to assist you with this experiment, and you can obtain them by downloading the files from the Blackboard.'''

We recommend you open the extracted folder with JupyterLab, which has an integrated terminal where you can run the simulation scripts.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. This class will be the basis for our model of the system, and is restricted to modelling the two dimensional case. It contains a number of stub methods that you will complete as you work through the experiment.

===A note about Python Classes===

The file IsingLattice.py makes use of a Python feature called ''class''. You should be familiar with the idea that all Python ''objects'' have an associated ''type'' -- you have already encountered objects such as integers, string, lists, and NumPy arrays. Classes allow us to define new types of object. In this case, we define a type called <code>IsingLattice</code>:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. We can see, for example, that each <code>IsingLattice</code> object has the methods <code>energy()</code>, <code>magnetisation()</code>, <code>montecarlostep(T)</code>, and <code>statistics()</code> (the method <code>__init__</code> will be explained shortly). Each of these methods takes a special first argument called self, which is a variable pointing to the current copy of IsingLattice. To use the new IsingLattice object that we have defined, we use code like:

<pre>
import IsingLattice

il = IsingLattice(5, 5) #create an IsingLattice object with 5 row and 5 columns

energy = il.energy() #call the energy() method of our IsingLattice object
</pre>

When writing code "within" the class definition, all attributes and methods belonging to the class must be accessed using the "self." notation. You can see an example of this in the <code>montecarlostep()</code> method:

<pre>
energy = self.energy()
</pre>

This line creates a variable called energy which contains whatever value was returned by the method <code>energy()</code> which is part of the <code>IsingLattice</code> class. If you just want to define a local variable instead, like an iterator in a loop, the <code>self.</code> notation is not required.

<pre>
def magnetisation(self):
# we have to use self.lattice, because the lattice variable is part of the object, not local to this function
for i in self.lattice:
...
# we use i, not self.i
# we don't need to access i outside this function
</pre>

Have a look at the other files that you have been given to see how this IsingLattice object is used.

===The Constructor===

The only method currently filled in is <code>__init__</code>. This is a special method called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this method is called automatically. It takes two arguments (excluding "self"), which are simply the number of rows and columns that our lattice will have. The attribute <code>self.lattice</code> is then created, which is a NumPy array with the requested number of rows and columns. The function <code>np.random.choice</code> is used to create a random collection of spins for our initial configuration. You can see the documentation for the function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html numpy.random.choice].

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

==Testing the files==

When you have completed the energy() and magnetisation() methods, we need to test that they work correctly.
We have provided you with a script, '''ILcheck.py''', which will create three IsingLattice objects and check their energies and magnetisations, displaying these on a graph. One configuration corresponds to the energy minimum, one to the energy maximum, and one to an intermediate state.

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

There are also some unit tests available in the <code>test_energy.py</code> file. These tests are designed to be run by simply running the [https://docs.pytest.org pytest] command in your terminal,
<pre>pytest</pre>
You should make sure your code passes the tests for the <code>energy()</code> and <code>magnetisation()</code> methods at this time. The file contains tests for methods that will only be implemented later on, so some failures are expected for now.
These sorts of test are useful because they are easier to write and faster to run and can be made to run automatically.

You can also debug and test your program on a Jupyter notebook or IPython console by importing your IsingLattice class, creating one object, and inspecting its properties. If you do this, we recommend that you start your session running the commands:

<pre>
%load_ext autoreload
%autoreload 2
</pre>

This will ensure that when you try executing your code, you always run the most recent version (please note this isn't fool proof and if you find problems you may need to restart the kernel)!

'''<big>This is the second section of the Ising model experiment. You can return to the previous page, [[Programming a 2D Ising Model/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

Programming a 2D Ising Model/Introduction to the Ising model

2026-05-03T09:27:20Z

Gmallia: /* 1. The Ising model */

'''<big>This is the first section of the Ising model experiment. You can return to the [[Programming_a_2D_Ising_Model|introduction page]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]].</big>'''

= 1. Introduction to the Ising model =
The Ising model is a simple physical model proposed in the 1920s to understand the behaviour of ferromagnets. In most materials, the magnetic dipoles of the atoms have no preferred orientation, and the bulk material has no net magnetic moment. In certain materials, however, there is preferential alignment of the magnetic dipoles along certain directions. These materials, such as iron, exhibit an overall magnetic moment, and are called ferromagnetic.

The phenomenon is based on competition between two physical effects — energy minimisation and entropy maximisation. On the one hand, parallel alignment of magnetic moments is energetically favourable. On the other hand, aligning the system in this way reduces the entropy.

In the Ising model, we consider a collection of spins, ''si'', each of which can either have the value "up" (+1) or "down" (-1), arranged on a lattice. In one dimension, this takes the form of a "chain" of adjacent spins, ''s1'' to ''sN''. In two dimensions we have a grid of cells, each of which contains one spin. In three dimensions, the system is a cuboid, and so on. The arrangements for 1, 2 and 3 dimensions are illustrated in '''figure 1'''.

[[File:ThirdYearCMPExpt-IsingSketch.png|300px|thumb|right|'''Figure 1''': Illustration of an Ising lattice in one (N=5), two (N=5×5), and three (N=5×5×5) dimensions. Red cells indicate the "up" spin", and blue cells indicate the "down" spin.]]

The rules of the model are simple:

* When no magnetic field is applied, the only energetic contribution comes from the interaction between spins in ''adjacent'' lattice cells
* The total interaction energy is defined as <math>- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>, where ''J'' is a constant which controls the strength of interaction and the notation ''j'' &in; neighbours(''i'') indicates that spin ''j'' must lie in a lattice cell adjacent to spin ''i''. The factor of one half is included to make sure that we do not count each spin-spin interaction twice.
* Periodic boundary conditions are applied. In other words, the "final" spin in a particular direction is considered to be adjacent to the "first" spin in that direction. In '''figure 2''', for example, cell 1 has neighbour interactions with cells 2, 3, 4, and 7, while cell 9 interacts with cells 3, 6, 7, and 8.

[[File:ThirdYearCMPExpt-IsingNumberedCell.png|300px|thumb|right|'''Figure 2''': Illustration of a 2D Ising lattice, with numbered cells to indicate the neighbour interactions.]]

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

==Phase Transitions==

The Ising model is interesting because it is one of the simplest physical models to display a phase transition (so long as we are in two dimensions or more). The temperature of the system controls the balance between the energetic driving force, which wants the system to adopt the lowest energy configuration that you imagined in the previous section, and the entropic driving force, which seeks to maximise the number of configurations available to the system.

'''<big>TASK 1b</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 (''D''=3, ''N''=1000)? How much entropy does the system gain by doing so?'''

At low temperatures, little thermal energy is available to overcome such energy barriers, and we expect the system to adopt the lowest energy (the energetic driving force will dominate). Under these conditions, we expect a large number of parallel spins, and we can characterise this by considering the total ''magnetisation'' of the system: <math>M = \sum_i s_i</math>. At higher temperatures, however, enough thermal energy will be available that spins can flip relatively freely. Under these conditions, the entropic driving force should dominate. The transition temperature between the two regimes is called the ''Curie temperature'', ''TC''. Below ''TC'', we observe the highly ordered magnetised phase, while above ''TC'' a disordered, demagnetised, phase emerges.

'''<big>TASK 1c</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 ''D''=3, ''N''=1000 at absolute zero?'''

We are interested, therefore, in determining the magnetisation ''M'', and energy ''E'', of the system as functions of temperature ''T''. We can write statistical mechanical equations for these quantities:

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

where ''α'' is used as a shorthand to represent all spins in the system, and ''Z'' is the partition function,
<math>Z = \sum_\alpha \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>.

It was discovered relatively quickly after the model was first proposed that the partition function in the one dimensional case can be solved analytically (though doing so is beyond the scope of this experiment):

<math>Z\left(T, N\right) = \left[2\cosh\left(\frac{J}{k_BT}\right)\right]^N</math>

Much later, [http://en.wikipedia.org/wiki/Lars_Onsager Lars Onsager] established that the two dimensional partition function can also be found analytically, though the mathematics is significantly more involved and the resulting expression is correspondingly much more complicated. In dimensions greater than two, no analytical solution is known.

Luckily for us, the problem can instead be tackled by numerical methods, and it is this that will be the focus of the experiment. We will discuss how the equations for ''M'' and ''E'' can be solved, but first we must set up a Python script to model the lattice.

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Previous section: [[Programming_a_2D_Ising_Model|Introduction page]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]].</big>'''

Programming a 2D Ising Model

2026-05-03T09:27:04Z

Gmallia: /* Structure of this Experiment */

'''<big>This is the Ising Model experiment (aka programming for simple simulations). This is a core experiment for year 3 students on the Chemistry with Molecular Physics degree.</big>'''

==Introduction==

In this exercise, you are going to use the Python that you learned in the first year to write a code to perform Monte-Carlo simulations of the 2D [http://en.wikipedia.org/wiki/Ising_model Ising model], a set of spins on a lattice which is used to model ferromagnetic behaviour, and also to analyse the results of the simulation to find the heat capacity of the system and the Curie temperature — the temperature below which the system is able to maintain a spontaneous magnetisation.

A more extensive introduction to the Ising Model will be given in the next page. This page gives some practical information about how the experiment.

==Assessment==
At the end of this experiment you must submit a report, writing up your findings. Each section of the experiment has a number of tasks that you should complete, labelled '''<big>TASK</big>'''. If this is a mathematical exercise, your report should contain a short summary of the solution. If it is a graphical exercise, you should include the relevant image. Your report should explain briefly what you did in each stage of the experiment, and what your findings were.

The report should be submitted as a PDF file to Turnitin, created with your favourite word processor or LaTeX, and should contain the annotated source code of the Python scripts that your write during the experiment. For clarity, code included in the report should be written in a monospaced font (e.g. Courier), which in LaTeX can be easily achieved by placing your code inside a <code>\texttt{}</code> command. '''The code, except where relevant to discussion, can be placed in an appendix.'''

As well as including relevant source code in your report, you should also submit a bundle with all your scripts and analysis code as a separate Blackboard submission.

===On the use of AI tools===
The code you write for this experiment is assessed, and as per the [https://bb.imperial.ac.uk/bbcswebdav/pid-3344481-dt-content-rid-21115105_1/xid-21115105_1 department's guidelines on the use of AI], it is not acceptable to use generative AI tools to write computer code for this assignment.

Writing computer programs yourself enhances your logical thinking and problem solving skills. Using generative AI tools to write your code makes you dumb. Don't choose to be dumb, don't use AI to write code for you.

==Getting Help==

The Graduate Teaching Assistants demonstrating this experiment are Fionn Carman and Lei Li. Demonstrators will be available in room 232A at the following times:

{| class="wikitable"
|+ Demonstrator timetable
! !! Monday !! Tuesday !! Wednesday !! Thusrsday !! Friday
|-
| Morning || 10:00 - Experiment intro 10:00-13:00 - Li || 10:00-13:00 - Fionn|| || 10:00-12:00 - João || 10:00-13:00 - Fionn
|-
| Afternoon || 14:00-17:00 - Li|| 14:00-17:00 - Fionn|| || || 14:00-17:00 - Fionn
|}

If you have questions outside these times, please use the Discussion Board on Blackboard.

João and Giuseppe are the members of academic staff responsible for this experiment. You are welcome to address your questions and comments to them.

==Structure of this Experiment==

This experimental manual has been broken up into a number of subsections. Direct links to each of them may be found below. You should attempt them in order, and you should complete all of them to finish the experiment.

# [[Programming_a_2D_Ising_Model/Introduction_to_the_Ising_model|Introduction to the Ising model]]
# [[Programming_a_2D_Ising_Model/Calculating the energy and magnetisation|Calculating the energy and magnetisation]]
# [[Programming_a_2D_Ising_Model/Introduction to Monte Carlo simulation|Introduction to the Monte Carlo simulation]]
# [[Programming_a_2D_Ising_Model/Accelerating the code|Accelerating the code]]
# [[Programming_a_2D_Ising_Model/The effect of temperature|The effect of temperature]]
# [[Programming_a_2D_Ising_Model/The effect of system size|The effect of system size]]
# [[Programming_a_2D_Ising_Model/Determining the heat capacity|Determining the heat capacity]]
# [[Programming_a_2D_Ising_Model/Locating the Curie temperature|Locating the Curie temperature]]

Programming a 2D Ising Model/Introduction to the Ising model

2026-05-03T09:25:35Z

Gmallia: /* Phase Transitions */

'''<big>This is the first section of the Ising model experiment. You can return to the [[Programming_a_2D_Ising_Model|introduction page]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]].</big>'''

= 1. The Ising model =
The Ising model is a simple physical model proposed in the 1920s to understand the behaviour of ferromagnets. In most materials, the magnetic dipoles of the atoms have no preferred orientation, and the bulk material has no net magnetic moment. In certain materials, however, there is preferential alignment of the magnetic dipoles along certain directions. These materials, such as iron, exhibit an overall magnetic moment, and are called ferromagnetic.

The phenomenon is based on competition between two physical effects — energy minimisation and entropy maximisation. On the one hand, parallel alignment of magnetic moments is energetically favourable. On the other hand, aligning the system in this way reduces the entropy.

In the Ising model, we consider a collection of spins, ''si'', each of which can either have the value "up" (+1) or "down" (-1), arranged on a lattice. In one dimension, this takes the form of a "chain" of adjacent spins, ''s1'' to ''sN''. In two dimensions we have a grid of cells, each of which contains one spin. In three dimensions, the system is a cuboid, and so on. The arrangements for 1, 2 and 3 dimensions are illustrated in '''figure 1'''.

[[File:ThirdYearCMPExpt-IsingSketch.png|300px|thumb|right|'''Figure 1''': Illustration of an Ising lattice in one (N=5), two (N=5×5), and three (N=5×5×5) dimensions. Red cells indicate the "up" spin", and blue cells indicate the "down" spin.]]

The rules of the model are simple:

* When no magnetic field is applied, the only energetic contribution comes from the interaction between spins in ''adjacent'' lattice cells
* The total interaction energy is defined as <math>- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>, where ''J'' is a constant which controls the strength of interaction and the notation ''j'' &in; neighbours(''i'') indicates that spin ''j'' must lie in a lattice cell adjacent to spin ''i''. The factor of one half is included to make sure that we do not count each spin-spin interaction twice.
* Periodic boundary conditions are applied. In other words, the "final" spin in a particular direction is considered to be adjacent to the "first" spin in that direction. In '''figure 2''', for example, cell 1 has neighbour interactions with cells 2, 3, 4, and 7, while cell 9 interacts with cells 3, 6, 7, and 8.

[[File:ThirdYearCMPExpt-IsingNumberedCell.png|300px|thumb|right|'''Figure 2''': Illustration of a 2D Ising lattice, with numbered cells to indicate the neighbour interactions.]]

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

==Phase Transitions==

The Ising model is interesting because it is one of the simplest physical models to display a phase transition (so long as we are in two dimensions or more). The temperature of the system controls the balance between the energetic driving force, which wants the system to adopt the lowest energy configuration that you imagined in the previous section, and the entropic driving force, which seeks to maximise the number of configurations available to the system.

'''<big>TASK 1b</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 (''D''=3, ''N''=1000)? How much entropy does the system gain by doing so?'''

At low temperatures, little thermal energy is available to overcome such energy barriers, and we expect the system to adopt the lowest energy (the energetic driving force will dominate). Under these conditions, we expect a large number of parallel spins, and we can characterise this by considering the total ''magnetisation'' of the system: <math>M = \sum_i s_i</math>. At higher temperatures, however, enough thermal energy will be available that spins can flip relatively freely. Under these conditions, the entropic driving force should dominate. The transition temperature between the two regimes is called the ''Curie temperature'', ''TC''. Below ''TC'', we observe the highly ordered magnetised phase, while above ''TC'' a disordered, demagnetised, phase emerges.

'''<big>TASK 1c</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 ''D''=3, ''N''=1000 at absolute zero?'''

We are interested, therefore, in determining the magnetisation ''M'', and energy ''E'', of the system as functions of temperature ''T''. We can write statistical mechanical equations for these quantities:

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

where ''α'' is used as a shorthand to represent all spins in the system, and ''Z'' is the partition function,
<math>Z = \sum_\alpha \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>.

It was discovered relatively quickly after the model was first proposed that the partition function in the one dimensional case can be solved analytically (though doing so is beyond the scope of this experiment):

<math>Z\left(T, N\right) = \left[2\cosh\left(\frac{J}{k_BT}\right)\right]^N</math>

Much later, [http://en.wikipedia.org/wiki/Lars_Onsager Lars Onsager] established that the two dimensional partition function can also be found analytically, though the mathematics is significantly more involved and the resulting expression is correspondingly much more complicated. In dimensions greater than two, no analytical solution is known.

Luckily for us, the problem can instead be tackled by numerical methods, and it is this that will be the focus of the experiment. We will discuss how the equations for ''M'' and ''E'' can be solved, but first we must set up a Python script to model the lattice.

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Previous section: [[Programming_a_2D_Ising_Model|Introduction page]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]].</big>'''

Programming a 2D Ising Model

2026-05-03T09:23:34Z

Gmallia: /* Structure of this Experiment */

'''<big>This is the Ising Model experiment (aka programming for simple simulations). This is a core experiment for year 3 students on the Chemistry with Molecular Physics degree.</big>'''

==Introduction==

In this exercise, you are going to use the Python that you learned in the first year to write a code to perform Monte-Carlo simulations of the 2D [http://en.wikipedia.org/wiki/Ising_model Ising model], a set of spins on a lattice which is used to model ferromagnetic behaviour, and also to analyse the results of the simulation to find the heat capacity of the system and the Curie temperature — the temperature below which the system is able to maintain a spontaneous magnetisation.

A more extensive introduction to the Ising Model will be given in the next page. This page gives some practical information about how the experiment.

==Assessment==
At the end of this experiment you must submit a report, writing up your findings. Each section of the experiment has a number of tasks that you should complete, labelled '''<big>TASK</big>'''. If this is a mathematical exercise, your report should contain a short summary of the solution. If it is a graphical exercise, you should include the relevant image. Your report should explain briefly what you did in each stage of the experiment, and what your findings were.

The report should be submitted as a PDF file to Turnitin, created with your favourite word processor or LaTeX, and should contain the annotated source code of the Python scripts that your write during the experiment. For clarity, code included in the report should be written in a monospaced font (e.g. Courier), which in LaTeX can be easily achieved by placing your code inside a <code>\texttt{}</code> command. '''The code, except where relevant to discussion, can be placed in an appendix.'''

As well as including relevant source code in your report, you should also submit a bundle with all your scripts and analysis code as a separate Blackboard submission.

===On the use of AI tools===
The code you write for this experiment is assessed, and as per the [https://bb.imperial.ac.uk/bbcswebdav/pid-3344481-dt-content-rid-21115105_1/xid-21115105_1 department's guidelines on the use of AI], it is not acceptable to use generative AI tools to write computer code for this assignment.

Writing computer programs yourself enhances your logical thinking and problem solving skills. Using generative AI tools to write your code makes you dumb. Don't choose to be dumb, don't use AI to write code for you.

==Getting Help==

The Graduate Teaching Assistants demonstrating this experiment are Fionn Carman and Lei Li. Demonstrators will be available in room 232A at the following times:

{| class="wikitable"
|+ Demonstrator timetable
! !! Monday !! Tuesday !! Wednesday !! Thusrsday !! Friday
|-
| Morning || 10:00 - Experiment intro 10:00-13:00 - Li || 10:00-13:00 - Fionn|| || 10:00-12:00 - João || 10:00-13:00 - Fionn
|-
| Afternoon || 14:00-17:00 - Li|| 14:00-17:00 - Fionn|| || || 14:00-17:00 - Fionn
|}

If you have questions outside these times, please use the Discussion Board on Blackboard.

João and Giuseppe are the members of academic staff responsible for this experiment. You are welcome to address your questions and comments to them.

==Structure of this Experiment==

This experimental manual has been broken up into a number of subsections. Direct links to each of them may be found below. You should attempt them in order, and you should complete all of them to finish the experiment.

# [[Programming_a_2D_Ising_Model/Introduction_to_the_Ising_model|The Ising model]]
# [[Programming_a_2D_Ising_Model/Calculating the energy and magnetisation|Calculating the energy and magnetisation]]
# [[Programming_a_2D_Ising_Model/Introduction to Monte Carlo simulation|Introduction to the Monte Carlo simulation]]
# [[Programming_a_2D_Ising_Model/Accelerating the code|Accelerating the code]]
# [[Programming_a_2D_Ising_Model/The effect of temperature|The effect of temperature]]
# [[Programming_a_2D_Ising_Model/The effect of system size|The effect of system size]]
# [[Programming_a_2D_Ising_Model/Determining the heat capacity|Determining the heat capacity]]
# [[Programming_a_2D_Ising_Model/Locating the Curie temperature|Locating the Curie temperature]]

Programming a 2D Ising Model/Introduction to the Ising model

2026-05-03T09:23:07Z

Gmallia: /* The Ising model */

'''<big>This is the first section of the Ising model experiment. You can return to the [[Programming_a_2D_Ising_Model|introduction page]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]].</big>'''

==1. The Ising model==
The Ising model is a simple physical model proposed in the 1920s to understand the behaviour of ferromagnets. In most materials, the magnetic dipoles of the atoms have no preferred orientation, and the bulk material has no net magnetic moment. In certain materials, however, there is preferential alignment of the magnetic dipoles along certain directions. These materials, such as iron, exhibit an overall magnetic moment, and are called ferromagnetic.

The phenomenon is based on competition between two physical effects — energy minimisation and entropy maximisation. On the one hand, parallel alignment of magnetic moments is energetically favourable. On the other hand, aligning the system in this way reduces the entropy.

In the Ising model, we consider a collection of spins, ''si'', each of which can either have the value "up" (+1) or "down" (-1), arranged on a lattice. In one dimension, this takes the form of a "chain" of adjacent spins, ''s1'' to ''sN''. In two dimensions we have a grid of cells, each of which contains one spin. In three dimensions, the system is a cuboid, and so on. The arrangements for 1, 2 and 3 dimensions are illustrated in '''figure 1'''.

[[File:ThirdYearCMPExpt-IsingSketch.png|300px|thumb|right|'''Figure 1''': Illustration of an Ising lattice in one (N=5), two (N=5×5), and three (N=5×5×5) dimensions. Red cells indicate the "up" spin", and blue cells indicate the "down" spin.]]

The rules of the model are simple:

* When no magnetic field is applied, the only energetic contribution comes from the interaction between spins in ''adjacent'' lattice cells
* The total interaction energy is defined as <math>- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>, where ''J'' is a constant which controls the strength of interaction and the notation ''j'' &in; neighbours(''i'') indicates that spin ''j'' must lie in a lattice cell adjacent to spin ''i''. The factor of one half is included to make sure that we do not count each spin-spin interaction twice.
* Periodic boundary conditions are applied. In other words, the "final" spin in a particular direction is considered to be adjacent to the "first" spin in that direction. In '''figure 2''', for example, cell 1 has neighbour interactions with cells 2, 3, 4, and 7, while cell 9 interacts with cells 3, 6, 7, and 8.

[[File:ThirdYearCMPExpt-IsingNumberedCell.png|300px|thumb|right|'''Figure 2''': Illustration of a 2D Ising lattice, with numbered cells to indicate the neighbour interactions.]]

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

==Phase Transitions==

The Ising model is interesting because it is one of the simplest physical models to display a phase transition (so long as we are in two dimensions or more). The temperature of the system controls the balance between the energetic driving force, which wants the system to adopt the lowest energy configuration that you imagined in the previous section, and the entropic driving force, which seeks to maximise the number of configurations available to the system.

'''<big>TASK 1b</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 (''D''=3, ''N''=1000)? How much entropy does the system gain by doing so?'''

At low temperatures, little thermal energy is available to overcome such energy barriers, and we expect the system to adopt the lowest energy (the energetic driving force will dominate). Under these conditions, we expect a large number of parallel spins, and we can characterise this by considering the total ''magnetisation'' of the system: <math>M = \sum_i s_i</math>. At higher temperatures, however, enough thermal energy will be available that spins can flip relatively freely. Under these conditions, the entropic driving force should dominate. The transition temperature between the two regimes is called the ''Curie temperature'', ''TC''. Below ''TC'', we observe the highly ordered magnetised phase, while above ''TC'' a disordered, demagnetised, phase emerges.

'''<big>TASK 1c</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 ''D''=3, ''N''=1000 at absolute zero?'''

We are interested, therefore, in determining the magnetisation ''M'', and energy ''E'', of the system as functions of temperature ''T''. We can write statistical mechanical equations for these quantities:

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

where ''α'' is used as a shorthand to represent all spins in the system, and ''Z'' is the partition function,
<math>Z = \sum_\alpha \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>.

It was discovered relatively quickly after the model was first proposed that the partition function in the one dimensional case can be solved analytically (though doing so is beyond the scope of this experiment):

<math>Z\left(T, N\right) = \left[2\cosh\left(\frac{J}{k_BT}\right)\right]^N</math>

Much later, [http://en.wikipedia.org/wiki/Lars_Onsager Lars Onsager] established that the two dimensional partition function can also be found analytically, though the mathematics is significantly more involved and the resulting expression is correspondingly much more complicated. In dimensions greater than two, no analytical solution is known.

Luckily for us, the problem can instead be tackled by numerical methods, and it is this that will be the focus of the experiment. We will discuss how the equations for ''M'' and ''E'' can be solved, but first we must set up a Python script to model the lattice.

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Previous section: [[Programming_a_2D_Ising_Model|Introduction page]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]].</big>'''

Programming a 2D Ising Model/Introduction to the Ising model

2026-05-03T09:22:52Z

Gmallia: /* The model */

'''<big>This is the first section of the Ising model experiment. You can return to the [[Programming_a_2D_Ising_Model|introduction page]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]].</big>'''

==The Ising model==
The Ising model is a simple physical model proposed in the 1920s to understand the behaviour of ferromagnets. In most materials, the magnetic dipoles of the atoms have no preferred orientation, and the bulk material has no net magnetic moment. In certain materials, however, there is preferential alignment of the magnetic dipoles along certain directions. These materials, such as iron, exhibit an overall magnetic moment, and are called ferromagnetic.

The phenomenon is based on competition between two physical effects — energy minimisation and entropy maximisation. On the one hand, parallel alignment of magnetic moments is energetically favourable. On the other hand, aligning the system in this way reduces the entropy.

In the Ising model, we consider a collection of spins, ''si'', each of which can either have the value "up" (+1) or "down" (-1), arranged on a lattice. In one dimension, this takes the form of a "chain" of adjacent spins, ''s1'' to ''sN''. In two dimensions we have a grid of cells, each of which contains one spin. In three dimensions, the system is a cuboid, and so on. The arrangements for 1, 2 and 3 dimensions are illustrated in '''figure 1'''.

[[File:ThirdYearCMPExpt-IsingSketch.png|300px|thumb|right|'''Figure 1''': Illustration of an Ising lattice in one (N=5), two (N=5×5), and three (N=5×5×5) dimensions. Red cells indicate the "up" spin", and blue cells indicate the "down" spin.]]

The rules of the model are simple:

* When no magnetic field is applied, the only energetic contribution comes from the interaction between spins in ''adjacent'' lattice cells
* The total interaction energy is defined as <math>- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>, where ''J'' is a constant which controls the strength of interaction and the notation ''j'' &in; neighbours(''i'') indicates that spin ''j'' must lie in a lattice cell adjacent to spin ''i''. The factor of one half is included to make sure that we do not count each spin-spin interaction twice.
* Periodic boundary conditions are applied. In other words, the "final" spin in a particular direction is considered to be adjacent to the "first" spin in that direction. In '''figure 2''', for example, cell 1 has neighbour interactions with cells 2, 3, 4, and 7, while cell 9 interacts with cells 3, 6, 7, and 8.

[[File:ThirdYearCMPExpt-IsingNumberedCell.png|300px|thumb|right|'''Figure 2''': Illustration of a 2D Ising lattice, with numbered cells to indicate the neighbour interactions.]]

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

==Phase Transitions==

The Ising model is interesting because it is one of the simplest physical models to display a phase transition (so long as we are in two dimensions or more). The temperature of the system controls the balance between the energetic driving force, which wants the system to adopt the lowest energy configuration that you imagined in the previous section, and the entropic driving force, which seeks to maximise the number of configurations available to the system.

'''<big>TASK 1b</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 (''D''=3, ''N''=1000)? How much entropy does the system gain by doing so?'''

At low temperatures, little thermal energy is available to overcome such energy barriers, and we expect the system to adopt the lowest energy (the energetic driving force will dominate). Under these conditions, we expect a large number of parallel spins, and we can characterise this by considering the total ''magnetisation'' of the system: <math>M = \sum_i s_i</math>. At higher temperatures, however, enough thermal energy will be available that spins can flip relatively freely. Under these conditions, the entropic driving force should dominate. The transition temperature between the two regimes is called the ''Curie temperature'', ''TC''. Below ''TC'', we observe the highly ordered magnetised phase, while above ''TC'' a disordered, demagnetised, phase emerges.

'''<big>TASK 1c</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 ''D''=3, ''N''=1000 at absolute zero?'''

We are interested, therefore, in determining the magnetisation ''M'', and energy ''E'', of the system as functions of temperature ''T''. We can write statistical mechanical equations for these quantities:

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

where ''α'' is used as a shorthand to represent all spins in the system, and ''Z'' is the partition function,
<math>Z = \sum_\alpha \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>.

It was discovered relatively quickly after the model was first proposed that the partition function in the one dimensional case can be solved analytically (though doing so is beyond the scope of this experiment):

<math>Z\left(T, N\right) = \left[2\cosh\left(\frac{J}{k_BT}\right)\right]^N</math>

Much later, [http://en.wikipedia.org/wiki/Lars_Onsager Lars Onsager] established that the two dimensional partition function can also be found analytically, though the mathematics is significantly more involved and the resulting expression is correspondingly much more complicated. In dimensions greater than two, no analytical solution is known.

Luckily for us, the problem can instead be tackled by numerical methods, and it is this that will be the focus of the experiment. We will discuss how the equations for ''M'' and ''E'' can be solved, but first we must set up a Python script to model the lattice.

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Previous section: [[Programming_a_2D_Ising_Model|Introduction page]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]].</big>'''

Programming a 2D Ising Model/Introduction to the Ising model

2026-05-03T09:21:57Z

Gmallia: /* Phase Transitions */

'''<big>This is the first section of the Ising model experiment. You can return to the [[Programming_a_2D_Ising_Model|introduction page]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]].</big>'''

==The model==
The Ising model is a simple physical model proposed in the 1920s to understand the behaviour of ferromagnets. In most materials, the magnetic dipoles of the atoms have no preferred orientation, and the bulk material has no net magnetic moment. In certain materials, however, there is preferential alignment of the magnetic dipoles along certain directions. These materials, such as iron, exhibit an overall magnetic moment, and are called ferromagnetic.

The phenomenon is based on competition between two physical effects — energy minimisation and entropy maximisation. On the one hand, parallel alignment of magnetic moments is energetically favourable. On the other hand, aligning the system in this way reduces the entropy.

In the Ising model, we consider a collection of spins, ''si'', each of which can either have the value "up" (+1) or "down" (-1), arranged on a lattice. In one dimension, this takes the form of a "chain" of adjacent spins, ''s1'' to ''sN''. In two dimensions we have a grid of cells, each of which contains one spin. In three dimensions, the system is a cuboid, and so on. The arrangements for 1, 2 and 3 dimensions are illustrated in '''figure 1'''.

[[File:ThirdYearCMPExpt-IsingSketch.png|300px|thumb|right|'''Figure 1''': Illustration of an Ising lattice in one (N=5), two (N=5×5), and three (N=5×5×5) dimensions. Red cells indicate the "up" spin", and blue cells indicate the "down" spin.]]

The rules of the model are simple:

* When no magnetic field is applied, the only energetic contribution comes from the interaction between spins in ''adjacent'' lattice cells
* The total interaction energy is defined as <math>- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>, where ''J'' is a constant which controls the strength of interaction and the notation ''j'' &in; neighbours(''i'') indicates that spin ''j'' must lie in a lattice cell adjacent to spin ''i''. The factor of one half is included to make sure that we do not count each spin-spin interaction twice.
* Periodic boundary conditions are applied. In other words, the "final" spin in a particular direction is considered to be adjacent to the "first" spin in that direction. In '''figure 2''', for example, cell 1 has neighbour interactions with cells 2, 3, 4, and 7, while cell 9 interacts with cells 3, 6, 7, and 8.

[[File:ThirdYearCMPExpt-IsingNumberedCell.png|300px|thumb|right|'''Figure 2''': Illustration of a 2D Ising lattice, with numbered cells to indicate the neighbour interactions.]]

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

==Phase Transitions==

The Ising model is interesting because it is one of the simplest physical models to display a phase transition (so long as we are in two dimensions or more). The temperature of the system controls the balance between the energetic driving force, which wants the system to adopt the lowest energy configuration that you imagined in the previous section, and the entropic driving force, which seeks to maximise the number of configurations available to the system.

'''<big>TASK 1b</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 (''D''=3, ''N''=1000)? How much entropy does the system gain by doing so?'''

At low temperatures, little thermal energy is available to overcome such energy barriers, and we expect the system to adopt the lowest energy (the energetic driving force will dominate). Under these conditions, we expect a large number of parallel spins, and we can characterise this by considering the total ''magnetisation'' of the system: <math>M = \sum_i s_i</math>. At higher temperatures, however, enough thermal energy will be available that spins can flip relatively freely. Under these conditions, the entropic driving force should dominate. The transition temperature between the two regimes is called the ''Curie temperature'', ''TC''. Below ''TC'', we observe the highly ordered magnetised phase, while above ''TC'' a disordered, demagnetised, phase emerges.

'''<big>TASK 1c</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 ''D''=3, ''N''=1000 at absolute zero?'''

We are interested, therefore, in determining the magnetisation ''M'', and energy ''E'', of the system as functions of temperature ''T''. We can write statistical mechanical equations for these quantities:

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

where ''α'' is used as a shorthand to represent all spins in the system, and ''Z'' is the partition function,
<math>Z = \sum_\alpha \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>.

It was discovered relatively quickly after the model was first proposed that the partition function in the one dimensional case can be solved analytically (though doing so is beyond the scope of this experiment):

<math>Z\left(T, N\right) = \left[2\cosh\left(\frac{J}{k_BT}\right)\right]^N</math>

Much later, [http://en.wikipedia.org/wiki/Lars_Onsager Lars Onsager] established that the two dimensional partition function can also be found analytically, though the mathematics is significantly more involved and the resulting expression is correspondingly much more complicated. In dimensions greater than two, no analytical solution is known.

Luckily for us, the problem can instead be tackled by numerical methods, and it is this that will be the focus of the experiment. We will discuss how the equations for ''M'' and ''E'' can be solved, but first we must set up a Python script to model the lattice.

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Previous section: [[Programming_a_2D_Ising_Model|Introduction page]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|2. Calculating the energy and magnetisation]].</big>'''

Programming a 2D Ising Model/Introduction to the Ising model

2026-05-03T09:19:04Z

Gmallia: /* Phase Transitions */

'''<big>This is the first section of the Ising model experiment. You can return to the [[Programming_a_2D_Ising_Model|introduction page]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]].</big>'''

==The model==
The Ising model is a simple physical model proposed in the 1920s to understand the behaviour of ferromagnets. In most materials, the magnetic dipoles of the atoms have no preferred orientation, and the bulk material has no net magnetic moment. In certain materials, however, there is preferential alignment of the magnetic dipoles along certain directions. These materials, such as iron, exhibit an overall magnetic moment, and are called ferromagnetic.

The phenomenon is based on competition between two physical effects — energy minimisation and entropy maximisation. On the one hand, parallel alignment of magnetic moments is energetically favourable. On the other hand, aligning the system in this way reduces the entropy.

In the Ising model, we consider a collection of spins, ''si'', each of which can either have the value "up" (+1) or "down" (-1), arranged on a lattice. In one dimension, this takes the form of a "chain" of adjacent spins, ''s1'' to ''sN''. In two dimensions we have a grid of cells, each of which contains one spin. In three dimensions, the system is a cuboid, and so on. The arrangements for 1, 2 and 3 dimensions are illustrated in '''figure 1'''.

[[File:ThirdYearCMPExpt-IsingSketch.png|300px|thumb|right|'''Figure 1''': Illustration of an Ising lattice in one (N=5), two (N=5×5), and three (N=5×5×5) dimensions. Red cells indicate the "up" spin", and blue cells indicate the "down" spin.]]

The rules of the model are simple:

* When no magnetic field is applied, the only energetic contribution comes from the interaction between spins in ''adjacent'' lattice cells
* The total interaction energy is defined as <math>- \frac{1}{2} J \sum_i^N \sum_{j\ \in\ \mathrm{neighbours}\left(i\right)} s_i s_j</math>, where ''J'' is a constant which controls the strength of interaction and the notation ''j'' &in; neighbours(''i'') indicates that spin ''j'' must lie in a lattice cell adjacent to spin ''i''. The factor of one half is included to make sure that we do not count each spin-spin interaction twice.
* Periodic boundary conditions are applied. In other words, the "final" spin in a particular direction is considered to be adjacent to the "first" spin in that direction. In '''figure 2''', for example, cell 1 has neighbour interactions with cells 2, 3, 4, and 7, while cell 9 interacts with cells 3, 6, 7, and 8.

[[File:ThirdYearCMPExpt-IsingNumberedCell.png|300px|thumb|right|'''Figure 2''': Illustration of a 2D Ising lattice, with numbered cells to indicate the neighbour interactions.]]

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

==Phase Transitions==

The Ising model is interesting because it is one of the simplest physical models to display a phase transition (so long as we are in two dimensions or more). The temperature of the system controls the balance between the energetic driving force, which wants the system to adopt the lowest energy configuration that you imagined in the previous section, and the entropic driving force, which seeks to maximise the number of configurations available to the system.

'''<big>TASK 1b</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 (''D''=3, ''N''=1000)? How much entropy does the system gain by doing so?'''

At low temperatures, little thermal energy is available to overcome such energy barriers, and we expect the system to adopt the lowest energy (the energetic driving force will dominate). Under these conditions, we expect a large number of parallel spins, and we can characterise this by considering the total ''magnetisation'' of the system: <math>M = \sum_i s_i</math>. At higher temperatures, however, enough thermal energy will be available that spins can flip relatively freely. Under these conditions, the entropic driving force should dominate. The transition temperature between the two regimes is called the ''Curie temperature'', ''TC''. Below ''TC'', we observe the highly ordered magnetised phase, while above ''TC'' a disordered, demagnetised, phase emerges.

'''<big>TASK 1c</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 ''D''=3, ''N''=1000 at absolute zero?'''

We are interested, therefore, in determining the magnetisation ''M'', and energy ''E'', of the system as functions of temperature ''T''. We can write statistical mechanical equations for these quantities:

<math>\left\langle M\right\rangle_T = \frac{1}{Z}\sum_\alpha M\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

<math>\left\langle E\right\rangle_T = \frac{1}{Z}\sum_\alpha E\left(\alpha\right) \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>

where ''α'' is used as a shorthand to represent all spins in the system, and ''Z'' is the partition function,
<math>Z = \sum_\alpha \exp \left\{-\frac{E\left(\alpha\right)}{k_BT}\right\} </math>.

It was discovered relatively quickly after the model was first proposed that the partition function in the one dimensional case can be solved analytically (though doing so is beyond the scope of this experiment):

<math>Z\left(T, N\right) = \left[2\cosh\left(\frac{J}{k_BT}\right)\right]^N</math>

Much later, [http://en.wikipedia.org/wiki/Lars_Onsager Lars Onsager] established that the two dimensional partition function can also be found analytically, though the mathematics is significantly more involved and the resulting expression is correspondingly much more complicated. In dimensions greater than two, no analytical solution is known.

Luckily for us, the problem can instead be tackled by numerical methods, and it is this that will be the focus of the experiment. We will discuss how the equations for ''M'' and ''E'' can be solved, but first we must set up a Python script to model the lattice.

'''<big>This is the first section of the Ising model experiment.</big>'''

'''<big>Previous section: [[Programming_a_2D_Ising_Model|Introduction page]]</big>'''

'''<big>Next section: [[Third year CMP compulsory experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]].</big>'''

Programming a 2D Ising Model

2026-02-12T12:20:11Z

Gmallia: /* Getting Help */

'''<big>This is the Ising Model experiment (aka programming for simple simulations). This is a core experiment for year 3 students on the Chemistry with Molecular Physics degree.</big>'''

==Introduction==

In this exercise, you are going to use the Python that you learned in the first year to write a code to perform Monte-Carlo simulations of the 2D [http://en.wikipedia.org/wiki/Ising_model Ising model], a set of spins on a lattice which is used to model ferromagnetic behaviour, and also to analyse the results of the simulation to find the heat capacity of the system and the Curie temperature — the temperature below which the system is able to maintain a spontaneous magnetisation.

A more extensive introduction to the Ising Model will be given in the next page. This page gives some practical information about how the experiment.

==Assessment==

At the end of this experiment you must submit a report, writing up your findings. Each section of the experiment has a number of tasks that you should complete, labelled '''<big>TASK</big>'''. If this is a mathematical exercise, your report should contain a short summary of the solution. If it is a graphical exercise, you should include the relevant image. Your report should explain briefly what you did in each stage of the experiment, and what your findings were.

The report should be submitted as a PDF file to Turnitin, created with your favourite word processor or LaTeX, and should contain the annotated source code of the Python scripts that your write during the experiment. For clarity, code included in the report should be written in a monospaced font (e.g. Courier), which in LaTeX can be easily achieved by placing your code inside a <code>\texttt{}</code> command. '''The code, except where relevant to discussion, can be placed in an appendix.'''

As well as including relevant source code in your report, you should also submit a bundle with all your scripts and analysis code as a separate Blackboard submission.

===On the use of AI tools===

The code you write for this experiment is assessed, and as per the [https://bb.imperial.ac.uk/bbcswebdav/pid-3344481-dt-content-rid-21115105_1/xid-21115105_1 department's guidelines on the use of AI], it is not acceptable to use generative AI tools to write computer code for this assignment.

Writing computer programs yourself enhances your logical thinking and problem solving skills. Using generative AI tools to write your code makes you dumb. Don't choose to be dumb, don't use AI to write code for you.

==Getting Help==

The Graduate Teaching Assistants demonstrating this experiment are Fionn Carman and Lei Li. Demonstrators will be available in room 232A at the following times:

{| class="wikitable"
|+ Demonstrator timetable
! !! Monday !! Tuesday !! Wednesday !! Thusrsday !! Friday
|-
| Morning || 10:00 - Experiment intro 10:00-13:00 - Fionn and Aidan || 10:00-13:00 - Ayse and Aidan|| || 10:00-13:00 - Fionn || 10:00-13:00 - Ayse and Aidan
|-
| Afternoon || 14:00-17:00 - Fionn|| 14:00-17:00 - Ayse and Aidan|| || 14:00-16:00 - Fionn and Aidan || 14:00-17:00 - Ayse
|}

If you have questions outside these times, please use the Discussion Board on Blackboard.

João and Giuseppe are the members of academic staff responsible for this experiment. You are welcome to address your questions and comments to them.

==Structure of this Experiment==

This experimental manual has been broken up into a number of subsections. Direct links to each of them may be found below. You should attempt them in order, and you should complete all of them to finish the experiment.

# [[Programming_a_2D_Ising_Model/Introduction_to_the_Ising_model|Introduction to the Ising model]]
# [[Programming_a_2D_Ising_Model/Calculating the energy and magnetisation|Calculating the energy and magnetisation]]
# [[Programming_a_2D_Ising_Model/Introduction to Monte Carlo simulation|Introduction to the Monte Carlo simulation]]
# [[Programming_a_2D_Ising_Model/Accelerating the code|Accelerating the code]]
# [[Programming_a_2D_Ising_Model/The effect of temperature|The effect of temperature]]
# [[Programming_a_2D_Ising_Model/The effect of system size|The effect of system size]]
# [[Programming_a_2D_Ising_Model/Determining the heat capacity|Determining the heat capacity]]
# [[Programming_a_2D_Ising_Model/Locating the Curie temperature|Locating the Curie temperature]]

Programming a 2D Ising Model

2026-02-12T12:19:45Z

Gmallia: /* Assessment */

'''<big>This is the Ising Model experiment (aka programming for simple simulations). This is a core experiment for year 3 students on the Chemistry with Molecular Physics degree.</big>'''

==Introduction==

In this exercise, you are going to use the Python that you learned in the first year to write a code to perform Monte-Carlo simulations of the 2D [http://en.wikipedia.org/wiki/Ising_model Ising model], a set of spins on a lattice which is used to model ferromagnetic behaviour, and also to analyse the results of the simulation to find the heat capacity of the system and the Curie temperature — the temperature below which the system is able to maintain a spontaneous magnetisation.

A more extensive introduction to the Ising Model will be given in the next page. This page gives some practical information about how the experiment.

==Assessment==

At the end of this experiment you must submit a report, writing up your findings. Each section of the experiment has a number of tasks that you should complete, labelled '''<big>TASK</big>'''. If this is a mathematical exercise, your report should contain a short summary of the solution. If it is a graphical exercise, you should include the relevant image. Your report should explain briefly what you did in each stage of the experiment, and what your findings were.

The report should be submitted as a PDF file to Turnitin, created with your favourite word processor or LaTeX, and should contain the annotated source code of the Python scripts that your write during the experiment. For clarity, code included in the report should be written in a monospaced font (e.g. Courier), which in LaTeX can be easily achieved by placing your code inside a <code>\texttt{}</code> command. '''The code, except where relevant to discussion, can be placed in an appendix.'''

As well as including relevant source code in your report, you should also submit a bundle with all your scripts and analysis code as a separate Blackboard submission.

===On the use of AI tools===

The code you write for this experiment is assessed, and as per the [https://bb.imperial.ac.uk/bbcswebdav/pid-3344481-dt-content-rid-21115105_1/xid-21115105_1 department's guidelines on the use of AI], it is not acceptable to use generative AI tools to write computer code for this assignment.

Writing computer programs yourself enhances your logical thinking and problem solving skills. Using generative AI tools to write your code makes you dumb. Don't choose to be dumb, don't use AI to write code for you.

==Getting Help==

The Graduate Teaching Assistants demonstrating this experiment are Fionn Carman, Ayse Halkseven and Aidan Chapman. Demonstrators will be available in room 232A at the following times:

{| class="wikitable"
|+ Demonstrator timetable
! !! Monday !! Tuesday !! Wednesday !! Thusrsday !! Friday
|-
| Morning || 10:00 - Experiment intro 10:00-13:00 - Fionn and Aidan || 10:00-13:00 - Ayse and Aidan|| || 10:00-13:00 - Fionn || 10:00-13:00 - Ayse and Aidan
|-
| Afternoon || 14:00-17:00 - Fionn|| 14:00-17:00 - Ayse and Aidan|| || 14:00-16:00 - Fionn and Aidan || 14:00-17:00 - Ayse
|}

If you have questions outside these times, please use the Discussion Board on Blackboard.

João and Giuseppe are the members of academic staff responsible for this experiment. You are welcome to address your questions and comments to them.

==Structure of this Experiment==

This experimental manual has been broken up into a number of subsections. Direct links to each of them may be found below. You should attempt them in order, and you should complete all of them to finish the experiment.

# [[Programming_a_2D_Ising_Model/Introduction_to_the_Ising_model|Introduction to the Ising model]]
# [[Programming_a_2D_Ising_Model/Calculating the energy and magnetisation|Calculating the energy and magnetisation]]
# [[Programming_a_2D_Ising_Model/Introduction to Monte Carlo simulation|Introduction to the Monte Carlo simulation]]
# [[Programming_a_2D_Ising_Model/Accelerating the code|Accelerating the code]]
# [[Programming_a_2D_Ising_Model/The effect of temperature|The effect of temperature]]
# [[Programming_a_2D_Ising_Model/The effect of system size|The effect of system size]]
# [[Programming_a_2D_Ising_Model/Determining the heat capacity|Determining the heat capacity]]
# [[Programming_a_2D_Ising_Model/Locating the Curie temperature|Locating the Curie temperature]]