View source for Talk:Mod:thk213

'''<span style="color:red">JC: General comments: Good report which demonstrates a good understanding of the experiment and your explanations of your results and graphs are clear and well thought out.  Look into the background theory a bit more to get a deeper understanding of what each step involves.</span>'''

=== Introduction ===

In this assignment, simulation of liquid is conducted using the Visual Molecular Dynamics (VMD) to understand the molecular motion described by equation of motions in Newton classical mechanics and properties of liquid i.e. heat capacity and compared with different models including radial distribution functions (RDF), mean squared velocity displacement (MSD) and velocity autocorrelation function (VACF). 

'''<span style="color:red">JC: You use LAMMPS to perform the simulation, VMD is only used to visualise the results.</span>'''

== Introduction to Molecular Dynamics Simulation ==

'''<big>TASK 1</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).

Using the classical harmonic oscillator equation, the positions of the atoms are calculated using <math>A=1</math>, <math>\omega=1</math>, and <math>\phi=0</math>. The "Error" column is completed by finding the difference between the position values obtained from velocity-Verlet and the classical harmonic oscillator equation. The "Energy" is obtained by substituting <math>m=1</math> and corresponding velocities using <math> (1) </math>.

<center><math>KE=\frac{1}{2}\ m_i{v_i}^2 \   (1)</math></center>

By inspection from Figure 1, maxima occur at a specific range. Using "Conditional Formatting" function in Excel, the positions of the maxima are determined accurately. The positions of the maxima are plotted as a function of time to generate a straight line with <math> y=0.004x - 7\times 10^{-5}</math> as shown in Figure 2. We can see the error is getting larger and larger over time. Maximum at each range increases. The error reaches a maximum point and fall back to zero and this repeats continuously. When the simulated and velocity-Verlet oscillate at the same time, the error is maximised and this is the peak at each range. 

'''<span style="color:red">JC: You should calculate the total energy, not just the kinetic energy.</span>'''

[[File:Figure 1 - A plot showing the difference of solutions obtained between the velocity-Verlet and classical harmonic oscillator equation.png|thumb|right|500 px|'''Figure 1 - A plot showing the difference of solutions obtained between the velocity-Verlet and classical harmonic oscillator equation''']]

'''<big>TASK 2</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:Figure 1 - A plot showing the maximum difference between the classical harmonic oscillator and the velocity Verlet solutions against time.png|thumb|right|500 px|'''Figure 2 - A plot showing the maximum difference between the classical harmonic oscillator and the velocity Verlet solutions against time''']]

'''Table 1'''

{| {{table}}border="1" cellpadding="4" cellspacing="0" style="border:#c9c9c9 1px solid; margin: 1em 1em 1em 0; border-collapse: collapse;" 
| align="center" style="background:#f0f0f0;"|'''Time'''
| align="center" style="background:#f0f0f0;"|'''Maxima'''
|-
| 2.00||7.58E-4
|-
| 4.90||2.01E-3
|-
| 8.00||3.30E-3
|-
| 11.10||4.60E-3
|-
| 14.20||5.91E-3
|-
| 
|}

From this plot, we can see the maximum of error increases linearly over time. 

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

As the timesteps increases from 0.1, the change in magnitude of total energy decreases significantly. A small change in timestep results a big change in total energy. At timestep = 1.98, the total energy change over the course of the simulation is within 1%,shown in Figure 3. However, when the timestep further raises to 2, the total energy is effectively constant which shows as a straight line and the oscillating behavior observed disappears. It is crucial to monitor and ensure the total energy of a system to approach a nearly constant based on the Law of Conservation of Energy, a fundamental principle built on the Noether's theorem. Total energy of the system must remain constant as time changes.

'''<span style="color:red">JC: Good choice of timestep.  Why did the oscillatory behaviour disappear when timestep is increased to 2?</span>'''

[[File:The total energy at timestep 1.98.png|thumb|right|500 px|'''Figure 3 - A plot showing the total energy of the system at timestep 1.98''']]

'''<big>TASK 4:</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 <math>r_{0}</math> and the force at this separation'''

Potential energy is zero at <math> r = r_{0} </math>

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

<math>\frac{\sigma^{12}}{r_{0}^{12}}=\frac{\sigma^6}{r_{0}^6}</math>

<math>r_{0}^{6}= \sigma^{6}</math>

<math>r_{0}= \sigma</math>

To find the force exerted at this separation, we need to understand that force can be expressed as followings in terms of potential:

<math>\mathbf{F}_i = - \frac{\mathrm{d}U\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_i}</math>

So we differentiate the <math>\phi\left(r_{0}\right)</math> function with respect to <math>r</math> and we substitute <math>r_{0}= \sigma</math> to obtain the force:

<math>\mathbf{F} = - \frac{\mathrm{d}U}{\mathrm{d}\mathbf{r}}_{(r_{0}=\sigma)} = - 4 \epsilon\left(\frac{12\sigma^{12}}{\sigma^{13}} + \frac{6\sigma^6}{\sigma^7}\right) = - 4 \epsilon\left(\frac{12}{\sigma} + \frac{6}{\sigma}\right)</math>

<math>\mathbf{F} = 4\epsilon\left(\frac{-6}{\sigma}\right)= \frac{-24\epsilon}{\sigma}</math>

'''Find the equilibrium separation, <math>r_{eq}</math>, and work out the well depth (<math>\phi\left(r_{eq}\right)</math>)'''

Equilibrium separation is achieved when the force is at minimum (i.e. <math>\mathbf{F}_i = - \frac{\mathrm{d}U}{\mathrm{d}\mathbf{r}} = 0 </math>)

<math>\frac{-12\sigma^{12}}{r_{eq}^{13}} + \frac{6\sigma^6}{r_{eq}^7} = 0 </math>

<math>\frac{12\sigma^{12}}{r_{eq}^{13}} = \frac{6\sigma^6}{r_{eq}^7} </math>

<math> r_{eq}^{6} = 2\sigma^{-6} </math>

<math> r_{eq} = 2^{1/6}\sigma </math>


At <math> r_{eq} = 2^{1/6}\sigma </math>, <math>\phi\left(r_{eq}\right) = - \epsilon </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>'''

<math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r = -0.0248 </math>

<math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r = -0.00818 </math>

<math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r = -0.00329 </math>

<math> \sigma</math> is the distance at which the intermolecular potential between the two particles is zero. When the distance increases, the values become less negative.

'''<span style="color:red">JC: The maths is correct, except that the force at r<sub>0</sub> which has the right magnitude, but should be positive rather than negative.</span>'''

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

Density of water = <math> 1 \ g \ cm^{-1} </math>

Mass of water in 1 mL of water = <math> 1 \ g </math>

No. of moles of water molecule in 1 mL of water = <math> 0.0556 \ mol </math>

No. of water molecules in 1 mL of water = <math>3.34\times 10^{22}</math>

Using the ratio, the volume of 10000 water molecules can be found directly.

<math>\frac{10000}{V} = \frac{3.34\times 10^{22}}{1} </math>

<math> V = \frac{10000}{3.34\times 10^{22}} = 2.99\times 10^{-19} mL </math>

'''<span style="color:red">JC: Both calculations are correct, remember density should have units of g cm<sup>-3</sup>.</span>'''

'''<big>TASK 6</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''?'''.

The new position of the atom is <math>\left(1.2, 1.1, 0.7\right)</math>. This is beyond the volume of a cubic box. After the periodic boundary conditions have been applied, the atom ends up at <math>\left(0.2, 0.1, 0.7\right)</math>.

'''<span style="color:red">JC: Correct.</span>'''

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

'''Find <math>r</math> in real units and the corresponding well depth in <math>\mathrm{kJ\ mol}^{-1}</math>'''

<math> 3.2 = \frac{r}{0.34}</math>

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

<math> \epsilon = k_{B}\times 120 = 1.38\times 10^{-23}\times 120 = 1.66\times 10^{-21} J </math>

<math>\phi\left(r_{0}\right) = 4\epsilon \left(\frac{\sigma^{12}}{r_{0}^{12}} - \frac{\sigma^6}{r_{0}^6}\right)\times\N_{A} </math>

<math>\phi\left(r_{0}\right) = -3.09\times 10^{-5} \ kJ \ mol^{-1} </math>

'''Find the reduced temperature <math>T^* = 1.5</math> in real units'''

<math> 1.5 = \frac{k_{B}\times T}{\epsilon}</math>

<math> 1.5 = 8.33\times 10^{-3}\times T</math>

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

'''<span style="color:red">JC: Values of r and T are correct.  You have shown above that the well depth is given by <math>-\epsilon</math>, what value of r have you used here as <math>r_{0}</math>?</span>'''

== Equilibration ==

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

Atoms in liquid are close together with no regular arrangement. To ensure we are simulating atoms in liquid, we need to specify the starting coordinates of the atoms. If the coordinates of two atoms generated are already too close together, the electrons in the orbital repel each other substantially and push the two atoms apart at a constant distance. This models the atom arrangement similar to a solid,which atoms are locked into place and has a long range order. The accumulation of large repulsive forces can crash down the simulation later on. If the coordinates of two atoms generated are too far, the atoms are well separated, behaving like a gas without any order.

'''<span style="color:red">JC: Correct, however, these simulations are purely classical so the repulsion is just due to the potential, not due to any quantum mechanical effects.</span>'''


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

Number of lattice point in simple cubic structure = <math>  1 </math>

<math> Number density = \frac{Number \ of \ lattice \ point}{Volume \ of \ the \ cubic \ unit \ cell} = \frac{1}{1.0072^{3}} = 0.8 </math>

Number of lattice points in FCC structure = <math> 4 </math>

Assume side length of the cubic unit cell = <math> a </math>

<math> 1.2 = \frac{4}{a^{3}} </math>

<math> a = 1.494 \ (in \ reduced \ units) </math>

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

As mentioned before, there are 4 lattice points in face-centred cubic lattice.

There are 1000 unit cells therefore 4000 atoms are created using the create_atoms command.

'''<span style="color:red">JC: Correct.</span>'''


'''<big>TASK 11</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>

<pre>mass 1 1.0</pre>
This is a mass command, indicating the mass of atom type 1 in the lattice box is equal to 1.0. The 1 in the middle specifies the atom type while 1.0 specifies the mass value of this atom type. It does not specify the style in this line. But we are dealing with lj style in this case so I assume the style is lj and the quantity is unitless.

<pre> pair_style lj/cut 3.0 </pre>
This command states that we only consider Lennard Jones Potential (LJP) of the atom pairs within 3 distance units. In this case, lj/cut 3.0 computes the global cut off of the LJP is 3. The potential of atom pairs greater than 3 distance units is zero. 

<pre> pair_coeff * * 1.0 1.0 </pre>
This command computes the pairwise force field coefficients of pairs of atom types. As the style we are working is LJP, the leading 1.0 represents the <math> \epsilon </math> in LJP is equal to 1 while the latter represents the <math> \sigma </math> in LJP is equal to 1 as well for all atom types from 1 to n (inclusive) and from n to N (inclusive) where n is an atom type and N is the number of atom types.<math> \epsilon </math> is the well depth which measures the extent of two atoms attracting each other; <math> \sigma </math> is the distance at which the potential between two atoms is zero.

'''<span style="color:red">JC: Good description of these commands, what's the advantage of using a Lennard Jones potential with a cut off?</span>'''

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

We are going to use the velocity-Verlet algorithm as shown in <math> (2) </math>.

<center><math>\mathbf{v}_i\left(t + \delta t\right) = \mathbf{v}_i\left(t + \frac{1}{2}\delta t\right) + \frac{1}{2}\mathbf{a}_i\left(t + \delta t\right)\delta t \ \ (2)</math></center>

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

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

Both yield the same outcome but the first method provides a more time-saving way if we need to change the number of timestep or the amount of time we want to run for different simulations. By doing so, instead of changing the timesteps or the length of simulation time line by line which incurs mistakes, we just need to change the value of the timestep for the timestep in line 2 and every lines with variable depending on this value will be changed accordingly, including the n_step variable as we see in line 3. This is also a convenient way to check if there is any problem in the script.

'''<span style="color:red">JC: Good explanation.</span>'''

'''<big>TASK 13</big>: make plots of the energy, temperature, and pressure, against time for the 0.001 timestep experiment (attach a picture to your report). Does the simulation reach equilibrium? How long does this take? When you have done this, make a single plot which shows the energy versus time for all of the timesteps (again, attach a picture to your report). Choosing a timestep is a balancing act: the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality; short timesteps, however, mean that the same number of simulation steps cover a shorter amount of actual time, and this is very unhelpful if the process you want to study requires observation over a long time. Of the five timesteps that you used, which is the largest to give acceptable results? Which one of the five is a particularly bad choice? Why?

As soon as the energy levels off at a fairly stable value approximately at t = 0.38, equilibrium is attained. As we see from Figure 5 and 6, average temperature and pressure also settle down at similar times. A plot of energy as a function of time for all the timesteps is plotted as Figure 7. As shown in Figure 7, timestep 0.001 and 0.0025 yield similar outcomes with minimized fluctuations of total energy, indicating the total energy is well conserved. For timesteps = 0.0075 and 0.01, the energy still levels off but at a value higher than the ones conducted with timestep 0.001 and 0.0025 with greater fluctuations. Timestep 0.015 is a bad choice as the simulation has completely destablized and the average energy does not level off. Equilibration state is not reached using this timestep. The total energy reaches at a minimum value at instant time but then oscillates with much bigger amplitude over time. The dynamic breaks down after reaching the maximum limit of timestep perhaps at t = 0.01. The two atoms are being too close together with this timestep, experiencing significant repulsive force therefore inducing significant errors in the simulation. The maximum timestep value is restricted by the higher order time derivatives neglected from the Taylor expansion in the velocity-Verlet algorithm used in this simulation.

[[File:Figure 4 - Plot showing the total energy of the system as a function of time for timestep 0.001 experiment.png|thumb|500 px|right|'''Figure 4 - A plot showing the total energy of the system as a function of time for timestep 0.001 experiment''']]

[[File:Figure 5 - Plot showing the temperature of the system as a function of time for timestep 0.001 experiment.png|thumb|500 px|right|'''Figure 5 - A plot showing the temperature of the system as a function of time for timestep 0.001 experiment''']]

[[File:Figure 6 - Plot showing the pressure as a function of time for timestep 0.001 experiment.png|thumb|500 px|right|'''Figure 6 - A plot showing the pressure as a function of time for timestep 0.001 experiment''']]

[[File:Figure 7 - A plot of total energy as a function of time with different timesteps.png|thumb|700 px|centre|'''Figure 7 - A plot of total energy as a function of time with different timesteps''']]

'''<span style="color:red">JC: Energy should be independent of timestep, so the timestep should be chosen to ensure that this is the case, choosing 0.0025 is a good idea.</span>'''

== Running simulations under specific conditions ==

'''<big>TASK 14</big>: 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 &mdash; 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. Submit these ten jobs to the HPC portal. While you wait for them to finish, you should read the next section.'''

Timestep 0.0025 is chosen as this is the largest acceptable timestep according to the results in the previous section which gives minimal fluctuations over time. Temperatures chosen are 1.7,1.9,2.1,2.3 and 2.5 while the pressure are 2.5 and 2.7. 

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

<math>T = \mathfrak{T}</math>

After rearranging the terms, 

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

<math> \gamma = (\frac {\mathfrak{T}}{T})^{1/2} </math>

'''<span style="color:red">JC: The result is correct, but the steps in your derivation are not very clear, especially writing <math>T = \mathfrak{T}</math> which is not correct.</span>'''

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

The three number specifies on what timesteps the input values are used to find the final average. 100 refers to Nevery, which means the input values are used every 100 timesteps; 1000 refers to Nrepeat, which means 1000 steps are used to find the final average on 100000 timestep; 100000 refers to Nfreq, which means the system will calculate the final average on 100000 timestep. So for our case, the values on timesteps 100,200,300,400,500,600,700,800,900,1000 will be used to compute the final average on timestep 100000.

'''<span style="color:red">JC: You are not just using timesteps up to 1000 for the averaging.  The average will be calculated on 1000 points each 100 timesteps apart: 100, 200, 300, ..., 99900, 100000.</span>'''

'''<big>TASK 17</big>: 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>. Use software of your choice to plot the density as a function of temperature for both of the pressures that you simulated.  Your graph(s) should include error bars in both the x and y directions. You should also include a line corresponding to the density predicted by the ideal gas law at that pressure. Is your simulated density lower or higher? Justify this. Does the discrepancy increase or decrease with pressure?'''

Equation <math> (3) </math> and <math> (4) </math> show the ideal gas equation and van der Waals equation respectively.

<center><math> PV = Nk_BT \ \ (3)</math></center>

<center><math>\left(P + \frac{an^2}{v^2}\right)\left(V-nb\right) = nRT \ \ (4)</math></center>

We are dealing with reduced units in the simulation therefore <math> k_B = 1 </math>

<math> \frac {P}{T} = \frac {N}{V} </math>

By substituting the corresponding temperature and pressure, pressure divided by the temperature is equal to the density.

Using the equations of state plotted in Figure 8, we can see how well our simulated results fit the ideal gas law equation. Ideal gas law model assumes all collisions between atoms (point masses) are perfectly elastic and neglects any attractive or repulsive interactions between atoms. As we can see, at constant pressure, density decreases as temperature increases. At higher temperature, atoms have higher kinetic energy and they can spread out more and move master, resulting in a higher volume and a lower density. The simulated density is half of the density predicted by ideal gas law. Error bars are fit into the simulated graphs only for both x and y directions. Under the ideal gas law conditions, the increase in volume is much less as the volume occupied by particles is negligible in this case compared to the volume of the box, resulting a higher density. The discrepancy increases with higher pressure and lower temperature. At higher pressure, atoms are closer together and volume occupied is larger than expected value predicted from the ideal gas equation.

van der Waals equation gives a more realistic and better prediction than ideal gas equation of how the density vary as temperature changes because it includes the two correction terms in the equation. The <math> \frac {an^2}{v^2} </math> term corrects the attraction forces between atoms while the <math> -nb </math> takes into account of the excluded volume. a and b are van der Waals coefficients subject to a specific type of atom. If the attraction is stronger, a will be larger. 

[[File:Figure 8 - A plot showing the density as a function of temperature at 2 different pressures at timestep = 0.0025.png|thumb|center|700 px|'''Figure 8 - A plot showing the density as a function of temperature at 2 different pressures at timestep = 0.0025''']]

'''<span style="color:red">JC: Good explanation, the van der Waal's equation would give a better fit to the simulated data than the ideal gas law.  Talking about collisions in the ideal gas could be a bit misleading though as there are no interactions between particles.</span>'''

==Calculating heat capacities using statistical physics==

'''<big>TASK 18</big>: As in the last section, you need to run simulations at ten phase points. In this section, we will be in density-temperature <math>\left(\rho^*, T^*\right)</math> phase space, rather than pressure-temperature phase space. The two densities required at <math>0.2</math> and <math>0.8</math>, and the temperature range is <math>2.0, 2.2, 2.4, 2.6, 2.8</math>. Plot <math>C_V/V</math> as a function of temperature, where <math>V</math> is the volume of the simulation cell, for both of your densities (on the same graph). Is the trend the one you would expect? Attach an example of one of your input scripts to your report.'''

Heat capacity is given by <math> (5) </math>, where <math>N</math> is the number of atoms,

<center><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} \ \ (5) </math></center>

10 simulations are run in canonical ensemble (NVT). Before generating the average temperature and the heat capacity, we need to prepare a system at equilibrium by melting a crystal. This is done by modifying the npt.in file we used in the previous section and change it to canonical ensemble. An example of the script is attached in the report. As we have set the constraint at constant volume, pressure is not allowed to change. After switching off the thermostat using the LAMMPS command, we run the simulations in 100000 steps. The results of the simulation are shown in Figure 9.

Here is the input script for the simulation of <math> \rho </math> = 0.2 and <math> T = 2.8 </math>:

<pre>

### DEFINE SIMULATION BOX GEOMETRY ###
lattice sc 0.2
region box block 0 15 0 15 0 15
create_box 1 box
create_atoms 1 box

### DEFINE PHYSICAL PROPERTIES OF ATOMS ###
mass 1 1.0
pair_style lj/cut/opt 3.0
pair_coeff 1 1 1.0 1.0
neighbor 2.0 bin

### SPECIFY THE REQUIRED THERMODYNAMIC STATE ###
variable T equal 2.8
variable timestep equal 0.0025

### ASSIGN ATOMIC VELOCITIES ###
velocity all create ${T} 12345 dist gaussian rot yes mom yes

### SPECIFY ENSEMBLE ###
timestep ${timestep}
fix nve all nve

### THERMODYNAMIC OUTPUT CONTROL ###
thermo_style custom time etotal temp press
thermo 10

### RECORD TRAJECTORY ###
dump traj all custom 1000 output-1 id x y z

### SPECIFY TIMESTEP ###

### RUN SIMULATION TO MELT CRYSTAL ###
run 10000
unfix nve
reset_timestep 0

### BRING SYSTEM TO REQUIRED STATE ###
variable tdamp equal ${timestep}*100
fix nvt all nvt temp ${T} ${T} ${tdamp}
run 10000
reset_timestep 0
unfix nvt
fix nve all nve

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp press density atoms vol
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
variable N equal atoms
variable N2 equal atoms*atoms
variable etotal equal etotal
variable etotal2 equal etotal*etotal
variable volume equal vol

fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_etotal v_dens2 v_temp2 v_press2 v_etotal2
run 100000

variable avedens equal f_aves[1]
variable avetemp equal f_aves[2]
variable avepress equal f_aves[3]
variable aveetotal equal f_aves[4]
variable aveetotal2 equal f_aves[8]

variable Cv equal ${N2}*(${aveetotal2}-${aveetotal}*${aveetotal})/(${avetemp}*${avetemp})

print "Averages"
print "--------"
print "Density: ${avedens}"
print "Temperature: ${avetemp}"
print "Pressure: ${avepress}" 
print "Heat Capacity: ${Cv}"
print "Total Energy: ${aveetotal}"
print "Total Energy2: ${aveetotal2}"
print "Volume: ${volume}"
print "No of atoms square: ${N2}"

</pre>

'''<span style="color:red">JC: Heat capacity correctly calculated in script.</span>'''

The temperature and heat capacity at constant volume are generated and shown in Table 2.

'''Table 2'''

{| {{table}}
| align="center" style="background:#f0f0f0;"|'''Density'''
| align="center" style="background:#f0f0f0;"|'''T(T*)'''
| align="center" style="background:#f0f0f0;"|'''Cv/V'''
|-
| 0.2||1.97||1.47E-05
|-
| 0.2||2.17||1.05E-05
|-
| 0.2||2.39||7.49E-06
|-
| 0.2||2.65||6.00E-06
|-
| 0.2||2.76||5.25E-06
|-
| 0.8||1.99||2.56E-05
|-
| 0.8||2.17||1.99E-05
|-
| 0.8||2.36||1.96E-05
|-
| 0.8||2.65||1.45E-05
|-
| 0.8||2.79||1.36E-05
|-
| 
|-
| 
|}

[[File:Figure 9 - A plot showing the heat capacity as a function of temperature at 2 different densities.png|thumb|center|500 px|Figure 9 - A plot showing the heat capacity as a function of temperature at 2 different densities]]

Heat capacity is related to a fluctuation in the canonical ensemble but not ensemble average. It is a physical property measuring how easy it is to excite the system to higher energy state. From a thermodynamics point of view, it is the amount of heat needed to raise the system temperature by one degree in Celsius or Kelvin. As we can see from Figure 9, as temperature increases, heat capacity decreases for both density. The drop is more significant when temperature is raised from 2.2 to 2.4 for <math> \rho </math> = 0.2 and from 2.4 to 2.6 for <math> \rho </math> = 0.8. This is an expected trend for a liquid. As temperature increases, atoms have more kinetic energy and more atoms are populated to higher energy state. The energy levels at higher state have a lower separation and get closer together. There is a higher density of states at higher energy level. Less energy is needed to populate another atom at higher energy level. Therefore the heat capacity at constant volume drops. Also, as mentioned in the paper by Trachenko [2], the decrease of liquid heat capacity as temperature increases is attributed to the increasing loss of two transverse modes with frequency <math> \omega < \frac{1}{\tau}</math>, <math> \tau </math> is the liquid relaxation time.

'''<span style="color:red">JC: Nice, clear explanation of the effect of temperature on heat capacity.  What about the change in heat capacity with density?</span>'''

== Structural properties and the radial distribution function ==

Radial distribution function is given in <math> (6) </math>

<center><math> g (r) = \frac {\rho (r)}{\rho} \ \ (6) </math></center>

<math> \rho (r) </math> is the local density; <math> \rho </math> is the overall density of atoms
   
In this section, the radial distribution function (RDF) of solid, liquid and vapour phases of the Lennard-Jones fluid are calculated using VMD. 

'''<big>TASK 19</big>: 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>\int g(r)\mathrm{d}r</math>. Plot the RDFs for the three systems on the same axes, and attach a copy to your report. Discuss qualitatively the differences between the three RDFs, and what this tells you about the structure of the system in each phase. In the solid case, illustrate which lattice sites the first three peaks correspond to. What is the lattice spacing? What is the coordination number for each of the first three peaks?'''

Use the literature as a reference to modify the density and temperature parameters in the input script for the RDF. (Table 3) The files are uploaded to the HPC portal and allow it to run. Download the trajectory file after the simulation finished and load it using VMD. The RDF of solid, liquid and vapour are plotted in Figure 10 and the integrals are plotted in Figure 11.

'''Table 3'''

{| {{table}}
| align="center" style="background:#f0f0f0;"|'''Phases'''
| align="center" style="background:#f0f0f0;"|'''Temperature'''
| align="center" style="background:#f0f0f0;"|'''Density'''
|-
| Solid||1||1.2
|-
| Liquid||1.2||0.8
|-
| Vapour||1.2||0.1
|-
| 
|}

[[File:Figure 10 - Radial distribution functions of solid, liquid and vapour.png|thumb|center|600 px|'''Figure 10 - Radial distribution functions of solid, liquid and vapour''']]

[[File:Figure 11 - A plot showing the number of integral over g (r) as a function of distance for solid, liquid and vapour.png|thumb|center|600 px|'''Figure 11 - A plot showing the number of integral over g (r) as a function of distance for solid, liquid and vapour''']]

RDF gives the probability of finding an atom in distance r from another atom (reference atom). It indicates how atoms are radially arranged among each other.  A peak shows a favourable distance of separation for the reference atom away from its neighbour. Three common features of RDFs in solid, liquid and vapour are observed. One is all g (r) is equal to 0 below r = 0.825 due to the strong repulsive forces experienced between two positively charged nuclei and the electron clouds. Secondly, g(r) go through a maximum at around r = 1 before they plunge, indicating the chance of finding a pair of atoms at this separation is maximized at r = 1 and decreases afterwards. Thirdly, at longer distances, the g(r) approach one because the overall density is equal to the local density. 

Solid has a ordered structure with regular spacing so neighboring atoms can be found at a defined distances as shown by the peaks in Figure 10. Atoms are closely packed so there is a higher chance to find nearest neighbor so there are more peaks observed in solid than liquid and gas. 

Unlike solid, liquid has an irregular spacing between atoms and does not have long range order but short range order only. The RDF of liquid appears to be the intermediate form between the RDF of solid and vapour. 4 peaks are seen in the RDF at 4r distances, indicating there are 4 favourable separations for the reference atom. Atoms are packed less closely and there is a lower chance of finding the nearest neighbour than solid. The maximum of g(r) is half of that in solid. For gas, it is a smooth RDF with only 1 peak is seen at r = 1.125, showing gas has the least ordered structure with even more irregular spacing between atoms. The short range order disappears completely. The g(r) is similar to that of liquid.

'''<span style="color:red">JC: Graphs look nice and good description of the RDFs and what they show.</span>'''

[[File:Figure 12 - Illustration of the first 3 nearest neighbour on a FCC lattice.png|thumb|right|500 px|'''Figure 12 - Illustration of the first 3 nearest neighbours on a FCC lattice''']]

<math> a </math> is the length of the fcc lattice,

The r values for the first three peaks for solid are 1.03, 1.53 and 1.83. They correspond to the 1st, 2nd and 3rd nearest atoms respectively. (Indicated as 1, 2 and 3 in Figure 12)

<math> r = 1.03 </math>, the distance of the 1st nearest atom is <math> a (\frac {1}{2})^{1/2} = 1.03 </math>

<math> a = 1.45 </math>

<math> r = 1.53 </math>, the distance of the 2nd nearest atom is <math> a = 1.53 </math>

<math> a = 1.53 </math>

<math> r = 1.83 </math>, the distance of the 3rd nearest atom is <math> a (\frac {3}{2})^{1/2} = 1.83 </math>

<math> a = 1.49 </math>

Average lattice spacing, <math> a_{avg} </math> = <math> \frac {1.45 + 1.53 + 1.49}{3} = 1.5 </math>.

Figure 13 shows a snapshot of the RDF of the solid with r ranging from 0 to 2. Using this figure, we can find the coordination number for each of the first 3 peaks. 

[[File:Figure illustrating the coordination number = 12 for the first nearest neighbour atom in FCC.png|thumb|left|'''Figure 14 - Figure illustrating the coordination number = 12 for the first nearest neighbour atom in FCC [6]''']]

[[File:Fcc first 12 neighbours.png|thumb|left|'''Figure 15 - Figure illustrating the position of the first 12 neighbours in FCC [5]''']]


[[File:Figure illustrating the coordination number = 6 for 2nd nearest neighbour atom in FCC.png|thumb|right|'''Figure 16 - Figure illustrating the coordination number = 6 for 2nd nearest neighbour atom in FCC [6]''']]

[[File:Figure illustrating the coordination number = 24 for the third nearest neighbour atom in FCC.png|thumb|right|'''Figure 17 - Figure illustrating the coordination number = 24 for the third nearest neighbour atom in FCC [6]''']]

Coordination number of the 2nd peak = 12 (Figure 14-15)

As we know, fcc is a cubic closed packed structure. There are 6 nearest neighbours in one layer, 3 atoms in the above and 3 in the bottom. 

Coordination number of the 2nd peak = 18-12 = 6 (Figure 16)

Coordination number of the 3rd peak = 42 - 18 = 24 (Figure 17)

[[File:Figure 13 - Figure showing the integral of g (r) of solid from r=0 to r = 2.png||500 px|center|thumb|'''Figure 13 - Figure showing the integral of g (r) of solid from r=0 to r = 2''']]

'''<span style="color:red">JC: Good use of pictures to illustrate which atoms contribute to each coordination number.</span>'''



== Dynamical properties and the diffusion coefficient  ==

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

'''<big>TASK 20</big>: In the D subfolder, there is a file ''liq.in'' that will run a simulation at specified density and temperature to calculate the mean squared displacement and velocity autocorrelation function of your system. Run one of these simulations for a vapour, liquid, and solid. You have also been given some simulated data from much larger systems (approximately one million atoms). You will need these files later.'''

Use the density and temperature determined in the early section to modify the input file.

'''<big>TASK 21</big>: make a plot for each of your simulations (solid, liquid, and gas), showing the mean squared displacement (the "total" MSD) as a function of timestep. Are these as you would expect? Estimate <math>D</math> in each case. Be careful with the units! Repeat this procedure for the MSD data that you were given from the one million atom simulations. '''

The MSD as a function of timestep for solid, liquid and vapour and the simulated results for one million atoms for each phases are plotted in Figure 18 - 20 respectively. To find the diffusion coefficient, as we know from <math> (7) </math>, the diffusion coefficient can be found from the plot of total MSD as a function of time. To convert timestep into time, I multiplied each by 0.002 for all the simulations. The MSD as a function of time for solid, liquid and vapour are plotted in Figure 21 - 23 respectively. 

[[File:Figure 14 - MSD for solid as a function of timestep.png|thumb|right|500 px|'''Figure 18 - MSD for solid as a function of timestep''']]

[[File:Figure 15 - MSD for liquid as a function of timestep.png|thumb|right|500 px|'''Figure 19 - MSD for liquid as a function of timestep''']]

[[File:Figure 16 - MSD for vapour as a function of timestep.png|thumb|right|500 px|'''Figure 20 - MSD for vapour as a function of timestep''']]

[[File:Figure 17 - MSD for solid as a function of time.png|thumb|right|500 px|'''Figure 21 - MSD for solid as a function of time''']]

[[File:Figure 18 - MSD for liquid as a function of time.png|thumb|right|500 px|'''Figure 22 - MSD for liquid as a function of time''']]

[[File:Figure 19 - MSD for vapour as a function of time.png|thumb|right|500 px|'''Figure 23 - MSD for vapour as a function of time]]

Diffusion coefficient, <math> D = \frac {1}{6} \times gradient </math>

Best fit line equations are fitted for all cases. For the solid and vapour case, in order to get a more accurate gradient, the gradient given in the Excel are not taken. I have tried to fit the equation within the linear range for both phases and it in turn gives a significant fluctuation and a best fit line in Excel is not suitable for both cases. Therefore, the gradient of the linear part of the MSD for solid and vapour are measured manually.

'''<span style="color:red">JC: This is correct, only the linear part of the MSD should be used for the fit.</span>'''

Results are shown in Table 4.

'''Table 4'''

{| {{table}}
| align="center" style="background:#f0f0f0;"|'''Phases'''
| align="center" style="background:#f0f0f0;"|'''Gradient of simulated result'''
| align="center" style="background:#f0f0f0;"|'''D (σ^{2}/unit time)'''
| align="center" style="background:#f0f0f0;"|'''Gradient of one million atoms'''
| align="center" style="background:#f0f0f0;"|'''D  (σ^{2}/unit time)'''
|-
| Solid||3.00E-05||5.00E-06||3.00E-05||5.00E-06
|-
| Liquid||0.509||8.48E-02||0.524||8.73E-02
|-
| Vapour||9.10||1.52E+00||15.2||2.53E+00
|-
| 
|}

In MSD, we model the displacement of an atom in a finite volume. As we can see from the graphs, the MSD for solid, liquid and vapour differs significantly within x = 0 to 500. This can be explained by the molecular dynamics through collision between atoms and the force of restoration on each atom to remain in its position in a lattice depending on the phases. For solid phase, the MSD attains a maximum at t = 0.152 and minimum at t = 0.318. Minimum is the point where atoms collide with another atom and lead to a displacement in an opposite direction to its original direction of collision. In MSD model, we cannot anticipate the direction of atoms after collision as the probability of getting both directions are the same. Atoms in solid are closely packed and they cannot diffuse very far from its lattice site. Unlike liquid, the atomic motions in solid are vibratory. There is much less collisions occur in solid compared to liquid and therefore the displacement of solid levels off quickly at around 0.002. This gives solid the lowest number of D.

'''<span style="color:red">JC: Particles in solids are fixed to their lattice points so MSD should be horizontal.</span>'''

[[File:Figure - A graph showing the linear part of the MSD for solid.png|thumb|center|500 px|'''Figure 24 - A graph showing the linear part of the MSD for solid''']]

Liquid is less closely packed and there are more space between atoms. The atomic motions are Brownian motion like, where atoms can continuously collide with the neighbouring atoms and diffusive behaviour is observed in liquid. As we can see, MSD increases linearly with time. As time proceeds, there is a greater displacement of atoms. This gives a greater D than solid.

For vapour, the atoms are not closely packed and more randomized with a lot of spacing between atoms. After a collision, there is a significant change in displacement. At non zero pressure, the MSD will increase drastically but finitely as time proceeds. There is the greatest diffusive behaviour observed in gas phase and thus gives the highest D among all three phases.

Veolocity autocorrection function (VACF) is shown in <math> (8) </math>:

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

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

'''<big>TASK 22</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!) Be sure to show your working in your writeup. On the same graph, with x range 0 to 500, plot <math>C\left(\tau\right)</math> with <math>\omega = 1/2\pi</math> and the VACFs from your liquid and solid simulations. What do the minima in the VACFs for the liquid and solid system represent? Discuss the origin of the differences between the liquid and solid VACFs. The harmonic oscillator VACF is very different to the Lennard Jones solid and liquid. Why is this? Attach a copy of your plot to your writeup.'''

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


Position of a classical harmonic oscillator is given by <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>


Differentiate it and we can get the velocity of the atom.


<math> v = \frac{\mathrm{d}\mathbf{x\left(t\right)}}{\mathrm{d}t} = -A\omega\sin(\omega t+\phi) </math>


Substitute the result obtained into the VACF,


<math>C(\tau) = \frac{\int_{-\infty}^{\infty}(-A\omega\sin(\omega t+\phi))(-A\omega\sin(\omega(t+\tau)+\phi))\mathrm{d}t}{\int_{-\infty}^{\infty}(-A\omega\sin(\omega t+ \phi))^2\mathrm{d}t}

</math>

<math> A^{2} </math> and <math> \omega^{2} </math> are cancelled out


<math>
C(\tau) = \frac{\int_{-\infty}^{\infty}\sin(\omega t+\phi)\sin(\omega(t+\tau)+\phi)\mathrm{d}t}{\int_{-\infty}^{\infty}\sin^2(\omega t+\phi) \mathrm{d}t}
</math>


<math>
C(\tau) = \frac{\int_{-\infty}^{\infty}\sin(\omega t+\phi)\sin((\omega t+\omega \tau)+\phi)\mathrm{d}t}{\int_{-\infty}^{\infty}\sin^2(\omega t+\phi) \mathrm{d}t}

</math>


Using trigonometry identity, we know <math>\sin(A+B)=\sin(A)\cos(B)+\cos(A)\sin(B)</math>


<math>
C(\tau) = \frac{\int_{-\infty}^{\infty}\sin(\omega t+\phi)[\sin(\omega t+\phi)cos(\omega \tau)+\cos(\omega t+\phi)\sin(\omega \tau)]\mathrm{d}t}{\int_{-\infty}^{\infty}\sin^2(\omega t+\phi)\mathrm{d}t}
</math>


<math>
C(\tau) = \frac{\int_{-\infty}^{\infty}\sin^2(\omega t+\phi)\cos(\omega \tau)\mathrm{d}t+\int_{-\infty}^{\infty}\sin(\omega t+\phi)cos(\omega t+\phi)\sin(w\tau) \mathrm{d}t}{\int_{-\infty}^{\infty}\sin^2(\omega t+\phi)\mathrm{d}t}
</math>


<math>
C(\tau)=\frac{\cos(\omega \tau)\int_{-\infty}^{\infty}\sin^2(\omega t+\phi)\mathrm{d}t}{\int_{-\infty}^{\infty}\sin^2(\omega t+\phi)\mathrm{d}t}+\frac{\sin(\omega \tau)\int_{-\infty}^{\infty}\sin(\omega t+\phi)\cos(\omega t+\phi)\mathrm{d}t}{\int_{-\infty}^{\infty}\sin^2(\omega t+\phi)\mathrm{d}t}
</math>


sin is an odd function and cosine is an even function. When an odd function is multiplied by an even function, the result is an odd function. 


The integral of an odd function with a symmetrical limit becomes 0, so 


<math>
\frac{\sin(\omega \tau)\int_{-\infty}^{\infty}\sin(\omega t+\phi)\cos(\omega t+\phi)\mathrm{d}t}{\int_{-\infty}^{\infty}\sin^2(\omega t+\theta)\mathrm{d}t} = 0
</math>


<math>C(\tau)=\cos(\omega \tau)</math>

'''<span style="color:red">JC: Derivation is correct and justifying each step makes it very clear to follow.</span>'''


The VACF of solid and liquid and the normalized VACF of 1D harmonic oscillator are plotted in Figure 25. 

[[File:Figure 20 - VACF for solid and liquid and 1D harmonic oscillator.png|thumb|right|500 px|'''Figure 25 - VACF for solid and liquid and 1D harmonic oscillator''']]

For 1D harmonic oscillator, there is a specific velocity at an instant time. There is no collision with other atoms and hence the velocity of the atoms will not change unless it hits the boundary and causes a change in velocity and direction. It results a periodic function as we see in Figure 25. However, atoms can collide with each other in solid, liquid and vapour. Any collision will change the velocity of an atoms. The velocity of the atom before and after collision are different and we say this is decorrelated. As we can see, the VACF decorrelates with time for solid and liquid but in different extent. In solid, there is a fixed regular pattern of atoms arrangement and they cannot move significantly away from their lattice site and collide with another atom. More oscillations are seen in the VACF for solid, indicating there are positive and negative correlations for the velocity as time proceeds before reaching 0. There is a reverse velocity after each oscillation. We can see the magnitude of correlation is different for each turn for solid because there are perturbative forces on atoms in solid to alter the oscillating behaviour. Unlike solid, liquid does not have fixed pattern. We can only see one minimum in the VACF for liquid before it drops to 0. Atoms in liquid can diffuse more and collide with neighbour atoms, which removes the oscillating behaviour we see in solid. The velocity decays rapidly with time. This also explains the value of D we obtain.

'''<big>TASK 23</big>: Use the trapezium rule to approximate the integral under the velocity autocorrelation function for the solid, liquid, and gas, and use these values to estimate <math>D</math> in each case. You should make a plot of the running integral in each case. Are they as you expect? Repeat this procedure for the VACF data that you were given from the one million atom simulations. What do you think is the largest source of error in your estimates of D from the VACF?'''

Trapezium rule is given as <math> (10) </math>

<math> y = f (x) </math>

[[File:Eq 10 - Trapezium rule.png|thumb|right|'''Eq 10 - Trapezium rule''']]

Using the trapezium rule, the integral for solid, liquid and vapour are obtained and the diffusion coefficients are found by multiplying the integral by <math> frac {1}{3} </math>.

The corresponding diffusion coefficients are calculated (Table 5)

'''Table 5'''
{| {{table}}
| align="center" style="background:#f0f0f0;"|'''Phases'''
| align="center" style="background:#f0f0f0;"|'''Integral'''
| align="center" style="background:#f0f0f0;"|'''D'''
| align="center" style="background:#f0f0f0;"|'''Integral'''
| align="center" style="background:#f0f0f0;"|'''D'''
|-
| Solid||5.64E-04||1.88E-04||1.37E-04||4.57E-05
|-
| Liquid||0.293||0.0977||0.27||0.0901
|-
| Vapour||5.23||1.74||9.81||3.27
|-
| 
|}

Table 6 shows a summary of D obtained from VACF and MSD.

'''Table 6'''
{| {{table}}
| align="center" style="background:#f0f0f0;"|'''Phases'''
| align="center" style="background:#f0f0f0;"|'''D obtained from MSD'''
| align="center" style="background:#f0f0f0;"|'''D obtained from MSD (1 million atoms)'''
| align="center" style="background:#f0f0f0;"|'''D obtained from VACF'''
| align="center" style="background:#f0f0f0;"|'''D obtained from VACF (1 million atoms)'''
|-
| Solid||6.05E-06||2.58E-06||1.88E-04||4.57E-05
|-
| Liquid||0.0848 ||0.0873 ||0.0977||0.0901
|-
| Vapour||1.49 ||2.68 ||1.74||3.27
|-
| 
|}

The integral and D found are what we expect for each phases. Both VACF and MSD model give the same trend, where vapour gives the largest D, followed by liquid and then solid. This match with the characteristics of each phases. Running integrals match the VACF we get in the previous task. For example, liquid has a single peak and then plateau while solid matches the fluctuations seen in the VACF .  The 1 million atom simulation systems give the same trend but different values for both model. As the simulation system size increases, less errors are encountered in the simulation as more atoms are taken into account. In MSD, greater fluctuations are seen in solid for a small system than a large system. The MSD does not effectively level off and becomes constant as time proceeds. Also, the D for vapour calculated from MSD is half of that of one million system for both model. 

To calculate D in VACF, we just use trapezium rule is only an approximation of the area under the curve for an integral. The integrals are not exact so this incurs a significant error in VACF. To enhance the accuracy, a smaller timestep can be used. The largest source of error in VACF is we only modify the temperature and density parameters to model a particular phase but not specify the positions of the atoms. Velocity is the only factor we take into account in VACF, which can be significantly affected by the position of the atoms. Position of atoms can vary significantly for each run and the velocity is also affected. For all phases and simulating systems, VACF has larger value than that of MSD. The trapezium rule may overestimate the area under the curve of the integral. Liquid gives the closest overall match between VACF and MSD. 

'''<span style="color:red">JC: Graphs look good, refer to them directly in the text to help illustrate your points.  The position of atoms isn't important as we average over all atoms when we calculate the VACF.</span>'''

[[File:Solid running.png|thumb||'''Figure 26 - Figure showing the running integral for solid''']]

[[File:Solid running one million.png|thumb|'''Figure 27 - Figure showing the running integral for solid with one million atoms''']]

[[File:Liquid running.png|thumb|'''Figure 28 - Figure showing the running integral for liquid''']]

[[File:Liquid one million atoms.png|thumb|'''Figure 29 - Figure showing the running integral for liquid with one million atoms]]

[[File:Running integral vapour.png|thumb|'''Figure 30 - Figure showing the running integral for vapour''']]

[[File:Running integral vapour one millon atoms.png|thumb|'''Figure 31 - Figure showing the running integral for vapour with one millon atoms''']]

===Reference===

[1] JI.Choe and B.Kim, Bull. Korean Chem.Soc. 2000, Vol.21(4) 419 - 424

[2] K. Trachenko, Phys. Rev B, 2008, 78, 104201

[3] D. Fincham, Computer Physics Communications 1986, 40, 263-269

[4] M. Winger, D. Trzesniak, R. Baron, W.F. van Gunsteren, Phys. Chem. Chem. Phys., 2009, 11, 1934-1941

[5] A Brief Note: Why Fcc Has 12 Nearest Neighbours, http://ph.qmul.ac.uk/~anthony/spfm/fcc.html,18 Dec. 2015

'''<span style="color:red">JC: Good use of literature.</span>'''

[6] FCC And Its Neighbours, http://www.ifmpan.poznan.pl/~urbaniak/fcc%20and%20its%20neigbours01.pdf, 18 Dec. 2015.