'''<big>This is the optional experiment which may be chosen by any third year student. If you are looking for the compulsory simulation experiment for students studying chemistry with molecular physics, you will find it [[Third_year_CMP_compulsory_experiment|here]].</big>'''

==Introduction==

Computer simulation is widely used to study a huge variety of chemical phenomena, from the behaviour of exotic materials under extreme conditions to protein folding and the properties of biological systems such as lipid membranes. In this experiment, we hope to give you a gentle introduction to one of the most powerful methods for the simulation of chemical systems, '''molecular dynamics simulation'''.

This course closely follows some of the ideas introduced in Professor Bresme's statistical thermodynamics lecture course. Do not worry if you are doing this experiment before the lectures begin, everything that you need to know to complete the experiment will be explained in these instructions. We will begin with a brief overview of the fundamental theory behind the method, before you start running your own simulations of a simple liquid using the college's [http://www.imperial.ac.uk/admin-services/ict/self-service/research-support/hpc/ high performance computing] facilities.

At the end of this experiment, you will have performed your own simulations using state of the art software packages used by researchers all around the world, and used those simulations to calculate both structural and dynamic properties of a simple liquid. You will have seen how the concepts of statistical physics introduced by Professor Bresme are needed to calculate thermodynamic quantities such as temperature and pressure in computer simulations, and you will see how computers can be used to validate those concepts.

All of the information that you need to complete the experiment is provided in these wiki pages. We have also tried to provide links to external resources and relevant textbooks where possible — unless explicitly stated, reading these resources '''is not required'''; they are provided only as further information for those interested in the subject.

==Assessment==

At the end of this experiment you must submit a report (pdf format via turnitin). The report should be structured:

* Introduction Questions (20% of the total mark)
* Results and Discussion (60% of total mark)
* Conclusion Questions (20% of total mark)
* References

Relevant supplementary material can be added at the end of the report so long as it supports your discussion. Please limit the introduction and conclusion sections (which are now just answers to questions) to a maximum of ~1.5 pages at size 12 font Times New Roman, or equivalent (e.g. if you are using latex or simply prefer another font). It is certainly possible to answer the introduction and conclusion questions with under this length limit, and as with most scientific writing, shorter concise but well thought out answers are preferable. We have not imposed a hard word limit because, equally, you should spend the majority of this lab learning, not editing and rewriting answers to just satisfy the word limit.

Please note that five "floating" marks have been reserved in the results and discussion section. These can be used to reward particularly insightful comments and/or explanations, on the basis of good/bad presentation, use of scientific language, or for other reasons.

The final task of the lab (task 10 - MSD diffusion simulations) is NOT part of the lab this year. The section has been left on the wiki for any interested student to give it a go. However, it should NOT be submitted as part of your lab report, and will not be marked.

==Getting Help==

Please feel free to ask any demonstrators for help during the lab sessions - that is what we are here for. Questions can also be asked via email.

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

==Structure of this Experiment==

This experimental manual has been broken up into a number of subsections. To help you plan your time it is suggested you complete the following at these times:

Monday (morning session): Theory - Introduction to molecular dynamics simulations

Monday (afternoon session): Theory + Equilibration (submit your files for running)

Tuesday (morning session): Equilibration (analyse your files)

Tuesday (afternoon session): Running simulations under specific conditions (submit and read)

Thursday (morning session): Checkpoint for progress. Submit your input files for the radial distribution function and analyse your equation of state from the previous section

Thursday (afternoon session).

Friday (morning session): Analyse MSD diffusion simulations

Friday (afternoon session): Analyse MSD diffusion simulations. Report write-up.

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 simulation experiment/Files to download|Downloading the Files]]
# [[Third year simulation experiment/Introduction to molecular dynamics simulation|Introduction to molecular dynamics simulation]]
# [[Third_year_simulation_experiment/Equilibration|Equilibration]]
# [[Third_year_simulation_experiment/Running_simulations_under_specific_conditions|Running simulations under specific conditions]]
# [[Third year simulation experiment/Structural properties and the radial distribution function|Structural properties and the radial distribution function]]
# [[Third year simulation experiment/Dynamical properties and the diffusion coefficient|Dynamical properties and the diffusion coefficient]]

== Introduction Questions ==
# Give a brief account on what molecular dynamics is, and (broadly speaking) what it aims to achieve. [5]
# Give two examples of instances where molecular dynamics simulations have/can be used (in a helpful way) instead of or in synergy with experimental techniques. In each case explain why. [10] (''Please consider only chemically relevant examples. E.g. To study geochemical processes which occur under extreme temperature/pressure conditions, and therefore cannot be studied experimentally.'')
# Explain what is meant by a thermodynamic ensemble, and the term "conserved quantity". [2]
# What is the Ergodic hypothesis? Explain its relevance to molecular dynamics simulations. [3]

2018-04-16T13:25:11Z

Org12: /* Conclusion */

=Tasks=
==Section 2: Introduction==

'''<big>TASK</big>: Open the file HO.xls. In it, the velocity-Verlet algorithm is used to model the behaviour of a classical harmonic oscillator. Complete the three columns "ANALYTICAL", "ERROR", and "ENERGY": "ANALYTICAL" should contain the value of the classical solution for the position at time <math>t</math>, "ERROR" should contain the ''absolute'' difference between "ANALYTICAL" and the velocity-Verlet solution (i.e. ERROR should always be positive -- make sure you leave the half step rows blank!), and "ENERGY" should contain the total energy of the oscillator for the velocity-Verlet solution. Remember that the position of a classical harmonic oscillator is given by <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math> (the values of <math>A</math>, <math>\omega</math>, and <math>\phi</math> are worked out for you in the sheet).'''

Excel file attached [[:File:AD5215_HO.xls|here]].

'''<big>TASK</big>: For the default timestep value, 0.1, estimate the positions of the maxima in the ERROR column as a function of time. Make a plot showing these values as a function of time, and fit an appropriate function to the data.'''

[[File:Ad5215 error vs time.png]]

'''<big>TASK</big>: Experiment with different values of the timestep. What sort of a timestep do you need to use to ensure that the total energy does not change by more than 1% over the course of your "simulation"? Why do you think it is important to monitor the total energy of a physical system when modelling its behaviour numerically?'''

The change in energy goes down as the timestep value becomes smaller. For a timestep of 0.25 the change in energy is 1.58% while for a timestep of 0.05 the change in energy is 0.06% . The energy change is 1.01% for a timestep of 0.2 . A timestep that is too large could lead to the simulation effectively "missing" any changes in the system that happen on a shorter timescale than that of the timestep. Therefore, it is important to monitor the energy to ensure that the change is not too drastic and we are observing the behaviour of the system closely.

'''<big>TASK:</big> For a single Lennard-Jones interaction, <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>, find the separation, <math>r_0</math>, at which the potential energy is zero. What is the force at this separation? Find the equilibrium separation, <math>r_{eq}</math>, and work out the well depth (<math>\phi\left(r_{eq}\right)</math>). Evaluate the integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math> when <math>\sigma = \epsilon = 1.0</math>.

''Find the separation at which the potential energy is zero''

[[File:Ad5215 lj zero pot.JPG]]

''FInd the force at <math>r=\sigma</math>''

The force is the derivative of potential wrt to distance:

[[File:Ad5215 lj Force.JPG]]

At separation <math>r=\sigma</math> this will be

[[File:Ad5215 force(R).JPG]]

''Equilibrium separation and well depth''<\span>

The equilibrium separation is the separation when <math> \frac{d \phi}{dr} = 0.</math>

[[File:Ad5215 lj equilibrium req.JPG]]

The well depth at this separation is

[[File:Ad5215 ls lj epsilon.JPG]]

 negative epsilon 

''Integrals''

[[File:Ad5215 lj int1.JPG]]

[[File:Ad5215 int2.JPG]]

[[File:Ad5215 int3.JPG]]

'''<big>TASK</big>: Estimate the number of water molecules in 1ml of water under standard conditions. Estimate the volume of <math>10000</math> water molecules under standard conditions.'''

<math> V = 1\ mL, \ \rho = 1\ g/mL, \ M = 18\ g/mol</math>

<math>
m = V \rho = 1\ g
</math>

<math>
N = nN_a = \frac{m}{M} \times N_a = \frac{6.022 \times 10^{23}}{18} = 3.35 \times 10{22}
</math>

N molecules occupy 1 mL. Therefore, the volume of 1 molecule of water will be <math>V_0 = \frac{1}{N} = \frac{1}{3.35 \times 10{22}} = 2.99 \times10^{-23}</math>

The volume of 1000 molecules will be <math>1000 \times V_0 = 2.99 \times10^{-20} </math>

'''<big>TASK</big>: Consider an atom at position <math>\left(0.5, 0.5, 0.5\right)</math> in a cubic simulation box which runs from <math>\left(0, 0, 0\right)</math> to <math>\left(1, 1, 1\right)</math>. In a single timestep, it moves along the vector <math>\left(0.7, 0.6, 0.2\right)</math>. At what point does it end up, ''after the periodic boundary conditions have been applied''?'''

It ends up at position (0.2, 0.1, 0.7).

'''<big>TASK</big>: The Lennard-Jones parameters for argon are <math>\sigma = 0.34\mathrm{nm}, \epsilon\ /\ k_B= 120 \mathrm{K}</math>. If the LJ cutoff is <math>r^* = 3.2</math>, what is it in real units? What is the well depth in <math>\mathrm{kJ\ mol}^{-1}</math>? What is the reduced temperature <math>T^* = 1.5</math> in real units?

<math>r^* = \frac{r}{\sigma}</math>

<math>r = \sigma \times r^*</math>

<math>r = 0.34 \times 3.2 = 1.088\ nm </math>

<math>\frac{\epsilon}{k_B} = 120\ K </math>

<math>\epsilon = 120 \times k_B </math> for 1 particle

For a mole of particles: <math>\epsilon = 120 \times k_B \times N_A </math>

<math> \epsilon = 0.997\, kJ mol^{-1} </math>

<math>T^* = \frac{k_BT}{\epsilon}</math>

<math>T = T^* \times \frac{\epsilon}{k_BT}</math>

<math>T = 1.5 \times 120 </math>

<math>T = 180\ K</math>

==Section 3: Equilibration==
'''<big>TASK</big>: Why do you think giving atoms random starting coordinates causes problems in simulations? Hint: what happens if two atoms happen to be generated close together?'''

Randomly generated positions can lead to two atoms being very close together, which would result in a large repulsive potential. This would then affect any propagation in time of the system, which would lead to undesirable behaviour, especially when using larger timesteps. The system would most likely behave "appropriately" for small enough timesteps, but this would require running longer simulations. This would be less effective; a larger timestep that still results in an accurate simulation is ideal. good! 

'''<big>TASK</big>: Satisfy yourself that this lattice spacing corresponds to a number density of lattice points of <math>0.8</math>. Consider instead a face-centred cubic lattice with a lattice point number density of 1.2. What is the side length of the cubic unit cell?'''

Density: <math> \rho = \frac{N}{V} = \frac{N}{l^3} </math> ---> Length: <math> l=\sqrt[3]{\frac{N}{\rho}} </math>

For a simple cubic lattice with <math> \rho = 0.8 </math>, the length, l, is:

<math> \sqrt[3]{\frac{1}{0.8}}=1.07722 </math>

For a face-centered cubic lattice with <math> \rho = 1.2 </math>, the length, l, is:

<math> \sqrt[3]{\frac{4}{1.2}}=1.4938 </math>

'''<big>TASK</big>: Consider again the face-centred cubic lattice from the previous task. How many atoms would be created by the create_atoms command if you had defined that lattice instead?'''

The box would still contain <math> 10 \times 10 \times 10 = 1000 </math> lattice units. For an FCC there are 4 atoms per lattice unit. Therefore the total number of atoms would be <math> 4 \times 1000 = 4000 </math> atoms.

'''<big>TASK</big>: Using the [http://lammps.sandia.gov/doc/Section_commands.html#cmd_5 LAMMPS manual], find the purpose of the following commands in the input script:'''

<pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

"Mass" sets the mass of a particular type of atom (in this case, the mass of type 1 atoms is 1). "Pair" refers to pair potentials. The lj/cut command computes the 12/6 Lennard-Jones potential, cut sets the cut-off for r. Pair_coeff sets the values for the 2 parameters, sigma and epsilon.

'''<big>TASK</big>: Given that we are specifying <math>\mathbf{x}_i\left(0\right)</math> and <math>\mathbf{v}_i\left(0\right)</math>, which integration algorithm are we going to use?'''

 Velocity Verlet. A simple Verlet algorithm wouldn't require the initial velocity, but would instead require <math>\mathbf{x}_i\left(-\delta t\right)</math>.

'''<big>TASK</big>: 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}
run 100000
</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. In light of this, what do you think the purpose of these lines is? Why not just write:'''
<pre>
timestep 0.001
run 100000
</pre>

The line "variable timestep equal 0.001" defines a variable timestep which is the assigned a value. This allows for the variable to be called later on if needed. This is convenient for the user as it means that if the same variable is required multiple times (in this case, the variable timestep is called twice) changing its value is easier, as this only needs to be done once (in the line defining the variable).

==Section 4: Running simulations under specific conditions==

'''<big>TASK</big>: We need to choose <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>

'''Solve these to determine <math>\gamma</math>.'''

[[File:Ad5215 ljls gamma.JPG]]

'''<big>TASK</big>: Use the [http://lammps.sandia.gov/doc/fix_ave_time.html manual page] to find out the importance of the three numbers ''100 1000 100000''. How often will values of the temperature, etc., be sampled for the average? How many measurements contribute to the average? Looking to the following line, how much time will you simulate?'''

From the manual, the structure of the command is "ave/time Nevery Nrepeat Nfreq".

*Nevery (100) = use input values every this many timesteps (total no. of data points is 100,000 so this will give 1000 values to be averaged in the next step; records an average of 100 values)

*Nrepeat (1000) = no. of times to use input values for calculating averages (i.e. average over 1000 values)

*Nfreq (100,000) = calculate averages every this many timesteps (same no. specified in the "run" command)

The timestep is 0.0025 and the simulation runs for 100,000 steps. Therefore we are simulating a total time of 250 seconds .

==Section 7: Dynamical properties and the diffusion coefficient==

'''<big>TASK</big>: In the theoretical section at the beginning, the equation for the evolution of the position of a 1D harmonic oscillator as a function of time was given. Using this, evaluate the normalised velocity autocorrelation function for a 1D harmonic oscillator (it is analytic!):

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} v\left(t\right)v\left(t + \tau\right)\mathrm{d}t}{\int_{-\infty}^{\infty} v^2\left(t\right)\mathrm{d}t}</math>

For a HO model:
[[File:Ad5215 Vacf HO model.JPG]]

The VACF is given by

[[File:Ad5215 VACF deriv.JPG]]

We evaluate the two integrals separately. Integral 2 is

[[File:Ad5215 Vacf i2.JPG]]

Using [[File:Ad5215 Vacf cos(2x).JPG]] we can write

[[File:Ad5215 Vacf i2 2.JPG]]

Integral 1 is

[[File:Ad5215 Vacf i1 1.JPG]]

Using [[File:Ad5215 Vacf sin(A+B).JPG]] we can write

[[File:Ad5215 Vacf i1 2.JPG]]

We evaluate the last integral separately.

[[File:Ad5215 Vacf i3.JPG]]

Substituting this back we obtain:

[[File:Ad5215 Vacf i1 3.JPG]]

Substituting I1 and I2 into the formula for C we obtain

[[File:Ad5215 Vacf c1.JPG]]

Sin is an odd function, i.e. [[File:Ad5215 Vacf sinodd.JPG]]. Thus

[[File:Ad5215 Vacf c2.JPG]]

 There is definitely a shorter way of doing this, but yes. 

=Report=

==Abstract==

Solid, liquid and gaseous systems were modeled using LAMMPS and a 12-6 Lennard-Jones forcefield. A suitable timestep of 0.0025 was determined. Simulations were run to obtain thermodynamic data (temperatures, pressures, densities, heat capacities). Calculated densities were found to be lower than those predicted by the ideal gas law. The simulated heat capacities showed a trend, decreasing with increasing temperature. The RDF was calculated for systems in all three phases. As expected, the RDF for a solid showed peaks decaying in amplitude. Lattice spacings and coordination numbers for the solid FCC lattice were calculated. The MSD and the VACF were plotted for the same three systems and the diffusion coefficient was calculated for both measurements. The two methods did not result in identical values; still, the difference was only 0.67% for the diffusion coefficients in the gas phase.
 Good, concise abstract that summaries main results. Perhaps 1 sentence of motivation would have been nice. 

==Introduction==

Food constitutes a huge part of our lives, whether we actively think about it or not. Cooking, in some form or other, has been around ever since the first human realised that fire makes raw meat tastier and more convenient to eat. Nowadays, the food industry is enormous and cooking has become a successful blend of art and science. Perhaps the most pertinent example of this is molecular gastronomy, a subdiscipline of food science that seeks to investigate the physical and chemical transformations of ingredients during the cooking process. In their review, Balham et al refer to molecular gastronomy as an "emerging scientific discipline".<ref>P. Barham, L. H. Skibsted, W. L. P. Bredie and J. Risbo, ''Symp.
A Q. J. Mod. Foreign Lit.'', 2010, 2313–2365.</ref>

In fact, molecular gastronomy relies heavily on processes such as gelification and infusion and on materials such as gels, foams and powders. It also uses equipment that is heavily reminiscent of a laboratory - some kitchens are fitted with butane burners, syringes and dehydrators.

Clever presentations and unusual sensations and are piquing the interest of many people; in some places, molecular gastronomy restaurants have become tourist attractions.<ref>D. Tüzünkan and A. Albayrak, ''Procedia
- Soc. Behav. Sci.'', 2015, '''195''', 446–452.</ref> Simulating the thermodynamic properties of systems can provide a better understanding of physical systems and can lead to the growth and development of this industry.

 An interesting topic and motivation! An explicit connection between gastronomy and the results you intend to present/discuss would be nice. For example how is the diffusion coefficient relevant? The introduction of a scientific paper usually includes the background theory, such as in your case, the equations for diffusion coefficient. You have included this in the methodology, which while logical, is not standard practice for most papers/journals. 

==Aims and Objectives==

* become familiar with LAMMPS and how simulating a physical system works
* model the behaviour of systems under the 12-6 Lennard-Jones potential
* calculate the physical properties (temperature, pressure, density, etc) of a system using such simulations
* comparing the results of the simulations to theory (e.g. simulated density vs density given by the ideal gas law)
* calculating the diffusion coefficient for systems in different phases (gas, liquid, solid) by two different methods (from the MSD and from the VACF)

==Methods==
'''General methods'''

A 12-6 Lennard-Jones system was modeled usings LAMMPS and all simulations were run using Imperial's High Performance Computer.
The potential for such a system is given by<ref name=":0">P. Atkins, J. De Paula, ''Physical
Chemistry'', OUP Oxford, 9th edn., 2009.</ref>:

<math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>

All calculations used reduced units:<math>r^* = r/{\sigma}</math>, <math>E^* = E/{\epsilon}</math> and <math>T^* = {k_BT}/{\epsilon}.</math>
The Lennard-Jones parameters <math>\sigma</math> and <math>\epsilon</math> were set to 1.0 for all simulations. The cutoff for Lennard-Jones interactions was set at r*=3 unless otherwise stated. The mass of the atoms was set to 1.0. The temperature, pressure, lattice density and timestep values were varied. For all calculations the velocity Verlet algorithm was employed. The atoms were assigned random velocities within the simulation, while ensuring that the Boltzmann distribution of states is followed.

'''Determining a suitable timestep'''

A simple cubic lattice with a density of 0.8 was defined and the simulation was populated with 1000 atoms (10 x 10 x 10 dimensions). The ensemble was defined as the microcanonical (NVE) ensemble. Five values for the timestep were tested: 0.015, 0.01, 0.0075, 0.0025, 0.001 and each simulation was run for a total time of 100 seconds. Values for the energy, temperature and pressure of the system were recorded at each step.

'''Variation of density with temperature and pressure'''

A simple cubic lattice with a density of 0.8 was defined and the simulation was populated with 3375 atoms (15 x 15 x 15 dimensions). The timestep for all simulations was set to 0.0025. The ensemble was defined as NPT and 10 different thermodynamic states were simulated (two pressures, 2.5 and 3.5, each associated with five temperatures: 3.0, 6.0, 9.0, 12.0 and 15.0). Values for the energy, temperature and pressure of the system were recorded at each step, as well as average values for the density, temperature and pressure of the system at the end of the simulation. Plots of density vs time were obtained, both for the simulated data and for densities predicted by the ideal gas law<ref name=":0" />, <math> PM = \rho RT.</math>

'''Heat capacity calculations'''

A simple cubic lattice with a density of 0.2 was defined and populated with 3375 atoms. The timestep for all simulations was set to 0.0025. An NVT ensemble was simulated for five temperatures: 2.0, 2.2, 2.4, 2.6 and 2.8. Then, a subsequent simulation was run to establish an NVE ensemble and measure the properties of the system. Average values for temperature, energy, volume and heat capacity were calculated. This procedure was repeated for a simple cubic lattice with a density of 0.8. An example input script can be found [[:File: Example_script_heatcap_ad5215.in|here]].

The heat capacity of a system is given by<ref name=":0" />:

<math>C_V = \frac{\partial E}{\partial T} = \frac{\mathrm{Var}\left[E\right]}{k_B T^2} = N^2\frac{\left\langle E^2\right\rangle - \left\langle E\right\rangle^2}{k_B T^2}</math>

where E is the internal energy and

<math>{\mathrm{Var}\left[E\right]}={\left\langle E^2\right\rangle - \left\langle E\right\rangle^2}</math>

is the variance in internal energy. The <math>N^2</math> term is required because LAMMPS automatically outputs the energy '''per atom''', not the '''total''' energy.

[[File:Ad5215 LJfluid phase diag.JPG|thumb|Fig. 1: Phase diagram for the Lennard-Jones fluid]]
'''RDF calculations'''

Three systems (a liquid, a solid and a gas) were modeled and populated with 3375 atoms each. Temperature and density values were taken from the Lennard-Jones fluid phase diagram<ref>J. P. Hansen and L.
Verlet, ''Phys. Rev.'', 1969, '''184''', 151–161.</ref> reproduced in Fig. 1. These were defined as:

*solid: fcc lattice, temperature 1.2, density 1.2;
*liquid: sc lattice, temperature 1.2 , density 0.8;
*vapour: sc lattice, temperature 1.2, density 0.05.

The ensemble was defined as NVT. The timestep for all simulations was set to 0.002. The trajectories of the atoms were recorded and VMD was used to calculate the radial distribution function and its integral from these trajectories. The data was then analysed using Python.

'''Diffusion coefficient calculations'''

The same three systems as above were modeled, this time with 8000 atoms each. The timestep was set to 0.002 and each simulation was run for 5000 steps. The Lennard-Jones cutoff was set to 3.2. The ensemble was defined as NVT. The mean squared displacement (MSD) and the velocity autocorrelation function (VACF) at each step were calculated for all systems. The data was analysed using Python. The MSD plots were fitted to a straight line and the gradient was used to calculate the diffusion coefficient. The VACF integrals were plotted as a function of time and, again, used to calculate the diffusion coefficient. The same data analysis was conducted using supplied data which modeled larger systems.

The '''MSD''' is a measure of the deviation of the position of a particle with respect to a reference position over time. It can be thought of as a measure of how much the system moves over time. The MSD is given by:

:<math>{\rm MSD}\equiv\langle (r-r_0)^2\rangle=\frac{1}{N}\sum_{n=1}^N (r_n(t) - r_n(0))^2</math>

The diffusion coefficient (D) can be calculated from the MSD, using<ref name=":1">O. J. Eder, ''J. Chem. Phys.'', 1977, '''66''', 3866–3870.</ref>:

<math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

The '''VACF''' is given by

<math>C\left(\delta t\right) = \left\langle \mathbf{v}\left(t\right) \cdot \mathbf{v}\left(t+\delta t\right)\right\rangle</math>

and is effectively a measure of how closely related the velocity of a particle is (at time t) to its initial velocity (at time t=0). This correlation is "perturbed" by collisions; at very long times (i.e. when t tends to infinity) we expect the VACF to be zero, as all particles will have collided at least once and their velocities will be uncorrelated.

The diffusion coefficient is proportional to the integral of the VACF<ref name=":1" />:

<math>D = \frac{1}{3}\int_0^\infty \mathrm{d}\delta t \left\langle\mathbf{v}\left(0\right)\cdot\mathbf{v}\left(\delta t\right)\right\rangle=\frac{1}{3}\int_0^\infty C\left(\delta t\right)</math>

 Good, I could reproduce your results with this information. 

==Results & Discussion==
===Equilibration===

Plots of energy, temperature and pressure vs time were obtained for a timestep value of 0.001. These are reproduced in Fig. 2-4. It can be seen that all three plots reach a "plateau" very quickly; equilibration is almost instantaneous. The values oscillate slightly, as a result of the approximations required by the simulation. These oscillations are however very small (note the scale of the y-axis).
{|

|[[File:Ad5215 Ts001 Eng.png|thumb|left|Fig. 2: Energy vs time (ts 0.001)]]
|[[File:Ad5215 Ts001 Temp.png|thumb|left|Fig. 3: Temperature vs time (ts 0.001)]]
|[[File:Ad5215 Ts001 Press.png|thumb|right|Fig. 4: Pressure vs time (ts 0.001)]]

|}

An important feature of these simulations is the timestep. The shorter the timestep the more "accurate" the simulation, but the more computational power this will require. In this case, plots of the total energy vs time were obtained for all five timestep values (Fig. 5).

[[File:Ad515 Allts energy.png|frame|center|Fig. 5: Energy vs time for all timesteps]]

The lowest energy is given by timestep values of 0.001 and 0.0025. These energies are almost identical; the 0.001 energy is lower, but the difference in energies is only 0.005%. In addition to this, for simulating a total time of e.g. 100s, ts = 0.0025 requires 40,000 steps, while ts = 0.001 requires 100,000 steps. Therefore, the 0.0025 timestep is the better choice, as the difference in energies is not large enough to warrant the use of more computational power (as required by the 0.001 timestep). The 0.015 timestep is a poor choice. Not only is the energy the highest of the five, but, unlike in the other four cases, the system does not reach equilibrium and the energy keeps increasing.

Based on this data, further simulations were run using a 0.0025 timestep (unless a different value was required by the lab script).

 Equilibration is not typically discussed in the results section of a scientific paper. Simply "systems were equilibrium for X timesteps/unit with a timestep of Y" would be sufficient. You get the marks for the task however. 

===Densities and the ideal gas law===

Plots of density vs temperature for two different pressures are reproduced in Fig. 6. The density given by the ideal gas law was also calculated and plotted on the same graph.

[[File:Ad5215 densvstemp.png|frame|center|Fig. 6: Plots of density vs temperature comparing experimental and theoretical data]]

It can be seen that the results of the simulation do not match those given by the theoretical approach. This is because the ideal gas law does not take into account any interaction between particles, i.e. it assumes that the gas behaves ideally. In the Lennard-Jones model the particles experience attractive and repulsive forces; the repulsive forces dominate be careful to convey precise meaning: when do repulsive forces dominate? and cause the system to be more diffuse and thus have a lower simulated density.

The discrepancy between theory and simulation increases with increasing pressure because this "pushes" the particles closer together and increases the effect of Lennard-Jones forces. It also increases with decreasing temperature; at low temperatures, the Lennard-Jones forces dominate, while at high temperatures thermal motion is more significant. I understand what you are trying to say here, however a more precise/succinct explanation would have been helpful. For example, rationalising your results in terms of the potential energy surface and available thermal energy ... "at the high temperature limit, LJ particles have enough kinetic energy to easily surmount all kinetic energy barriers in the PES", or something more elegantly worded. 

===Heat capacities===

The variation in heat capacity with temperature is shown in Fig. 7. The system is under the NVE ensemble so we are dealing with the isochoric heat capacity, CV.

[[File:Ad5215 Heatcaps.png|frame|center|Fig. 7: Heat capacity variation with temperature]]

The heat capacity decreases with increasing temperature. This agrees with the formula for heat capacity, which shows inverse proportionality to temperature. However, this is not the full extent of the explanation.

Heat capacity is a measure of how much energy (heat) is required to increase the temperature of a system. At higher temperatures more energetic states become available and the spacing between them decreases - this makes populating higher states easier, leading to a decrease in heat capacity.

In addition to this, increasing the temperature can lead to a phase change and thus to an increase in the degrees of freedom available to the system (e.g. melting causes a solid - rigid, fewer degrees of freedom - to change into a liquid).

===The Radial Distribution Function (RDF)===

[[File:Ad5215 RDFs 3phase.png|frame|right|Fig. 8: RDFs for systems in different phases]]

The radial distribution function for a solid, liquid and gas is reproduced in Fig. 8. The RDF shows how a system is arranged, relative to the position of one particle in the system; effectively, it is a measure of long-range order. A peak corresponds to a shell of atoms around the central particle. The intensity of the peak (effectively its integral) is proportional to the number of atoms within this shell.

All three RDFs (vapour, liquid, solid) show an initial peak, but differ in their behaviour at longer distances. The RDF for a system in the gas phase rapidly reaches a value of 1 and plateaus. This is because a gas is, by its very nature, disordered. Atoms are free to move and they tend to disperse, not arrange themselves in shells. The RDF for a liquid oscillates slightly after the initial peak but also plateaus at 1 after a short distance. The initial peak corresponds to a solvation shell around the central particle. The subsequent smaller peaks show that a liquid has some degree of order - the forces between the particles are strong enough to restrict their movement to a degree.

[[File:Ad5215 small lat.png|thumb|FCC lattice showing first three neighbouring lattice sites for a central atom (light pink)]]

The RDF for the solid system is different to the other two, as it shows long-range order. It does not plateau but instead shows peaks of decreasing intensity. This can be explained by looking at the structure of the solid crystal. This was defined in the simulation as a face-centred cubic (FCC) lattice, shown in Figure 9. The particles are arranged in shells, at distances which depend on the lattice spacing of the crystal. This can be calculated from the lattice density (1.2).

The first three peaks in the RDF plot correspond to the first three neighbouring sites of the central particle, coloured in blue, purple and green respectively. The lattice spacing and the coordination number of each site can be calculated by considering the geometry of the crystal:

*Shell 1 is found at <math> r_1 = \frac{\sqrt{2}}{2}a = 1.056</math> and holds <math>12</math> atoms.
*Shell 2 is found at <math> r_2 = a = 1.494</math> and holds <math>6</math>atoms.
*Shell 3 is found at <math> r_3 = \frac{\sqrt{6}}{2}a = 1.830</math> and holds <math>24</math> atoms.

These values agree with those found by the RDF. The coordination numbers match those calculated from the running integral, which has values of 12.15, 17.98 and 42.3 respectively.

===The diffusion coefficient, D===

Plots of the total mean squared displacement vs time are reproduced in Appendix A. Plots of the VACF vs time are reproduced in Appendix B. Plots of the VACF running integral vs time are reproduced in Appendix B.

Figure 10 below shows the time evolution of the VACF for a solid, a liquid and a gas, as well as for an ideal harmonic oscillator.

[[File:Ad5215 Velocity Autocorrelation Functions small.png|frame|center|Fig. 10: VACF plots for small scale simulations]]

The VACF is effectively a measure of how closely related the velocity of a particle is (at time t) to its initial velocity (at time t=0). This correlation is "perturbed" by collisions or by interactions with other particles; at very long times (i.e. when t tends to infinity) we expect the VACF to be zero and the velocities to be uncorrelated. What does the first minimum indicate? 

The harmonic oscillator shows perfectly oscillatory behaviour, with constant amplitude in time: the velocity goes from an initial state to an uncorrelated one and then back to the initial state. The solid shows similar behaviour: the VACF oscillates about 0 but dampens with time. This is because in a solid the atoms have fixed positions in a lattice; the forces between the particles are strong and these will oscillate in place for a while. The VACF takes much longer to reach zero than in the case of a liquid or a gas. The gas VACF tends slowly to zero; the interactions between particles in a gas are minimal, which means that the velocity at time is not very different from an initial velocity. A liquid is somewhere in-between these two phases: the particles have more freedom of movement than they do in a solid, but the attractive forces are strong enough to cause a perturbation in the velocities; the VACF shows a very slight oscillation, but then quickly dampens to zero.

The diffusion coefficients can be calculated in two different ways, either from the gradient of an MSD plot or from the integral of the VACF. These calculations were performed and the D values are given in Table 1.

{|class="wikitable"
|+ Table 1: Diffusion coefficients <math> D / m^2 s^{-1}</math>
|-
! rowspan="2" | Phase
! colspan="2" | MSD data
! colspan="2" | VACF data
|-
! Small simulation
! Large simulation
! Small simulation
! Large simulation
|-
! Gas
| 2.536
| 2.542
| 3.294
| 3.268
|-
! Gas, linear region
| 3.317
| 3.217
| ---
| ---
|-
! Liquid
| 0.085
| 0.087
| 0.098
| 0.090
|-
! Solid
| 5.825 x 10-7
| 4.391 x 10-4
| -1.845 x 10^-4
| 4.558 x 10^-5
|-
|}

The largest error in the case of the MSD measurements comes from the fact that the gas phase gives rise to a curved plot, which cannot be feasibly fitted to a straight line. This is because particles in a gas will diffuse readily and thus the system will take longer to reach equilibrium (the linear region) than, say, a liquid. A longer simulation that would allow the system to reach equilibrium and then collect a larger amount of data would provide a more accurate result, as in this case the linear region that was used only comprised about 20% of the total data. In the case of a liquid, the MSD function is linear and provides the best fit and thus the most accurate result. The MSD for a solid establishes linear behaviour quickly, as the particles are "fixed". The diffusion coefficient values are very small; this shows that in a solid no diffusion (or almost none) takes place.

In the case of the VACF measurements the largest error comes from using the trapezium rule to compute the integral. The smaller the timestep, the more accurate the measurement - in this case the timestep is relatively small but some error still remains.

The errors in both of these measurements cause the diffusion coefficients to differ slightly. The difference between the MSD- and VACF-calculated diffusion coefficients are 0.67%, 15.3% and 31780% for the gas, liquid and solid phases respectively. The difference in the case of the solid phase is incredibly large but not significant as we have established diffusion is not a significant process for solids. The difference for the liquid phase is quite small and likely comes from the poor fit of the MSD plot and the short duration of the simulation.

===Appendix A: MSD plots===
<gallery>
File:Ad5215 Gas phase MSD, large atom count.png|Gas phase MSD (large scale)
File:Ad5215 Gas phase MSD, large atom count - linear region.png|Gas phase MSD, linear region (large scale)
File:Ad5215 Gas phase MSD, small atom count.png|Gas phase MSD (small scale)
File:Ad5215 as phase MSD, small atom count - linear region.png|Gas phase MSD, linear region (small scale)
File:Ad5215 Liquid phase MSD, large atom count.png|Liquid phase MSD (large scale)
File:Ad5215 Liquid phase MSD, small atom count.png|Liquid phase MSD (small scale
File:Ad5215 Solid phase MSD, large atom count.png|Solid phase MSD (large scale)
File:Ad5215 Solid phase MSD, small atom count.png|Solid phase MSD (small scale)
</gallery>

===Appendix B: VACF running integral plots===
<gallery>
File:Ad5215 Gas phase VACF Integral, large scale.png|Gas phase (large scale)
File:Ad5215 Gas phase VACF Integral, small scale.png|Gas phase (small scale)
File:Ad5215 Liquid phase VACF Integral, large scale.png|Liquid phase (large scale)
File:Ad5215 Liquid phase VACF Integral, small scale.png|Liquid phase (small scale)
File:AD5215 Solid phase VACF Integral, large scale.png|Solid phase (large scale)
File:Ad5215 Solid phase VACF Integral, small scale.png|Solid phase (small scale)
</gallery>

==Conclusion==

This lab provided insight into the thermodynamic properties of systems and how these change with phase, temperature, pressure. Of course, the systems investigated were small, but modeling larger, more complicated systems is possible and could prove useful. A particular domain where this kind of research would be invaluable is, as previously mentioned, molecular gastronomy: understanding phase changes and the properties of liquids, solids and gels can lead to the advancement of this (pseudo)-scientific discipline.

 This conclusion is a bit short: it should include a summary of the results, and perhaps a short outlook. It would also be nice to convey the importance/meaning of your results and conclusions in the context of molecular gastronomy, which you motivated your study with, but have not mentioned since the introduction.

2018-04-16T13:02:59Z

Org12: /* Equilibration */

=Tasks=
==Section 2: Introduction==

'''<big>TASK</big>: Open the file HO.xls. In it, the velocity-Verlet algorithm is used to model the behaviour of a classical harmonic oscillator. Complete the three columns "ANALYTICAL", "ERROR", and "ENERGY": "ANALYTICAL" should contain the value of the classical solution for the position at time <math>t</math>, "ERROR" should contain the ''absolute'' difference between "ANALYTICAL" and the velocity-Verlet solution (i.e. ERROR should always be positive -- make sure you leave the half step rows blank!), and "ENERGY" should contain the total energy of the oscillator for the velocity-Verlet solution. Remember that the position of a classical harmonic oscillator is given by <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math> (the values of <math>A</math>, <math>\omega</math>, and <math>\phi</math> are worked out for you in the sheet).'''

Excel file attached [[:File:AD5215_HO.xls|here]].

'''<big>TASK</big>: For the default timestep value, 0.1, estimate the positions of the maxima in the ERROR column as a function of time. Make a plot showing these values as a function of time, and fit an appropriate function to the data.'''

[[File:Ad5215 error vs time.png]]

'''<big>TASK</big>: Experiment with different values of the timestep. What sort of a timestep do you need to use to ensure that the total energy does not change by more than 1% over the course of your "simulation"? Why do you think it is important to monitor the total energy of a physical system when modelling its behaviour numerically?'''

The change in energy goes down as the timestep value becomes smaller. For a timestep of 0.25 the change in energy is 1.58% while for a timestep of 0.05 the change in energy is 0.06% . The energy change is 1.01% for a timestep of 0.2 . A timestep that is too large could lead to the simulation effectively "missing" any changes in the system that happen on a shorter timescale than that of the timestep. Therefore, it is important to monitor the energy to ensure that the change is not too drastic and we are observing the behaviour of the system closely.

'''<big>TASK:</big> For a single Lennard-Jones interaction, <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>, find the separation, <math>r_0</math>, at which the potential energy is zero. What is the force at this separation? Find the equilibrium separation, <math>r_{eq}</math>, and work out the well depth (<math>\phi\left(r_{eq}\right)</math>). Evaluate the integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math> when <math>\sigma = \epsilon = 1.0</math>.

''Find the separation at which the potential energy is zero''

[[File:Ad5215 lj zero pot.JPG]]

''FInd the force at <math>r=\sigma</math>''

The force is the derivative of potential wrt to distance:

[[File:Ad5215 lj Force.JPG]]

At separation <math>r=\sigma</math> this will be

[[File:Ad5215 force(R).JPG]]

''Equilibrium separation and well depth''<\span>

The equilibrium separation is the separation when <math> \frac{d \phi}{dr} = 0.</math>

[[File:Ad5215 lj equilibrium req.JPG]]

The well depth at this separation is

[[File:Ad5215 ls lj epsilon.JPG]]

 negative epsilon 

''Integrals''

[[File:Ad5215 lj int1.JPG]]

[[File:Ad5215 int2.JPG]]

[[File:Ad5215 int3.JPG]]

'''<big>TASK</big>: Estimate the number of water molecules in 1ml of water under standard conditions. Estimate the volume of <math>10000</math> water molecules under standard conditions.'''

<math> V = 1\ mL, \ \rho = 1\ g/mL, \ M = 18\ g/mol</math>

<math>
m = V \rho = 1\ g
</math>

<math>
N = nN_a = \frac{m}{M} \times N_a = \frac{6.022 \times 10^{23}}{18} = 3.35 \times 10{22}
</math>

N molecules occupy 1 mL. Therefore, the volume of 1 molecule of water will be <math>V_0 = \frac{1}{N} = \frac{1}{3.35 \times 10{22}} = 2.99 \times10^{-23}</math>

The volume of 1000 molecules will be <math>1000 \times V_0 = 2.99 \times10^{-20} </math>

'''<big>TASK</big>: Consider an atom at position <math>\left(0.5, 0.5, 0.5\right)</math> in a cubic simulation box which runs from <math>\left(0, 0, 0\right)</math> to <math>\left(1, 1, 1\right)</math>. In a single timestep, it moves along the vector <math>\left(0.7, 0.6, 0.2\right)</math>. At what point does it end up, ''after the periodic boundary conditions have been applied''?'''

It ends up at position (0.2, 0.1, 0.7).

'''<big>TASK</big>: The Lennard-Jones parameters for argon are <math>\sigma = 0.34\mathrm{nm}, \epsilon\ /\ k_B= 120 \mathrm{K}</math>. If the LJ cutoff is <math>r^* = 3.2</math>, what is it in real units? What is the well depth in <math>\mathrm{kJ\ mol}^{-1}</math>? What is the reduced temperature <math>T^* = 1.5</math> in real units?

<math>r^* = \frac{r}{\sigma}</math>

<math>r = \sigma \times r^*</math>

<math>r = 0.34 \times 3.2 = 1.088\ nm </math>

<math>\frac{\epsilon}{k_B} = 120\ K </math>

<math>\epsilon = 120 \times k_B </math> for 1 particle

For a mole of particles: <math>\epsilon = 120 \times k_B \times N_A </math>

<math> \epsilon = 0.997\, kJ mol^{-1} </math>

<math>T^* = \frac{k_BT}{\epsilon}</math>

<math>T = T^* \times \frac{\epsilon}{k_BT}</math>

<math>T = 1.5 \times 120 </math>

<math>T = 180\ K</math>

==Section 3: Equilibration==
'''<big>TASK</big>: Why do you think giving atoms random starting coordinates causes problems in simulations? Hint: what happens if two atoms happen to be generated close together?'''

Randomly generated positions can lead to two atoms being very close together, which would result in a large repulsive potential. This would then affect any propagation in time of the system, which would lead to undesirable behaviour, especially when using larger timesteps. The system would most likely behave "appropriately" for small enough timesteps, but this would require running longer simulations. This would be less effective; a larger timestep that still results in an accurate simulation is ideal. good! 

'''<big>TASK</big>: Satisfy yourself that this lattice spacing corresponds to a number density of lattice points of <math>0.8</math>. Consider instead a face-centred cubic lattice with a lattice point number density of 1.2. What is the side length of the cubic unit cell?'''

Density: <math> \rho = \frac{N}{V} = \frac{N}{l^3} </math> ---> Length: <math> l=\sqrt[3]{\frac{N}{\rho}} </math>

For a simple cubic lattice with <math> \rho = 0.8 </math>, the length, l, is:

<math> \sqrt[3]{\frac{1}{0.8}}=1.07722 </math>

For a face-centered cubic lattice with <math> \rho = 1.2 </math>, the length, l, is:

<math> \sqrt[3]{\frac{4}{1.2}}=1.4938 </math>

'''<big>TASK</big>: Consider again the face-centred cubic lattice from the previous task. How many atoms would be created by the create_atoms command if you had defined that lattice instead?'''

The box would still contain <math> 10 \times 10 \times 10 = 1000 </math> lattice units. For an FCC there are 4 atoms per lattice unit. Therefore the total number of atoms would be <math> 4 \times 1000 = 4000 </math> atoms.

'''<big>TASK</big>: Using the [http://lammps.sandia.gov/doc/Section_commands.html#cmd_5 LAMMPS manual], find the purpose of the following commands in the input script:'''

<pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

"Mass" sets the mass of a particular type of atom (in this case, the mass of type 1 atoms is 1). "Pair" refers to pair potentials. The lj/cut command computes the 12/6 Lennard-Jones potential, cut sets the cut-off for r. Pair_coeff sets the values for the 2 parameters, sigma and epsilon.

'''<big>TASK</big>: Given that we are specifying <math>\mathbf{x}_i\left(0\right)</math> and <math>\mathbf{v}_i\left(0\right)</math>, which integration algorithm are we going to use?'''

 Velocity Verlet. A simple Verlet algorithm wouldn't require the initial velocity, but would instead require <math>\mathbf{x}_i\left(-\delta t\right)</math>.

'''<big>TASK</big>: 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}
run 100000
</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. In light of this, what do you think the purpose of these lines is? Why not just write:'''
<pre>
timestep 0.001
run 100000
</pre>

The line "variable timestep equal 0.001" defines a variable timestep which is the assigned a value. This allows for the variable to be called later on if needed. This is convenient for the user as it means that if the same variable is required multiple times (in this case, the variable timestep is called twice) changing its value is easier, as this only needs to be done once (in the line defining the variable).

==Section 4: Running simulations under specific conditions==

'''<big>TASK</big>: We need to choose <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>

'''Solve these to determine <math>\gamma</math>.'''

[[File:Ad5215 ljls gamma.JPG]]

'''<big>TASK</big>: Use the [http://lammps.sandia.gov/doc/fix_ave_time.html manual page] to find out the importance of the three numbers ''100 1000 100000''. How often will values of the temperature, etc., be sampled for the average? How many measurements contribute to the average? Looking to the following line, how much time will you simulate?'''

From the manual, the structure of the command is "ave/time Nevery Nrepeat Nfreq".

*Nevery (100) = use input values every this many timesteps (total no. of data points is 100,000 so this will give 1000 values to be averaged in the next step; records an average of 100 values)

*Nrepeat (1000) = no. of times to use input values for calculating averages (i.e. average over 1000 values)

*Nfreq (100,000) = calculate averages every this many timesteps (same no. specified in the "run" command)

The timestep is 0.0025 and the simulation runs for 100,000 steps. Therefore we are simulating a total time of 250 seconds .

==Section 7: Dynamical properties and the diffusion coefficient==

'''<big>TASK</big>: In the theoretical section at the beginning, the equation for the evolution of the position of a 1D harmonic oscillator as a function of time was given. Using this, evaluate the normalised velocity autocorrelation function for a 1D harmonic oscillator (it is analytic!):

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} v\left(t\right)v\left(t + \tau\right)\mathrm{d}t}{\int_{-\infty}^{\infty} v^2\left(t\right)\mathrm{d}t}</math>

For a HO model:
[[File:Ad5215 Vacf HO model.JPG]]

The VACF is given by

[[File:Ad5215 VACF deriv.JPG]]

We evaluate the two integrals separately. Integral 2 is

[[File:Ad5215 Vacf i2.JPG]]

Using [[File:Ad5215 Vacf cos(2x).JPG]] we can write

[[File:Ad5215 Vacf i2 2.JPG]]

Integral 1 is

[[File:Ad5215 Vacf i1 1.JPG]]

Using [[File:Ad5215 Vacf sin(A+B).JPG]] we can write

[[File:Ad5215 Vacf i1 2.JPG]]

We evaluate the last integral separately.

[[File:Ad5215 Vacf i3.JPG]]

Substituting this back we obtain:

[[File:Ad5215 Vacf i1 3.JPG]]

Substituting I1 and I2 into the formula for C we obtain

[[File:Ad5215 Vacf c1.JPG]]

Sin is an odd function, i.e. [[File:Ad5215 Vacf sinodd.JPG]]. Thus

[[File:Ad5215 Vacf c2.JPG]]

 There is definitely a shorter way of doing this, but yes. 

=Report=

==Abstract==

Solid, liquid and gaseous systems were modeled using LAMMPS and a 12-6 Lennard-Jones forcefield. A suitable timestep of 0.0025 was determined. Simulations were run to obtain thermodynamic data (temperatures, pressures, densities, heat capacities). Calculated densities were found to be lower than those predicted by the ideal gas law. The simulated heat capacities showed a trend, decreasing with increasing temperature. The RDF was calculated for systems in all three phases. As expected, the RDF for a solid showed peaks decaying in amplitude. Lattice spacings and coordination numbers for the solid FCC lattice were calculated. The MSD and the VACF were plotted for the same three systems and the diffusion coefficient was calculated for both measurements. The two methods did not result in identical values; still, the difference was only 0.67% for the diffusion coefficients in the gas phase.
 Good, concise abstract that summaries main results. Perhaps 1 sentence of motivation would have been nice. 

==Introduction==

Food constitutes a huge part of our lives, whether we actively think about it or not. Cooking, in some form or other, has been around ever since the first human realised that fire makes raw meat tastier and more convenient to eat. Nowadays, the food industry is enormous and cooking has become a successful blend of art and science. Perhaps the most pertinent example of this is molecular gastronomy, a subdiscipline of food science that seeks to investigate the physical and chemical transformations of ingredients during the cooking process. In their review, Balham et al refer to molecular gastronomy as an "emerging scientific discipline".<ref>P. Barham, L. H. Skibsted, W. L. P. Bredie and J. Risbo, ''Symp.
A Q. J. Mod. Foreign Lit.'', 2010, 2313–2365.</ref>

In fact, molecular gastronomy relies heavily on processes such as gelification and infusion and on materials such as gels, foams and powders. It also uses equipment that is heavily reminiscent of a laboratory - some kitchens are fitted with butane burners, syringes and dehydrators.

Clever presentations and unusual sensations and are piquing the interest of many people; in some places, molecular gastronomy restaurants have become tourist attractions.<ref>D. Tüzünkan and A. Albayrak, ''Procedia
- Soc. Behav. Sci.'', 2015, '''195''', 446–452.</ref> Simulating the thermodynamic properties of systems can provide a better understanding of physical systems and can lead to the growth and development of this industry.

 An interesting topic and motivation! An explicit connection between gastronomy and the results you intend to present/discuss would be nice. For example how is the diffusion coefficient relevant? The introduction of a scientific paper usually includes the background theory, such as in your case, the equations for diffusion coefficient. You have included this in the methodology, which while logical, is not standard practice for most papers/journals. 

==Aims and Objectives==

* become familiar with LAMMPS and how simulating a physical system works
* model the behaviour of systems under the 12-6 Lennard-Jones potential
* calculate the physical properties (temperature, pressure, density, etc) of a system using such simulations
* comparing the results of the simulations to theory (e.g. simulated density vs density given by the ideal gas law)
* calculating the diffusion coefficient for systems in different phases (gas, liquid, solid) by two different methods (from the MSD and from the VACF)

==Methods==
'''General methods'''

A 12-6 Lennard-Jones system was modeled usings LAMMPS and all simulations were run using Imperial's High Performance Computer.
The potential for such a system is given by<ref name=":0">P. Atkins, J. De Paula, ''Physical
Chemistry'', OUP Oxford, 9th edn., 2009.</ref>:

<math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>

All calculations used reduced units:<math>r^* = r/{\sigma}</math>, <math>E^* = E/{\epsilon}</math> and <math>T^* = {k_BT}/{\epsilon}.</math>
The Lennard-Jones parameters <math>\sigma</math> and <math>\epsilon</math> were set to 1.0 for all simulations. The cutoff for Lennard-Jones interactions was set at r*=3 unless otherwise stated. The mass of the atoms was set to 1.0. The temperature, pressure, lattice density and timestep values were varied. For all calculations the velocity Verlet algorithm was employed. The atoms were assigned random velocities within the simulation, while ensuring that the Boltzmann distribution of states is followed.

'''Determining a suitable timestep'''

A simple cubic lattice with a density of 0.8 was defined and the simulation was populated with 1000 atoms (10 x 10 x 10 dimensions). The ensemble was defined as the microcanonical (NVE) ensemble. Five values for the timestep were tested: 0.015, 0.01, 0.0075, 0.0025, 0.001 and each simulation was run for a total time of 100 seconds. Values for the energy, temperature and pressure of the system were recorded at each step.

'''Variation of density with temperature and pressure'''

A simple cubic lattice with a density of 0.8 was defined and the simulation was populated with 3375 atoms (15 x 15 x 15 dimensions). The timestep for all simulations was set to 0.0025. The ensemble was defined as NPT and 10 different thermodynamic states were simulated (two pressures, 2.5 and 3.5, each associated with five temperatures: 3.0, 6.0, 9.0, 12.0 and 15.0). Values for the energy, temperature and pressure of the system were recorded at each step, as well as average values for the density, temperature and pressure of the system at the end of the simulation. Plots of density vs time were obtained, both for the simulated data and for densities predicted by the ideal gas law<ref name=":0" />, <math> PM = \rho RT.</math>

'''Heat capacity calculations'''

A simple cubic lattice with a density of 0.2 was defined and populated with 3375 atoms. The timestep for all simulations was set to 0.0025. An NVT ensemble was simulated for five temperatures: 2.0, 2.2, 2.4, 2.6 and 2.8. Then, a subsequent simulation was run to establish an NVE ensemble and measure the properties of the system. Average values for temperature, energy, volume and heat capacity were calculated. This procedure was repeated for a simple cubic lattice with a density of 0.8. An example input script can be found [[:File: Example_script_heatcap_ad5215.in|here]].

The heat capacity of a system is given by<ref name=":0" />:

<math>C_V = \frac{\partial E}{\partial T} = \frac{\mathrm{Var}\left[E\right]}{k_B T^2} = N^2\frac{\left\langle E^2\right\rangle - \left\langle E\right\rangle^2}{k_B T^2}</math>

where E is the internal energy and

<math>{\mathrm{Var}\left[E\right]}={\left\langle E^2\right\rangle - \left\langle E\right\rangle^2}</math>

is the variance in internal energy. The <math>N^2</math> term is required because LAMMPS automatically outputs the energy '''per atom''', not the '''total''' energy.

[[File:Ad5215 LJfluid phase diag.JPG|thumb|Fig. 1: Phase diagram for the Lennard-Jones fluid]]
'''RDF calculations'''

Three systems (a liquid, a solid and a gas) were modeled and populated with 3375 atoms each. Temperature and density values were taken from the Lennard-Jones fluid phase diagram<ref>J. P. Hansen and L.
Verlet, ''Phys. Rev.'', 1969, '''184''', 151–161.</ref> reproduced in Fig. 1. These were defined as:

*solid: fcc lattice, temperature 1.2, density 1.2;
*liquid: sc lattice, temperature 1.2 , density 0.8;
*vapour: sc lattice, temperature 1.2, density 0.05.

The ensemble was defined as NVT. The timestep for all simulations was set to 0.002. The trajectories of the atoms were recorded and VMD was used to calculate the radial distribution function and its integral from these trajectories. The data was then analysed using Python.

'''Diffusion coefficient calculations'''

The same three systems as above were modeled, this time with 8000 atoms each. The timestep was set to 0.002 and each simulation was run for 5000 steps. The Lennard-Jones cutoff was set to 3.2. The ensemble was defined as NVT. The mean squared displacement (MSD) and the velocity autocorrelation function (VACF) at each step were calculated for all systems. The data was analysed using Python. The MSD plots were fitted to a straight line and the gradient was used to calculate the diffusion coefficient. The VACF integrals were plotted as a function of time and, again, used to calculate the diffusion coefficient. The same data analysis was conducted using supplied data which modeled larger systems.

The '''MSD''' is a measure of the deviation of the position of a particle with respect to a reference position over time. It can be thought of as a measure of how much the system moves over time. The MSD is given by:

:<math>{\rm MSD}\equiv\langle (r-r_0)^2\rangle=\frac{1}{N}\sum_{n=1}^N (r_n(t) - r_n(0))^2</math>

The diffusion coefficient (D) can be calculated from the MSD, using<ref name=":1">O. J. Eder, ''J. Chem. Phys.'', 1977, '''66''', 3866–3870.</ref>:

<math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

The '''VACF''' is given by

<math>C\left(\delta t\right) = \left\langle \mathbf{v}\left(t\right) \cdot \mathbf{v}\left(t+\delta t\right)\right\rangle</math>

and is effectively a measure of how closely related the velocity of a particle is (at time t) to its initial velocity (at time t=0). This correlation is "perturbed" by collisions; at very long times (i.e. when t tends to infinity) we expect the VACF to be zero, as all particles will have collided at least once and their velocities will be uncorrelated.

The diffusion coefficient is proportional to the integral of the VACF<ref name=":1" />:

<math>D = \frac{1}{3}\int_0^\infty \mathrm{d}\delta t \left\langle\mathbf{v}\left(0\right)\cdot\mathbf{v}\left(\delta t\right)\right\rangle=\frac{1}{3}\int_0^\infty C\left(\delta t\right)</math>

 Good, I could reproduce your results with this information. 

==Results & Discussion==
===Equilibration===

Plots of energy, temperature and pressure vs time were obtained for a timestep value of 0.001. These are reproduced in Fig. 2-4. It can be seen that all three plots reach a "plateau" very quickly; equilibration is almost instantaneous. The values oscillate slightly, as a result of the approximations required by the simulation. These oscillations are however very small (note the scale of the y-axis).
{|

|[[File:Ad5215 Ts001 Eng.png|thumb|left|Fig. 2: Energy vs time (ts 0.001)]]
|[[File:Ad5215 Ts001 Temp.png|thumb|left|Fig. 3: Temperature vs time (ts 0.001)]]
|[[File:Ad5215 Ts001 Press.png|thumb|right|Fig. 4: Pressure vs time (ts 0.001)]]

|}

An important feature of these simulations is the timestep. The shorter the timestep the more "accurate" the simulation, but the more computational power this will require. In this case, plots of the total energy vs time were obtained for all five timestep values (Fig. 5).

[[File:Ad515 Allts energy.png|frame|center|Fig. 5: Energy vs time for all timesteps]]

The lowest energy is given by timestep values of 0.001 and 0.0025. These energies are almost identical; the 0.001 energy is lower, but the difference in energies is only 0.005%. In addition to this, for simulating a total time of e.g. 100s, ts = 0.0025 requires 40,000 steps, while ts = 0.001 requires 100,000 steps. Therefore, the 0.0025 timestep is the better choice, as the difference in energies is not large enough to warrant the use of more computational power (as required by the 0.001 timestep). The 0.015 timestep is a poor choice. Not only is the energy the highest of the five, but, unlike in the other four cases, the system does not reach equilibrium and the energy keeps increasing.

Based on this data, further simulations were run using a 0.0025 timestep (unless a different value was required by the lab script).

 Equilibration is not typically discussed in the results section of a scientific paper. Simply "systems were equilibrium for X timesteps/unit with a timestep of Y" would be sufficient. You get the marks for the task however. 

===Densities and the ideal gas law===

Plots of density vs temperature for two different pressures are reproduced in Fig. 6. The density given by the ideal gas law was also calculated and plotted on the same graph.

[[File:Ad5215 densvstemp.png|frame|center|Fig. 6: Plots of density vs temperature comparing experimental and theoretical data]]

It can be seen that the results of the simulation do not match those given by the theoretical approach. This is because the ideal gas law does not take into account any interaction between particles, i.e. it assumes that the gas behaves ideally. In the Lennard-Jones model the particles experience attractive and repulsive forces; the repulsive forces dominate and cause the system to be more diffuse and thus have a lower simulated density.

The discrepancy between theory and simulation increases with increasing pressure because this "pushes" the particles closer together and increases the effect of Lennard-Jones forces. It also increases with decreasing temperature; at low temperatures, the Lennard-Jones forces dominate, while at high temperatures thermal motion is more significant.

===Heat capacities===

The variation in heat capacity with temperature is shown in Fig. 7. The system is under the NVE ensemble so we are dealing with the isochoric heat capacity, CV.

[[File:Ad5215 Heatcaps.png|frame|center|Fig. 7: Heat capacity variation with temperature]]

The heat capacity decreases with increasing temperature. This agrees with the formula for heat capacity, which shows inverse proportionality to temperature. However, this is not the full extent of the explanation.

Heat capacity is a measure of how much energy (heat) is required to increase the temperature of a system. At higher temperatures more energetic states become available and the spacing between them decreases - this makes populating higher states easier, leading to a decrease in heat capacity.

In addition to this, increasing the temperature can lead to a phase change and thus to an increase in the degrees of freedom available to the system (e.g. melting causes a solid - rigid, fewer degrees of freedom - to change into a liquid).

===The Radial Distribution Function (RDF)===

[[File:Ad5215 RDFs 3phase.png|frame|right|Fig. 8: RDFs for systems in different phases]]

The radial distribution function for a solid, liquid and gas is reproduced in Fig. 8. The RDF shows how a system is arranged, relative to the position of one particle in the system; effectively, it is a measure of long-range order. A peak corresponds to a shell of atoms around the central particle. The intensity of the peak (effectively its integral) is proportional to the number of atoms within this shell.

All three RDFs (vapour, liquid, solid) show an initial peak, but differ in their behaviour at longer distances. The RDF for a system in the gas phase rapidly reaches a value of 1 and plateaus. This is because a gas is, by its very nature, disordered. Atoms are free to move and they tend to disperse, not arrange themselves in shells. The RDF for a liquid oscillates slightly after the initial peak but also plateaus at 1 after a short distance. The initial peak corresponds to a solvation shell around the central particle. The subsequent smaller peaks show that a liquid has some degree of order - the forces between the particles are strong enough to restrict their movement to a degree.

[[File:Ad5215 small lat.png|thumb|FCC lattice showing first three neighbouring lattice sites for a central atom (light pink)]]

The RDF for the solid system is different to the other two, as it shows long-range order. It does not plateau but instead shows peaks of decreasing intensity. This can be explained by looking at the structure of the solid crystal. This was defined in the simulation as a face-centred cubic (FCC) lattice, shown in Figure 9. The particles are arranged in shells, at distances which depend on the lattice spacing of the crystal. This can be calculated from the lattice density (1.2).

The first three peaks in the RDF plot correspond to the first three neighbouring sites of the central particle, coloured in blue, purple and green respectively. The lattice spacing and the coordination number of each site can be calculated by considering the geometry of the crystal:

*Shell 1 is found at <math> r_1 = \frac{\sqrt{2}}{2}a = 1.056</math> and holds <math>12</math> atoms.
*Shell 2 is found at <math> r_2 = a = 1.494</math> and holds <math>6</math>atoms.
*Shell 3 is found at <math> r_3 = \frac{\sqrt{6}}{2}a = 1.830</math> and holds <math>24</math> atoms.

These values agree with those found by the RDF. The coordination numbers match those calculated from the running integral, which has values of 12.15, 17.98 and 42.3 respectively.

===The diffusion coefficient, D===

Plots of the total mean squared displacement vs time are reproduced in Appendix A. Plots of the VACF vs time are reproduced in Appendix B. Plots of the VACF running integral vs time are reproduced in Appendix B.

Figure 10 below shows the time evolution of the VACF for a solid, a liquid and a gas, as well as for an ideal harmonic oscillator.

[[File:Ad5215 Velocity Autocorrelation Functions small.png|frame|center|Fig. 10: VACF plots for small scale simulations]]

The VACF is effectively a measure of how closely related the velocity of a particle is (at time t) to its initial velocity (at time t=0). This correlation is "perturbed" by collisions or by interactions with other particles; at very long times (i.e. when t tends to infinity) we expect the VACF to be zero and the velocities to be uncorrelated.

The harmonic oscillator shows perfectly oscillatory behaviour, with constant amplitude in time: the velocity goes from an initial state to an uncorrelated one and then back to the initial state. The solid shows similar behaviour: the VACF oscillates about 0 but dampens with time. This is because in a solid the atoms have fixed positions in a lattice; the forces between the particles are strong and these will oscillate in place for a while. The VACF takes much longer to reach zero than in the case of a liquid or a gas. The gas VACF tends slowly to zero; the interactions between particles in a gas are minimal, which means that the velocity at time is not very different from an initial velocity. A liquid is somewhere in-between these two phases: the particles have more freedom of movement than they do in a solid, but the attractive forces are strong enough to cause a perturbation in the velocities; the VACF shows a very slight oscillation, but then quickly dampens to zero.

The diffusion coefficients can be calculated in two different ways, either from the gradient of an MSD plot or from the integral of the VACF. These calculations were performed and the D values are given in Table 1.

{|class="wikitable"
|+ Table 1: Diffusion coefficients <math> D / m^2 s^{-1}</math>
|-
! rowspan="2" | Phase
! colspan="2" | MSD data
! colspan="2" | VACF data
|-
! Small simulation
! Large simulation
! Small simulation
! Large simulation
|-
! Gas
| 2.536
| 2.542
| 3.294
| 3.268
|-
! Gas, linear region
| 3.317
| 3.217
| ---
| ---
|-
! Liquid
| 0.085
| 0.087
| 0.098
| 0.090
|-
! Solid
| 5.825 x 10-7
| 4.391 x 10-4
| -1.845 x 10^-4
| 4.558 x 10^-5
|-
|}

The largest error in the case of the MSD measurements comes from the fact that the gas phase gives rise to a curved plot, which cannot be feasibly fitted to a straight line. This is because particles in a gas will diffuse readily and thus the system will take longer to reach equilibrium (the linear region) than, say, a liquid. A longer simulation that would allow the system to reach equilibrium and then collect a larger amount of data would provide a more accurate result, as in this case the linear region that was used only comprised about 20% of the total data. In the case of a liquid, the MSD function is linear and provides the best fit and thus the most accurate result. The MSD for a solid establishes linear behaviour quickly, as the particles are "fixed". The diffusion coefficient values are very small; this shows that in a solid no diffusion (or almost none) takes place.

In the case of the VACF measurements the largest error comes from using the trapezium rule to compute the integral. The smaller the timestep, the more accurate the measurement - in this case the timestep is relatively small but some error still remains.

The errors in both of these measurements cause the diffusion coefficients to differ slightly. The difference between the MSD- and VACF-calculated diffusion coefficients are 0.67%, 15.3% and 31780% for the gas, liquid and solid phases respectively. The difference in the case of the solid phase is incredibly large but not significant as we have established diffusion is not a significant process for solids. The difference for the liquid phase is quite small and likely comes from the poor fit of the MSD plot and the short duration of the simulation.

===Appendix A: MSD plots===
<gallery>
File:Ad5215 Gas phase MSD, large atom count.png|Gas phase MSD (large scale)
File:Ad5215 Gas phase MSD, large atom count - linear region.png|Gas phase MSD, linear region (large scale)
File:Ad5215 Gas phase MSD, small atom count.png|Gas phase MSD (small scale)
File:Ad5215 as phase MSD, small atom count - linear region.png|Gas phase MSD, linear region (small scale)
File:Ad5215 Liquid phase MSD, large atom count.png|Liquid phase MSD (large scale)
File:Ad5215 Liquid phase MSD, small atom count.png|Liquid phase MSD (small scale
File:Ad5215 Solid phase MSD, large atom count.png|Solid phase MSD (large scale)
File:Ad5215 Solid phase MSD, small atom count.png|Solid phase MSD (small scale)
</gallery>

===Appendix B: VACF running integral plots===
<gallery>
File:Ad5215 Gas phase VACF Integral, large scale.png|Gas phase (large scale)
File:Ad5215 Gas phase VACF Integral, small scale.png|Gas phase (small scale)
File:Ad5215 Liquid phase VACF Integral, large scale.png|Liquid phase (large scale)
File:Ad5215 Liquid phase VACF Integral, small scale.png|Liquid phase (small scale)
File:AD5215 Solid phase VACF Integral, large scale.png|Solid phase (large scale)
File:Ad5215 Solid phase VACF Integral, small scale.png|Solid phase (small scale)
</gallery>

==Conclusion==

This lab provided insight into the thermodynamic properties of systems and how these change with phase, temperature, pressure. Of course, the systems investigated were small, but modeling larger, more complicated systems is possible and could prove useful. A particular domain where this kind of research would be invaluable is, as previously mentioned, molecular gastronomy: understanding phase changes and the properties of liquids, solids and gels can lead to the advancement of this (pseudo)-scientific discipline.