Talk:Mod:IOPJKL

Liquid Simulations Computational Experiment

Theory

Task 1: Calculations of a Harmonic Oscillator using both Classical and Velocity-Verlet solutions

Figure 1.a) Positions calculated via the two solutions vs time b) Difference between the positions calculated via different methods against time c) Total energy of the harmonic oscillator calculated using Velocity-Verlet solutions vs time

Position calculations from the two methods: The position $x$ at time $t$ was determined using $x(t)=Acos(\omega t+\varphi )$ where; $A=1,\omega =1$ and $\varphi =0$. These were plotted against time(Figure 1.a).

Error analysis of the two positions: The difference between the two positions was calculated and plotted against time(Figure 1.b). The error increases with time and position, this loss of precision could come from two sources. As the position gets further from the original position the $a(t)d{t}^2$ part of the equation has less importance also the errors could build up since the previous position is used in the preceding position calculation.

NJ: Yes and no. Your point about the acceleration being less significant isn't true (that depends on how large the forces are), but the other point is spot on. Well done.

Total Energy of the Harmonic-Oscillator vs time: Values for velocity $v$ and position $x$ were obtained from the Velocity-Verlet calculation, values for the force-constant $k$ and mass $m$ were both given as 1. Using $E=\frac{1}{2}m{v}^2+\frac{1}{2}k{x}^2$ the total energies were calculated and plotted against time(Figure 1.c).

Task 2: Plot of maxima error positions vs time

The coordinates of the 5 maximas in figure 2 were estimated and plotted on a new graph. This yielded a linear function where

$x(t)=0.0004x-7*10^{-5}$

.

Figure 2. Linear plot of the maxima error positions(figure 1.b) against time

Task 3: Timesteps to minimise variation in total energy

Multiple timesteps were used to calculate the total energy of the harmonic oscillator. The percentage difference was taken between the initial energy assumed to be correct and the energy of the first minima, these results were then plotted(figure 3). In order to keep the change in total energy below 1% a timestep equal to or below 0.2 should be used. The change in total energy of the simulation should be kept to a minimum as the total energy should not vary over time in the harmonic oscillator.

NJ: Can you say a bit more about this? Why shouldn't it vary?

Figure 3. Plot of percentage energy difference and timestep

Task 4: Lennard-Jones interaction

1. Calculation of distance at zero-potential, r = r₀: $\varphi \left(r_0\right)=4ϵ(\frac{{\sigma }^{12}}{r_0^{12}}-\frac{{\sigma }^6}{r_0^6})$

where φ = potential energy r = the distance ε = well depth and σ = distance from zero potential distance
                                $\varphi \left(r_0\right)=4ϵ({\sigma }^{12}-{\sigma }^6{r}_0^6)$

$=4ϵ(1-\frac{r_0^6}{{\sigma }^6})=0$

$=1-\frac{r_0^6}{{\sigma }^6}=0$

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

$\sigma =r_0$

2. Force at r₀:

${\mathbf{F}}_0=-\frac{\mathrm{d}\varphi \left({\mathbf{r}}^N\right)}{\mathrm{d}{\mathbf{r}}_0}=-\frac{\mathrm{d}}{\mathrm{d}{\mathbf{r}}_0}4ϵ(\frac{{\sigma }^{12}}{r_0^{12}}-\frac{{\sigma }^6}{r_0^6})=4ϵ(-12\frac{{\sigma }^{12}}{r_0^{13}}+6\frac{{\sigma }^6}{r_0^7})$

Since $\sigma =r_0$

$=4ϵ(-12\frac{{\sigma }^{12}}{{\sigma }^{13}}+6\frac{{\sigma }^6}{{\sigma }^7})=4ϵ\frac{-6}{\sigma }=-\frac{24ϵ}{\sigma }$

NJ: You've dropped a negative here.

3. Finding the equilibrium separation r_eq:

${\mathbf{F}}_{eq}=-\frac{\mathrm{d}\varphi \left({\mathbf{r}}^N\right)}{\mathrm{d}{\mathbf{r}}_{eq}}=-\frac{\mathrm{d}}{\mathrm{d}{\mathbf{r}}_{eq}}4ϵ(\frac{{\sigma }^{12}}{r_{eq}^{12}}-\frac{{\sigma }^6}{r_{eq}^6})=4ϵ(-12\frac{{\sigma }^{12}}{r_{eq}^{13}}+6\frac{{\sigma }^6}{r_{eq}^7})=-48ϵ\left(\frac{{\sigma }^{12}}{r_{eq}^{13}}\right)+24ϵ\left(\frac{{\sigma }^6}{r_{eq}^7}\right)=0$

$48ϵ\left(\frac{{\sigma }^{12}}{r_{eq}^{13}}\right)=24ϵ\left(\frac{{\sigma }^6}{r_{eq}^7}\right)$

$2ϵ\left(\frac{{\sigma }^{12}}{r_{eq}^{13}}\right)=ϵ\left(\frac{{\sigma }^6}{r_{eq}^7}\right)$

$2ϵ\left({\sigma }^{12}\right)=ϵ\left({\sigma }^6{r}_{eq}^6\right)$

$$\begin{gathered} 2ϵ\left({\sigma }^6\right)=ϵ\left( \\ {r}_{eq}^6\right) \end{gathered}$$

$$\begin{gathered} 2\left({\sigma }^6\right)=\left( \\ {r}_{eq}^6\right) \end{gathered}$$

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

4. Finding the well-depth at equilibrium:

$\varphi \left(r_{eq}\right)=4ϵ(\frac{{\sigma }^{12}}{r_{eq}^{12}}-\frac{{\sigma }^6}{r_{eq}^6})=4ϵ(\frac{{\sigma }^{12}}{(\sqrt[6]{2}\sigma {)}^{12}}-\frac{{\sigma }^6}{(\sqrt[6]{2}\sigma {)}^6})=4ϵ(\frac{{\sigma }^{12}}{2^2{\sigma }^{12}}-\frac{{\sigma }^6}{2{\sigma }^6})=4ϵ(\frac{1}{4}-\frac{1}{2})=4ϵ\left(\frac{-1}{4}\right)=-ϵ$

5. Evaluation of integrals:

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

$\underset{2}{\overset{\mathrm{\infty }}{\int }}4ϵ(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})\mathrm{d}r=\underset{2}{\overset{\mathrm{\infty }}{\int }}4ϵ(\frac{1^{12}}{r^{12}}-\frac{1}{r^6})\mathrm{d}r=\underset{2}{\overset{\mathrm{\infty }}{\int }}4ϵ(\frac{1}{r^{12}}-\frac{1}{r^6})\mathrm{d}r$

$=[4ϵ(\frac{1}{-11r^{11}}-\frac{1}{-5r^5}){]}_2^{\mathrm{\infty }}$

$=4ϵ(\frac{1}{-11(\mathrm{\infty }{)}^{11}}-\frac{1}{-5(\mathrm{\infty }{)}^5})-4ϵ(\frac{1}{-11(2{)}^{11}}-\frac{1}{-5(2{)}^5})$

$=-0+4ϵ(\frac{1}{-11(2{)}^{11}}-\frac{1}{-5(2{)}^5})=-24.82KJ/mol$

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

$\underset{2.5}{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=\underset{2.5}{\overset{\mathrm{\infty }}{\int }}4ϵ(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})\mathrm{d}r=\underset{2.5}{\overset{\mathrm{\infty }}{\int }}4ϵ(\frac{1^{12}}{r^{12}}-\frac{1}{r^6})\mathrm{d}r=\underset{2.5}{\overset{\mathrm{\infty }}{\int }}4ϵ(\frac{1}{r^{12}}-\frac{1}{r^6})\mathrm{d}r$

$=[4ϵ(\frac{1}{-11r^{11}}-\frac{1}{-5r^5}){]}_{2.5}^{\mathrm{\infty }}$

$=4ϵ(\frac{1}{-11(\mathrm{\infty }{)}^{11}}-\frac{1}{-5(\mathrm{\infty }{)}^5})-4ϵ(\frac{1}{-11(2.5{)}^{11}}-\frac{1}{-5(2.5{)}^5})$

$=-0+4ϵ(\frac{1}{-11(2.5{)}^{11}}-\frac{1}{-5(2.5{)}^5})=-8.18KJ/mol$

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

$\underset{3}{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=\underset{3}{\overset{\mathrm{\infty }}{\int }}4ϵ(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})\mathrm{d}r=\underset{3}{\overset{\mathrm{\infty }}{\int }}4ϵ(\frac{1^{12}}{r^{12}}-\frac{1}{r^6})\mathrm{d}r=\underset{3}{\overset{\mathrm{\infty }}{\int }}4ϵ(\frac{1}{r^{12}}-\frac{1}{r^6})\mathrm{d}r$

$[4ϵ(\frac{1}{-11r^{11}}-\frac{1}{-5r^5}){]}_3^{\mathrm{\infty }}$

$4ϵ(\frac{1}{-11(\mathrm{\infty }{)}^{11}}-\frac{1}{-5(\mathrm{\infty }{)}^5})-4ϵ(\frac{1}{-11(3{)}^{11}}-\frac{1}{-5(3{)}^5})$

$=-0+4ϵ(\frac{1}{-11(3{)}^{11}}-\frac{1}{-5(3{)}^5})=-3.29KJ/mol$

NJ: Be careful with units, these should be reduced.

Task 5: Estimating the number of molecules in volumes of H₂O

1. Number of molecules in 1mL:

$Mr(H_{2}O)=18.02g/mol,\rho =1g/c{m}^3$

$n=\frac{m}{Mr}=\frac{1}{18.02}=5.55x10^{-2}$ $Numberofmolecules=5.55x10^{-2}x6.022x10^{23}=3.34192x10^{22}$

2. Volume of 1000 molecules:

$n=\frac{1000}{6.022x10^{23}}=1.66x10^{-21}mol$ $m=1.66x10^{-21}x18.02=2.99x10^{-20}g$ $V=\frac{2.99x10^{-20}}{1}=2.99x10^{-20}c{m}^3$ NJ: This is out by a factor of 10, should be ${𝒪}(10^{-19})$

Task 6: Periodic boundary conditions

When the atom at $(0.5,0.5,0.5)$ in the cubic simulation box $(0,0,0,)$ to $(1,1,1)$ moves along the vector $(0.7,0.6,0.2)$ it ends up at $(0.2,0.1,0.7)$ (when the periodic boundry condition is taken into consideration).

Task 7: Converting reduced values to real values:

1. Reduced distance

$r*=3.2=\frac{r}{\sigma }$

$r=3.2x0.34x10^{-9}=1.09x10^{-}9=1.09nm$

2. Well-depth

$\frac{ϵ}{K_B}=120$

$ϵ=120K_B=\frac{8.69x10^{24}}{6.022x10^{23}}=997.71Jmo{l}^{-1}=9.98x10^{-1}KJmo{l}^{-}1$

3. Temperature

$T^{*}=1.5=\frac{K_{B}T}{ϵ}$

$T=\frac{T^{*}ϵ}{K_B}=\frac{1.5x1.657x10^{-21}}{K_B}=180K$

Equilibration

Task 8: Problems with random starting coordinates in non-solid simulations

If atoms were given random starting coordinates then they could end up within a contact proximity, they could bump into each other and then they would not only experience the potential of each but also the kinetic energy.

NJ: Not really. They would be closer than they "should be", as you suggest, and this leads to very large forces between atoms. Consequently, atoms move apart from each other very quickly, and a very short timestep would be needed to properly calculate the trajectory.

Task 9: The number density of lattice points calculations:

The number density n is given by the following equation

$n=\frac{N}{V}$ where N is the number of atoms/lattice points per unit cell and V is volumes of the unit cell.

1. The number density of a simple cubic lattice with inter lattice point distance of 1.07722

For a simple cubic lattice there is one atom/lattice point, therefore;

$n=\frac{N}{V}=\frac{N}{r^{*3}}=\frac{1}{1.07722^3}=0.8$

2. Side length of a face-centred cubic lattice with a number density of 1.2

A face-centred cubic lattice has 4 atoms/lattice points, therfore;

$r^{*}=\sqrt[3]{V}=\sqrt[3]{\frac{N}{n}}=\sqrt[3]{\frac{4}{1.2}}=1.49380$

Task 10: How many atoms are created by create_atoms 1 box command for a face-centred cubic lattice?

Since there are 4 atoms per lattice and the box is ten lattice spacings larger than the single lattice spacing then there would be 4000 atoms.

Task 11: LAMMPS commands for setting the properties of atoms:

1. Mass command

mass I value = mass 1 1.0  

This command is used for setting the mass of the atom, in the general command I is the atom type and value is the mass(in reduced units). In the example the atom is type 1 and the mass is 1.0.

2. Pair_style command

pair_style style args = pair_style lj/cut 3.0

The pair_style command is used to describe paired interactions, style is the interaction which is being looked at and args is the argument used for that style. In the example above the style is defining the cutoff distance, that is the distance at which interactions become negligible. The argument for the cutoff distance is 3.0, after which the simulation will not compute.

3. pair_coeff command

pair_coeff * * args = pair_coeff * * 1.0 1.0 This command is used to define the force field coefficients between pairs of atoms, the asterix terms are to define the coefficients for all atoms, the second part of this command, args, defines the coeffecients.

NJ: What are the coefficients in this case?

Task 12: Specifying the velocity of each atom

The Velocity-Verlet algorithm is used to specify velocity and position $$\begin{gathered} {\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 \\ (9) \end{gathered}$$

Task 13: Running the simulation:

The command $(100/timestep) command allows for the immediate evaluation of the the formula with the pre-defined variable i.e. the timestep. The following line variable n_steps equal floor(100/$({timestep}) will be broken down into individual contributions. The variable component is used to allows for the return to re-run a command, this allows creates a loop. The n_steps part describes the number of timesteps you wish to calculate which is defined in the 'Run Simulation' command. Equal is used to create a mathematical argument. Floor(11/$({timestep}) is useful in setting the maximum timestep bigger than the original one. In the run simulation part first we define the number of timesteps to be considered and then report the floor value i.e. the total time it takes. These lines are essential in creating a loop where each new timestep may be calculated subsequently. Make it easier to change the timestep, all values dependedent on the timestep will be changed with it and so it reduces possible human errors when changing all variables and also makes it easier and faster.

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 NJ: This section makes sure that you always simulate the same time (nsteps*dt), no matter what timestep is chosen.

Task 14: Checking equilibration

1. Pressure, Temperature and energy against time for simulation ran at 0.001 timestep

Figure 5. Temperature, pressure and total energy against time for the simulation using a timestep of 0.001

2. Equilibration

The simulation reaches equilibration at approximately 0.28, this is the timestep at which the simulation becomes relatively constant.

3. Maximising results accuracy whilst keeping the timestep within a reasonable timescale

The best timestep for maximum accuracy and sufficient information is 0.0025. It overlaps with 0.001 and so although it is more than twice the time as 0.001 it maintains a satisfactory level of accuracy. The others give equilibrium energies which are too far off the most accurate value, with the worst being for 0.01 as it diverges quite dramatically.

NJ: I think you mean 0.015

Figure 5. The energy vs time plots of simulations with varying timesteps

Running simulations under specific conditions

Task 15: Finding the expression for factor $\gamma$

$\gamma$ is a factor by which the velocity can be corrected.

$\frac{1}{2}\sum m_i{v}_i^2=\frac{3}{2}N{K}_{B}T$ (1) $\frac{1}{2}\sum m_i{v}_i^2=\frac{3}{2}N{K}_B{T}^{\prime }$ (2)

Where T' is the target temperature and T is the instantaneous temperature

$\frac{1}{2}\sum m_i{v}_i^2=\frac{3}{2}N{K}_{B}T$ $\frac{1}{2}{m}_i{v}_i^2=\frac{1}{2}N{K}_{B}T$ $m_i{v}_i^2=N{K}_{B}T$ $v_i^2=\frac{N{K}_{B}T}{m_i}$ $v_i=\sqrt[2]{\frac{N{K}_{B}T}{m_i}}$

Insert v_i into equation x

$\frac{1}{2}\sum m_i{\gamma }^2\frac{N{K}_{B}T}{m_i}=\frac{3}{2}N{K}_B{T}^{\prime }$ $\frac{1}{2}{m}_i{\gamma }^2\frac{N{K}_{B}T}{m_i}=\frac{1}{2}N{K}_B{T}^{\prime }$ ${\gamma }^{2}N{K}_{B}T=N{K}_B{T}^{\prime }$ ${\gamma }^2=\frac{N{K}_B{T}^{\prime }}{N{K}_{B}T}$ ${\gamma }^2=\frac{T^{\prime }}{T}$ $\gamma =\sqrt[2]{\frac{T^{\prime }}{T}}$

NJ: Good effort, but you can't take a single atomic velocity out of the sum like this.

Task 16: The importance of number, repetition and frequency of average calculations with timesteps

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

100 sets the number of timesteps after which another thermodynamic value is used to calculate the average, 1000 is the number of thermodynamic values which are used for the average and 1000000 sets the frequency at which averages are calculated.

100000 timesteps will be run.

Task 17: Equation of state plots

5 temperatures above the critical temeperature were chosen(1.6, 1.8, 2.0, 3.0 and 5.0) these were all simulated with a pressure of 2.6(n-1) and also 3.3(n-2). The timestep chose was 0.01 as in previous experiments it was the largest timesteps yielding the most accuracy.

NJ: What do you mean by (n-1) and (n-2)? 0.01 is too large a timestep, for the reasons you gave earlier!

The pressures were plotted against the temperatures along with their error bars.

The ideal gas densities were calculated with the following equation: $\rho =\frac{N}{V}=\frac{P}{K_{B}T}=\frac{P}{T}$ where K_B = 1 in reduced terms

Figure 6. Densities n-1 and n-2 calculated via simulation and Ideal Gas methods plotted against temperature

The simulated density is lower than that of the ideal gas method, this can be explained by an assumption to satisy an ideal gas. The atoms must not have any forces upon them and also they mustn't take up any space. In the simulation the atoms do take up space and therefore they have a lower density. This discrepancy increases with pressure. The general trend of the two plots is a decrease in heat capacity as temperature is increased.

NJ: Why does the discrepancy increase with pressure? What is the effect of changing the temperature? Discuss more about why your simulation isn't an ideal gas - why can't it be more dense than the ideal gas?

NJ: You shouldn't just join up the points in plots like this. You could fit a line or curve instead, if you wanted to.

Heat Capacity Calculation

Task 18: The temperature dependence of the isochoric heat capacity

Figure 7. Cv/V as a function of T* at P= 0.2 and 0.8

A LAAMPS script which could calculate the isochoric heat capacity from calculated thermodynamic values was used figure 7. The Cv/V was calculated with the following equation; $C_v=N^2\frac{<E^2>-<E{>}^2}{K_B{T}^2}$

Figure 8. Cv/V as a function of T* at P= 0.2 and 0.8

The effect on heat capacity at different temperatures and pressures

The decrease in heat capacity at higher temperatures is a result of the non-classical behaviour of the system at low temperatures. The energy levels behave as if they are discrete, the lower the temperature the more discrete they are. As illustrated here, the energy levels are closer together at higher temperatures e.g.$T_2$ and further apart at lower temperatures e.g.$T_1$. When a heat increment, $dq$, is introduced into the two systems we see that it has a different effect in both. At $T_2$, $dq$ will give a greater temperature increase than $T_1$, this causes the denominator in the equation below to become more important thus decreasing the isochoric heat capacity.

$C_V=\frac{dq}{dT}$

NJ: This is a good explanation, but you have to be careful with it. This is only an analogy - your simulated system is strictly classical, and so doesn't *actually* have discrete energy levels. Thinking about it in the way you've described, however, allows us to explain the trend.

The difference between the two pressured systems in figure 7 is what would be expected, the higher pressure plot has a higher isochoric heat capacity. At a higher pressure the atoms are closer together, they can thus transfer heat more effectively with a given temperature change than a system with a lower pressure, since the atoms are further apart.

NJ: This isn't about transferring heat though, it's about storing it. The answer is that heat capacity is extensive; the more atoms you have per volume, the more degrees of freedom there are in which to store heat.

Structural properties and the radial distribution function

Task 18: The Radial Distribution Functions of a solid, liquid and gas

1. Comparisons of the solid, liquid and gas radial distributions functions as well as the integrated function

figure 9

figure 10

The distribution plot(figure 9) of the gas has one short, broad peak which is close to the reference atom; as you move further away the distribution becomes constant. This characteristic shows the small probability of finding an atom close to the reference atom; also when it becomes constant this shows how there is negligible probability of finding another atom at the distance studied. When looking at the integrated function(figure 10), it emphasizes how few atoms you will find in the distance, this is characteristic of a gas' dispersed atoms.

NJ: No! There is only a negligible probability if g(r)=0

The liquid's distribution has a taller first peak than the gas' and then multiple peaks tapering off as you move further away; it is evident that there is a higher probability of finding an atom throughout the space than in a gas. Additionally the integrated distribution(figure 10) shows a larger number of atoms as you move through the space. This demonstrates the higher proximity of the atoms in the liquid state.

Finally the solid's distribution has a very sharp and tall peak close to the reference atom but then as you move further away the peaks become smaller but continuously fluctuate. The sharpness of the peaks and their continuous distribution arise from the highly ordered structure of the solid. The peaks are sharp unlike the gas and liquid because all of the atoms are in fixed positions and so at a certain position there will be a number of atoms in an equivalent position around the sphere of the distribution function. Furthermore the solid's first peak is considerably larger than the other phases, the probability is much greater due to the closeness of the atoms. The integrated function shows a very large number of atoms in the space, this again illustrates the closeness of the atoms in the solid state.

NJ: Your description of the solid is quite good. If the atoms are fixed, though, why aren't the peaks just vertical lines? Why do they still have a finite width? It's not surprising that the solid has more neighbours than the liquid, and then the gas - the densities decrease in that order. The liquid has much more structure than the gas, though.

2. Using the solid RDF to extract information on the lattice structure

figure 11

The three first peaks of the solid's radial distribution function are distances 1.0, 1.5 and 1.8. The lattice spacing is 1.5 and the calculated distances $r_a$,$r_b$ and $r_c$(where $r$ is the distance from the reference atom to the atoms a, b and c) were found to be 1.0, 1.5 and 1.8(respectively). Therefore the first three peaks 1, 2 and 3 correspond to the lattice points a, b and c(respectively).

NJ: How did you calculate the lattice distances? You should get three expressions for the distances, all of which involve a. You can then get three different estimates of a.

The coordination of atom a, b and c is 12, 6 and 24(respectively). These values were interpolated from the magnification of the integrated solid distribution function.

NJ: Nice plot. Is this consistent with an FCC lattice?

Dynamical properties and the diffusion coefficient

Task 19: The mean squared displacement(MSD) of large and small simulations

figure 12

$MSD(t)=<((t_0+t)-(t_0)){>}^2$ where; $t_0=$ the initial timestep and $t=$the new timestep.

The mean squared displacement(equation above) is a measure of how much random motion there is within the system. In this simulation the ability to move randomly is determined by the amount of free space and therefore the density. In the gaseous phase there is a lot of free space, the atoms move freely and encounter only a small number of atoms to collide with thus impeding it's motion; the linear increase of the plot reflects this behavior. The liquid has a lower degree of random motion than the gas, due closer proximity of the atoms, forces such as drag will now effect the atom's motion as it moves past other atoms. The two first plots show that the atom is able to move in drift, however the liquid experiences more drag and has a smaller gradient. The solid displays the least amount of translational freedom, the close packed structure prevents any random translational motion; it does however show some random movement, arising from vibrations of the constrained atoms. Therefore the plots in figure 12 are a reasonable representation of the random movement of the three phases.

NJ: Why is the vapour MSD curved?

Finding the diffusion coefficient $D$

$D=\frac{1}{6}\frac{\delta <r^2(t)>}{\delta t}>$ where $D=$ the diffusion coefficient in $m{s}^{-2}$,$<r^2(t)=$ MSD and $t=$ timestep The gradient of the slopes give the diffusion coefficient, however this uses a dimensionless timestep denominator, therefore we must convert it to time units by dividing by the timestep(0.002).

The gradient was found by fitting a linear trendline through the function.

Solid

Conversion of the gradient to give time units gave $D=4.17x10^{-6}m{s}^{-2}$ for the million atom simulation and $D=5.83x10^{-7}m{s}^{-2}$ for the 3200 atom simulation.

Liquid

Conversion of the gradient to give time units gave $D=8.33x10^{-2}m{s}^{-2}$ for the million atom simulation and $D=8.33x10^{-2}m{s}^{-2}$ for the 3200 atom simulation.

Gas

A trendline was fitted to the function, since the simulation took some time to produce the linear behavior it was fitted to be linear with the latter part of the function. Conversion of the gradient to give time units gave $D=3.09m{s}^{-2}$ for the million atom simulation and $D=2.43m{s}^{-2}$ for the 3200 atom simulation.

NJ: Be careful with units, these should be reduced. What do you think is the origin of the differences between your results and the million atom results?

Task 20: The velocity auto-correlation function(VACF)

Evaluating the 1D harmonic oscillator VACF using the expression for the position as a function of time

The normalized VACF is given by;

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

We must simplify this expression in order for it to be useful, to do this we must first find expressions for $v(t)$ and $v(t+\tau )$:

The expression for position as a function of time for the HO is defined as; $x(t)=A\cos (\omega t+\varphi )$

This can be substituted into the expression of velocity as a function of time;

$v(t)=\frac{dx(t)}{dt}=\frac{dA\cos (\omega t+\varphi )}{dt}$

This is then differentiated using the chain rule;

$\frac{dx(t)}{dt}=\frac{dx(t)}{du}x\frac{du}{dt}$

Where $u=\omega t+\varphi$

$v(t)=-\omega A\sin (\omega t+\varphi )$

We also need our velocity as a function of time after $\tau$, this is acheived by substituting in $t=t+\tau$;

$v(t+\tau )=-\omega A\sin (\omega (t+\tau )+\varphi )$ $=-\omega A\sin (\omega t+\omega \tau +\varphi )$

The two velocity expressions are substituted into the VACF;

$C(\tau )=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\omega A\sin (\omega t+\varphi )\omega A\sin (\omega t+\omega \tau +\varphi )dt}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\omega A\sin (\omega t+\varphi {)}^{2}dt}$

$=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi )\sin (\omega t+\omega \tau +\varphi )dt}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi {)}^{2}dt}$

Using the trigonometric identity $\sin (\alpha +\beta )=\sin (\alpha )\cos (\beta )+\cos \alpha \sin \beta$ we can re-write our equation;

where $\omega t+\varphi =\alpha$ and $\omega \tau =\beta$

$C(\tau )=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi )\sin (\omega t+\varphi )\cos (\omega \tau )+\cos (\omega t+\varphi )\sin (\omega \tau )dt}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi {)}^{2}dt}$

$C(\tau )=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi {)}^2\cos (\omega \tau )+\cos (\omega t+\varphi )\sin (\omega \tau )dt}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi {)}^{2}dt}$

The fraction can be expanded;

$C(\tau )=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi {)}^2\cos (\omega \tau )dt}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi {)}^{2}dt}+\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\cos (\omega t+\varphi )\sin (\omega \tau )dt}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi {)}^{2}dt}$

Then since $\cos (\omega \tau )$ does not contain a t term it can be moved outside of the integral;

$=\frac{\cos (\omega \tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi {)}^{2}dt}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi {)}^{2}dt}+\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\cos (\omega t+\varphi )\sin (\omega \tau )dt}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi {)}^{2}dt}$

The alike terms cancel to give;

$=\cos (\omega \tau )+\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\cos (\omega t+\varphi )\sin (\omega \tau )dt}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi {)}^{2}dt}$

Finally $\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\cos (\omega t+\varphi )\sin (\omega \tau )dt}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi {)}^{2}dt}=0$, this is because $sin$ and $\frac{\cos }{\sin }$ are odd functions. When odd functions are integrated between symmetrical boundaries the product is zero.

Therefore;

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

NJ: Nicely presented.

Analysis of the Lennard-Jones(LJ) solid, liquid and harmonic osillator(HO) VACFs against position

figure 13 VACF plots for LJ solid and liquid simulations and the HO

The minima in the VACFs represent a change in the direction of the velocity, this caused by collisions with other atoms.

The liquid and solid VACFs are similar but have some important differences. The solid VACF reaches its minima much quicker than the liquid; the atoms in the solid are much closer together and so have more surrounding atoms, it therefore takes less time for atoms to collide. The decorrelated solid function shows some oscillation about $C(\tau )=0$ this arises from the jostling about with other atoms in the fixed lattice and so we see a sort of harmonic oscillator motion. The liquid does not display this behaviour because they are free to move and so will be able to "rapidly destroy any oscillatory motion." by the "diffusive motion".

The LJ VACFs become decorrelated after colliding with other atoms, after collision the subsequent velocities no longer depend on the initial velocity. The Harmonic Oscillator’s VACF is fully correlated with the initial velocity. It behaves periodically, always returning through the same points after a certain amount of time. The differences seen between these plots are rooted in the essence of the representations. The HO models two sphere’s connected by a spring, the spring may expand and contract with perfect elasticity. In addition to this there are no external forces which could cause the velocities to decorelate. The LJ representions are effected by surrounding atoms, the velocities loose their correlation after a collision.

NJ: Excellent. What aspect of the liquid structure do you think causes it to have the single broad minimum?

Task 19: Integrating using the trapezium rule to get D

The plots of $C(\tau )vstime$ of the three phases in a small and large system were plotted along with their running integral.

$D=\frac{1}{3}\underset{0}{\overset{\mathrm{\infty }}{\int }}d\tau <v(0)xv(\tau )>=\frac{\frac{V_n+V_{n+1}}{2}x({\mathrm{T}}_{n+1}-{\mathrm{T}}_n)+r_{n-1}}{3}$

The first part of this equation is the relationship between $D$, the diffusion coefficient and the VACF. The second part is the equation used in excel to calculate the running integral. Where $V$ is the velocity, $\mathrm{T}$ is the timestep and $r$ is the last integral calculated. The last value of the running integral was divided by 3 to give $D$.

Solid(figures 14-15)

figure 14 VACF and running integral vs time for 3200 atom solid simulation

$D=1.8x10^{-4}m{s}^{-2}$ for the small solid simulation and $4.55x10^{-5}m{s}^{-2}$ for the million atom simulation. 

figure 15 VACF and running integral vs time for 1000000 atom solid simulation

Liquid(figures 16-17)

[[File:liquid8000.png|200px|thumb|left|figure 16 VACF and running integral vs time for 8000 atom liquid simulation]$D=9.8x10^{-2}m{s}^{-2}$ for the small liquid simulation and $9.0x10^{-2}s{m}^{-2}$ for the million atom simulation.

figure 17 VACF and running integral vs time for 1000000 atom liquid simulation

Gas(figures 18-19)

figure 18 VACF and running integral vs time for 8000 atom gas simulation

$D=2.34m{s}^{-2}$

for the small gas simulation and

$3.27m{s}^{-2}$

for the million atom simulation.

figure 19 VACF and running integral vs time for 1000000 atom gas simulation

The diffusion coefficients are as would be expected for the respective states. An atom in a tight lattice is not free to diffuse and is locked in position, this is reflected in its very small value. The liquid has a more important diffusion coefficient since the atoms are free to move, however they do experience more impeding forces such as drag while passing other moving atoms. Finally the gas is also a reasonable estimation, it is much greater than the solid and liquid since the atoms are a lot more disperse and so move very quickly and encountering relatively few restrictive forces.

Sources of error in the estimates

A larger number of atoms included in the simulation gives a more accurate result, therefore the million atom simulation and the smaller one do have a small difference in values which can be accounted for by this. Another important source of error comes from the running integral calculation, each new area will have the previous area added to it, this results in a build-up of errors. The two solid diffusion coefficients calculated from the running integral are an order larger than the coefficient calculated using the mean squared displacement method.

NJ: How do the MSD and VACF results compare? Again, be careful with units - these should all be reduced.

References: