File:0.001tpe.png

2015-12-04T18:45:57Z

Mcd13:

Rep:Mod:IOPJKL

2015-12-04T18:29:24Z

Mcd13:

'''Liquid Simulations Computational Experiment'''

==Theory==

===Task 1: Calculations of a Harmonic Oscillator using both Classical and Velocity-Verlet solutions===
[[File:All3graphs.png|800px|thumb|right|''Figure 1.a) Positions calculated via the two solutions vs time b) Difference between the positions calculated via different methods against time c) Total energy of the harmonic oscillator calculated using Velocity-Verlet solutions vs time '']]

'''Position calculations from the two methods''': The position <math>x</math> at time <math>t</math> was determined using <math>x(t) = A cos(\omega t + \phi)</math> where; <math>A = 1, \omega = 1</math> and <math> \phi= 0 </math>. These were plotted against time(Figure 1.a).

'''Error analysis of the two positions''': The difference between the two positions was calculated and plotted against time(Figure 1.b). The error increases with time and position, this loss of precision could come from two sources. As the position gets further from the original position the <math>a(t)dt^2</math> part of the equation has less importance also the errors could build up since the previous position is used in the preceding position calculation.

'''Total Energy of the Harmonic-Oscillator vs time''': Values for velocity <math>v</math> and position <math>x</math> were obtained from the Velocity-Verlet calculation, values for the force-constant <math>k</math> and mass <math>m</math> were both given as 1. Using <math> E= \frac{1}{2}mv^2 + \frac{1}{2}kx^2</math> the total energies were calculated and plotted against time(Figure 1.c).

===Task 2: Plot of maxima error positions vs time===
The coordinates of the 5 maximas in figure 2 were estimated and plotted on a new graph. This yielded a linear function where <math>x(t) = 0.0004x - 7*10^{-5}</math>. [[File:MaximaDiff.png|200px|thumb|right|''Figure 2. Linear plot of the maxima error positions(figure 1.b) against time '']]

===Task 3: Timesteps to minimise variation in total energy===

Multiple timesteps were used to calculate the total energy of the harmonic oscillator. The percentage difference was taken between the initial energy assumed to be correct and the energy of the first minima, these results were then plotted(figure 3). In order to keep the change in total energy below 1% a timestep equal to or below 0.2 should be used. The change in total energy of the simulation should be kept to a minimum as the total energy should not vary over time in the harmonic oscillator.[[File:%vstimestep.png|200px|thumb|right|''Figure 3. Plot of percentage energy difference and timestep'']]

===Task 4: Lennard-Jones interaction===

'''1. Calculation of distance at zero-potential, r = r0:'''
<math>\phi\left(r_0\right) = 4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right)</math>
where φ = potential energy r = the distance ε = well depth and σ = distance from zero potential distance
<math>\phi\left(r_0\right)= 4\epsilon \left( \sigma^{12} - \sigma^6 r_0^6 \right)</math>
:<math>= 4\epsilon \left( 1 - \frac{r_0^6}{\sigma^6} \right) = 0</math>
:<math>= 1 - \frac{r_0^6}{\sigma^6} = 0 </math>
:<math> \sigma^6 = r_0^6 </math>
:<math> \sigma = r_0 </math>

2. '''Force at r0:'''

:<math>\mathbf{F}_0 = - \frac{\mathrm{d}\phi\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_0} = - \frac{\mathrm{d}}{\mathrm{d}\mathbf{r}_0} 4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right) = 4\epsilon \left( -12\frac{\sigma^{12}}{r_0^{13}} + 6\frac{\sigma^6}{r_0^7} \right)</math>
Since <math>\sigma = r_0 </math>
:<math>= 4\epsilon\left( -12\frac{\sigma^{12}}{\sigma^{13}} + 6\frac{\sigma^6}{\sigma^7} \right) = 4\epsilon\frac{-6}{\sigma} = -\frac{24\epsilon}{\sigma}</math>

'''3. Finding the equilibrium separation ''req:'' '''

:<math>\mathbf{F}_{eq} = - \frac{\mathrm{d}\phi\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_{eq}} = - \frac{\mathrm{d}}{\mathrm{d}\mathbf{r}_{eq}} 4\epsilon \left( \frac{\sigma^{12}}{r_{eq}^{12}} - \frac{\sigma^6}{r_{eq}^6} \right) = 4\epsilon \left( -12\frac{\sigma^{12}}{r_{eq}^{13}} + 6\frac{\sigma^6}{r_{eq}^7} \right) = -48\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) + 24\epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right) = 0</math>
:<math>48\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) = 24\epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right)</math>
:<math>2\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) = \epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right)</math>
:<math>2\epsilon \left(\sigma^{12}\right) = \epsilon \left(\sigma^6 r_{eq}^6 \right)</math>
:<math>2\epsilon \left(\sigma^{6}\right) = \epsilon \left(\ r_{eq}^6 \right)</math>
:<math>2\left(\sigma^{6}\right) = \left(\ r_{eq}^6 \right)</math>
:<math>\sqrt[6]{2}\sigma = r_{eq}</math>

'''4. Finding the well-depth at equilibrium:'''

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

'''5. Evaluation of integrals:'''

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

:<math>\int_{2}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{2}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{2}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>= [4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_2^{\infty}</math>
:<math>= 4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(2)^{11}} - \frac{1}{-5(2)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(2)^{11}} - \frac{1}{-5(2)^5} \right) = -24.82 KJ/mol </math>

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

:<math>\int_{2.5}^\infty \phi\left(r\right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>= [4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_{2.5}^{\infty}</math>
:<math>= 4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(2.5)^{11}} - \frac{1}{-5(2.5)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(2.5)^{11}} - \frac{1}{-5(2.5)^5} \right) = -8.18 KJ/mol </math>

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

:<math>\int_{3}^\infty \phi\left(r\right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>[4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_{3}^{\infty}</math>
:<math>4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(3)^{11}} - \frac{1}{-5(3)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(3)^{11}} - \frac{1}{-5(3)^5} \right) = -3.29 KJ/mol </math>

===Task 5: Estimating the number of molecules in volumes of H2O===

1. '''Number of molecules in 1mL:'''
:<math> Mr(H_2O) = 18.02g/mol,\rho = 1g/cm^3 </math>
<math>n = \frac{m}{Mr} = \frac{1}{18.02} = 5.55x10^{-2}</math>
<math>Number of molecules = 5.55x10^{-2} x 6.022x10^{23} = 3.34192x10^{22}</math>

2. '''Volume of 1000 molecules:'''

<math>n = \frac{1000}{6.022x10^{23}} = 1.66x10^{-21} mol</math>
<math>m = 1.66x10^{-21} x 18.02 = 2.99x10^{-20}g</math>
<math>V = \frac{2.99x10^{-20}}{1} = 2.99x10^{-20}cm^3</math>

===Task 6: Periodic boundary conditions===

When the atom at <math>(0.5,0.5,0.5)</math> in the cubic simulation box <math>(0,0,0,)</math> to <math>(1,1,1)</math> moves along the vector <math>(0.7,0.6,0.2)</math> it ends up at <math>(0.2,0.1,0.7)</math> (when the periodic boundry condition is taken into consideration).

===Task 7: Converting reduced values to real values:===

1. '''Reduced distance'''

:<math> r* = 3.2 = \frac {r}{\sigma} </math>
:<math> r = 3.2 x 0.34x10^{-9} = 1.09x10^-9 = 1.09nm </math>

2. '''Well-depth'''

:<math> \frac {\epsilon}{K_B} = 120 </math>
:<math> \epsilon = 120K_B = \frac {8.69x10^{24}}{6.022x10^{23}} = 997.71 J mol^{-1} = 9.98x10^{-1} KJ mol^-1</math>

3. '''Temperature'''

:<math>T^* = 1.5 = \frac {K_BT}{\epsilon} </math>
:<math>T = \frac {T^* \epsilon}{K_B} = \frac {1.5 x 1.657x10^{-21}}{K_B} = 180K</math>

==Equilibration==

===Task 8: Problems with random starting coordinates in non-solid simulations===

If atoms were given random starting coordinates then they could end up within a contact proximity, they could bump into each other and then they would not only experience the potential of each but also the kinetic energy.

===Task 9: The number density of lattice points calculations:===

The number density ''n'' is given by the following equation

<math> n = \frac {N}{V}</math> where N is the number of atoms/lattice points per unit cell and V is volumes of the unit cell.

1. The number density of a simple cubic lattice with inter lattice point distance of 1.07722

For a simple cubic lattice there is one atom/lattice point, therefore;

:<math> n = \frac {N}{V} = \frac {N}{r^{*3}} = \frac {1}{1.07722^3} = 0.8 </math>

2. '''Side length of a face-centred cubic lattice with a number density of 1.2'''

A face-centred cubic lattice has 4 atoms/lattice points, therfore;

<math> r^* = \sqrt[3]{V} = \sqrt[3]{\frac {N}{n}} = \sqrt[3]{\frac {4}{1.2}} = 1.49380 </math>

===Task 10: How many atoms are created by create_atoms 1 box command for a face-centred cubic lattice?===

Since there are 4 atoms per lattice and the box is ten lattice spacings larger than the single lattice spacing then there would be 4000 atoms.

===Task 11: LAMMPS commands for setting the properties of atoms:===

1. Mass command

mass ''I'' ''value'' = mass 1 1.0
This command is used for setting the mass of the atom, in the general command I is the atom type and value is the mass(in reduced units). In the example the atom is type 1 and the mass is 1.0.

2. Pair_style command

pair_style style args = pair_style lj/cut 3.0

The pair_style command is used to describe paired interactions, style is the interaction which is being looked at and args is the argument used for that style. In the example above the style is defining the cutoff distance, that is the distance at which interactions become negligible. The argument for the cutoff distance is 3.0, after which the simulation will not compute.

3. pair_coeff command

pair_coeff ''* * args'' = pair_coeff * * 1.0 1.0
This command is used to define the force field coefficients between pairs of atoms, the asterix terms are to define the coefficients for all atoms, the second part of this command, args, defines the coeffecients.

===Task 12: Specifying the velocity of each atom===

The Velocity-Verlet algorithm is used to specify velocity and position
<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 \ \ (9)</math>

===Task 13: Running the simulation:===
The command $(100/timestep) command allows for the immediate evaluation of the the formula with the pre-defined variable i.e. the timestep. The following line variable n_steps equal floor(100/$({timestep}) will be broken down into individual contributions. The variable component is used to allows for the return to re-run a command, this allows creates a loop. The n_steps part describes the number of timesteps you wish to calculate which is defined in the 'Run Simulation' command. Equal is used to create a mathematical argument. Floor(11/$({timestep}) is useful in setting the maximum timestep bigger than the original one. In the run simulation part first we define the number of timesteps to be considered and then report the floor value i.e. the total time it takes. These lines are essential in creating a loop where each new timestep may be calculated subsequently.
Make it easier to change the timestep, all values dependedent on the timestep will be changed with it and so it reduces possible human errors when changing all variables and also makes it easier and faster.

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

===Task 14: Checking equilibration===

1. Pressure, Temperature and energy against time for simulation ran at 0.001 timestep

PIC fig 4

2. Equilibration

The simulation reaches equilibration at approximately 0.28, this is the timestep at which the simulation becomes relatively constant.

3. Maximising results accuracy whilst keeping the timestep within a reasonable timescale

The best timestep for maximum accuracy and sufficient information is 0.0025. It overlaps with 0.001 and so although it is more than twice the time as 0.001 it maintains a satisfactory level of accuracy. The others give equilibrium energies which are too far off the most accurate value, with the worst being for 0.01 as it diverges quite dramatically.

[[File:Energyplotsvstime.png|400px|thumb|right|''Figure 5. The energy vs time plots of simulations with varying timesteps '']]

==Running simulations under specific conditions==

===Task 15: Finding the expression for factor <math> \gamma </math>===

<math> \gamma</math> is a factor by which the velocity can be corrected.

<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT </math> (1)
<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT' </math> (2)

Where ''T''' is the target temperature and ''T'' is the instantaneous temperature

<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT </math>
<math> \frac {1}{2} m_iv_i^2 = \frac {1}{2}NK_BT </math>
<math> m_iv_i^2 = NK_BT </math>
<math> v_i^2 = \frac {NK_BT}{m_i} </math>
<math> v_i = \sqrt[2]{\frac {NK_BT}{m_i}} </math>

Insert v_i into equation x

<math> \frac {1}{2} \sum m_i \gamma^2 \frac {NK_BT}{m_i} = \frac {3}{2}NK_BT' </math>
<math> \frac {1}{2} m_i \gamma^2 \frac {NK_BT}{m_i} = \frac {1}{2}NK_BT' </math>
<math> \gamma^2NK_BT = NK_BT' </math>
<math> \gamma^2 = \frac {NK_BT'}{NK_BT} </math>
<math> \gamma^2 = \frac {T'}{T} </math>
<math> \gamma = \sqrt[2]{\frac {T'}{T}} </math>

===Task 16: The importance of number, repetition and frequency of average calculations with timesteps===

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

100 sets the number of timesteps after which another thermodynamic value is used to calculate the average, 1000 is the number of thermodynamic values which are used for the average and 1000000 sets the frequency at which averages are calculated.

100000 timesteps will be run.

===Task 17: Equation of state plots===

5 temperatures above the critical temeperature were chosen(1.6, 1.8, 2.0, 3.0 and 5.0) these were all simulated with a pressure of 2.6(''n-1'') and also 3.3(''n-2''). The timestep chose was 0.01 as in previous experiments it was the largest timesteps yielding the most accuracy.
The pressures were plotted against the temperatures along with their error bars.

The ideal gas densities were calculated with the following equation:
<math> \rho = \frac {N}{V} = \frac {P}{K_BT} = \frac {P}{T} </math> where K_B = 1 in reduced terms
[[File:DensityvsT.png|400px|thumb|right|''Figure 6. Densities n-1 and n-2 calculated via simulation and Ideal Gas methods plotted against temperature'']]

The simulated density is lower than that of the ideal gas method, this can be explained by an assumption to satisy an ideal gas. The atoms must not have any forces upon them and also they mustn't take up any space. In the simulation the atoms do take up space and therefore they have a lower density. This discrepancy increases with pressure. The general trend of the two plots is a decrease in heat capacity as temperature is increased.

==Heat Capacity Calculation==

===Task 18: The temperature dependence of the isochoric heat capacity===
[[File:CvVcommand.png|400px|thumb|left|''Figure 7. Cv/V as a function of T* at P= 0.2 and 0.8'']]

A LAAMPS script which could calculate the isochoric heat capacity from calculated thermodynamic values was used figure 7. The Cv/V calculation was done using equation x.

[[File:CvVvsT.png|400px|thumb|right|''Figure 8. Cv/V as a function of T* at P= 0.2 and 0.8'']]
The heat capacity represents the ability to transfer heat for a given temperature or how much energy will effect a given temperature.

The decrease in heat capacity at higher temperatures is a result of the non-classical behaviour of the system at low temperatures. The energy levels behave as if they are discrete, the lower the temperature the more discrete they are. As illustrated here, the energy levels are closer together at higher temperatures(<math>T_2</math>) and further apart at lower temperatures(<math>T_1</math>). When a heat increment <math>dq</math> is introduced into the two systems we see that it has a different effect in both. At T2 the dq will give a greater temperature increase than T1, this causes the denominator in the equation below to become more important thus decreasing the isochoric heat capacity.

[[Media:Cvexp.png]]

The difference between the two pressured systems in ''figure 7'' is what would be expected, the higher pressure plot has a higher isochoric heat capacity. At a higher pressure the atoms are closer together, they can thus transfer heat more effectively with a given temperature change than a system with a lower pressure since the atoms are further apart.

==Structural properties and the radial distribution function==
===Task 18: The Radial Distribution Functions of a solid, liquid and gas===

'''1. Comparisons of the solid, liquid and gas radial distributions functions as well as the integrated function '''
[[File:rdfslg.png|400px|thumb|left|figure 9]]
[[File:rdfint.png|400px|thumb|left|figure 10]]

The distribution plot(figure 9) of the gas has one short, broad peak which is close to the reference atom; as you move further away the distribution becomes constant. This characteristic shows the small probability of finding an atom close to the reference atom; also when it becomes constant this shows how there is negligible probability of finding another atom at the distance studied. When looking at the integrated function(figure 10), it emphasizes how few atoms you will find in the distance, this is characteristic of a gas' dispersed atoms.

The liquid's distribution has a taller first peak than the gas' and then multiple peaks tapering off as you move further away; it is evident that there is a higher probability of finding an atom throughout the space than in a gas. Additionally the integrated distribution(figure 10) shows a larger number of atoms as you move through the space. This demonstrates the higher proximity of the atoms in the liquid state.

Finally the solid's distribution has a very sharp and tall peak close to the reference atom but then as you move further away the peaks become smaller but continuously fluctuate. The sharpness of the peaks and their continuous distribution arise from the highly ordered structure of the solid. The peaks are sharp unlike the gas and liquid because all of the atoms are in fixed positions and so at a certain position there will be a number of atoms in an equivalent position around the sphere of the distribution function. Futhermore the solid's first peak is considerably larger than the other phases, the probibility is much greater sdue to the closeness of the atoms. The integrated function shows a very large number of atoms in the space, this again illustrates the closeness of the atoms in the solid state.

'''2. Using the solid RDF to extract information on the lattice structure'''

[[File:rdf3d.png|200px|thumb|right|figure 11]]

The three first peaks of the solid's radial distribution function are distances 1.0, 1.5 and 1.8. The lattice spacing is 1.5 and the calculated distances ra, rb and rc(where r is the distance from the reference atom to the atoms a, b and c) were found to be 1.0, 1.5 and 1.8(respectively). Therefore the first three peaks 1, 2 and 3 correspond to the lattice points a, b and c(respectively).

The coordination of atom a, b and c is 12, 6 and 24(respectively). These values were interpolated from the magnification of the integrated solid distribution function.

==Dynamical properties and the diffusion coefficient==

===Task 19: The mean squared displacement(MSD) of large and small simulations===

[[File:msdslv.png|800px|thumb|right|figure 12]]
<math> MSD(t) = <((t_0 + t)-(t_0))>^2 </math> where; <math>t_0 =</math> the initial timestep and <math>t =</math>the new timestep Equation x

The mean squared displacement(equation x) is a measure of how much random motion there is within the system. In this simulation the ability to move randomly is determined by the amount of free space and therefore the density. In the gaseous phase there is a lot of free space, the atoms move freely and encounter only a small number of atoms to collide with thus impeding it's motion; the linear increase of the plot reflects this behaviour. The liquid has a lower degree of random motion than the gas, due closer proximity of the atoms, forces such as drag will now effect the atom's motion as it moves past other atoms. The two first plots show that the atom is able to move in drift, however the liquid experiences more drag and has a smaller gradient. The solid displays the least amount of translational freedom, the close packed structure prevents any random translational motion; it does however show some random movement, arising from vibrations of the constrained atoms. Therefore the plots in figure 12 are a reasonable representation of the random movement of the three phases.

Finding the diffusion coefficient <math>D</math>

<math>D = \frac{1}{6} \frac{\delta <r^2(t)>}{\delta t}> </math> where <math>D =</math> the diffusion coefficient in <math>m s^{-2}</math>,<math><r^2(t) =</math> MSD and <math> t =</math> timestep
The gradient of the slopes give the diffusion coefficent, however this uses a dimensionless timestep denominator, therefore we must convert it to time units by dividing by the timestep(0.002).

The gradient was found by fitting a linear trendline through the function.

Solid

Conversion of the gradient to give time units gave <math>D = 4.17 x 10^{-6} ms^{-2}</math> for the million atom simulation and <math>D = 5.83E-08 x 10^{-7} ms^{-2}</math> for the 3200 atom simulation.

Liquid

Conversion of the gradient to give time units gave <math>D = 8.33 x 10^{-2} ms^{-2}</math> for the million atom simulation and <math>D = 8.33 x 10^{-2} ms^{-2}</math> for the 3200 atom simulation.

Gas

A trendline was fitted to the function, since the simulation took some time to produce the linear behavior it was fitted to be linear with the latter part of the function. Conversion of the gradient to give time units gave <math>D = 3.09 ms^{-2}</math> for the million atom simulation and <math>D = 2.43 ms^{-2}</math> for the 3200 atom simulation.

===Task 20: The velocity auto-correlation function(VACF)===

==Evaluating the 1D harmonic oscillator VACF using the expression for the position as a function of time==

The normalized VACF is given by;

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

We must simplify this expression in order for it to be useful, to do this we must first find expressions for <math>v(t)</math> and <math>v(t + \tau)</math>:

The expression for position as a function of time for the HO is defined as;
<math>x(t) = A\cos(\omega t + \phi)</math>

This can be substituted into the expression of velocity as a function of time;

<math>v(t) = \frac{dx(t)}{dt} = \frac{dA\cos(\omega t + \phi)}{dt}</math>

This is then differentiated using the chain rule;

<math> \frac{dx(t)}{dt} = \frac{dx(t)}{du} x \frac{du}{dt} </math>

Where <math>u = \omega t + \phi</math>

<math>v(t) = -\omega A\sin(\omega t + \phi) </math>

We also need our velocity as a function of time after <math>\tau</math>, this is acheived by substituting in <math>t = t + \tau</math>;

<math> v(t + \tau) = -\omega A\sin(\omega (t + \tau) + \phi)</math>
<math> = -\omega A\sin(\omega t + \omega\tau + \phi) </math>

The two velocity expressions are substituted into the VACF;

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

::::<math> = \frac{\int_{-\infty}^{\infty} \sin(\omega t + \phi) \sin(\omega t + \omega\tau + \phi)\,dt}{\int_{-\infty}^{\infty} \sin(\omega t + \phi)^2\,dt} </math>

Using the trigonometric identity <math> \sin(\alpha + \beta) = \sin(\alpha) \cos(\beta) + \cos \alpha \sin \beta </math> we can re-write our equation;

where <math> \omega t + \phi = \alpha </math> and <math> \omega \tau = \beta </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)\,dt}{ \int_{-\infty}^{\infty} \sin(\omega t + \phi)^2\,dt}</math>

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

The fraction can be expanded;

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

Then since <math>\cos(\omega \tau)</math> does not contain a t term it can be moved outside of the integral;

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

The alike terms cancel to give;

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

Finally <math>\frac{\int_{-\infty}^{\infty}\cos(\omega t + \phi) \sin(\omega \tau)\,dt}{\int_{-\infty}^{\infty}\sin(\omega t + \phi)^2\,dt} = 0</math>, this is because <math>sin</math> and <math>\frac{\cos}{\sin}</math> are odd functions. When odd functions are integrated between symmetrical boundaries the product is zero.

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

Analysis of the solid, liquid and harmonic osillator VACFs against position
[[File:VACFplots.png|800px|thumb|right|figure 13]]
The minima in the VACFs represent a change in the direction of the velocity, this caused by collisions with other atoms.

The liquid and solid VACFs are similar but have some important differences. The solid atoms correlation with their initial velocity decreases faster than the liquids. The minima for the solid is much earlier and also much more negative. This happens because the solid atoms are much closer to each other and there are a lot more of them, resulting in an interaction after less time. When the correlation becomes negative it represents and inversed correlation which is where the velocity is no longer correlated to it's initial velocity but rather external velocities of other atoms. At the minima the hanges direction after colliding fully with other atoms

cked structure, the correlation becomes more negative beause the atoms are tightly constrained and so the interactions are very strong with other atoms. The atoms in the liquid are free moving and so when they collide they can move to stabilise The decorrelated solid function shows some oscillation about <math>C(\tau)=0</math> this arises from the jostling with other atoms in the fixed lattice and so we see a sort of harmonic oscillator motion. The liquid does not display this behavior because they are free to move and so will be able to "rapidly destroy any oscillatory motion."[x] by the "diffusive motion".

===Task 19: Integrating using the trapezium rule to get D===

The plots of <math> C(\tau) vs time </math> of the three phases in a small and large system were plotted along with their running integral.
[[File:solid3200.png|800px|thumb|right|figure 14]] [[File:solidM.png|800px|thumb|right|figure 15]] [[File:liquid8000.png|800px|thumb|right|figure 16]] [[File:liquidM.png|800px|thumb|right|figure 17]] [[File:gas8000.png|800px|thumb|right|figure 18]] [[File:gasM.png|800px|thumb|right|figure 19]]

Rep:Mod:IOPJKL

2015-12-04T18:09:02Z

Mcd13: /* Heat Capacity Calculation */

'''Liquid Simulations Computational Experiment'''

==Theory==

===Task 1: Calculations of a Harmonic Oscillator using both Classical and Velocity-Verlet solutions===
[[File:All3graphs.png|800px|thumb|right|''Figure 1.a) Positions calculated via the two solutions vs time b) Difference between the positions calculated via different methods against time c) Total energy of the harmonic oscillator calculated using Velocity-Verlet solutions vs time '']]

'''Position calculations from the two methods''': The position <math>x</math> at time <math>t</math> was determined using <math>x(t) = A cos(\omega t + \phi)</math> where; <math>A = 1, \omega = 1</math> and <math> \phi= 0 </math>. These were plotted against time(Figure 1.a).

'''Error analysis of the two positions''': The difference between the two positions was calculated and plotted against time(Figure 1.b). The error increases with time and position, this loss of precision could come from two sources. As the position gets further from the original position the <math>a(t)dt^2</math> part of the equation has less importance [x] also the errors could build up since the previous position is used in the preceding position calculation.

'''Total Energy of the Harmonic-Oscillator vs time''': Values for velocity <math>v</math> and position <math>x</math> were obtained from the Velocity-Verlet calculation, values for the force-constant <math>k</math> and mass <math>m</math> were both given as 1. Using ,math> E= \frac{1}{2}mv^2 + \frac{1}{2}kx^2'' the total energies were calculated and plotted against time(Figure 1.c).

===Task 2: Plot of maxima error positions vs time===
The coordinates of the 5 maximas in figure 2 were estimated and plotted on a new graph. This yielded a linear function where <math>x(t) = 0.0004x - 7*10^{-5}</math>. [[File:MaximaDiff.png|200px|thumb|right|''Figure 2. Linear plot of the maxima error positions(figure 1.b) against time '']]

===Task 3: Timesteps to minimise variation in total energy===

Multiple timesteps were used to calculate the total energy of the harmonic oscillator. The percentage difference was taken between the initial energy assumed to be correct and the energy of the first minima, these results were then plotted(figure 3). In order to keep the change in total energy below 1% a timestep equal to or below 0.2 should be used. The change in total energy of the simulation should be kept to a minimum as the total energy should not vary over time in the harmonic oscillator.[[File:%vstimestep.png|200px|thumb|right|''Figure 3. Plot of percentage energy difference and timestep'']]

===Task 4: Lennard-Jones interaction===

'''1. Calculation of distance at zero-potential, r = r0:'''
<math>\phi\left(r_0\right) = 4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right)</math>
where φ = potential energy r = the distance ε = well depth and σ = distance from zero potential distance
<math>\phi\left(r_0\right)= 4\epsilon \left( \sigma^{12} - \sigma^6 r_0^6 \right)</math>
:<math>= 4\epsilon \left( 1 - \frac{r_0^6}{\sigma^6} \right) = 0</math>
:<math>= 1 - \frac{r_0^6}{\sigma^6} = 0 </math>
:<math> \sigma^6 = r_0^6 </math>
:<math> \sigma = r_0 </math>

2. '''Force at r0:'''

:<math>\mathbf{F}_0 = - \frac{\mathrm{d}\phi\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_0} = - \frac{\mathrm{d}}{\mathrm{d}\mathbf{r}_0} 4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right) = 4\epsilon \left( -12\frac{\sigma^{12}}{r_0^{13}} + 6\frac{\sigma^6}{r_0^7} \right)</math>
Since <math>\sigma = r_0 </math>
:<math>= 4\epsilon\left( -12\frac{\sigma^{12}}{\sigma^{13}} + 6\frac{\sigma^6}{\sigma^7} \right) = 4\epsilon\frac{-6}{\sigma} = -\frac{24\epsilon}{\sigma}</math>

'''3. Finding the equilibrium separation ''req:'' '''

:<math>\mathbf{F}_{eq} = - \frac{\mathrm{d}\phi\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_{eq}} = - \frac{\mathrm{d}}{\mathrm{d}\mathbf{r}_{eq}} 4\epsilon \left( \frac{\sigma^{12}}{r_{eq}^{12}} - \frac{\sigma^6}{r_{eq}^6} \right) = 4\epsilon \left( -12\frac{\sigma^{12}}{r_{eq}^{13}} + 6\frac{\sigma^6}{r_{eq}^7} \right) = -48\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) + 24\epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right) = 0</math>
:<math>48\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) = 24\epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right)</math>
:<math>2\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) = \epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right)</math>
:<math>2\epsilon \left(\sigma^{12}\right) = \epsilon \left(\sigma^6 r_{eq}^6 \right)</math>
:<math>2\epsilon \left(\sigma^{6}\right) = \epsilon \left(\ r_{eq}^6 \right)</math>
:<math>2\left(\sigma^{6}\right) = \left(\ r_{eq}^6 \right)</math>
:<math>\sqrt[6]{2}\sigma = r_{eq}</math>

'''4. Finding the well-depth at equilibrium:'''

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

'''5. Evaluation of integrals:'''

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

:<math>\int_{2}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{2}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{2}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>= [4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_2^{\infty}</math>
:<math>= 4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(2)^{11}} - \frac{1}{-5(2)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(2)^{11}} - \frac{1}{-5(2)^5} \right) = -24.82 KJ/mol </math>

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

:<math>\int_{2.5}^\infty \phi\left(r\right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>= [4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_{2.5}^{\infty}</math>
:<math>= 4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(2.5)^{11}} - \frac{1}{-5(2.5)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(2.5)^{11}} - \frac{1}{-5(2.5)^5} \right) = -8.18 KJ/mol </math>

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

:<math>\int_{3}^\infty \phi\left(r\right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>[4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_{3}^{\infty}</math>
:<math>4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(3)^{11}} - \frac{1}{-5(3)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(3)^{11}} - \frac{1}{-5(3)^5} \right) = -3.29 KJ/mol </math>

===Task 5: Estimating the number of molecules in volumes of H2O===

1. '''Number of molecules in 1mL:'''
:<math> Mr(H^2O) = 18.02g/mol,\rho = 1g/cm^3 </math>
<math>n = \frac{m}{Mr} = \frac{1}{18.02} = 5.55x10^{-2}</math>
<math>Number of molecules = 5.55x10^{-2} x 6.022x10^{23} = 3.34192x10^{22}</math>

2. '''Volume of 1000 molecules:'''

<math>n = \frac{1000}{6.022x10^{23}} = 1.66x10^{-21} mol</math>
<math>m = 1.66x10^{-21} x 18.02 = 2.99x10^{-20}g</math>
<math>V = \frac{2.99x10^{-20}}{1} = 2.99x10^{-20}cm^3</math>

===Task 6: Periodic boundary conditions===

When the atom at <math>(0.5,0.5,0.5)</math> in the cubic simulation box(<math>(0,0,0,) to (1,1,1)</math>) moves along the vector <math>(0.7,0.6,0.2)</math> it ends up at <math>(0.2,0.1,0.7)</math> (when the periodic boundry condition is taken into consideration).

===Task 7: Converting reduced values to real values:===

1. '''Reduced distance'''

:<math> r* = 3.2 = \frac {r}{\sigma} </math>
:<math> r = 3.2 x 0.34x10^{-9} = 1.09x10^-9 = 1.09nm </math>

2. '''Well-depth'''

:<math> \frac {\epsilon}{K_B} = 120K </math>
:<math>\epsilon = 120K_B = \frac {8.69x10^{24}}{6.022x10^{23}} = 997.71 J/mol = 9.98x10^{-1}</math>

3. '''Temperature'''

:<math> T^* = 1.5 = \frac {K_BT}{\epsilon} </math>
:<math>T = \frac {T^* \epsilon}{K_B} = \frac {1.5 x 1.657x10^{-21}}{K_B} = 180K</math>

==Equilibration==

===Task 8: Problems with random starting coordinates in non-solid simulations===

If atoms were given random starting coordinates then they could end up within a contact proximity, they could bump into each other and then they would not only experience the potential of each but also the kinetic energy.

===Task 9: The number density of lattice points calculations:===

The number density ''n'' is given by the following equation
<math> n = \frac {N}{V}</math> where N is the number of atoms/lattice points per unit cell and V is volumes of the unit cell.

1. The number density of a simple cubic lattice with inter lattice point distance of 1.07722

For a simple cubic lattice there is one atom/lattice point, using equation x
:<math> n = \frac {N}{V} = \frac {N}{r^{*3}} = \frac {1}{1.07722^3} = 0.8 </math>

2. '''Side length of a face-centred cubic lattice with a number density of 1.2'''

A face-centred cubic lattice has 4 atoms/lattice points.
<math> r^* = \sqrt[3]{V} = \sqrt[3]{\frac {N}{n}} = \sqrt[3]{\frac {4}{1.2}} = 1.49380 </math>

===Task 10: How many atoms are created by create_atoms 1 box command for a face-centred cubic lattice?===
Since there are 4 atoms per lattice and the box is ten lattice spacings larger than the single lattice spacing then there would be 4000 atoms.

===Task 11: LAMMPS commands for setting the properties of atoms:===

1. Mass command

mass ''I'' ''value'' = mass 1 1.0
This command is used for setting the mass of the atom, in the general command I is the atom type and value is the mass(in reduced units). In the example the atom is type 1 and the mass is 1.0.

2. Pair_style command

pair_style style args = pair_style lj/cut 3.0

The pair_style command is used to describe paired interactions, style is the interaction which is being looked at and args is the argument used for that style. In the example above the style is defining the cutoff distance, that is the distance at which interactions become negligible. The argument for the cutoff distance is 3.0, after which the simulation will not compute.

3. pair_coeff command

pair_coeff ''* * args'' = pair_coeff * * 1.0 1.0
This command is used to define the force field coefficients between pairs of atoms, the asterix terms are to define the coefficients for all atoms, the second part of this command, args, defines the coeffecients.

===Task 12: Specifying the velocity of each atom===

The Velocity-Verlet algorithm is used to specify velocity and position
<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 \ \ (9)</math>

===Task 13: Running the simulation:===
The command $(100/timestep) command allows for the immediate evaluation of the the formula with the pre-defined variable i.e. the timestep. The following line variable n_steps equal floor(100/$({timestep}) will be broken down into individual contributions. The variable component is used to allows for the return to re-run a command, this allows creates a loop. The n_steps part describes the number of timesteps you wish to calculate which is defined in the 'Run Simulation' command. Equal is used to create a mathematical argument. Floor(11/$({timestep}) is useful in setting the maximum timestep bigger than the original one. In the run simulation part first we define the number of timesteps to be considered and then report the floor value i.e. the total time it takes. These lines are essential in creating a loop where each new timestep may be calculated subsequently.
Make it easier to change the timestep, all values dependedent on the timestep will be changed with it and so it reduces possible human errors when changing all variables and also makes it easier and faster.

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

===Task 14: Checking equilibration===

1. Pressure, Temperature and energy against time for simulation ran at 0.001 timestep

PIC fig 4

2. Equilibration

The simulation reaches equilibration at approximately 0.28, this is the timestep at which the simulation becomes relatively constant.

3. Maximising results accuracy whilst keeping the timestep within a reasonable timescale
The best timestep for maximum accuracy and sufficient information is 0.0025. It overlaps with 0.001 and so although it is more than twice the time as 0.001 it maintains a satisfactory level of accuracy. The others give equilibrium energies which are too far off the most accurate value, with the worst being for 0.01 as it diverges quite dramatically.

[[File:Energyplotsvstime.png|400px|thumb|right|''Figure 5. The energy vs time plots of simulations with varying timesteps '']]

==Running simulations under specific conditions==

===Task 15: Finding the expression for factor <math> \gamma </math>===

<math> \gamma</math> is a factor by which the velocity can be corrected.

<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT </math> (1)
<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT' </math> (2)

Where ''T''' is the target temperature and ''T'' is the instantaneous temperature

<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT </math>
<math> \frac {1}{2} m_iv_i^2 = \frac {1}{2}NK_BT </math>
<math> m_iv_i^2 = NK_BT </math>
<math> v_i^2 = \frac {NK_BT}{m_i} </math>
<math> v_i = \sqrt[2]{\frac {NK_BT}{m_i}} </math>

Insert v_i into equation x

<math> \frac {1}{2} \sum m_i \gamma^2 \frac {NK_BT}{m_i} = \frac {3}{2}NK_BT' </math>
<math> \frac {1}{2} m_i \gamma^2 \frac {NK_BT}{m_i} = \frac {1}{2}NK_BT' </math>
<math> \gamma^2NK_BT = NK_BT' </math>
<math> \gamma^2 = \frac {NK_BT'}{NK_BT} </math>
<math> \gamma^2 = \frac {T'}{T} </math>
<math> \gamma = \sqrt[2]{\frac {T'}{T}} </math>

===Task 16: The importance of number, repetition and frequency of average calculations with timesteps===

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

100 sets the number of timesteps after which another thermodynamic value is used to calculate the average, 1000 is the number of thermodynamic values which are used for the average and 1000000 sets the frequency at which averages are calculated.

100000 timesteps will be run.

===Task 17: Equation of state plots===

5 temperatures above the critical temeperature were chosen(1.6, 1.8, 2.0, 3.0 and 5.0) these were all simulated with a pressure of 2.6(''n-1'') and also 3.3(''n-2''). The timestep chose was 0.01 as in previous experiments it was the largest timesteps yielding the most accuracy.
The pressures were plotted against the temperatures along with their error bars.

The ideal gas densities were calculated with the following equation:
<math> \rho = \frac {N}{V} = \frac {P}{K_BT} = \frac {P}{T} </math> where K_B = 1 in reduced terms
[[File:DensityvsT.png|400px|thumb|right|''Figure 6. Densities n-1 and n-2 calculated via simulation and Ideal Gas methods plotted against temperature'']]

The simulated density is lower than that of the ideal gas method, this can be explained by an assumption to satisy an ideal gas. The atoms must not have any forces upon them and also they mustn't take up any space. In the simulation the atoms do take up space and therefore they have a lower density. This discrepancy increases with pressure. The general trend of the two plots is a decrease in heat capacity as temperature is increased.

==Heat Capacity Calculation==

===Task 18: The temperature dependence of the isochoric heat capacity===
[[File:CvVcommand.png|400px|thumb|left|''Figure 7. Cv/V as a function of T* at P= 0.2 and 0.8'']]

A LAAMPS script which could calculate the isochoric heat capacity from calculated thermodynamic values was used figure 7. The Cv/V calculation was done using equation x.

[[File:CvVvsT.png|400px|thumb|right|''Figure 8. Cv/V as a function of T* at P= 0.2 and 0.8'']]
The heat capacity represents the ability to transfer heat for a given temperature or how much energy will effect a given temperature.

The decrease in heat capacity at higher temperatures is a result of the non-classical behaviour of the system at low temperatures. The energy levels behave as if they are discrete, the lower the temperature the more discrete they are. As illustrated here, the energy levels are closer together at higher temperatures(<math>T_2</math>) and further apart at lower temperatures(<math>T_1</math>). When a heat increment <math>dq</math> is introduced into the two systems we see that it has a different effect in both. At T2 the dq will give a greater temperature increase than T1, this causes the denominator in the equation below to become more important thus decreasing the isochoric heat capacity.

[[Media:Cvexp.png]]

The difference between the two pressured systems in ''figure 7'' is what would be expected, the higher pressure plot has a higher isochoric heat capacity. At a higher pressure the atoms are closer together, they can thus transfer heat more effectively with a given temperature change than a system with a lower pressure since the atoms are further apart.

==Structural properties and the radial distribution function==
===Task 18: The Radial Distribution Functions of a solid, liquid and gas===

'''1. Comparisons of the solid, liquid and gas radial distributions functions as well as the integrated function '''
[[File:rdfslg.png|400px|thumb|left|figure 9]]
[[File:rdfint.png|400px|thumb|left|figure 10]]

The distribution plot(figure 9) of the gas has one short, broad peak which is close to the reference atom; as you move further away the distribution becomes constant. This characteristic shows the small probability of finding an atom close to the reference atom; also when it becomes constant this shows how there is negligible probability of finding another atom at the distance studied. When looking at the integrated function(figure 10), it emphasizes how few atoms you will find in the distance, this is characteristic of a gas' dispersed atoms.

The liquid's distribution has a taller first peak than the gas' and then multiple peaks tapering off as you move further away; it is evident that there is a higher probability of finding an atom throughout the space than in a gas. Additionally the integrated distribution(figure 10) shows a larger number of atoms as you move through the space. This demonstrates the higher proximity of the atoms in the liquid state.

Finally the solid's distribution has a very sharp and tall peak close to the reference atom but then as you move further away the peaks become smaller but continuously fluctuate. The sharpness of the peaks and their continuous distribution arise from the highly ordered structure of the solid. The peaks are sharp unlike the gas and liquid because all of the atoms are in fixed positions and so at a certain position there will be a number of atoms in an equivalent position around the sphere of the distribution function. Futhermore the solid's first peak is considerably larger than the other phases, the probibility is much greater sdue to the closeness of the atoms. The integrated function shows a very large number of atoms in the space, this again illustrates the closeness of the atoms in the solid state.

'''2. Using the solid RDF to extract information on the lattice structure'''

[[File:rdf3d.png|200px|thumb|right|figure 11]]

The three first peaks of the solid's radial distribution function are distances 1.0, 1.5 and 1.8. The lattice spacing is 1.5 and the calculated distances ra, rb and rc(where r is the distance from the reference atom to the atoms a, b and c) were found to be 1.0, 1.5 and 1.8(respectively). Therefore the first three peaks 1, 2 and 3 correspond to the lattice points a, b and c(respectively).

The coordination of atom a, b and c is 12, 6 and 24(respectively). These values were interpolated from the magnification of the integrated solid distribution function.

==Dynamical properties and the diffusion coefficient==

===Task 19: The mean squared displacement(MSD) of large and small simulations===

[[File:msdslv.png|800px|thumb|right|figure 12]]
<math> MSD(t) = <((t_0 + t)-(t_0))>^2 </math> where; <math>t_0 =</math> the initial timestep and <math>t =</math>the new timestep Equation x

The mean squared displacement(equation x) is a measure of how much random motion there is within the system. In this simulation the ability to move randomly is determined by the amount of free space and therefore the density. In the gaseous phase there is a lot of free space, the atoms move freely and encounter only a small number of atoms to collide with thus impeding it's motion; the linear increase of the plot reflects this behaviour. The liquid has a lower degree of random motion than the gas, due closer proximity of the atoms, forces such as drag will now effect the atom's motion as it moves past other atoms. The two first plots show that the atom is able to move in drift, however the liquid experiences more drag and has a smaller gradient. The solid displays the least amount of translational freedom, the close packed structure prevents any random translational motion; it does however show some random movement, arising from vibrations of the constrained atoms. Therefore the plots in figure 12 are a reasonable representation of the random movement of the three phases.

Finding the diffusion coefficient <math>D</math>

<math>D = \frac{1}{6} \frac{\delta <r^2(t)>}{\delta t}> </math> where <math>D =</math> the diffusion coefficient in <math>m s^{-2}</math>,<math><r^2(t) =</math> MSD and <math> t =</math> timestep
The gradient of the slopes give the diffusion coefficent, however this uses a dimensionless timestep denominator, therefore we must convert it to time units by dividing by the timestep(0.002).

The gradient was found by fitting a linear trendline through the function.

Solid

Conversion of the gradient to give time units gave <math>D = 4.17 x 10^{-6} ms^{-2}</math> for the million atom simulation and <math>D = 5.83E-08 x 10^{-7} ms^{-2}</math> for the 3200 atom simulation.

Liquid

Conversion of the gradient to give time units gave <math>D = 8.33 x 10^{-2} ms^{-2}</math> for the million atom simulation and <math>D = 8.33 x 10^{-2} ms^{-2}</math> for the 3200 atom simulation.

Gas

A trendline was fitted to the function, since the simulation took some time to produce the linear behavior it was fitted to be linear with the latter part of the function. Conversion of the gradient to give time units gave <math>D = 3.09 ms^{-2}</math> for the million atom simulation and <math>D = 2.43 ms^{-2}</math> for the 3200 atom simulation.

===Task 20: The velocity auto-correlation function(VACF)===

==Evaluating the 1D harmonic oscillator VACF using the expression for the position as a function of time==

The normalized VACF is given by;

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

We must simplify this expression in order for it to be useful, to do this we must first find expressions for <math>v(t)</math> and <math>v(t + \tau)</math>:

The expression for position as a function of time for the HO is defined as;
<math>x(t) = A\cos(\omega t + \phi)</math>

This can be substituted into the expression of velocity as a function of time;

<math>v(t) = \frac{dx(t)}{dt} = \frac{dA\cos(\omega t + \phi)}{dt}</math>

This is then differentiated using the chain rule;

<math> \frac{dx(t)}{dt} = \frac{dx(t)}{du} x \frac{du}{dt} </math>

Where <math>u = \omega t + \phi</math>

<math>v(t) = -\omega A\sin(\omega t + \phi) </math>

We also need our velocity as a function of time after <math>\tau</math>, this is acheived by substituting in <math>t = t + \tau</math>;

<math> v(t + \tau) = -\omega A\sin(\omega (t + \tau) + \phi)</math>
<math> = -\omega A\sin(\omega t + \omega\tau + \phi) </math>

The two velocity expressions are substituted into the VACF;

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

::::<math> = \frac{\int_{-\infty}^{\infty} \sin(\omega t + \phi) \sin(\omega t + \omega\tau + \phi)\,dt}{\int_{-\infty}^{\infty} \sin(\omega t + \phi)^2\,dt} </math>

Using the trigonometric identity <math> \sin(\alpha + \beta) = \sin(\alpha) \cos(\beta) + \cos \alpha \sin \beta </math> we can re-write our equation;

where <math> \omega t + \phi = \alpha </math> and <math> \omega \tau = \beta </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)\,dt}{ \int_{-\infty}^{\infty} \sin(\omega t + \phi)^2\,dt}</math>

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

The fraction can be expanded;

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

Then since <math>\cos(\omega \tau)</math> does not contain a t term it can be moved outside of the integral;

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

The alike terms cancel to give;

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

Finally <math>\frac{\int_{-\infty}^{\infty}\cos(\omega t + \phi) \sin(\omega \tau)\,dt}{\int_{-\infty}^{\infty}\sin(\omega t + \phi)^2\,dt} = 0</math>, this is because <math>sin</math> and <math>\frac{\cos}{\sin}</math> are odd functions. When odd functions are integrated between symmetrical boundaries the product is zero.

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

Analysis of the solid, liquid and harmonic osillator VACFs against position
[[File:VACFplots.png|800px|thumb|right|figure 13]]
The minima in the VACFs represent a change in the direction of the velocity, this caused by collisions with other atoms.

The liquid and solid VACFs are similar but have some important differences. The solid atoms correlation with their initial velocity decreases faster than the liquids. The minima for the solid is much earlier and also much more negative. This happens because the solid atoms are much closer to each other and there are a lot more of them, resulting in an interaction after less time. When the correlation becomes negative it represents and inversed correlation which is where the velocity is no longer correlated to it's initial velocity but rather external velocities of other atoms. At the minima the hanges direction after colliding fully with other atoms

cked structure, the correlation becomes more negative beause the atoms are tightly constrained and so the interactions are very strong with other atoms. The atoms in the liquid are free moving and so when they collide they can move to stabilise The decorrelated solid function shows some oscillation about <math>C(\tau)=0</math> this arises from the jostling with other atoms in the fixed lattice and so we see a sort of harmonic oscillator motion. The liquid does not display this behavior because they are free to move and so will be able to "rapidly destroy any oscillatory motion."[x] by the "diffusive motion".

===Task 19: Integrating using the trapezium rule to get D===

The plots of <math> C(\tau) vs time </math> of the three phases in a small and large system were plotted along with their running integral.
[[File:solid3200.png|800px|thumb|right|figure 14]] [[File:solidM.png|800px|thumb|right|figure 15]] [[File:liquid8000.png|800px|thumb|right|figure 16]] [[File:liquidM.png|800px|thumb|right|figure 17]] [[File:gas8000.png|800px|thumb|right|figure 18]] [[File:gasM.png|800px|thumb|right|figure 19]]

File:CvExp.png

2015-12-04T18:03:56Z

Mcd13:

Rep:Mod:IOPJKL

2015-12-04T14:01:05Z

Mcd13: /* Evaluating the 1D harmonic oscillator VACF using the expression for the position as a function of time */

'''Liquid Simulations Computational Experiment'''

==Theory==

===Task 1: Calculations of a Harmonic Oscillator using both Classical and Velocity-Verlet solutions===
[[File:All3graphs.png|800px|thumb|right|''Figure 1.a) Positions calculated via the two solutions vs time b) Difference between the positions calculated via different methods against time c) Total energy of the harmonic oscillator calculated using Velocity-Verlet solutions vs time '']]

'''Position calculations from the two methods''': The position <math>x</math> at time <math>t</math> was determined using <math>x(t) = A cos(\omega t + \phi)</math> where; <math>A = 1, \omega = 1</math> and <math> \phi= 0 </math>. These were plotted against time(Figure 1.a).

'''Error analysis of the two positions''': The difference between the two positions was calculated and plotted against time(Figure 1.b). The error increases with time and position, this loss of precision could come from two sources. As the position gets further from the original position the <math>a(t)dt^2</math> part of the equation has less importance [x] also the errors could build up since the previous position is used in the preceding position calculation.

'''Total Energy of the Harmonic-Oscillator vs time''': Values for velocity <math>v</math> and position <math>x</math> were obtained from the Velocity-Verlet calculation, values for the force-constant <math>k</math> and mass <math>m</math> were both given as 1. Using ,math> E= \frac{1}{2}mv^2 + \frac{1}{2}kx^2'' the total energies were calculated and plotted against time(Figure 1.c).

===Task 2: Plot of maxima error positions vs time===
The coordinates of the 5 maximas in figure 2 were estimated and plotted on a new graph. This yielded a linear function where <math>x(t) = 0.0004x - 7*10^{-5}</math>. [[File:MaximaDiff.png|200px|thumb|right|''Figure 2. Linear plot of the maxima error positions(figure 1.b) against time '']]

===Task 3: Timesteps to minimise variation in total energy===

Multiple timesteps were used to calculate the total energy of the harmonic oscillator. The percentage difference was taken between the initial energy assumed to be correct and the energy of the first minima, these results were then plotted(figure 3). In order to keep the change in total energy below 1% a timestep equal to or below 0.2 should be used. The change in total energy of the simulation should be kept to a minimum as the total energy should not vary over time in the harmonic oscillator.[[File:%vstimestep.png|200px|thumb|right|''Figure 3. Plot of percentage energy difference and timestep'']]

===Task 4: Lennard-Jones interaction===

'''1. Calculation of distance at zero-potential, r = r0:'''
<math>\phi\left(r_0\right) = 4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right)</math>
where φ = potential energy r = the distance ε = well depth and σ = distance from zero potential distance
<math>\phi\left(r_0\right)= 4\epsilon \left( \sigma^{12} - \sigma^6 r_0^6 \right)</math>
:<math>= 4\epsilon \left( 1 - \frac{r_0^6}{\sigma^6} \right) = 0</math>
:<math>= 1 - \frac{r_0^6}{\sigma^6} = 0 </math>
:<math> \sigma^6 = r_0^6 </math>
:<math> \sigma = r_0 </math>

2. '''Force at r0:'''

:<math>\mathbf{F}_0 = - \frac{\mathrm{d}\phi\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_0} = - \frac{\mathrm{d}}{\mathrm{d}\mathbf{r}_0} 4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right) = 4\epsilon \left( -12\frac{\sigma^{12}}{r_0^{13}} + 6\frac{\sigma^6}{r_0^7} \right)</math>
Since <math>\sigma = r_0 </math>
:<math>= 4\epsilon\left( -12\frac{\sigma^{12}}{\sigma^{13}} + 6\frac{\sigma^6}{\sigma^7} \right) = 4\epsilon\frac{-6}{\sigma} = -\frac{24\epsilon}{\sigma}</math>

'''3. Finding the equilibrium separation ''req:'' '''

:<math>\mathbf{F}_{eq} = - \frac{\mathrm{d}\phi\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_{eq}} = - \frac{\mathrm{d}}{\mathrm{d}\mathbf{r}_{eq}} 4\epsilon \left( \frac{\sigma^{12}}{r_{eq}^{12}} - \frac{\sigma^6}{r_{eq}^6} \right) = 4\epsilon \left( -12\frac{\sigma^{12}}{r_{eq}^{13}} + 6\frac{\sigma^6}{r_{eq}^7} \right) = -48\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) + 24\epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right) = 0</math>
:<math>48\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) = 24\epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right)</math>
:<math>2\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) = \epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right)</math>
:<math>2\epsilon \left(\sigma^{12}\right) = \epsilon \left(\sigma^6 r_{eq}^6 \right)</math>
:<math>2\epsilon \left(\sigma^{6}\right) = \epsilon \left(\ r_{eq}^6 \right)</math>
:<math>2\left(\sigma^{6}\right) = \left(\ r_{eq}^6 \right)</math>
:<math>\sqrt[6]{2}\sigma = r_{eq}</math>

'''4. Finding the well-depth at equilibrium:'''

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

'''5. Evaluation of integrals:'''

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

:<math>\int_{2}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{2}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{2}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>= [4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_2^{\infty}</math>
:<math>= 4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(2)^{11}} - \frac{1}{-5(2)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(2)^{11}} - \frac{1}{-5(2)^5} \right) = -24.82 KJ/mol </math>

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

:<math>\int_{2.5}^\infty \phi\left(r\right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>= [4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_{2.5}^{\infty}</math>
:<math>= 4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(2.5)^{11}} - \frac{1}{-5(2.5)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(2.5)^{11}} - \frac{1}{-5(2.5)^5} \right) = -8.18 KJ/mol </math>

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

:<math>\int_{3}^\infty \phi\left(r\right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>[4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_{3}^{\infty}</math>
:<math>4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(3)^{11}} - \frac{1}{-5(3)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(3)^{11}} - \frac{1}{-5(3)^5} \right) = -3.29 KJ/mol </math>

===Task 5: Estimating the number of molecules in volumes of H2O===

1. '''Number of molecules in 1mL:'''
:<math> Mr(H^2O) = 18.02g/mol,\rho = 1g/cm^3 </math>
<math>n = \frac{m}{Mr} = \frac{1}{18.02} = 5.55x10^{-2}</math>
<math>Number of molecules = 5.55x10^{-2} x 6.022x10^{23} = 3.34192x10^{22}</math>

2. '''Volume of 1000 molecules:'''

<math>n = \frac{1000}{6.022x10^{23}} = 1.66x10^{-21} mol</math>
<math>m = 1.66x10^{-21} x 18.02 = 2.99x10^{-20}g</math>
<math>V = \frac{2.99x10^{-20}}{1} = 2.99x10^{-20}cm^3</math>

===Task 6: Periodic boundary conditions===

When the atom at <math>(0.5,0.5,0.5)</math> in the cubic simulation box(<math>(0,0,0,) to (1,1,1)</math>) moves along the vector <math>(0.7,0.6,0.2)</math> it ends up at <math>(0.2,0.1,0.7)</math> (when the periodic boundry condition is taken into consideration).

===Task 7: Converting reduced values to real values:===

1. '''Reduced distance'''

:<math> r* = 3.2 = \frac {r}{\sigma} </math>
:<math> r = 3.2 x 0.34x10^{-9} = 1.09x10^-9 = 1.09nm </math>

2. '''Well-depth'''

:<math> \frac {\epsilon}{K_B} = 120K </math>
:<math>\epsilon = 120K_B = \frac {8.69x10^{24}}{6.022x10^{23}} = 997.71 J/mol = 9.98x10^{-1}</math>

3. '''Temperature'''

:<math> T^* = 1.5 = \frac {K_BT}{\epsilon} </math>
:<math>T = \frac {T^* \epsilon}{K_B} = \frac {1.5 x 1.657x10^{-21}}{K_B} = 180K</math>

==Equilibration==

===Task 8: Problems with random starting coordinates in non-solid simulations===

If atoms were given random starting coordinates then they could end up within a contact proximity, they could bump into each other and then they would not only experience the potential of each but also the kinetic energy.

===Task 9: The number density of lattice points calculations:===

The number density ''n'' is given by the following equation
<math> n = \frac {N}{V}</math> where N is the number of atoms/lattice points per unit cell and V is volumes of the unit cell.

1. The number density of a simple cubic lattice with inter lattice point distance of 1.07722

For a simple cubic lattice there is one atom/lattice point, using equation x
:<math> n = \frac {N}{V} = \frac {N}{r^{*3}} = \frac {1}{1.07722^3} = 0.8 </math>

2. '''Side length of a face-centred cubic lattice with a number density of 1.2'''

A face-centred cubic lattice has 4 atoms/lattice points.
<math> r^* = \sqrt[3]{V} = \sqrt[3]{\frac {N}{n}} = \sqrt[3]{\frac {4}{1.2}} = 1.49380 </math>

===Task 10: How many atoms are created by create_atoms 1 box command for a face-centred cubic lattice?===
Since there are 4 atoms per lattice and the box is ten lattice spacings larger than the single lattice spacing then there would be 4000 atoms.

===Task 11: LAMMPS commands for setting the properties of atoms:===

1. Mass command

mass ''I'' ''value'' = mass 1 1.0
This command is used for setting the mass of the atom, in the general command I is the atom type and value is the mass(in reduced units). In the example the atom is type 1 and the mass is 1.0.

2. Pair_style command

pair_style style args = pair_style lj/cut 3.0

The pair_style command is used to describe paired interactions, style is the interaction which is being looked at and args is the argument used for that style. In the example above the style is defining the cutoff distance, that is the distance at which interactions become negligible. The argument for the cutoff distance is 3.0, after which the simulation will not compute.

3. pair_coeff command

pair_coeff ''* * args'' = pair_coeff * * 1.0 1.0
This command is used to define the force field coefficients between pairs of atoms, the asterix terms are to define the coefficients for all atoms, the second part of this command, args, defines the coeffecients.

===Task 12: Specifying the velocity of each atom===

The Velocity-Verlet algorithm is used to specify velocity and position
<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 \ \ (9)</math>

===Task 13: Running the simulation:===
The command $(100/timestep) command allows for the immediate evaluation of the the formula with the pre-defined variable i.e. the timestep. The following line variable n_steps equal floor(100/$({timestep}) will be broken down into individual contributions. The variable component is used to allows for the return to re-run a command, this allows creates a loop. The n_steps part describes the number of timesteps you wish to calculate which is defined in the 'Run Simulation' command. Equal is used to create a mathematical argument. Floor(11/$({timestep}) is useful in setting the maximum timestep bigger than the original one. In the run simulation part first we define the number of timesteps to be considered and then report the floor value i.e. the total time it takes. These lines are essential in creating a loop where each new timestep may be calculated subsequently.
Make it easier to change the timestep, all values dependedent on the timestep will be changed with it and so it reduces possible human errors when changing all variables and also makes it easier and faster.

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

===Task 14: Checking equilibration===

1. Pressure, Temperature and energy against time for simulation ran at 0.001 timestep

PIC fig 4

2. Equilibration

The simulation reaches equilibration at approximately 0.28, this is the timestep at which the simulation becomes relatively constant.

3. Maximising results accuracy whilst keeping the timestep within a reasonable timescale
The best timestep for maximum accuracy and sufficient information is 0.0025. It overlaps with 0.001 and so although it is more than twice the time as 0.001 it maintains a satisfactory level of accuracy. The others give equilibrium energies which are too far off the most accurate value, with the worst being for 0.01 as it diverges quite dramatically.

[[File:Energyplotsvstime.png|400px|thumb|right|''Figure 5. The energy vs time plots of simulations with varying timesteps '']]

==Running simulations under specific conditions==

===Task 15: Finding the expression for factor <math> \gamma </math>===

<math> \gamma</math> is a factor by which the velocity can be corrected.

<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT </math> (1)
<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT' </math> (2)

Where ''T''' is the target temperature and ''T'' is the instantaneous temperature

<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT </math>
<math> \frac {1}{2} m_iv_i^2 = \frac {1}{2}NK_BT </math>
<math> m_iv_i^2 = NK_BT </math>
<math> v_i^2 = \frac {NK_BT}{m_i} </math>
<math> v_i = \sqrt[2]{\frac {NK_BT}{m_i}} </math>

Insert v_i into equation x

<math> \frac {1}{2} \sum m_i \gamma^2 \frac {NK_BT}{m_i} = \frac {3}{2}NK_BT' </math>
<math> \frac {1}{2} m_i \gamma^2 \frac {NK_BT}{m_i} = \frac {1}{2}NK_BT' </math>
<math> \gamma^2NK_BT = NK_BT' </math>
<math> \gamma^2 = \frac {NK_BT'}{NK_BT} </math>
<math> \gamma^2 = \frac {T'}{T} </math>
<math> \gamma = \sqrt[2]{\frac {T'}{T}} </math>

===Task 16: The importance of number, repetition and frequency of average calculations with timesteps===

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

100 sets the number of timesteps after which another thermodynamic value is used to calculate the average, 1000 is the number of thermodynamic values which are used for the average and 1000000 sets the frequency at which averages are calculated.

100000 timesteps will be run.

===Task 17: Equation of state plots===

5 temperatures above the critical temeperature were chosen(1.6, 1.8, 2.0, 3.0 and 5.0) these were all simulated with a pressure of 2.6(''n-1'') and also 3.3(''n-2''). The timestep chose was 0.01 as in previous experiments it was the largest timesteps yielding the most accuracy.
The pressures were plotted against the temperatures along with their error bars.

The ideal gas densities were calculated with the following equation:
<math> \rho = \frac {N}{V} = \frac {P}{K_BT} = \frac {P}{T} </math> where K_B = 1 in reduced terms
[[File:DensityvsT.png|400px|thumb|right|''Figure 6. Densities n-1 and n-2 calculated via simulation and Ideal Gas methods plotted against temperature'']]

The simulated density is lower than that of the ideal gas method, this can be explained by an assumption to satisy an ideal gas. The atoms must not have any forces upon them and also they mustn't take up any space. In the simulation the atoms do take up space and therefore they have a lower density. This discrepancy increases with pressure. The general trend of the two plots is a decrease in heat capacity as temperature is increased.

==Heat Capacity Calculation==

===Task 18: The temperature dependence of the isochoric heat capacity===
[[File:CvVcommand.png|400px|thumb|left|''Figure 7. Cv/V as a function of T* at P= 0.2 and 0.8'']]

A LAAMPS script which could calculate the isochoric heat capacity from calculated thermodynamic values was used figure 7. The Cv/V calculation was done using equation x.

[[File:CvVvsT.png|400px|thumb|right|''Figure 8. Cv/V as a function of T* at P= 0.2 and 0.8'']]
The heat capacity represents the ability to transfer heat for a given temperature or how much energy will effect a given temperature.

The general decrease of Cv/V vs T is what would be expected. At a higher temperature the energy levels are closer together and at a lower temperature they are further apart. A heat increment will have a greater effect on the higher temperature system, it will cause a greater temperature change since more energy levels are accessible, resulting in a low heat capacity. At a lower temperature the same heat increment that caused a large temperature change at higher temperature will not have the same effect, less energy levels are accessible and so the heat capacity will be higher.

The difference between the two pressured systems in ''figure 7'' is what would be expected, the higher pressure plot has a higher isochoric heat capacity. At a higher pressure the atoms are closer together, they can thus transfer heat more efficiently with a given temperature change than a system with a lower temperature since the atoms are further away.

==Structural properties and the radial distribution function==
===Task 18: The Radial Distribution Functions of a solid, liquid and gas===

'''1. Comparisons of the solid, liquid and gas radial distributions functions as well as the integrated function '''
[[File:rdfslg.png|400px|thumb|left|figure 9]]
[[File:rdfint.png|400px|thumb|left|figure 10]]

The distribution plot(figure 9) of the gas has one short, broad peak which is close to the reference atom; as you move further away the distribution becomes constant. This characteristic shows the small probability of finding an atom close to the reference atom; also when it becomes constant this shows how there is negligible probability of finding another atom at the distance studied. When looking at the integrated function(figure 10), it emphasizes how few atoms you will find in the distance, this is characteristic of a gas' dispersed atoms.

The liquid's distribution has a taller first peak than the gas' and then multiple peaks tapering off as you move further away; it is evident that there is a higher probability of finding an atom throughout the space than in a gas. Additionally the integrated distribution(figure 10) shows a larger number of atoms as you move through the space. This demonstrates the higher proximity of the atoms in the liquid state.

Finally the solid's distribution has a very sharp and tall peak close to the reference atom but then as you move further away the peaks become smaller but continuously fluctuate. The sharpness of the peaks and their continuous distribution arise from the highly ordered structure of the solid. The peaks are sharp unlike the gas and liquid because all of the atoms are in fixed positions and so at a certain position there will be a number of atoms in an equivalent position around the sphere of the distribution function. Futhermore the solid's first peak is considerably larger than the other phases, the probibility is much greater sdue to the closeness of the atoms. The integrated function shows a very large number of atoms in the space, this again illustrates the closeness of the atoms in the solid state.

'''2. Using the solid RDF to extract information on the lattice structure'''

[[File:rdf3d.png|200px|thumb|right|figure 11]]

The three first peaks of the solid's radial distribution function are distances 1.0, 1.5 and 1.8. The lattice spacing is 1.5 and the calculated distances ra, rb and rc(where r is the distance from the reference atom to the atoms a, b and c) were found to be 1.0, 1.5 and 1.8(respectively). Therefore the first three peaks 1, 2 and 3 correspond to the lattice points a, b and c(respectively).

The coordination of atom a, b and c is 12, 6 and 24(respectively). These values were interpolated from the magnification of the integrated solid distribution function.

==Dynamical properties and the diffusion coefficient==

===Task 19: The mean squared displacement(MSD) of large and small simulations===

[[File:msdslv.png|800px|thumb|right|figure 12]]
<math> MSD(t) = <((t_0 + t)-(t_0))>^2 </math> where; <math>t_0 =</math> the initial timestep and <math>t =</math>the new timestep Equation x

The mean squared displacement(equation x) is a measure of how much random motion there is within the system. In this simulation the ability to move randomly is determined by the amount of free space and therefore the density. In the gaseous phase there is a lot of free space, the atoms move freely and encounter only a small number of atoms to collide with thus impeding it's motion; the linear increase of the plot reflects this behaviour. The liquid has a lower degree of random motion than the gas, due closer proximity of the atoms, forces such as drag will now effect the atom's motion as it moves past other atoms. The two first plots show that the atom is able to move in drift, however the liquid experiences more drag and has a smaller gradient. The solid displays the least amount of translational freedom, the close packed structure prevents any random translational motion; it does however show some random movement, arising from vibrations of the constrained atoms. Therefore the plots in figure 12 are a reasonable representation of the random movement of the three phases.

Finding the diffusion coefficient <math>D</math>

<math>D = \frac{1}{6} \frac{\delta <r^2(t)>}{\delta t}> </math> where <math>D =</math> the diffusion coefficient in <math>m s^{-2}</math>,<math><r^2(t) =</math> MSD and <math> t =</math> timestep
The gradient of the slopes give the diffusion coefficent, however this uses a dimensionless timestep denominator, therefore we must convert it to time units by dividing by the timestep(0.002).

The gradient was found by fitting a linear trendline through the function.

Solid

Conversion of the gradient to give time units gave <math>D = 4.17 x 10^{-6} ms^{-2}</math> for the million atom simulation and <math>D = 5.83E-08 x 10^{-7} ms^{-2}</math> for the 3200 atom simulation.

Liquid

Conversion of the gradient to give time units gave <math>D = 8.33 x 10^{-2} ms^{-2}</math> for the million atom simulation and <math>D = 8.33 x 10^{-2} ms^{-2}</math> for the 3200 atom simulation.

Gas

A trendline was fitted to the function, since the simulation took some time to produce the linear behavior it was fitted to be linear with the latter part of the function. Conversion of the gradient to give time units gave <math>D = 3.09 ms^{-2}</math> for the million atom simulation and <math>D = 2.43 ms^{-2}</math> for the 3200 atom simulation.

===Task 20: The velocity auto-correlation function(VACF)===

==Evaluating the 1D harmonic oscillator VACF using the expression for the position as a function of time==

The normalized VACF is given by;

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

We must simplify this expression in order for it to be useful, to do this we must first find expressions for <math>v(t)</math> and <math>v(t + \tau)</math>:

The expression for position as a function of time for the HO is defined as;
<math>x(t) = A\cos(\omega t + \phi)</math>

This can be substituted into the expression of velocity as a function of time;

<math>v(t) = \frac{dx(t)}{dt} = \frac{dA\cos(\omega t + \phi)}{dt}</math>

This is then differentiated using the chain rule;

<math> \frac{dx(t)}{dt} = \frac{dx(t)}{du} x \frac{du}{dt} </math>

Where <math>u = \omega t + \phi</math>

<math>v(t) = -\omega A\sin(\omega t + \phi) </math>

We also need our velocity as a function of time after <math>\tau</math>, this is acheived by substituting in <math>t = t + \tau</math>;

<math> v(t + \tau) = -\omega A\sin(\omega (t + \tau) + \phi)</math>
<math> = -\omega A\sin(\omega t + \omega\tau + \phi) </math>

The two velocity expressions are substituted into the VACF;

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

::::<math> = \frac{\int_{-\infty}^{\infty} \sin(\omega t + \phi) \sin(\omega t + \omega\tau + \phi)\,dt}{\int_{-\infty}^{\infty} \sin(\omega t + \phi)^2\,dt} </math>

Using the trigonometric identity <math> \sin(\alpha + \beta) = \sin(\alpha) \cos(\beta) + \cos \alpha \sin \beta </math> we can re-write our equation;

where <math> \omega t + \phi = \alpha </math> and <math> \omega \tau = \beta </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)\,dt}{ \int_{-\infty}^{\infty} \sin(\omega t + \phi)^2\,dt}</math>

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

The fraction can be expanded;

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

Then since <math>\cos(\omega \tau)</math> does not contain a t term it can be moved outside of the integral;

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

The alike terms cancel to give;

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

Finally <math>\frac{\int_{-\infty}^{\infty}\cos(\omega t + \phi) \sin(\omega \tau)\,dt}{\int_{-\infty}^{\infty}\sin(\omega t + \phi)^2\,dt} = 0</math>, this is because <math>sin</math> and <math>\frac{\cos}{\sin}</math> are odd functions. When odd functions are integrated between symmetrical boundaries the product is zero.

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

Analysis of the solid, liquid and harmonic osillator VACFs against position
[[File:VACFplots.png|800px|thumb|right|figure 13]]
The minima in the VACFs represent a change in the direction of the velocity, this caused by collisions with other atoms.

The liquid and solid VACFs are similar but have some important differences. The solid atoms correlation with their initial velocity decreases faster than the liquids. The minima for the solid is much earlier and also much more negative. This happens because the solid atoms are much closer to each other and there are a lot more of them, resulting in an interaction after less time. When the correlation becomes negative it represents and inversed correlation which is where the velocity is no longer correlated to it's initial velocity but rather external velocities of other atoms. At the minima the hanges direction after colliding fully with other atoms

cked structure, the correlation becomes more negative beause the atoms are tightly constrained and so the interactions are very strong with other atoms. The atoms in the liquid are free moving and so when they collide they can move to stabilise The decorrelated solid function shows some oscillation about <math>C(\tau)=0</math> this arises from the jostling with other atoms in the fixed lattice and so we see a sort of harmonic oscillator motion. The liquid does not display this behavior because they are free to move and so will be able to "rapidly destroy any oscillatory motion."[x] by the "diffusive motion".

===Task 19: Integrating using the trapezium rule to get D===

The plots of <math> C(\tau) vs time </math> of the three phases in a small and large system were plotted along with their running integral.
[[File:solid3200.png|800px|thumb|right|figure 14]] [[File:solidM.png|800px|thumb|right|figure 15]] [[File:liquid8000.png|800px|thumb|right|figure 16]] [[File:liquidM.png|800px|thumb|right|figure 17]] [[File:gas8000.png|800px|thumb|right|figure 18]] [[File:gasM.png|800px|thumb|right|figure 19]]

Rep:Mod:IOPJKL

2015-12-04T13:55:22Z

Mcd13: /* Evaluating the 1D harmonic oscillator VACF using the expression for the position as a function of time */

'''Liquid Simulations Computational Experiment'''

==Theory==

===Task 1: Calculations of a Harmonic Oscillator using both Classical and Velocity-Verlet solutions===
[[File:All3graphs.png|800px|thumb|right|''Figure 1.a) Positions calculated via the two solutions vs time b) Difference between the positions calculated via different methods against time c) Total energy of the harmonic oscillator calculated using Velocity-Verlet solutions vs time '']]

'''Position calculations from the two methods''': The position <math>x</math> at time <math>t</math> was determined using <math>x(t) = A cos(\omega t + \phi)</math> where; <math>A = 1, \omega = 1</math> and <math> \phi= 0 </math>. These were plotted against time(Figure 1.a).

'''Error analysis of the two positions''': The difference between the two positions was calculated and plotted against time(Figure 1.b). The error increases with time and position, this loss of precision could come from two sources. As the position gets further from the original position the <math>a(t)dt^2</math> part of the equation has less importance [x] also the errors could build up since the previous position is used in the preceding position calculation.

'''Total Energy of the Harmonic-Oscillator vs time''': Values for velocity <math>v</math> and position <math>x</math> were obtained from the Velocity-Verlet calculation, values for the force-constant <math>k</math> and mass <math>m</math> were both given as 1. Using ,math> E= \frac{1}{2}mv^2 + \frac{1}{2}kx^2'' the total energies were calculated and plotted against time(Figure 1.c).

===Task 2: Plot of maxima error positions vs time===
The coordinates of the 5 maximas in figure 2 were estimated and plotted on a new graph. This yielded a linear function where <math>x(t) = 0.0004x - 7*10^{-5}</math>. [[File:MaximaDiff.png|200px|thumb|right|''Figure 2. Linear plot of the maxima error positions(figure 1.b) against time '']]

===Task 3: Timesteps to minimise variation in total energy===

Multiple timesteps were used to calculate the total energy of the harmonic oscillator. The percentage difference was taken between the initial energy assumed to be correct and the energy of the first minima, these results were then plotted(figure 3). In order to keep the change in total energy below 1% a timestep equal to or below 0.2 should be used. The change in total energy of the simulation should be kept to a minimum as the total energy should not vary over time in the harmonic oscillator.[[File:%vstimestep.png|200px|thumb|right|''Figure 3. Plot of percentage energy difference and timestep'']]

===Task 4: Lennard-Jones interaction===

'''1. Calculation of distance at zero-potential, r = r0:'''
<math>\phi\left(r_0\right) = 4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right)</math>
where φ = potential energy r = the distance ε = well depth and σ = distance from zero potential distance
<math>\phi\left(r_0\right)= 4\epsilon \left( \sigma^{12} - \sigma^6 r_0^6 \right)</math>
:<math>= 4\epsilon \left( 1 - \frac{r_0^6}{\sigma^6} \right) = 0</math>
:<math>= 1 - \frac{r_0^6}{\sigma^6} = 0 </math>
:<math> \sigma^6 = r_0^6 </math>
:<math> \sigma = r_0 </math>

2. '''Force at r0:'''

:<math>\mathbf{F}_0 = - \frac{\mathrm{d}\phi\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_0} = - \frac{\mathrm{d}}{\mathrm{d}\mathbf{r}_0} 4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right) = 4\epsilon \left( -12\frac{\sigma^{12}}{r_0^{13}} + 6\frac{\sigma^6}{r_0^7} \right)</math>
Since <math>\sigma = r_0 </math>
:<math>= 4\epsilon\left( -12\frac{\sigma^{12}}{\sigma^{13}} + 6\frac{\sigma^6}{\sigma^7} \right) = 4\epsilon\frac{-6}{\sigma} = -\frac{24\epsilon}{\sigma}</math>

'''3. Finding the equilibrium separation ''req:'' '''

:<math>\mathbf{F}_{eq} = - \frac{\mathrm{d}\phi\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_{eq}} = - \frac{\mathrm{d}}{\mathrm{d}\mathbf{r}_{eq}} 4\epsilon \left( \frac{\sigma^{12}}{r_{eq}^{12}} - \frac{\sigma^6}{r_{eq}^6} \right) = 4\epsilon \left( -12\frac{\sigma^{12}}{r_{eq}^{13}} + 6\frac{\sigma^6}{r_{eq}^7} \right) = -48\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) + 24\epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right) = 0</math>
:<math>48\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) = 24\epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right)</math>
:<math>2\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) = \epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right)</math>
:<math>2\epsilon \left(\sigma^{12}\right) = \epsilon \left(\sigma^6 r_{eq}^6 \right)</math>
:<math>2\epsilon \left(\sigma^{6}\right) = \epsilon \left(\ r_{eq}^6 \right)</math>
:<math>2\left(\sigma^{6}\right) = \left(\ r_{eq}^6 \right)</math>
:<math>\sqrt[6]{2}\sigma = r_{eq}</math>

'''4. Finding the well-depth at equilibrium:'''

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

'''5. Evaluation of integrals:'''

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

:<math>\int_{2}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{2}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{2}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>= [4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_2^{\infty}</math>
:<math>= 4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(2)^{11}} - \frac{1}{-5(2)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(2)^{11}} - \frac{1}{-5(2)^5} \right) = -24.82 KJ/mol </math>

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

:<math>\int_{2.5}^\infty \phi\left(r\right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>= [4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_{2.5}^{\infty}</math>
:<math>= 4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(2.5)^{11}} - \frac{1}{-5(2.5)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(2.5)^{11}} - \frac{1}{-5(2.5)^5} \right) = -8.18 KJ/mol </math>

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

:<math>\int_{3}^\infty \phi\left(r\right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>[4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_{3}^{\infty}</math>
:<math>4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(3)^{11}} - \frac{1}{-5(3)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(3)^{11}} - \frac{1}{-5(3)^5} \right) = -3.29 KJ/mol </math>

===Task 5: Estimating the number of molecules in volumes of H2O===

1. '''Number of molecules in 1mL:'''
:<math> Mr(H^2O) = 18.02g/mol,\rho = 1g/cm^3 </math>
<math>n = \frac{m}{Mr} = \frac{1}{18.02} = 5.55x10^{-2}</math>
<math>Number of molecules = 5.55x10^{-2} x 6.022x10^{23} = 3.34192x10^{22}</math>

2. '''Volume of 1000 molecules:'''

<math>n = \frac{1000}{6.022x10^{23}} = 1.66x10^{-21} mol</math>
<math>m = 1.66x10^{-21} x 18.02 = 2.99x10^{-20}g</math>
<math>V = \frac{2.99x10^{-20}}{1} = 2.99x10^{-20}cm^3</math>

===Task 6: Periodic boundary conditions===

When the atom at <math>(0.5,0.5,0.5)</math> in the cubic simulation box(<math>(0,0,0,) to (1,1,1)</math>) moves along the vector <math>(0.7,0.6,0.2)</math> it ends up at <math>(0.2,0.1,0.7)</math> (when the periodic boundry condition is taken into consideration).

===Task 7: Converting reduced values to real values:===

1. '''Reduced distance'''

:<math> r* = 3.2 = \frac {r}{\sigma} </math>
:<math> r = 3.2 x 0.34x10^{-9} = 1.09x10^-9 = 1.09nm </math>

2. '''Well-depth'''

:<math> \frac {\epsilon}{K_B} = 120K </math>
:<math>\epsilon = 120K_B = \frac {8.69x10^{24}}{6.022x10^{23}} = 997.71 J/mol = 9.98x10^{-1}</math>

3. '''Temperature'''

:<math> T^* = 1.5 = \frac {K_BT}{\epsilon} </math>
:<math>T = \frac {T^* \epsilon}{K_B} = \frac {1.5 x 1.657x10^{-21}}{K_B} = 180K</math>

==Equilibration==

===Task 8: Problems with random starting coordinates in non-solid simulations===

If atoms were given random starting coordinates then they could end up within a contact proximity, they could bump into each other and then they would not only experience the potential of each but also the kinetic energy.

===Task 9: The number density of lattice points calculations:===

The number density ''n'' is given by the following equation
<math> n = \frac {N}{V}</math> where N is the number of atoms/lattice points per unit cell and V is volumes of the unit cell.

1. The number density of a simple cubic lattice with inter lattice point distance of 1.07722

For a simple cubic lattice there is one atom/lattice point, using equation x
:<math> n = \frac {N}{V} = \frac {N}{r^{*3}} = \frac {1}{1.07722^3} = 0.8 </math>

2. '''Side length of a face-centred cubic lattice with a number density of 1.2'''

A face-centred cubic lattice has 4 atoms/lattice points.
<math> r^* = \sqrt[3]{V} = \sqrt[3]{\frac {N}{n}} = \sqrt[3]{\frac {4}{1.2}} = 1.49380 </math>

===Task 10: How many atoms are created by create_atoms 1 box command for a face-centred cubic lattice?===
Since there are 4 atoms per lattice and the box is ten lattice spacings larger than the single lattice spacing then there would be 4000 atoms.

===Task 11: LAMMPS commands for setting the properties of atoms:===

1. Mass command

mass ''I'' ''value'' = mass 1 1.0
This command is used for setting the mass of the atom, in the general command I is the atom type and value is the mass(in reduced units). In the example the atom is type 1 and the mass is 1.0.

2. Pair_style command

pair_style style args = pair_style lj/cut 3.0

The pair_style command is used to describe paired interactions, style is the interaction which is being looked at and args is the argument used for that style. In the example above the style is defining the cutoff distance, that is the distance at which interactions become negligible. The argument for the cutoff distance is 3.0, after which the simulation will not compute.

3. pair_coeff command

pair_coeff ''* * args'' = pair_coeff * * 1.0 1.0
This command is used to define the force field coefficients between pairs of atoms, the asterix terms are to define the coefficients for all atoms, the second part of this command, args, defines the coeffecients.

===Task 12: Specifying the velocity of each atom===

The Velocity-Verlet algorithm is used to specify velocity and position
<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 \ \ (9)</math>

===Task 13: Running the simulation:===
The command $(100/timestep) command allows for the immediate evaluation of the the formula with the pre-defined variable i.e. the timestep. The following line variable n_steps equal floor(100/$({timestep}) will be broken down into individual contributions. The variable component is used to allows for the return to re-run a command, this allows creates a loop. The n_steps part describes the number of timesteps you wish to calculate which is defined in the 'Run Simulation' command. Equal is used to create a mathematical argument. Floor(11/$({timestep}) is useful in setting the maximum timestep bigger than the original one. In the run simulation part first we define the number of timesteps to be considered and then report the floor value i.e. the total time it takes. These lines are essential in creating a loop where each new timestep may be calculated subsequently.
Make it easier to change the timestep, all values dependedent on the timestep will be changed with it and so it reduces possible human errors when changing all variables and also makes it easier and faster.

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

===Task 14: Checking equilibration===

1. Pressure, Temperature and energy against time for simulation ran at 0.001 timestep

PIC fig 4

2. Equilibration

The simulation reaches equilibration at approximately 0.28, this is the timestep at which the simulation becomes relatively constant.

3. Maximising results accuracy whilst keeping the timestep within a reasonable timescale
The best timestep for maximum accuracy and sufficient information is 0.0025. It overlaps with 0.001 and so although it is more than twice the time as 0.001 it maintains a satisfactory level of accuracy. The others give equilibrium energies which are too far off the most accurate value, with the worst being for 0.01 as it diverges quite dramatically.

[[File:Energyplotsvstime.png|400px|thumb|right|''Figure 5. The energy vs time plots of simulations with varying timesteps '']]

==Running simulations under specific conditions==

===Task 15: Finding the expression for factor <math> \gamma </math>===

<math> \gamma</math> is a factor by which the velocity can be corrected.

<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT </math> (1)
<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT' </math> (2)

Where ''T''' is the target temperature and ''T'' is the instantaneous temperature

<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT </math>
<math> \frac {1}{2} m_iv_i^2 = \frac {1}{2}NK_BT </math>
<math> m_iv_i^2 = NK_BT </math>
<math> v_i^2 = \frac {NK_BT}{m_i} </math>
<math> v_i = \sqrt[2]{\frac {NK_BT}{m_i}} </math>

Insert v_i into equation x

<math> \frac {1}{2} \sum m_i \gamma^2 \frac {NK_BT}{m_i} = \frac {3}{2}NK_BT' </math>
<math> \frac {1}{2} m_i \gamma^2 \frac {NK_BT}{m_i} = \frac {1}{2}NK_BT' </math>
<math> \gamma^2NK_BT = NK_BT' </math>
<math> \gamma^2 = \frac {NK_BT'}{NK_BT} </math>
<math> \gamma^2 = \frac {T'}{T} </math>
<math> \gamma = \sqrt[2]{\frac {T'}{T}} </math>

===Task 16: The importance of number, repetition and frequency of average calculations with timesteps===

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

100 sets the number of timesteps after which another thermodynamic value is used to calculate the average, 1000 is the number of thermodynamic values which are used for the average and 1000000 sets the frequency at which averages are calculated.

100000 timesteps will be run.

===Task 17: Equation of state plots===

5 temperatures above the critical temeperature were chosen(1.6, 1.8, 2.0, 3.0 and 5.0) these were all simulated with a pressure of 2.6(''n-1'') and also 3.3(''n-2''). The timestep chose was 0.01 as in previous experiments it was the largest timesteps yielding the most accuracy.
The pressures were plotted against the temperatures along with their error bars.

The ideal gas densities were calculated with the following equation:
<math> \rho = \frac {N}{V} = \frac {P}{K_BT} = \frac {P}{T} </math> where K_B = 1 in reduced terms
[[File:DensityvsT.png|400px|thumb|right|''Figure 6. Densities n-1 and n-2 calculated via simulation and Ideal Gas methods plotted against temperature'']]

The simulated density is lower than that of the ideal gas method, this can be explained by an assumption to satisy an ideal gas. The atoms must not have any forces upon them and also they mustn't take up any space. In the simulation the atoms do take up space and therefore they have a lower density. This discrepancy increases with pressure. The general trend of the two plots is a decrease in heat capacity as temperature is increased.

==Heat Capacity Calculation==

===Task 18: The temperature dependence of the isochoric heat capacity===
[[File:CvVcommand.png|400px|thumb|left|''Figure 7. Cv/V as a function of T* at P= 0.2 and 0.8'']]

A LAAMPS script which could calculate the isochoric heat capacity from calculated thermodynamic values was used figure 7. The Cv/V calculation was done using equation x.

[[File:CvVvsT.png|400px|thumb|right|''Figure 8. Cv/V as a function of T* at P= 0.2 and 0.8'']]
The heat capacity represents the ability to transfer heat for a given temperature or how much energy will effect a given temperature.

The general decrease of Cv/V vs T is what would be expected. At a higher temperature the energy levels are closer together and at a lower temperature they are further apart. A heat increment will have a greater effect on the higher temperature system, it will cause a greater temperature change since more energy levels are accessible, resulting in a low heat capacity. At a lower temperature the same heat increment that caused a large temperature change at higher temperature will not have the same effect, less energy levels are accessible and so the heat capacity will be higher.

The difference between the two pressured systems in ''figure 7'' is what would be expected, the higher pressure plot has a higher isochoric heat capacity. At a higher pressure the atoms are closer together, they can thus transfer heat more efficiently with a given temperature change than a system with a lower temperature since the atoms are further away.

==Structural properties and the radial distribution function==
===Task 18: The Radial Distribution Functions of a solid, liquid and gas===

'''1. Comparisons of the solid, liquid and gas radial distributions functions as well as the integrated function '''
[[File:rdfslg.png|400px|thumb|left|figure 9]]
[[File:rdfint.png|400px|thumb|left|figure 10]]

The distribution plot(figure 9) of the gas has one short, broad peak which is close to the reference atom; as you move further away the distribution becomes constant. This characteristic shows the small probability of finding an atom close to the reference atom; also when it becomes constant this shows how there is negligible probability of finding another atom at the distance studied. When looking at the integrated function(figure 10), it emphasizes how few atoms you will find in the distance, this is characteristic of a gas' dispersed atoms.

The liquid's distribution has a taller first peak than the gas' and then multiple peaks tapering off as you move further away; it is evident that there is a higher probability of finding an atom throughout the space than in a gas. Additionally the integrated distribution(figure 10) shows a larger number of atoms as you move through the space. This demonstrates the higher proximity of the atoms in the liquid state.

Finally the solid's distribution has a very sharp and tall peak close to the reference atom but then as you move further away the peaks become smaller but continuously fluctuate. The sharpness of the peaks and their continuous distribution arise from the highly ordered structure of the solid. The peaks are sharp unlike the gas and liquid because all of the atoms are in fixed positions and so at a certain position there will be a number of atoms in an equivalent position around the sphere of the distribution function. Futhermore the solid's first peak is considerably larger than the other phases, the probibility is much greater sdue to the closeness of the atoms. The integrated function shows a very large number of atoms in the space, this again illustrates the closeness of the atoms in the solid state.

'''2. Using the solid RDF to extract information on the lattice structure'''

[[File:rdf3d.png|200px|thumb|right|figure 11]]

The three first peaks of the solid's radial distribution function are distances 1.0, 1.5 and 1.8. The lattice spacing is 1.5 and the calculated distances ra, rb and rc(where r is the distance from the reference atom to the atoms a, b and c) were found to be 1.0, 1.5 and 1.8(respectively). Therefore the first three peaks 1, 2 and 3 correspond to the lattice points a, b and c(respectively).

The coordination of atom a, b and c is 12, 6 and 24(respectively). These values were interpolated from the magnification of the integrated solid distribution function.

==Dynamical properties and the diffusion coefficient==

===Task 19: The mean squared displacement(MSD) of large and small simulations===

[[File:msdslv.png|800px|thumb|right|figure 12]]
<math> MSD(t) = <((t_0 + t)-(t_0))>^2 </math> where; <math>t_0 =</math> the initial timestep and <math>t =</math>the new timestep Equation x

The mean squared displacement(equation x) is a measure of how much random motion there is within the system. In this simulation the ability to move randomly is determined by the amount of free space and therefore the density. In the gaseous phase there is a lot of free space, the atoms move freely and encounter only a small number of atoms to collide with thus impeding it's motion; the linear increase of the plot reflects this behaviour. The liquid has a lower degree of random motion than the gas, due closer proximity of the atoms, forces such as drag will now effect the atom's motion as it moves past other atoms. The two first plots show that the atom is able to move in drift, however the liquid experiences more drag and has a smaller gradient. The solid displays the least amount of translational freedom, the close packed structure prevents any random translational motion; it does however show some random movement, arising from vibrations of the constrained atoms. Therefore the plots in figure 12 are a reasonable representation of the random movement of the three phases.

Finding the diffusion coefficient <math>D</math>

<math>D = \frac{1}{6} \frac{\delta <r^2(t)>}{\delta t}> </math> where <math>D =</math> the diffusion coefficient in <math>m s^{-2}</math>,<math><r^2(t) =</math> MSD and <math> t =</math> timestep
The gradient of the slopes give the diffusion coefficent, however this uses a dimensionless timestep denominator, therefore we must convert it to time units by dividing by the timestep(0.002).

The gradient was found by fitting a linear trendline through the function.

Solid

Conversion of the gradient to give time units gave <math>D = 4.17 x 10^{-6} ms^{-2}</math> for the million atom simulation and <math>D = 5.83E-08 x 10^{-7} ms^{-2}</math> for the 3200 atom simulation.

Liquid

Conversion of the gradient to give time units gave <math>D = 8.33 x 10^{-2} ms^{-2}</math> for the million atom simulation and <math>D = 8.33 x 10^{-2} ms^{-2}</math> for the 3200 atom simulation.

Gas

A trendline was fitted to the function, since the simulation took some time to produce the linear behavior it was fitted to be linear with the latter part of the function. Conversion of the gradient to give time units gave <math>D = 3.09 ms^{-2}</math> for the million atom simulation and <math>D = 2.43 ms^{-2}</math> for the 3200 atom simulation.

===Task 20: The velocity auto-correlation function(VACF)===

==Evaluating the 1D harmonic oscillator VACF using the expression for the position as a function of time==

The normalized VACF is given by;

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

We must simplify this expression in order for it to be useful, to do this we must first find expressions for <math>v(t)</math> and <math>v(t + \tau)</math>:

The expression for position as a function of time for the HO is defined as;
<math>x(t) = A\cos(\omega t + \phi)</math>

This can be substituted into the expression of velocity as a function of time;

<math>v(t) = \frac{dx(t)}{dt} = \frac{dA\cos(\omega t + \phi)}{dt}</math>

This is then differentiated using the chain rule;

<math> \frac{dx(t)}{dt} = \frac{dx(t)}{du} x \frac{du}{dt} </math>

Where <math>u = \omega t + \phi</math>

<math>v(t) = -\omega A\sin(\omega t + \phi) </math>

We also need our velocity as a function of time after <math>\tau</math>, this is acheived by substituting in <math>t = t + \tau</math>;

<math> v(t + \tau) = -\omega A\sin(\omega (t + \tau) + \phi)</math>
<math> = -\omega A\sin(\omega t + \omega\tau + \phi) </math>

The two velocity expressions are substituted into the VACF;

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

::::<math> = \frac{\int_{-\infty}^{\infty} \sin(\omega t + \phi) \sin(\omega t + \omega\tau + \phi)\,dt}{\int_{-\infty}^{\infty} \sin(\omega t + \phi)^2\,dt} </math>

Using the trigonometric identity <math> \sin(\alpha + \beta) = \sin(\alpha) \cos(\beta) + \cos \alpha \sin \beta </math> we can re-write our equation;

where <math> \omega t + \phi = \alpha </math> and <math> \omega \tau = \beta </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)\,dt}{ \int_{-\infty}^{\infty} \sin(\omega t + \phi)^2\,dt}</math>

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

The fraction can be expanded;

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

Then since <math>\cos(\omega \tau)</math> does not contain a t term it can be moved outside of the integral;

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

The alike terms cancel to give;

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

Finally <math>\frac{\int_{-\infty}^{\infty}\cos(\omega t + \phi) \sin(\omega \tau)\,dt}{\int_{-\infty}^{\infty}\sin(\omega t + \phi)^2\,dt} = 0</math>, this is because <math>sin</math> and <math>\frac{\cos}{\sin}</math> are odd functions. When odd functions are integrated between symmetrical boundaries the product is zero.

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

Analysis of the solid, liquid and harmonic osillator VACFs against position
[[File:VACFplots.png|800px|thumb|right|figure 13]]
The minima in the VACFs represent a change in the direction of the velocity, this caused by collisions with other atoms.

The liquid and solid VACFs are similar but have some important differences. The solid atoms correlation with their initial velocity decreases faster than the liquids. The minima for the solid is much earlier and also much more negative. This happens because the solid atoms are much closer to each other and there are a lot more of them, resulting in an interaction after less time. When the correlation becomes negative it represents and inversed correlation i.e. the velocity is no longer correlated to it's initial velocity but rather external velocities of other atoms.

cked structure, the correlation becomes more negative beause the atoms are tightly constrained and so the interactions are very strong with other atoms. The atoms in the liquid are free moving and so when they collide they can move to stabilise The decorrelated solid function shows some oscillation about <math>C(\tau)=0</math> this arises from the jostling with other atoms in the fixed lattice and so we see a sort of harmonic oscillator motion. The liquid does not display this behavior because they are free to move and so will be able to "rapidly destroy any oscillatory motion."[x] by the "diffusive motion".

===Task 19: Integrating using the trapezium rule to get D===

The plots of <math> C(\tau) vs time </math> of the three phases in a small and large system were plotted along with their running integral.
[[File:solid3200.png|800px|thumb|right|figure 14]] [[File:solidM.png|800px|thumb|right|figure 15]] [[File:liquid8000.png|800px|thumb|right|figure 16]] [[File:liquidM.png|800px|thumb|right|figure 17]] [[File:gas8000.png|800px|thumb|right|figure 18]] [[File:gasM.png|800px|thumb|right|figure 19]]

Rep:Mod:IOPJKL

2015-12-04T13:21:18Z

Mcd13:

'''Liquid Simulations Computational Experiment'''

==Theory==

===Task 1: Calculations of a Harmonic Oscillator using both Classical and Velocity-Verlet solutions===
[[File:All3graphs.png|800px|thumb|right|''Figure 1.a) Positions calculated via the two solutions vs time b) Difference between the positions calculated via different methods against time c) Total energy of the harmonic oscillator calculated using Velocity-Verlet solutions vs time '']]

'''Position calculations from the two methods''': The position <math>x</math> at time <math>t</math> was determined using <math>x(t) = A cos(\omega t + \phi)</math> where; <math>A = 1, \omega = 1</math> and <math> \phi= 0 </math>. These were plotted against time(Figure 1.a).

'''Error analysis of the two positions''': The difference between the two positions was calculated and plotted against time(Figure 1.b). The error increases with time and position, this loss of precision could come from two sources. As the position gets further from the original position the <math>a(t)dt^2</math> part of the equation has less importance [x] also the errors could build up since the previous position is used in the preceding position calculation.

'''Total Energy of the Harmonic-Oscillator vs time''': Values for velocity <math>v</math> and position <math>x</math> were obtained from the Velocity-Verlet calculation, values for the force-constant <math>k</math> and mass <math>m</math> were both given as 1. Using ,math> E= \frac{1}{2}mv^2 + \frac{1}{2}kx^2'' the total energies were calculated and plotted against time(Figure 1.c).

===Task 2: Plot of maxima error positions vs time===
The coordinates of the 5 maximas in figure 2 were estimated and plotted on a new graph. This yielded a linear function where <math>x(t) = 0.0004x - 7*10^{-5}</math>. [[File:MaximaDiff.png|200px|thumb|right|''Figure 2. Linear plot of the maxima error positions(figure 1.b) against time '']]

===Task 3: Timesteps to minimise variation in total energy===

Multiple timesteps were used to calculate the total energy of the harmonic oscillator. The percentage difference was taken between the initial energy assumed to be correct and the energy of the first minima, these results were then plotted(figure 3). In order to keep the change in total energy below 1% a timestep equal to or below 0.2 should be used. The change in total energy of the simulation should be kept to a minimum as the total energy should not vary over time in the harmonic oscillator.[[File:%vstimestep.png|200px|thumb|right|''Figure 3. Plot of percentage energy difference and timestep'']]

===Task 4: Lennard-Jones interaction===

'''1. Calculation of distance at zero-potential, r = r0:'''
<math>\phi\left(r_0\right) = 4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right)</math>
where φ = potential energy r = the distance ε = well depth and σ = distance from zero potential distance
<math>\phi\left(r_0\right)= 4\epsilon \left( \sigma^{12} - \sigma^6 r_0^6 \right)</math>
:<math>= 4\epsilon \left( 1 - \frac{r_0^6}{\sigma^6} \right) = 0</math>
:<math>= 1 - \frac{r_0^6}{\sigma^6} = 0 </math>
:<math> \sigma^6 = r_0^6 </math>
:<math> \sigma = r_0 </math>

2. '''Force at r0:'''

:<math>\mathbf{F}_0 = - \frac{\mathrm{d}\phi\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_0} = - \frac{\mathrm{d}}{\mathrm{d}\mathbf{r}_0} 4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right) = 4\epsilon \left( -12\frac{\sigma^{12}}{r_0^{13}} + 6\frac{\sigma^6}{r_0^7} \right)</math>
Since <math>\sigma = r_0 </math>
:<math>= 4\epsilon\left( -12\frac{\sigma^{12}}{\sigma^{13}} + 6\frac{\sigma^6}{\sigma^7} \right) = 4\epsilon\frac{-6}{\sigma} = -\frac{24\epsilon}{\sigma}</math>

'''3. Finding the equilibrium separation ''req:'' '''

:<math>\mathbf{F}_{eq} = - \frac{\mathrm{d}\phi\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_{eq}} = - \frac{\mathrm{d}}{\mathrm{d}\mathbf{r}_{eq}} 4\epsilon \left( \frac{\sigma^{12}}{r_{eq}^{12}} - \frac{\sigma^6}{r_{eq}^6} \right) = 4\epsilon \left( -12\frac{\sigma^{12}}{r_{eq}^{13}} + 6\frac{\sigma^6}{r_{eq}^7} \right) = -48\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) + 24\epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right) = 0</math>
:<math>48\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) = 24\epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right)</math>
:<math>2\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) = \epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right)</math>
:<math>2\epsilon \left(\sigma^{12}\right) = \epsilon \left(\sigma^6 r_{eq}^6 \right)</math>
:<math>2\epsilon \left(\sigma^{6}\right) = \epsilon \left(\ r_{eq}^6 \right)</math>
:<math>2\left(\sigma^{6}\right) = \left(\ r_{eq}^6 \right)</math>
:<math>\sqrt[6]{2}\sigma = r_{eq}</math>

'''4. Finding the well-depth at equilibrium:'''

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

'''5. Evaluation of integrals:'''

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

:<math>\int_{2}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{2}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{2}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>= [4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_2^{\infty}</math>
:<math>= 4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(2)^{11}} - \frac{1}{-5(2)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(2)^{11}} - \frac{1}{-5(2)^5} \right) = -24.82 KJ/mol </math>

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

:<math>\int_{2.5}^\infty \phi\left(r\right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>= [4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_{2.5}^{\infty}</math>
:<math>= 4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(2.5)^{11}} - \frac{1}{-5(2.5)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(2.5)^{11}} - \frac{1}{-5(2.5)^5} \right) = -8.18 KJ/mol </math>

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

:<math>\int_{3}^\infty \phi\left(r\right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>[4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_{3}^{\infty}</math>
:<math>4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(3)^{11}} - \frac{1}{-5(3)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(3)^{11}} - \frac{1}{-5(3)^5} \right) = -3.29 KJ/mol </math>

===Task 5: Estimating the number of molecules in volumes of H2O===

1. '''Number of molecules in 1mL:'''
:<math> Mr(H^2O) = 18.02g/mol,\rho = 1g/cm^3 </math>
<math>n = \frac{m}{Mr} = \frac{1}{18.02} = 5.55x10^{-2}</math>
<math>Number of molecules = 5.55x10^{-2} x 6.022x10^{23} = 3.34192x10^{22}</math>

2. '''Volume of 1000 molecules:'''

<math>n = \frac{1000}{6.022x10^{23}} = 1.66x10^{-21} mol</math>
<math>m = 1.66x10^{-21} x 18.02 = 2.99x10^{-20}g</math>
<math>V = \frac{2.99x10^{-20}}{1} = 2.99x10^{-20}cm^3</math>

===Task 6: Periodic boundary conditions===

When the atom at <math>(0.5,0.5,0.5)</math> in the cubic simulation box(<math>(0,0,0,) to (1,1,1)</math>) moves along the vector <math>(0.7,0.6,0.2)</math> it ends up at <math>(0.2,0.1,0.7)</math> (when the periodic boundry condition is taken into consideration).

===Task 7: Converting reduced values to real values:===

1. '''Reduced distance'''

:<math> r* = 3.2 = \frac {r}{\sigma} </math>
:<math> r = 3.2 x 0.34x10^{-9} = 1.09x10^-9 = 1.09nm </math>

2. '''Well-depth'''

:<math> \frac {\epsilon}{K_B} = 120K </math>
:<math>\epsilon = 120K_B = \frac {8.69x10^{24}}{6.022x10^{23}} = 997.71 J/mol = 9.98x10^{-1}</math>

3. '''Temperature'''

:<math> T^* = 1.5 = \frac {K_BT}{\epsilon} </math>
:<math>T = \frac {T^* \epsilon}{K_B} = \frac {1.5 x 1.657x10^{-21}}{K_B} = 180K</math>

==Equilibration==

===Task 8: Problems with random starting coordinates in non-solid simulations===

If atoms were given random starting coordinates then they could end up within a contact proximity, they could bump into each other and then they would not only experience the potential of each but also the kinetic energy.

===Task 9: The number density of lattice points calculations:===

The number density ''n'' is given by the following equation
<math> n = \frac {N}{V}</math> where N is the number of atoms/lattice points per unit cell and V is volumes of the unit cell.

1. The number density of a simple cubic lattice with inter lattice point distance of 1.07722

For a simple cubic lattice there is one atom/lattice point, using equation x
:<math> n = \frac {N}{V} = \frac {N}{r^{*3}} = \frac {1}{1.07722^3} = 0.8 </math>

2. '''Side length of a face-centred cubic lattice with a number density of 1.2'''

A face-centred cubic lattice has 4 atoms/lattice points.
<math> r^* = \sqrt[3]{V} = \sqrt[3]{\frac {N}{n}} = \sqrt[3]{\frac {4}{1.2}} = 1.49380 </math>

===Task 10: How many atoms are created by create_atoms 1 box command for a face-centred cubic lattice?===
Since there are 4 atoms per lattice and the box is ten lattice spacings larger than the single lattice spacing then there would be 4000 atoms.

===Task 11: LAMMPS commands for setting the properties of atoms:===

1. Mass command

mass ''I'' ''value'' = mass 1 1.0
This command is used for setting the mass of the atom, in the general command I is the atom type and value is the mass(in reduced units). In the example the atom is type 1 and the mass is 1.0.

2. Pair_style command

pair_style style args = pair_style lj/cut 3.0

The pair_style command is used to describe paired interactions, style is the interaction which is being looked at and args is the argument used for that style. In the example above the style is defining the cutoff distance, that is the distance at which interactions become negligible. The argument for the cutoff distance is 3.0, after which the simulation will not compute.

3. pair_coeff command

pair_coeff ''* * args'' = pair_coeff * * 1.0 1.0
This command is used to define the force field coefficients between pairs of atoms, the asterix terms are to define the coefficients for all atoms, the second part of this command, args, defines the coeffecients.

===Task 12: Specifying the velocity of each atom===

The Velocity-Verlet algorithm is used to specify velocity and position
<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 \ \ (9)</math>

===Task 13: Running the simulation:===
The command $(100/timestep) command allows for the immediate evaluation of the the formula with the pre-defined variable i.e. the timestep. The following line variable n_steps equal floor(100/$({timestep}) will be broken down into individual contributions. The variable component is used to allows for the return to re-run a command, this allows creates a loop. The n_steps part describes the number of timesteps you wish to calculate which is defined in the 'Run Simulation' command. Equal is used to create a mathematical argument. Floor(11/$({timestep}) is useful in setting the maximum timestep bigger than the original one. In the run simulation part first we define the number of timesteps to be considered and then report the floor value i.e. the total time it takes. These lines are essential in creating a loop where each new timestep may be calculated subsequently.
Make it easier to change the timestep, all values dependedent on the timestep will be changed with it and so it reduces possible human errors when changing all variables and also makes it easier and faster.

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

===Task 14: Checking equilibration===

1. Pressure, Temperature and energy against time for simulation ran at 0.001 timestep

PIC fig 4

2. Equilibration

The simulation reaches equilibration at approximately 0.28, this is the timestep at which the simulation becomes relatively constant.

3. Maximising results accuracy whilst keeping the timestep within a reasonable timescale
The best timestep for maximum accuracy and sufficient information is 0.0025. It overlaps with 0.001 and so although it is more than twice the time as 0.001 it maintains a satisfactory level of accuracy. The others give equilibrium energies which are too far off the most accurate value, with the worst being for 0.01 as it diverges quite dramatically.

[[File:Energyplotsvstime.png|400px|thumb|right|''Figure 5. The energy vs time plots of simulations with varying timesteps '']]

==Running simulations under specific conditions==

===Task 15: Finding the expression for factor <math> \gamma </math>===

<math> \gamma</math> is a factor by which the velocity can be corrected.

<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT </math> (1)
<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT' </math> (2)

Where ''T''' is the target temperature and ''T'' is the instantaneous temperature

<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT </math>
<math> \frac {1}{2} m_iv_i^2 = \frac {1}{2}NK_BT </math>
<math> m_iv_i^2 = NK_BT </math>
<math> v_i^2 = \frac {NK_BT}{m_i} </math>
<math> v_i = \sqrt[2]{\frac {NK_BT}{m_i}} </math>

Insert v_i into equation x

<math> \frac {1}{2} \sum m_i \gamma^2 \frac {NK_BT}{m_i} = \frac {3}{2}NK_BT' </math>
<math> \frac {1}{2} m_i \gamma^2 \frac {NK_BT}{m_i} = \frac {1}{2}NK_BT' </math>
<math> \gamma^2NK_BT = NK_BT' </math>
<math> \gamma^2 = \frac {NK_BT'}{NK_BT} </math>
<math> \gamma^2 = \frac {T'}{T} </math>
<math> \gamma = \sqrt[2]{\frac {T'}{T}} </math>

===Task 16: The importance of number, repetition and frequency of average calculations with timesteps===

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

100 sets the number of timesteps after which another thermodynamic value is used to calculate the average, 1000 is the number of thermodynamic values which are used for the average and 1000000 sets the frequency at which averages are calculated.

100000 timesteps will be run.

===Task 17: Equation of state plots===

5 temperatures above the critical temeperature were chosen(1.6, 1.8, 2.0, 3.0 and 5.0) these were all simulated with a pressure of 2.6(''n-1'') and also 3.3(''n-2''). The timestep chose was 0.01 as in previous experiments it was the largest timesteps yielding the most accuracy.
The pressures were plotted against the temperatures along with their error bars.

The ideal gas densities were calculated with the following equation:
<math> \rho = \frac {N}{V} = \frac {P}{K_BT} = \frac {P}{T} </math> where K_B = 1 in reduced terms
[[File:DensityvsT.png|400px|thumb|right|''Figure 6. Densities n-1 and n-2 calculated via simulation and Ideal Gas methods plotted against temperature'']]

The simulated density is lower than that of the ideal gas method, this can be explained by an assumption to satisy an ideal gas. The atoms must not have any forces upon them and also they mustn't take up any space. In the simulation the atoms do take up space and therefore they have a lower density. This discrepancy increases with pressure. The general trend of the two plots is a decrease in heat capacity as temperature is increased.

==Heat Capacity Calculation==

===Task 18: The temperature dependence of the isochoric heat capacity===
[[File:CvVcommand.png|400px|thumb|left|''Figure 7. Cv/V as a function of T* at P= 0.2 and 0.8'']]

A LAAMPS script which could calculate the isochoric heat capacity from calculated thermodynamic values was used figure 7. The Cv/V calculation was done using equation x.

[[File:CvVvsT.png|400px|thumb|right|''Figure 8. Cv/V as a function of T* at P= 0.2 and 0.8'']]
The heat capacity represents the ability to transfer heat for a given temperature or how much energy will effect a given temperature.

The general decrease of Cv/V vs T is what would be expected. At a higher temperature the energy levels are closer together and at a lower temperature they are further apart. A heat increment will have a greater effect on the higher temperature system, it will cause a greater temperature change since more energy levels are accessible, resulting in a low heat capacity. At a lower temperature the same heat increment that caused a large temperature change at higher temperature will not have the same effect, less energy levels are accessible and so the heat capacity will be higher.

The difference between the two pressured systems in ''figure 7'' is what would be expected, the higher pressure plot has a higher isochoric heat capacity. At a higher pressure the atoms are closer together, they can thus transfer heat more efficiently with a given temperature change than a system with a lower temperature since the atoms are further away.

==Structural properties and the radial distribution function==
===Task 18: The Radial Distribution Functions of a solid, liquid and gas===

'''1. Comparisons of the solid, liquid and gas radial distributions functions as well as the integrated function '''
[[File:rdfslg.png|400px|thumb|left|figure 9]]
[[File:rdfint.png|400px|thumb|left|figure 10]]

The distribution plot(figure 9) of the gas has one short, broad peak which is close to the reference atom; as you move further away the distribution becomes constant. This characteristic shows the small probability of finding an atom close to the reference atom; also when it becomes constant this shows how there is negligible probability of finding another atom at the distance studied. When looking at the integrated function(figure 10), it emphasizes how few atoms you will find in the distance, this is characteristic of a gas' dispersed atoms.

The liquid's distribution has a taller first peak than the gas' and then multiple peaks tapering off as you move further away; it is evident that there is a higher probability of finding an atom throughout the space than in a gas. Additionally the integrated distribution(figure 10) shows a larger number of atoms as you move through the space. This demonstrates the higher proximity of the atoms in the liquid state.

Finally the solid's distribution has a very sharp and tall peak close to the reference atom but then as you move further away the peaks become smaller but continuously fluctuate. The sharpness of the peaks and their continuous distribution arise from the highly ordered structure of the solid. The peaks are sharp unlike the gas and liquid because all of the atoms are in fixed positions and so at a certain position there will be a number of atoms in an equivalent position around the sphere of the distribution function. Futhermore the solid's first peak is considerably larger than the other phases, the probibility is much greater sdue to the closeness of the atoms. The integrated function shows a very large number of atoms in the space, this again illustrates the closeness of the atoms in the solid state.

'''2. Using the solid RDF to extract information on the lattice structure'''

[[File:rdf3d.png|200px|thumb|right|figure 11]]

The three first peaks of the solid's radial distribution function are distances 1.0, 1.5 and 1.8. The lattice spacing is 1.5 and the calculated distances ra, rb and rc(where r is the distance from the reference atom to the atoms a, b and c) were found to be 1.0, 1.5 and 1.8(respectively). Therefore the first three peaks 1, 2 and 3 correspond to the lattice points a, b and c(respectively).

The coordination of atom a, b and c is 12, 6 and 24(respectively). These values were interpolated from the magnification of the integrated solid distribution function.

==Dynamical properties and the diffusion coefficient==

===Task 19: The mean squared displacement(MSD) of large and small simulations===

[[File:msdslv.png|800px|thumb|right|figure 12]]
<math> MSD(t) = <((t_0 + t)-(t_0))>^2 </math> where; <math>t_0 =</math> the initial timestep and <math>t =</math>the new timestep Equation x

The mean squared displacement(equation x) is a measure of how much random motion there is within the system. In this simulation the ability to move randomly is determined by the amount of free space and therefore the density. In the gaseous phase there is a lot of free space, the atoms move freely and encounter only a small number of atoms to collide with thus impeding it's motion; the linear increase of the plot reflects this behaviour. The liquid has a lower degree of random motion than the gas, due closer proximity of the atoms, forces such as drag will now effect the atom's motion as it moves past other atoms. The two first plots show that the atom is able to move in drift, however the liquid experiences more drag and has a smaller gradient. The solid displays the least amount of translational freedom, the close packed structure prevents any random translational motion; it does however show some random movement, arising from vibrations of the constrained atoms. Therefore the plots in figure 12 are a reasonable representation of the random movement of the three phases.

Finding the diffusion coefficient <math>D</math>

<math>D = \frac{1}{6} \frac{\delta <r^2(t)>}{\delta t}> </math> where <math>D =</math> the diffusion coefficient in <math>m s^{-2}</math>,<math><r^2(t) =</math> MSD and <math> t =</math> timestep
The gradient of the slopes give the diffusion coefficent, however this uses a dimensionless timestep denominator, therefore we must convert it to time units by dividing by the timestep(0.002).

The gradient was found by fitting a linear trendline through the function.

Solid

Conversion of the gradient to give time units gave <math>D = 4.17 x 10^{-6} ms^{-2}</math> for the million atom simulation and <math>D = 5.83E-08 x 10^{-7} ms^{-2}</math> for the 3200 atom simulation.

Liquid

Conversion of the gradient to give time units gave <math>D = 8.33 x 10^{-2} ms^{-2}</math> for the million atom simulation and <math>D = 8.33 x 10^{-2} ms^{-2}</math> for the 3200 atom simulation.

Gas

A trendline was fitted to the function, since the simulation took some time to produce the linear behavior it was fitted to be linear with the latter part of the function. Conversion of the gradient to give time units gave <math>D = 3.09 ms^{-2}</math> for the million atom simulation and <math>D = 2.43 ms^{-2}</math> for the 3200 atom simulation.

===Task 20: The velocity auto-correlation function(VACF)===

==Evaluating the 1D harmonic oscillator VACF using the expression for the position as a function of time==

The normalized VACF is given by;

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

We must simplify this expression in order for it to be useful, to do this we must first find expressions for <math>v(t)</math> and <math>v(t + \tau)</math>:

The expression for position as a function of time for the HO is defined as;
<math>x(t) = A\cos(\omega t + \phi)</math>

This can be substituted into the expression of velocity as a function of time;

<math>v(t) = \frac{dx(t)}{dt} = \frac{dA\cos(\omega t + \phi)}{dt}</math>

This is then differentiated using the chain rule;

<math> \frac{dx(t)}{dt} = \frac{dx(t)}{du} x \frac{du}{dt} </math>

Where <math>u = \omega t + \phi</math>

<math>v(t) = -\omega A\sin(\omega t + \phi) </math>

We also need our velocity as a function of time after <math>\tau</math>, this is acheived by substituting in <math>t = t + \tau</math>;

<math> v(t + \tau) = -\omega A\sin(\omega (t + \tau) + \phi)</math>
<math> = -\omega A\sin(\omega t + \omega\tau + \phi) </math>

The two velocity expressions are substituted into the VACF;

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

::::<math> = \frac{\int_{-\infty}^{\infty} \sin(\omega t + \phi) \sin(\omega t + \omega\tau + \phi)\,dt}{\int_{-\infty}^{\infty} \sin(\omega t + \phi)^2\,dt} </math>

Using the trigonometric identity <math> \sin(\alpha + \beta) = \sin(\alpha) \cos(\beta) + \cos \alpha \sin \beta </math> we can re-write our equation;

where <math> \omega t + \phi = \alpha </math> and <math> \omega \tau = \beta </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)\,dt}{ \int_{-\infty}^{\infty} \sin(\omega t + \phi)^2\,dt}</math>

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

The fraction can be expanded;

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

Then since <math>\cos(\omega \tau)</math> does not contain a t term it can be moved outside of the integral;

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

The alike terms cancel to give;

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

Finally <math>\frac{\int_{-\infty}^{\infty}\cos(\omega t + \phi) \sin(\omega \tau)\,dt}{\int_{-\infty}^{\infty}\sin(\omega t + \phi)^2\,dt} = 0</math>, this is because <math>sin</math> and <math>\frac{\cos}{\sin}</math> are odd functions. When odd functions are integrated between symmetrical boundaries the product is zero.

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

Analysis of the solid, liquid and harmonic osillator VACFs against position
[[File:VACFplots.png|800px|thumb|right|figure 13]]
The minima in the VACFs represent a change in the direction of the velocity, this caused by collisions with other atoms.

The liquid and solid plots are similar in general however there are some important differences arising from their physical structures. The solid has a steeper decrease in correlation and so it has an earlier minima than the liquid. This is caused by the tight ordered structure of the solid, other atoms are encountered a lot more quickly since they are closer. Another difference is the deepness of the minima, the solid's is significantly more negative than the liquid. Again this is due to the solid's close packed structure, the correlation becomes more negative beause the atoms are tightly constrained and so the interactions are very strong with other atoms. The atoms in the liquid are free moving and so when they collide they can move to stabilise The decorrelated solid function shows some oscillation about <math>C(\tau)=0</math> this arises from the jostling with other atoms in the fixed lattice and so we see a sort of harmonic oscillator motion. The liquid does not display this behavior because they are free to move and so will be able to "rapidly destroy any oscillatory motion."[x] by the "diffusive motion".

===Task 19: Integrating using the trapezium rule to get D===

The plots of <math> C(\tau) vs time </math> of the three phases in a small and large system were plotted along with their running integral.
[[File:solid3200.png|800px|thumb|right|figure 14]] [[File:solidM.png|800px|thumb|right|figure 15]] [[File:liquid8000.png|800px|thumb|right|figure 16]] [[File:liquidM.png|800px|thumb|right|figure 17]] [[File:gas8000.png|800px|thumb|right|figure 18]] [[File:gasM.png|800px|thumb|right|figure 19]]

Rep:Mod:IOPJKL

2015-12-04T12:54:30Z

2015-12-03T21:15:53Z

Mcd13: /* Evaluating the 1D harmonic oscillator VACF using the expression for the position as a function of time */

'''Liquid Simulations Computational Experiment'''

==Theory==

===Task 1: Calculations of a Harmonic Oscillator using both Classical and Velocity-Verlet solutions===
[[File:All3graphs.png|800px|thumb|right|''Figure 1.a) Positions calculated via the two solutions vs time b) Difference between the positions calculated via different methods against time c) Total energy of the harmonic oscillator calculated using Velocity-Verlet solutions vs time '']]

'''Position calculations from the two methods''': The position <math>x</math> at time <math>t</math> was determined using <math>x(t) = A cos(\omega t + \phi)</math> where; <math>A = 1, \omega = 1</math> and <math> \phi= 0 </math>. These were plotted against time(Figure 1.a).

'''Error analysis of the two positions''': The difference between the two positions was calculated and plotted against time(Figure 1.b). The error increases with time and position, this loss of precision could come from two sources. As the position gets further from the original position the <math>a(t)dt^2</math> part of the equation has less importance [x] also the errors could build up since the previous position is used in the preceding position calculation.

'''Total Energy of the Harmonic-Oscillator vs time''': Values for velocity <math>v</math> and position <math>x</math> were obtained from the Velocity-Verlet calculation, values for the force-constant <math>k</math> and mass <math>m</math> were both given as 1. Using ,math> E= \frac{1}{2}mv^2 + \frac{1}{2}kx^2'' the total energies were calculated and plotted against time(Figure 1.c).

===Task 2: Plot of maxima error positions vs time===
The coordinates of the 5 maximas in figure 2 were estimated and plotted on a new graph. This yielded a linear function where <math>x(t) = 0.0004x - 7*10^{-5}</math>. [[File:MaximaDiff.png|200px|thumb|right|''Figure 2. Linear plot of the maxima error positions(figure 1.b) against time '']]

===Task 3: Timesteps to minimise variation in total energy===

Multiple timesteps were used to calculate the total energy of the harmonic oscillator. The percentage difference was taken between the initial energy assumed to be correct and the energy of the first minima, these results were then plotted(figure 3). In order to keep the change in total energy below 1% a timestep equal to or below 0.2 should be used. The change in total energy of the simulation should be kept to a minimum as the total energy should not vary over time in the harmonic oscillator.[[File:%vstimestep.png|200px|thumb|right|''Figure 3. Plot of percentage energy difference and timestep'']]

===Task 4: Lennard-Jones interaction===

'''1. Calculation of distance at zero-potential, r = r0:'''
<math>\phi\left(r_0\right) = 4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right)</math>
where φ = potential energy r = the distance ε = well depth and σ = distance from zero potential distance
<math>\phi\left(r_0\right)= 4\epsilon \left( \sigma^{12} - \sigma^6 r_0^6 \right)</math>
:<math>= 4\epsilon \left( 1 - \frac{r_0^6}{\sigma^6} \right) = 0</math>
:<math>= 1 - \frac{r_0^6}{\sigma^6} = 0 </math>
:<math> \sigma^6 = r_0^6 </math>
:<math> \sigma = r_0 </math>

2. '''Force at r0:'''

:<math>\mathbf{F}_0 = - \frac{\mathrm{d}\phi\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_0} = - \frac{\mathrm{d}}{\mathrm{d}\mathbf{r}_0} 4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right) = 4\epsilon \left( -12\frac{\sigma^{12}}{r_0^{13}} + 6\frac{\sigma^6}{r_0^7} \right)</math>
Since <math>\sigma = r_0 </math>
:<math>= 4\epsilon\left( -12\frac{\sigma^{12}}{\sigma^{13}} + 6\frac{\sigma^6}{\sigma^7} \right) = 4\epsilon\frac{-6}{\sigma} = -\frac{24\epsilon}{\sigma}</math>

'''3. Finding the equilibrium separation ''req:'' '''

:<math>\mathbf{F}_{eq} = - \frac{\mathrm{d}\phi\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_{eq}} = - \frac{\mathrm{d}}{\mathrm{d}\mathbf{r}_{eq}} 4\epsilon \left( \frac{\sigma^{12}}{r_{eq}^{12}} - \frac{\sigma^6}{r_{eq}^6} \right) = 4\epsilon \left( -12\frac{\sigma^{12}}{r_{eq}^{13}} + 6\frac{\sigma^6}{r_{eq}^7} \right) = -48\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) + 24\epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right) = 0</math>
:<math>48\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) = 24\epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right)</math>
:<math>2\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) = \epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right)</math>
:<math>2\epsilon \left(\sigma^{12}\right) = \epsilon \left(\sigma^6 r_{eq}^6 \right)</math>
:<math>2\epsilon \left(\sigma^{6}\right) = \epsilon \left(\ r_{eq}^6 \right)</math>
:<math>2\left(\sigma^{6}\right) = \left(\ r_{eq}^6 \right)</math>
:<math>\sqrt[6]{2}\sigma = r_{eq}</math>

'''4. Finding the well-depth at equilibrium:'''

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

'''5. Evaluation of integrals:'''

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

:<math>\int_{2}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{2}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{2}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>= [4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_2^{\infty}</math>
:<math>= 4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(2)^{11}} - \frac{1}{-5(2)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(2)^{11}} - \frac{1}{-5(2)^5} \right) = -24.82 KJ/mol </math>

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

:<math>\int_{2.5}^\infty \phi\left(r\right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>= [4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_{2.5}^{\infty}</math>
:<math>= 4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(2.5)^{11}} - \frac{1}{-5(2.5)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(2.5)^{11}} - \frac{1}{-5(2.5)^5} \right) = -8.18 KJ/mol </math>

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

:<math>\int_{3}^\infty \phi\left(r\right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>[4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_{3}^{\infty}</math>
:<math>4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(3)^{11}} - \frac{1}{-5(3)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(3)^{11}} - \frac{1}{-5(3)^5} \right) = -3.29 KJ/mol </math>

===Task 5: Estimating the number of molecules in volumes of H2O===

1. '''Number of molecules in 1mL:'''
:<math> Mr(H^2O) = 18.02g/mol,\rho = 1g/cm^3 </math>
<math>n = \frac{m}{Mr} = \frac{1}{18.02} = 5.55x10^{-2}</math>
<math>Number of molecules = 5.55x10^{-2} x 6.022x10^{23} = 3.34192x10^{22}</math>

2. '''Volume of 1000 molecules:'''

<math>n = \frac{1000}{6.022x10^{23}} = 1.66x10^{-21} mol</math>
<math>m = 1.66x10^{-21} x 18.02 = 2.99x10^{-20}g</math>
<math>V = \frac{2.99x10^{-20}}{1} = 2.99x10^{-20}cm^3</math>

===Task 6: Periodic boundary conditions===

When the atom at <math>(0.5,0.5,0.5)</math> in the cubic simulation box(<math>(0,0,0,) to (1,1,1)</math>) moves along the vector <math>(0.7,0.6,0.2)</math> it ends up at <math>(0.2,0.1,0.7)</math> (when the periodic boundry condition is taken into consideration).

===Task 7: Converting reduced values to real values:===

1. '''Reduced distance'''

:<math> r* = 3.2 = \frac {r}{\sigma} </math>
:<math> r = 3.2 x 0.34x10^{-9} = 1.09x10^-9 = 1.09nm </math>

2. '''Well-depth'''

:<math> \frac {\epsilon}{K_B} = 120K </math>
:<math>\epsilon = 120K_B = \frac {8.69x10^{24}}{6.022x10^{23}} = 997.71 J/mol = 9.98x10^{-1}</math>

3. '''Temperature'''

:<math> T^* = 1.5 = \frac {K_BT}{\epsilon} </math>
:<math>T = \frac {T^* \epsilon}{K_B} = \frac {1.5 x 1.657x10^{-21}}{K_B} = 180K</math>

==Equilibration==

===Task 8: Problems with random starting coordinates in non-solid simulations===

If atoms were given random starting coordinates then they could end up within a contact proximity, they could bump into each other and then they would not only experience the potential of each but also the kinetic energy.

===Task 9: The number density of lattice points calculations:===

The number density ''n'' is given by the following equation
<math> n = \frac {N}{V}</math> where N is the number of atoms/lattice points per unit cell and V is volumes of the unit cell.

1. The number density of a simple cubic lattice with inter lattice point distance of 1.07722

For a simple cubic lattice there is one atom/lattice point, using equation x
:<math> n = \frac {N}{V} = \frac {N}{r^{*3}} = \frac {1}{1.07722^3} = 0.8 </math>

2. '''Side length of a face-centred cubic lattice with a number density of 1.2'''

A face-centred cubic lattice has 4 atoms/lattice points.
<math> r^* = \sqrt[3]{V} = \sqrt[3]{\frac {N}{n}} = \sqrt[3]{\frac {4}{1.2}} = 1.49380 </math>

===Task 10: How many atoms are created by create_atoms 1 box command for a face-centred cubic lattice?===
Since there are 4 atoms per lattice and the box is ten lattice spacings larger than the single lattice spacing then there would be 4000 atoms.

===Task 11: LAMMPS commands for setting the properties of atoms:===

1. Mass command

mass ''I'' ''value'' = mass 1 1.0
This command is used for setting the mass of the atom, in the general command I is the atom type and value is the mass(in reduced units). In the example the atom is type 1 and the mass is 1.0.

2. Pair_style command

pair_style style args = pair_style lj/cut 3.0

The pair_style command is used to describe paired interactions, style is the interaction which is being looked at and args is the argument used for that style. In the example above the style is defining the cutoff distance, that is the distance at which interactions become negligible. The argument for the cutoff distance is 3.0, after which the simulation will not compute.

3. pair_coeff command

pair_coeff ''* * args'' = pair_coeff * * 1.0 1.0
This command is used to define the force field coefficients between pairs of atoms, the asterix terms are to define the coefficients for all atoms, the second part of this command, args, defines the coeffecients.

===Task 12: Specifying the velocity of each atom===

The Velocity-Verlet algorithm is used to specify velocity and position
<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 \ \ (9)</math>

===Task 13: Running the simulation:===
The command $(100/timestep) command allows for the immediate evaluation of the the formula with the pre-defined variable i.e. the timestep. The following line variable n_steps equal floor(100/$({timestep}) will be broken down into individual contributions. The variable component is used to allows for the return to re-run a command, this allows creates a loop. The n_steps part describes the number of timesteps you wish to calculate which is defined in the 'Run Simulation' command. Equal is used to create a mathematical argument. Floor(11/$({timestep}) is useful in setting the maximum timestep bigger than the original one. In the run simulation part first we define the number of timesteps to be considered and then report the floor value i.e. the total time it takes. These lines are essential in creating a loop where each new timestep may be calculated subsequently.
Make it easier to change the timestep, all values dependedent on the timestep will be changed with it and so it reduces possible human errors when changing all variables and also makes it easier and faster.

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

===Task 14: Checking equilibration===

1. Pressure, Temperature and energy against time for simulation ran at 0.001 timestep

PIC fig 4

2. Equilibration

The simulation reaches equilibration at approximately 0.28, this is the timestep at which the simulation becomes relatively constant.

3. Maximising results accuracy whilst keeping the timestep within a reasonable timescale
The best timestep for maximum accuracy and sufficient information is 0.0025. It overlaps with 0.001 and so although it is more than twice the time as 0.001 it maintains a satisfactory level of accuracy. The others give equilibrium energies which are too far off the most accurate value, with the worst being for 0.01 as it diverges quite dramatically.

[[File:Energyplotsvstime.png|400px|thumb|right|''Figure 5. The energy vs time plots of simulations with varying timesteps '']]

==Running simulations under specific conditions==

===Task 15: Finding the expression for factor <math> \gamma </math>===

<math> \gamma</math> is a factor by which the velocity can be corrected.

<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT </math> (1)
<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT' </math> (2)

Where ''T''' is the target temperature and ''T'' is the instantaneous temperature

<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT </math>
<math> \frac {1}{2} m_iv_i^2 = \frac {1}{2}NK_BT </math>
<math> m_iv_i^2 = NK_BT </math>
<math> v_i^2 = \frac {NK_BT}{m_i} </math>
<math> v_i = \sqrt[2]{\frac {NK_BT}{m_i}} </math>

Insert v_i into equation x

<math> \frac {1}{2} \sum m_i \gamma^2 \frac {NK_BT}{m_i} = \frac {3}{2}NK_BT' </math>
<math> \frac {1}{2} m_i \gamma^2 \frac {NK_BT}{m_i} = \frac {1}{2}NK_BT' </math>
<math> \gamma^2NK_BT = NK_BT' </math>
<math> \gamma^2 = \frac {NK_BT'}{NK_BT} </math>
<math> \gamma^2 = \frac {T'}{T} </math>
<math> \gamma = \sqrt[2]{\frac {T'}{T}} </math>

===Task 16: The importance of number, repetition and frequency of average calculations with timesteps===

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

100 sets the number of timesteps after which another thermodynamic value is used to calculate the average, 1000 is the number of thermodynamic values which are used for the average and 1000000 sets the frequency at which averages are calculated.

100000 timesteps will be run.

===Task 17: Equation of state plots===

5 temperatures above the critical temeperature were chosen(1.6, 1.8, 2.0, 3.0 and 5.0) these were all simulated with a pressure of 2.6(''n-1'') and also 3.3(''n-2''). The timestep chose was 0.01 as in previous experiments it was the largest timesteps yielding the most accuracy.
The pressures were plotted against the temperatures along with their error bars.

The ideal gas densities were calculated with the following equation:
<math> \rho = \frac {N}{V} = \frac {P}{K_BT} = \frac {P}{T} </math> where K_B = 1 in reduced terms
[[File:DensityvsT.png|400px|thumb|right|''Figure 6. Densities n-1 and n-2 calculated via simulation and Ideal Gas methods plotted against temperature'']]

The simulated density is lower than that of the ideal gas method, this can be explained by an assumption to satisy an ideal gas. The atoms must not have any forces upon them and also they mustn't take up any space. In the simulation the atoms do take up space and therefore they have a lower density. This discrepancy increases with pressure. The general trend of the two plots is a decrease in heat capacity as temperature is increased.

==Heat Capacity Calculation==

===Task 18: The temperature dependence of the isochoric heat capacity===
[[File:CvVcommand.png|400px|thumb|left|''Figure 7. Cv/V as a function of T* at P= 0.2 and 0.8'']]

A LAAMPS script which could calculate the isochoric heat capacity from calculated thermodynamic values was used figure 7. The Cv/V calculation was done using equation x.

[[File:CvVvsT.png|400px|thumb|right|''Figure 8. Cv/V as a function of T* at P= 0.2 and 0.8'']]
The heat capacity represents the ability to transfer heat for a given temperature or how much energy will effect a given temperature.

The general decrease of Cv/V vs T is what would be expected. At a higher temperature the energy levels are closer together and at a lower temperature they are further apart. A heat increment will have a greater effect on the higher temperature system, it will cause a greater temperature change since more energy levels are accessible, resulting in a low heat capacity. At a lower temperature the same heat increment that caused a large temperature change at higher temperature will not have the same effect, less energy levels are accessible and so the heat capacity will be higher.

The difference between the two pressured systems in ''figure 7'' is what would be expected, the higher pressure plot has a higher isochoric heat capacity. At a higher pressure the atoms are closer together, they can thus transfer heat more efficiently with a given temperature change than a system with a lower temperature since the atoms are further away.

==Structural properties and the radial distribution function==
===Task 18: The Radial Distribution Functions of a solid, liquid and gas===

'''1. Comparisons of the solid, liquid and gas radial distributions functions as well as the integrated function '''
[[File:rdfslg.png|400px|thumb|left|figure 9]]
[[File:rdfint.png|400px|thumb|left|figure 10]]

The distribution plot(figure 9) of the gas has one short, broad peak which is close to the reference atom; as you move further away the distribution becomes constant. This characteristic shows the small probability of finding an atom close to the reference atom; also when it becomes constant this shows how there is negligible probability of finding another atom at the distance studied. When looking at the integrated function(figure 10), it emphasizes how few atoms you will find in the distance, this is characteristic of a gas' dispersed atoms.

The liquid's distribution has a taller first peak than the gas' and then multiple peaks tapering off as you move further away; it is evident that there is a higher probability of finding an atom throughout the space than in a gas. Additionally the integrated distribution(figure 10) shows a larger number of atoms as you move through the space. This demonstrates the higher proximity of the atoms in the liquid state.

Finally the solid's distribution has a very sharp and tall peak close to the reference atom but then as you move further away the peaks become smaller but continuously fluctuate. The sharpness of the peaks and their continuous distribution arise from the highly ordered structure of the solid. The peaks are sharp unlike the gas and liquid because all of the atoms are in fixed positions and so at a certain position there will be a number of atoms in an equivalent position around the sphere of the distribution function. Futhermore the solid's first peak is considerably larger than the other phases, the probibility is much greater sdue to the closeness of the atoms. The integrated function shows a very large number of atoms in the space, this again illustrates the closeness of the atoms in the solid state.

'''2. Using the solid RDF to extract information on the lattice structure'''

[[File:rdf3d.png|200px|thumb|right|figure 11]]

The three first peaks of the solid's radial distribution function are distances 1.0, 1.5 and 1.8. The lattice spacing is 1.5 and the calculated distances ra, rb and rc(where r is the distance from the reference atom to the atoms a, b and c) were found to be 1.0, 1.5 and 1.8(respectively). Therefore the first three peaks 1, 2 and 3 correspond to the lattice points a, b and c(respectively).

The coordination of atom a, b and c is 12, 6 and 24(respectively). These values were interpolated from the magnification of the integrated solid distribution function.

==Dynamical properties and the diffusion coefficient==

===Task 19: The mean squared displacement(MSD) of large and small simulations===

[[File:msdslv.png|800px|thumb|right|figure 12]]
<math> MSD(t) = <((t_0 + t)-(t_0))>^2 </math> where; <math>t_0 =</math> the initial timestep and <math>t =</math>the new timestep Equation x

The mean squared displacement(equation x) is a measure of how much random motion there is within the system. In this simulation the ability to move randomly is determined by the amount of free space and therefore the density. In the gaseous phase there is a lot of free space, the atoms move freely and encounter only a small number of atoms to collide with thus impeding it's motion; the linear increase of the plot reflects this behaviour. The liquid has a lower degree of random motion than the gas, due closer proximity of the atoms, forces such as drag will now effect the atom's motion as it moves past other atoms. The two first plots show that the atom is able to move in drift, however the liquid experiences more drag and has a smaller gradient. The solid displays the least amount of translational freedom, the close packed structure prevents any random translational motion; it does however show some random movement, arising from vibrations of the constrained atoms. Therefore the plots in figure 12 are a reasonable representation of the random movement of the three phases.

Finding the diffusion coefficient <math>D</math>

<math>D = \frac{1}{6} \frac{\delta <r^2(t)>}{\delta t}> </math> where <math>D =</math> the diffusion coefficient in <math>m s^{-2}</math>,<math><r^2(t) =</math> MSD and <math> t =</math> timestep
The gradient of the slopes give the diffusion coefficent, however this uses a dimensionless timestep denominator, therefore we must convert it to time units by dividing by the timestep(0.002).

The gradient was found by fitting a linear trendline through the function.

Solid

Conversion of the gradient to give time units gave <math>D = 4.17 x 10^{-6} ms^{-2}</math> for the million atom simulation and <math>D = 5.83E-08 x 10^{-7} ms^{-2}</math> for the 3200 atom simulation.

Liquid

Conversion of the gradient to give time units gave <math>D = 8.33 x 10^{-2} ms^{-2}</math> for the million atom simulation and <math>D = 8.33 x 10^{-2} ms^{-2}</math> for the 3200 atom simulation.

Gas

A trendline was fitted to the function, since the simulation took some time to produce the linear behavior it was fitted to be linear with the latter part of the function. Conversion of the gradient to give time units gave <math>D = 3.09 ms^{-2}</math> for the million atom simulation and <math>D = 2.43 ms^{-2}</math> for the 3200 atom simulation.

===Task 20: The velocity auto-correlation function(VACF)===

==Evaluating the 1D harmonic oscillator VACF using the expression for the position as a function of time==

The normalized VACF is given by;

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

From task 2 we define <math>x(t),/math>

<math>x(t) = A\cos(\omega t + \phi)</math>

<math>v(t) = \frac{dx(t)}{dt} = \frac{dA\cos(\omega t + \phi}{dt}</math>

Using the chain rule
<math> \frac{dx(t)}{dt} = \frac{dx(t)}{du} x \frac{du}{dt} </math>
Where <math>u = \omega t + \phi</math>
<math>v(t) = -\omega A\sin(\omega t + \phi) </math>
For <math>v(t + \tau </math> we let <math>t = t + \tau </math>
Therefore;
<math> v(t + \tau) = -\omega A\sin(\omega (t + \tau) + \phi)</math>
<math> = -\omega A\sin(\omega t + \omega\tau + \phi) </math>

The two velocity expressions are substituted into the velocity autocorrelation function;

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

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

Since <math> \sin(\alpha + \beta) = \sin(\alpha) \cos(\beta) + \cos \alpha \sin \beta </math>

where <math> \omega t + \phi = \alpha </math> and <math> \omega \tau = \beta </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)\,dt}{ \int_{-\infty}^{\infty}\sin(\omega t + \phi)^2\,dt}</math>
<math> = \frac{\sin(\omega t + \phi)^2 \cos(\omega\tau) + \sin(\omega t + \phi) \cos(\omega t + \phi) \sin(\omega \tau)\,dt}{\int_{\-infty}^{\infty}(\sin(\omega t + \phi)^2dt} </math>

= \frac{(\sin(\omega t + \phi)^2 \cos(\omega\tau))dt}{\int_{\-infty}^{\infty}(\sin(\omega t + \phi)^2dt} + \frac{\sin(\omega t + \phi) \cos(\omega t + \phi) \sin(\omega \tau)dt}{\int_{\-infty}^{\infty}(\sin(\omega t + \phi)^2dt}

<math> = </math>

Plotting and analysis of the solid, liquid and harmonic osillator VACFs against position

Comparison and analysis of the solid and liquid VACFs

The VACFs for the two simulations in figure 13 have three main features of interest at different distances; the initial part at a minuscule distance is the ''free flight'' part where the function is at 0, this is where there are no restrictive forces this then starts to decay when the atom encounters surrounding atoms decreasing the function. In the solid we can see it starts decaying faster than the liquid, this is due to the higher number of surrounding atoms. The next significant behaviour is at a moderate distance, where the function reaches zero. At zero there are no residual correlations and so at this point the function starts to increase to the constant uncorrelated function since C(t) = 0. We can see that it takes the solid a lot less time for the atom to become uncorrelated than the liquid. This is because there is much less random motion in the solid's ordered structure than the liquid. Another interesting part of the function is the decay into a negative function, this behavior is called back-scattering, where the atom moves about from interactions with neighboring atoms. The atom in the solid experiences a larger degree of back-scattering, since it is confined to a tight face-centered lattice and jostles about with it's nearest neighbors. Where the functions become constant we can see that the solid has fluctuations around it's uncorrelated function, this arises from atomic vibrations of the close packed atoms; a small fluctuation may also be seen for the liquid, however it experiences less vibrational motion since the atoms are not locked in a tight configuration.

Comparison and analysis of the solid and liquid VACFs

The harmonic oscillator models a diatomic atom connected by a spring, the spring expands and contracts in a periodic fashion, resulting in a repetitive back and forth motion which is fixed. The Lennard-Jones potential models the interaction between two atoms, there is an attractive and negative term which are dependent on distance and which are not fixed. This results in a periodic behaviour, where the atom goes from free flight, to zero correlation with a lot of backscattering and then back up to free flight. This periodic nature arises from the fixed springs constrained motion. The lennard-jones atom has a more random behavior, since the interatomic distance is not fixed the correlation is highly dependent on distance.

===Task 19: Integrating using the trapezium rule to get D===

The plots of <math> C(\tau) vs time </math> of the three phases in a small and large system were plotted along with their running integral.

File:Solid3200.png

2015-12-03T20:16:05Z

Mcd13:

File:Liquid8000.png

2015-12-03T20:15:32Z

Mcd13:

2015-12-02T17:17:43Z

Mcd13:

'''Liquid Simulations Computational Experiment'''

=='''Theory'''==

===Task 1: Calculations of a Harmonic Oscillator using both Classical and Velocity-Verlet solutions===
[[File:All3graphs.png|800px|thumb|right|''Figure 1.a) Positions calculated via the two solutions vs time b) Difference between the positions calculated via different methods against time c) Total energy of the harmonic oscillator calculated using Velocity-Verlet solutions vs time '']]

'''Position calculations from the two methods''': The position <math>x</math> at time <math>t</math> was determined using <math>x(t) = A cos(\omega t + \phi)</math> where; <math>A = 1, \omega = 1</math> and <math> \phi= 0 </math>. These were plotted against time(Figure 1.a).

'''Error analysis of the two positions''': The difference between the two positions was calculated and plotted against time(Figure 1.b). The error increases with time and position, this loss of precision could come from two sources. As the position gets further from the original position the <math>a(t)dt^2</math> part of the equation has less importance [x] also the errors could build up since the previous position is used in the preceding position calculation.

'''Total Energy of the Harmonic-Oscillator vs time''': Values for velocity <math>v</math> and position <math>x</math> were obtained from the Velocity-Verlet calculation, values for the force-constant <math>k</math> and mass <math>m</math> were both given as 1. Using ,math> E= \frac{1}{2}mv^2 + \frac{1}{2}kx^2'' the total energies were calculated and plotted against time(Figure 1.c).

===Task 2: Plot of maxima error positions vs time===
The coordinates of the 5 maximas in figure 2 were estimated and plotted on a new graph. This yielded a linear function where <math>x(t) = 0.0004x - 7*10^{-5}</math>. [[File:MaximaDiff.png|200px|thumb|right|''Figure 2. Linear plot of the maxima error positions(figure 1.b) against time '']]

===Task 3: Timesteps to minimise variation in total energy===

Multiple timesteps were used to calculate the total energy of the harmonic oscillator. The percentage difference was taken between the initial energy assumed to be correct and the energy of the first minima, these results were then plotted(figure 3). In order to keep the change in total energy below 1% a timestep equal to or below 0.2 should be used. The change in total energy of the simulation should be kept to a minimum as the total energy should not vary over time in the harmonic oscillator.[[File:%vstimestep.png|200px|thumb|right|''Figure 3. Plot of percentage energy difference and timestep'']]

===Task 4: Lennard-Jones interaction===

'''1. Calculation of distance at zero-potential, r = r0:'''
<math>\phi\left(r_0\right) = 4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right)</math>
where φ = potential energy r = the distance ε = well depth and σ = distance from zero potential distance
<math>\phi\left(r_0\right)= 4\epsilon \left( \sigma^{12} - \sigma^6 r_0^6 \right)</math>
:<math>= 4\epsilon \left( 1 - \frac{r_0^6}{\sigma^6} \right) = 0</math>
:<math>= 1 - \frac{r_0^6}{\sigma^6} = 0 </math>
:<math> \sigma^6 = r_0^6 </math>
:<math> \sigma = r_0 </math>

2. '''Force at r0:'''

:<math>\mathbf{F}_0 = - \frac{\mathrm{d}\phi\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_0} = - \frac{\mathrm{d}}{\mathrm{d}\mathbf{r}_0} 4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right) = 4\epsilon \left( -12\frac{\sigma^{12}}{r_0^{13}} + 6\frac{\sigma^6}{r_0^7} \right)</math>
Since <math>\sigma = r_0 </math>
:<math>= 4\epsilon\left( -12\frac{\sigma^{12}}{\sigma^{13}} + 6\frac{\sigma^6}{\sigma^7} \right) = 4\epsilon\frac{-6}{\sigma} = -\frac{24\epsilon}{\sigma}</math>

'''3. Finding the equilibrium separation ''req:'' '''

:<math>\mathbf{F}_{eq} = - \frac{\mathrm{d}\phi\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_{eq}} = - \frac{\mathrm{d}}{\mathrm{d}\mathbf{r}_{eq}} 4\epsilon \left( \frac{\sigma^{12}}{r_{eq}^{12}} - \frac{\sigma^6}{r_{eq}^6} \right) = 4\epsilon \left( -12\frac{\sigma^{12}}{r_{eq}^{13}} + 6\frac{\sigma^6}{r_{eq}^7} \right) = -48\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) + 24\epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right) = 0</math>
:<math>48\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) = 24\epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right)</math>
:<math>2\epsilon \left(\frac{\sigma^{12}}{r_{eq}^{13}}\right) = \epsilon \left(\frac{\sigma^6}{r_{eq}^7} \right)</math>
:<math>2\epsilon \left(\sigma^{12}\right) = \epsilon \left(\sigma^6 r_{eq}^6 \right)</math>
:<math>2\epsilon \left(\sigma^{6}\right) = \epsilon \left(\ r_{eq}^6 \right)</math>
:<math>2\left(\sigma^{6}\right) = \left(\ r_{eq}^6 \right)</math>
:<math>\sqrt[6]{2}\sigma = r_{eq}</math>

'''4. Finding the well-depth at equilibrium:'''

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

'''5. Evaluation of integrals:'''

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

:<math>\int_{2}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{2}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{2}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>= [4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_2^{\infty}</math>
:<math>= 4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(2)^{11}} - \frac{1}{-5(2)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(2)^{11}} - \frac{1}{-5(2)^5} \right) = -24.82 KJ/mol </math>

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

:<math>\int_{2.5}^\infty \phi\left(r\right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{2.5}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>= [4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_{2.5}^{\infty}</math>
:<math>= 4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(2.5)^{11}} - \frac{1}{-5(2.5)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(2.5)^{11}} - \frac{1}{-5(2.5)^5} \right) = -8.18 KJ/mol </math>

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

:<math>\int_{3}^\infty \phi\left(r\right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{1^{12}}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r = \int_{3}^\infty 4\epsilon \left( \frac{1}{r^{12}} - \frac{1}{r^6} \right)\mathrm{d}r</math>
:<math>[4\epsilon \left( \frac{1}{-11r^{11}} - \frac{1}{-5r^5} \right)]_{3}^{\infty}</math>
:<math>4\epsilon \left( \frac{1}{-11(\infty)^{11}} - \frac{1}{-5(\infty)^5} \right) - 4\epsilon \left( \frac{1}{-11(3)^{11}} - \frac{1}{-5(3)^5} \right)</math>
:<math>= -0 + 4\epsilon \left( \frac{1}{-11(3)^{11}} - \frac{1}{-5(3)^5} \right) = -3.29 KJ/mol </math>

===Task 5: Estimating the number of molecules in volumes of H2O===

1. '''Number of molecules in 1mL:'''
:<math> Mr(H^2O) = 18.02g/mol,\rho = 1g/cm^3 </math>
<math>n = \frac{m}{Mr} = \frac{1}{18.02} = 5.55x10^{-2}</math>
<math>Number of molecules = 5.55x10^{-2} x 6.022x10^{23} = 3.34192x10^{22}</math>

2. '''Volume of 1000 molecules:'''

<math>n = \frac{1000}{6.022x10^{23}} = 1.66x10^{-21} mol</math>
<math>m = 1.66x10^{-21} x 18.02 = 2.99x10^{-20}g</math>
<math>V = \frac{2.99x10^{-20}}{1} = 2.99x10^{-20}cm^3</math>

===Task 6: Periodic boundary conditions===

When the atom at <math>(0.5,0.5,0.5)</math> in the cubic simulation box(<math>(0,0,0,) to (1,1,1)</math>) moves along the vector <math>(0.7,0.6,0.2)</math> it ends up at <math>(0.2,0.1,0.7)</math> (when the periodic boundry condition is taken into consideration).

===Task 7: Converting reduced values to real values:===

1. '''Reduced distance'''

:<math> r* = 3.2 = \frac {r}{\sigma} </math>
:<math> r = 3.2 x 0.34x10^{-9} = 1.09x10^-9 = 1.09nm </math>

2. '''Well-depth'''

:<math> \frac {\epsilon}{K_B} = 120K </math>
:<math>\epsilon = 120K_B = \frac {8.69x10^{24}}{6.022x10^{23}} = 997.71 J/mol = 9.98x10^{-1}</math>

3. '''Temperature'''

:<math> T^* = 1.5 = \frac {K_BT}{\epsilon} </math>
:<math>T = \frac {T^* \epsilon}{K_B} = \frac {1.5 x 1.657x10^{-21}}{K_B} = 180K</math>

==Equilibration==

===Task 8: Problems with random starting coordinates in non-solid simulations===

If atoms were given random starting coordinates then they could end up within a contact proximity, they could bump into each other and then they would not only experience the potential of each but also the kinetic energy.

===Task 9: The number density of lattice points calculations:===

The number density ''n'' is given by the following equation
<math> n = \frac {N}{V}</math> where N is the number of atoms/lattice points per unit cell and V is volumes of the unit cell.

1. The number density of a simple cubic lattice with inter lattice point distance of 1.07722

For a simple cubic lattice there is one atom/lattice point, using equation x
:<math> n = \frac {N}{V} = \frac {N}{r^{*3}} = \frac {1}{1.07722^3} = 0.8 </math>

2. '''Side length of a face-centred cubic lattice with a number density of 1.2'''

A face-centred cubic lattice has 4 atoms/lattice points.
<math> r^* = \sqrt[3]{V} = \sqrt[3]{\frac {N}{n}} = \sqrt[3]{\frac {4}{1.2}} = 1.49380 </math>

===Task 10: How many atoms are created by create_atoms 1 box command for a face-centred cubic lattice?===
Since there are 4 atoms per lattice and the box is ten lattice spacings larger than the single lattice spacing then there would be 4000 atoms.

===Task 11: LAMMPS commands for setting the properties of atoms:===

1. Mass command

mass ''I'' ''value'' = mass 1 1.0
This command is used for setting the mass of the atom, in the general command I is the atom type and value is the mass(in reduced units). In the example the atom is type 1 and the mass is 1.0.

2. Pair_style command

pair_style style args = pair_style lj/cut 3.0

The pair_style command is used to describe paired interactions, style is the interaction which is being looked at and args is the argument used for that style. In the example above the style is defining the cutoff distance, that is the distance at which interactions become negligible. The argument for the cutoff distance is 3.0, after which the simulation will not compute.

3. pair_coeff command

pair_coeff ''* * args'' = pair_coeff * * 1.0 1.0
This command is used to define the force field coefficients between pairs of atoms, the asterix terms are to define the coefficients for all atoms, the second part of this command, args, defines the coeffecients.

===Task 12: Specifying the velocity of each atom===

The Velocity-Verlet algorithm is used to specify velocity and position
<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 \ \ (9)</math>

===Task 13: Running the simulation:===
The command $(100/timestep) command allows for the immediate evaluation of the the formula with the pre-defined variable i.e. the timestep. The following line variable n_steps equal floor(100/$({timestep}) will be broken down into individual contributions. The variable component is used to allows for the return to re-run a command, this allows creates a loop. The n_steps part describes the number of timesteps you wish to calculate which is defined in the 'Run Simulation' command. Equal is used to create a mathematical argument. Floor(11/$({timestep}) is useful in setting the maximum timestep bigger than the original one. In the run simulation part first we define the number of timesteps to be considered and then report the floor value i.e. the total time it takes. These lines are essential in creating a loop where each new timestep may be calculated subsequently.
Make it easier to change the timestep, all values dependedent on the timestep will be changed with it and so it reduces possible human errors when changing all variables and also makes it easier and faster.

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

===Task 14: Checking equilibration===

1. Pressure, Temperature and energy against time for simulation ran at 0.001 timestep

PIC fig 4

2. Equilibration

The simulation reaches equilibration at approximately 0.28, this is the timestep at which the simulation becomes relatively constant.

3. Maximising results accuracy whilst keeping the timestep within a reasonable timescale
The best timestep for maximum accuracy and sufficient information is 0.0025. It overlaps with 0.001 and so although it is more than twice the time as 0.001 it maintains a satisfactory level of accuracy. The others give equilibrium energies which are too far off the most accurate value, with the worst being for 0.01 as it diverges quite dramatically.

[[File:Energyplotsvstime.png|400px|thumb|right|''Figure 5. The energy vs time plots of simulations with varying timesteps '']]

==Running simulations under specific conditions==

===Task 15: Finding the expression for factor <math> \gamma </math>===

<math> \gamma</math> is a factor by which the velocity can be corrected.

<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT </math> (1)
<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT' </math> (2)

Where ''T''' is the target temperature and ''T'' is the instantaneous temperature

<math> \frac {1}{2} \sum m_i v_i^2 = \frac {3}{2}NK_BT </math>
<math> \frac {1}{2} m_iv_i^2 = \frac {1}{2}NK_BT </math>
<math> m_iv_i^2 = NK_BT </math>
<math> v_i^2 = \frac {NK_BT}{m_i} </math>
<math> v_i = \sqrt[2]{\frac {NK_BT}{m_i}} </math>

Insert v_i into equation x

<math> \frac {1}{2} \sum m_i \gamma^2 \frac {NK_BT}{m_i} = \frac {3}{2}NK_BT' </math>
<math> \frac {1}{2} m_i \gamma^2 \frac {NK_BT}{m_i} = \frac {1}{2}NK_BT' </math>
<math> \gamma^2NK_BT = NK_BT' </math>
<math> \gamma^2 = \frac {NK_BT'}{NK_BT} </math>
<math> \gamma^2 = \frac {T'}{T} </math>
<math> \gamma = \sqrt[2]{\frac {T'}{T}} </math>

===Task 16: The importance of number, repetition and frequency of average calculations with timesteps===

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

100 sets the number of timesteps after which another thermodynamic value is used to calculate the average, 1000 is the number of thermodynamic values which are used for the average and 1000000 sets the frequency at which averages are calculated.

100000 timesteps will be run.

===Task 17: Equation of state plots===

5 temperatures above the critical temeperature were chosen(1.6, 1.8, 2.0, 3.0 and 5.0) these were all simulated with a pressure of 2.6(''n-1'') and also 3.3(''n-2''). The timestep chose was 0.01 as in previous experiments it was the largest timesteps yielding the most accuracy.
The pressures were plotted against the temperatures along with their error bars.

The ideal gas densities were calculated with the following equation:
<math> \rho = \frac {N}{V} = \frac {P}{K_BT} = \frac {P}{T} </math> where K_B = 1 in reduced terms
[[File:DensityvsT.png|400px|thumb|right|''Figure 6. Densities n-1 and n-2 calculated via simulation and Ideal Gas methods plotted against temperature'']]

The simulated density is lower than that of the ideal gas method, this can be explained by an assumption to satisy an ideal gas. The atoms must not have any forces upon them and also they mustn't take up any space. In the simulation the atoms do take up space and therefore they have a lower density. This discrepancy increases with pressure. The general trend of the two plots is a decrease in heat capacity as temperature is increased.

==Heat Capacity Calculation==

===Task 18: The temperature dependence of the isochoric heat capacity===
[[File:CvVcommand.png|400px|thumb|left|''Figure 7. Cv/V as a function of T* at P= 0.2 and 0.8'']]

A LAAMPS script which could calculate the isochoric heat capacity from calculated thermodynamic values was used figure 7. The Cv/V calculation was done using equation x.

[[File:CvVvsT.png|400px|thumb|right|''Figure 8. Cv/V as a function of T* at P= 0.2 and 0.8'']]
The heat capacity represents the ability to transfer heat for a given temperature or how much energy will effect a given temperature.

The general decrease of Cv/V vs T is what would be expected. At a higher temperature the energy levels are closer together and at a lower temperature they are further apart. A heat increment will have a greater effect on the higher temperature system, it will cause a greater temperature change since more energy levels are accessible, resulting in a low heat capacity. At a lower temperature the same heat increment that caused a large temperature change at higher temperature will not have the same effect, less energy levels are accessible and so the heat capacity will be higher.

The difference between the two pressured systems in ''figure 7'' is what would be expected, the higher pressure plot has a higher isochoric heat capacity. At a higher pressure the atoms are closer together, they can thus transfer heat more efficiently with a given temperature change than a system with a lower temperature since the atoms are further away.

==Structural properties and the radial distribution function==
===Task 18: The Radial Distribution Functions of a solid, liquid and gas===

'''1. Comparisons of the solid, liquid and gas radial distributions functions as well as the integrated function '''
[[File:rdfslg.png|400px|thumb|left|figure 9]]
[[File:rdfint.png|400px|thumb|left|figure 10]]

The distribution plot(figure 9) of the gas has one short, broad peak which is close to the reference atom; as you move further away the distribution becomes constant. This characteristic shows the small probability of finding an atom close to the reference atom; also when it becomes constant this shows how there is negligible probability of finding another atom at the distance studied. When looking at the integrated function(figure 10), it emphasizes how few atoms you will find in the distance, this is characteristic of a gas' dispersed atoms.

The liquid's distribution has a taller first peak than the gas' and then multiple peaks tapering off as you move further away; it is evident that there is a higher probability of finding an atom throughout the space than in a gas. Additionally the integrated distribution(figure 10) shows a larger number of atoms as you move through the space. This demonstrates the higher proximity of the atoms in the liquid state.

Finally the solid's distribution has a very sharp and tall peak close to the reference atom but then as you move further away the peaks become smaller but continuously fluctuate. The sharpness of the peaks and their continuous distribution arise from the highly ordered structure of the solid. The peaks are sharp unlike the gas and liquid because all of the atoms are in fixed positions and so at a certain position there will be a number of atoms in an equivalent position around the sphere of the distribution function. Futhermore the solid's first peak is considerably larger than the other phases, the probibility is much greater sdue to the closeness of the atoms. The integrated function shows a very large number of atoms in the space, this again illustrates the closeness of the atoms in the solid state.

'''2. Using the solid RDF to extract information on the lattice structure'''

[[File:rdf3d.png|200px|thumb|right|figure 11]]

The three first peaks of the solid's radial distribution function are distances 1.0, 1.5 and 1.8. The lattice spacing is 1.5 and the calculated distances ra, rb and rc(where r is the distance from the reference atom to the atoms a, b and c) were found to be 1.0, 1.5 and 1.8(respectively). Therefore the first three peaks 1, 2 and 3 correspond to the lattice points a, b and c(respectively).

The coordination of atom a, b and c is 12, 6 and 24(respectively). These values were interpolated from the magnification of the integrated solid distribution function.

==Dynamical properties and the diffusion coefficient==

===Task 19: The mean squared displacement of large and small simulations===

[[File:msdslv.png|800px|thumb|right|figure 12]]

The mean squared displacement(equation x) is a measure of how much random motion there is within the system. In this simulation the ability to move randomly is determined by the amount of free space and therefore the density. In the gaseous phase there is a lot of free space, the atoms move freely and encounter only a small number of atoms to collide with thus impeding it's motion; the linear increase of the plot reflects this behaviour. The liquid has a lower degree of random motion than the gas, due closer proximity of the atoms, forces such as drag will now effect the atom's motion as it moves past other atoms. The two first plots show that the atom is able to move in drift, however the liquid experiences more drag and has a smaller gradient. The solid displays the least amount of translational freedom, the close packed structure prevents any random translational motion; it does however show some random movement, arising from vibrations of the constrained atoms. Therefore the plots in figure 12 are a reasonable representation of the random movement of the three phases.