View source for Talk:Mod:P0TAT0

'''<span style="color:red">JC: General comments: Most tasks answered, but a few parts missing from some questions.  Make sure that you fully understand the theory behind this experiment and what you are being asked to calculate and plot in each task.</span>'''

= Molecular Dynamics Simulations of Simple Liquids =
[[User:Ab9314|Ab9314]] ([[User talk:Ab9314|talk]]) 11:11, 24 March 2017 (UTC)

== Introduction ==

Molecular dynamics (MD) is a method to simulate the motion of many particles by numerically integrating Newton's Laws of Motion over time. The simulation can be conducted in an ensemble, e.g. NVT or NVE, by setting limits on temperatures and pressures (thermostatting and barostatting respectively).<ref name=sutmann>{{cite journal | author1=Godehard Suttman | year=2002 | title=Classical Molecular Dynamics| journal=NIC Series | volume=10 | ISBN=3-00-009057-6 | pages=211-254 | doi=|10.1016/0364-5916(91)90030-N}}></ref> MD can be used to simulate the behaviour of protein folding, lipid bilayers, and other biological systems, thereby making this a very practical and convenient way to experiment with complex systems.<ref name=biological_apps>{{cite journal | author1=Karplus, Martin | year=2002 | title=Molecular dynamics simulations of biomolecules| journal=Nature Structural Biology | volume=9 | pages=646-652 | doi=10.1038/nsb0902-646}}</ref> The potential between particles were calculated using pairwise-additive 12-6 Lennard-Jones potentials. All simulations were performed on the HPC system using the [http://lammps.sandia.gov LAMMPS software package].

=== Velocity Verlet Algorithm ===

To begin the simulation, a lattice of atoms can be initiated with atoms a user-defined length apart (inputted in reduced units). All particles were given an initial velocity (according to Maxwell-Boltzmann distribution). Integrating each particle's velocity over time iteratively gave the new positions of the atoms in the system. This was done using the velocity-Verlet algorithm.

At first, the timestep of the simulation was varied to observe the errors in the positions with passing time when using the velocity-Verlet algorithm. This was done by calculating the exact position of a simple harmonic oscillator, and comparing against the position given by the velocity-Verlet calculation.

For a simple harmonic oscillator, <math>x(t) = A \cos(\omega t + \phi)</math>. Using an Excel spreadsheet, the error in the Velocity-Verlet algorithm was calculated. In this case, <math> A = \omega = 1</math> and <math> \phi = 0 </math>. The error is plotted in Figure 1, and the maxima in the error is plotted against time in Figure 2.

<gallery mode=packed heights=205px class="center">
image:sho_error.PNG|Figure 1: Error in displacement (when calculated by velocity-Verlet algorithm) against time for a simple harmonic oscillator.
image:sho_max_error.PNG|Figure 2: Maximum error in displacement against time for the harmonic oscillator.
</gallery>

By changing the timestep in the spreadsheet, it was found that the error in <math>x(t)</math> as calculated by the velocity-Verlet algorithm is less than <math>1%</math> over the course of the simulation when the timestep is less than <math>0.025</math>. 

The total energy of the system should be monitored to make sure that the energy is conserved in the system, and to check if the system is in thermal equilibrium. If the total energy is not stable, then the thermodynamic values derived from the system will not be equilibrium values, and cannot be compared to other MD systems.

'''<span style="color:red">JC: Good choice of timestep, the total energy should be approximately constant to ensure that the simulations give physically meaningful results.  Why does the error oscillate over time?</span>'''

===Atomic Forces===

In order to make the calculations easier, the LJ-potentials were given a 'cut-off' distance - any atoms separated by a distance larger than this truncation distance will not be included in the LJ force calculations.<ref name=RDF>{{cite journal | author1=Allen, M.P. | author2=Tildesley, D.J. | year=1987 | title=Computer Simulation of Liquids | volume=1 | journal=Oxford: Clarenden Press | ISBN=0-19-855375-7 | pages=27-29}}</ref>

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

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

<math>\int \phi \left ( r \right ) {d}r = 4 \left ( - \frac{1}{12r^{11}} + \frac{1}{6r^5} \right ) + c = \frac{2}{3r^5} - \frac{1}{3r^{11}} + c</math>

Thus, when limits are introduced:

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

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

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

From the above calculations, the force of attraction reduces very quickly when the distance from the Lennard-Jones centre is greater than <math>2\sigma</math>. Hence, truncation at such a distance will introduce only a very small error in the system.

'''<span style="color:red">JC: You are missing half of the maths questions - what is r0, the force at r0, req and the potential at req?  The integral is not correct, the integral of 1/(r^12) is -12/(r^11), not -1/(12 r^11).</span>'''

===Periodic Boundary Conditions===

As it is very heavy computationally to simulate even millimole amounts in LAMMPS, Periodic Boundary Conditions (PBCs) are used to cut down the calculation times.<ref name=RDF>{{cite journal | author1=Allen, M.P. | author2=Tildesley, D.J. | year=1987 | title=Computer Simulation of Liquids | volume=1| journal=Oxford: Clarenden Press | ISBN=0-19-855375-7 | pages=24-27}}</ref> The reason for this difficulty can be seen by considering the number of particles in a small volume of water.

For example, 1 mL of water contains:
<math>\frac{1}{18} \cdot 6.022 \cdot 10^{23} = 3.346 \cdot 10^{22}</math> molecules.

In most cases, the number of particles used will be between 10<sup>3</sup> and 10<sup>4</sup>.

10<sup>4</sup> water molecules take up a volume of: <math>\frac{10^4}{6.022 \cdot 10^{23}} \cdot 18 = 2.99 \cdot 10^{-19} \, mL</math>

PBCs allow for the relatively small number of particles in our simulation to act as though it's part of a bulk solution. To this end, the particles are simulated such that when moves off one edge of the box, it reappears on the opposite face with the same velocity. Hence, the number of particles in the box, N, is always conserved.

For example, consider an atom at <math>(0.5, 0.5, 0.5)</math> in a cubic box with vertices at <math>(0, 0, 0)</math> and <math>(1, 1, 1)</math>. If this atom moves by vector <math>(0.7, 0.6, 0.2)</math> in a single timestep, then the final position of the particle will be: 
<math>(0.2, 0.1, 0.7)</math>.

===Reduced Units===

Reduced units are used to make the numbers more manageable.

For Argon, <math>\sigma = 0.34 \, nm</math>. Given that the LJ-cutoff in reduced units is at <math>r^*=3.2</math>, the cutoff in nanometers is: <math>r = r^* \, \sigma = 3.2 \cdot 0.34 = 1.088 \, nm</math>.

Given <math>\frac{\epsilon}{k_B} = 120</math>,

well depth, <math> \epsilon = 120 k_B = 1.66 \cdot 10^{-21} \, J = \frac{1.66 \cdot 10^{-21} \cdot N_A}{1000} kJ \, mol^{-1} = 9.98 \cdot 10^{-1} \, kJ \, mol^{-1}</math>.

Therefore, when <math>T^* = 1.5</math>, and given <math>\frac{\epsilon}{k_B} = 120</math>,

<math>T = \frac{\epsilon}{k_B} \cdot T^* = 120 \cdot 1.5 = 180 \, K</math>.

'''<span style="color:red">JC: Calculations correct.  PBCs are used to remove the effects of the edges of the simulation box and make simulations a small systems behave like a bulk systems.</span>'''

==Equilibriation==

===The Simulation Box===

The particles for the simulation were initiated in a simple cubic lattice structure. A random set of coordinates were not used. This was to avoid two atoms being generated too close together - this would lead to a very large potential, and cause the pair to split violently apart. This simulation would be unlikely to be a good representation of a real life scenario of a crystal melting.

'''<span style="color:red">JC: High repulsive forces make the simulation unstable and can cause it to crash unless a much smaller timestep is used.</span>'''

In this simulation, a lattice density of 0.8 was used (0.8 lattice points per unit volume). If a face-centred cubic (fcc) lattice had been used instead, with a lattice density of 1.2, then the side length of the cubic unit cell, <math>a</math>, can be found thus:

<math>\rho = \frac{N}{V}</math>

<math>1.2 = \frac{4 \cdot 1000}{(10 \cdot a)^3}</math>

<math>a = \frac{\sqrt[3]{\frac{4000}{1.2}}}{10} = 1.49380</math>

'''<span style="color:red">JC: Correct.  How many more atoms would be created in an fcc lattice compared to a simple cubic lattice?</span>'''

===Properties of atoms===

The reaction conditions such as the lattice density are passed into LAMMPS using text commands in the input file. Part of the commands are shown below:

<pre>mass 1 1.0 </pre>

<code> mass</code> defines the property of the atom that is being defined. <code>1</code> indicates the type of particle being defined; all particles in LAMMPS are given a type, in form of an integer number. For example, if considering a mixture with 2 types of atoms, then the particles in LAMMPS would be differentiated by having type <code>1</code> and type <code>2</code> particles, and so on.

<code>1.0</code> indicates the [http://http://lammps.sandia.gov/doc/mass.html  mass] of the particle in g mol<sup>-1</sup>.


<pre>pair_style lj/cut 3.0</pre>

This line defines the [http://lammps.sandia.gov/doc/pair_style.html potential between two interacting particles]. Once again, the <code>pair_style</code> keyword indicates the property being defined. <code>lj/cut</code> defines the potential as a truncated Lennard-Jones potential, while the <code>3.0</code> defines the potential to be truncated at <math>3.0 \, \sigma</math>. Overall, the potential is being defined as:

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


<pre>pair_coeff * * 1.0 1.0</pre>

This line is defining the [http://lammps.sandia.gov/doc/pair_coeff.html force that each particle exerts] on all other particles is the same as that defined by the negative of the derivative of the Lennard-Jones potential.

'''<span style="color:red">JC: What do the asterisks mean in this command and what do the number refer to for the Lennard-Jones potential?  Why is a cutoff used for the Lennard-Jones potential?</span>'''

The generated particles are then given initial coordinates, and their initial velocities. As <math> x_i(0)</math> and <math>v_i(0)</math> are both specified, the velocity-Verlet integration is used for the simulations.

===Running the Simulation===

The final lines of the input file are as follows:
<pre>variable timestep equal 0.001

variable n_steps equal floor(100/${timestep})

### RUN SIMULATION ###
run ${n_steps}</pre>

<code>variable timestep equal 0.001 </code> creates a variable so that any subsequent occurence of the string <code>${timestep}</code> is replaced by <code>0.001</code>. The number of steps the simulation is run for is also determined by the value of the timestep. By assigning a variable name to the timestep, editing the timestep value will automatically adjust the number of steps, avoiding tedious editing of the value in two different places each time the timestep is changed.

'''<span style="color:red">JC: Correct, using variables makes it easier to run different simulations with different valued parameters using the same script.</span>'''

===Checking Equilibriation===

Several different simulations were run using different timesteps, and the total energy of the system was monitored at each step. Figure 3 shows the results of the simulations.

[[File:Engvstime.jpg|center|frame|none|alt=Energy vs. Time|Figure 3: Energy vs. Time for differing timesteps]]

All systems reach equilibrium, except for the system where the timestep is 0.015. The equilibrium was reached by around <math> time \approx 0.6</math> - in quite a short amount of time.

Of the five timesteps used, the largest acceptable timestep is 0.01. A particularly bad choice would be a timestep of 0.015. This is because the total energy of the system seems to be increasing indefinitely with increasing time with large fluctuations. This is a poor representation of a real system, hence this simulation timestep length is too long for this system.

'''<span style="color:red">JC: The average total energy should not depend on the choice of timestep so timesteps of 0.01 and 0.0075 are not suitable.  Both 0.0025 and 0.001 give roughly the same average total energy so the best choice would be the larger of these two (0.0025) to enable longer simulations to be run more efficiently.</span>'''

==NPT Simulation==

This section looks at running simulation under specific conditions, using an NPT ensemble. This is done using barostats and thermostats.

===Temperature and Pressure Control===
Temperature, <math>T</math>, fluctuates in the simulation, and the target temperature is <math>\mathfrak{T}</math>. As temperature is related to the total kinetic energy of particles within the system, the temperature is kept roughly constant by multiplying the velocities of the particles by a factor <math>\gamma</math>. Hence, these two equations can be surmised from the above, from which an expression of <math>\gamma</math> is derived:

<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} 
\implies \frac{1}{2} \gamma^2 \sum_i m_i v_i^2 = \frac{3}{2} N k_B \mathfrak{T} </math>

<math>\frac{1}{2} \gamma^2 \cdot \frac{3}{2} N k_B T = \frac{3}{2} N k_B \mathfrak{T} </math>

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

'''<span style="color:red">JC: Nearly, there shouldn't be a factor of 2 in the square root, check your working again.</span>'''

===The Input Script===

The thermodynamic property values of the system were periodically computed and saved to file using the following command in the input script:

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

These values are [http://lammps.sandia.gov/doc/fix_ave_time.html time-averaged] in a specific way. The arguments <code>100 1000 100000</code> instruct LAMMPS to compute the properties every 100000 steps using preceding 1000 values in intervals of 100 values. For example, at the 200000th step, the properties are averaged using values at steps number 100100, 100200, 100300... 199000, 200000.

'''<span style="color:red">JC: Correct.</span>'''

===Plotting the Equations of State===

Simulations of a crystal melting were run at two different pressures, P = 2.5 and P = 3.0, with 8 temperatures for each pressure. A graph of density against temperature was then plotted, along with the density predicted by the Ideal Gas Law with respect to temperature at these two pressures. The results are shown in Figure 4. It can be seen that in general, the MD simulation predicts a lower density than the ideal gas law. This is because the ideal gas law does not take the forces between particles into account. In the Lennard-Jones potential model, the repulsive forces are very strong (compared to the attractive forces). Hence, the MD simulations show the particles will prefer to be more spread out (i.e. less dense) so as to minimise the unfavourable repulsive forces. The discrepancy decreases with increasing temperature, because the particles in the MD systems will move further apart. This in turn decreases the magnitude of the forces between particles, bringing it closer to the ideal system of no forces between particles.


[[Image:Pressure vs temperature.png|thumb|750px|center|none|alt=Density vs. Temperature|Figure 4: Density vs. Temperature for varying pressure as found by MD simulations, compared to results predicted by the Ideal Gas Law.]]

'''<span style="color:red">JC: Correct, how does the discrepancy between the simulation and ideal gas results change with pressure?  The ideal gas is a good approximation to a dilute (high temperature, low pressure) gas.</span>'''

==Heat Capacity Calculations using Statistical Physics==

To investigate the effect of changing temperature on the heat capacity, a simulation was run in the NVT ensemble for two different densities (<math>\rho = 0.2</math> and <math>\rho = 0.8</math>). The results of the simulations are summarised in Figure 5.


[[File:Cv_v_vs_temp.png|thumb|750px|center|none|alt=Heat Capacity per Unit Volume vs. Reduced Temperature|Figure 5: Heat Capacity per Unit Volume vs. Reduced Temperature for differing lattice densities.]]

The expected value for <math>\frac{C_V}{V}</math> in an ideal system is:

<math>\frac{C_V}{V} = \frac{3}{2} \cdot N \cdot \frac{1}{V} = \frac{3}{2} \cdot \rho</math>

Hence, for an ideal gas, when <math>\rho = 0.3, \, \frac{C_V}{V} = 0.1</math> and when <math>\rho = 0.8, \, \frac{C_V}{V} = 1.2</math>.

The value of <math>C_V</math> found in the simulations are around 4 orders of magnitude smaller than what is predicted the ideal system, even at low densities (more ideal system). Also, the simulation shows variation of <math>C_V</math> with increasing temperature (approximately <math>C_V \propto \frac{1}{T^n}, \, \, \, n < 2</math>).<ref name=heat_capacity_data>{{cite journal | author1=Dinsdale, A.T. | year=1991 | title= SGTE data for pure elements| journal=CALPHAD | volume=15 | issue=4 | pages=317-425 | doi=|10.1016/0364-5916(91)90030-N}}</ref> The value of n is unknown, and must be found by fitting different curves to the data. We know that <math>C_V \propto \frac{Var(E)}{T^2}</math> and that <math>Var(E) \propto T^m, \, m > 0</math>. Hence, we can estimate that <math>n = 2 - m \implies n < 2</math>. This is in contrast to the ideal system, which predicts no variation of heat capacity with temperature. 

However, qualitatively, the trends seen are reasonable. As the temperature is raised, the particles have more and more energy, yet the number of available degrees of freedom fall. As <math>C_V</math> is a measure of how easily the system is excited to higher energy states - an already excited system requires even more energy to put the particles in a higher vibrational energy states. Increasing the number of particles, <math>N</math>, also leads to a larger number of degrees of freedom, leading to an increase in the heat capacity, <math>C_V</math>. Consequently, increasing <math>\rho = \frac{N}{V}</math> also leads to an increase in <math>\frac{C_V}{V}</math>. This trend is clearly seen in Figure 3.

Finally, the apparent reversal in trends from <math>T^* = 2.3</math> to <math>T^* = 2.4</math> is showing the change from a liquid to a gaseous phase.

'''<span style="color:red">JC: Some good suggestions.  Why do you think the number of degrees of freedom reduces as energy is increased?  Remember these simulations are classical simulations of spherical particles, there are no vibrational energy levels.</span>'''

===The Input Script===

One of the scripts used for running the NVT simulations is shown below:

<pre>
### CREATING VARIABLE TO DEFINE LATTICE DENSITY
variable rho equal 0.8

### DEFINE SIMULATION BOX GEOMETRY ###
lattice sc ${rho}
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.2
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
variable pdamp equal ${timestep}*1000
fix nvt all nvt temp ${T} ${T} ${tdamp}
run 10000
reset_timestep 0

### ENDING NVT ONCE SYSTEM IN REQUIRED STATE. NVE RESTATED.
unfix nvt
fix nve all nve

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp press density
variable temp equal temp
variable temp2 equal temp*temp
variable eng equal etotal
variable eng2 equal etotal*etotal
fix aves all ave/time 100 1000 100000 v_temp v_temp2 v_eng v_eng2
run 100000

variable avetemp equal f_aves[2]
variable errtemp equal sqrt(f_aves[5]-f_aves[2]*f_aves[2])

### CALCULATION OF C_V
variable vareng equal (f_aves[8]-f_aves[7]*f_aves[7])
variable cv equal (3375*3375*${vareng})/(${avetemp}*${avetemp})

### PRINTING CALCULATED VALUES TO OUTPUT FILE
print "Averages"
print "--------"
print "Temperature: ${avetemp}"
print "Stderr: ${errtemp}"
print "Heat Capacity: ${cv}"
</pre>

==Structural Properties and the Radial Distribution Function==

The Radial Distribution Function (RDF, or <math>g(r)</math>) can provide some insight into the structure of the particles in the simulation. The three phases, solid, liquid and gas, have distinctive RDF plots. These RDF plots can be compared directly to X-Ray diffraction plots produced experimentally, so that the identity of the systems can be verified by this property.<ref name=RDF>{{cite journal | author1=Allen, M.P. | author2=Tildesley, D.J. | year=1987 | title=Computer Simulation of Liquids | volume=1| journal=Oxford: Clarenden Press | ISBN=0-19-855375-7 | pages=54-55}}</ref> Figure 6 shows the RDF plots acquired by varying the density and temperature of the NVT simulation.

<gallery mode=packed heights=395px class="center">
image:RDF_allStates.png|Figure 6: The RDF of the same particles in solid, liquid and gaseous phase. The initial part of the graph where <math>RDF=0</math> indicates space between adjacent particles. This figure also shows that <math>g(r)</math> tends to <math>1</math> as <math>r</math> becomes large.
image:RDF_integrated_allStates.png|Figure 7: The integral of the RDF of the same particles in solid, liquid and gaseous phase.
</gallery>

[[File:fcc_lattice.png|left|frame|alt=3 closest neighbours of any atom on an fcc lattice.|Figure 8: Three closest neighbours of each atom on an face-centred lattice shown by vectors '''1''', '''2''' and '''3''' in increasing order of magnitude.]]

As can be seen from Figure 6, from <math>\sigma=0</math> to <math>\sigma \approx 1 </math> is where the repulsive forces between the adjacent atoms are too great to be overcome, and hence is empty space.

The RDF of solids have sharp peaks at well defined spacing. The atoms responsible for each of the first three main peaks are outlined in Table 1 using Figure 8 as a reference. The lattice parameter is the distance between two adjacent atoms on the fcc lattice along a straight edge (i.e. not diagonally adjacent). This appears to be approximately <math>1.775 \sigma</math>. The RDF decreases slowly as the radius increases.

By using the plot of <math>\int g(r) \, \mathrm{d}V</math> against <math>r</math> (see Figure 7), the coordination number around a single reference atom as the distance, <math>r</math>, increases. The coordination numbers at symmetry related lattice vectors '''1''', '''2''' and '''3''' for this fcc lattice are tabulated in Table 1.

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: 20px;"
|+ Table 1: Summary of data inferred from RDF data.
! Vector !! Distance from central atom/<math>\sigma</math> !! Coordination Number
|-
| 1 || 1.025 || 8
|-
| 2 || 1.775 || 8
|-
| 3 || 2.725 || 16
|}

The liquid has a smooth RDF, with the peaks similar distances apart, but its decay is faster than that of the solid's. The progressively smaller peaks indicate the coordination shells around each atom. 

Finally, the RDF of the gas phase tends to 1 very quickly, with no visible oscillations about <math> g(r) = 1 </math>. This indicates a uniform distribution of particles in the simulation box, as is expected of a gas in equilibrium.

'''<span style="color:red">JC: Solid has long range order, liquid only has short range order.  How did you calculate the lattice parameter, it should be roughly the distance to the second peak which is less than 1.775.  Could you have also calculated it from the first and third peaks and taken an average?  The coordination numbers for the first 3 peaks should be 12, 6 and 24, did you zoom in to the integral graph to get these values?</span>'''

==Dynamical properties and the Diffusion Coefficient==

===Mean Square Displacement (MSD)===
[[File: MSD illustration.PNG|left|thumb|200px|Figure 9: A particle moving from initial position '''r(0)''' at time <math>0</math> to another position '''r(t)''' at time <math>t</math>.]]

The motion of the particles in the molecular simulation were investigated by finding several properties of the systems.

The [http://isaacs.sourceforge.net/phys/msd.html Mean Square Deviation (MSD)] of the system is a measure of how freely the particles in the system can move to a new position (see Figure 9). The MSD is calculated thus:

<math>MSD = [\mathbf{r(0)} - \mathbf{r(t)}]^2</math>

Hence, as time progresses in a system with diffusion, the sum of MSD over all time intervals will be increasing. As <math>t \rightarrow \infty</math>, the rate of increase of MSD is approximately constant. The diffusion coefficient, <math>D</math>, is defined as:

<math> D = \frac{1}{6} \, \frac{\partial \langle \mathbf{r}^2\mathbf{(t)} \rangle }{\partial t} </math>

This equation is true as <math>t \rightarrow \infty</math>; at small values of <math>t</math>, there are few or no collisions, and the increase in <math>MSD \propto t^2</math>. <math>D</math> can thus be estimated by fitting a straight line to the linear portion of a plot of <math>MSD</math> against time. Figures 10, 11, 12 and 13 show the calculation of the diffusion coefficient from the <math>MSD</math> plots. Table 2 records the values of <math>D</math> - a clear trend of increasing values of diffusion coefficients with decreasing density can be seen from the data. Atoms in gases can therefore diffuse 20 times as fast as atoms in solids.

<gallery mode=packed heights=360px class="center">
image:allStates_MSD.png|Figure 10: The MSD of the particles at at all states. The values increases with relation to time - the fastest increase is for the gas phase, followed by the liquid phase, and the solid phase is the slowest.
image:Solid msd fitted.png|Figure 11: The MSD of the particles in the solid phase. As the particles' movements are restricted and not averaged, the plot is not smooth.
image:liquid_msd_fitted.png|Figure 12: The MSD of the particles in the liquid phase. The plot is linear, indicating <math>MSD \propto t</math>.
image:gas_msd_fitted.png|Figure 13: The MSD of the particles in the gas phase. Initially, <math>MSD \propto t^2</math>, but as <math>t \rightarrow \inf, \, MSD \propto t</math>.
</gallery>

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: 20px;"1
|+ Table 2: Diffusion coefficient data for solid, liquid and gaseous states.
! State !! <math>\rho^*</math> !! Diffusion Coefficients, <math>D/\sigma^2 \, timestep^{-1}</math> 
|-
| Solid || 1.3 || style="text-align: center;"| <math>2.74 \cdot 10^{-6}</math>
|-
| Liquid || 0.8 || style="text-align: center;"| <math>3.47 \cdot 10^{-4}</math>
|-
| Gas || 0.1 || style="text-align: center;"| <math>5.23 \cdot 10^{-3}</math>
|}

'''<span style="color:red">JC: Would you expect the MSD to increase for a solid?  Did you remember to change the lattice type to fcc in the script?  Good explanation that MSD is proportional to t^2 at short times, this is especially clear in the gas phase plot.</span>'''

The above simulations were performed with thousands of particles, and this can be considered as a relatively small amount - not even a fraction of a femtomole of atoms. To simulate a larger system, a set of simulations were run with 1 million atoms. The difference in the MSD plots (see Figures 14, 15, 16 and 17) is most marked in the solid state - whereas in the smaller system, there was a definite increase in MSD over time for all timesteps, the larger system oscillates slowly oscillates around <math>MSD = 0.0066</math>. Physically, this observation shows that most of the particles in the solid are vibrating and colliding into each other very often. In the smaller system, the collisions are much less frequent, so that <math>|\mathbf{r(0)} - \mathbf{r(t)}|</math> is larger for each atom. The diffusion coefficients of the million atoms in the three states are tabulated in Table 3.

<gallery mode=packed heights=360px class="center">
image:msd_mill_allStates.png|Figure 14: The MSD of the 1 million particles at at all states.
image:msd_mill_sol.png|Figure 15: The MSD of the particles in the solid phase. This plot is jagged, and the gradient is almost <math>0</math>: indicating virtually no diffusion in solid state.
image:msd_mill_liq.png|Figure 16: The MSD of the particles in the liquid phase. The plot is linear from <math>t = 1000</math>.
image:msd_mill_vap.png|Figure 17: The MSD of the particles in the gaseous phase. The plot is linear from <math>t = 3000</math>.
</gallery>

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: 20px;"1
|+ Table 3: Diffusion coefficient data for 1 million atoms in solid, liquid and vapour states.
! State !! Diffusion Coefficients, <math>D/\sigma^2 \, timestep^{-1}</math> 
|-
| Solid || style="text-align: center;"| <math>2.37 \cdot 10^{-8}</math>
|-
| Liquid || style="text-align: center;"| <math>3.50 \cdot 10^{-4}</math>
|-
| Gas || style="text-align: center;"| <math>1.25 \cdot 10^{-2}</math>
|}

'''<span style="color:red">JC: Figure 15 is what you would expect for a solid phase.</span>'''

===Velocity Autocorrelation Function===

The diffusion coefficient, <math>D</math>, can also be estimated from the [http://www.compsoc.man.ac.uk/~lucky/Democritus/Theory/vaf.html velocity autocorrelation function (VACF)], <math>C(\tau)</math>, if <math>C(\tau) \rightarrow 0</math> at large <math>\tau</math><ref name=RDF>{{cite journal | author1=Allen, M.P. | author2=Tildesley, D.J. | year=1987 | title=Computer Simulation of Liquids | volume=1| journal=Oxford: Clarenden Press | ISBN=0-19-855375-7 | pages=59-60}}</ref>:

<math>C(\tau) = \langle \mathbf{v(t)} \cdot \mathbf{v(t + \tau)} \rangle</math>

<math>D = \frac{1}{3} \, \int_0^\infty \langle \mathbf{v(0)} \cdot \mathbf{v(\tau)} \rangle \, \mathrm{d} \tau</math>

Below, the normalised velocity autocorrelation function is evaluated for a 1D harmonic oscillator.

Given <math>C(\tau) = \frac{\int_{-\infty}^{\infty} v\left(t\right)v\left(t + \tau\right)\mathrm{d}t}{\int_{-\infty}^{\infty} v^2(t)\mathrm{d}t}</math>, and <math>x(t)=A \cos (\omega t + \phi) \implies v(t)=-A \omega \sin (\omega t + \phi)</math>:

<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^2 \omega^2 \sin^2 (\omega t + \phi) \mathrm{d} t}</math>

Simplifying,

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

Using the trigonometric identity <math>\sin(a + b) \equiv \sin(a) \cos(b) + \sin(b) \cos(a) </math>, the following can be surmised:

<math>\sin [\omega(t + \tau) + \phi)] = \sin(\omega t + \omega \tau + \phi) = \sin(\omega t + \phi) \cos(\omega \tau) + \cos(\omega t + \phi) \sin(\omega \tau)</math>

Hence,

<math>C(\tau) = \frac{\int_{-\infty}^{\infty} \sin^2 (\omega t + \phi) \cos(\omega \tau) + \sin(\omega t + \phi) \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{\cos(\omega \tau)\int_{-\infty}^{\infty} \sin^2 (\omega t + \phi) + \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>

Rearranging <math>\sin^2(x) \equiv \frac{ 1 - \cos(2x)}{2}</math>,

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

Simplifying,

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

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

However, <math>\cos(-a) = \cos(a)</math> as this is an even function.

Thus, <math>\cos(\infty) - \cos(-\infty) = 0</math>:

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

'''<span style="color:red">JC: Correct result, the derivation could have been simplified by realising that the function sin(x)*cos(x) is an odd function and so the integral is zero.</span>'''

This result shows that for a simple harmonic oscillator, the VACF is a sinusoidal function that depends on the time passed, as well as the frequency of the oscillator. 

The VACF of the solid and liquid states were calculated during the simulations. In addition, the analytical solution to the VACF of a particle doing simple harmonic motion with frequency, <math>\omega = \frac{1}{2 \pi}</math> was also calculated. These are all plotted in Figure 18. 

<gallery mode=packed heights=395px class="center">
image:vacf_allStates.png|Figure 18: The plots of velocity autocorrelation functions with for atoms in a solid, a liquid, and for a simple harmonic oscillator. At large <math>t</math>, <math>C(\tau) \rightarrow 0</math> for the solid and liquid systems.
image:vacf_gas.png|Figure 19: VACF of atoms in the gaseous phase decay much more slowly owing to the much weaker interatomic forces.
</gallery>


From Figure 18, it can be observed that the VACF of the simple harmonic oscillator is sinusoidal, whereas the liquid and solid VACF are best described as damped sinusoidal. The solid's VACF tends to 0 faster that that of the liquid. 

The global minima in the plots for the solid and liquid indicate when all the particles have had at least one collision. By this point, each of the atoms have a different velocity vector associated with it. Hence, the VACF rebounds towards 0. As the simple harmonic oscillator is a lone particle's motion, with no interactions with any other atoms, hence its VACF is very different to the Lennard-Jones solid and liquid. Its velocity vector also changes in a cyclic manner, and does not change from one oscillation to another. Hence, the amplitude of the VACF does not change as time increases.

In the solid, the particles are in close contact due to the high density - the atoms are arranged so that the repulsive forces are minimised, and the attractive forces maximised. This is the distance at which the LJ-potential is at a minimum. Here, the atoms are very stable, and are reluctant to leave the locations they are at.Their motion is therefore a vibration around their eqiulibrium positions, with a reversal of velocity at the end of each oscillation. This explains the oscillatory behaviour of the VACF. However, the damping is explained by perturbative forces acting on the atoms that disrupt the oscillations. Liquids are similar to solids - however, the density here is lower, and diffusion is faster. Here, the diffusive forces are the perturbative forces, and hence the VACF decays much more rapidly than the solid's VACF.

In contrast, the VACF from a system in the gaseous state (see Figure 19) decays much more slowly. The large separation distances between particles in the gaseous phase mean the forces are weaker. Hence, the inter=particle forces between them only affect their velocities very slowly, over a long time. For this reason, the VACF decays much more slowly than the liquid and the solid VACF.

'''<span style="color:red">JC: What are the perturbative forces in the solid?  Collisions between particles randomise particle velocities and cause the VACF to decorrelate, but there are no collisions in the harmonic oscillator.</span>'''

<gallery mode=packed heights=395px class="center">
image:vacf_trapz.png|Figure 20: Cumulative integral of VACF against time for all phases as an approximation to the diffusion coefficient.
image:Vacf_trapz_liq_sol.png|Figure 21: Cumulative integral of VACF against time for solid and liquid phases.
</gallery>

The area under the VACF, and thus the diffusion coefficient, <math>D</math>, was estimated by using the trapezium rule, by using the following formula:

<math>D = \frac{1}{3} \, \int_0^\infty \langle \mathbf{v(0)} \cdot \mathbf{v(\tau)} \rangle \, \mathrm{d} \tau</math>

The cumulative integral is plotted in Figures 18 and 19. The gas has a very large <math>D</math>, compared to the liquid and the solid; the diffusion coefficient of the liquid is approximately 100 times as large as that of the solid. The trend of lower <math>D</math> with higher density is as found earlier in Tables 2 and 3. 

<gallery mode=packed heights=395px class="center">
image:int_vacf_mill_allStates.png|Figure 22: Cumulative integral of VACF against time for all phases for 1 million atoms.
image:Int_vacf_mill_liq_sol.png|Figure 23: Cumulative integral of VACF against time for solid and liquid phases for 1 million atoms.
</gallery>

{| class="wikitable" style="margin-left: auto; margin-right: auto; border: 20px;"1
|+ Table 4: Diffusion coefficient data derived from MSD and VACD compared for solid, liquid and gaseous states.
! State !! Diffusion Coefficients (MSD), <math>D/\sigma^2 \, timestep^{-1}</math> !! Diffusion Coefficients (VACD), <math>D/\sigma^2 \, timestep^{-1}</math>
|-
| Solid || style="text-align: center;"| <math>2.37 \cdot 10^{-8}</math> || style="text-align: center;"| <math>0.0487</math>
|-
| Liquid || style="text-align: center;"| <math>3.50 \cdot 10^{-4}</math> || style="text-align: center;"| <math>15.1 </math>
|-
| Gas || style="text-align: center;"| <math>1.25 \cdot 10^{-2}</math> || style="text-align: center;"| <math>544</math>
|}

A larger system with a million atoms was also simulated to simulate a system more akin to a real-life system. These plots (see Figures 22 and 23) are much smoother and give a more accurate value of the diffusion coefficient, <math>D</math>. The values of <math>D</math> acquired for the 1 million atom simulation from the VACD is compared to that of the values acquired from the MSD in Table 4. The table shows <math>D</math> values gotten from the VACD values are much higher than those derived from the MSD values. The largest source of errors of calculation of <math>D</math> from VACF is from the use of trapezium rule to approximate the area under the graph.

'''<span style="color:red">JC: Calculating D from the VACF relies on the integral of the VACF reaching a plateau within the simulation time, is this the case here?  Do you think the MSD or VACF derived values for the diffusion coefficient are more accurate?</span>'''

==References==
{{reflist}}