File:MSD Solid RS7913.png

2018-03-27T21:12:49Z

Rs7913: Rs7913 uploaded a new version of File:MSD Solid RS7913.png

File:MSD Solid RS7913.png

2018-03-27T21:12:13Z

Rs7913: Rs7913 uploaded a new version of File:MSD Solid RS7913.png

Rep:Molecular Dynamics RS7913

2018-03-27T21:04:10Z

Rs7913:

== Tasks ==

'''<big>TASK</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).'''

See "Analytical" tab in [https://wiki.ch.ic.ac.uk/wiki/images/8/80/HO_RS7913.xls HO_RS7913.xls]

'''<big>TASK</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.'''

See "Analytical" tab in [https://wiki.ch.ic.ac.uk/wiki/images/8/80/HO_RS7913.xls HO_RS7913.xls]

'''<big>TASK</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?'''

A timestep of 0.2 is just about short enough to ensure that the total energy doesn't change by more than 1% over the course of the simulation, since it results in the total energy oscillating between 0.5 and 0.495. The total energy of this simple harmonic oscillator was calculated by summing the kinetic energy (calculated using the velocity from the velocity-Verlet solution) and the potential energy (calculated using the classical simple harmonic motion equation to calculate the force, F = -kx). In a simple harmonic oscillator, the total energy of the system remains constant as kinetic and potential energy are interchanged. It is important to monitor the total energy of a physical system such as this because the amount by which the total energy varies corresponds to the accuracy of the simulation. Of course, the more accurate the approximation, the smaller the amount by which the energy will vary, and a smaller timestep will result in a more accurate approximation.

See "Timestep" tab in [https://wiki.ch.ic.ac.uk/wiki/images/8/80/HO_RS7913.xls HO_RS7913.xls]

'''<big>TASK:</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. What is the force at this separation? Find the equilibrium separation, <math>r_{eq}</math>, and work out the well depth (<math>\phi\left(r_{eq}\right)</math>). 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>.

[[File:TASK2 RS7913.png|frameless]]

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

[[File:TASK3.png RS7913|frameless|640x640px]]

'''<big>TASK</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''?'''

(0.5, 0.5, 0.5) + (0.7, 0.6, 0.2) = (0.2, 0.1, 0.7)

'''<big>TASK</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? 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?

[[File:TASK4 RS7913.png|frameless|600x600px]]

'''<big>TASK</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?'''

If random starting coordinates were given to the atoms, two or more atoms may be generated very close to one another. This would mean that these atoms would experience a very strong repulsive force between one another, as given by the Lennard-Jones potential, since the '1/r12' term dominates at shorter distance. Hence, the system may not reach equilibrium by the time the simulation has run - although increasing the length of time a simulation is run would increase the chance of the system reaching equilibrium, this would be of computational cost.

'''<big>TASK</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?'''

p = m/v

Volume of a lattice point = (1.07722)3

Therefore, number density of lattice points = 1/(1.07722)3 = 0.8

Face-centered cubic (FCC) lattice contains 4 lattice points within the unit cell

Therefore, the volume of one lattice point = 4/1.2

Side length of FCC unit cell = (4/1.2)1/3 = 1.49380

'''<big>TASK</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?'''

The create_atoms command created 1000 atoms for the simple cubic lattice, for which each unit cell contains 1 atom. Since the FCC unit cell contains 4 atoms, the create_atoms command would have created 1000 x 4 = 4000 atoms if this lattice had been defined instead.

'''<big>TASK</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
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

"mass 1 1.0" sets the mass of atom type 1 to 1.0

"pair_style lj/cut 3.0" sets the formula used to calculate pairwise interactions between atoms as that defined by the Lennard-Jones potential, with a cutoff of r* = 3.0

"pair_coeff * * 1.0 1.0" sets the forcefield coefficients - in this case the Lennard-Jones parameters σ and ε - to 1.0

'''<big>TASK</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 the velocity-Verlet algorithm, since the starting positions of the atoms''' '''('''<math>\mathbf{x}_i\left(0\right)</math>''') and their velocities at the same time''' '''('''<math>\mathbf{v}_i\left(0\right)</math>''') must be defined when using the velocity-Varlet. When using the Verlet algorithm, the starting positions of the atoms ('''<math>\mathbf{x}_i\left(0\right)</math>''') and their positions one timestep in the past (<math>\mathbf{x}_i\left(-\delta t\right)</math>) must be defined.

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

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

This allows for the timestep to be easily modified by modifying the value "${timestep}", and other sections of the script can also be linked to "${timestep}" to ensure that the timestep for all these sections of the script is set to the desired "${timestep}" value.

'''<big>TASK</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?'''

Plots of the total energy, temperature, and pressure, against time for the 0.001 timestep experiment are shown below. It can be seen that the simulation reaches equilibrium at around 0.5 s, after which the total energy oscillates about the average equilibrium value.
[[File:EnergyVsTime RS7913.png|left|frameless|538x538px]]
[[File:TempVsTime RS7913.png|frameless|538x538px]]
[[File:PressureVsTime_RS7913.png|left|frameless|538x538px]]
[[File:TotEnergyVsTAllTimesteps RS7913.png|frameless|538x538px]]

The total energy against time for all timesteps has also been plotted above (bottom right graph). It is clear that the average total energy for the experiments carried out using time steps of 0.001 s and 0.0025 s were very similar; the total energy for both of these experiments oscillated about virtually the same average equilibrium value. The 0.0075 s and 0.01 s timestep experiments both reached equilibrium, however the average energy at equilibrium was considerably greater for the 0.0075 s experiment and even greater for the 0.01 s experiment, and this is indicative that using these experiments are not as accurate. In addition, the amount by which the total energy oscillated once equilibrium was reached was greater for the 0.0075 s experiment and even greater for the 0.01 s experiment, further testifying that these simulations were not as accurate. In the case of the 0.015 s timestep, the simulation did not even reach equilibrium over the duration of the simulation; the energy just continued increasing further and further. This 0.015 s timestep simulation is clearly not accurate, and using such a large timestep is a particularly bad choice. Whilst using the 0.001 s timestep would give the most accurate results, the 0.0025 s timestep is sufficiently short to give acceptable results since the average total energies from both of these calculations are virtually the same. Using a timestep greater than 0.0025 s would likely give inaccurate results.

'''<big>TASK</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.'''

Timestep of 0.0025 chosen

T = 2, 3.5, 5, 6.5, 8

p = 2.5, 3

'''<big>TASK</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>.'''

[[File:TASK 5 RS7913.png|frameless|263x263px]]

'''<big>TASK</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?'''

The number 100 sets the Nevery variable, which is the number of timesteps between which each sample of the (thermodynamic) data is taken.

The number 1000 sets the Nrepeat variable, which is the number of timesteps between which averages of this sampled data (set by Nevery) are taken.

The number 100000 sets the Nfreq variable, which is the number of timesteps over which an average of the previous averages (determined by Nrepeat) are taken.

The following line "run 100000" instructs LAMMPS to run the simulation over 100000 timesteps. Hence, setting the Nfreq variable to 100000 as done for this simulation essentially instructs LAMMPS to take an average over the whole course of the simulation.

'''<big>TASK</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:TASK 6 RS7913.png|frameless|800x800px]]

The simulated density, at both reduced pressures of 2.5 and 3.0, was lower than that predicted by the ideal gas law. This is because the ideal gas law models the atoms as simply hard spheres, meaning the atoms cannot get as close to one another as in the molecular dynamics simulations carried out, in which the Lennard-Jones potential is used to model intermolecular interactions. This function is a more accurate representation of the true intermolecular interactions, and accounts for the more gradual onset of repulsion between two atoms as they are brought closer together, caused by electrostatic repulsion between the outer electrons - the deformability of atoms is not accounted for by the ideal gas law. This discrepancy increases with increasing pressure since, at higher pressure, the particles spend more time closer to one another, and bump into one another more often, therefore the deformability of atoms accounted for in the molecular dynamics simulation becomes more significant relative to the ideal gas law hard sphere model.

'''<big>TASK</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:TASK 7 RS7913.png|frameless|740x740px]]

For both densities of 0.2 and 0.8, the heat capacity at constant volume, [[C]]v, decreases with increasing temperature. This trend would not be expected; the heat capacity per unit volume should remain constant with respect to temperature because the inherent chemical properties of the system have not changed. Since both the density and volume - because we are considering the heat capacity of some unit volume V (Cv/V) - are fixed, the number of particles being considered does not change as the temperature is varied, hence the internal energy remains constant. Since the heat capacity at constant volume can be defined as Cv = dU/dT, the fact that the internal energy remains constant as the temperature is varied inherently means that the heat capacity also remains constant as the temperature is varied.

The heat capacity, [[C]]v, also increases with increasing density, which would be the expected trend. Since we are considering the heat capacity of some unit volume (V) of the system (Cv/V), a higher density of atoms means a greater number of atoms in that unit volume (V) being considered. It makes sense that a greater number of atoms would require more heat energy to raise the average kinetic energy of the particles by some particular amount than that required by that of fewer atoms in that same volume. Since temperature corresponds to the average kinetic energy of particles in a system, it would require more heat energy to raise the temperature of a given volume of some material of higher density by a given amount compared to that required by the same system at lower density - i.e. the more dense system has a higher heat capacity per unit volume.

The input script for the calculation carried out at a density of 0.2 and temperature of 2.0 can be found [https://wiki.ch.ic.ac.uk/wiki/images/7/7b/NvT_D_0.2_T_2.0_RS7913.in here].

'''<big>TASK</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 RDFs for the solid, liquid and gaseous systems have been plotted below:

[[File:RDF RS7913.png|frameless|740x740px]]

The RDF of the system in the gas phase has one broad peak around r = 1 - the average density value - which then tapers off to a constant value of g(r) = 1. This means that these molecules in the gas phase essentially have only one coordination sphere, and no further order. This is because the density of the gas is so much lower that longer range interactions (r > 2) are essentially non-existent, and the distribution of gas molecules at a distance r > 2.5 is essentially random.

In the liquid phase, there are 3 discernible peaks in the RDF, corresponding to the first, second and third solvation spheres. Since the density of the liquid is considerably higher than that of the gas - 10 times higher than the density of the gas - longer-range interactions become more important, and this gives rise to the short-range order that can be seen in the RDF.

In the solid phase, the atoms (have been commanded, by the software, to) occupy an ordered FCC lattice structure, in which the atoms can only move by very small amounts as the structure vibrates – the atoms are essentially held still. Hence, there are far more peaks in the RDF corresponding to the long-range order in the crystalline structure. The first, second and third peaks in the RDF, at around r = 1.03, 1.49, and 1.83 correspond to the first, second and third nearest neighbours in the FCC lattice. The second nearest neighbour is the length of one side of the cubic unit cell, which was earlier calculated to be equal to 1.49380 (for a FCC lattice with a point number density of 1.2). The value of 1.49 from the RDF agrees well with the calculated value of 1.49380. In the FCC lattice, there are 12 nearest neighbours, 6 second nearest neighbours, and 24 third nearest neighbours, and the fact that the third peak in the RDF is of considerably greater intensity (g(r)) than the second reflects the fact that there are many more (4 times as many) second nearest neighbours in comparison to first nearest neighbours. The first, second and third nearest neighbours in the FCC lattice are illustrated in the diagram below:

[[File:FCC Nearest Neighbours RS7913.png|frameless|500px]]

It is worth noting that using the LJ cutoff of r = 3.2, neglecting interactions where r > 3.2, may have given rise to some error in the case of the liquid and gas simulations. It can be seen that the RDF of both the liquid and gas simulations have decayed to a constant value of g(r) = 1 when r > 3.2. In the case of the liquid simulation, it would not be surprising if a more accurate simulation with a higher LJ cutoff gave rise to further peak(s) in the RDF when r > 3.2, as well as better defined peaks that are closer to one another for r < 3.2. In the case of the gas simulation, using a higher LJ cutoff would likely result in a broader peak than that seen in the RDF above, that decays to the constant value of g(r) = 1 at a radius greater than 3.2.

'''<big>TASK</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.'''

Simulations were run for a vapour (density = 0.05)<ref name=":0">J. P. Hansen, L. Verlet, ''Phys. Rev.'', 1969, '''184''', 151.</ref>, liquid (density = 0.8)<ref name=":0" /> and solid (density = 1.2<ref name=":0" />, with lattice type set to fcc), at a reduced temperature of 1.0. The system contains 8000 atoms (20x20x20 grid).

'''<big>TASK</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.'''

The total MSD as a function of timestep for the solid, liquid and gas simulations carried out with N = 8000 and N = 1000000 have been plotted below, and a linear fit has been applied to each of these functions:
[[File:MSD Gas RS7913.png|centre|frameless|700x700px]]
[[File:MSD Liquid RS7913.png|centre|frameless|700x700px]]
[[File:MSD Solid RS7913.png|centre|frameless|700x700px]]

Each plot of the MSD vs time (timestep) is essentially a measure of how far particles in the system have moved with respect to time. In all cases we would expect the function to be linear, since any motion of particles corresponds to random Brownian motion, or random vibrations in the case of the solid system, rather than some external force acting on the particles. An upwards curve would indicate that there is some sort of drift of the molecules, due to some net force acting on the particles, and a plateau in the function would indicate that the particle is confined, or constrained, in some sort of way. In the case of the liquid, excellent linear behaviour is maintained throughout the simulation, however, in the solid and gas phases, it does require some time for linear behaviour to be established, and this could be attributed to the time required for the simulation to equilibrate following the initial random assignment of velocities to all of the particles. In the solid phase, oscillations can be seen that fade away, which correspond to anisotropic vibrations in the system spreading over the whole solid as equilibrium is reached, and the system does equilibrate relatively quickly. The MSD function for the gas phase, however, seems to be an upwards curve, indicating some sort of external force acting on the particles. This could be attributed to the system not having reached equilibrium within the timeframe of the simulation - it may require more time to reach equilibrium in the gas simulation due to the decreased density of the system meaning fewer collisions, and hence a longer time required to distribute all the kinetic energy (roughly) evenly between all particles in the system. This may be resolved by using a smaller timestep, or by running the simulation for a longer period of time. I would not expect that increasing the number of particles in the system would result in equilibrium being reached, since both the 8000-atom and 1000000-atom simulations result in similar curves with a positive second derivative (curvature) throughout the simulation.

As expected, the gradient was greatest for the gas, followed by the liquid, whereas the gradient for the solid phase was essentially zero. This is because particles in the gas phase move the most (due to the lower density meaning movement of particles is not constricted by collisions), followed by those in the liquid. The fact that the solid function is essentially flat reflects the fact that the particles do not have net movement in the solid, but only some vibrations of course.

In the case of the liquid and gas simulations, increasing the size of the system from 8000 to 1000000 atoms doesn't seem to result in any great qualitative difference in terms of the MSD. However, in the case of the solid simulation it can clearly be seen that the MSD flattens out to a virtually straight line in the 1000000-atom simulation, whereas for the 8000-atom simulation, the function does not really completely flatten out over the course of the simulation, indicating that this smaller system size may not be quite large enough to gain accurate information in terms of the MSD.

The diffusion coefficient, D, was estimated by dividing the gradient of the linear fit by 6:

{| class="wikitable"
!
! colspan="3" |Diffusion Coefficient
|-
!Phase
!N = 8000
!N = 1000000
!Literature Value
|-
|Solid
|5.5533 x 10-3
|5.0833 x 10-3
|
|-
|Liquid
|1.6667 x 10-4
|1.6667 x 10-4
|
|-
|Gas
|8.3333 x 10-9
|8.3333 x 10-9
|
|}

As expected....

'''<big>TASK</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>

'''Be sure to show your working in your writeup. 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.'''

The normalised velocity autocorrelation function, C(τ), for a 1D harmonic oscillator was evaluated as follows:

[[File:VACF SHM Working RS7913.png|frameless|815x815px]]

C(τ) was plotted as a function of time (timestep (τ), between 0 and 500) for the 1D harmonic oscillator, as well as for the Lennard-Jones solid and liquid simulations carried out with N = 8000:

[[File:VACF Vs Timestep Harmonic RS7913.png|frameless|740x740px]]

'''<big>TASK</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?'''

The trapezium rule was used to approximate the integral under the VACF, using a function that created trapezia of 1 timestep 'width'.
[[File:VACF Running Integral N=8000 RS7913.png|left|frameless|538x538px]]
[[File:VACF Running Integral N=1000000 RS7913.png|frameless|538x538px]]

D was calculated using the by dividing the integral of C(τ) over all τ (0 to 5000) by 3.
{| class="wikitable"
!
! colspan="3" |Diffusion Coefficient
|-
!Phase
!N = 8000
!N = 1000000
!Literature Value
|-
|Solid
|0.060364158
|0.022764742
|
|-
|Liquid
|34.47327454
|45.04573814
|
|-
|Gas
|1089.389257
|1634.232899
|
|}

2018-03-27T01:45:58Z

Rs7913: Rs7913 uploaded a new version of File:FCC Nearest Neighbours RS7913.png

File:FCC Nearest Neighbours RS7913.png

2018-03-27T01:44:11Z

Rs7913: Rs7913 uploaded a new version of File:FCC Nearest Neighbours RS7913.png

2018-03-27T00:39:07Z

Rs7913: