View source for Rep:DSB113

== Completing Harmonic Oscillator ==
The Velocity-Verlet algorithm is used to model the behaviour of a classical harmonic oscillator. 

Energy in this case was just kinetic energy not total energy. Potential energy was added to this using the formula  Φ= 0.5kx<sup>2,</sup> where k is the force constant and x is distance. [[File:HO.PNG|options|file=]]

=== '''Calculation Parameters''' ===
{| class="wikitable"
|Φ
|0.00
|-
|A
|1.00
|-
|ω
|1.00
|-
|mass
|1.00
|-
|k
|1.00
|-
|original timestep
|0.1
|-
|updated timestep
|0.2
|}
The position of maxima in error were estimated and plotted against time generating the formula; y = <span lang="EN-US">0.0004x- 7E-05. The timestep was varied to provide an error within 1%. The timestep used was 0.2, as increasing the timestep increased the error. The original timestep used through the Harmonic Oscillator model was 0.1. This reduced the error. It is important to monitor the total energy of a physical system to ensure energy is conserved and to see any changes. </span>

[[File:analytical01.PNG|options|file=]]
 Graph showing the Analytical solution against timestep

[[File:KE01.PNG|options|file=]]
 Graph showing the Kinetic Energy against timestep
[[File:error01.PNG|options|file=]]
 Graph showing the error against timestep

[[File:errormax.PNG|options|file=]]
 Graph showing the maximum error against timestep

[[File:totenergy.PNG|options|file=]]
 Graph showing the total energy against timestep

The graphs above are for a timestep of 0.1. The graphs below show an increase in timestep to 0.2 providing an error within 1%. 

[[File:analytical02.PNG|options|file=]]
 Graph showing the Analytical solution against timestep when timestep has been doubled from 0.1 to 0.2

[[File:KE02.PNG|options|file=]]
 Graph showing the Kinetic Energy against timestep when timestep has been doubled from 0.1 to 0.2
[[File:error02.PNG|options|file=]]
 Graph showing the error against timestep when timestep has been doubled from 0.1 to 0.2 to ensure error is only 1%

=== Lennard-Jones Interaction ===
[[File:LJ interaction phi.PNG|link=|left]]When the potential energy is zero, the formula rearranges to show that r<sub>0</sub>=σ.

The force at this separation is calculated by differentiating [[File:U.PNG|options|file=]] to give [[File:Force.PNG|400x90px|file=]]

[[File:Equi.PNG|file=|525x525px]]

=== Periodic Boundary Conditions ===
Realistic volumes of liquid cannot be simulated due to the number of water particles involved as shown below. 

The number of water molecules in 1ml of water under standard conditions is calculated using density of water and Avagadros Number. 1ml=1g, 1/18 (Molar mass of water) x 6.022 x 10<sup>23 </sup>=3.34 x 10<sup>22</sup>.

The volume of 10000 water molecules is calculated by doing the reverse of above. 10000/(6.022 x 10<sup>23</sup>) = 1.66 x 10<sup>-20</sup>. This is the number of Moles of 10000 water molecules. Multiply this by 18 to find the mass of 10000 water molecules which is 3 x 10<sup>-19</sup> g converted to cm<sup>3 </sup>is the same as density of water is 1 g/cm<sup>3</sup>.

If an atom at position (0.5, 0.5, 0.5) moves along vector (0.7, 0.6, 0.2) it ends up at (0.2, 0.1, 0.7) after periodic boundary conditions have been applied.

The Lennard-Jones parameter for argon are σ=0.34nm ε/K<sub>B</sub>=120K.

If the cutoff is 3.2, r is 1.088 in real units as r*= r/σ.

The well depth is 1.656 x 10<sup>-21</sup> as ε/K<sub>B</sub>=120K where ε is the well depth and K<sub>B </sub>is the Boltzmann constant. 

Similarly if the reduced temperature (T*) is 1.5, T is 180K in real units as T* = (k<sub>B</sub>T)/ ε

== Creating the Simulation Box ==
As there is no order a single point of reference cannot be used to find positions of all other atoms in liquid and vapour phases. 

Giving atoms random starting coordinates causes problems in simulations as atoms may land on top of each other and repel causing minor explosions. 

If the distance between points of a lattice are 1.07722 in reduced units, the number density of a 0.8

<span lang="EN">1 lattice point. </span>

<span lang="EN">For a face-centred cubic lattice, with a number density of 1.2, the side length of the cubic unit cell is 1.36 in reduced units. Through the create_atoms command, 3000 atoms would be created if this lattice had been defined.</span>

=== <span lang="EN">Setting the Properties of the''' '''atoms'''<nowiki/>''' ===
<span lang="EN">The purpose of:</span>
 <span lang="EN">mass 1 1.0</span>
<span lang="EN">is to define the mass of atom type 1 which is 1.0</span>
 <span lang="EN">pair_style lj/cut 3.0</span>
<span lang="EN">is defining the Lennard-Jones cutoff as 3.0 with no coulomb interaction.</span>
 <span lang="EN">pair_coeff ** 1.0 1.0</span>
<span lang="EN">defines that for each type of atom and each other type of atom, the interaction remains the same at 1.0.</span>

<span lang="EN">Since both x<sub>i</sub> (0) and v<sub>i'' ''</sub>(0) are being specified, the Velocity-Verlet algorithm will be used.</span>

<span lang="EN">The reason just </span>
 <span lang="EN">timestep 0.001</span>

 <span lang="EN">run 100000</span>
<span lang="EN">is not written is due to the fact that the timestep being used numerous times. Using the other method makes it easier to change timesteps and hence easier to compare the results. </span>

=== '''<span lang="EN">Checking Equilibration </span>''' ===
It is of paramount importance that the system reaches equilibrium state.

[[File:totenergy1.PNG|options|file=]]
 Graph showing the Total Energy against timestep
[[File:pressure.PNG|options|file=]]
 Graph showing the Pressure against timestep
[[File:temp.PNG|options|file=]]
 Graph showing the Temperature against timestep
[[File:totents.PNG|options|file=]]
 Graph showing the Total Energy against timestep at individual timesteps
<span lang="EN">The simulation reaches equilibrium at 0.41. A timestep of 0.01 is the largest to give acceptable results. The 0.015 timestep is a particularly bad choice as the simulation doesn't reach equilibrium. The energy of this time
step decreases initially before increasing rapidly. </span>

== '''<span lang="EN">Temperature and Pressure Control</span>''' ==
The systme upto this point was run under a microconical ensemble. Now the system will be run under isobaric-isothermal ensemble. <span lang="EN">The timestep used for this is 0.0075</span>

{| class="MsoNormalTable"
 
  |
<nowiki>  </nowiki>'''Temperature'''
  
  |
<nowiki>  </nowiki>'''Pressure 1'''
  
  |
<nowiki>  </nowiki>'''Pressure 2'''  
 
 |-
  |
<nowiki>  </nowiki>1.6
  
  |
<nowiki>  </nowiki>2.6
  
  |
<nowiki>  </nowiki>3.5  
 
 |-
  |
<nowiki>  </nowiki>2.0
  
  |
<nowiki>  </nowiki>2.6
  
  |
<nowiki>  </nowiki>3.5  
 
 |-
  |
<nowiki>  </nowiki>2.4
  
  |
<nowiki>  </nowiki>2.6
  
  |
<nowiki>  </nowiki>3.5  
 
 |-
  |
<nowiki>  </nowiki>2.8
  
  |
<nowiki>  </nowiki>2.6
  
  |
<nowiki>  </nowiki>3.5  
 
 |-
  |
<nowiki>  </nowiki>3.2
  
  |
<nowiki>  </nowiki>2.6
  
  |
<nowiki>  </nowiki>3.5  
 
|}

As T fluctuates, a constant factor is used to ensure the target temperature is achieved. If temperature is higher than target, the constant is less than 1 however, if the target is higher than the temperature the constant is greater than 1. U<span lang="EN">sing these 2 equations:</span>

[[File:Thermostat.PNG|file=|281x281px]]

by dividing the first equation from the second equation, it can be seen that <math>\gamma = \sqrt{\left (\frac{\mathfrak{T}}{T}  \right )}</math>

If the <em>ave</em> setting is <em>window</em>, then the values produced on timesteps that are multiples of <em>Nfreq</em> are summed and averaged within a moving “window” of time, so that the last M values are used to produce the output. E.g. if M = 3 and Nfreq = 1000, then the output on step 10000 will be the average of the individual values on steps 8000,9000,10000. Outputs on early steps will average over less than M values if they are not available.

Values of temperature will be sampled every 100 units for the average. 1000 measurements contribute to the average. Since the timestep used in this case is 0.0075 and there are 1000 measurements, the time simulated will be 7.5

=== Plotting Equations of State ===

[[File:DVT26.PNG|file=]]
 Graph showing density against temperature at pressure of 2.6 with error bars
[[File:IGL26.PNG|file=]]
 Graph showing density against temperature using the ideal gas law at pressure of 2.6
[[File:DVT35.PNG|file=]]
 Graph showing density against temperature at pressure of 3.5 with error bars
[[File:IGL35.PNG|file=]]
 Graph showing density against temperature using the ideal gas law at pressure of 3.5
The simulated density using the ideal gas law simulates higher densities at both pressures due to assumptions made under the law. The ideal gas law doesn't account for molecular size or intermolecular interactions. As temperature increases, the effect of intermolecular interactions decreases. As pressure increases the discrepancy increases.

== Calculating Heat Capacity ==
In thermodynamics, many uantities can be calculated by considering how much the system can fluctuate from its average equilibrium state.

[[File:heatcap.PNG|options|file=]]
 Graph showing Heat capacity/volume against temperature at density of 0.2 and 0.8

A copy of an example input script can be found here [[Media:d8t28.ogg]]. This is an example when density is 0.8 and Temperature is 2.8. 

At both densities, the heat capacity/volume decreases as temperature increases.  Lower density has lower heat capacity as heat capacity is an intensive property.

== Radial Distribution Function ==
The RDF provides the distance between an atom and its nearest neighbour, second nearest, etc. It can be found experimentally hence ensures the simulation is reproducing correctly. The atomic trajectory is recorded to generate the RDF.

[[File:g(r).PNG|options|file=]]
 Graph showing the Radial Distribution Function denoted by g(r) of the 3 phases against timestep
[[File:Ig(r).PNG|options|file=]]
 Graph showing the integrated radial distribution function of the 3 phases against timestep
[[File:Lattice1.PNG|options|file=]]
 FCC lattice representation
The RDF shows the distance between x from the centre of the atom and the probability. A gas phase has very few neighboring atoms as the RDF shows the fewest peaks while solid has the most neighboring atoms and hence the most peaks. The first 3 peaks correspond to the 4 atoms closest to the atom being observed, followed by the atom in a similar position but 1 up and finally another 4 atoms which are 1.5 up. Peaks 1 and 3 correspond to 4 atoms as the atoms are equidistant from the atom being observed.  The lattice spacing is 0.8 in reduced units calculated by finding the difference between 2 maximas which correspond to 0.5 of the hypotenuse.The coordination number for each of the first 3 peaks is 12, 18 and 42 respectively. This is found by looking at the first 3 minimas after a peak on g(r) and finding the equivalent timestep and hence the integrated g(r) value.

== Mean Squared Displacement ==
The MSD provides information on how much atoms in the simulation move around, characterised by the diffucsion coefficient. 

[[File:LiqMSD.PNG|options|file=]]
 Graph showing the mean squared displacement against timestep for a liquid
[[File:SolidMSD.PNG|options|file=]]
 Graph showing the mean squared displacement against timestep for a solid
[[File:GasMSD.PNG|options|file=]]
 Graph showing the mean squared displacement against timestep for a vapour
[[File:LiqMSD1M.PNG|options|file=]]
 Graph showing the mean squared displacement against timestep for a liquid for a million particles
[[File:SolidMSD1M.PNG|options|file=]]
 Graph showing the mean squared displacement against timestep for a solid for a million particles
[[File:GasMSD1M.PNG|options|file=]]
 Graph showing the mean squared displacement against timestep for a vapour for a million particles
The estimated D for liquid is 0.0849 reduced units, for solid its 0.0000005 reduced units and for vapour, its 0.39215 reduced units. The general trend observed is what is expected. A solid is tightly packed so will not move much however, a vapour is loosely packed and hence has the highest D. Increasing the number of particles to a million has very little effect on the Solid and liquid where D remains near the same, 0.000005 and 0.08726. However there is a large difference for the vapour phase which has D of 2.933 reduced units.

== Velocity Autocorrelation Function ==
The diffusion coefficient can also be calculated using VACF which provides the difference in velocity between 2 different times. Velocities should not be correlated at long times indicating velocity is random. 

[[File:integration.PNG|file=|849x849px]]

[[File:Harmonic.PNG|options|file=]]
 Graph showing the 1D harmonic oscillator function of the Velocity Autocorrelation Function denoted by C(t)
[[File:LJVACF.PNG|options|file=]]
 Lennard-Jones liquid and solid harmonic oscillator VACF.
The minima in the VACFs for liquid and solid system represents the repulsion between atoms as the repulsive forces are stronger when the atoms are closer together. At the maxima, the atoms are moving in the same direction as there is a strong correlation. The difference between the liquid and solid is that liquid particles are able to move more compared to solid particles, this dampens the VACF as liquid particles have weaker inter-particle interactions. The harmonic oscillator is different to the Lennard Jones solid and liquid as the harmonic oscillator is not driven or damped. LJ solid and liquid account for friction/damping which slows the motion of the system by decreasing the velocity.
{| class="wikitable"
!Phase
!Running Integral
!Diffusion Coefficient
|-
|Gas
|3.382857
|1.127619
|-
|Liquid
|0.587281
|0.19576
|-
|Solid
|<nowiki>-0.00111</nowiki>
|<nowiki>-0.00037</nowiki>
|}

[[File:ILiq.PNG|options|file=]]
 Integral for a liquid against timestep
[[File:RILiq.PNG|options|file=]]
 Running Integral for a liquid against timestep
[[File:ISol.PNG|options|file=]]
 Integral for a solid against timestep
[[File:RISol.PNG|options|file=]]
 Running Integral for a solid against timestep
[[File:IVap.PNG|options|file=]]
 Integral for a vapour against timestep
[[File:RIVap.PNG|options|file=]]
 Running Integral for a vapour against timestep
The Integral was calculated using the trapezium rule. The largest source of error is decreasing the size of each trapezium to increase the accuracy.