ChemWiki - User contributions [en]

File:Integral of VACF of 1 million solid .png

2016-02-26T00:16:17Z

Nhl13: Nhl13 uploaded a new version of File:Integral of VACF of 1 million solid .png

File:Integral of VACF of 1 million liquid.png

2016-02-26T00:15:51Z

Nhl13: Nhl13 uploaded a new version of File:Integral of VACF of 1 million liquid.png

File:Integral of VACF of 1 million gas.png

2016-02-26T00:15:27Z

Nhl13: Nhl13 uploaded a new version of File:Integral of VACF of 1 million gas.png

File:Integral of VACF of solid.png

2016-02-26T00:15:06Z

Nhl13: Nhl13 uploaded a new version of File:Integral of VACF of solid.png

File:Integral of VACF of liquid.png

2016-02-26T00:14:43Z

Nhl13: Nhl13 uploaded a new version of File:Integral of VACF of liquid.png

2016-02-25T21:50:19Z

Nhl13: /* Equilibration */

==='''Simulation of a simple liquid'''===

==Introduction to Molecular Dynamics Simulation==

'''<big>TASK 1</big>: Open the file HO.xls. In it, the velocity-Verlet algorithm is used to model the behaviour of a classical harmonic oscillator. Complete the three columns "ANALYTICAL", "ERROR", and "ENERGY": "ANALYTICAL" should contain the value of the classical solution for the position at time <math>t</math>, "ERROR" should contain the ''absolute'' difference between "ANALYTICAL" and the velocity-Verlet solution (i.e. ERROR should always be positive -- make sure you leave the half step rows blank!), and "ENERGY" should contain the total energy of the oscillator for the velocity-Verlet solution. Remember that the position of a classical harmonic oscillator is given by <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math> (the values of <math>A</math>, <math>\omega</math>, and <math>\phi</math> are worked out for you in the sheet).'''

[[File:position against time.PNG|thumbnail|thumb|350px|right|Fig.1 Analytical position vs. Time.]]

[[File:Energy against time.png|thumbnail|thumb|350px|right|Fig.2 Total Energy from velocity-Verlet solution vs. Time.]]

[[File:Error against time.PNG|thumbnail|thumb|350px|right|Fig.3 Error vs. Time.]]

Analytical position of a classical harmonic oscillator was plotted in figure 1, where the analytical position was calculated by equation:

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

where <math>A=1, \omega = 1, \phi=0</math>

The total energy of the classical harmonic oscillator for the velocity-Verlet solution was calculated by equation:

<math>Total Energy = Kinetic Energy + Potential Energy</math>,

For a classical harmonic oscillator, <math>KE=\frac{1}{2}\ mv^2</math> and <math>PE=\frac{1}{2}\ kx^2</math>,

and x(t) from velocity-Verlet solution was used instead of analytical x(t).

The total energy of the oscillator was plotted in figure 2. The energy fluctuated between 0.500 and 0.49875, which was relatively small when compared to the maximum total energy, i.e. 0.25% of the total energy.

The error in position between classical solution and velocity-Verlet solution was plotted in figure 3. It showed that the error gradually increased with time.

'''<big>TASK 2</big>: For the default timestep value, 0.1, estimate the positions of the maxima in the ERROR column as a function of time. Make a plot showing these values as a function of time, and fit an appropriate function to the data.'''

Figure 4 showed the maxima of error in position between the classical and velocity-Verlet solution respectively was plotted against time. It gave a straight line with a positive slope that indicated the error was directly proportional to time. The corresponding function was y = 0.0004x, when the line was set to intercept at origin.

[[File:Position of error maxima.PNG|thumbnail|thumb|350px|right|Fig 4.Position of error maxima vs. Time.]]

'''<big>TASK 3</big>: Experiment with different values of the timestep. What sort of a timestep do you need to use to ensure that the total energy does not change by more than 1% over the course of your "simulation"? Why do you think it is important to monitor the total energy of a physical system when modelling its behaviour numerically?'''

[[File:Different timestep.PNG|thumbnail|thumb|350px|right|Fig.5 Fluctuation in energy with different timestep.]]

The maximum total energy of the harmonic oscillator was 0.500. From figure 5, the degree of fluctuation increased when time step increased. The total energy had a fluctuation of 5x10^-3, i.e. 1% of the maximum energy, when the timestep reached 0.20 s. Hence, the timestep had to be equal or smaller than 0.20 s in order to obtain an accurate simulation.

'''<big>TASK 4.1:</big> For a single Lennard-Jones interaction, <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>, find the separation, <math>r_0</math>, at which the potential energy is zero.'''

When potential energy <math>\phi = 0</math>, <math>r=r_{0}</math>

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

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

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

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

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

'''<big>TASK 4.2</big>: What is the force at this separation?

F (r) = - dU/dr

F (r) = - d [4 ε (σ12 / r12) + 4 ε (σ6 / r06)] / dr

F (r) = d [ - 4 ε (σ12 / r12) + 4 ε (σ6 / r06)] / dr

F (r) = d [ - 4 ε (σ12 r-12) + 4 ε (σ6 r-6)] / dr

F (r) = 48 ε σ12 r-13 - 24 ε σ6 r-7

When r = r0, r0 = σ

F (r0) = 48 ε σ12 σ-13 - 24 ε σ6 σ-7

F (r0) = 48 ε σ-1 - 24 ε σ-1

F (r0) = 24 ε / σ

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

When F = 0, r = req

F (r) = 48 ε σ12 r-13 - 24 ε σ6 r-7

0 = 48 ε σ12 req-13 - 24 ε σ6 req-7

48 ε σ12 req-13 = 24 ε σ6 req-7

2 σ12 / σ6 = req13 / req7

2 σ6 = req6

req = 21/6 σ

When req = 21/6 σ,

Φ = 4 ε [(σ12 /r12 ) - (σ6 / r6)]

Φ = 4 ε {[σ12 /(21/6 σ)12] - [σ6 / (21/6 σ)6]}

Φ = 4 ε [(1/4) - (1/2)]

Φ = - ε

'''<big>TASK 4.4</big>: Evaluate the integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math> when <math>\sigma = \epsilon = 1.0</math>.

When <math>\sigma=\epsilon=1.0</math>

= ʃ 4 ε [(σ12 / r12 ) - (σ6 / r6)] dr

= [-4 / 11 r11 + 4 / 5 r5]

Therefore,

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

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

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

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

Density of water = 1 g cm-3

Mass of water in 1 ml of water = 1 g

No. of moles in 1 ml of water = 1 g ÷ 18.0 g mol-1 = 0.0556 mol

No. of water molecules in 1 ml of water = 0.0556 mol × 6.022×1023 = 3.34×1022 molecules

No. of moles of 10000 water molecules = 10000 ÷ 6.022×1022 = 1.66×10-20 mol

Mass of 10000 water molecules = 1.66×10-20 mol × 18.0 g mol-1 = 2.99×10-19 g

volume of 10000 water molecules = 2.99×10-19 g ÷ 1 g cm-1 = 2.99×10-19 ml

'''<big>TASK 6</big>: Consider an atom at position <math>\left(0.5, 0.5, 0.5\right)</math> in a cubic simulation box which runs from <math>\left(0, 0, 0\right)</math> to <math>\left(1, 1, 1\right)</math>. In a single timestep, it moves along the vector <math>\left(0.7, 0.6, 0.2\right)</math>. At what point does it end up, ''after the periodic boundary conditions have been applied''?'''.

Periodic boundary condition states that boxes are repeated infinitely in all direction, hence, if one molecule exits the boundary of the original box, one molecule from the replicated box enters the original box through the opposite face to retain the number of molecules. As a result, we could work out the final position of the atom simply by adding the movement vector to its initial position vector, i.e. (0.5, 0.5, 0.5) + (0.7, 0.6, 0.2) = (1.2, 1.1, 0.7). As the box had a dimension of (1.0, 1.0, 1.0), so the resulting position vector was (0.2, 0.1, 0.7).

'''<big>TASK 7.1</big>: The Lennard-Jones parameters for argon are <math>\sigma = 0.34\mathrm{nm}, \epsilon\ /\ k_B= 120 \mathrm{K}</math>. If the LJ cutoff is <math>r^* = 3.2</math>, what is it in real units? '''

r* = r / σ

r = r* σ

r = (3.2)(0.34×10-9)

r = 1.088 nm

'''<big>TASK 7.2</big>: What is the well depth in <math>\mathrm{kJ\ mol}^{-1}</math>? '''

Φ = - ε, as calculated above, and ε / kB = 120K

Φ = - 120 × kB

Φ = - 1.66×10-21 J ÷ 1000 × NA

Φ = 0.998 kJmol-1

'''<big>TASK 7.3</big>: What is the reduced temperature <math>T^* = 1.5</math> in real units?'''

T* = kBT ÷ ε

T = (T*)(ε) ÷ kB

T = (1.5)(120)(kB) ÷ kB

T = (1.5)(120)

T = 180 K

==Equilibration==

'''<big>TASK 8</big>: Why do you think giving atoms random starting coordinates causes problems in simulations? Hint: what happens if two atoms happen to be generated close together?'''

In liquid, molecules are held together by intermolecular bonds. They are not fixed and can flow around and have a certain average intermolecular distance. If the molecules are too far apart, the attractive interaction between them is not as strong as it used to and behaves as gas instead of liquid. In contrast, if molecules are too close together, they experience a repulsion force between them. The liquid simulation would be affected in both extremes.

'''<big>TASK 9</big>: Satisfy yourself that this lattice spacing corresponds to a number density of lattice points of <math>0.8</math>. Consider instead a face-centred cubic lattice with a lattice point number density of 1.2. What is the side length of the cubic unit cell?'''

<math>number density = \frac{number of atom}{volume}</math> can be used to rationalise the relationship between the number density, and the lattice spacing, a.

Case 1: In a simple cubic lattice, there is only 1 atom inside a lattice unit cell.

number density = number of atom / volume

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

<math>a=1.0772 m</math>

This result agreed with the value obtained from LAMMPS.

Case 2: In a face-centred cubic lattice, there are 1 atom at each corner and 1 atom on each face. Atom at the corner and on the face contribute to 1/4 and 1/2 atom to the lattice. Hence, number of atom in a fcc lattice = 1/8 × 8 + 1/2 × 6 = 4

number density = number of atom / volume

1.2 = 4 / a3

a = 1.494 m

Therefore, a fcc lattice with a number density of 1.2 will have a side length of 1.494 m.

'''<big>TASK 10</big>: Consider again the face-centred cubic lattice from the previous task. How many atoms would be created by the create_atoms command if you had defined that lattice instead?'''

Each fcc unit cell consists of 4 atoms. With the same 10×10×10 box, we are having 1000 unit cells, hence, fcc lattice will consist of 4×1000 = 4000 atoms.

'''<big>TASK 11</big>: Using the [http://lammps.sandia.gov/doc/Section_commands.html#cmd_5 LAMMPS manual], find the purpose of the following commands in the input script:'''

<pre>mass 1 1.0</pre>Syntax of the command<pre>mass I value</pre>Mass command defines the mass of atoms, where I is the number of atom type and value is the mass of atom.

In our case, there is only 1 type of atom and its mass is 1.0 in reduced unit.<pre>pair_style lj/cut 3.0</pre>Syntax of the command<pre>pair_style style args</pre>pair_style command compute different styles of pairwise interaction. The second <nowiki>''style''</nowiki> defines a particular type of pairwise interaction, and the <nowiki>''args''</nowiki> is the arguments used by that particular style. In our case, lj/cut is used to compute the cutoff Lennard-Jones potential with no Coulomb, and the cutoff potential is set at 3.0 reduced unit.<pre>pair_coeff * * 1.0 1.0</pre>Syntax of the command

pair_coeff command specifies the pairwise force field coefficients for one or more pairs of atom types. <nowiki>''I,J'' is the atom types and ''args''</nowiki> is the coefficients for one or more pairs of atom types. In out case, I and J equal * which means there is N types of atom and means all types from 1 to N, and both coefficients are set at 1.0.

'''<big>TASK 12</big>: Given that we are specifying <math>\mathbf{x}_i\left(0\right)</math> and <math>\mathbf{v}_i\left(0\right)</math>, which integration algorithm are we going to use?'''

We are using Velocity Verlet Algorithm, see below.

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

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

### RUN SIMULATION ###
run ${n_steps}
run 100000
</pre>
'''The second line (starting "variable timestep...") tells LAMMPS that if it encounters the text ${timestep} on a subsequent line, it should replace it by the value given. In this case, the value ${timestep} is always replaced by 0.001. In light of this, what do you think the purpose of these lines is? Why not just write:'''
<pre>
timestep 0.001
run 100000</pre>
Those extra lines would allow the variable to change automatically when timestep in the second line is changed.

'''<big>TASK 14</big>: make plots of the energy, temperature, and pressure, against time for the 0.001 timestep experiment (attach a picture to your report). Does the simulation reach equilibrium? How long does this take? When you have done this, make a single plot which shows the energy versus time for all of the timesteps (again, attach a picture to your report). Choosing a timestep is a balancing act: the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality; short timesteps, however, mean that the same number of simulation steps cover a shorter amount of actual time, and this is very unhelpful if the process you want to study requires observation over a long time. Of the five timesteps that you used, which is the largest to give acceptable results? Which one of the five is a ''particularly'' bad choice? Why?'''[[File:total energy equilibration.PNG|thumbnail|thumb|350px|right|Fig6. Equilibration of total energy.]][[File:temperature equilibration.PNG|thumbnail|thumb|350px|right|Fig.7 Equilibration of Temperature.]][[File:pressure equilibration.PNG|thumbnail|thumb|350px|right|Fig.8 Equilibration of Pressure.]][[File:total energy at different timestep.PNG|thumbnail|thumb|350px|right|Fig.9 Total energy at different timestep.]]

All three properties reached equilibrium after short period of time when timestep was 0.001 s. The total energy, temperature and pressure reached the equilibrium after approximately 0.42, 0.32 and 0.25 s respectively. When timestep was as large as 0.015 s, the total energy did not reach equilibrium. This was particularly bad for simulation. When the timestep was equal to or lower than 0.01 s, the energy did reach equilibrium. However, the fluctuation of the energy level was quite large when timestep equal 0.01 and 0.075 s. Also, the timestep cannot be too short, as it will lower the efficiency of the simulation. As a result, the most acceptable timestep for simulation to compensate both disadvantages was 0.0025 s.

==Running simulations under specific conditions==

'''<big>TASK 15</big>: Choose 5 temperatures (above the critical temperature <math>T^* = 1.5</math>), and two pressures (you can get a good idea of what a reasonable pressure is in Lennard-Jones units by looking at the average pressure of your simulations from the last section). This gives ten phase points — five temperatures at each pressure. Create 10 copies of npt.in, and modify each to run a simulation at one of your chosen <math>\left(p, T\right)</math> points. You should be able to use the results of the previous section to choose a timestep. Submit these ten jobs to the HPC portal. While you wait for them to finish, you should read the next section.'''

In our simulation, timestep of 0.0025 s was used, as it allowed the oscillator to reach equilibrium and had a small fluctuation in energy level. The two pressures used were 2.55 and 2.65 atm, as there were inside the range of the pressure equilibrium from the last section. Finally, the temperature used were 1.7, 1.9, 2.1, 2.3 and 2.5, where all of them were above the critical temperature.

'''<big>TASK 16</big>: We need to choose <math>\gamma</math> so that the temperature is correct <math>T = \mathfrak{T}</math> if we multiply every velocity <math>\gamma</math>. We can write two equations:'''

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

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

'''Solve these to determine <math>\gamma</math>.'''

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

<math>\frac{\sum_i m_i v_i^2}{3T} = N k_B</math>

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

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

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

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

'''<big>TASK 17</big>: Use the [http://lammps.sandia.gov/doc/fix_ave_time.html manual page] to find out the importance of the three numbers ''100 1000 100000''. How often will values of the temperature, etc., be sampled for the average? How many measurements contribute to the average? Looking to the following line, how much time will you simulate?'''

100 is the Nevery, 1000 is the Nrepeat and 100000 is the Nfreq. Nevery defines how often are the input values is used for calculation, Nrepeat defines the number of times that the input values are used to calculate the average, and Nfreq defines the number of timestep for calculating an average. In our case, the input values are used every 100 timesteps, average calculated from input values are repeated for 1000 times, and the averages are calculated every 100000 timesteps. We used timestep of 0.0025 s, so our simulation time is 100000 x 0.0025 = 250 (reduced unit)

'''<big>TASK 18</big>: When your simulations have finished, download the log files as before. At the end of the log file, LAMMPS will output the values and errors for the pressure, temperature, and density <math>\left(\frac{N}{V}\right)</math>. Use software of your choice to plot the density as a function of temperature for both of the pressures that you simulated. Your graph(s) should include error bars in both the x and y directions. You should also include a line corresponding to the density predicted by the ideal gas law at that pressure. Is your simulated density lower or higher? Justify this. Does the discrepancy increase or decrease with pressure?'''

[[File:density vs temp.PNG|thumbnail|thumb|350px|right|Fig.10 Density vs. Temp at different pressure.]]

Figure 10 showed that the simulated density decreases with temperature and increases with pressure respectively. This could be rationalised easily. As temperature increases, the molecules have higher kinetic energy and overcome the attractive interaction between molecules, thus a lower density. On the other hand, as pressure increases, the molecules are experiencing a greater force and are compressed together, thus a higher density. The error in density was much smaller than the error in temperature, shown in figure 10. The simulated density was much lower than the density calculated from the ideal gas law. The ideal gas law does not take into account of the repulsive force between molecules, while the simulation was able to mimic the repulsive force. Therefore, the simulated density was lower than the calculated density as repulsive force between molecules pushed each other away. The discrepancy between simulated and calculated density increases with pressure. This could be rationalised by knowing that the gas behaves more as a real gas instead of ideal gas when pressure increase, hence the discrepancy increases.

==Calculating heat capacities using statistical physics==

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

[[File:heat capacity against temperature.png|350px|right|thumb|fig.11 heat capacity against temperature at different density]]

Cv/V was plotted against temperature in figure 11. For both density, the simulation indicated the volume of the liquid remain unchanged, i.e. 16875 and 4218.75. The graph showed that the Cv/V decrease with temperature. The extent of decrease was greater when temperature was small and started to decrease slower as temperature increase gradually. However, there was an anomalous when T=2.4, D=0.8, which possibly arose due Brownian motion of the liquid.

Example of the script is shown below, where density is 0.2 and temperature is 2.0.

<pre>

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

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

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

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

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

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

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

### SPECIFY TIMESTEP ###

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

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

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp press density atoms vol
variable dens equal density
variable dens2 equal density*density
variable temp equal temp
variable temp2 equal temp*temp
variable press equal press
variable press2 equal press*press
variable N equal atoms
variable N2 equal atoms*atoms
variable etotal equal etotal
variable etotal2 equal etotal*etotal
variable volume equal vol

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

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

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

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

</pre>

==Structural properties and the radial distribution function==

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

The following table shows the conditions under which the simulations of the Lennard-Jones system were performed.
{| class="wikitable"
|-
! !! Temperature !! Density
|-
|Solid ||1.25 || 1.25
|-
| Liquid || 1.1 || 0.7
|-
|Gas || 1.1 || 0.05
|}

[[File:RDF for solid, liquid and gas.PNG|thumbnail|thumb|350px|right|Fig.12 RDF for solid, liquid and gas.]]

[[File:Integral of RDF for solid, liquid and gas.PNG|thumbnail|thumb|350px|right|Fig.13 Integral of RDF for solid, liquid and gas.]]

RDF provides the density of atoms at particular radius to a reference atom. Each peak in RDF indicates a particularly favoured separation distance for the neighbouring atoms to a reference atom. Thus, RDF reveals the atomic structure of the system being simulated. Figure 12 showed that RDF of solid had lots of sharp peaks, while there were only a few broad peaks in liquid, and just one broad peak for gas. It reveals atoms in solid are highly ordered and stays at particular distance to the reference atom. Sharp peak also indicates the atoms are not spread out randomly but are concentrated at favoured particular distances. The coordination number at favoured particular distance can be determined from the integral of RDF. In contrast, a few peaks in liquid RDF shows its system is less ordered and atoms are able to move around freely and rarely stay at a particular distance to the reference atom. Peaks still shows up illustrates that the atoms are still held together by some sort of attraction forces. On the other hand, broad peaks are observed in RDF of liquid and gas, where peak for gas is broader than those for liquid. This can be rationalised by the Brownian motion of the atoms in the system. The atoms in gas have higher energy, therefore, can move more randomly and resulted in a broad distribution at particular distance, where the only peak observed indicates the gas atoms are still interacting with each other.

[[File:1st and 2nd coordinating neighbours of fcc.PNG|thumbnail|thumb|350px|right|Fig.14 1st and 2nd coordinating neighbours of fcc.]]

[[File:3rd coordinating neighbours of fcc.PNG|thumbnail|thumb|350px|right|Fig.15 3rd coordinating neighbours of fcc.]]

The first three peaks at the RDF of soild were at a, b and c. These corresponded to the first three nearest coordinating neighbours to a reference atom. The lattice spacing between the reference atom and the nearest neighbours can be calculated by Pythagoras theorem. It was found that the peaks in RDF were comparable to the calculated results. On the other hand, the integral of RDF, which is accumulative, showed the coordination number of the neighbouring atoms, and this agreed with the theoretical fcc model. See figure 14 and 15. Both peak and integral of RDF gave evidence that RDF can be used to determine lattice structure. From the simulation, we know the lattice spacing, a, is 1.47361. By visualising the fcc lattice model, we know the nearest neighbour has a separation of <math>\frac{\sqrt{2}}{2}a</math>, the second nearest neighbour has a separation of <math>a</math>, and the third nearest neighbour has a separation of <math>\sqrt{\frac{3}{2}}a</math>. The results were summarised in the table below.

{| class="wikitable"
|-
! no of peak !! RDF Peak, Lattice spacing !! Calculated lattice spacing !! Coordination Number
|-
| First || a || 1.042 || 12
|-
| Second || b || 1.474 || 6
|-
| Third || c || 1. 805|| 24
|}

==Dynamical properties and the diffusion coefficient==

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

The following table shows the conditions under which the simulations were performed.
{| class="wikitable"
|-
! !! Temperature !! Density
|-
|Solid ||1.25 || 1.25
|-
| Liquid || 1.1 || 0.7
|-
|Gas || 1.1 || 0.05
|}

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

[[File:MSD of gas.png|350px|right|thumb|Fig.16 MSD of gas against timestep.]]

[[File:MSD of liquid.png|350px|right|thumb|Fig.17 MSD of liquid against timestep.]]

[[File:MSD of solid.png|350px|right|thumb|Fig.18 MSD of solid against timestep.]]

From equation <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>,
we can calculate D in reduced unit, D*, from the MSD-timestep graph, i.e. <math>D^{*} = \frac{1}{6} gradient</math>. The timestep of the simulation was 0.002 s, therefore, <math>D=\frac{D^{*}}{0.002}</math>. The results are illustrated in the table below.

{| class="wikitable"
|-
! '''D from MSD''' !! Solid !! Liquid !! Gas
|-
| Simulated Data || <math> 4.28* 10^{-6} </math>|| <math> 0.121 </math>||<math> 2.92 </math>
|-
| 1 Million Atom Data || <math> 4.39* 10^{-6} </math>|| <math> 0.0873 </math>|| <math> 3.00 </math>
|}

The D values obtained were expected. Gas had the highest D among the three, as gas molecules had the highest kinetic energy, and can move around freely with only a little restriction. Liquid had a smaller D value, while liquid molecules can move around freely, but have lower kinetic energy, pack closer together and have less space to move around. In contrast, solid had a very low D value of <math> 10^{-6} </math>. It is because solid molecules have the least kinetic energy. Also, they are ordered in a rigid shape, and can only oscillate about its equilibrium position. On the other hand, Simulated data was similar to the 1 million atom data given, which showed consistency of the the simulation.

'''<big>TASK 23</big>: In the theoretical section at the beginning, the equation for the evolution of the position of a 1D harmonic oscillator as a function of time was given. Using this, evaluate the normalised velocity autocorrelation function for a 1D harmonic oscillator (it is analytic!):

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

The position of the oscillator is defined as <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

The velocity is the first derivative of the position equation: <math> v = \frac{\mathrm{d}\mathbf{x\left(t\right)}}{\mathrm{d}t} = -A\omega\sin(\omega t+\phi) </math>

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

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

And<math> - A^{2} \omega^{2} </math> can be cancelled out.

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

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

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

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

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

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

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

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

We know that Sine is an odd function and cosine is an even function. Also, odd function multiplies even function gives in an odd function. The integral of an odd function with a symmetrical limit becomes 0.

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

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

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

Figure 19 compared the VACF of the solid and liquid simulation and the normalised harmonic oscillator. Solid simulation showed a minima at timestep = 47, while liquid simulatoin did not give a minima, its VACF flucatuated about 0 after teimstep = 109. Minima was observed because there are atom collisions in the solid and liquid simulation, and caused change in direction. Moreover, collisions in Lennard-Jones system will lead to loss of velocity and loss of energy, hence, the VACF is no longer correlated with time after a period of time and remain at 0. In contrast, harmonic oscillator is shown as a periodic function because there is no atom collision within the system. Therefore, no lost of energy and the oscillator will oscillate forever.

[[File:VACF.png|350px|right|thumb|Fig.19 VACF of solid, liquid and harmonic oscillator against timestep.]]

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

{| class="wikitable"
|-
! '''D from VACF''' !! Solid !! Liquid !! Gas !! D from MSD!! Solid !! Liquid !! Gas
|-
| Simulated Data || <math> 5.51 * 10^{-4}</math>|| <math> 0.113 </math> || <math> 0.885 </math> || Simulated Data || <math>4.28*10^{-4}</math> || <math>0.121</math>|| <math>2.92</math>
|-
| 1 Million Atom Data || <math> 2.78 * 10^{-4} </math>|| <math> 0.0825 </math> || <math> 0.978 </math> ||1 Million Atom Data || <math>4.39*10^{-4}</math> || <math>0.0873 </math>|| <math>3.00</math>
|}
[[File:Integral of VACF of gas.png|350px|right|thumb|Fig.20 Integral of VACF of gas against timestep.]]

[[File:Integral of VACF of liquid.png|350px|right|thumb|Fig.21 Integral of VACF of liquid against timestep.]]

[[File:Integral of VACF of solid.png|350px|right|thumb|Fig.22 Integral of VACF of solid against timestep.]]

[[File:Integral of VACF of 1 million gas.png|350px|right|thumb|Fig.23 Integral of VACF of 1 million gas molecules against timestep.]]

[[File:Integral of VACF of 1 million liquid.png|350px|right|thumb|Fig.24 Integral of VACF of 1 million liquid molecules against timestep.]]

[[File:Integral of VACF of 1 million solid .png|350px|right|thumb|Fig.25 Integral of VACF of 1 million solid against timestep.]]

Rep:Mod:nhl13

2016-02-25T21:48:17Z

Nhl13: /* Equilibration */

==='''Simulation of a simple liquid'''===

==Introduction to Molecular Dynamics Simulation==

'''<big>TASK 1</big>: Open the file HO.xls. In it, the velocity-Verlet algorithm is used to model the behaviour of a classical harmonic oscillator. Complete the three columns "ANALYTICAL", "ERROR", and "ENERGY": "ANALYTICAL" should contain the value of the classical solution for the position at time <math>t</math>, "ERROR" should contain the ''absolute'' difference between "ANALYTICAL" and the velocity-Verlet solution (i.e. ERROR should always be positive -- make sure you leave the half step rows blank!), and "ENERGY" should contain the total energy of the oscillator for the velocity-Verlet solution. Remember that the position of a classical harmonic oscillator is given by <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math> (the values of <math>A</math>, <math>\omega</math>, and <math>\phi</math> are worked out for you in the sheet).'''

[[File:position against time.PNG|thumbnail|thumb|350px|right|Fig.1 Analytical position vs. Time.]]

[[File:Energy against time.png|thumbnail|thumb|350px|right|Fig.2 Total Energy from velocity-Verlet solution vs. Time.]]

[[File:Error against time.PNG|thumbnail|thumb|350px|right|Fig.3 Error vs. Time.]]

Analytical position of a classical harmonic oscillator was plotted in figure 1, where the analytical position was calculated by equation:

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

where <math>A=1, \omega = 1, \phi=0</math>

The total energy of the classical harmonic oscillator for the velocity-Verlet solution was calculated by equation:

<math>Total Energy = Kinetic Energy + Potential Energy</math>,

For a classical harmonic oscillator, <math>KE=\frac{1}{2}\ mv^2</math> and <math>PE=\frac{1}{2}\ kx^2</math>,

and x(t) from velocity-Verlet solution was used instead of analytical x(t).

The total energy of the oscillator was plotted in figure 2. The energy fluctuated between 0.500 and 0.49875, which was relatively small when compared to the maximum total energy, i.e. 0.25% of the total energy.

The error in position between classical solution and velocity-Verlet solution was plotted in figure 3. It showed that the error gradually increased with time.

'''<big>TASK 2</big>: For the default timestep value, 0.1, estimate the positions of the maxima in the ERROR column as a function of time. Make a plot showing these values as a function of time, and fit an appropriate function to the data.'''

Figure 4 showed the maxima of error in position between the classical and velocity-Verlet solution respectively was plotted against time. It gave a straight line with a positive slope that indicated the error was directly proportional to time. The corresponding function was y = 0.0004x, when the line was set to intercept at origin.

[[File:Position of error maxima.PNG|thumbnail|thumb|350px|right|Fig 4.Position of error maxima vs. Time.]]

'''<big>TASK 3</big>: Experiment with different values of the timestep. What sort of a timestep do you need to use to ensure that the total energy does not change by more than 1% over the course of your "simulation"? Why do you think it is important to monitor the total energy of a physical system when modelling its behaviour numerically?'''

[[File:Different timestep.PNG|thumbnail|thumb|350px|right|Fig.5 Fluctuation in energy with different timestep.]]

The maximum total energy of the harmonic oscillator was 0.500. From figure 5, the degree of fluctuation increased when time step increased. The total energy had a fluctuation of 5x10^-3, i.e. 1% of the maximum energy, when the timestep reached 0.20 s. Hence, the timestep had to be equal or smaller than 0.20 s in order to obtain an accurate simulation.

'''<big>TASK 4.1:</big> For a single Lennard-Jones interaction, <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>, find the separation, <math>r_0</math>, at which the potential energy is zero.'''

When potential energy <math>\phi = 0</math>, <math>r=r_{0}</math>

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

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

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

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

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

'''<big>TASK 4.2</big>: What is the force at this separation?

F (r) = - dU/dr

F (r) = - d [4 ε (σ12 / r12) + 4 ε (σ6 / r06)] / dr

F (r) = d [ - 4 ε (σ12 / r12) + 4 ε (σ6 / r06)] / dr

F (r) = d [ - 4 ε (σ12 r-12) + 4 ε (σ6 r-6)] / dr

F (r) = 48 ε σ12 r-13 - 24 ε σ6 r-7

When r = r0, r0 = σ

F (r0) = 48 ε σ12 σ-13 - 24 ε σ6 σ-7

F (r0) = 48 ε σ-1 - 24 ε σ-1

F (r0) = 24 ε / σ

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

When F = 0, r = req

F (r) = 48 ε σ12 r-13 - 24 ε σ6 r-7

0 = 48 ε σ12 req-13 - 24 ε σ6 req-7

48 ε σ12 req-13 = 24 ε σ6 req-7

2 σ12 / σ6 = req13 / req7

2 σ6 = req6

req = 21/6 σ

When req = 21/6 σ,

Φ = 4 ε [(σ12 /r12 ) - (σ6 / r6)]

Φ = 4 ε {[σ12 /(21/6 σ)12] - [σ6 / (21/6 σ)6]}

Φ = 4 ε [(1/4) - (1/2)]

Φ = - ε

'''<big>TASK 4.4</big>: Evaluate the integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math> when <math>\sigma = \epsilon = 1.0</math>.

When <math>\sigma=\epsilon=1.0</math>

= ʃ 4 ε [(σ12 / r12 ) - (σ6 / r6)] dr

= [-4 / 11 r11 + 4 / 5 r5]

Therefore,

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

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

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

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

Density of water = 1 g cm-3

Mass of water in 1 ml of water = 1 g

No. of moles in 1 ml of water = 1 g ÷ 18.0 g mol-1 = 0.0556 mol

No. of water molecules in 1 ml of water = 0.0556 mol × 6.022×1023 = 3.34×1022 molecules

No. of moles of 10000 water molecules = 10000 ÷ 6.022×1022 = 1.66×10-20 mol

Mass of 10000 water molecules = 1.66×10-20 mol × 18.0 g mol-1 = 2.99×10-19 g

volume of 10000 water molecules = 2.99×10-19 g ÷ 1 g cm-1 = 2.99×10-19 ml

'''<big>TASK 6</big>: Consider an atom at position <math>\left(0.5, 0.5, 0.5\right)</math> in a cubic simulation box which runs from <math>\left(0, 0, 0\right)</math> to <math>\left(1, 1, 1\right)</math>. In a single timestep, it moves along the vector <math>\left(0.7, 0.6, 0.2\right)</math>. At what point does it end up, ''after the periodic boundary conditions have been applied''?'''.

Periodic boundary condition states that boxes are repeated infinitely in all direction, hence, if one molecule exits the boundary of the original box, one molecule from the replicated box enters the original box through the opposite face to retain the number of molecules. As a result, we could work out the final position of the atom simply by adding the movement vector to its initial position vector, i.e. (0.5, 0.5, 0.5) + (0.7, 0.6, 0.2) = (1.2, 1.1, 0.7). As the box had a dimension of (1.0, 1.0, 1.0), so the resulting position vector was (0.2, 0.1, 0.7).

'''<big>TASK 7.1</big>: The Lennard-Jones parameters for argon are <math>\sigma = 0.34\mathrm{nm}, \epsilon\ /\ k_B= 120 \mathrm{K}</math>. If the LJ cutoff is <math>r^* = 3.2</math>, what is it in real units? '''

r* = r / σ

r = r* σ

r = (3.2)(0.34×10-9)

r = 1.088 nm

'''<big>TASK 7.2</big>: What is the well depth in <math>\mathrm{kJ\ mol}^{-1}</math>? '''

Φ = - ε, as calculated above, and ε / kB = 120K

Φ = - 120 × kB

Φ = - 1.66×10-21 J ÷ 1000 × NA

Φ = 0.998 kJmol-1

'''<big>TASK 7.3</big>: What is the reduced temperature <math>T^* = 1.5</math> in real units?'''

T* = kBT ÷ ε

T = (T*)(ε) ÷ kB

T = (1.5)(120)(kB) ÷ kB

T = (1.5)(120)

T = 180 K

==Equilibration==

'''<big>TASK 8</big>: Why do you think giving atoms random starting coordinates causes problems in simulations? Hint: what happens if two atoms happen to be generated close together?'''

In liquid, molecules are held together by intermolecular bonds. They are not fixed and can flow around and have a certain average intermolecular distance. If the molecules are too far apart, the attractive interaction between them is not as strong as it used to and behaves as gas instead of liquid. In contrast, if molecules are too close together, they experience a repulsion force between them. The liquid simulation would be affected in both extremes.

'''<big>TASK 9</big>: Satisfy yourself that this lattice spacing corresponds to a number density of lattice points of <math>0.8</math>. Consider instead a face-centred cubic lattice with a lattice point number density of 1.2. What is the side length of the cubic unit cell?'''

Equation of number density can be used to rationalise the relationship between the number density, and the lattice spacing, a.

Case 1: In a simple cubic lattice, there is only 1 atom inside the lattice.

number density = number of atom / volume

0.8 = 1 / a3

a = 1.0772 m

This result agreed with the value obtained from LAMMPS.

Case 2: In a face-centred cubic lattice, there are 1 atom at each corner and 1 atom on each face. Atom at the corner and on the face contribute to 1/4 and 1/2 atom to the lattice. Hence, number of atom in a fcc lattice = 1/8 × 8 + 1/2 × 6 = 4

number density = number of atom / volume

1.2 = 4 / a3

a = 1.494 m

Therefore, a fcc lattice with a number density of 1.2 will have a side length of 1.494 m.

'''<big>TASK 10</big>: Consider again the face-centred cubic lattice from the previous task. How many atoms would be created by the create_atoms command if you had defined that lattice instead?'''

Each fcc unit cell consists of 4 atoms. With the same 10×10×10 box, we are having 1000 unit cells, hence, fcc lattice will consist of 4×1000 = 4000 atoms.

'''<big>TASK 11</big>: Using the [http://lammps.sandia.gov/doc/Section_commands.html#cmd_5 LAMMPS manual], find the purpose of the following commands in the input script:'''

<pre>mass 1 1.0</pre>Syntax of the command<pre>mass I value</pre>Mass command defines the mass of atoms, where I is the number of atom type and value is the mass of atom.

In our case, there is only 1 type of atom and its mass is 1.0 in reduced unit.<pre>pair_style lj/cut 3.0</pre>Syntax of the command<pre>pair_style style args</pre>pair_style command compute different styles of pairwise interaction. The second <nowiki>''style''</nowiki> defines a particular type of pairwise interaction, and the <nowiki>''args''</nowiki> is the arguments used by that particular style. In our case, lj/cut is used to compute the cutoff Lennard-Jones potential with no Coulomb, and the cutoff potential is set at 3.0 reduced unit.<pre>pair_coeff * * 1.0 1.0</pre>Syntax of the command

pair_coeff command specifies the pairwise force field coefficients for one or more pairs of atom types. <nowiki>''I,J'' is the atom types and ''args''</nowiki> is the coefficients for one or more pairs of atom types. In out case, I and J equal * which means there is N types of atom and means all types from 1 to N, and both coefficients are set at 1.0.

'''<big>TASK 12</big>: Given that we are specifying <math>\mathbf{x}_i\left(0\right)</math> and <math>\mathbf{v}_i\left(0\right)</math>, which integration algorithm are we going to use?'''

We are using Velocity Verlet Algorithm, see below.

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

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

### RUN SIMULATION ###
run ${n_steps}
run 100000
</pre>
'''The second line (starting "variable timestep...") tells LAMMPS that if it encounters the text ${timestep} on a subsequent line, it should replace it by the value given. In this case, the value ${timestep} is always replaced by 0.001. In light of this, what do you think the purpose of these lines is? Why not just write:'''
<pre>
timestep 0.001
run 100000</pre>
Those extra lines would allow the variable to change automatically when timestep in the second line is changed.

'''<big>TASK 14</big>: make plots of the energy, temperature, and pressure, against time for the 0.001 timestep experiment (attach a picture to your report). Does the simulation reach equilibrium? How long does this take? When you have done this, make a single plot which shows the energy versus time for all of the timesteps (again, attach a picture to your report). Choosing a timestep is a balancing act: the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality; short timesteps, however, mean that the same number of simulation steps cover a shorter amount of actual time, and this is very unhelpful if the process you want to study requires observation over a long time. Of the five timesteps that you used, which is the largest to give acceptable results? Which one of the five is a ''particularly'' bad choice? Why?'''[[File:total energy equilibration.PNG|thumbnail|thumb|350px|right|Fig6. Equilibration of total energy.]][[File:temperature equilibration.PNG|thumbnail|thumb|350px|right|Fig.7 Equilibration of Temperature.]][[File:pressure equilibration.PNG|thumbnail|thumb|350px|right|Fig.8 Equilibration of Pressure.]][[File:total energy at different timestep.PNG|thumbnail|thumb|350px|right|Fig.9 Total energy at different timestep.]]

All three properties reached equilibrium after short period of time when timestep was 0.001 s. The total energy, temperature and pressure reached the equilibrium after approximately 0.42, 0.32 and 0.25 s respectively. When timestep was as large as 0.015 s, the total energy did not reach equilibrium. This was particularly bad for simulation. When the timestep was equal to or lower than 0.01 s, the energy did reach equilibrium. However, the fluctuation of the energy level was quite large when timestep equal 0.01 and 0.075 s. Also, the timestep cannot be too short, as it will lower the efficiency of the simulation. As a result, the most acceptable timestep for simulation to compensate both disadvantages was 0.0025 s.

==Running simulations under specific conditions==

'''<big>TASK 15</big>: Choose 5 temperatures (above the critical temperature <math>T^* = 1.5</math>), and two pressures (you can get a good idea of what a reasonable pressure is in Lennard-Jones units by looking at the average pressure of your simulations from the last section). This gives ten phase points — five temperatures at each pressure. Create 10 copies of npt.in, and modify each to run a simulation at one of your chosen <math>\left(p, T\right)</math> points. You should be able to use the results of the previous section to choose a timestep. Submit these ten jobs to the HPC portal. While you wait for them to finish, you should read the next section.'''

In our simulation, timestep of 0.0025 s was used, as it allowed the oscillator to reach equilibrium and had a small fluctuation in energy level. The two pressures used were 2.55 and 2.65 atm, as there were inside the range of the pressure equilibrium from the last section. Finally, the temperature used were 1.7, 1.9, 2.1, 2.3 and 2.5, where all of them were above the critical temperature.

'''<big>TASK 16</big>: We need to choose <math>\gamma</math> so that the temperature is correct <math>T = \mathfrak{T}</math> if we multiply every velocity <math>\gamma</math>. We can write two equations:'''

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

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

'''Solve these to determine <math>\gamma</math>.'''

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

<math>\frac{\sum_i m_i v_i^2}{3T} = N k_B</math>

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

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

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

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

'''<big>TASK 17</big>: Use the [http://lammps.sandia.gov/doc/fix_ave_time.html manual page] to find out the importance of the three numbers ''100 1000 100000''. How often will values of the temperature, etc., be sampled for the average? How many measurements contribute to the average? Looking to the following line, how much time will you simulate?'''

100 is the Nevery, 1000 is the Nrepeat and 100000 is the Nfreq. Nevery defines how often are the input values is used for calculation, Nrepeat defines the number of times that the input values are used to calculate the average, and Nfreq defines the number of timestep for calculating an average. In our case, the input values are used every 100 timesteps, average calculated from input values are repeated for 1000 times, and the averages are calculated every 100000 timesteps. We used timestep of 0.0025 s, so our simulation time is 100000 x 0.0025 = 250 (reduced unit)

'''<big>TASK 18</big>: When your simulations have finished, download the log files as before. At the end of the log file, LAMMPS will output the values and errors for the pressure, temperature, and density <math>\left(\frac{N}{V}\right)</math>. Use software of your choice to plot the density as a function of temperature for both of the pressures that you simulated. Your graph(s) should include error bars in both the x and y directions. You should also include a line corresponding to the density predicted by the ideal gas law at that pressure. Is your simulated density lower or higher? Justify this. Does the discrepancy increase or decrease with pressure?'''

[[File:density vs temp.PNG|thumbnail|thumb|350px|right|Fig.10 Density vs. Temp at different pressure.]]

Figure 10 showed that the simulated density decreases with temperature and increases with pressure respectively. This could be rationalised easily. As temperature increases, the molecules have higher kinetic energy and overcome the attractive interaction between molecules, thus a lower density. On the other hand, as pressure increases, the molecules are experiencing a greater force and are compressed together, thus a higher density. The error in density was much smaller than the error in temperature, shown in figure 10. The simulated density was much lower than the density calculated from the ideal gas law. The ideal gas law does not take into account of the repulsive force between molecules, while the simulation was able to mimic the repulsive force. Therefore, the simulated density was lower than the calculated density as repulsive force between molecules pushed each other away. The discrepancy between simulated and calculated density increases with pressure. This could be rationalised by knowing that the gas behaves more as a real gas instead of ideal gas when pressure increase, hence the discrepancy increases.

==Calculating heat capacities using statistical physics==

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

[[File:heat capacity against temperature.png|350px|right|thumb|fig.11 heat capacity against temperature at different density]]

Cv/V was plotted against temperature in figure 11. For both density, the simulation indicated the volume of the liquid remain unchanged, i.e. 16875 and 4218.75. The graph showed that the Cv/V decrease with temperature. The extent of decrease was greater when temperature was small and started to decrease slower as temperature increase gradually. However, there was an anomalous when T=2.4, D=0.8, which possibly arose due Brownian motion of the liquid.

Example of the script is shown below, where density is 0.2 and temperature is 2.0.

<pre>

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

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

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

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

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

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

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

### SPECIFY TIMESTEP ###

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

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

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp press density atoms vol
variable dens equal density
variable dens2 equal density*density
variable temp equal temp
variable temp2 equal temp*temp
variable press equal press
variable press2 equal press*press
variable N equal atoms
variable N2 equal atoms*atoms
variable etotal equal etotal
variable etotal2 equal etotal*etotal
variable volume equal vol

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

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

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

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

</pre>

==Structural properties and the radial distribution function==

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

The following table shows the conditions under which the simulations of the Lennard-Jones system were performed.
{| class="wikitable"
|-
! !! Temperature !! Density
|-
|Solid ||1.25 || 1.25
|-
| Liquid || 1.1 || 0.7
|-
|Gas || 1.1 || 0.05
|}

[[File:RDF for solid, liquid and gas.PNG|thumbnail|thumb|350px|right|Fig.12 RDF for solid, liquid and gas.]]

[[File:Integral of RDF for solid, liquid and gas.PNG|thumbnail|thumb|350px|right|Fig.13 Integral of RDF for solid, liquid and gas.]]

RDF provides the density of atoms at particular radius to a reference atom. Each peak in RDF indicates a particularly favoured separation distance for the neighbouring atoms to a reference atom. Thus, RDF reveals the atomic structure of the system being simulated. Figure 12 showed that RDF of solid had lots of sharp peaks, while there were only a few broad peaks in liquid, and just one broad peak for gas. It reveals atoms in solid are highly ordered and stays at particular distance to the reference atom. Sharp peak also indicates the atoms are not spread out randomly but are concentrated at favoured particular distances. The coordination number at favoured particular distance can be determined from the integral of RDF. In contrast, a few peaks in liquid RDF shows its system is less ordered and atoms are able to move around freely and rarely stay at a particular distance to the reference atom. Peaks still shows up illustrates that the atoms are still held together by some sort of attraction forces. On the other hand, broad peaks are observed in RDF of liquid and gas, where peak for gas is broader than those for liquid. This can be rationalised by the Brownian motion of the atoms in the system. The atoms in gas have higher energy, therefore, can move more randomly and resulted in a broad distribution at particular distance, where the only peak observed indicates the gas atoms are still interacting with each other.

[[File:1st and 2nd coordinating neighbours of fcc.PNG|thumbnail|thumb|350px|right|Fig.14 1st and 2nd coordinating neighbours of fcc.]]

[[File:3rd coordinating neighbours of fcc.PNG|thumbnail|thumb|350px|right|Fig.15 3rd coordinating neighbours of fcc.]]

The first three peaks at the RDF of soild were at a, b and c. These corresponded to the first three nearest coordinating neighbours to a reference atom. The lattice spacing between the reference atom and the nearest neighbours can be calculated by Pythagoras theorem. It was found that the peaks in RDF were comparable to the calculated results. On the other hand, the integral of RDF, which is accumulative, showed the coordination number of the neighbouring atoms, and this agreed with the theoretical fcc model. See figure 14 and 15. Both peak and integral of RDF gave evidence that RDF can be used to determine lattice structure. From the simulation, we know the lattice spacing, a, is 1.47361. By visualising the fcc lattice model, we know the nearest neighbour has a separation of <math>\frac{\sqrt{2}}{2}a</math>, the second nearest neighbour has a separation of <math>a</math>, and the third nearest neighbour has a separation of <math>\sqrt{\frac{3}{2}}a</math>. The results were summarised in the table below.

{| class="wikitable"
|-
! no of peak !! RDF Peak, Lattice spacing !! Calculated lattice spacing !! Coordination Number
|-
| First || a || 1.042 || 12
|-
| Second || b || 1.474 || 6
|-
| Third || c || 1. 805|| 24
|}

==Dynamical properties and the diffusion coefficient==

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

The following table shows the conditions under which the simulations were performed.
{| class="wikitable"
|-
! !! Temperature !! Density
|-
|Solid ||1.25 || 1.25
|-
| Liquid || 1.1 || 0.7
|-
|Gas || 1.1 || 0.05
|}

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

[[File:MSD of gas.png|350px|right|thumb|Fig.16 MSD of gas against timestep.]]

[[File:MSD of liquid.png|350px|right|thumb|Fig.17 MSD of liquid against timestep.]]

[[File:MSD of solid.png|350px|right|thumb|Fig.18 MSD of solid against timestep.]]

From equation <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>,
we can calculate D in reduced unit, D*, from the MSD-timestep graph, i.e. <math>D^{*} = \frac{1}{6} gradient</math>. The timestep of the simulation was 0.002 s, therefore, <math>D=\frac{D^{*}}{0.002}</math>. The results are illustrated in the table below.

{| class="wikitable"
|-
! '''D from MSD''' !! Solid !! Liquid !! Gas
|-
| Simulated Data || <math> 4.28* 10^{-6} </math>|| <math> 0.121 </math>||<math> 2.92 </math>
|-
| 1 Million Atom Data || <math> 4.39* 10^{-6} </math>|| <math> 0.0873 </math>|| <math> 3.00 </math>
|}

The D values obtained were expected. Gas had the highest D among the three, as gas molecules had the highest kinetic energy, and can move around freely with only a little restriction. Liquid had a smaller D value, while liquid molecules can move around freely, but have lower kinetic energy, pack closer together and have less space to move around. In contrast, solid had a very low D value of <math> 10^{-6} </math>. It is because solid molecules have the least kinetic energy. Also, they are ordered in a rigid shape, and can only oscillate about its equilibrium position. On the other hand, Simulated data was similar to the 1 million atom data given, which showed consistency of the the simulation.

'''<big>TASK 23</big>: In the theoretical section at the beginning, the equation for the evolution of the position of a 1D harmonic oscillator as a function of time was given. Using this, evaluate the normalised velocity autocorrelation function for a 1D harmonic oscillator (it is analytic!):

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

The position of the oscillator is defined as <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

The velocity is the first derivative of the position equation: <math> v = \frac{\mathrm{d}\mathbf{x\left(t\right)}}{\mathrm{d}t} = -A\omega\sin(\omega t+\phi) </math>

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

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

And<math> - A^{2} \omega^{2} </math> can be cancelled out.

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

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

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

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

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

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

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

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

We know that Sine is an odd function and cosine is an even function. Also, odd function multiplies even function gives in an odd function. The integral of an odd function with a symmetrical limit becomes 0.

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

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

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

Figure 19 compared the VACF of the solid and liquid simulation and the normalised harmonic oscillator. Solid simulation showed a minima at timestep = 47, while liquid simulatoin did not give a minima, its VACF flucatuated about 0 after teimstep = 109. Minima was observed because there are atom collisions in the solid and liquid simulation, and caused change in direction. Moreover, collisions in Lennard-Jones system will lead to loss of velocity and loss of energy, hence, the VACF is no longer correlated with time after a period of time and remain at 0. In contrast, harmonic oscillator is shown as a periodic function because there is no atom collision within the system. Therefore, no lost of energy and the oscillator will oscillate forever.

[[File:VACF.png|350px|right|thumb|Fig.19 VACF of solid, liquid and harmonic oscillator against timestep.]]

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

{| class="wikitable"
|-
! '''D from VACF''' !! Solid !! Liquid !! Gas !! D from MSD!! Solid !! Liquid !! Gas
|-
| Simulated Data || <math> 5.51 * 10^{-4}</math>|| <math> 0.113 </math> || <math> 0.885 </math> || Simulated Data || <math>4.28*10^{-4}</math> || <math>0.121</math>|| <math>2.92</math>
|-
| 1 Million Atom Data || <math> 2.78 * 10^{-4} </math>|| <math> 0.0825 </math> || <math> 0.978 </math> ||1 Million Atom Data || <math>4.39*10^{-4}</math> || <math>0.0873 </math>|| <math>3.00</math>
|}
[[File:Integral of VACF of gas.png|350px|right|thumb|Fig.20 Integral of VACF of gas against timestep.]]

[[File:Integral of VACF of liquid.png|350px|right|thumb|Fig.21 Integral of VACF of liquid against timestep.]]

[[File:Integral of VACF of solid.png|350px|right|thumb|Fig.22 Integral of VACF of solid against timestep.]]

[[File:Integral of VACF of 1 million gas.png|350px|right|thumb|Fig.23 Integral of VACF of 1 million gas molecules against timestep.]]

[[File:Integral of VACF of 1 million liquid.png|350px|right|thumb|Fig.24 Integral of VACF of 1 million liquid molecules against timestep.]]

[[File:Integral of VACF of 1 million solid .png|350px|right|thumb|Fig.25 Integral of VACF of 1 million solid against timestep.]]

File:Position of error maxima.PNG

2016-02-25T21:43:06Z

Nhl13: Nhl13 uploaded a new version of File:Position of error maxima.PNG

2016-02-25T20:06:04Z

Nhl13:

File:Integral of VACF of 1 million liquid.png

2016-02-25T19:57:31Z

Nhl13:

File:Integral of VACF of 1 million gas.png

2016-02-25T19:57:07Z

Nhl13:

File:Integral of VACF of solid.png

2016-02-25T19:56:44Z

Nhl13:

File:Integral of VACF of liquid.png

2016-02-25T19:56:23Z

Nhl13: