View source for Talk:Y3C Liquid Simulation:hjt14

'''<span style="color:red">JC: General comments: Very good report, results are well presented and accompanied by good explanations which demonstrate a good understanding of the background theory behind this experiment, well done.</span>'''

==Section 1: Running First Simulation==
No task.

==Section 2: Introduction to Molecular Dynamics Simulation==
===Numerical Integration===
====Task====
[[Media:HO hjt14.xls| Link]] to the completed HO.xls file. The followings are the graph constructed using data from the excel worksheet. The first graph is an overlap of two graph of position of simple harmonic oscillator, <math> x(t) </math>, as a function of time, <math> t</math>; one calculated using velocity verlet algorithm, and one calculated using <math>  x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>. The second graph is the graph of total energy of the simple harmonic oscillator, <math> E_{total}(t) </math>, as a function of time, <math> t</math>. The third graph is the graph of absolute error (difference between in position, <math> x(t) </math>, calculated using classical harmonic oscillator solution and calculated using velocity verlet algorithm) as a function of time, <math> t</math>. The maximas were located from each parabola and were fitted with a linear function. The equation of the linear function is displayed on the graph.

{| class="wikitable" border="1"
! Graphs from excel file
|-
| [[File:Graph of Position of Simple Harmonic Oscillator Against Time (Time Step= 0.100).PNG|x500px|thumb|Graph 1]] 
|-
| [[File:Graph of Total Energy of Simple Harmonic Oscillator Against Time (Time Step= 0.100).PNG|x500px|thumb|Graph 2]]
|-
| [[File:Graph of Absolute Error Against Time (Time Step= 0.100).PNG|x500px|thumb|Graph 3]]
|}

'''<span style="color:red">JC: Why does the error oscillate?</span>'''

Since we are working in a close system, i.e., no energy exchange with surrounding, the total energy of the system should be a constant. What we observed from graph 2 is that the total energy of the system oscillates about an average value. The deviation from the average energy is treated as the error introduced by the approximation we made in velocity verlet algorithm. The percentage error of total energy at different timestep were computed using following equation and the results are summarized in graph 4:

Percentage Error <math> = \frac{(E_{max}-E_{min})}{2\times E_{average}} \times100 % </math>

[[File:Graph of Percentage Error of Total Energy Against Timestep.PNG|x500px|thumb|centre|Graph 4]]

It can be seen from graph 4 that for timestep beyond 0.28, the percentage error in total energy will exceed 1%.

'''<span style="color:red">JC: Good, thorough analysis and choice of timestep.</span>'''

===Atomic Forces===
====Task====
Lennard-Jones Equation: <math>\phi (r)=4\epsilon \left( \frac{\sigma ^{12}}{r^{12}}-\frac{\sigma ^6}{r^6} \right) </math>

When potential energy is set to zero, i.e.<math>\phi =0 </math>:

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

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

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

<math>\sigma  = r = r_0</math>

The value of <math>r</math> when potential energy <math>\phi</math> is set to zero,<math> r_0 </math>, is equal to <math>\sigma </math>.

The force acting on a particle a distance <math> r </math> away from another particle can be calculated using following equation:

<math> F= - \frac {d\phi(r)}{dr} </math>

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

When <math> r= r_0 =\sigma </math>:

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

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

<math>F=\frac{24\epsilon}{\sigma} </math>

At distance <math> r_0 </math>, where the potential energy, <math>\phi</math>, is zero, the force acting on a particle is equal to <math>\frac{24\epsilon}{\sigma} </math>. 

When <math>F=0</math>, defines <math> r=r_{eq}</math>:

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

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

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

<math> r_{eq}= \sqrt[6]{2} \sigma </math>


The equilibrium distance,<math> r_{eq} </math> , between 2 particles, which is defined as separation of 2 particles when the force acting on each particle is zero, is equal to <math>\sqrt[6]{2} \sigma </math>, <math> \sigma </math> is a constant and is specific to particle under study. At equalibrium distance, <math> r_{eq} </math>, the potential energy, <math> \phi(r_{eq}) </math>, is at minimum and can be caculated as follow:

<math> \phi(r_{eq}) = 4 \epsilon \left(\frac{\sigma^{12}}{(\sqrt[6]{2}\sigma)^{12}}-\frac{\sigma^6}{(\sqrt[6]{2}\sigma)^6} \right) </math>

<math> \phi(r_{eq}) = 4 \epsilon \left(\frac{1}{4}-\frac{1}{2} \right) </math>

<math> \phi(r_{eq}) = -\epsilon </math>


The following part shows the working of the integration of potential energy, <math> \phi </math>, with respect to <math> r </math>. The integral is equivalent to the area under the potential energy graph encompassed by the integration range, <math> x </math> to <math> y</math>.

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

<math> \int_x^y \phi(r) dr = 4\epsilon \left[\frac{-\sigma^{12}}{11r^{11}}+\frac{\sigma^6}{5r^5} \right]_x^y </math>

When <math> \epsilon =1</math> and <math> \sigma=1 </math>:

<math> \int_x^y \phi(r) dr = 4 \left[\frac{-1^{12}}{11r^{11}}+\frac{1^6}{5r^5} \right]_x^y </math>


{| class="wikitable" border="1"
|When <math> x=2\sigma=2 </math> and <math> y=\infty </math>: 
|When <math> x=2.5\sigma=2.5 </math> and <math> y=\infty </math>:
|When <math> x=3\sigma=3 </math> and <math> y=\infty </math>:
|-
|<math> \int_2^\infty \phi(r) dr = 4 \left[\frac{-1^{12}}{11r^{11}}+\frac{1^6}{5r^5} \right]_2^\infty </math>

<math> \int_2^\infty \phi(r) dr = 4 \left[0-\left(\frac{-1^{12}}{11(2^{11})}+\frac{1^6}{5(2^5)} \right)\right] </math>

<math> \int_2^\infty \phi(r) dr = -\frac{699}{28160} = -0.0248 </math> 
|<math> \int_{2.5}^\infty \phi(r) dr = 4 \left[\frac{-1^{12}}{11r^{11}}+\frac{1^6}{5r^5} \right]_{2.5}^\infty </math>

<math> \int_{2.5}^\infty \phi(r) dr = 4 \left[0-\left(\frac{-1^{12}}{11(2.5^{11})}+\frac{1^6}{5(2.5^5)} \right)\right] </math>

<math> \int_{2.5}^\infty \phi(r) dr = -0.00818 </math>
|<math> \int_3^\infty \phi(r) dr = 4 \left[\frac{-1^{12}}{11r^{11}}+\frac{1^6}{5r^5} \right]_3^\infty </math>

<math> \int_3^\infty \phi(r) dr = 4 \left[0-\left(\frac{-1^{12}}{11(3^{11})}+\frac{1^6}{5(3^5)} \right)\right] </math>

<math> \int_3^\infty \phi(r) dr = -0.00329 </math>
|}

'''<span style="color:red">JC: All maths correct and well laid out.</span>'''

===Periodic Boundary Conditions===
====Task====
The number of water molecules in 1 ml of water under standard conditions (298K and 1 atm) is calculated as follow:

Density of water under standard conditions<math> = 0.9970 g ml^{-1} </math>

Mass of 1 ml of water <math>=0.9970 g </math>

Molar mass of water <math> = 2(1.008)+15.999 = 18.015 gmol^{-1} </math>

Number of water molecule in 1 ml of water <math> = \frac{0.9970(6.022\times 10^{23})}{18.015} = 3.333\times10^{22} </math>

The volume occupied by 1000 water molecules under standard conditions is calculated as follow:

Mass of 1000 water molecules <math>=\frac{1000}{6.022\times10^{23}} (18.015) = 2.992\times 10^{-20}g  </math>

Volume occupied by 1000 water molecules <math> = \frac{2.992\times 10^{-20}}{0.9970}=3.001\times 10^{-20}ml </math>

====Task====
Under periodic boundary condition, the final position of an atom, after it move from 
<math>
\begin{pmatrix}
0.5\\
0.5\\
0.5
\end{pmatrix} 
</math>
by 
<math>
\begin{pmatrix}
0.7\\
0.6\\
0.2
\end{pmatrix} 
</math>
in a cubic simulation box which runs from <math> \left(0, 0, 0\right)</math> to <math>\left(1, 1, 1\right) </math>, can be calculated as follow:

New position of atom <math> =
\begin{pmatrix}
0.5\\
0.5\\
0.5
\end{pmatrix}
+
\begin{pmatrix}
0.7\\
0.6\\
0.2
\end{pmatrix}
=
\begin{pmatrix}
1.2\\
1.1\\
0.7
\end{pmatrix}
</math>

Apply periodic boundary condition:

New position of atom
<math>=
\begin{pmatrix}
1.2-1\\
1.1-1\\
0.2
\end{pmatrix}
=
\begin{pmatrix}
0.2\\
0.1\\
0.2
\end{pmatrix}
</math>

===Reduced Units===
====Task====
{| class="wikitable" border="1"
|The real units, <math> r </math>, is calculated as follow:
|The well depth in <math> kJ mol^{-1} </math> is calculated as follow:
|The reduced temperature <math> T^* = 1.5 </math> in real units is calculated as follow:
|-
|
<math> r^*= \frac{r}{\sigma} </math>

<math> r^*= 3.2 </math> and <math> \sigma= 0.34 nm </math>

<math> r= 3.2(0.34) nm </math>

<math> r= 1.088 nm </math>
|
<math> \frac{\epsilon}{N_Ak_B} =120 K </math>

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

<math> \epsilon = 120(8.31) J mol^{-1} </math>

<math> \epsilon = 997 J mol^{-1}= 0.997 kJ mol^{-1} </math>
|
<math> T^*=\frac{k_BT}{\epsilon} </math>

<math> 1.5= \frac{1}{120}T </math>

<math> T= 180 K </math>
|}

'''<span style="color:red">JC: All calculations correct and clearly laid out.</span>'''

==Section 3: Equilibration==
===Creating the Simulation Box===
====Task====
Within an input file for molecular simulation, the initial conditions stipulated may not be the equilibrium conditions. During a simulation, LAMMPS will try to bring the system down to the equilibrium state by iteration within the allowed timesteps specified in the input file. If random coordinate is assigned to each atom, there is a chance that two or more atoms will be assigned coordinates that are too close together. This will mean that the atoms overlap in space. According to Lennard-Jones equation, this condition will give rise to the extremely high potential energy (approach infinity) and large repulsion between atoms. As a result of this situation, calculations done using LAMMPS sumulator will run into errors and the system will not be brought to equilibrium state.

====Task====
A unit cell of a simple cubic lattice is a cube dotted with lattice point at each vertice. Each vertice is shared between 4 unit cells, hence each vertice will contribute a quarter of a lattice point to a unit cell. Consequently, each unit cell will have lattice point, 4 corners each contributes a quarter of a lattice point. The number density can be calculated by dividing the number of lattice point in a unit cell by the volume of a unit cell.

If we have a unit cell of size <math> 1.07722 \times 1.07722 \times 1.07722 </math> (in redued units):

Number density of a simple cubic lattice <math> = \frac{1}{1.07722 \times 1.07722 \times 1.07722} =0.8 </math>

If we have a face-centered cubic lattice with lattice point number density of <math>1.2</math>, we can calculate the size of unit cell as follow:

Number of lattice point in a face-centered cubic unit cell <math> =4</math>

Volume of unit cell <math>= l\times l\times l = l^3 </math>

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

<math> l=\sqrt[3]{3\frac{1}{3}} =1.49 </math>

The lattice parameter of a face-centered cubic lattice with lattice point number density of <math>1.2</math> is <math> 1.49 </math> (in reduced units).


====Task====
The number of unit cell created in the defined region is <math> 1000 </math>. In each face-centered cubic unit cell, <math> 4 </math> atoms will be generated. Consequently, the total number of atom that will be generated if a face-centered cubic lattice is used instead of a simple cubic lattice is <math>4000</math>.

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

===Setting the Properties of the Atoms===
====Task====
{| width=50%
| <pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>
|}

The first line of the command set the mass of the atom of type 1 to 1.0 (reduced units). The second line of the command defines the equation used to model the pairwise interaction energy between atoms. In this example, we are using Lennard-Jones Equation with a cutoff distance of 3.0 (reduced units). The third line of the command defines the coefficients of Lennard-Jones Equation, i.e <math> \epsilon =\sigma = 1.0 </math>. The asterisks (* *) mean that the Lennard-Jones model will be imposed on every pair of atoms.

'''<span style="color:red">JC: Good, why do we use a cutoff with a Lennard-Jones potential?</span>'''

====Task====
Given that we are specifying <math>x_i\left(0\right) </math> and <math>v_i\left(0\right) </math>, the integration algorithm that we will be using is velocity verlet  algorithm.

===Running the Simulation===
====Task====
<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>

By assigning a variable (timestep) to timestep (the value), it allows us to call 'timestep' at other part of the input script, meaning that we do not have to retype the exact value of timestep repeatly and allow us to call it simply by typing ${timestep}. Another advantage of assigning a variable (timestep) to timestep (the value), is that if we need to change the timestep of a simulation, we only need to change the value assigned to the variable 'timestep' and timestep in other part of the input file will be updated accordingly.

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

===Checking Equilibrium===
====Task====
{| class="wikitable" border="1"
|+ Graphs Obtained at 0.001 Timestep
|-
| [[File:Intro 0.001 Timestep total energy hjt14.png|thumb|500px|Graph 5]] 
| [[File:Intro 0.001 Timestep Temperature hjt14.png|thumb|500px|Graph 6]] 
| [[File:Intro 0.001 Timestep Pressure hjt14.png|thumb|500px|Graph 7]] 
|}

Above show three graphs obtained from LAMMPS simulation ran at 0.001 timestep. It can be seen that the system did achieved equilibrium and the inset in each graph show that system achieve equilibrium in approximately 0.4 fs. LAMMPS simulations were also ran at other timestep and the results for total energy were summarised in graph 8.   

[[File:Intro Total Energy Diff Timestep hjt14.png|thumb|centre|500px|Graph 8]]

From graph 8, it can be seen that the most suitable timestep is 0.0025 as it is the maximum value which still allow the system under study to achieve equilibrium (a compromise between accuracy and time). 0.015 is a particularly bad choice because the system simulated at that timestep did not converge.

'''<span style="color:red">JC: Good choice of timestep, but why are 0.01 and 0.0075 no good - should the average total energy depend on the timestep?</span>'''

==Section 4: Running Simulations under Specific Conditions==
===Thermostats and Barostats===
====Task====
In order to bring the system to target temperature, <math> \mathfrak{T} </math>, we multiply the velocity of each particle with a constant factor, <math>\gamma </math>, which will help to modulate the velocity of each atoms. The factor,<math>\gamma </math>, can be derived as follow:

<math> E_K=\frac{1}{2} \sum_i m_i v_i^2 = \frac{3}{2}NK_BT </math>  

Let <math> T= \mathfrak{T} </math> and multiply velocity of each atom with a constant factor, <math>\gamma </math>

<math> E_K=\frac{1}{2} \sum_i m_i \gamma^2v_i^2 = \frac{3}{2}NK_B\mathfrak{T} </math>  

<math> \gamma^2 \frac{1}{2} \sum_i m_i v_i^2 = \frac{3}{2}NK_B\mathfrak{T} </math>

<math> \gamma^2 \frac{3}{2}NK_BT = \frac{3}{2}NK_B\mathfrak{T} </math>

<math> \gamma = \sqrt{\frac{\mathfrak{T}}{T}} </math>

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

===Examining the Input Script===
====Task====
<pre>sampled
fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_dens2 v_temp2 v_press2
</pre>

The numbers and its position in the series will affect how the averages are calculated. The roles that each number play are summarised as follow:

{| class="wikitable"
! Position
! Value
! Role
|-
| First
| 100
| Sample data every 100 timestep
|-
| Second
| 1000
| Perform sampling 1000 times
|-
|Third
|100000
|Average all the sampled data collected within 100000 timestep
|}

===Plotting the Equations of State===
====Task====
{| class="wikitable"
! colspan='2'|Conditions used
! rowspan='2'|Graph
|-
!Pressure
!Temperature
|-
| rowspan='5'| 2.0 
| 1.8
| rowspan='10'| [[File:Section 4 NPT graph hjt14.png|thumb|500px|Graph 9]]
|-
|2.0
|-
|2.2
|-
|2.4
|-
|2.6
|-
| rowspan='5'| 3.0
| 1.8
|-
|2.0
|-
|2.2
|-
|2.4
|-
|2.6
|}

Shown in the table, are 10 phase points used in LAMMPS molecular dynamic simulation to compute density under these conditions. The theoretical value of density at these phase points were also calculated using reduced ideal gas equation. The simulated density and calculated density are summarised in graph 9. The reduced ideal gas equation used in calculation is shown below:

<math>

\rho =\frac{P}{T^*};   \rho=\frac{N}{V}

</math>

In graph 9, it can be seen that:

a) at all temperature, the simulated density are lower than the theoretical density.

Explanation: In ideal gas, we assume there are no inter-molecular interactions between particles, i.e., no attraction and no repulsion. However, in our LAMMPS simulation, we model the inter-molecular interactions using Lennord-Jones equation and this breaks the fundamental assumption made in the derivation of ideal gas equation. According to Lennard-Jones equation, below a certain bond distance, a particles pair will experience repulsive forces which will keep them at a distance away from each other. In this case, the particles are further apart from each other than they would be if they behave ideally. Consequently, the number of particle per unit volume, which is the density, is lowered.   

b) the simulated density deviate considerably from ideal behavior at low temperature. At higher temperature, the magnitude of deviation decreases and the theoretical density and simulated density start to converge to a single value.

Explanation: At low temperature, the thermal energy, <math> k_BT </math>, is small relative to the interaction energy. Therefore the interaction energy can not be neglected and contribute significantly to the thermodynamic properties of our system. As temperature starts to increase, the thermal energy increases and at one point the magnitude of the thermal energy will become comparable to the interaction energy. At that point, the interactions between particles can be easily perturbed by thermal fluctuation and hence become irrelevant in determining the thermodynamic properties of our system. In other words, at higher temperature, interactions between particles become insignificant and system approaches ideal behavior. Therefore, simulated density and theoretical density converge to a single value.

'''<span style="color:red">JC: Good explanations.  Why are the ideal gas and simulation results closer at lower pressure?</span>'''

==Section 5: Calculating Heat Capacities Using Statistical Physics==
====Task====
[[Media:Nvt(0.2, 2.0, 0.0025).in|Link]] to one of the input file. 
{| class="wikitable"
! colspan='2'|Conditions used
! rowspan='2'|Graph
|-
!Density
!Temperature
|-
| rowspan='5'| 0.2 
| 2.0
| rowspan='10'| [[File:Section 5 CvV vs T.png|thumb|500px|Graph 10]]
|-
|2.2
|-
|2.4
|-
|2.6
|-
|2.8
|-
| rowspan='5'| 0.8
| 2.0
|-
|2.2
|-
|2.4
|-
|2.6
|-
|2.8
|}

Pressure is defined as <math> \frac{N}{V} </math>, at higher pressure there will be more particle per unit volume. In section 4, we encounter following equation which correlates the instantaneous temperature of a system to total kinetic energy:

<math>


E_K = \frac{3}{2} N k_B T = \frac{1}{2}\sum_i m_i v_i^2

</math>

Occupying the same volume, a system at higher pressure will have a greater number of particles. Therefore, for the same increment in temperature, the energy we need to provide to a system at higher density is greater than that of a system at lower pressure simply because there are more particles to which we have to provide the energy. Consequently, the heat capacity per unit volume of a system at higher pressure is greater than that of system at lower pressure simply because we need to provide more energy for the same increment in temperature.

According to equipartition theorem, each degree of freedom of an atom will contribute <math> \frac{1}{2} k_B </math> to the total heat capacity of the system. As we are modeling a constant number of atoms, we expect the heat capacity of the system to be a constant. For both graphs, there is no clear reason to why the heat capacity reduces with temperature. One possible reason is that the density of state decreases with temperature. Higher density of state means that the energy levels are closer together and are easy to populate the energy states and hence correlate to lower heat capacity; Lower density of state means that the energy gap between energy levels is larger and and are difficult to populate the energy states. For same increment in temperature, greater energy need to be provided, and hence this correlate to higher heat capacity. One way to verify this is to obtain a plot of density of state versus energy.

'''<span style="color:red">JC: Good explanation of the trend with pressure and good ideas to explain the trend with temperature - more analysis, beyond the scope of this experiment, would be needed to confirm this trend.</span>'''

==Section 6: Structural Properties and the Radial Distribution Function==
===Calculating <math>g(r) </math> in VMD===
====Task====
<math> g(r)</math> can be understood as the number of particle per unit volume of shell of size <math> dr </math> divided by the bulk density.

{|class="wikitable"
!Graphs
!Observed trends in graphs when going from solid to liquid and to gas
!Remark
|-
|rowspan='3'| [[File:Part 6 gr hjt14.png|thumb|600px|Graph 11]]

[[File:Part 6 int g hjt14.png|thumb|600px|Graph 12]]

|The peak intensity decreases and has greater width
|In solid, the atoms are closely packed and hence are in higher coordination environment and this translates into high <math> g(r) </math>. From solid to gas, the coordination environment drop as the atoms are more loosely packed in this case and  hence <math>g(r)</math> decreases from solid to gas. The width of the peak is related to the distribution of <math> r </math>of atoms in the first coordination shell, then second coordination shell and etc.. In solid, the peaks are more confined than that of liquid and gas. This is because solid has low degree of freedom in motion, i.e., no translation, hence <math> r </math> of atoms in each coordination shell has smaller distribution and this gives rise to relatively sharp peaks. From liquid to gas, the translational degree of freedom increases, and hence the position of atoms in coordination shell becomes less and less well defined. This uncertainty in <math>r </math>gives rise to wide peaks.  
|-
|The rate of decay of <math> g(r) </math> increases
|This is because solid has long range order in the spatial arrangement of atoms. In liquid and gas, the ordered structure is short range and over <math> r</math> the peaks flattened out as atoms distribution is getting more and more random. 
|-
|Coordination environment decreases
|'''REFERRING TO SOLID''' The first coordination shell of solid consisted of 12 particles (cuboctahedron) and corresponds to the first peak of solid <math>g(r)</math> plot. The second coordination shell of solid consists of 6 particles (octahedron) and corresponds to the second peak of solid <math>g(r)</math> plot. The third coordination shell of solid consists of 24 particles and corresponds to the third peak of solid <math>g(r)</math> plot. The position of atoms in these shells are shown below relative to a reference point (coloured black):

[[File:FCC neighbour.jpg|thumb|400px|centre|FCC lattice made using crystal maker. Black atom is the reference atom, yellow atoms are in the first coordination shell, red atoms are in the second coordination shell and grey atoms are in the third coordination shell.]]

The lattice spacing of FCC lattice of solid is the same as the radius of the second coordination shell and is equal to 1.475 (refer graph 11)
 
|}

'''<span style="color:red">JC: The main result from the RDF is that the solid has short and long range order, the liquid has only short range order and the gas is completely disordered.  Good diagram of an fcc lattice to show the positions on the nearest neighbours.  Could you calculate the lattice spacing from the first and third peaks as well using the geometry of an fcc lattice and then average over them?</span>'''

==Section 7: Dynamic Properties and the Diffusion Coefficient ==
===Mean Squared Displacement===
====Task====
{| class="wikitable"
|
[[File:Part 7 Diffusion hjt146.png|thumb|centre|500px|Graph 13]]
|
[[File:Part 7 Diffusion hjt14.png|thumb|centre|500px|Graph 14]]
|
{| class='wikitable'
!State
!Gradient (w.r.t. time)
!Diffusion Coefficient (w.r.t. time)
|-
|Gas (million)
|18.19 (linear part)
|3.032
|-
|Gas
|4.915
|0.8192
|-
|Liquid (million)
|0.5236
|0.08727
|-
|Liquid
|0.5094
|0.0849
|-
|Solid (million)
|0.00002635
|0.000004392
|-
|Solid
|0.00535
|0.000892
|}

|}

From graph 13, it can be seen that the diffusion coefficients increase from solid to gas and these are expected. In solid, the particles are closely packed, to diffuse, solid particles have to overcome large repulsion when sliding pass other particles and this contribute to the large activation energy for diffusive motion. Hence solid has negligible diffusion coefficient due to large energy barrier.  From liquid to gas, the packing decreases meaning that the particles are further apart, making it easier for particle to move pass each other. The is reflected in the increase in diffusion coefficient when going from liquid to gas.

'''<span style="color:red">JC: Plot lines of best fit on the graph as well, did you fit to the entire data range or just to the linear part?</span>'''

===Velocity Autocorrelation Function===
====Task====
<math>
x(t)= A \cos \left(\omega t+\phi \right) ;  v(t)=\frac{dx(t)}{dt}=-A \omega \sin\left(\omega t+ \phi \right)
</math>

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


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


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




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


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


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


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


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

The second term goes to zero and then we are left with:


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


'''<span style="color:red">JC: Why does the second integral go to zero?  An even function multiplied by an odd function is an odd function, and the integral of an odd function with limits equally spaced about zero is zero.</span>'''

====Task====
{|
|[[File:Section 7 VACF graph hjt14.png|thumb|500px|Graph 15]]
|[[File:Section 7 VACF and SHO.png|thumb|500px|Graph 16]]
|}


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

Equation shown above projects velocity at <math> \left(t+\tau\right) </math> onto velocity at <math> \left(t\right) </math> and what we obtain is <math> C\left(\tau\right) </math>, a number which quantifies the correlation between the two velocities. From gas to solid, the rate of decay of <math> C\left(\tau\right) </math> increases and this is related to the rate of collision. In gas and liquid, particles collide with one another, however, in liquid, the frequency of collision is higher. Rapid exchange of energy between liquid particle causes a liquid particle to lost the initial velocity component at a faster rate than that of gas. From liquid to solid, the situation changes slightly. In a solid, unlike gas and liquid, the particles only vibrate about the equilibrium position. Because the motion of solid particle is confined within a small space, it correlate to initial velocity to a longer period of time. The minima in solid and liquid VACF vs time graph (see Graph 15 and Graph 16) is related to back scattering of particle after collision, hence the <math> C(\tau) </math> is negative in sign (inverse of the initial course). 

In simple harmonic oscillator the particles are vibrating about an equilibrium position and there is no collision with another particle, therefore it retain the memory of the initial velocity throughout the simulation. The simple harmonic oscillator aims to illustrate that the velocity correlation in ideal system will continue for infinite amount of time. 

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

====Task====
{| class='wikitable'
| [[File:Section 7 Running Int hjt14.png|thumb|500px|Graph 17]]
|
{|class='wikitable'
|Phase
|Average value of plateau region of graph
|Diffusion Coefficient from VACF
|Diffusion Coefficient from MSD
|-
|Solid (million)
|0.000192943
|0.000064314
|0.000004392
|-
|Solid
|0.00345
|0.00115
|0.000892
|-
|Liquid (million)
|0.25621
|0.08540
|0.08727
|-
|Liquid
|0.2316
|0.0772
|0.0849
|-
|Gas (million)
|9.8054
|3.2685
|3.032
|-
|Gas
|2.44197
|0.81399
|0.8192
|}
|}


The diffusion coefficients calculated using two different methods are in very close agreement and this is expected. 


'''Source of error'''

According to following equation:

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

The integration range is <math>0</math> to <math>\infty</math>. However in our calculation the integration was stopped at time<math> = 10 </math>, and this could contribute to the error of data.  Another source of error is the integration method. By using trapezium rule, we are dividing the graph into different rectangular strips. However, the graph itself can not be be describe perfectly using rectangular strips and will eventually result in overestimation (strips bigger than graph) or underestimation (strips smaller graph).

'''<span style="color:red">JC: Approximating the integral by reducing the range from infinity to 10 will cause a more significant error if the running integral hasn't fully plateaued by 10.</span>'''