Rep:Mod:CEW complab 2

Liquid Simulations

Introduction to molecular dynamics simulation

Numerical Integration

The Verlet algorithm and the modified velocity-Verlet algorithm can be used to numerically calculate the positions of atoms in a molecular dynamics simulation. These numerical methods require the simulation to be discretised into a series of timesteps, rather than treating the atomic positions, velocities and forces as continuous functions of time. The velocity-Verlet algorithm is:

${\mathbf{v}}_i(t+\delta t)={\mathbf{v}}_i(t+\frac{1}{2}\delta t)+\frac{1}{2}{\mathbf{a}}_i(t+\delta t)\delta t$

${\mathbf{v}}_i=\text{velocity of atom }i$

$\delta t=\text{timestep}$

${\mathbf{a}}_i=\text{accelaration of atom }i$

The plot below in figure 1 shows the atomic positions as a function of time as calculated by the velocity-Verlet algorithm, and the classical harmonic oscillator, where:

$$\begin{gathered} x(t)=Acos(\omega t+\varphi ), \\ A=1.00, \\ \omega =1.00, \\ \varphi =0.00. \end{gathered}$$

Figure 2 plots the energy as a function of time, which was calculated by summing the kinetic energy term, $\frac{1}{2}m{v}^2$, and the potential energy term, $\frac{1}{2}k{x}^2$, and figure 3 plots the error, which was calculated as the difference in the positions found by the velocity-Verlet algorithm and the classical harmonic oscillator, as a function of time. Figure 4 plots the error maxima from figure 3 as a function of time.

Figure 1: Plot of time vs position for the positions given by the velocity-Verlet algorithm "x(t)", and by the classical harmonic oscillator "ANALYTICAL".	Figure 2: Plot of the time vs total energy (kinetic and potential energy).
Figure 3: Plot of time vs error (difference in positions).	Figure 4: Plot of time vs error for the error maxima from Figure 3.

The choice of timestep can influence the error of the calculation, as a small timestep is desired to most accurately simulate the system but calculations with a smaller timestep take longer to run than those with a larger timestep. By the harmonic oscillator the total energy should be a constant over the course of the simulation, and it was found that a timestep of $0.21$ is required to ensure the total energy does not change by more than 1% over the course of the simulation. This can be determined by varying the timestep and calculating the size of the fluctuations of the total energy for the simulation, compared to the average constant energy value that would arise from the harmonic oscillator, so monitoring the total energy of of the system when modelling it numerically is important as it allows for the error of the calculation to be determined.

Atomic Forces

The Lennard-Jones potential describes molecular interactions, and is made up of a repulsive and an attractive part. A Lennard-Jones potential is shown in figure 5 and the equation that governs it is given below:

Figure 5: Lennard-Jones Potential

$\varphi (r)=4ϵ[\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6}]$

$\varphi =\text{intermolecular potential}$

$ϵ=\text{well depth}$

$\sigma =\text{Van der Waals radius}$

$r=\text{separation distance}$

Setting this to zero enables the separation at zero potential, $r_o$, to be found:

$\varphi (r)=4ϵ[\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6}]=0$

$\frac{{\sigma }^{12}}{r^{12}}=\frac{{\sigma }^6}{r^6}$

$r^6={\sigma }^6$

$r_o=\sigma$

The force is the derivative of the potential with respect to the separation and is shown for the Lennard-Jones potential below:

$F=-\frac{d\varphi (r)}{dr}=\frac{48ϵ{\sigma }^{12}}{r^{13}}-\frac{24ϵ{\sigma }^6}{r^7}$

When $r=r_o$ the force is given by:

$F=\frac{24ϵ}{\sigma }$

The equilibrium separation,$r_{eq}$, occurs when the force is zero so is found by:

$\frac{d\varphi (r)}{dr}=\frac{48ϵ{\sigma }^{12}}{r^{13}}-\frac{24ϵ{\sigma }^6}{r^7}=0$

$\frac{48ϵ{\sigma }^{12}}{r^{13}}=\frac{24ϵ{\sigma }^6}{r^7}$

$2{\sigma }^6=r^6$

$r_{eq}=\sigma \sqrt[6]{2}=1.12\sigma$

At $r=r_{eq}$ the depth of the potential well is given by:

$\varphi (r)=4ϵ[\frac{{\sigma }^{12}}{4{\sigma }^{12}}-\frac{{\sigma }^6}{2{\sigma }^6}]=4ϵ\times -\frac{1}{4}=-ϵ$

Taking $\sigma =ϵ=1.0$, the integral below can be expressed as:

$\int \varphi \left(r\right)\mathrm{d}r=\frac{4}{5r^5}-\frac{4}{11r^{11}}+C$

This result can be used to evaluate the integrals below:

$\underset{2\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=-0.0248$

$\underset{2.5\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=-0.0082$

$\underset{3\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=-0.0033$

Periodic Boundary Conditions

For simulations, realistic volumes of particles cannot be used as this leads to a huge number of atoms that need to be simulated. This can be shown by considering a system of water molecules:

Taking the concentration of water as $55.5{\text{ moldm}}^{-3}=0.0555{\text{ molml}}^{-1}$, under standard conditions, the number of molecules of water in $1\text{ ml}$ is the concentration of water multiplied by Avogadro's number ($6.02\times 10^{23}$):

$0.0555\times N_A=3.34\times 10^{22}\text{ molecules}$

The volume of $10000$ water molecules under standard conditions can be found by dividing the number of water molecules by Avogadro's number to convert to the number of moles of water, and by the concentration of water:

$\frac{10000}{0.0555N_A}=2.99\times 10^{-19}\text{ ml}$

For the simulations run it would not be possible to simulate $1\text{ ml}$ of water due to the large number of particles, however, applying periodic boundary conditions allows for bulk systems to be simulated with a small system volume. Applying periodic boundary conditions ensures that the number of particles is kept constant, and an example of applying these conditions is described below:

After an atom at position $(0.5,0.5,0.5)$ in a cubic simulation box which runs from $(0,0,0)$ to $(1,1,1)$ has been moved along the vector $(0.7,0.6,0.2)$, it will end up in the position $(0.2,0.1,0.7)$, due to the application of periodic boundary conditions, not outside the simulation box.

Reduced Units

The simulations run are carried out in reduced units. The example for argon below demonstrates how reduced units can be converted into real units:

The Lennard-Jones parameters for argon are $$\begin{gathered} \sigma =0.34\text{ nm},ϵ \\ / \\ {k}_B=120\text{ K} \end{gathered}$$, and the cutoff separation is $r^{*}=3.2$. These values are given in reduced units and can be converted into real units by:

$r=r^{*}\sigma =1.088\text{ nm}$

The well depth is given by $ϵ$, so can be found as:

$ϵ=120k_B=1.656\times 10^{-21}\text{ J}=0.99{\text{ kJmol}}^{-1}$

The reduced temperature is $T^{*}=1.5$, and can be converted into real units by:

$T=\frac{T^{*}ϵ}{k_B}=180\text{ K}$

Equilibration

Creating the simulation box

In these simulations, when particles are too close together they will have a high, repulsive force. Randomly generating the starting coordinates can lead to some atoms being very close to each other, which results in very large repulsive forces between them, and this can cause the calculation to fail due to the size of the force. Instead simulations start from a lattice, which will equilibrate over time. For a simple cubic lattice unit cell with lattice spacing $1.0772$, the number density of lattice points is found by:

$\frac{N}{V}=\frac{1}{1.0772^3}=0.800$

For a face centred cubic (FCC) lattice unit cell with the number density of lattice points $1.2$, the lattice spacing can be found using:

$\text{lattice spacing}=\sqrt[3]{\frac{\text{number of lattice points}}{\text{number density of lattice points}}}=\sqrt[3]{\frac{4}{1.2}}=1.4938$

A simulation for the simple cubic lattice with the input file command below leads to the formation of $1000$ atoms, as there is one atom per unit cell:

create_atoms 1 box

and this is acknowledged in the ouput file by the line:

Created 1000 atoms

For an FCC lattice the input command would lead to the formation of $4000$ atoms, as there are four atoms per unit cell in the FCC lattice.

Setting the properties of the atoms

The properties of the atoms in the simulation are defined by the lines below:

mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0

The first line of the script means the mass of particle $1$ is $1.0$, the second line means the global cutoff for the Lennard-Jones interactions is at a distance of $3.0$, and the third line means the pairwise force field coefficients for all atoms, from atoms $1$ to $N$, are $1.0$ and $1.0$. For these simulations the velocity-verlet algorithm is being used, as ${\mathbf{x}}_i\left(0\right)$ and ${\mathbf{v}}_i\left(0\right)$ have been specified.

Running the simulation

The lines from an input file below:

 ### SPECIFY TIMESTEP ###
 variable timestep equal 0.001
 variable n_steps equal floor(100/${timestep})
 variable n_steps equal floor(100/0.001)
 timestep ${timestep}
 timestep 0.001
 
 ### RUN SIMULATION ###
 run ${n_steps}
 run 100000

could be replaced by:

 timestep 0.001
 run 100000

The advantages of the first method are that a variable "timestep" is defined, so every time:

${timestep}

is used in the input file, the amount defined by the line:

variable timestep equal 0.001

is used. This means the simulation will run for the same amount of time, irrespective of the timestep used as the variable "n_steps" is defined as:

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

and this value is then used to determine the number of timesteps the simulation is run for in the line:

run ${n_steps}

Using the second method would require the number of timesteps needed to a run a simulation of a certain length to be calculated manually for each timestep used, which would take longer and could lead to errors.

Checking equilibration

It is important to check that the system reaches equilibrium over the course of the simulation. For the experiment with the timestep $0.001$ the simulation does reach equilibrium, at time $0.5$, as can be seen in Figures 6, 7, and 8. Figure 9 shows a plot of the energy of all five of the experiments, which were each run with a different timestep. It can be seen that the experiment run with timestep $0.015$ gave a very poor result, as the energy does not reach equilibrium. The largest timestep used to give a useful result is $0.01$ as it reaches equilibrium. However, for timesteps above $0.0025$ the energy is dependent on the timestep chosen, which is seen by the energies averaging at increasingly higher values for timesteps $0.0075$ and $0.01$, so the timestep $0.0025$ has been chosen to carry out further calculations.

Figure 6: Plot of time vs energy.	Figure 7: Plot of time vs temperature.
Figure 8: Plot of time vs pressure.	Figure 9: Plot of time vs energy for all of the timesteps.

Running simulations under specific conditions

Thermostats and Barostats

$\gamma$ is a constant factor that is required to keep the instantaneous temperature, $T$, and the target temperature, $\mathfrak{𝔗}$, equal. This is required to ensure the kinetic energy of the system remains at the correct value. It can be found using equipartition theory, where each degree of freedom contributes $\frac{1}{2}{k}_{B}T$, on average, to the energy. This gives equations one and two, which are divided by each other to give $\gamma$ in terms of $T$ and $\mathfrak{𝔗}$.

Equation one: $\frac{1}{2}\sum _i{m}_i{\left(\gamma v_i\right)}^2=\frac{{\gamma }^2}{2}\sum _i{m}_i{v}_i^2=\frac{3}{2}N{k}_B\mathfrak{𝔗}$

Equation two: $\frac{1}{2}\sum _i{m}_i{v}_i^2=\frac{3}{2}N{k}_{B}T$

${\gamma }^2=\frac{\mathfrak{𝔗}}{T}$

$\gamma =\sqrt{\frac{\mathfrak{𝔗}}{T}}$

Examining the Input Script

The input script below describes how average values will be determined.

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp press density
variable dens equal density
variable dens2 equal density*density
variable temp equal temp
variable temp2 equal temp*temp
variable press equal press
variable press2 equal press*press
fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_dens2 v_temp2 v_press2
run 100000

IN the penultimate line,

$100$

is the how often input values will be taken,

$1000$

is the number of times to use input values for calculating averages, and

$100000$

is how often averages are calculated. In this case averages will be calculated every

$100000$

timesteps, using

$1000$

measurements from the simulation, which are found by sampling the values every

$100$

timesteps before the average is calculated. The final line is the number of timesteps that the simulation will run for, so in this case

$100000$

timesteps of

$0.0025$

will be carried out, so the simulation will run for time

$250$

.

Temperature and Pressure Control

Simulations using the velocity-Verlet algorithm on the Lennard-Jones system were carried out at pressures $2.5$ and $3.0$, and temperatures $0.5$, $1.0$, $1.5$, $2.0$ and $2.5$ (values in reduced units), with timestep $0.0025$. The pressures and temperatures were chosen as they are close to the equilibrium values that were previously calculated, and the timestep was chosen at it was the largest that gave valid results. The plots in figures 10 and 11show both the computed values for the density using the velocity-Verlet algorithm and the predicted values, found using the perfect gas law with $k_B=1$ as the simulations are run in reduced units:

$\frac{N}{V}=\frac{P}{k_{B}T}$

Figure 10: Plot of density versus temperature for $2.5$ pressure.	Figure 11: Plot of density vs temperature for $3.0$ pressure.

The perfect gas law assumes that the volume of the particles is negligible and that there are no intermolecular interactions between the particles, so is best applied to dilute gas systems. The difference between the computed and predicted values increases with pressure because the system becomes less dilute, so less ideal. The computed values are higher than the predicted values as they were found considering intermolecular interactions, as is instructed in the script by the lines below (purpose of commands discussed previously):

pair_style lj/cut/opt 3.0
pair_coeff 1 1 1.0 1.0

Calculating heat capacities using statistical physics

The heat capacity of a system is the amount of energy needed to increase the temperature of the system by $1\text{K}$, so is a measure of the amount of thermal energy that can be absorbed. Generally this increases with temperature, as more degrees of freedom are possible (rotational and electronic, in addition to translational) so the system can absorb more thermal energy, but for these simulations the particles are taken as hard spheres so no rotations are possible, and since the simulations are classical no electronic transitions are considered. In the canonical ensemble (NVT) the heat capacity can be calculated using:

$C_V=\frac{\partial E}{\partial T}=\frac{\mathrm{V}\mathrm{a}\mathrm{r}\left[E\right]}{k_B{T}^2}=N^2\frac{⟨E^2⟩-{⟨E⟩}^2}{k_B{T}^2}$

The heat capacity was found using this equation for simulations of a Lennard-Jones system, with densities $0.2$ and $0.8$, at temperatures of $2.0$, $2.2$, $2.4$, $2.6$ and $2.8$ (all values in reduced units), with timestep $0.0025$. Figure 12 shows plots of heat capacity over volume vs temperature for each of the densities.

Figure 12: Plot of heat capacity over volume, vs temperature for a Lennard-Jones system at densities $0.2$ and $0.8$.

The plot in Figure 12 doesn't follow the expected increasing heat capacity with temperature, but instead the heat capacity decreases with temperature. This can be explained by considering that, at higher energies, the energy levels are closer together so for a given energy level there is a higher degeneracy. This means that in order to achieve a specific population of energy levels at a higher temperature, less energy is required than would be needed for the equivalent density of states at a lower temperature. Also, the heat capacity of the system with density $0.2$ is lower than that of the system with density $0.8$. This is due to there being more particles per unit volume at the higher density, so to increase the temperature by $1\text{K}$ there are more particles to absorb the energy before the temperature of the system is raised, at the higher density.

The input file for this simulation can be seen here: File:Cew 41.in.

Structural properties and the radial distribution function

The solid, liquid and vapour phases of a Lennard-Jones system were simulated using the densities and temperatures given below (in reduced units) :

Phase	Density	Temperature
Solid	1.20	1.40
Liquid	0.80	1.20
Vapour	0.01	1.11

The plots of the radial distribution function (RDF) and its integral from these simulations are shown in figures 13 and 14.

Figure 13: Plots of the RDF for the solid, liquid and vapour phases of the Lennard-Jones system.	Figure 14: Plots of the integral of the RDF for the solid, liquid and vapour phases of the Lennard-Jones system.

The peaks in the RDFs (figure 13) correspond to the nearest neighbours, so the RDF for the solid phase Lennard-Jones system has many clear peaks. However, those for the liquid and vapour phases do not due to the absence of long range order so the peaks become too small to be observed as the distance between nearest neighbour is too long. For the solid phase, the first three peaks in the RDF correspond to the first three nearest neighbours, which are illustrated in figure 15. The coordination numbers for these peaks can be found by comparing the peak positions in the RDF and the integration of the RDF (figure 14) at the at these positions. This analysis gives the coordination numbers $5.5$, $8.1$ and $18.6$ for the first, second and third peaks respectively. The lattice spacing, $a$, can be determined using trigonometry from the first nearest neighbour separation, $2R=1.025$ (determined from figure 13):

$a=4Rcos(45)=2R\sqrt{2}=1.450\text{ (3 d.p.)}$

Alternatively the lattice spacing can be taken as the distance to the second nearest neighbour, which results in a lattice spacing of $1.425$. This is good agreement with the calculated result above.

Figure 15: FCC lattice unit cell showing the three nearest neighbours (N.B.: not all atoms in unit cell shown)

Dynamical properties and the diffusion coefficient

Mean Squared Displacement

The solid, liquid and vapour phases of a Lennard-Jones system were simulated using the densities and temperatures used previously given, and timestep $0.002$. From these simulations the mean squared displacement (MSD) was calculated. Figures 16 to 21 below show plots of the MSD vs the timestep for a Lennard-Jones solid, liquid and gas system, with $8000$ and $1000000$ atoms. The gradient of the line increases on moving from the solid to the liquid to the vapour phase, which was expected, as the atoms are able to move most easily in the vapour phase, so will have a greater MSD.

Figure 16: Plot of MSD vs timestep for a Lennard-Jones solid, with $8000$ atoms.	Figure 17: Plot of MSD vs timestep for a Lennard-Jones liquid, with $8000$ atoms.	Figure 18: Plot of MSD vs timestep for a Lennard-Jones vapour, with $8000$ atoms.
Figure 19: Plot of MSD vs timestep for a Lennard-Jones solid, with $1000000$ atoms.	Figure 20: Plot of MSD vs timestep for a Lennard-Jones liquid, with $1000000$ atoms.	Figure 21: Plot of MSD vs timestep for a Lennard-Jones vapour, with $1000000$ atoms.

The diffusion coefficient can be found from the mean squared displacement by the equation below:

$D=\frac{1}{6}\frac{\partial ⟨r^2\left(t\right)⟩}{\partial t}$

The gradient of the line, once it has established linear behaviour, can be taken and converted to a function of time (instead of timestep) by dividing the gradient by the timestep, $0.002$. This can then be divided by $6$ to give the diffusion coefficient. The results are summarised below:

Type of System	with $8000$ atoms	with $1000000$ atoms
Lennard-Jones Solid	Gradient$\approx 0$, $D\approx 0$	Gradient$\approx 0$, $D\approx 0$
Lennard-Jones Liquid	Gradient$=0.001$, $D=\frac{0.50}{6}=0.083$	Gradient$=0.001$, $D=\frac{0.50}{6}=0.083$
Lennard-Jones Vapour	Gradient$=0.080$, $D=\frac{40}{6}=6.667$	Gradient$=0.016$, $D=\frac{8}{6}=1.333$

Velocity Autocorrelation Function

The velocity autocorrelation function (VACF), given by $C\left(\tau \right)$, is another method that can be used to calculate the diffusion coefficient, as:

$D=\frac{1}{3}\underset{0}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}\tau ⟨\mathbf{v}\left(0\right)\cdot \mathbf{v}\left(\tau \right)⟩$

The VACF can be found by evaluating :

$C\left(\tau \right)=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}v\left(t\right)v(t+\tau )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{v}^2\left(t\right)\mathrm{d}t}$

For the 1D harmonic oscillator:

$v\left(t\right)=\frac{\mathrm{d}x(t)}{\mathrm{d}t}$ and $x(t)=Acos(\omega t+\varphi )$

The VACF for the 1D harmonic oscillator can be evaluated to give the result shown below:

$C\left(\tau \right)=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}v\left(t\right)v(t+\tau )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{v}^2\left(t\right)\mathrm{d}t}$

$$\begin{gathered} C\left(\tau \right)=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}[-\omega A \\ sin(\omega t+\varphi )][-\omega A \\ sin(\omega (t+\tau )+\varphi )]\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}[-\omega A \\ sin(\omega t+\varphi )][-\omega A \\ sin(\omega t+\varphi )]\mathrm{d}t} \end{gathered}$$

$$\begin{gathered} C\left(\tau \right)=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}[-\omega A \\ sin(a)][-\omega A \\ sin(a+b)]\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}[-\omega A \\ sin(a)][-\omega A \\ sin(a)]\mathrm{d}t} \end{gathered}$$ ($a=\omega t+\varphi$, $b=t+\tau$)

$$\begin{gathered} C\left(\tau \right)=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{A}^2{\omega }^2 \\ si{n}^{2}a \\ cosb \\ +sina \\ cosa \\ sinb \\ \mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{A}^2{\omega }^2 \\ si{n}^{2}a \\ \mathrm{d}t} \end{gathered}$$

$$\begin{gathered} C\left(\tau \right)=\frac{A^2{\omega }^2 \\ cosb \\ \underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^{2}a \\ \mathrm{d}t}{A^2{\omega }^2 \\ \underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^{2}a \\ \mathrm{d}t}+\frac{A^2{\omega }^2 \\ sinb \\ \underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}sina \\ cosa \\ \mathrm{d}t}{A^2{\omega }^2 \\ \underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^{2}a \\ \mathrm{d}t} \end{gathered}$$

$C\left(\tau \right)=cos(\omega \tau )$

Figure 22 shows the VACF for a Lennard-Jones solid and liquid with $8000$atoms, which both show fluctuations due to changes in velocity of the particles. These are caused by collisions with other particles in the system, which cause a change in the direction of the motion of the particle, hence the change in velocity. The differences between the fluctuations observed in the solid and liquid VACFs is due to the distances between the particles, so in the solid the particles are closer together so collide more frequently than in the liquid, which leads to more fluctuations in the VACF for the solid. Furthermore, for both the solid and liquid the VACF decays to zero, as the energy of the particles is dispersed randomly throughout the system upon collisions between particles. The differences between the harmonic oscillator VACF ("analytical") and the Lennard-Jones solid and liquid system are that there are regular fluctuations in the harmonic oscillator, and that the system doesn't decay to zero. The regular fluctuations are caused by changes of velocity each time the spring reaches its fully extended state, as is governed by Hooke's law:

$F=-kx$

The system doesn't decay to zero because there are no collisions in the harmonic oscillator, so the energy of the particles remains constant and isn't randomly dispersed among the particles.

Figure 22: Plot of the velocity autocorrelation function vs timestep for a Lennard-Jones solid and liquid, and for the harmonic oscillator ("analytical").

The integral under the VACF can be estimated using the trapezium rule, and this can be used to estimate the diffusion coefficient, as described above. Figures 23, 24 and 25 show the running integrals for each of the Lennard-Jones solid, liquid and vapour phases with $8000$ atoms and figures 26, 27 and 28 show the running integrals for them with $1000000$ atoms. The running integrals for the solid systems show that the VACF reaches equilibrium, where the gradient decreases to close to zero. This is also true for the liquid simulation with $1000000$ atoms, but not for the other simulations of the liquid and vapour phases. The solid reaches equilibrium the most rapidly as the atoms are able to move the least, but this occurs most slowly in the vapour systems as the particles have more energy so are able to move around more rapidly. This means it takes a longer amount of time for the velocities to reach an average, equilibrium value.

Figure 23: Running integral for the VACF for the Lennard-Jones solid, with $8000$ atoms.	Figure 24: Running integral for the VACF for the Lennard-Jones liquid, with $8000$ atoms.	Figure 25: Running integral for the VACF for the Lennard-Jones vapour, with $8000$ atoms.
Figure 26: Running integral for the VACF for the Lennard-Jones solid, with $1000000$ atoms.	Figure 27: Running integral for the VACF for the Lennard-Jones liquid, with $1000000$ atoms.	Figure 28: Running integral for the VACF for the Lennard-Jones vapour, with $1000000$ atoms.

For the Lennard-Jones solid, liquid and vapour, with $8000$ and $1000000$ atoms, the diffusion coefficients were predicted by the method described above to give the results in the table below. The largest source of error in the estimates of the diffusion coefficient from the VACF is that it is impossible to calculate the integral for infinite time, so this introduces error into calculating the diffusion coefficient, especially when the system doesn't reach an equilibrium state.

Type of System	with $8000$ atoms	with $1000000$ atoms
Lennard-Jones Solid	$D=\frac{1}{3}(8.758\times 10^{-5})=2.919\times 10^{-5}\approx 0$	$D=\frac{1}{3}(2.742\times 10^{-4})=9.124\times 10^{-5}\approx 0$
Lennard-Jones Liquid	$D=\frac{1}{3}(0.587)=0.200$	$D=\frac{1}{3}(0.270)=0.090$
Lennard-Jones Vapour	$D=\frac{1}{3}(47.1)=15.7$	$D=\frac{1}{3}(19.6)=6.54$

Summary

Molecular dynamics simulations can be used to determine a lot of thermodynamic data about a system, and good agreement of simulated data and classically calculated data is seen. However, it is important to note that error can arise due to simulations being run with a small number of atoms, and over a limited period of time.