View source for Talk:Efr114ls

'''<span style="color:red">JC: General comments: Tasks which were answered were answered well with good explanations, but the questions on the RDF and VACF were missing.  Make sure you attempt all of the tasks before spending too long on your answers to the earlier tasks as the later sections are more useful and more interesting.</span>'''

__TOC__
== Introduction to Molecular Dynamics Simulation==
[[File:Initial graphsefr114.PNG|thumb|Figure 1a) A graph showing how position varies with time for both the Verlocity-Verlet model (blue points) of the harmonic oscillator and the analytical model of the harmonic oscillator (red points)
Figure 1b) A plot displaying how total energy of the system varies with time

Figure 1c) A plot displaying how absolute error varies with time 
]]
The file 'HO.xls' contains a Velocity-Verlet model for the behaviour of a classical mechanic oscillator, this file was then expanded to include:
* Values for position at time <math> t </math> using the classical mechanics model 
* Values for absolute difference between the classical mechanics model and the Velocity-Verlet model
* Values for total energy of the system 
This processed file can be accessed here: [[File:HOprocessedefr114.xls]]

[[File:Bloopblood.PNG|centre|thumb|Figure 2: A plot showing the relationship between maxima of absolute error against time.]]

The maximum deviation from the mean total energy, measured in percentage, was calculated for 20 timestep values ranging from 0.025 to 0.4. A graph and trend line were then plotted,seen in Fig.3 to produce a quadratic function and extrapolated for 1% deviation. The value of timestep used in further simulations must not exceed 0.282 - as above this value the total energy of the system fluctuations exceed the 1% limit.

The total energy of a system is one way of monitoring if the method of simulation is valid. For a simulation to be valid it must obey the 1st law of thermodynamics. Total energy is also a useful property to track the process of equilibration and can be compared with experimental data and literature to determine the accury of the simulation method.

[[File:Errorefr114section1.PNG|centre|thumb|Figure 3 A plot showing how maximum deviation from average varies with time ]]

'''<span style="color:red">JC: Why does the error oscillate over time?  Good, thorough analysis of the choice of timestep.</span>'''

=== The Lennard-Jones Model  ===
The Lennard-Jones model is widely used to describe how potential varies with distance between two non-ionic particles. It is not a completely accurate representation for a potential energy surface but its simplicity makes it an essential tool in computational chemistry.

The potential between the two particles is zero when <math> r_0=\sigma </math>, (calculations for this can be seen below). Sigma is the equivalent of half the internuclear distance between the non bonding particles-when the particles are the same type of atom, this value is known as the Van Der Waal's radius.  

The equation for the Lennard-Jones Potential for an interaction between two non ionic bodies is : 
<div style="text-align: center;">
<math> \phi(r) = 4\epsilon\bigg(\frac {\sigma^{12}}{r^{12}} -  \frac {\sigma^{6}}{r^{6}}\bigg) </math>

<math> 0 = 4\epsilon\bigg(\frac {\sigma^{12}}{r^{12}} -  \frac {\sigma^{6}}{r^{6}}\bigg) </math>
<math> r_0=\sigma</math>
</div>
The force due to this interaction is given by the negative of the derivative of potential in respect to to time, resulting in: 
<div style="text-align: center;">
<math> F(r)= -\frac{d\phi(r) }{dr }= 24\epsilon\bigg(\bigg(\frac{\sigma^{6}}{r^{7}}\bigg) - 2\bigg(\frac{\sigma^{12}}{r^{13}}\bigg)\bigg)
</math>

Equilibruim separation, req, is achieved when the sum of the attractive forces F(r) is equal to zero. From this equation, the equilibrium separation can be determined.

<math> F(r)=0 </math>

<math> \frac{\sigma^6}{r_{eq}^{7}}=2\frac{\sigma^12}{r_{eq}^{13}} </math>

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

<math>r_{eq} =\sqrt[6]2</math>
</div>

By substituting the req value back into the Lennard-Jones potential equation, the depth of the potential and hence the energy required to fully dissociate the two atoms can be determined: phi=epsilon

<math> \int_{\infty}^{2\sigma} \phi(r) \, dr = 4\epsilon\sigma\bigg(\frac{1}{5}\bigg(\frac{\sigma^{5}}{r^{5}}\bigg) - \frac{1}{11} \bigg( \frac{\sigma^{11}}{r^{11}}\bigg)\bigg) \bigg |_{\infty}^{2\sigma} = - 2.48*10^{-2}</math> 

<math> \int_{\infty}^{2.5\sigma} \phi(r) \, dr = 4\epsilon\sigma\bigg(\frac{1}{5}\bigg(\frac{\sigma^{5}}{r^{5}}\bigg) - \frac{1}{11} \bigg( \frac{\sigma^{11}}{r^{11}}\bigg)\bigg) \bigg |_{\infty}^{2.5\sigma} = - 8.18*10^{-3} </math>

<math> \int_{\infty}^{3\sigma} \phi(r) \, dr = 4\epsilon\sigma\bigg(\frac{1}{5}\bigg(\frac{\sigma^{5}}{r^{5}}\bigg) - \frac{1}{11} \bigg( \frac{\sigma^{11}}{r^{11}}\bigg)\bigg) \bigg |_{\infty}^{3\sigma} = = - 3.29*10^{-3} </math>

'''<span style="color:red">JC: Correct, the force is at r = sigma so you can simplify the expression a bit more.  Sigma can be thought of as the particle diameter, not the particle radius.</span>'''

=== Periodic Boundary Conditions ===
Under standard conditions, the molar volume of water is <math> 18 mL mol^{-1}</math>, multiplication of this value by <math>N_A</math> produces an estimate for the number of water molecules in a millilitre<math> n = 3.34*10^{22} mL^{-1} </math>. The division of this value by 10,000 gives an approximate value for the volume of 10,000 molecules of water <math> = 2.99*10^{-19} </math>

The purpose of these calculations is to demostrate, one millimetre would corresponds to an extremely large system for calculations. Periodic Boundary conditions are essential for the simulations involved in this experiment, they allow unit cells to simulate a bulk liquid. The unit cells are projected in each direction of the axis and if an atom's trajectory leaves the unit cell, its trajectory will continue through the opposite edge of the trajectory. For example, an atom in a simple cubic unit cell that as the initial coordinates of <math> (0.5, 0.5, 0.5) </math> and trajectory <math>(0.7, 0.6, 0.2) </math> will have the final coordinates of <math> (0.7,0.6,0.2) </math>

=== Reduced units ===
Most software and literature involving Lennard-Jones potential uses reduced units. It reduces values to more managable orders of magnitude.

<math> r*=\frac{r}{\sigma} </math>
<math> E*=\frac{E}{\epsilon} </math>
<math> T*=\frac{K_BT}{\epsilon} </math>

Argon has the properties of <math> \sigma= 0.34 nm </math>, <math> \frac{\epsilon}{k_B}=120 K </math>, a Lennard-Jones cut-off value of <math> r*= 3.2 </math> and in the system <math> T*=1.5 </math>.
In real units, the Lennard jones cut-off value, <math> r= 1.1 nm </math> and the well depth, <math> E_w= 1.0 kJ mol^{-1} </math> with a temperature of <math> 180 K</math>.

'''<span style="color:red">JC: All calculations correct.</span>'''

==Equilibration==

=== Creating the Simulation Box ===
'''Atom postion''' Within the input file for a simulation, the initial
positions for the atoms within the simulation box must be stated. A liquid, by definition, has no long range order so it would be a logical presumption to use randomly generated atomic coordinates, however this method can be problematic as it allows the creation of two atoms in close proximity. The issue involved with two atoms with separation of less than <math> \sigma </math> is that, as seen in the Lennard-Jones plot, the potential between the atoms rapidly increases and tends to infinity. The creation of an intial state containing these high potentials and hence high energies would result in LAMMPS running to error and crashing when the simulation was run, unless an extremely small timestep is used- which limits the range of calculations possible.

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

The set up and conditions of the simulation don’t necessarily have to correspond to  exact equilibrium conditions to occur, the conditions must be in a realistic range to equilibrium conditions, allowing the use of ordered particle coordinates. A primitive cubic lattice and other basic lattice structures may be used to control initial inter-atomic distances, preventing the issues that arise with randomization of coordinates.

In following simulations, a primitive cubic lattice is uses to order the atoms, this lattice contains one lattice point per unit cell. Simple calculations for lattice constant, <math>a</math> (in reduced units) if lattice type and number density  are known.  
<div style="text-align: center;"><math>Number Density=\frac{Number\ of\ lattice\ points}{Volume}</math> <math>0.8=\frac{1}{a^{3}}</math>
</div>
<div style="text-align: center;">
<math>a=\sqrt[3]{1.25} = 1.07722 </math>
</div>

A face centred cubic lattice, however has 4 lattice points per unit cell as seen in Figure BLAH, but the same calculation can still be used to calculate lattive constant
<div style="text-align: center;"><math>Number Density=\frac{Number\ of\ lattice\ points}{Volume}</math> <math>1.2=\frac{4}{a^{3}}</math></div>
<div style="text-align: center;">
<math>a=\sqrt[3]\frac{4}{1.2} = 1.49380 </math>
</div>
'''Number of atoms''' For a box with parameters defined as:
 region box block 0 10 0 10 0 10
 create_box 1 box
This creates a simulation box of the dimensions of 10 lattice spacing in each axis direction which corresponds to 10 unit cells in each direction and hence 1000 unit cells overall. If the lattice type used is primitive cubic with one lattice point per unit cell, 1000 lattice points and hence 1000 atoms would be created in the simulation box. If the lattice type used is face-centred cubic(FCC) with four lattice points per unit cell, 4000 lattice points and hence 4000 atoms would be created in the simulation box.

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

'''Other parameters'''

The command <code> mass [type of atom] [mass of atom] </code> is used for setting the mass for each type of atom. The first value indicates the type of atom and the second value indicates the mass of that type of the atom. In the simulations run in this simulation only one type of atom is used and hence the code <code> mass 1 1.0 </code> is used, however if there were more than one type of atom being used a code such as <code> mass 1 1.0 2 1.5 </code> would be used.

'''<span style="color:red">JC: Correct, for more than one type of atom I think you would have to use two separate mass commands.</span>'''

The command <code> pair_style [type of interaction] [parameter]</code> is used for setting the type of interactions between atoms. The first value defines the type of interaction between atoms such a Lennard-Jones <code>lj</code>, Embedded Atom Model <code>eam</code> and Stillinger-Weber<code>sw</code>and <code>/cut</code> shows that there is a cutoff <code>/coul</code> or the lack of it shows that the presence or the absence of coloumbic interactions. The second value is the cutoff parameter in reduced distance units this value must be large than σ. The command used in the simulations in this experiment is <code> pair_style lj/cut 3.0 </code> determines that the type of interactions are Lennard-Jones with a cutoff of 3.0 reduced distance units and with no columbic interactions or any other interactions.

'''<span style="color:red">JC: Here the cutoff is applied to the Lennard-Jones potential, not a Coloumbic potential, why is a cutoff used?</span>'''

The command <code> pair_coeff [type of atom][type of atom] [coefficient] [coefficient]</code> is used to set the force field coefficients between the type of atoms. The first two values represent the atoms that are interaction and can be either numbers to represent to specify the types of atoms interaction or two asterisks to define all of the types of atoms are interacting. The last two values define the coefficients of these interactions representing <math>\sigma</math> and <math>\epsilon</math> respectively. The code used in this experiment is <code> pair_coeff * * 1.0 1.0 </code> determining that all types of atoms are interacting equally with each other with coefficients of 1.0 resulting in symmetrical interactions.

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

Algorithm In the simulations in this experiment, both the initial position <math>\mathbf{x}_i\left(0\right)</math> and intial velocity<math>\mathbf{v}_i\left(0\right)</math> are specified and hence the Velocity-Verlet intergration algorithm is used. The atoms are positioned in and ordered lattice - the primitive cubic lattice  and the intial velocities are randomized and generated using the Maxwell-Boltzmann distribution

The plots in Figures 4,5 and 6 are plots of reduced energy, pressure and temperature against time for a 0.001 timestep and the system can be seen to equilibrate by time 0.25. After which, the values fluctuate around the mean with deviation of less that 1%.
<div><ul> 
<li style="display: inline-block;"> [[File:Energyefr114.PNG|thumb|none|350px|Figure 4: A plot displaying how the total energy of the system varies with time. dt=0.001]] </li>
<li style="display: inline-block;">[[File:Tempefr114.PNG|thumb|none|350px|Figure 5: A plot displaying how the temperature of the system varies with time. dt=0.001]] </li>
<li style="display: inline-block;"> [[File:Pressefr114.PNG|thumb|none|350px|Figure 6: A plot displaying how the pressure of the system varies with time. dt=0.01]] </li>

</ul></div>Figure 7 shows the plot of variation of total energy against time for 5 different timesteps: 0.001, 0.0025, 0.0075, 0.01 and 0.015. The timestep of 0.015 is two large and causes the total energy to diverge away from equilibrium. The remaining timesteps all converge to equilibrium successfully, however a compromise must be made between accuracy of the equilibrium total energy and length of simulation. Total energy of the system should be independant to timestep and both 0.01 and 0.0075 timesteps result in a equilibrium total energy that is significantly higher than the lower timesteps, suggesting that these timesteps are less accurate. The most suitable timestep would be 0.0025 as it has similar accuracy to 0.001 but has a longer simulation.[[File:Energiesefr114.PNG|none|thumb|350px|Figure 7: A plot that displays the relationship between total energy and time for 5 different timesteps]]

'''<span style="color:red">JC: Good choice of timestep.</span>'''

==Simulations Under Specific Conditions==

'''Changes in Input File'''   

The average kinetic energy of the system, as stated by the equipartition theorem and due to the system having three degrees of freedom (x, y, z), can be given by: 
<div style="text-align: center;">
<math>E_{K} =\frac{1}{2}\sum_{i}m_{i}v_{i}^{2}=\frac{3}{2}Nk_{B}T</math>
</div>


For an NpT system, the temperature must be maintained  at target temperature T. In MD, the kinetic energy and hence the temperature of the system fluctuates throughout the simulation, to prevent this fluctuation from occuring, a factor of GAMMA can be introduced to correct the fluctuation which multiplies the velocity of each of the atoms by a factor to maintain a constant temperature. 

<div style="text-align: center;">
<math>E_{K} =\frac{1}{2}\sum_{i}m_{i}(\gamma v_{i})^{2}=\frac{3}{2}Nk_{B}\mathfrak{T}</math>
</div>
These two equations can then be treated simulaniously to produce the desired scaling constant- allowing the isothermal system:

<div style="text-align: center;">
<math>\gamma^{2} =\frac{\mathfrak{T}}{T}</math>  <math>\gamma =\sqrt\frac{\mathfrak{T}}{T}</math>
</div>

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

The input file for NpT simulations contains a command in the format of <code> fix ID group-ID ave/time Nevery Nrepeat Nfreq value1 value2 ... keyword args ... </code>, this determines the calculations for the thermodynamic property averages.
The <code>Nevery</code> specifies to use input values every N timesteps,.

The <code>Nrepeat</code> specifies that an N number of times is used to input values to calculate averages.

The <code>Nfreq</code> specifies to calculate averages every N timesteps.

For this simulation, Nevery=100, Nrepeat=1000 and Nfreq=10000

In summary,  1000 input values result in the simulation of <math>1*10^5</math> time units in this simulation.

'''<span style="color:red">JC: Make sure that you understand these three values.</span>'''

=== The NpT Simulation ===
In simulations previous to this section, the ensemble being investigated was canonical (N, V, T kept constant) however in this simulation, the isothermal-isobaric ensemble is used (N, p, T kept constant). The timestep used for these simulations was 0.0025 as determined in the previous section. The pressures used were 2.5 and 3.0 in reduced units - decided due to the average temperature of  2.62 reduced unitswhen the system was in equilibrium in the previous simulations. The temperatures used were 1.60, 1.85, 2.10, 2.35 and 2.6 in reduced units - values must be greater than 1.5 as this is the critical temperature.<div>
The ideal gas plots for pressures 2.5 and 3.0, seen in Fig.8 were calculated using the formula <math> \rho*=\frac{p*}{T*} </math>. For both pressures tested in these simulations, the ideal gas plot results in substantially higher values for the simulation for the Lennard-Jones liquid. These differences arise from the particle modeled, both ideal gases and Lennard-Jones liquids are based on broad and contrasting assumptions. An ideal gas is composed of identical particles of negligible volume and have no intermolecular forces between particles- resulting in the total energy of the system is completely from the kinetic energy of the particles and there are no potential energy contributions. In contrast to our simulations which have Lennard-Jones interactions between particles in the system resulting in a non-zero potential.  The discrepancy between the two models decreases and the temperature values increase, this could be due to the fact that at higher temperatures, average kinetic energy of the particles increases  and contributes more to the total energy of the system and hence Lennard-Jones fluid behaves more alike to an ideal gas with an increase of temperature. 

<div><ul> 
<li style="display: inline-block;">[[File:Zoomyzoomefr114.png|thumb|none|450px|Figure 8: A figure displaying the relationship between density and temperature for an LJ fluid including error bars]]</li>
<li style="display: inline-block;">[[File:Nptefr114.png|thumb|none|450px|Figure 9: A figure displaying the relationship between density and temperature for both the LJ fluids and ideal gases.]] </li>
</ul></div>

'''<span style="color:red">JC: Specifically it is the lack of repulsive forces in the ideal gas that gives it a higher density than the Lennard-Jones simulations.  The ideal gas is a good approximation to a dilute gas (high temperature and low density) as in this situation interactions between particles are less significant.</span>'''

==Heat Capacity Calculations==
Heat capacity is the amount of energy required to heat a system by a unit temperature. It can be given using the equation:

<math>C_V=N^{2}\frac{\langle E^{2}\rangle-\langle E\rangle^{2}}{k_BT^{2}} </math>

Heat capacity per Volume can be calculated using the equation:

<math>C_V=N*\rho\frac{\langle E^{2}\rangle-\langle E\rangle^{2}}{k_BT^{2}} </math>

The <math> N^{2} </math> variable here is a correction factor that is required due to the LAMMPS software's method of calculating heat capacity, it also provides an easy conversion of heat capacity to heat capacity per volume due to the <math> \rho=\frac{N}{V} </math> relationship.

The ensemble used in these simulations is the canonical ensemble. In total ten simulations were run with densities of 0.8 or 0.2 in reduced units and temperatures of 2.0, 2.2, 2.4, 2.6 and 2.8 in reduced units. <code>variable heatcap</code> and <code> variable heatcapperv</code> were built into the code to calculate the heat capacity and heat capacity per volume directly. Pressure calculations were removed from the script to reduce simulation times.

<math>C_V</math> was plotted as a function of temperature to produce Fig.10. Increasing the temperature of a system usually increases the specific heat capacity of a system, however as seen in Fig. 10,in this simulation and increase of the temperature of a system results in a decrease of Heat capacity for both densities. In comparison to literature, the trend is most likely due to the Frenkel line being crossed resulting in an dynamic transition from a rigid liquid to a non-rigid gas-like fluid, in this state an increase in temperature  results in lower vibrational energy but higher diffusion rates which results in a reduction of internal energy. The physical definition of heat capacity <math> C_V=\frac{\delta U}{\delta T} </math>, therefore a reduction in internal energy results in a reduction of specific heat capacity. 
 file://icnas4.cc.ic.ac.uk/efr114/ncomms3331.pdf

'''<span style="color:red">JC: Good suggestion for the trend in heat capacity with temperature, more analysis beyond the scope of this experiment would be needed to confirm this trend, why do you think heat capacity usually increases with temperature?  What about the trend with density?  What function did you fit to your results and why?</span>'''

[[File:Spec heatcapefr114.PNG|none|thumb|450px|Figure 10: A graph to show how heat capacity varies with temperature.]]
Below is an intro file used for the heat capacity simulations.
{| class="wikitable"
|### DEFINE SIMULATION BOX GEOMETRY ###
variable dens equal 0.8

lattice sc ${dens}

region box block 0 15 0 15 0 15

create_box 1 box

create_atoms 1 box



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



<nowiki>###</nowiki> SPECIFY THE REQUIRED THERMODYNAMIC STATE ###

variable T equal 2.0

variable timestep equal 0.0025



<nowiki>###</nowiki> ASSIGN ATOMIC VELOCITIES ###

velocity all create ${T} 12345 dist gaussian rot yes mom yes



<nowiki>###</nowiki> SPECIFY ENSEMBLE ###

timestep ${timestep}

fix nve all nve



<nowiki>###</nowiki> THERMODYNAMIC OUTPUT CONTROL ###

thermo_style custom time etotal temp

thermo 10



<nowiki>###</nowiki> RECORD TRAJECTORY ###

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



<nowiki>###</nowiki> RUN SIMULATION TO MELT CRYSTAL ###

run 10000

unfix nve

reset_timestep 0



<nowiki>###</nowiki> BRING SYSTEM TO REQUIRED STATE ###

variable tdamp equal ${timestep}*100

fix nvt all nvt temp ${T} ${T} ${tdamp} 

run 10000

reset_timestep 0



<nowiki>###</nowiki> MEASURE SYSTEM STATE ###

thermo_style custom step etotal temp density

variable etot equal etotal

variable etot2 equal etotal*etotal

variable temp equal temp

fix aves all ave/time 100 1000 100000 v_temp v_etot v_etot2

unfix nvt 

fix nve all nve  

run 100000  

variable averagetemp equal f_aves[1]

variable heatcap equal (atoms*atoms*(f_aves[3]-f_aves[2]*f_aves[2])/(f_aves[1]*f_aves[1]))

variable heatcapperv equal ((atoms*${dens})*(f_aves[3]-f_aves[2]*f_aves[2])/(f_aves[1]*f_aves[1]))





print "Temp: ${averagetemp}"

print "Heat Capacity: ${heatcap}"

print "Heat Capacity per unit volume: ${heatcapperv}"
|}

==Structural Properties and the Radial Distribution Function==
{| class="wikitable"
!Phase
!Temperature (reduced units)
!Density (reduced units)
|-
|Vapour
|1.5
|0.01
|-
|Liquid
|1.25
|0.8
|-
|Solid
|1
|1.2
|}

==Dynamical Properties and the Diffusion coefficient==
The MSD measures the deviation of the particles postion in respect to a standard reference position over time. Simulations for 8,000 particles were run and compared to simulations of 1,000,000 atoms.
As seen below, for both 8,000 and 1,000,000 particles, the solid-state simulations have a clear trend, the MSD rapidly increases until equilibrium is reached, after which there is very minimal variation- this is because solids are not capable of undergoing brownian motion and particle motion is confined. Both the liquid-phase and the vapour-phase increase steadily with time as both fluids and vapours experience brownian motion. The system with 1,000,000 particles gives a more accurate equilibrium with less fluctuation and deviation from the main value, this is due to the large number of particles in the system, resulting in a more narrow gaussian distribution and a better defined value.

{| class="wikitable"
!
!Solid 
!Liquid
!Vapour
|-
|8,000 particles
|[[File:Smallsoledefr114.PNG|frameless]]
|[[File:Smallliqefr114.PNG|frameless]]
|[[File:Smallvapefr114.PNG|frameless]]
|-
|1,000,000 particles
|[[File:Efr114largesolid.PNG|frameless]]
|[[File:Efr114largeliq.PNG|frameless]]
|[[File:Efr114argeMSD.PNG|frameless]]
|}

Diffusion Coefficient(D) 

The diffusive behaviour of the system's that were tested can be simplified into a diffusive coeffictiant, the formula for which is:
<math> D=\frac{1}{6} \frac{\delta \langle r^{2} (t)\rangle}{\delta t} </math>

This can be easily determined from the linear plots to produce the values seen below.

{| class="wikitable"
!
! colspan="2" |Diffusion Coefficient <math>m^2 s^{-1}</math>
|-
|
|8,000 particle simulation
|1,000,000 particle simulation
|-
|Vapour
|1.74
|3.12
|-
|Liquid
|9.26*10-2
|7.9810-2
|-
|Solid
|<nowiki>-2.31*10-7</nowiki>
|<nowiki>-8.75*10-7</nowiki>
|}
These values are in agreement with both the MSD against time plots and with theory.

'''VACF's'''
The velocity autocorrelation function can also be used to calculate the diffusion coefficient:
<math> C(T)=\langle v(t)*V(t+T)\rangle </math>

'''<span style="color:red">JC: Plots look good, include the line of best fit that you used to calculate D on your plots.  Why does the gas phase MSD begin as a curved line (ballistic regime)?  Did you calculate D from the VACF?</span>'''