ChemWiki - User contributions [en]

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

2025-03-04T15:26:14Z

Ac5116: Add the general expression of the partition function to avoid confusion.

'''<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, <math>s_i</math>, 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, <math>s_1</math> to <math>s_N</math>. 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=5x5), and three (N=5x5x5) 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 <math>J</math> is a constant which controls the strength of interaction and the notation <math>j\ \in\ \mathrm{neighbours}\left(i\right)</math> indicates that spin <math>j</math> must lie in a lattice cell adjacent to spin <math>i</math>. 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 <math>J>0</math>) is <math>E\ =\ -DNJ</math>, where <math>D</math> is the number of dimensions and <math>N</math> is the total number of spins. What is the multiplicity of this state? Calculate its entropy.'''

==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 (<math>D=3,\ N=1000</math>)? 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'', <math>T_C</math>. Below <math>T_C</math>, we observe the highly ordered magnetised phase, while above <math>T_C</math> 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 <math>D = 3,\ N=1000</math> at absolute zero?'''

We are interested, therefore, in determining the magnetisation <math>M</math>, and energy <math>E</math>, of the system as functions of temperature <math>T</math>. 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 <math>\alpha</math> is used as a shorthand to represent all spins in the system, and <math>Z</math> 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 <math>M</math> and <math>E</math> 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. 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>'''

Programming a 2D Ising Model

2025-03-03T10:00:49Z

Ac5116: Update hours

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

The mark breakdown for the tasks is given below. You will also be assessed on the quality of the code and the report.

'''Task 1a-c''': 6 marks

'''Task 2a-b''': 10 marks

'''Task 3a-c''': 12 marks

'''Task 4a-c''': 10 marks

'''Task 5a-b''': 12 marks

'''Task 6''': 10 marks

'''Task 7a-b''': 10 marks

'''Task 8a-d''': 20 marks

'''Report + Code''': 10 marks

'''Total:''' 100 marks

This mark will be divided by 5 and rounded to the nearest integer to obtain a final mark out of 20.

'''The deadline for report submission is friday the 14th of March at 18:00.'''

===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]]

Programming a 2D Ising Model/Locating the Curie temperature

2025-03-02T23:09:37Z

Ac5116: more <code> blocks. Added a note about the origin of the scaling equation.

'''<big>This is the eighth (and final) section of the Ising model experiment. You can return to the previous page, [[Programming a 2D Ising Model/Determining the heat capacity|Determining the heat capacity]], or go back to the [[Programming a 2D Ising Model|Introduction]].</big>'''

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 <math>T = T_C</math>. 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 <math>T_{C, L}</math> is the Curie temperature of an <math>L\times L</math>lattice, <math>T_{C,\infty}</math> is the Curie temperature of an infinite lattice, and <math>A</math> is a constant in which we are not especially interested. This scaling equation comes from expanding in the limit of large <math>L</math>.

'''<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: <math>T, E, E^2, M, M^2, C</math> (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 <math>T_{\mathrm{min}}</math> and <math>T_{\mathrm{max}}</math>, 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: peak_T_range, containing 1000 equally spaced temperature values between <math>T_{\mathrm{min}}</math> and <math>T_{\mathrm{max}}</math> (whatever you chose those values to be), and fitted_C_values, containing the fitted heat capacity at each of those points. Use the NumPy max function to find the maximum in C. If you store the maximum value of C in the variable Cmax, 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 colums: the lattice side length (2,4,8, etc.), and the temperature at which C is a maximum. This is your estimate of <math>T_C</math> for that side length. Make a plot that uses the scaling relation given above to determine <math>T_{C,\infty}</math>. By doing a little research, 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. You can return to the previous page, [[Programming a 2D Ising Model/Determining the heat capacity|Determining the heat capacity]], or go back to the [[Programming a 2D Ising Model|Introduction]].</big>'''

Programming a 2D Ising Model/Determining the heat capacity

2025-03-02T23:07:13Z

Ac5116: /* Calculating the heat capacity */

'''<big>This is the seventh section of the Ising model experiment. You can return to the previous page, [[Programming a 2D Ising Model/The effect of system size|The effect of system size]], or jump ahead to the next section, [[Programming a 2D Ising Model/Locating the Curie temperature|Locating the Curie temperature]].</big>'''

==Calculating the heat capacity==

As we have seen, increasing the temperature above <math>T_C</math> 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 <math>T = T_C</math> exactly, <math>C</math> 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, <math>\chi</math>, 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, <math>\mathrm{Var}[X]</math>, the mean of its square <math>\left\langle X^2\right\rangle</math>, and its squared mean <math>\left\langle X\right\rangle^2</math>. You may find that the data around the peak is very noisy — this is normal, and is a result of being in the critical region. As before, use the plot controls to save 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. You can return to the previous page, [[Programming a 2D Ising Model/The effect of system size|The effect of system size]], or jump ahead to the next section, [[Programming a 2D Ising Model/Locating the Curie temperature|Locating the Curie temperature]].</big>'''

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

2025-03-02T23:05:23Z

Ac5116: <code> blocks

'''<big>This is the sixth section of the Ising model experiment. You can return to the previous page, [[Programming a 2D Ising Model/The effect of temperature|The effect of temperature]], or jump ahead to the next section, [[Programming a 2D Ising Model/Determining the heat capacity|Determining the heat capacity]].</big>'''

==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. You can return to the previous page, [[Programming a 2D Ising Model/The effect of temperature|The effect of temperature]], or jump ahead to the next section, [[Programming a 2D Ising Model/Determining the heat capacity|Determining the heat capacity]].</big>'''

Programming a 2D Ising Model/The effect of temperature

2025-03-02T23:03:22Z

Ac5116: add code blocks. Warn about off-by-one errors.

'''<big>This is the fifth section of the Ising model experiment. You can return to the previous page, [[Programming a 2D Ising Model/Accelerating the code|Accelerating the code]], or jump ahead to the next section, [[Programming a 2D Ising Model/The effect of system size|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 <math>N_\mathrm{spins}</math> 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 <math>N</math> 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 ''np.arange'' function near the top of the file. Starting from a uniformly aligned lattice, a simulation of length n_steps will be performed at temperature T_cold, and the average energy and magnetisation will be recorded. A new simulation of length n_steps will then be performed at temperature T_cold + dT, 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: ave_energies, ave_sqenergies (squared energies), ave_magnetisations, ave_sqmagnetisations (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 savetxt() stores your array of output data on disk — you will need it later. Save the file as ''8x8.dat'' so that you know which lattice size it came from.'''

'''<big>This is the fifth section of the Ising model experiment. You can return to the previous page, [[Programming a 2D Ising Model/Accelerating the code|Accelerating the code]], or jump ahead to the next section, [[Programming a 2D Ising Model/The effect of system size|The effect of system size]].</big>'''

Programming a 2D Ising Model/Accelerating the code

2025-03-02T23:01:09Z

Ac5116: typos and grammar

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

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. You can return to the previous page, [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]], or jump ahead to the next section, [[Programming a 2D Ising Model/The effect of temperature|The effect of temperature]].</big>'''

Programming a 2D Ising Model/Accelerating the code

2025-03-02T23:00:23Z

Ac5116: Highlight where code should be added to the report. Point out difference between profile and timetrial

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

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 mesaures 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 ther 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. You can return to the previous page, [[Programming a 2D Ising Model/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]], or jump ahead to the next section, [[Programming a 2D Ising Model/The effect of temperature|The effect of temperature]].</big>'''

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

2025-03-02T22:56:06Z

Ac5116: Emphasise code adding code to report. Highlight code with <code> blocks. Correct reduced units defnition

'''<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 <math>\left\langle M\right\rangle_T</math>?'''

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 <math>\alpha</math> 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, <math>\alpha_0</math>, with energy <math>E_0</math>.
# Choose a single spin '''at random''', and "flip" it, to generate a new configuration <math>\alpha_1</math>
# Calculate the energy of this new configuration, <math>E_1</math>
# Calculate the energy difference between the states, <math>\Delta E = E_1 - E_0</math>
## If the <math>\Delta E \leq 0</math> (the spin flipping decreased the energy), then we '''accept''' the new configuration.
##* We set <math>\alpha_0 = \alpha_1</math>, and <math>E_0 = E_1</math>, and then '''go to step 5'''
## If <math>\Delta E > 0</math>, the spin flipping increased the energy. By considering the probability of observing the starting and final states, <math>\alpha_0</math> and <math>\alpha_1</math>, 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, <math>R</math>, in the interval <math>[0,1)</math>
### If <math>R \leq \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''accept''' the new configuration.
###* We set <math>\alpha_0 = \alpha_1</math>, and <math>E_0 = E_1</math>, and then '''go to step 5'''
### If <math>R > \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''reject''' the new configuration.
###* <math>\alpha_0</math> and <math>E_0</math> 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, <math>\Delta E</math> 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 <math>N_\mathrm{spins}</math>-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. You may assume that the energy returned by your <code>energy()</code> function is in units of <math>J</math> and the temperature is in units <math>T^* = k_BT/J</math>! Complete the <code>statistics()</code> function. This should return the following quantities whenever it is called: <math><E>, <E^2>, <M>, <M^2></math>, and <math>N_\mathrm{steps}</math>, the number of Monte Carlo steps that have elapsed. Put your code for these functions in the report.'''

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 <math>T=1.0</math>, 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 <math>T < T_C</math>, do you expect a spontaneous magnetisation (i.e. do you expect <math>\left\langle M\right\rangle \neq 0</math>)? When the state of the simulation appears to stop changing (when you have reached an equilibrium state), use the controls to export the output to 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 next section, [[Programming a 2D Ising Model/Accelerating the code|Accelerating the code]].</big>'''

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

2025-03-02T22:45:08Z

Ac5116: Added link to the Metropolis Paper

'''<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 <math>\left\langle M\right\rangle_T</math>?'''

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 <math>\alpha</math> 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, <math>\alpha_0</math>, with energy <math>E_0</math>.
# Choose a single spin '''at random''', and "flip" it, to generate a new configuration <math>\alpha_1</math>
# Calculate the energy of this new configuration, <math>E_1</math>
# Calculate the energy difference between the states, <math>\Delta E = E_1 - E_0</math>
## If the <math>\Delta E \leq 0</math> (the spin flipping decreased the energy), then we '''accept''' the new configuration.
##* We set <math>\alpha_0 = \alpha_1</math>, and <math>E_0 = E_1</math>, and then '''go to step 5'''
## If <math>\Delta E > 0</math>, the spin flipping increased the energy. By considering the probability of observing the starting and final states, <math>\alpha_0</math> and <math>\alpha_1</math>, 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, <math>R</math>, in the interval <math>[0,1)</math>
### If <math>R \leq \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''accept''' the new configuration.
###* We set <math>\alpha_0 = \alpha_1</math>, and <math>E_0 = E_1</math>, and then '''go to step 5'''
### If <math>R > \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''reject''' the new configuration.
###* <math>\alpha_0</math> and <math>E_0</math> 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, <math>\Delta E</math> 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 <math>N_\mathrm{spins}</math>-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, '''montecarlostep(T)''', 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: '''self.E_tally''', '''self.E2_tally''', '''self.M_tally''', '''self.M2_tally''', and '''self.n_steps'''. You should use these values to keep a running sum of all of 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 '''statistics()''', that takes no arguments.

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

You have been provided with a script, '''ILanim.py''', which will use your IsingLattice 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 <math>T=1.0</math>, 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 '''statistics()''' function to print the average energy and magnetisation from the simulation.

'''<big>TASK 3c</big>: If <math>T < T_C</math>, do you expect a spontaneous magnetisation (i.e. do you expect <math>\left\langle M\right\rangle \neq 0</math>)? When the state of the simulation appears to stop changing (when you have reached an equilibrium state), use the controls to export the output to 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 next section, [[Programming a 2D Ising Model/Accelerating the code|Accelerating the code]].</big>'''

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

2025-03-02T22:35:33Z

Ac5116: Emphasise code adding code to report. Highlight code with <code> blocks. Correct reduced units defnition

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

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

To perform some of the analysis tasks later in the experiment, you may want to use a Jupyter notebook. If running commands from a Jupyter notebook or IPython console, you should start the session by 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 you may need to restart the kernel)!

==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 <math>J=1.0</math> at all times (in fact, we are working in ''reduced units'' in which all energies are a multiple of <math>J</math>, but there will be more information about this in later sections). Do not worry about the efficiency of the code at the moment — we will address the speed in a later part of the experiment. 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. You are welcome to do this yourself — from JupyterLab you can open an IPython console and simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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>
These sorts of test are useful because they are easier to write and faster to run and can be run automatically.

'''<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/Locating the Curie temperature

2024-03-08T13:26:40Z

Ac5116:

'''<big>This is the eighth (and final) section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Determining the heat capacity|Determining the heat capacity]], or go back to the [[Third year CMP compulsory experiment|Introduction]].</big>'''

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 <math>T = T_C</math>. 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 <math>T_{C, L}</math> is the Curie temperature of an <math>L\times L</math>lattice, <math>T_{C,\infty}</math> is the Curie temperature of an infinite lattice, and <math>A</math> is a constant in which we are not especially interested.

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

==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 polyfit] and [http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyval.html#numpy.polyval polyval] 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 <math>\chi</math> vs T, as well as a fitted polynomial. Try changing the degree of the polynomial to improve the fit — in general, it might be difficult to get a good fit! Attach a PNG of an example fit to your report.'''

===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 <math>T_{\mathrm{min}}</math> and <math>T_{\mathrm{max}}</math>, 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: peak_T_range, containing 1000 equally spaced temperature values between <math>T_{\mathrm{min}}</math> and <math>T_{\mathrm{max}}</math> (whatever you chose those values to be), and fitted_C_values, containing the fitted heat capacity at each of those points. Use the NumPy max function to find the maximum in C. If you store the maximum value of C in the variable Cmax, 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 colums: the lattice side length (2,4,8, etc.), and the temperature at which C is a maximum. This is your estimate of <math>T_C</math> for that side length. Make a plot that uses the scaling relation given above to determine <math>T_{C,\infty}</math>. By doing a little research online, you should be able to find the theoretical exact Curie temperature for the infinite 2D Ising lattice. How does your value compare to this? Are you surprised by how good/bad the agreement is? Attach a PNG of this final graph to your report, and discuss briefly what you think the major sources of error are in your estimate.'''

'''<big>This is the eighth (and final) section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Determining the heat capacity|Determining the heat capacity]], or go back to the [[Third year CMP compulsory experiment|Introduction]].</big>'''

Programming a 2D Ising Model/Accelerating the code

2024-03-04T13:09:51Z

Ac5116:

'''<big>This is the fourth section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]], or jump ahead to the next section, [[Third year CMP compulsory experiment/The effect of temperature|The effect of temperature]].</big>'''

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 in C++, 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 C code 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>ILProfiling.py</code> script 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. The data is outputted into a file which can be visualised with [https://jiffyclub.github.io/snakeviz/ snakeviz].

If you run the script and open the file with snakeviz (you may first need to install it via anaconda),

<code>
snakeviz ising.prof
</code>

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.

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

'''<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 third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]], or jump ahead to the next section, [[Third year CMP compulsory experiment/The effect of temperature|The effect of temperature]].</big>'''

Programming a 2D Ising Model/Determining the heat capacity

2024-03-04T13:04:12Z

Ac5116: /* Calculating the heat capacity */

'''<big>This is the seventh section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/The effect of system size|The effect of system size]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Locating the Curie temperature|Locating the Curie temperature]].</big>'''

==Calculating the heat capacity==

As we have seen, increasing the temperature above <math>T_C</math> 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 <math>T = T_C</math> exactly, <math>C</math> 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, <math>\chi</math>, 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.

'''<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, <math>\mathrm{Var}[X]</math>, the mean of its square <math>\left\langle X^2\right\rangle</math>, and its squared mean <math>\left\langle X\right\rangle^2</math>. You may find that the data around the peak is very noisy — this is normal, and is a result of being in the critical region. As before, use the plot controls to save your a PNG image of your plot and attach this to the report.'''

'''<big>This is the seventh section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/The effect of system size|The effect of system size]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Locating the Curie temperature|Locating the Curie temperature]].</big>'''

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

2024-03-03T23:12:19Z

Ac5116: /* Getting the files for the experiment */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 Box [https://imperialcollegelondon.box.com/v/ising-lab-files here].''' If this link is broken you can alternatively try an older(slightly broken) version [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend you open the extracted folder with VSCode. This makes it easy to edit your code and run simulations all in one place.

We recommend that you use Conda Prompt (Windows) or Terminal(Mac/Linux) to run your simulation scripts. This should be available on all the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''Conda Prompt'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment.

'''NOTE: QTConsole does not work for running the ILanim.py script you run later in the lab. The best solution is to just run the script from the command line using '''python ILanim.py''' (if you're in the directory where files are stored). On a Mac/Linux computer you run this from '''Terminal''' or on windows you can use the app 'Conda Prompt''''

If running commands from a Jupyter Notebook or ipython terminal, you shuold start the session by 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 you may need to restart the console)!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your OneDrive drive, then in your terminal or console use cd H:\... to change to the correct directory.

Alternatively, open the folder in VS Code and create a new built-in terminal. Doing so opens the terminal in the correct folder.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<pre>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_steps = 0
</pre>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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() function of our IsingLattice object
</pre>

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

This line creates a variable called energy which contains whatever value was returned by the function energy() which is part of the IsingLattice class. If you just want to define a local variable instead, like an iterator in a loop, the "self." 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 function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice 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 this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 VSCode 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 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/en/8.0.x/ pytest] command in your terminal,
<pre>pytest</pre>
These sorts of test are useful because they are easier to write and faster to run and can be run automatically.

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

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

2024-03-03T23:07:14Z

Ac5116: /* Testing the files */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 Box [https://imperialcollegelondon.box.com/v/ising-lab-files here].''' If this link is broken you can alternatively try an older(slightly broken) version [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use Conda Prompt (Windows) or Terminal(Mac/Linux) to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''Conda Prompt'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment.

'''NOTE: QTConsole does not work for running the ILanim.py script you run later in the lab. The best solution is to just run the script from the command line using '''python ILanim.py''' (if you're in the directory where files are stored). On a Mac/Linux computer you run this from '''Terminal''' or on windows you can use the app 'Conda Prompt''''

When using the QTConsole, you should start the session by 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 you may need to restart the console)!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<pre>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_steps = 0
</pre>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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() function of our IsingLattice object
</pre>

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

This line creates a variable called energy which contains whatever value was returned by the function energy() which is part of the IsingLattice class. If you just want to define a local variable instead, like an iterator in a loop, the "self." 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 function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice 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 this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 VSCode 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 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/en/8.0.x/ pytest] command in your terminal,
<pre>pytest</pre>
These sorts of test are useful because they are easier to write and faster to run and can be run automatically.

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

Programming a 2D Ising Model/The effect of temperature

2024-03-03T22:58:19Z

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

'''<big>This is the fifth section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Accelerating the code|Accelerating the code]], or jump ahead to the next section, [[Third year CMP compulsory experiment/The effect of system size|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. How many cycles (remember a cycle is <math>N_\mathrm{spins}</math> steps) are ''typically'' needed for the system to go from its random starting position to the equilibrium state? Modify your statistics() and montecarlostep() functions so that the first N cycles of the simulation are ignored when calculating the averages. You should state in your report what period you chose to ignore, and include graphs from ILfinalframe.py to illustrate your motivation in choosing this 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 ''np.arange'' function near the top of the file. Starting from a random lattice, a simulation of length n_steps will be performed at temperature T_cold, and the average energy and magnetisation will be recorded. A new simulation of length n_steps will then be performed at temperature T_cold + dT, 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: ave_energies, ave_sqenergies (squared energies), ave_magnetisations, ave_sqmagnetisations (squared magnetisations).

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

'''<big>This is the fifth section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Accelerating the code|Accelerating the code]], or jump ahead to the next section, [[Third year CMP compulsory experiment/The effect of system size|The effect of system size]].</big>'''

Programming a 2D Ising Model/The effect of temperature

2024-03-03T22:56:14Z

Ac5116: /* Correcting the Averaging Code */

'''<big>This is the fifth section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Accelerating the code|Accelerating the code]], or jump ahead to the next section, [[Third year CMP compulsory experiment/The effect of system size|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. How many cycles (remember a cycle is <math>N_\mathrm{spins}</math> steps) are ''typically'' needed for the system to go from its random starting position to the equilibrium state? Modify your statistics() and montecarlostep() functions so that the first N cycles of the simulation are ignored when calculating the averages. You should state in your report what period you chose to ignore, and include graphs from ILfinalframe.py to illustrate your motivation in choosing this 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 ''np.arange'' function near the top of the file. Starting from a random lattice, a simulation of length n_steps will be performed at temperature T_cold, and the average energy and magnetisation will be recorded. A new simulation of length n_steps will then be performed at temperature T_cold + dT, 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: ave_energies, ave_sqenergies (squared energies), ave_magnetisations, ave_sqmagnetisations (squared magnetisations).

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

'''<big>This is the fifth section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Accelerating the code|Accelerating the code]], or jump ahead to the next section, [[Third year CMP compulsory experiment/The effect of system size|The effect of system size]].</big>'''

Programming a 2D Ising Model/Accelerating the code

2024-03-03T22:51:06Z

Ac5116: Added extra details about why numpy has better performance and extension for delta energy method

'''<big>This is the fourth section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]], or jump ahead to the next section, [[Third year CMP compulsory experiment/The effect of temperature|The effect of temperature]].</big>'''

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 in C++, 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 C code 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>ILProfiling.py</code> script 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. The data is outputted into a file which can be visualised with [https://jiffyclub.github.io/snakeviz/ snakeviz].

If you run the script and open the file with snakeviz (you may first need to install it via anaconda),

<code>
snakeviz ising.prof
</code>

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.

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

'''<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> Extension </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 third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]], or jump ahead to the next section, [[Third year CMP compulsory experiment/The effect of temperature|The effect of temperature]].</big>'''

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

2024-03-03T22:43:03Z

Ac5116: /* Modifying IsingLattice.py */

'''<big>This is the third section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]], or jump ahead to the next section, [[Third year CMP compulsory experiment/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 <math>1\times 10^9</math> configurations per second with our computer. How long will it take to evaluate a single value of <math>\left\langle M\right\rangle_T</math>?'''

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 <math>\alpha</math> 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 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 it first in its list of the [http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=814652 "Top 10 Algorithms of the 20th Century"].

The algorithm is as follows:

# Start from a given configuration of spins, <math>\alpha_0</math>, with energy <math>E_0</math>.
# Choose a single spin '''at random''', and "flip" it, to generate a new configuration <math>\alpha_1</math>
# Calculate the energy of this new configuration, <math>E_1</math>
# Calculate the energy difference between the states, <math>\Delta E = E_1 - E_0</math>
## If the <math>\Delta E \leq 0</math> (the spin flipping decreased the energy), then we '''accept''' the new configuration.
##* We set <math>\alpha_0 = \alpha_1</math>, and <math>E_0 = E_1</math>, and then '''go to step 5'''
## If <math>\Delta E > 0</math>, the spin flipping increased the energy. By considering the probability of observing the starting and final states, <math>\alpha_0</math> and <math>\alpha_1</math>, 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, <math>R</math>, in the interval <math>[0,1)</math>
### If <math>R \leq \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''accept''' the new configuration.
###* We set <math>\alpha_0 = \alpha_1</math>, and <math>E_0 = E_1</math>, and then '''go to step 5'''
### If <math>R > \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''reject''' the new configuration.
###* <math>\alpha_0</math> and <math>E_0</math> 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, <math>\Delta E</math> 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 <math>N_\mathrm{spins}</math>-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 "Understanding Molecular Simulation", by Frenkel and Smit, but a discussion of this is not required for the experiment. You can find this as an ebook in the library [https://library-search.imperial.ac.uk/permalink/f/o1297h/TN_cdi_askewsholts_vlebooks_9780080519982 here]. Chapter 7 of 'Statistical Mechanics: Theory and Molecular Simulation' by Tuckerman also gives a good description.

==Modifying IsingLattice.py==

Our IsingLattice object contains a stub function, '''montecarlostep(T)''', which takes a single argument — the temperature at which we want to do the cycle. There are also attributes to measure the average energy and magnetisation: '''self.E''', '''self.E2''', '''self.M''', '''self.M2''', and '''self.n_steps'''. You should use these values to keep a running sum of all of 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 '''statistics()''', that takes no arguments.

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

You have been provided with a script, '''ILanim.py''', which will use your IsingLattice object to run a Monte Carlo simulation and display the output on the screen in real time. By default, this simulation will run an <math>8\times 8</math> lattice at <math>T=1.0</math>, 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 '''statistics()''' function to print the average energy and magnetisation from the simulation.

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

'''<big>This is the third section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Accelerating the code|Accelerating the code]].</big>'''

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

2024-03-03T22:41:02Z

Ac5116: /* A note about Python Classes */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 Box [https://imperialcollegelondon.box.com/v/ising-lab-files here].''' If this link is broken you can alternatively try an older(slightly broken) version [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use Conda Prompt (Windows) or Terminal(Mac/Linux) to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''Conda Prompt'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment.

'''NOTE: QTConsole does not work for running the ILanim.py script you run later in the lab. The best solution is to just run the script from the command line using '''python ILanim.py''' (if you're in the directory where files are stored). On a Mac/Linux computer you run this from '''Terminal''' or on windows you can use the app 'Conda Prompt''''

When using the QTConsole, you should start the session by 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 you may need to restart the console)!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<pre>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_steps = 0
</pre>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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() function of our IsingLattice object
</pre>

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

This line creates a variable called energy which contains whatever value was returned by the function energy() which is part of the IsingLattice class. If you just want to define a local variable instead, like an iterator in a loop, the "self." 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 function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice 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 this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 IPython Qt console using the command '''<pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

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/en/8.0.x/ pytest] command in your terminal,
<pre>pytest</pre>
These sorts of test are useful because they are easier to write and faster to run and can be run automatically.

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

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

2024-03-03T22:40:17Z

Ac5116: /* Testing the files */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 Box [https://imperialcollegelondon.box.com/v/ising-lab-files here].''' If this link is broken you can alternatively try an older(slightly broken) version [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use Conda Prompt (Windows) or Terminal(Mac/Linux) to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''Conda Prompt'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment.

'''NOTE: QTConsole does not work for running the ILanim.py script you run later in the lab. The best solution is to just run the script from the command line using '''python ILanim.py''' (if you're in the directory where files are stored). On a Mac/Linux computer you run this from '''Terminal''' or on windows you can use the app 'Conda Prompt''''

When using the QTConsole, you should start the session by 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 you may need to restart the console)!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<code>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_steps = 0
</code>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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:

<code>
import IsingLattice

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

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

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

This line creates a variable called energy which contains whatever value was returned by the function energy() which is part of the IsingLattice class. If you just want to define a local variable instead, like an iterator in a loop, the "self." 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 function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice 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 this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 IPython Qt console using the command '''<pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

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/en/8.0.x/ pytest] command in your terminal,
<pre>pytest</pre>
These sorts of test are useful because they are easier to write and faster to run and can be run automatically.

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

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

2024-03-03T22:34:48Z

Ac5116: /* The Constructor */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 Box [https://imperialcollegelondon.box.com/v/ising-lab-files here].''' If this link is broken you can alternatively try an older(slightly broken) version [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use Conda Prompt (Windows) or Terminal(Mac/Linux) to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''Conda Prompt'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment.

'''NOTE: QTConsole does not work for running the ILanim.py script you run later in the lab. The best solution is to just run the script from the command line using '''python ILanim.py''' (if you're in the directory where files are stored). On a Mac/Linux computer you run this from '''Terminal''' or on windows you can use the app 'Conda Prompt''''

When using the QTConsole, you should start the session by 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 you may need to restart the console)!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<code>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_steps = 0
</code>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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:

<code>
import IsingLattice

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

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

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

This line creates a variable called energy which contains whatever value was returned by the function energy() which is part of the IsingLattice class. If you just want to define a local variable instead, like an iterator in a loop, the "self." 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 function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice 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 this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 IPython Qt console using the command '''<pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

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

2024-03-03T22:33:51Z

Ac5116: /* A note about Python Classes */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 Box [https://imperialcollegelondon.box.com/v/ising-lab-files here].''' If this link is broken you can alternatively try an older(slightly broken) version [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use Conda Prompt (Windows) or Terminal(Mac/Linux) to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''Conda Prompt'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment.

'''NOTE: QTConsole does not work for running the ILanim.py script you run later in the lab. The best solution is to just run the script from the command line using '''python ILanim.py''' (if you're in the directory where files are stored). On a Mac/Linux computer you run this from '''Terminal''' or on windows you can use the app 'Conda Prompt''''

When using the QTConsole, you should start the session by 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 you may need to restart the console)!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<code>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_steps = 0
</code>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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:

<code>
import IsingLattice

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

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

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

This line creates a variable called energy which contains whatever value was returned by the function energy() which is part of the IsingLattice class. If you just want to define a local variable instead, like an iterator in a loop, the "self." 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 function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice is then created, which is a NumPy array with the requested number of rows and columns. The function <pre>np.random.choice</pre> is used to create a random collection of spins for our initial configuration. You can see the documentation for this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 IPython Qt console using the command '''<pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

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

2024-03-03T22:30:03Z

Ac5116: /* A note about Python Classes */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 Box [https://imperialcollegelondon.box.com/v/ising-lab-files here].''' If this link is broken you can alternatively try an older(slightly broken) version [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use Conda Prompt (Windows) or Terminal(Mac/Linux) to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''Conda Prompt'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment.

'''NOTE: QTConsole does not work for running the ILanim.py script you run later in the lab. The best solution is to just run the script from the command line using '''python ILanim.py''' (if you're in the directory where files are stored). On a Mac/Linux computer you run this from '''Terminal''' or on windows you can use the app 'Conda Prompt''''

When using the QTConsole, you should start the session by 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 you may need to restart the console)!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<code>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_steps = 0
</code>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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:

<code>
import IsingLattice

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

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

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

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

<code>

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

</code>

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

===The Constructor===

The only function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice is then created, which is a NumPy array with the requested number of rows and columns. The function <pre>np.random.choice</pre> is used to create a random collection of spins for our initial configuration. You can see the documentation for this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 IPython Qt console using the command '''<pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

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

2024-03-03T22:28:51Z

Ac5116: /* Importance Sampling */

'''<big>This is the third section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]], or jump ahead to the next section, [[Third year CMP compulsory experiment/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 <math>1\times 10^9</math> configurations per second with our computer. How long will it take to evaluate a single value of <math>\left\langle M\right\rangle_T</math>?'''

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 <math>\alpha</math> 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 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 it first in its list of the [http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=814652 "Top 10 Algorithms of the 20th Century"].

The algorithm is as follows:

# Start from a given configuration of spins, <math>\alpha_0</math>, with energy <math>E_0</math>.
# Choose a single spin '''at random''', and "flip" it, to generate a new configuration <math>\alpha_1</math>
# Calculate the energy of this new configuration, <math>E_1</math>
# Calculate the energy difference between the states, <math>\Delta E = E_1 - E_0</math>
## If the <math>\Delta E \leq 0</math> (the spin flipping decreased the energy), then we '''accept''' the new configuration.
##* We set <math>\alpha_0 = \alpha_1</math>, and <math>E_0 = E_1</math>, and then '''go to step 5'''
## If <math>\Delta E > 0</math>, the spin flipping increased the energy. By considering the probability of observing the starting and final states, <math>\alpha_0</math> and <math>\alpha_1</math>, 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, <math>R</math>, in the interval <math>[0,1)</math>
### If <math>R \leq \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''accept''' the new configuration.
###* We set <math>\alpha_0 = \alpha_1</math>, and <math>E_0 = E_1</math>, and then '''go to step 5'''
### If <math>R > \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''reject''' the new configuration.
###* <math>\alpha_0</math> and <math>E_0</math> 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, <math>\Delta E</math> 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 <math>N_\mathrm{spins}</math>-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 "Understanding Molecular Simulation", by Frenkel and Smit, but a discussion of this is not required for the experiment. You can find this as an ebook in the library [https://library-search.imperial.ac.uk/permalink/f/o1297h/TN_cdi_askewsholts_vlebooks_9780080519982 here]. Chapter 7 of 'Statistical Mechanics: Theory and Molecular Simulation' by Tuckerman also gives a good description.

==Modifying IsingLattice.py==

Our IsingLattice object contains a stub function, '''montecarlocycle(T)''', which takes a single argument — the temperature at which we want to do the cycle. There are also attributes to measure the average energy and magnetisation: '''self.E''', '''self.E2''', '''self.M''', '''self.M2''', and '''self.n_cycles'''. You should use these values to keep a running sum of all of the energies, squared energies, magnetisations, and squared magnetisations that your system has sampled, as well as the number of cycles that have been performed. Finally, there is also a function '''statistics()''', that takes no arguments.

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

You have been provided with a script, '''ILanim.py''', which will use your IsingLattice object to run a Monte Carlo simulation and display the output on the screen in real time. By default, this simulation will run an <math>8\times 8</math> lattice at <math>T=1.0</math>, 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 '''statistics()''' function to print the average energy and magnetisation from the simulation.

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

'''<big>This is the third section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Accelerating the code|Accelerating the code]].</big>'''

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

2023-03-13T15:09:50Z

Ac5116: /* Getting the files for the experiment */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 Box [https://imperialcollegelondon.box.com/v/ising-lab-files here].''' If this link is broken you can alternatively try an older(slightly broken) version [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use Conda Prompt (Windows) or Terminal(Mac/Linux) to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''Conda Prompt'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment.

'''NOTE: QTConsole does not work for running the ILanim.py script you run later in the lab. The best solution is to just run the script from the command line using '''python ILanim.py''' (if you're in the directory where files are stored). On a Mac/Linux computer you run this from '''Terminal''' or on windows you can use the app 'Conda Prompt''''

When using the QTConsole, you should start the session by 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 you may need to restart the console)!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<code>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_cycles = 0
</code>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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:

<code>
import IsingLattice

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

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

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

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

<code>

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

</code>

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

===The Constructor===

The only function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice is then created, which is a NumPy array with the requested number of rows and columns. The function <pre>np.random.choice</pre> is used to create a random collection of spins for our initial configuration. You can see the documentation for this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 IPython Qt console using the command '''<pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

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

2023-03-13T15:08:36Z

Ac5116: /* Getting the files for the experiment */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 Box [https://imperialcollegelondon.box.com/v/ising-lab-files here].''' If this link is broken you can alternatively try an older(slightly broken) version [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use the IPython Qt or the regular ipython console (which you can access via spyder) to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''IPython QTConsole'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment.

'''NOTE: QTConsole does not work for running the ILanim.py script you run later in the lab. The best solution is to just run the script from the command line using '''python ILanim.py''' (if you're in the directory where files are stored). On a Mac/Linux computer you run this from '''Terminal''' or on windows you can use the app 'Conda Prompt''''

When using the QTConsole, you should start the session by 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 you may need to restart the console)!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<code>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_cycles = 0
</code>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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:

<code>
import IsingLattice

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

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

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

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

<code>

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

</code>

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

===The Constructor===

The only function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice is then created, which is a NumPy array with the requested number of rows and columns. The function <pre>np.random.choice</pre> is used to create a random collection of spins for our initial configuration. You can see the documentation for this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 IPython Qt console using the command '''<pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

Programming a 2D Ising Model

2023-03-10T15:56:27Z

Ac5116: /* Getting Help */

'''<big>This is the Ising Model experiment (aka programming for simple simulations). This is a core experiment for students on the Chemistry with Molecular Physics degree. If you are not a CMP student, you might be looking for the 3rd year [[Third_year_simulation_experiment|Liquid simulation experiment]].</big>'''

==Introduction==

In the first year, you were introduced to [https://www.python.org/ Python], a scripting language which is rapidly becoming the de facto language for everyday scientific programming. Python is an interpreted, rather than compiled language, and is rather more forgiving than older languages such as C or FORTRAN. This reduces the amount of time that we need to spend programming and debugging. The downside to all this is that the execution of a Python program is much slower than a compiled equivalent. As a compromise, we usually let large codes written in a compiled language (usually C/C++) do the hefty numerical work for us, and then use scripting languages like Python to analyse the results.

These large codes for numerical work (you may have already used GAUSSIAN for electronic structure calculations) typically take arcane text files as input, and produce equally arcane text files as output. If, for the sake of example, you run twenty different molecular dynamics simulations, and each of them produces an output file which contains information about the density of the system, then extracting this information by hand would be very tedious (and if you run hundreds or thousands of simulations, virtually impossible), but this sort of task is the thing at which languages like Python really excel.

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.

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

'''YOU MUST SUBMIT YOUR REPORT ON BLACKBOARD BY 12 NOON ON THE WEDNESDAY FOLLOWING THE EXPERIMENT.'''

'''Your report should contain the source code to *any* Python scripts that you write during the experiment, annotated with comments. You can write your report in Word or LaTeX, but you should submit a **single PDF file**.'''

The mark breakdown for the tasks is given below. You will also be assessed on the quality of the code and the report.

'''Task 1a-c''': 6 marks

'''Task 2a-b''': 10 marks

'''Task 3a-c''': 12 marks

'''Task 4a-c''': 10 marks

'''Task 5a-b''': 12 marks

'''Task 6''': 10 marks

'''Task 7a-b''': 10 marks

'''Task 8a-d''': 20 marks

'''Report + Code''': 10 marks

'''Total:''' 100 marks

This mark will be divided by 5 and rounded to the nearest integer to obtain a final mark out of 20.

For clarity, code included in the report should be written in a monospaced font (e.g. Courier) or a pasted screenshot. The former can be easily achieved in LaTeX by placing your code inside a <code>\texttt{}</code> command. '''The code, except where relevant to discussion, can be placed in an appendix.'''

If you wish, you are welcome to do the experiment on your own computer. You will need the [https://www.anaconda.com/products/individual Anaconda] scientific Python distribution (if you are a Windows, Mac or Linux user). You may ask for help with installing this, but it is not part of the experiment — students who have scientific questions will take priority.

==Getting Help==

The demonstrators for this exercise are Aidan Chapman (aidan.chapman16@imperial.ac.uk) and Jonathan Hedley (jonathan.hedley16@imperial.ac.uk). The assessors will be Aidan Chapman, Jonathan Hedley and Prof. Fernando Bresme. The demonstrators will be available in person between 10-12 and 2-3 pm on Monday and Thursday and 11-12 and 2-4pm on Tuesday and Friday ('''Note: the experiment does not "run" on Wednesdays'''). If you have questions outside these times, you are welcome to send them via e-mail.

The member of academic staff responsible for this exercise is Prof. Fernando Bresme (f.bresme@imperial.ac.uk).

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

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

Third year simulation experiment/Software Instructions

2021-10-22T12:27:28Z

Ac5116:

==Running LAMMPS==

===Command line basics ===

==VMD==

==Reading LAMMPS Log Files==

==Installing LAMMPS On Your Own Machine (OPTIONAL)==

==Accessing the Virtual Machines==

Third year simulation experiment/Software Instructions

2021-10-22T12:26:19Z

Ac5116: Created page with "==Running LAMMPS== ===Command line basics === ==VMD== ==Reading LAMMPS Log Files== ==Installing LAMMPS On Your Own Machine (OPTIONAL)=="

==Running LAMMPS==

===Command line basics ===

==VMD==

==Reading LAMMPS Log Files==

==Installing LAMMPS On Your Own Machine (OPTIONAL)==

Third year simulation experiment/Structural properties and the radial distribution function

2021-10-20T13:31:52Z

Ac5116: /* Simulations in this section */

'''<big>This is the fifth section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Running_simulations_under_specific_conditions|Running simulations under specific conditions]], or jump ahead to the next section, [[Third year simulation experiment/Dynamical properties and the diffusion coefficient|Dynamical properties and the diffusion coefficient]].</big>'''

We can characterise the structure of systems that we simulate using [http://en.wikipedia.org/wiki/Radial_distribution_function radial distribution functions], which we denote <math>g(r)</math>. Calculating the RDF for a simulation is very useful — it can tell us the distances from an atom at which we will find it's nearest neighbour, second nearest neighbour, and so on; it is also a quantity that can be accessed experimentally, and so provides a good check that the forcefield in our simulation is correctly reproducing the structural features.

In this section, you are going to use VMD to calculate the radial distribution function for the solid, liquid, and vapour phases of the Lennard-Jones fluid.

===Simulations in this section===

The '''RDF''' subfolder contains an example input script that you can use to record an atomic trajectory to generate RDFs for the solid, liquid, and vapour phase Lennard Jones systems. Make three copies of that script (one for each phase), and modify the density and temperature parameters to give the phase that you want (a phase diagram for the Lennard-Jones system can be found [http://journals.aps.org/pr/abstract/10.1103/PhysRev.184.151 here]). <big>'''Note: when simulating the solid, you will need to change the lattice type in the lattice command to fcc, rather than sc. '''</big>

'''<big>TASK 9: perform simulations of the Lennard-Jones system in the three phases. When each is complete, open the trajectories (.dump files) and calculate <math>g(r)</math> and <math>4\pi \int g(r) r^{2}\mathrm{d}r</math>. See instructions at the end of this page.</big>'''

=== Briefly explain what a Radial Distribution Function is. What is the relationship between the coordination number and RDF? [2] ===

=== Plot the RDFs for the three systems. [2] ===

=== Discuss qualitatively the differences between the three RDFs, and what this tells you about the structure of the system in each phase. [5] ===

=== In the solid case, illustrate which lattice sites the first three peaks correspond to [2]. What is the lattice spacing [1]? What are the coordination number for each of the first three peaks [1]? ===
'''<big>This is the sixth section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Running simulations under specific conditions|Running simulations under specific conditions]], or jump ahead to the next section, [[Third year simulation experiment/Dynamical properties and the diffusion coefficient|Dynamical properties and the diffusion coefficient]].</big>'''

===Calculating <math>g(r)</math> in VMD===

# Start VMD as before and load the trajectory that you want to analyse.
# Select '''Extensions''' -> '''Analysis''' -> '''Radial Pair Distribution Function g(r)'''
# Set '''Selection 1''' to '''all''' and '''Selection 2''' to '''all'''
# Change '''delta r''' to '''0.05''' — this is the distance between points in the generated RDF.
# Ensure that '''Use PBC''', '''Display g(r)''', '''Display Int(g(r))''', and '''Save to File''' are checked
# Click '''Compute g(r)'''

After a short pause while it performs the calculation, VMD will display both the RDF, and its running integral. You will then be prompted to save this data — choose a location for the file that you will be able to find easily later.

Third year simulation experiment/Structural properties and the radial distribution function

2021-10-20T13:29:39Z

Ac5116: /* Simulations in this section */

'''<big>This is the fifth section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Running_simulations_under_specific_conditions|Running simulations under specific conditions]], or jump ahead to the next section, [[Third year simulation experiment/Dynamical properties and the diffusion coefficient|Dynamical properties and the diffusion coefficient]].</big>'''

We can characterise the structure of systems that we simulate using [http://en.wikipedia.org/wiki/Radial_distribution_function radial distribution functions], which we denote <math>g(r)</math>. Calculating the RDF for a simulation is very useful — it can tell us the distances from an atom at which we will find it's nearest neighbour, second nearest neighbour, and so on; it is also a quantity that can be accessed experimentally, and so provides a good check that the forcefield in our simulation is correctly reproducing the structural features.

In this section, you are going to use VMD to calculate the radial distribution function for the solid, liquid, and vapour phases of the Lennard-Jones fluid.

===Simulations in this section===

The '''RDF''' subfolder contains an example input script that you can use to record an atomic trajectory to generate RDFs for the solid, liquid, and vapour phase Lennard Jones systems. Make three copies of that script (one for each phase), and modify the density and temperature parameters to give the phase that you want (a phase diagram for the Lennard-Jones system can be found [http://journals.aps.org/pr/abstract/10.1103/PhysRev.184.151 here]). <big>'''Note: when simulating the solid, you will need to change the lattice type in the lattice command to fcc, rather than sc.'''</big>

'''<big>TASK 9: perform simulations of the Lennard-Jones system in the three phases. When each is complete, download the trajectory and calculate <math>g(r)</math> and <math>4\pi \int g(r) r^{2}\mathrm{d}r</math>. See instructions at the end of this page.</big>'''

=== Briefly explain what a Radial Distribution Function is. What is the relationship between the coordination number and RDF? [2] ===

=== Plot the RDFs for the three systems. [2] ===

=== Discuss qualitatively the differences between the three RDFs, and what this tells you about the structure of the system in each phase. [5] ===

=== In the solid case, illustrate which lattice sites the first three peaks correspond to [2]. What is the lattice spacing [1]? What are the coordination number for each of the first three peaks [1]? ===
'''<big>This is the sixth section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Running simulations under specific conditions|Running simulations under specific conditions]], or jump ahead to the next section, [[Third year simulation experiment/Dynamical properties and the diffusion coefficient|Dynamical properties and the diffusion coefficient]].</big>'''

===Calculating <math>g(r)</math> in VMD===

# Start VMD as before and load the trajectory that you want to analyse.
# Select '''Extensions''' -> '''Analysis''' -> '''Radial Pair Distribution Function g(r)'''
# Set '''Selection 1''' to '''all''' and '''Selection 2''' to '''all'''
# Change '''delta r''' to '''0.05''' — this is the distance between points in the generated RDF.
# Ensure that '''Use PBC''', '''Display g(r)''', '''Display Int(g(r))''', and '''Save to File''' are checked
# Click '''Compute g(r)'''

After a short pause while it performs the calculation, VMD will display both the RDF, and its running integral. You will then be prompted to save this data — choose a location for the file that you will be able to find easily later.

Third year simulation experiment/Running simulations under specific conditions

2021-10-20T13:19:41Z

Ac5116: /* Temperature and Pressure Control */

'''<big>This is the fourth section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Equilibration|Equilibration]], or jump ahead to the next section, [[Third_year_simulation_experiment/Structural_properties_and_the_radial_distribution_function| Structural Properties and the Radial Distribution Functions]].</big>'''

'''<big>THE FILES THAT YOU NEED FOR THIS SECTION ARE FOUND IN THE "NpT" SUBFOLDER.</big>'''

==Changing Ensemble==

So far, we have been able to do some simulations in which the number of particles and the volume of the simulation cell are held constant. The energy is also constant (within a certain degree of error, which is introduced by the approximations that we make to do the simulation). If the simulation is a working properly, then the pressure and temperature of the system should also reach a constant ''average'' value (although there will again be fluctuations). In the statistical thermodynamics lectures, you met the concept of ensembles, which are used in statistical mechanics to represent different sorts of experimental conditions. The simulations we have done so far are described by the ''microcanonical'', or NVE ensemble (the letters represent those thermodynamic quantities which are constant).

As chemists, we often want to understand what happens under particular experimental conditions — at 298K under 1 atmosphere of pressure, for example. These sorts of conditions are described by different ensembles in statistical mechanics, such as the NVT (''canonical'') or NpT (''isobaric-isothermal'') ensembles.

In this section, we are going to modify our simulations from the previous section to run under NpT conditions, and sketch an equation of state for our model fluid at atmospheric pressure.

==Temperature and Pressure Control==

The file npt.in can be used to perform a constant temperature/pressure simulation of our model fluid. It starts by melting a simple cubic crystal, just as before, so much of this file will look familiar to you. You will notice a new section near the top, however, called '''### SPECIFY THE REQUIRED THERMODYNAMIC STATE ###'''. It contains three ''variables'' — these are used by the script later on to define the desired temperature, pressure, and timestep. The ellipses need to be replaced by the actual temperature, pressure and timestep that you want to use, so

<pre>
variable T equal 0.5
variable p equal 1.0
variable timestep equal 0.75
</pre>

would run a simulation at <math>T=0.5,\ p=1.0,\ \delta t=0.75</math>. You should remember from the [[Third_year_simulation_experiment/Equilibration|Equilibration]] section that this is a poor choice of timestep!

'''<big>TASK 8: Choose 5 temperatures (above the critical temperature <math>T^* = 1.5</math>), and two pressures (you can get a good idea of what a reasonable pressure is in Lennard-Jones units by looking at the average pressure of your simulations from the last section). This gives ten phase points — five temperatures at each pressure. Create 10 copies of npt.in, and modify each to run a simulation at one of your chosen <math>\left(p, T\right)</math> points. You should be able to use the results of the previous section to choose a timestep. Run these 10 scripts on the virtual machine. When your simulations have finished, download the log files as before. At the end of the log file, LAMMPS will output the values and errors for the pressure, temperature, and density <math>\left(\frac{N}{V}\right)</math>. </big>'''

=== Plot the density as a function of temperature for both pressures that you simulated. Include a line corresponding to the predictions made by the ideal gas law. [3] ===

=== How do your results compare to the ideal gas law? Do deviations increase/decrease with temperature and pressure? Explain. [7] ===

=== Do you expect your simulation results to be in better or worse agreement with the Van der Waals equation of state? Why? [3] ===

===Thermostats and Barostats - controlling the thermodynamic properties===
The statistical thermodynamics lectures will have introduced you to the ''equipartition theorem'', which states that, on average, every degree of freedom in a system at equilibrium will have <math>\frac{1}{2}k_B T</math> of energy. In our system with <math>N</math> atoms, each with 3 degrees of freedom, we can write
<math>E_K = \frac{3}{2} N k_B T</math>

<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T</math>

At the end of every timestep, we use the left hand side of this equation to calculate the kinetic energy, then divide by <math>\frac{3}{2}Nk_B</math> to get the ''instantaneous'' temperature <math>T</math>. In general, <math>T</math> will fluctuate, and will be different to our ''target'' temperature, <math>\mathfrak{T}</math> (this is whatever value we specify in the input script). We can change the temperature by multiplying every velocity by a constant factor, <math>\gamma</math>.

* If <math> T > \mathfrak{T} </math>, then the kinetic energy of the system is too high, and we need to reduce it. <math>\gamma < 1</math>
* If <math> T < \mathfrak{T} </math>, then the kinetic energy of the system is too low, and we need to increase it. <math>\gamma > 1</math>

We need to choose a scaling parameter <math>\gamma</math> so that the temperature is correct <math>T = \mathfrak{T}</math> if we multiply every velocity <math>\gamma</math>. We can write two equations:'''
<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T</math>

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B \mathfrak{T}</math>

By combining these equations, one can see that <math> \gamma = \sqrt{\frac{\mathfrak{T}}{T}} </math> (satisfy yourself that this is true!). A target value of <math> \gamma </math> of 1 is required and thus, dependent on whether it's larger or smaller than 1 the simulation can target the desired temperature.

Controlling the pressure is a little more involved, but the principle is largely the same: at each timestep, the pressure of the system is calculated; if the pressure is too high, then the simulation box is made a little larger, while if the pressure is too low the box is made smaller. Simulations in which the pressure is controlled are thus in the NpT ensemble — the volume of the simulation box is not constant!

===Examining the Input Script===

Open one of your input scripts (it doesn't matter which), and look at the section '''### BRING SYSTEM TO REQUIRED STATE ###'''. The line <pre>fix npt all npt temp ${T} ${T} ${tdamp} iso ${p} ${p} ${pdamp}</pre> is the one responsible for switching on the temperature and pressure control. LAMMPS actually allows us to heat or cool the system over the course of a simulation, if we want to — this is the reason that the temperature appears twice in this line. The first ${T} is the desired starting temperature, and the second is the desired temperature at the end of the simulation. We want a constant average temperature, so we specify the same value twice. The same goes for the pressure.

Now look at the lines near the end of the file:

<pre>
### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp press density
variable dens equal density
variable dens2 equal density*density
variable temp equal temp
variable temp2 equal temp*temp
variable press equal press
variable press2 equal press*press
fix aves all ave/time 100 1000 30000 v_dens v_temp v_press v_dens2 v_temp2 v_press2
run 30000
</pre>

The first command, ''thermo_style'', controls which thermodynamic properties are recorded, as before. The next lines are used to measure ''average'' thermodynamic properties for the system. To draw our equations of state, we need to know the average temperature, pressure, and density, and the statistical errors in those quantities. The six variable lines link those quantities (and their squared values, needed for the errors), to variable names that we can use in the averaging command, which is the line starting ''fix aves...''. This command takes a number of input values and averages them every so many timesteps. Exactly how often this happens depends in the values of the three numbers which follow ''ave/time''.

'''<big>This is the fourth section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Equilibration|Equilibration]], or jump ahead to the next section, [[Third year simulation experiment/Structural properties and the radial distribution function|Structural properties and the radial distribution function]].</big>'''

Third year simulation experiment/Equilibration

2021-10-20T13:17:00Z

Ac5116: /* Checking equilibration */

'''<big>This is the third section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Introduction_to_molecular_dynamics_simulation|Introduction to molecular dynamics simulation]], or jump ahead to the next section, [[Third year simulation experiment/Running simulations under specific conditions|Running simulations under specific conditions]].</big>'''
 
We will be using the LAMMPS program to carry out our molecular dynamics simulations.
'''In several places in this section, we will ask you to consult the LAMMPS manual to find out things about how the software works. You can find the manual [https://lammps.sandia.gov/doc/Manual.html here].''' We appreciate that the format of this document can make it a little hard to navigate, but it is the definitive resource on how different commands in LAMMPS work, and is therefore invaluable. The files you will need for this section can be found in the intro folder downloaded previously.

===Creating the simulation box===
In the previous section, it was pointed out that before we can start a simulation, we need to know the initial states of all of the atoms in the system. Exactly what information we need about each atom depends on which method of numerical integration we need, but at the very least we need to specify the starting position of each atom. If we wanted to simulate a crystal, this information would be quite easy to come by — we could just look up the crystal structure, and use that to generate coordinates for however many unit cells we wanted. For this purpose, LAMMPS includes a command which generates crystal lattice structures.

Generating coordinates for atoms in a liquid is more difficult. There is no long range order, so we can't use a single point of reference to work out the positions of every other atom like we can in a solid. We could generate a random position for each atom. This would certainly create a disordered structure, but causes larger problems when we try to run the simulation.

Instead, we are going to place the atoms on the lattice points of a simple cubic lattice. This, of course, is not a situation in which the system is likely to be found physically. It turns out, though, that if we simulate for enough time we will find that the atoms rearrange themselves into more realistic configurations. We will discuss towards the end of this section exactly what is meant by "enough time"!

Consider the line in the input file <pre>lattice sc 0.8</pre> This command [https://lammps.sandia.gov/doc/lattice.html (further info)] creates a grid of points forming a simple cubic lattice (one lattice point per unit cell). The parameter <math>0.8</math> specifies the number density (number of lattice points per unit volume). In a corresponding output file, you will see the line <pre>Lattice spacing in x,y,z = 1.07722 1.07722 1.07722</pre> This indicates that the distance between the points of this lattice is <math>1.07722</math> (in reduced units, remember!).

The next lines in the input file are
<pre>
region box block 0 5 0 5 0 5
create_box 1 box
</pre>
The corresponding log file output is
<pre>
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (5.3860867 5.3860867 5.3860867)
</pre>

The region command [https://lammps.sandia.gov/doc/region.html (further info)] simply defines a geometrical region in space, which we call "box". In this case, "box" is a cube extending ten lattice spacings from the origin in all three dimensions. The subsequent create_box command [https://lammps.sandia.gov/doc/create_box.html (further info)] tells LAMMPS to use the geometrical region called "box" as a template for the simulation box. The number 1 between "create_box" and "box" indicates that our simulation will contain only one type (species) of atom.

So far we have defined a simulation box which is based around a virtual simple cubic lattice. Our box contains 125 (5x5x5) unit cells of this lattice, and so contains 125 lattice points. We now need to fill our simulation box with atoms. The input command is <pre>create_atoms 1 box</pre> [https://lammps.sandia.gov/doc/create_atoms.html (further info)] while the log file simply contains an acknowledgement of this <pre>Created 125 atoms</pre>

The create_atoms command has two arguments; the first tells LAMMPS that all of the atoms that we create will be of type 1. Every atom in the simulation has a type — because we will be simulating a pure fluid, containing only one chemical species, every atom will have the same type. The actual type that we assign to each atom is arbitrary — type 1 does not, for example, need to correspond to the element with atomic number 1 (hydrogen). If we wanted to simulate water, we might make the hydrogen atoms type 1 and the oxygen atoms type 2. We will specify the physical and chemical properties of each atom type later in the input script.

The remaining data in the log file isn't very instructive as it stands — it simply contains a list of the thermodynamic properties of the simulation at certain intervals. In a few sections time, we will plot this data, but for now you can close the log file. Keep the input script open.

===Setting the properties of the atoms===

In addition to their positions, we also need the physical properties of the atoms to be able to perform the simulation. We set these properties on a 'per-type' basis, so that every atom of the same type has the same mass and the same interactions.

So far we have created 125 atoms, and we know the starting (<math>t = 0</math>) position for each of them. We have also set their masses, and told LAMMPS what sort of forces to calculate between them. The final thing we need to specify to completely specify the initial conditions is the velocity of each atom.

Choosing initial velocities for the atoms is a little easier than choosing initial positions. From the statistical thermodynamics lectures, you should know that, at equilibrium, the velocities of atoms in any system must be distributed according to the [http://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution Maxwell-Boltzmann (MB) distribution]. If we know the masses of the atoms, and we know what temperature we want to simulate, then we can determine the relevant MB distribution function. LAMMPS is able to give every atom a random velocity whilst ensuring that overall the MB distribution is followed. This is the purpose of the line

<pre>
velocity all create 1.5 12345 dist gaussian rot yes mom yes
</pre>

You can see the manual page for this command [http://lammps.sandia.gov/doc/velocity.html here], but the key sections are:

* '''all''': the ''group'' of atoms on which the command acts. '''all''' simply specifies that we want every atom to have a velocity assigned to it.
* '''1.5''': the temperature, <math>T</math>, needed to calculate the MB distribution(in reduced units, as always)

===Monitoring thermodynamic properties===

We need to be sure that our simulation is correctly modelling whatever physical system we want to study. It is relatively easy to set up simulations, but how can we be sure that the "results" we get make sense? One of the best ways is to calculate from the simulation things that we can measure in experiment, and see if they agree. For example, we might want to simulate our system at a particular temperature and pressure, and measure the resulting density. If we repeat this over a range of temperatures at the same pressure, we will be able to plot an ''equation of state'', which we could compare to experimental measurements.

LAMMPS is able to calculate a great deal of thermodynamic information for us (you can see a full list of the properties it is able to calculate [http://lammps.sandia.gov/doc/thermo_style.html here]), but in these first simulations we are only interested in those properties specified in these commands:

<pre>
thermo_style custom time etotal temp press
thermo 10
</pre>

The first controls which properties will be printed out in the log file. In this case, we print how much time we have simulated so far (which is ''not'' the same as how long it has taken us to simulate it!), the total energy of the atoms, their temperature, and their pressure. The second line tells LAMMPS to print this information on every 10th timestep.

===Running the simulation===

'''Look at the lines below.'''
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(50/${timestep})
timestep ${timestep}

### RUN SIMULATION ###
run ${n_steps}
</pre>
'''The second line (starting "variable timestep...") tells LAMMPS that if it encounters the text ${timestep} on a subsequent line, it should replace it by the value given. In this case, the value ${timestep} is always replaced by 0.001.'''
 

'''<big> It is now time to run your first simulation, submit the input script with the data file in the intro folder of the files you have downloaded Try changing the timestep - what happens when you make the timestep larger?. </big>'''

===Visualising the trajectory===

The '''trajectory files''' contain the positions of all the atoms in the simulation, recorded at a set interval (for all of these simulations, this was every ten timesteps — this is controlled by the '''dump''' command in the input scripts). We use a programme called [http://www.ks.uiuc.edu/Research/vmd/ '''VMD'''] to view these trajectories, which you should find is already installed on both the desktop and laptop computers. You can run VMD from the start menu with '''Start''' -> '''All Programs''' -> '''University of Illinois''' -> '''VMD'''.

====Loading a Trajectory====

We'll start by looking at the output of the 0.02 timestep simulation. In the '''VMD Main''' window, select the menu option '''File''' -> '''New Molecule'''. Click the '''Browse''' button, then select the relevant trajectory file. In the '''Determine file type''' dropdown, select '''LAMMPS Trajectory'''. Then click '''Load'''.

You will see that the '''VMD 1.9.1 OpenGL Display''' window now shows a horrible mess. VMD's default behaviour is to draw lines between atoms which it thinks might be chemically bonded. Our system doesn't model chemical bonds, so we want to turn this off. In the '''VMD Main''' window, select the menu option '''Graphics''' -> '''Representations'''. This shows a list of "representations" of our atoms. You will see that at the moment, there is a single representation listed, and it is selected. It will have the ''Lines'' style, the ''Name'' colour, and the selection ''all''. "Selection" simply tells VMD which atoms we want it to draw. We want to show every atom, so the current selection is fine. The ''name'' colouring method just makes VMD give atoms colours according to their specified type. The colour isn't important to us, so we can leave this be too. The "style" tells VMD what we want it to display for each atom. Change the '''Drawing Method''' from ''Lines'' to ''VDW''. You will see that the mess of lines is replaced by a mess of low resolution, overlapping spheres. Change the '''Sphere Scale''' to 0.3, and the '''Sphere Resolution''' to 17. The result should look a little smoother. Close the '''Graphical Representations''' window. You will notice that in the bottom right of the '''VMD Main''' window, there is a small play button. Click this, and you will see the animated version of your simulation trajectory.

By clicking and dragging with the mouse, you can rotate the simulation box (though this may be sluggish). At any time, you can reset the view by pressing the equals key.

====Tracking a Single Particle====
To illustrate the periodic boundary conditions that we are using, we are going to draw almost all of the atoms as points, but we will pick a single atom at random to draw as a sphere. This will make it easy to see how a single atom moves through the box. Reset the display using the equals key, then use the '''VMD Main''' window controls to pause the trajectory and reset it to the first trajectory (play with the different buttons until you find the one that does this). You should see the perfect cubic lattice. Use the option '''Display''' -> '''Orthographic''' to change the drawing mode, then rotate the displayed crystal so that you are looking at one vertex (looking down the 111 direction, in crystallographic terms).

Open the '''Graphical Representations''' window again. Change the representation style from '''VDW''' to points, then click the '''Create Rep''' button. This creates a second representation, allowing a subset of the atoms to be drawn in a different way. The '''Selected Atoms''' box allows us to choose which atoms this representation applies to. We just want to pick two of them at random — VMD assigns every atom an index, from 0 to N-1. In our case, there are 125 atoms, so choose two numbers between 0 and 124. Changed the '''Selected Atoms''' field to <pre>index i or index j</pre> where i and j are your chosen numbers, press return, then change the '''Drawing Method''' to '''VDW'''. You should now see only two atoms represented by spheres, with the rest shown as small points. In the '''VMD Main''' window, click play. Try rotating the box, and changing the playback speed.

You will see that sometimes one of the spheres seems to change position across the box very rapidly — this occurs when it reaches one periodic boundary, and is reflected back across the other face. Try playing with some of the other representation types in VMD — it is a very powerful package, which is often used to render images of simulated proteins, so many of its options aren't relevant to our simple system!

===Checking equilibration===

When we first set up a simulation, it is very important to make sure that our system reaches an equilibrium state. We characterise equilibrium by the average values of thermodynamic quantities becoming constant (due to the approximations that we have made, there will always be fluctuations, but the average values will become constant).

In this section, we are going to plot the thermodynamic output of the simulation to see how long it takes to reach the equilibrium state (and indeed, whether this happens at all). Instructions are given below to import data from the LAMMPS log file into Microsoft Excel. Once you have the data in a spreadsheet, you can plot it. If you know how to use some of the other plotting software available on the chemistry computers (like Origin), you are welcome to use it.

# Open a blank Excel workbook
# Copy the data in the textfile into the first cell
# With these data highlighted, click the Data tab and "Text to Columns"
# Click "Delimited", continue and let it be space delimited
# Click finish
You can then export this and data as a .csv and analyse in Python or Excel, as you wish.

Challenge: Can you write a python script or function that extracts this data for a given file automatically?

'''<big>TASK 7: </big>'''

=== What does it mean for a simulation to "reach equilibrium"? Why is this important in terms of sampling from an ensemble using molecular dynamics? [2] ===

=== Plot the energy (potential, kinetic and total), temperature and pressure, against time for the 0.001 timestep experiment [2] ===

=== Does the simulation reach equilibrium? How can you tell? [2] ===

=== Make a single plot which shows the energy vs. time for the timesteps you have simulated [2]. ===

=== Of the timesteps that you used, which timestep will you use for subsequent simulations and why? [6] ===
''(Think about what is happening "physically" as you increase/decrease the timestep. Also, what features of each timeseries are indicative of the simulation's "health"?)''

'''<big>This is the third section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Introduction_to_molecular_dynamics_simulation|Introduction to molecular dynamics simulation]], or jump ahead to the next section, [[Third year simulation experiment/Running simulations under specific conditions|Running simulations under specific conditions]].</big>'''

Third year simulation experiment/Equilibration

2021-10-20T13:14:42Z

Ac5116: /* Plot the energy (potential, kinetic and total), temperature and pressure, against time for the 0,001 timestep experiment [2] */

'''<big>This is the third section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Introduction_to_molecular_dynamics_simulation|Introduction to molecular dynamics simulation]], or jump ahead to the next section, [[Third year simulation experiment/Running simulations under specific conditions|Running simulations under specific conditions]].</big>'''
 
We will be using the LAMMPS program to carry out our molecular dynamics simulations.
'''In several places in this section, we will ask you to consult the LAMMPS manual to find out things about how the software works. You can find the manual [https://lammps.sandia.gov/doc/Manual.html here].''' We appreciate that the format of this document can make it a little hard to navigate, but it is the definitive resource on how different commands in LAMMPS work, and is therefore invaluable. The files you will need for this section can be found in the intro folder downloaded previously.

===Creating the simulation box===
In the previous section, it was pointed out that before we can start a simulation, we need to know the initial states of all of the atoms in the system. Exactly what information we need about each atom depends on which method of numerical integration we need, but at the very least we need to specify the starting position of each atom. If we wanted to simulate a crystal, this information would be quite easy to come by — we could just look up the crystal structure, and use that to generate coordinates for however many unit cells we wanted. For this purpose, LAMMPS includes a command which generates crystal lattice structures.

Generating coordinates for atoms in a liquid is more difficult. There is no long range order, so we can't use a single point of reference to work out the positions of every other atom like we can in a solid. We could generate a random position for each atom. This would certainly create a disordered structure, but causes larger problems when we try to run the simulation.

Instead, we are going to place the atoms on the lattice points of a simple cubic lattice. This, of course, is not a situation in which the system is likely to be found physically. It turns out, though, that if we simulate for enough time we will find that the atoms rearrange themselves into more realistic configurations. We will discuss towards the end of this section exactly what is meant by "enough time"!

Consider the line in the input file <pre>lattice sc 0.8</pre> This command [https://lammps.sandia.gov/doc/lattice.html (further info)] creates a grid of points forming a simple cubic lattice (one lattice point per unit cell). The parameter <math>0.8</math> specifies the number density (number of lattice points per unit volume). In a corresponding output file, you will see the line <pre>Lattice spacing in x,y,z = 1.07722 1.07722 1.07722</pre> This indicates that the distance between the points of this lattice is <math>1.07722</math> (in reduced units, remember!).

The next lines in the input file are
<pre>
region box block 0 5 0 5 0 5
create_box 1 box
</pre>
The corresponding log file output is
<pre>
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (5.3860867 5.3860867 5.3860867)
</pre>

The region command [https://lammps.sandia.gov/doc/region.html (further info)] simply defines a geometrical region in space, which we call "box". In this case, "box" is a cube extending ten lattice spacings from the origin in all three dimensions. The subsequent create_box command [https://lammps.sandia.gov/doc/create_box.html (further info)] tells LAMMPS to use the geometrical region called "box" as a template for the simulation box. The number 1 between "create_box" and "box" indicates that our simulation will contain only one type (species) of atom.

So far we have defined a simulation box which is based around a virtual simple cubic lattice. Our box contains 125 (5x5x5) unit cells of this lattice, and so contains 125 lattice points. We now need to fill our simulation box with atoms. The input command is <pre>create_atoms 1 box</pre> [https://lammps.sandia.gov/doc/create_atoms.html (further info)] while the log file simply contains an acknowledgement of this <pre>Created 125 atoms</pre>

The create_atoms command has two arguments; the first tells LAMMPS that all of the atoms that we create will be of type 1. Every atom in the simulation has a type — because we will be simulating a pure fluid, containing only one chemical species, every atom will have the same type. The actual type that we assign to each atom is arbitrary — type 1 does not, for example, need to correspond to the element with atomic number 1 (hydrogen). If we wanted to simulate water, we might make the hydrogen atoms type 1 and the oxygen atoms type 2. We will specify the physical and chemical properties of each atom type later in the input script.

The remaining data in the log file isn't very instructive as it stands — it simply contains a list of the thermodynamic properties of the simulation at certain intervals. In a few sections time, we will plot this data, but for now you can close the log file. Keep the input script open.

===Setting the properties of the atoms===

In addition to their positions, we also need the physical properties of the atoms to be able to perform the simulation. We set these properties on a 'per-type' basis, so that every atom of the same type has the same mass and the same interactions.

So far we have created 125 atoms, and we know the starting (<math>t = 0</math>) position for each of them. We have also set their masses, and told LAMMPS what sort of forces to calculate between them. The final thing we need to specify to completely specify the initial conditions is the velocity of each atom.

Choosing initial velocities for the atoms is a little easier than choosing initial positions. From the statistical thermodynamics lectures, you should know that, at equilibrium, the velocities of atoms in any system must be distributed according to the [http://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution Maxwell-Boltzmann (MB) distribution]. If we know the masses of the atoms, and we know what temperature we want to simulate, then we can determine the relevant MB distribution function. LAMMPS is able to give every atom a random velocity whilst ensuring that overall the MB distribution is followed. This is the purpose of the line

<pre>
velocity all create 1.5 12345 dist gaussian rot yes mom yes
</pre>

You can see the manual page for this command [http://lammps.sandia.gov/doc/velocity.html here], but the key sections are:

* '''all''': the ''group'' of atoms on which the command acts. '''all''' simply specifies that we want every atom to have a velocity assigned to it.
* '''1.5''': the temperature, <math>T</math>, needed to calculate the MB distribution(in reduced units, as always)

===Monitoring thermodynamic properties===

We need to be sure that our simulation is correctly modelling whatever physical system we want to study. It is relatively easy to set up simulations, but how can we be sure that the "results" we get make sense? One of the best ways is to calculate from the simulation things that we can measure in experiment, and see if they agree. For example, we might want to simulate our system at a particular temperature and pressure, and measure the resulting density. If we repeat this over a range of temperatures at the same pressure, we will be able to plot an ''equation of state'', which we could compare to experimental measurements.

LAMMPS is able to calculate a great deal of thermodynamic information for us (you can see a full list of the properties it is able to calculate [http://lammps.sandia.gov/doc/thermo_style.html here]), but in these first simulations we are only interested in those properties specified in these commands:

<pre>
thermo_style custom time etotal temp press
thermo 10
</pre>

The first controls which properties will be printed out in the log file. In this case, we print how much time we have simulated so far (which is ''not'' the same as how long it has taken us to simulate it!), the total energy of the atoms, their temperature, and their pressure. The second line tells LAMMPS to print this information on every 10th timestep.

===Running the simulation===

'''Look at the lines below.'''
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(50/${timestep})
timestep ${timestep}

### RUN SIMULATION ###
run ${n_steps}
</pre>
'''The second line (starting "variable timestep...") tells LAMMPS that if it encounters the text ${timestep} on a subsequent line, it should replace it by the value given. In this case, the value ${timestep} is always replaced by 0.001.'''
 

'''<big> It is now time to run your first simulation, submit the input script with the data file in the intro folder of the files you have downloaded Try changing the timestep - what happens when you make the timestep larger?. </big>'''

===Visualising the trajectory===

The '''trajectory files''' contain the positions of all the atoms in the simulation, recorded at a set interval (for all of these simulations, this was every ten timesteps — this is controlled by the '''dump''' command in the input scripts). We use a programme called [http://www.ks.uiuc.edu/Research/vmd/ '''VMD'''] to view these trajectories, which you should find is already installed on both the desktop and laptop computers. You can run VMD from the start menu with '''Start''' -> '''All Programs''' -> '''University of Illinois''' -> '''VMD'''.

====Loading a Trajectory====

We'll start by looking at the output of the 0.02 timestep simulation. In the '''VMD Main''' window, select the menu option '''File''' -> '''New Molecule'''. Click the '''Browse''' button, then select the relevant trajectory file. In the '''Determine file type''' dropdown, select '''LAMMPS Trajectory'''. Then click '''Load'''.

You will see that the '''VMD 1.9.1 OpenGL Display''' window now shows a horrible mess. VMD's default behaviour is to draw lines between atoms which it thinks might be chemically bonded. Our system doesn't model chemical bonds, so we want to turn this off. In the '''VMD Main''' window, select the menu option '''Graphics''' -> '''Representations'''. This shows a list of "representations" of our atoms. You will see that at the moment, there is a single representation listed, and it is selected. It will have the ''Lines'' style, the ''Name'' colour, and the selection ''all''. "Selection" simply tells VMD which atoms we want it to draw. We want to show every atom, so the current selection is fine. The ''name'' colouring method just makes VMD give atoms colours according to their specified type. The colour isn't important to us, so we can leave this be too. The "style" tells VMD what we want it to display for each atom. Change the '''Drawing Method''' from ''Lines'' to ''VDW''. You will see that the mess of lines is replaced by a mess of low resolution, overlapping spheres. Change the '''Sphere Scale''' to 0.3, and the '''Sphere Resolution''' to 17. The result should look a little smoother. Close the '''Graphical Representations''' window. You will notice that in the bottom right of the '''VMD Main''' window, there is a small play button. Click this, and you will see the animated version of your simulation trajectory.

By clicking and dragging with the mouse, you can rotate the simulation box (though this may be sluggish). At any time, you can reset the view by pressing the equals key.

====Tracking a Single Particle====
To illustrate the periodic boundary conditions that we are using, we are going to draw almost all of the atoms as points, but we will pick a single atom at random to draw as a sphere. This will make it easy to see how a single atom moves through the box. Reset the display using the equals key, then use the '''VMD Main''' window controls to pause the trajectory and reset it to the first trajectory (play with the different buttons until you find the one that does this). You should see the perfect cubic lattice. Use the option '''Display''' -> '''Orthographic''' to change the drawing mode, then rotate the displayed crystal so that you are looking at one vertex (looking down the 111 direction, in crystallographic terms).

Open the '''Graphical Representations''' window again. Change the representation style from '''VDW''' to points, then click the '''Create Rep''' button. This creates a second representation, allowing a subset of the atoms to be drawn in a different way. The '''Selected Atoms''' box allows us to choose which atoms this representation applies to. We just want to pick two of them at random — VMD assigns every atom an index, from 0 to N-1. In our case, there are 125 atoms, so choose two numbers between 0 and 124. Changed the '''Selected Atoms''' field to <pre>index i or index j</pre> where i and j are your chosen numbers, press return, then change the '''Drawing Method''' to '''VDW'''. You should now see only two atoms represented by spheres, with the rest shown as small points. In the '''VMD Main''' window, click play. Try rotating the box, and changing the playback speed.

You will see that sometimes one of the spheres seems to change position across the box very rapidly — this occurs when it reaches one periodic boundary, and is reflected back across the other face. Try playing with some of the other representation types in VMD — it is a very powerful package, which is often used to render images of simulated proteins, so many of its options aren't relevant to our simple system!

===Checking equilibration===

When we first set up a simulation, it is very important to make sure that our system reaches an equilibrium state. We characterise equilibrium by the average values of thermodynamic quantities becoming constant (due to the approximations that we have made, there will always be fluctuations, but the average values will become constant).

In this section, we are going to plot the thermodynamic output of the simulation to see how long it takes to reach the equilibrium state (and indeed, whether this happens at all). Instructions are given below to import data from the LAMMPS log file into Microsoft Excel. Once you have the data in a spreadsheet, you can plot it. If you know how to use some of the other plotting software available on the chemistry computers (like Origin), you are welcome to use it.

# Open a blank Excel workbook
# Copy the data in the textfile into the first cell
# With these data highlighted, click the Data tab and "Text to Columns"
# Click "Delimited", continue and let it be space delimited
# Click finish

'''<big>TASK 7: </big>'''

=== What does it mean for a simulation to "reach equilibrium"? Why is this important in terms of sampling from an ensemble using molecular dynamics? [2] ===

=== Plot the energy (potential, kinetic and total), temperature and pressure, against time for the 0.001 timestep experiment [2] ===

=== Does the simulation reach equilibrium? How can you tell? [2] ===

=== Make a single plot which shows the energy vs. time for the timesteps you have simulated [2]. ===

=== Of the timesteps that you used, which timestep will you use for subsequent simulations and why? [6] ===
''(Think about what is happening "physically" as you increase/decrease the timestep. Also, what features of each timeseries are indicative of the simulation's "health"?)''

'''<big>This is the third section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Introduction_to_molecular_dynamics_simulation|Introduction to molecular dynamics simulation]], or jump ahead to the next section, [[Third year simulation experiment/Running simulations under specific conditions|Running simulations under specific conditions]].</big>'''

Third year simulation experiment/Equilibration

2021-10-20T13:10:35Z

Ac5116: /* Running the simulation */

'''<big>This is the third section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Introduction_to_molecular_dynamics_simulation|Introduction to molecular dynamics simulation]], or jump ahead to the next section, [[Third year simulation experiment/Running simulations under specific conditions|Running simulations under specific conditions]].</big>'''
 
We will be using the LAMMPS program to carry out our molecular dynamics simulations.
'''In several places in this section, we will ask you to consult the LAMMPS manual to find out things about how the software works. You can find the manual [https://lammps.sandia.gov/doc/Manual.html here].''' We appreciate that the format of this document can make it a little hard to navigate, but it is the definitive resource on how different commands in LAMMPS work, and is therefore invaluable. The files you will need for this section can be found in the intro folder downloaded previously.

===Creating the simulation box===
In the previous section, it was pointed out that before we can start a simulation, we need to know the initial states of all of the atoms in the system. Exactly what information we need about each atom depends on which method of numerical integration we need, but at the very least we need to specify the starting position of each atom. If we wanted to simulate a crystal, this information would be quite easy to come by — we could just look up the crystal structure, and use that to generate coordinates for however many unit cells we wanted. For this purpose, LAMMPS includes a command which generates crystal lattice structures.

Generating coordinates for atoms in a liquid is more difficult. There is no long range order, so we can't use a single point of reference to work out the positions of every other atom like we can in a solid. We could generate a random position for each atom. This would certainly create a disordered structure, but causes larger problems when we try to run the simulation.

Instead, we are going to place the atoms on the lattice points of a simple cubic lattice. This, of course, is not a situation in which the system is likely to be found physically. It turns out, though, that if we simulate for enough time we will find that the atoms rearrange themselves into more realistic configurations. We will discuss towards the end of this section exactly what is meant by "enough time"!

Consider the line in the input file <pre>lattice sc 0.8</pre> This command [https://lammps.sandia.gov/doc/lattice.html (further info)] creates a grid of points forming a simple cubic lattice (one lattice point per unit cell). The parameter <math>0.8</math> specifies the number density (number of lattice points per unit volume). In a corresponding output file, you will see the line <pre>Lattice spacing in x,y,z = 1.07722 1.07722 1.07722</pre> This indicates that the distance between the points of this lattice is <math>1.07722</math> (in reduced units, remember!).

The next lines in the input file are
<pre>
region box block 0 5 0 5 0 5
create_box 1 box
</pre>
The corresponding log file output is
<pre>
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (5.3860867 5.3860867 5.3860867)
</pre>

The region command [https://lammps.sandia.gov/doc/region.html (further info)] simply defines a geometrical region in space, which we call "box". In this case, "box" is a cube extending ten lattice spacings from the origin in all three dimensions. The subsequent create_box command [https://lammps.sandia.gov/doc/create_box.html (further info)] tells LAMMPS to use the geometrical region called "box" as a template for the simulation box. The number 1 between "create_box" and "box" indicates that our simulation will contain only one type (species) of atom.

So far we have defined a simulation box which is based around a virtual simple cubic lattice. Our box contains 125 (5x5x5) unit cells of this lattice, and so contains 125 lattice points. We now need to fill our simulation box with atoms. The input command is <pre>create_atoms 1 box</pre> [https://lammps.sandia.gov/doc/create_atoms.html (further info)] while the log file simply contains an acknowledgement of this <pre>Created 125 atoms</pre>

The create_atoms command has two arguments; the first tells LAMMPS that all of the atoms that we create will be of type 1. Every atom in the simulation has a type — because we will be simulating a pure fluid, containing only one chemical species, every atom will have the same type. The actual type that we assign to each atom is arbitrary — type 1 does not, for example, need to correspond to the element with atomic number 1 (hydrogen). If we wanted to simulate water, we might make the hydrogen atoms type 1 and the oxygen atoms type 2. We will specify the physical and chemical properties of each atom type later in the input script.

The remaining data in the log file isn't very instructive as it stands — it simply contains a list of the thermodynamic properties of the simulation at certain intervals. In a few sections time, we will plot this data, but for now you can close the log file. Keep the input script open.

===Setting the properties of the atoms===

In addition to their positions, we also need the physical properties of the atoms to be able to perform the simulation. We set these properties on a 'per-type' basis, so that every atom of the same type has the same mass and the same interactions.

So far we have created 125 atoms, and we know the starting (<math>t = 0</math>) position for each of them. We have also set their masses, and told LAMMPS what sort of forces to calculate between them. The final thing we need to specify to completely specify the initial conditions is the velocity of each atom.

Choosing initial velocities for the atoms is a little easier than choosing initial positions. From the statistical thermodynamics lectures, you should know that, at equilibrium, the velocities of atoms in any system must be distributed according to the [http://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution Maxwell-Boltzmann (MB) distribution]. If we know the masses of the atoms, and we know what temperature we want to simulate, then we can determine the relevant MB distribution function. LAMMPS is able to give every atom a random velocity whilst ensuring that overall the MB distribution is followed. This is the purpose of the line

<pre>
velocity all create 1.5 12345 dist gaussian rot yes mom yes
</pre>

You can see the manual page for this command [http://lammps.sandia.gov/doc/velocity.html here], but the key sections are:

* '''all''': the ''group'' of atoms on which the command acts. '''all''' simply specifies that we want every atom to have a velocity assigned to it.
* '''1.5''': the temperature, <math>T</math>, needed to calculate the MB distribution(in reduced units, as always)

===Monitoring thermodynamic properties===

We need to be sure that our simulation is correctly modelling whatever physical system we want to study. It is relatively easy to set up simulations, but how can we be sure that the "results" we get make sense? One of the best ways is to calculate from the simulation things that we can measure in experiment, and see if they agree. For example, we might want to simulate our system at a particular temperature and pressure, and measure the resulting density. If we repeat this over a range of temperatures at the same pressure, we will be able to plot an ''equation of state'', which we could compare to experimental measurements.

LAMMPS is able to calculate a great deal of thermodynamic information for us (you can see a full list of the properties it is able to calculate [http://lammps.sandia.gov/doc/thermo_style.html here]), but in these first simulations we are only interested in those properties specified in these commands:

<pre>
thermo_style custom time etotal temp press
thermo 10
</pre>

The first controls which properties will be printed out in the log file. In this case, we print how much time we have simulated so far (which is ''not'' the same as how long it has taken us to simulate it!), the total energy of the atoms, their temperature, and their pressure. The second line tells LAMMPS to print this information on every 10th timestep.

===Running the simulation===

'''Look at the lines below.'''
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(50/${timestep})
timestep ${timestep}

### RUN SIMULATION ###
run ${n_steps}
</pre>
'''The second line (starting "variable timestep...") tells LAMMPS that if it encounters the text ${timestep} on a subsequent line, it should replace it by the value given. In this case, the value ${timestep} is always replaced by 0.001.'''
 

'''<big> It is now time to run your first simulation, submit the input script with the data file in the intro folder of the files you have downloaded Try changing the timestep - what happens when you make the timestep larger?. </big>'''

===Visualising the trajectory===

The '''trajectory files''' contain the positions of all the atoms in the simulation, recorded at a set interval (for all of these simulations, this was every ten timesteps — this is controlled by the '''dump''' command in the input scripts). We use a programme called [http://www.ks.uiuc.edu/Research/vmd/ '''VMD'''] to view these trajectories, which you should find is already installed on both the desktop and laptop computers. You can run VMD from the start menu with '''Start''' -> '''All Programs''' -> '''University of Illinois''' -> '''VMD'''.

====Loading a Trajectory====

We'll start by looking at the output of the 0.02 timestep simulation. In the '''VMD Main''' window, select the menu option '''File''' -> '''New Molecule'''. Click the '''Browse''' button, then select the relevant trajectory file. In the '''Determine file type''' dropdown, select '''LAMMPS Trajectory'''. Then click '''Load'''.

You will see that the '''VMD 1.9.1 OpenGL Display''' window now shows a horrible mess. VMD's default behaviour is to draw lines between atoms which it thinks might be chemically bonded. Our system doesn't model chemical bonds, so we want to turn this off. In the '''VMD Main''' window, select the menu option '''Graphics''' -> '''Representations'''. This shows a list of "representations" of our atoms. You will see that at the moment, there is a single representation listed, and it is selected. It will have the ''Lines'' style, the ''Name'' colour, and the selection ''all''. "Selection" simply tells VMD which atoms we want it to draw. We want to show every atom, so the current selection is fine. The ''name'' colouring method just makes VMD give atoms colours according to their specified type. The colour isn't important to us, so we can leave this be too. The "style" tells VMD what we want it to display for each atom. Change the '''Drawing Method''' from ''Lines'' to ''VDW''. You will see that the mess of lines is replaced by a mess of low resolution, overlapping spheres. Change the '''Sphere Scale''' to 0.3, and the '''Sphere Resolution''' to 17. The result should look a little smoother. Close the '''Graphical Representations''' window. You will notice that in the bottom right of the '''VMD Main''' window, there is a small play button. Click this, and you will see the animated version of your simulation trajectory.

By clicking and dragging with the mouse, you can rotate the simulation box (though this may be sluggish). At any time, you can reset the view by pressing the equals key.

====Tracking a Single Particle====
To illustrate the periodic boundary conditions that we are using, we are going to draw almost all of the atoms as points, but we will pick a single atom at random to draw as a sphere. This will make it easy to see how a single atom moves through the box. Reset the display using the equals key, then use the '''VMD Main''' window controls to pause the trajectory and reset it to the first trajectory (play with the different buttons until you find the one that does this). You should see the perfect cubic lattice. Use the option '''Display''' -> '''Orthographic''' to change the drawing mode, then rotate the displayed crystal so that you are looking at one vertex (looking down the 111 direction, in crystallographic terms).

Open the '''Graphical Representations''' window again. Change the representation style from '''VDW''' to points, then click the '''Create Rep''' button. This creates a second representation, allowing a subset of the atoms to be drawn in a different way. The '''Selected Atoms''' box allows us to choose which atoms this representation applies to. We just want to pick two of them at random — VMD assigns every atom an index, from 0 to N-1. In our case, there are 125 atoms, so choose two numbers between 0 and 124. Changed the '''Selected Atoms''' field to <pre>index i or index j</pre> where i and j are your chosen numbers, press return, then change the '''Drawing Method''' to '''VDW'''. You should now see only two atoms represented by spheres, with the rest shown as small points. In the '''VMD Main''' window, click play. Try rotating the box, and changing the playback speed.

You will see that sometimes one of the spheres seems to change position across the box very rapidly — this occurs when it reaches one periodic boundary, and is reflected back across the other face. Try playing with some of the other representation types in VMD — it is a very powerful package, which is often used to render images of simulated proteins, so many of its options aren't relevant to our simple system!

===Checking equilibration===

When we first set up a simulation, it is very important to make sure that our system reaches an equilibrium state. We characterise equilibrium by the average values of thermodynamic quantities becoming constant (due to the approximations that we have made, there will always be fluctuations, but the average values will become constant).

In this section, we are going to plot the thermodynamic output of the simulation to see how long it takes to reach the equilibrium state (and indeed, whether this happens at all). Instructions are given below to import data from the LAMMPS log file into Microsoft Excel. Once you have the data in a spreadsheet, you can plot it. If you know how to use some of the other plotting software available on the chemistry computers (like Origin), you are welcome to use it.

# Open a blank Excel workbook
# Copy the data in the textfile into the first cell
# With these data highlighted, click the Data tab and "Text to Columns"
# Click "Delimited", continue and let it be space delimited
# Click finish

'''<big>TASK 7: </big>'''

=== What does it mean for a simulation to "reach equilibrium"? Why is this important in terms of sampling from an ensemble using molecular dynamics? [2] ===

=== Plot the energy (potential, kinetic and total), temperature and pressure, against time for the 0,001 timestep experiment [2] ===

=== Does the simulation reach equilibrium? How can you tell? [2] ===

=== Make a single plot which shows the energy vs. time for the timesteps you have simulated [2]. ===

=== Of the timesteps that you used, which timestep will you use for subsequent simulations and why? [6] ===
''(Think about what is happening "physically" as you increase/decrease the timestep. Also, what features of each timeseries are indicative of the simulation's "health"?)''

'''<big>This is the third section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Introduction_to_molecular_dynamics_simulation|Introduction to molecular dynamics simulation]], or jump ahead to the next section, [[Third year simulation experiment/Running simulations under specific conditions|Running simulations under specific conditions]].</big>'''

Third year simulation experiment/Equilibration

2021-10-20T13:08:37Z

Ac5116: /* Creating the simulation box */

'''<big>This is the third section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Introduction_to_molecular_dynamics_simulation|Introduction to molecular dynamics simulation]], or jump ahead to the next section, [[Third year simulation experiment/Running simulations under specific conditions|Running simulations under specific conditions]].</big>'''
 
We will be using the LAMMPS program to carry out our molecular dynamics simulations.
'''In several places in this section, we will ask you to consult the LAMMPS manual to find out things about how the software works. You can find the manual [https://lammps.sandia.gov/doc/Manual.html here].''' We appreciate that the format of this document can make it a little hard to navigate, but it is the definitive resource on how different commands in LAMMPS work, and is therefore invaluable. The files you will need for this section can be found in the intro folder downloaded previously.

===Creating the simulation box===
In the previous section, it was pointed out that before we can start a simulation, we need to know the initial states of all of the atoms in the system. Exactly what information we need about each atom depends on which method of numerical integration we need, but at the very least we need to specify the starting position of each atom. If we wanted to simulate a crystal, this information would be quite easy to come by — we could just look up the crystal structure, and use that to generate coordinates for however many unit cells we wanted. For this purpose, LAMMPS includes a command which generates crystal lattice structures.

Generating coordinates for atoms in a liquid is more difficult. There is no long range order, so we can't use a single point of reference to work out the positions of every other atom like we can in a solid. We could generate a random position for each atom. This would certainly create a disordered structure, but causes larger problems when we try to run the simulation.

Instead, we are going to place the atoms on the lattice points of a simple cubic lattice. This, of course, is not a situation in which the system is likely to be found physically. It turns out, though, that if we simulate for enough time we will find that the atoms rearrange themselves into more realistic configurations. We will discuss towards the end of this section exactly what is meant by "enough time"!

Consider the line in the input file <pre>lattice sc 0.8</pre> This command [https://lammps.sandia.gov/doc/lattice.html (further info)] creates a grid of points forming a simple cubic lattice (one lattice point per unit cell). The parameter <math>0.8</math> specifies the number density (number of lattice points per unit volume). In a corresponding output file, you will see the line <pre>Lattice spacing in x,y,z = 1.07722 1.07722 1.07722</pre> This indicates that the distance between the points of this lattice is <math>1.07722</math> (in reduced units, remember!).

The next lines in the input file are
<pre>
region box block 0 5 0 5 0 5
create_box 1 box
</pre>
The corresponding log file output is
<pre>
Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (5.3860867 5.3860867 5.3860867)
</pre>

The region command [https://lammps.sandia.gov/doc/region.html (further info)] simply defines a geometrical region in space, which we call "box". In this case, "box" is a cube extending ten lattice spacings from the origin in all three dimensions. The subsequent create_box command [https://lammps.sandia.gov/doc/create_box.html (further info)] tells LAMMPS to use the geometrical region called "box" as a template for the simulation box. The number 1 between "create_box" and "box" indicates that our simulation will contain only one type (species) of atom.

So far we have defined a simulation box which is based around a virtual simple cubic lattice. Our box contains 125 (5x5x5) unit cells of this lattice, and so contains 125 lattice points. We now need to fill our simulation box with atoms. The input command is <pre>create_atoms 1 box</pre> [https://lammps.sandia.gov/doc/create_atoms.html (further info)] while the log file simply contains an acknowledgement of this <pre>Created 125 atoms</pre>

The create_atoms command has two arguments; the first tells LAMMPS that all of the atoms that we create will be of type 1. Every atom in the simulation has a type — because we will be simulating a pure fluid, containing only one chemical species, every atom will have the same type. The actual type that we assign to each atom is arbitrary — type 1 does not, for example, need to correspond to the element with atomic number 1 (hydrogen). If we wanted to simulate water, we might make the hydrogen atoms type 1 and the oxygen atoms type 2. We will specify the physical and chemical properties of each atom type later in the input script.

The remaining data in the log file isn't very instructive as it stands — it simply contains a list of the thermodynamic properties of the simulation at certain intervals. In a few sections time, we will plot this data, but for now you can close the log file. Keep the input script open.

===Setting the properties of the atoms===

In addition to their positions, we also need the physical properties of the atoms to be able to perform the simulation. We set these properties on a 'per-type' basis, so that every atom of the same type has the same mass and the same interactions.

So far we have created 125 atoms, and we know the starting (<math>t = 0</math>) position for each of them. We have also set their masses, and told LAMMPS what sort of forces to calculate between them. The final thing we need to specify to completely specify the initial conditions is the velocity of each atom.

Choosing initial velocities for the atoms is a little easier than choosing initial positions. From the statistical thermodynamics lectures, you should know that, at equilibrium, the velocities of atoms in any system must be distributed according to the [http://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution Maxwell-Boltzmann (MB) distribution]. If we know the masses of the atoms, and we know what temperature we want to simulate, then we can determine the relevant MB distribution function. LAMMPS is able to give every atom a random velocity whilst ensuring that overall the MB distribution is followed. This is the purpose of the line

<pre>
velocity all create 1.5 12345 dist gaussian rot yes mom yes
</pre>

You can see the manual page for this command [http://lammps.sandia.gov/doc/velocity.html here], but the key sections are:

* '''all''': the ''group'' of atoms on which the command acts. '''all''' simply specifies that we want every atom to have a velocity assigned to it.
* '''1.5''': the temperature, <math>T</math>, needed to calculate the MB distribution(in reduced units, as always)

===Monitoring thermodynamic properties===

We need to be sure that our simulation is correctly modelling whatever physical system we want to study. It is relatively easy to set up simulations, but how can we be sure that the "results" we get make sense? One of the best ways is to calculate from the simulation things that we can measure in experiment, and see if they agree. For example, we might want to simulate our system at a particular temperature and pressure, and measure the resulting density. If we repeat this over a range of temperatures at the same pressure, we will be able to plot an ''equation of state'', which we could compare to experimental measurements.

LAMMPS is able to calculate a great deal of thermodynamic information for us (you can see a full list of the properties it is able to calculate [http://lammps.sandia.gov/doc/thermo_style.html here]), but in these first simulations we are only interested in those properties specified in these commands:

<pre>
thermo_style custom time etotal temp press
thermo 10
</pre>

The first controls which properties will be printed out in the log file. In this case, we print how much time we have simulated so far (which is ''not'' the same as how long it has taken us to simulate it!), the total energy of the atoms, their temperature, and their pressure. The second line tells LAMMPS to print this information on every 10th timestep.

===Running the simulation===

'''Look at the lines below.'''
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(100/${timestep})
timestep ${timestep}

### RUN SIMULATION ###
run ${n_steps}
</pre>
'''The second line (starting "variable timestep...") tells LAMMPS that if it encounters the text ${timestep} on a subsequent line, it should replace it by the value given. In this case, the value ${timestep} is always replaced by 0.001.'''
 

'''<big> It is now time to run your first simulation, submit the input script with the data file in the intro folder of the files you have downloaded Try changing the timestep - what happens when you make the timestep larger?. </big>'''

===Visualising the trajectory===

The '''trajectory files''' contain the positions of all the atoms in the simulation, recorded at a set interval (for all of these simulations, this was every ten timesteps — this is controlled by the '''dump''' command in the input scripts). We use a programme called [http://www.ks.uiuc.edu/Research/vmd/ '''VMD'''] to view these trajectories, which you should find is already installed on both the desktop and laptop computers. You can run VMD from the start menu with '''Start''' -> '''All Programs''' -> '''University of Illinois''' -> '''VMD'''.

====Loading a Trajectory====

We'll start by looking at the output of the 0.02 timestep simulation. In the '''VMD Main''' window, select the menu option '''File''' -> '''New Molecule'''. Click the '''Browse''' button, then select the relevant trajectory file. In the '''Determine file type''' dropdown, select '''LAMMPS Trajectory'''. Then click '''Load'''.

You will see that the '''VMD 1.9.1 OpenGL Display''' window now shows a horrible mess. VMD's default behaviour is to draw lines between atoms which it thinks might be chemically bonded. Our system doesn't model chemical bonds, so we want to turn this off. In the '''VMD Main''' window, select the menu option '''Graphics''' -> '''Representations'''. This shows a list of "representations" of our atoms. You will see that at the moment, there is a single representation listed, and it is selected. It will have the ''Lines'' style, the ''Name'' colour, and the selection ''all''. "Selection" simply tells VMD which atoms we want it to draw. We want to show every atom, so the current selection is fine. The ''name'' colouring method just makes VMD give atoms colours according to their specified type. The colour isn't important to us, so we can leave this be too. The "style" tells VMD what we want it to display for each atom. Change the '''Drawing Method''' from ''Lines'' to ''VDW''. You will see that the mess of lines is replaced by a mess of low resolution, overlapping spheres. Change the '''Sphere Scale''' to 0.3, and the '''Sphere Resolution''' to 17. The result should look a little smoother. Close the '''Graphical Representations''' window. You will notice that in the bottom right of the '''VMD Main''' window, there is a small play button. Click this, and you will see the animated version of your simulation trajectory.

By clicking and dragging with the mouse, you can rotate the simulation box (though this may be sluggish). At any time, you can reset the view by pressing the equals key.

====Tracking a Single Particle====
To illustrate the periodic boundary conditions that we are using, we are going to draw almost all of the atoms as points, but we will pick a single atom at random to draw as a sphere. This will make it easy to see how a single atom moves through the box. Reset the display using the equals key, then use the '''VMD Main''' window controls to pause the trajectory and reset it to the first trajectory (play with the different buttons until you find the one that does this). You should see the perfect cubic lattice. Use the option '''Display''' -> '''Orthographic''' to change the drawing mode, then rotate the displayed crystal so that you are looking at one vertex (looking down the 111 direction, in crystallographic terms).

Open the '''Graphical Representations''' window again. Change the representation style from '''VDW''' to points, then click the '''Create Rep''' button. This creates a second representation, allowing a subset of the atoms to be drawn in a different way. The '''Selected Atoms''' box allows us to choose which atoms this representation applies to. We just want to pick two of them at random — VMD assigns every atom an index, from 0 to N-1. In our case, there are 125 atoms, so choose two numbers between 0 and 124. Changed the '''Selected Atoms''' field to <pre>index i or index j</pre> where i and j are your chosen numbers, press return, then change the '''Drawing Method''' to '''VDW'''. You should now see only two atoms represented by spheres, with the rest shown as small points. In the '''VMD Main''' window, click play. Try rotating the box, and changing the playback speed.

You will see that sometimes one of the spheres seems to change position across the box very rapidly — this occurs when it reaches one periodic boundary, and is reflected back across the other face. Try playing with some of the other representation types in VMD — it is a very powerful package, which is often used to render images of simulated proteins, so many of its options aren't relevant to our simple system!

===Checking equilibration===

When we first set up a simulation, it is very important to make sure that our system reaches an equilibrium state. We characterise equilibrium by the average values of thermodynamic quantities becoming constant (due to the approximations that we have made, there will always be fluctuations, but the average values will become constant).

In this section, we are going to plot the thermodynamic output of the simulation to see how long it takes to reach the equilibrium state (and indeed, whether this happens at all). Instructions are given below to import data from the LAMMPS log file into Microsoft Excel. Once you have the data in a spreadsheet, you can plot it. If you know how to use some of the other plotting software available on the chemistry computers (like Origin), you are welcome to use it.

# Open a blank Excel workbook
# Copy the data in the textfile into the first cell
# With these data highlighted, click the Data tab and "Text to Columns"
# Click "Delimited", continue and let it be space delimited
# Click finish

'''<big>TASK 7: </big>'''

=== What does it mean for a simulation to "reach equilibrium"? Why is this important in terms of sampling from an ensemble using molecular dynamics? [2] ===

=== Plot the energy (potential, kinetic and total), temperature and pressure, against time for the 0,001 timestep experiment [2] ===

=== Does the simulation reach equilibrium? How can you tell? [2] ===

=== Make a single plot which shows the energy vs. time for the timesteps you have simulated [2]. ===

=== Of the timesteps that you used, which timestep will you use for subsequent simulations and why? [6] ===
''(Think about what is happening "physically" as you increase/decrease the timestep. Also, what features of each timeseries are indicative of the simulation's "health"?)''

'''<big>This is the third section of the third year simulation experiment. You can return to the previous page, [[Third_year_simulation_experiment/Introduction_to_molecular_dynamics_simulation|Introduction to molecular dynamics simulation]], or jump ahead to the next section, [[Third year simulation experiment/Running simulations under specific conditions|Running simulations under specific conditions]].</big>'''

Third year simulation experiment/Files to download

2021-10-20T13:03:55Z

Ac5116: /* Getting the files for the experiment */

'''<big>This is the first section of the third year simulation experiment. You can return to the introduction page, [[Third year simulation experiment]], or jump ahead to the next section, [[Third year simulation experiment/Introduction to molecular dynamics simulation|Introduction to molecular dynamics simulation]].</big>'''

All of the simulations that you run in this experiment are going to be performed in Virtual Machines (VM). Before the beginning of your lab session, you should have received an invitation email to register/connect to the VM.

'''The small-scale simulations that we will perform in this experiment should not be too long a few minutes. '''

In each section of the exercise, we have tried to provide a number of mathematical and/or research exercises that you should attempt while you are waiting for the simulations in that section to be completed. You can also use this time to write your report on the previous sections!

In this first section, we will teach you how to login to the Virtual Machines and submit an example simulation. While you wait for that example to complete, you can move on to the next section and read about the theory of molecular dynamics simulations

==Connecting to the Virtual Machine==

You will be using a Microsoft Azure Lab Virtual Machine (VM) to access the software needed for the lab.

: 1. You will receive an invitation email, before your session starts, with the subject Register for Lab - IC_Chemistry_UK.
: 2. Follow the link to 'Register for the lab' in the email which will direct you to the Azure Lab Services page.
: 3. On the web page, you should see a box corresponding to the VM shared with you. In the bottom left, slide the toggle from left to right to start the VM (It may take some time to start running).
: 4. Once running, at the bottom right of the box, click the computer-style icon next to the three vertical dots.
A file with the extension rdp will be downloaded: IC_Chemistry_UK.rdp. Follow the instructions for your operating service below to use the file:

'''Windows'''

:: a. Navigate to where the file has downloaded and double click on the file to open.

'''Linux'''

:: a. Open a terminal window and through the terminal, go to the location where the file has downloaded.
:: b. Type the command: <pre> remmina IC_Chemistry_UK.rdp </pre> to run the file.

'''Mac'''

:: a. Download and install the Microsoft Remote Desktop app for Mac OS.
:: b. Open the Microsoft Remote Desktop app
:: c. On the top toolbar bar, click on the cog icon and then select Import from RDP file..
:: d. Navigate and select the downloaded rdp file.
:: e. There should now be an 'IC_Chemistry_UK' PC showing, double click on this to initialise.

: 5. You should be asked to Accept Certificate?, select Yes.
: 6. You will be asked to Enter authentication credentials:
:: a. Change the username into "chemistry" by removing "~/".
:: b. Enter Imperial2021 as the password.
The VM should now launch and you will be taken to a Windows desktop where you should be able to see the LAMMPS icon on the desktop. Double click this to launch.

If you have any problems with accessing the VM then let a demonstrator know.



==Getting the files for the experiment==

You can download all of the files that you will need for this experiment [https://imperialcollegelondon.box.com/s/mjxn9zxl67y10pgn8wkwlin05hmyqdsn from this address]. You should copy the folder '''ImperialChem-Year3SimExpt1415-master''' to a location of your choice. It contains a number of subfolders — one for each section of the experiment. Every subsequent page of this lab manual will begin with a line telling you which folder contains the necessary files, like this one:

<big>'''THE FILES THAT YOU NEED FOR THIS SECTION ARE FOUND IN THE "Intro" SUBFOLDER'''.</big>

Have a look in the '''"Intro"''' folder now. It contains a file called '''melt_crystal.in''', which you should open with a text editor (like Notepad++ or VSCode which you can find on Software Hub. VSCode is also on the Virtual Machines). This file is called an "input script", and it controls how the simulation software operates. We will perform all of our simulations with a software package called [http://lammps.sandia.gov LAMMPS]. Over the course of the experiment, you will learn what all of the commands in this file mean. The [https://lammps.sandia.gov/doc/Commands.html LAMMPS manual] contains a lot of valuable information about each of the commands, if you don't understand one of them or want to look up what some parameter means you can look it up there. To make life easier, we put "further info" links in the wiki. For now, we are going to use this file to run a few trial simulations.

'''<big>This is the first section of the third year simulation experiment. You can return to the introduction page, [[Third year simulation experiment]], or jump ahead to the next section, [[Third year simulation experiment/Introduction to molecular dynamics simulation|Introduction to molecular dynamics simulation]].</big>'''

Third year simulation experiment/Files to download

2021-10-20T12:51:45Z

Ac5116: /* Getting the files for the experiment */

'''<big>This is the first section of the third year simulation experiment. You can return to the introduction page, [[Third year simulation experiment]], or jump ahead to the next section, [[Third year simulation experiment/Introduction to molecular dynamics simulation|Introduction to molecular dynamics simulation]].</big>'''

All of the simulations that you run in this experiment are going to be performed in Virtual Machines (VM). Before the beginning of your lab session, you should have received an invitation email to register/connect to the VM.

'''The small-scale simulations that we will perform in this experiment should not be too long a few minutes. '''

In each section of the exercise, we have tried to provide a number of mathematical and/or research exercises that you should attempt while you are waiting for the simulations in that section to be completed. You can also use this time to write your report on the previous sections!

In this first section, we will teach you how to login to the Virtual Machines and submit an example simulation. While you wait for that example to complete, you can move on to the next section and read about the theory of molecular dynamics simulations

==Connecting to the Virtual Machine==

You will be using a Microsoft Azure Lab Virtual Machine (VM) to access the software needed for the lab.

: 1. You will receive an invitation email, before your session starts, with the subject Register for Lab - IC_Chemistry_UK.
: 2. Follow the link to 'Register for the lab' in the email which will direct you to the Azure Lab Services page.
: 3. On the web page, you should see a box corresponding to the VM shared with you. In the bottom left, slide the toggle from left to right to start the VM (It may take some time to start running).
: 4. Once running, at the bottom right of the box, click the computer-style icon next to the three vertical dots.
A file with the extension rdp will be downloaded: IC_Chemistry_UK.rdp. Follow the instructions for your operating service below to use the file:

'''Windows'''

:: a. Navigate to where the file has downloaded and double click on the file to open.

'''Linux'''

:: a. Open a terminal window and through the terminal, go to the location where the file has downloaded.
:: b. Type the command: <pre> remmina IC_Chemistry_UK.rdp </pre> to run the file.

'''Mac'''

:: a. Download and install the Microsoft Remote Desktop app for Mac OS.
:: b. Open the Microsoft Remote Desktop app
:: c. On the top toolbar bar, click on the cog icon and then select Import from RDP file..
:: d. Navigate and select the downloaded rdp file.
:: e. There should now be an 'IC_Chemistry_UK' PC showing, double click on this to initialise.

: 5. You should be asked to Accept Certificate?, select Yes.
: 6. You will be asked to Enter authentication credentials:
:: a. Change the username into "chemistry" by removing "~/".
:: b. Enter Imperial2021 as the password.
The VM should now launch and you will be taken to a Windows desktop where you should be able to see the LAMMPS icon on the desktop. Double click this to launch.

If you have any problems with accessing the VM then let a demonstrator know.



==Getting the files for the experiment==

You can download all of the files that you will need for this experiment [https://imperialcollegelondon.box.com/s/mjxn9zxl67y10pgn8wkwlin05hmyqdsn from this address]. You should copy the folder '''ImperialChem-Year3SimExpt1415-master''' to a location of your choice. It contains a number of subfolders — one for each section of the experiment. Every subsequent page of this lab manual will begin with a line telling you which folder contains the necessary files, like this one:

<big>'''THE FILES THAT YOU NEED FOR THIS SECTION ARE FOUND IN THE "Intro" SUBFOLDER'''.</big>

Have a look in the '''"Intro"''' folder now. It contains a file called '''melt_crystal.in''', which you should open with a text editor (like Notepad++ which you can find on your Desktop). This file is called an "input script", and it controls how the simulation software operates. We will perform all of our simulations with a software package called [http://lammps.sandia.gov LAMMPS]. Over the course of the experiment, you will learn what all of the commands in this file mean. The [https://lammps.sandia.gov/doc/Commands.html LAMMPS manual] contains a lot of valuable information about each of the commands, if you don't understand one of them or want to look up what some parameter means you can look it up there. To make life easier, we put "further info" links in the wiki. For now, we are going to use this file to run a few trial simulations.

'''<big>This is the first section of the third year simulation experiment. You can return to the introduction page, [[Third year simulation experiment]], or jump ahead to the next section, [[Third year simulation experiment/Introduction to molecular dynamics simulation|Introduction to molecular dynamics simulation]].</big>'''

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

2021-02-16T14:54:23Z

Ac5116: /* Importance Sampling */

'''<big>This is the third section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]], or jump ahead to the next section, [[Third year CMP compulsory experiment/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 <math>1\times 10^9</math> configurations per second with our computer. How long will it take to evaluate a single value of <math>\left\langle M\right\rangle_T</math>?'''

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 <math>\alpha</math> 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 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 it first in its list of the [http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=814652 "Top 10 Algorithms of the 20th Century"].

The algorithm is as follows:

# Start from a given configuration of spins, <math>\alpha_0</math>, with energy <math>E_0</math>.
# Choose a single spin '''at random''', and "flip" it, to generate a new configuration <math>\alpha_1</math>
# Calculate the energy of this new configuration, <math>E_1</math>
# Calculate the energy difference between the states, <math>\Delta E = E_1 - E_0</math>
## If the <math>\Delta E \leq 0</math> (the spin flipping decreased the energy), then we '''accept''' the new configuration.
##* We set <math>\alpha_0 = \alpha_1</math>, and <math>E_0 = E_1</math>, and then '''go to step 5'''
## If <math>\Delta E > 0</math>, the spin flipping increased the energy. By considering the probability of observing the starting and final states, <math>\alpha_0</math> and <math>\alpha_1</math>, 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, <math>R</math>, in the interval <math>[0,1)</math>
### If <math>R \leq \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''accept''' the new configuration.
###* We set <math>\alpha_0 = \alpha_1</math>, and <math>E_0 = E_1</math>, and then '''go to step 5'''
### If <math>R > \exp \left\{-\frac{\Delta E}{k_BT}\right\}</math>, we '''reject''' the new configuration.
###* <math>\alpha_0</math> and <math>E_0</math> are left unchanged. '''Go to step 5'''
# Update the running averages of the energy and magnetisation.
# Monte Carlo "cycle" 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, <math>\Delta E</math> 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!

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 "Understanding Molecular Simulation", by Frenkel and Smit, but a discussion of this is not required for the experiment. You can find this as an ebook in the library [https://library-search.imperial.ac.uk/permalink/f/o1297h/TN_cdi_askewsholts_vlebooks_9780080519982 here]. Chapter 7 of 'Statistical Mechanics: Theory and Molecular Simulation' by Tuckerman also gives a good description.

==Modifying IsingLattice.py==

Our IsingLattice object contains a stub function, '''montecarlocycle(T)''', which takes a single argument — the temperature at which we want to do the cycle. There are also attributes to measure the average energy and magnetisation: '''self.E''', '''self.E2''', '''self.M''', '''self.M2''', and '''self.n_cycles'''. You should use these values to keep a running sum of all of the energies, squared energies, magnetisations, and squared magnetisations that your system has sampled, as well as the number of cycles that have been performed. Finally, there is also a function '''statistics()''', that takes no arguments.

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

You have been provided with a script, '''ILanim.py''', which will use your IsingLattice object to run a Monte Carlo simulation and display the output on the screen in real time. By default, this simulation will run an <math>8\times 8</math> lattice at <math>T=1.0</math>, 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 '''statistics()''' function to print the average energy and magnetisation from the simulation.

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

'''<big>This is the third section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Calculating the energy and magnetisation|Calculating the energy and magnetisation]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Accelerating the code|Accelerating the code]].</big>'''

Programming a 2D Ising Model/Locating the Curie temperature

2020-11-26T14:15:36Z

Ac5116: /* Finding the peak in C */

'''<big>This is the eighth (and final) section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Determining the heat capacity|Determining the heat capacity]], or go back to the [[Third year CMP compulsory experiment|Introduction]].</big>'''

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 <math>T = T_C</math>. 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 <math>T_{C, L}</math> is the Curie temperature of an <math>L\times L</math>lattice, <math>T_{C,\infty}</math> is the Curie temperature of an infinite lattice, and <math>A</math> is a constant in which we are not especially interested.

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

==Polynomial fitting==

To find the temperature at which the heat capacity has its maximum, 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 polyfit] and [http://docs.scipy.org/doc/numpy/reference/generated/numpy.polyval.html#numpy.polyval polyval] functions for this purpose. 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 vs T, as well as a fitted polynomial. Try changing the degree of the polynomial to improve the fit — in general, it might be difficult to get a good fit! Attach a PNG of an example fit to your report.'''

===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 <math>T_{\mathrm{min}}</math> and <math>T_{\mathrm{max}}</math>, 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 of the heat capacity! 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: peak_T_range, containing 1000 equally spaced temperature values between <math>T_{\mathrm{min}}</math> and <math>T_{\mathrm{max}}</math> (whatever you chose those values to be), and fitted_C_values, containing the fitted heat capacity at each of those points. Use the NumPy max function to find the maximum in C. If you store the maximum value of C in the variable Cmax, 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 colums: the lattice side length (2,4,8, etc.), and the temperature at which C is a maximum. This is your estimate of <math>T_C</math> for that side length. Make a plot that uses the scaling relation given above to determine <math>T_{C,\infty}</math>. By doing a little research online, you should be able to find the theoretical exact Curie temperature for the infinite 2D Ising lattice. How does your value compare to this? Are you surprised by how good/bad the agreement is? Attach a PNG of this final graph to your report, and discuss briefly what you think the major sources of error are in your estimate.'''

'''<big>This is the eighth (and final) section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Determining the heat capacity|Determining the heat capacity]], or go back to the [[Third year CMP compulsory experiment|Introduction]].</big>'''

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

2020-11-23T11:55:29Z

Ac5116: /* Getting the files for the experiment */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 Box [https://imperialcollegelondon.box.com/v/ising-lab-files here].''' If this link is broken you can alternatively try an older(slightly broken) version [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use the IPython Qt or the regular ipython console (which you can access via spyder) to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''IPython QTConsole'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment.

NOTE: QTConsole does not work for running the ILanim.py script you run latter in the lab. The best solution is to just run the script from the command line using '''python ILanim.py''' (if you're in the directory where files are stored). On a Mac/Linux computer you run this from '''Terminal''' or on windows you can use the app 'Conda Prompt'

When using the QTConsole, you should start the session by 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 you may need to restart the console)!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<code>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_cycles = 0
</code>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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:

<code>
import IsingLattice

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

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

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

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

<code>

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

</code>

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

===The Constructor===

The only function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice is then created, which is a NumPy array with the requested number of rows and columns. The function <pre>np.random.choice</pre> is used to create a random collection of spins for our initial configuration. You can see the documentation for this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 IPython Qt console using the command '''<pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

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

2020-11-23T11:13:51Z

Ac5116: /* The model */

'''<big>This is the first section of the third year CMP experiment. You can return to the introduction page, [[Third year CMP compulsory experiment]], 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, <math>s_i</math>, each of which can either have the value "up" (<math>+1</math>) or "down" (<math>-1</math>), arranged on a lattice. In one dimension, this takes the form of a "chain" of adjacent spins, <math>s_1</math> to <math>s_N</math>. 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=5x5), and three (N=5x5x5) 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 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 <math>J</math> is a constant which controls the strength of interaction and the notation <math>j\ \in\ \mathrm{neighbours}\left(i\right)</math> indicates that spin <math>j</math> must lie in a lattice cell adjacent to spin <math>i</math>. 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 <math>J>0</math>) is <math>E\ =\ -DNJ</math>, where <math>D</math> is the number of dimensions and <math>N</math> is the total number of spins. What is the multiplicity of this state? Calculate its entropy.'''

==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 (<math>D=3,\ N=1000</math>)? 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'', <math>T_C</math>. Below <math>T_C</math>, we observe the highly ordered magnetised phase, while above <math>T_C</math> 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 <math>D = 3,\ N=1000</math> at absolute zero?'''

We are interested, therefore, in determining the magnetisation <math>M</math>, and energy <math>E</math>, of the system as functions of temperature <math>T</math>. 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 <math>\alpha</math> is used as a shorthand to represent all spins in the system, and <math>Z</math> is the partition function. '''If you are doing this experiment before the statistical thermodynamics lecture course, you may not have seen expressions of this nature before. If so, please ask or e-mail the demonstrator, who will be happy to give you a brief introduction to the topic.''' 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 the 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 <math>M</math> and <math>E</math> can be solved, but first we must set up a Python script to model the lattice.

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

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

2020-10-26T09:35:34Z

Ac5116: /* Getting the files for the experiment */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 Box [https://imperialcollegelondon.box.com/v/ising-lab-files here].''' If this link is broken you can alternatively try an older(slightly broken) version [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use the IPython Qt console to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''IPython QTConsole'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment. When using the QTConsole, you should start the session by 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!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<code>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_cycles = 0
</code>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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:

<code>
import IsingLattice

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

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

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

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

<code>

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

</code>

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

===The Constructor===

The only function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice is then created, which is a NumPy array with the requested number of rows and columns. The function <pre>np.random.choice</pre> is used to create a random collection of spins for our initial configuration. You can see the documentation for this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 IPython Qt console using the command '''<pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

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

2020-10-20T14:16:17Z

Ac5116: /* Getting the files for the experiment */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 Box [https://imperialcollegelondon.box.com/v/ising-lab-files here]. If this link is broken you can alternatively try an older(slightly broken) version [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use the IPython Qt console to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''IPython QTConsole'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment. When using the QTConsole, you should start the session by 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!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<code>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_cycles = 0
</code>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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:

<code>
import IsingLattice

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

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

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

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

<code>

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

</code>

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

===The Constructor===

The only function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice is then created, which is a NumPy array with the requested number of rows and columns. The function <pre>np.random.choice</pre> is used to create a random collection of spins for our initial configuration. You can see the documentation for this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 IPython Qt console using the command '''<pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

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

2020-10-20T14:15:49Z

Ac5116: /* Getting the files for the experiment */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 Box [https://imperialcollegelondon.box.com/v/ising-lab-files here]. If this link is broken you can alternatively try an older(slightly broken) [here https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use the IPython Qt console to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''IPython QTConsole'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment. When using the QTConsole, you should start the session by 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!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<code>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_cycles = 0
</code>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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:

<code>
import IsingLattice

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

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

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

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

<code>

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

</code>

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

===The Constructor===

The only function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice is then created, which is a NumPy array with the requested number of rows and columns. The function <pre>np.random.choice</pre> is used to create a random collection of spins for our initial configuration. You can see the documentation for this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 IPython Qt console using the command '''<pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

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

2020-10-20T13:03:35Z

Ac5116: /* A note about Python Classes */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 file [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use the IPython Qt console to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''IPython QTConsole'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment. When using the QTConsole, you should start the session by 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!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<code>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_cycles = 0
</code>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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:

<code>
import IsingLattice

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

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

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

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

<code>

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

</code>

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

===The Constructor===

The only function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice is then created, which is a NumPy array with the requested number of rows and columns. The function <pre>np.random.choice</pre> is used to create a random collection of spins for our initial configuration. You can see the documentation for this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 IPython Qt console using the command '''<pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

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

2020-10-20T13:03:11Z

Ac5116: /* A note about Python Classes */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 file [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use the IPython Qt console to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''IPython QTConsole'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment. When using the QTConsole, you should start the session by 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!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<code>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_cycles = 0
</code>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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:

<code>
import IsingLattice

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

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

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

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

<code>
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
</code>

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

===The Constructor===

The only function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice is then created, which is a NumPy array with the requested number of rows and columns. The function <pre>np.random.choice</pre> is used to create a random collection of spins for our initial configuration. You can see the documentation for this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 IPython Qt console using the command '''<pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

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

2020-10-20T13:02:26Z

Ac5116: /* The Constructor */

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''

==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 file [https://github.com/niallj/ImperialChem-Year3CMPExpt/archive/master.zip here].

We recommend that you use the IPython Qt console to do your simulation work. This should be available on all of the college computers, where you can access it through '''Start''' -> '''All Programs''' -> '''Anaconda (64-bit)''' -> '''IPython QTConsole'''. On your own Mac or Linux computer you can simply type 'ipython' into the terminal. You are very welcome to use the IPython Notebook, that you were introduced to last year, to perform some of the analysis tasks later in the experiment. When using the QTConsole, you should start the session by 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!

For this section of the experiment, you will need the files '''IsingLattice.py''' and '''ILcheck.py'''. Extract them from the zip file to a convenient location in your H drive, then in the IPython QTConsole use cd H:\... to move to the correct directory.

==Modifying the files==

The file '''IsingLattice.py''' contains a Python class called IsingLattice. Open this file in your favourite text editor. 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 functions 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 that you may not have met before, by defining a ''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 IsingLattice:

<code>
class IsingLattice:
</code>

Each class is allowed to have certain attributes and functions (known as methods) associated with it. Lines 5-10 of IsingLattice.py create 5 variables that every object of type IsingLattice has:

<code>
E = 0.0

E2 = 0.0

M = 0.0

M2 = 0.0

n_cycles = 0
</code>

There are also functions defined, on lines 12, 17, 22, etc. We can see, for example, that each IsingLattice object has the functions energy(), magnetisation(), montecarlostep(T), and statistics() (the function __init__ will be explained shortly). Each of these functions 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:

<code>
import IsingLattice

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

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

When writing code "within" the class definition, all variables and functions belonging to the class must be accessed using the "self." notation. You can see an example of this in the montecarlostep() function:

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

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

<code>
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
</code>

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

===The Constructor===

The only function currently filled in is __init__. This is a special function called a ''constructor''. Whenever a new IsingLattice object is created, using code like <code>il = IsingLattice(5,5)</code>, this function 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 self.lattice is then created, which is a NumPy array with the requested number of rows and columns. The function <pre>np.random.choice</pre> is used to create a random collection of spins for our initial configuration. You can see the documentation for this function [https://numpy.org/doc/stable/reference/random/generated/numpy.random.choice.html here].

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

==Testing the files==

When you have completed the energy() and magnetisation() functions, we need to test that they work correctly. You are welcome to do this yourself — from the IPython Qt console (ensuring that you are in the correct directory), you can simply type <pre>import IsingLattice</pre> then create objects of your class with the command <pre>il = IsingLattice.IsingLattice(n_rows, n_columns)</pre>
We have also 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 IPython Qt console using the command '''<pre>%run ILcheck.py</pre> The displayed window has a series of control buttons in the bottom left, one of which will allow you to export the figure as a PNG image. Save an image of the ILcheck.py output, and include it in your report.'''

'''<big>This is the second section of the third year CMP experiment. You can return to the previous page, [[Third_year_CMP_compulsory_experiment/Introduction to the Ising model|Introduction to the Ising model]], or jump ahead to the next section, [[Third year CMP compulsory experiment/Introduction to Monte Carlo simulation|Introduction to Monte Carlo simulation]].</big>'''