= Molecular dynamics simulation =

=== Tasks 1,2 and 3: Velocity-Verlet algorithm: ===

:The variation of the position of an atom was calculated using the classical harmonic oscillator approach and the Velocity-Verlet algorithm. Image 1 shows the difference between the two methods.
The orange line traces de maximum error peaks. It can be seen that these maximum errors increase linearly (<math> y= 2*10^{-5}x + 3*10^{5}</math>) with time. The reason behind this is that the simulation of each time-step requires that of the previous one, so the error builds up.

[[File:Al3214error.jpg|800px||alt=Alt text|Caption text]]

For the following simulations, since the energy of an ensemble should be constant over time we needed to find the value of the time-step that would give the lowest standard deviation in order words, the least energy fluctuations. The following table shows that 0.03 gives a 1% STD.
{| class="wikitable"
|+ Energy Fluctuation with time.
!Time-step
!Standard deviation
|-
|0.1
|0.18
|-
|0.01
|0.17
|-
|0.001
|0.003
|-
|0.002
|0.012
|}

=== Task 4: Lennard Jones Potential:===
:*'''Finding <math>r_0</math> at which the potential is zero''':

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

:: Since <math>\phi\left(r_o\right)= 0 </math>:

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

::And:

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

:<math> \sigma^{6} r_0^{12} = \sigma^{12} r_0^{6}</math>

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

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

:<math> \sigma^{6}= r_0^{6}</math>

::'''Finally <math> \sigma= r_0</math>'''

:*'''Force at seperation at <math>r_0</math> equal to <math> \frac{d(\phi\left(r_0\right))}{dr_0}</math>'''

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

:<math> F = \frac{-4\epsilon\sigma^{12}12}{r_0^{13}} + \frac{4\epsilon\sigma^{6}6}{r_0^{7}}</math>

:<math>F=\frac{-48\epsilon\sigma^{12}}{r_0^{13}} + \frac{24\epsilon\sigma^{6}}{r_0^{7}}</math>

:'''Since <math> \sigma= r_0</math> at this point''',

:<math>F=\frac{-48\epsilon r_0^{12}}{r_0^{13}} + \frac{24\epsilon r_0^{6}}{r_0^{7}}</math>

:<math>F=\frac{-48\epsilon+24\epsilon r_0^{6}}{r_0^{7}}</math>

:<math>F=\frac{-24\epsilon}{r_0}</math> or <math>F=\frac{-24\epsilon}{\sigma}</math>

:*'''Finding the equilibrium separation:'''
: The equilibrium separation <math> r_{eq}</math> corresponds to the energy minimum of the Lennard Jones potential curve. As it is the minimum it will be a solution of equating the first derivative of <math>\phi\left(r\right)</math> to zero.
:Calculating <math> \frac{d(\phi\left(r\right))}{dr}</math> and knowing it has to equal to zero.

::<math> \frac{d(\phi\left(r\right))}{dr} = \frac{-48\epsilon \sigma^{12}}{r^{13}} + \frac{24\epsilon \sigma^{6}}{r^{7}}</math> = <math>\frac{-48\epsilon \sigma^{12} r^{7}}{r^{20}} + \frac{24\epsilon \sigma^{6}r^{13}}{r^{20}}</math> = <math>\frac{-48\epsilon \sigma^{12} r^{7} +24\epsilon \sigma^{6}r^{13}}{r^{20}}</math>
::<math>\frac{-24\epsilon \sigma^{6} r^{7}(-2\sigma^{6}+r^{6})}{r^{20}} = 0</math>
::<math>(-2\sigma^{6}+r^{6})=0</math>
::<math>2\sigma^{6}=r^{6}</math>
:: '''<math> r_{eq}=\sqrt[6]{2}\sigma</math>''''

:* Now the well depth can be found
::<math> \phi(r_{eq}) =4\epsilon( \frac{\sigma^{12}}{\sqrt[6]{2}^{12}\sigma^{12}} - \frac{\sigma^6}{\sqrt[6]{2}^{6}\sigma^{6}}) </math> = <math> 4\epsilon( \frac{1}{2^{2}} - \frac{1}{2}) </math> =<math> 4\epsilon( \frac{-1}{4}) </math>=<math> -\epsilon </math>

*Evaluating integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>. With values of <math>\sigma = \epsilon = 1.0</math>
<math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r = \int_{2\sigma}^\infty\mathrm{4\epsilon \left( \frac{\sigma^{12}}{r_0^{12}} - \frac{\sigma^6}{r_0^6} \right)d}r = [ \frac{1}{4}(-\frac{\sigma^{12}}{11*2^{11}\sigma^{11}}+\frac{\sigma^{6}}{5*2^{5}\sigma^{5}})]^{\infty}_{2\sigma} </math> = <math>\frac{1}{4}(-\frac{1^{12}}{11*2^{11}1^{11}}+\frac{1^{6}}{5*2^{5}1^{5}})= -0.0248</math>

:The same process was followed for the other integrals:
:<math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r = -\frac{1}{4}\epsilon(-\frac{\sigma^{12}}{11*2.5^{11}\sigma^{11}}+\frac{\sigma^{6}}{5*2.5^{5}\sigma^{5}})= -0.00818</math>
:<math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r = -\frac{1}{4}\epsilon(-\frac{\sigma^{12}}{11*3^{11}\sigma^{11}}+\frac{\sigma^{6}}{5*3^{5}\sigma^{5}})= -0.00329</math>

===Task 5: Number of water molecules in 1ml of water and volume of 1000 water molecules===
:1mL of water equals 1 gr of water. Dividing this gram by 18.01528 g/mol, the molar density of water we get 0.0555 mol of water. Finally by multiplying by Avogadro's number we find that there are <math>3.34228 10^{22}</math> molecules of water in 1 mL.
: Proportionally, 10000 molecules will occupy a volume of <math>2.991*10^{-19}</math> mL.

=== Task 6: ===
: An atom at position <math>\left(0.5, 0.5, 0.5\right)</math> sits in a cubic simulation box of dimensions from <math>\left(0, 0, 0\right)</math> to <math>\left(1, 1, 1\right)</math> with periodic boundary conditions. If this atom moves along the vector <math>\left(0.7, 0.6, 0.2\right)</math> it will end up at <math>\left(0.2, 0.1, 0.7\right)</math>.

=== Task 7: Reduced Units ===

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

''<math>r = 1.088</math>'''

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

'''<math>\epsilon\ = 0.014 {kJ\ mol}^{-1}</math> '''

===Task 8: Density and lattice spacing===
: Using the formula for density:
::<math> \rho = \frac{m}{V}</math><math>=\frac{1}{1.077^{3}}</math><math>= 0.79999\approx 0.8</math>
: If we consider a fcc lattice, (4 atoms per unit cell), with density of 1.2 and using the formula above:
::<math> 1.2 = \frac{4}{x^{3}}</math>
:Where x is the side length of a unit cell
::<math> x=\sqrt[3]{\frac{4}{1.2}}\approx 1.275</math>

===Task 9 Number of atoms with an fcc lattice===
With the established set-up, 1000 unit cells are created. If an fcc lattice had been chosen 1000 unit cells would have also been formed each containing 4 atoms which add up to a total of 4000 atoms.

===Task 10: Input script commands. ===
:The commands are used to determine the properties of the atoms.
<pre> mass 1 1.0</pre> Sets the mass for all atoms of type 1 to a value of 1.0. It is worth remembering that for this simulation only 1 atom type has been used.

<pre>pair_style lj/cut 3.0</pre> Calculates the coulombic pairwise interaction of two neighbouring atoms as the distance increases between them until a maximum of 0.3 in the Lennard-Jones potential.

<pre>pair_coeff * * 1.0 1.0</pre> Specifies the pairwise force field coefficients between atoms pairs. The asterisks indicate that all atoms types from 1 to N will have a value of 1.

===Task 11: Integration algorithm ===

By specifying <math>\mathbf{x}_i\left(0\right)</math> and <math>\mathbf{v}_i\left(0\right)</math>, we are using the Velocity-Verlet algorithm.

===Task 12: Understanding the code===
<pre>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</pre>

The code above is used to constrain the number of steps depending on the time-step.
This is easily visualized by choosing two different time-steps:

::timestep equal 0.001 <math>\implies</math> n_steps = 100000
::timestep equal 0.002 <math>\implies</math> n_steps = 50000

:The two time-steps have different number of steps but the sum of the time-steps is constant and its 100.
:In other words, the number of <math>\delta t</math>'s in the Velocity Verlet algorithm is different but their sum is constant.
:This sum of the time-steps is not dependent on the actual duration of the simulation.

===Task 13: Variation of Energy, Pressure and Temperature plots with the time-step.===
[[File:Stabilisation_at_0.001.png]]
[[File:Energy stabilistion timestep.png]]

It can be seen from Image that for a 0.001 times-step the simulation reaches equilibrium at around 0.02. Although pressure and temperature fluctuate throughout the simulation the energy remains constant with a standard deviation of___ that was previously calculated.This result was therefore expected.

If we know compare the effect of the time-step on the energy the following conclusions can be drawn:

* The larger the time-step, the less accurate the Velocity-Verlet algorithm, choosign a larger time-step causes the error to increase
* Both 0.001 and 0.0025 give accurate energy readings but choosing 0.0025 will allow the simulations to run for longer.
* Above 0.003 The standard deviation of the energy increases to
* 0.015 as a time step does not reach an energy stabilisation. 

In view of the above, the simulations will be run for both 0.001 and 0.0025 as time-steps.

===Task 14: Determining <math>\gamma</math>===
:<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T</math>
To achieve <math>T = \mathfrak{T}</math> every <math> v_i </math> has to be multiplied by <math>\gamma</math>

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

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

Cancelling and rearranging leads us to:

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

=== Task 15: Numbers 100 1000 100000 ===

*'''First''': Use the input values '''every''' this number of time steps.
*'''Second''': number of times to '''repeat''' the input values for calculating the averages.
*'''Third''': To calculate averages for this '''frequency''' of timesteps

===Task 16: Editing the equations of state===
[[File:0.0025uploadIdGasANA.png|caption|Variation of density with temperature at different pressures (blue). The orange line represents the behaviour of a perfect gas]]

The graphs above show the change in density with respect to temperature. These simulations of ensembles were run at a time-step of 0.0025 that as shown in Task 1 ensures that the energy of the system does not fluctuate more than 1%. From these graphs several conclusions can be drawn:

- An increase in temperature causes a decrease in density.

- An increase in pressure causes an increase in density.

===Task 17: Heat Capacities===
[[File:Heat_capacity_with_T_analo.png]]

For this simulation a file was created [[File:Log_file.rtf|input file]]

From the plot above it can be seen that a temperature increase causes a decrease in heat capacity.

The reasoning behind this trend is that because the simulation establishes a certain lattice structure there is a constraint on the positions of the atoms. This constraint stops the atoms from gaining free energy. At constant volume the heat capacity can be written as <math> C_V = (\frac{dU}{dT})_V</math> <ref name="definition"/> . Consequently if the increase in internal energy is not greater that the temperature increase the heat capacity will decrease.

A change in density changes the volume of the lattice cells, as calculated in Task . At higher densities the heat capacity of the system increases. This is due to the lattice cells being more effectively filled and hence more atoms per area can potentially increase their internal energy.

===Task 18: Radial distribution functions===

[[File:RDf+intRDFal3214analo.png|thumb|500px|]]

[[File:Nearest_neighboursal32143214.png|thumb|500px|]]

===Task 19: Mean Squared Displacement===
[[File:Allal3214.png|thumb|900px|]]

The MSD was plotted as seen on the images on the right. Fronm those plots, the gradient was calculated and divided by 6 to give D, the diffusion coefficient. the results are summarized in the following table:
{| class="wikitable"
|+ D from MSD plots.
!
!Solid
!Liquid
!Gas
|-
|Millions
|8.33E-9
|1.66E-4
|0.005
|-
|fcc
|1.16E-9
|0.0011
|0.00303
|-
|sq lat
|3.3E-6
|0.0011
|0.00302
|-
|}

There are two trends than can be inferred from the MDS plots.
*The diffusion coefficient increases when going from a solid phase to a liquid or a gas.
*The greater the density of the lattice, fcc with 4 atoms per simulation box, the lower the diffusion coefficient.

Under the same density conditions of the ensemble, the intrinsic density of an square lattice with 1 atom per simulation box, causes it to have an MSD plot which resembles after a liquid more than after a solid.

===References===

<ref name="definition"/>Peter Atkins, Julio de Paula, Atkins Physical Chemistry, Oxford 2008.