Talk:Mod:Simulation of a Simple Liquid by Leo Hodges

Introduction to molecular dynamics simulation

Open the file HO.xls. In it, the velocity-Verlet algorithm is used to model the behaviour of a classical harmonic oscillator. Complete the three columns "ANALYTICAL", "ERROR", and "ENERGY": "ANALYTICAL" should contain the value of the classical solution for the position at time $t$, "ERROR" should contain the absolute difference between "ANALYTICAL" and the velocity-Verlet solution (i.e. ERROR should always be positive -- make sure you leave the half step rows blank!), and "ENERGY" should contain the total energy of the oscillator for the velocity-Verlet solution. Remember that the position of a classical harmonic oscillator is given by $x\left(t\right)=A\cos (\omega t+\varphi )$ (the values of $A$, $\omega$, and $\varphi$ are worked out for you in the sheet).

If,

$A=1,\omega =1$ and $\varphi =0$

Then the position using the analytical method is given by:

$x\left(t\right)=\cos (t)$

This is plotted below alongside the Velocity-Verlet, in order to show the comparison:

Figure 1: Position as a function of time using analytical and Velocity-Verlet methods

They are shown to be in very good agreement. The absolute error between the two methods, i.e. the difference between the two, is plotted below:

Figure 2: Difference between analytical and Velocity-Verlet methods

The error between the two methods is ultimately very small but is periodic and increases as total time elapsed increases.

The total energy is given by the sum of the potential and kinetic energies:

$E=\frac{1}{2}m{v}^2+\frac{1}{2}k{x}^2$

where $k=1,m=1$

Figure 3: Total energy as a function of time from the Velocity-Verlet method

The plot shows a wave oscillating around an average value, as expected for an ideal harmonic oscillator. The amplitude of oscillating is small and the total energy remains constant as there is no source of energy loss. NJ: This is a bit confused. Does your graph show constant energy, or not?

For the default timestep value, 0.1, estimate the positions of the maxima in the ERROR column as a function of time. Make a plot showing these values as a function of time, and fit an appropriate function to the data.

Figure 4: Linear correlation of error maxima

NJ: Why do you think the error behaves like this?

Experiment with different values of the timestep. What sort of a timestep do you need to use to ensure that the total energy does not change by more than 1% over the course of your "simulation"? Why do you think it is important to monitor the total energy of a physical system when modelling its behaviour numerically?'

The largest timestep that can be used, such that the amplitude of oscillation does not cause the energy to vary from the average by more than 1%, is 0.285. NJ: This is a bit big, you should have found that 0.2 is the largest value.

Figure 5: Total energy as a function of time from the Velocity-Verlet method with 0.285 timestep

By monitoring the energy of a system, it can be shown that it is constant, as expected for a closed system. NJ: You should say something like, "to see if the energy remains constant", rather than "to show that it is constant". You can't know until you've checked.

For a single Lennard-Jones interaction, $\varphi \left(r\right)=4ϵ(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})$, find the separation, $r_0$, at which the potential energy is zero. What is the force at this separation? Find the equilibrium separation, $r_{eq}$, and work out the well depth ($\varphi \left(r_{eq}\right)$). Evaluate the integrals $\underset{2\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r$, $\underset{2.5\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r$, and $\underset{3\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r$ when $\sigma =ϵ=1.0$.

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

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

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

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

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

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

$r_0=\sigma$

The force is given by:

$F=-\frac{d\varphi \left(r\right)}{dr}$

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

As,

$r_0=\sigma$

Then,

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

$F=4ϵ(\frac{6}{\sigma })$

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

The equilibrium separation, r_e, is the distance where $F=-\frac{d\varphi \left(r\right)}{dr}=0$.

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

$\frac{12{\sigma }^{12}}{{\sigma }^{13}}=\frac{6{\sigma }^6}{{\sigma }^7}$

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

$r_e=2^{\frac{1}{6}}\sigma$

This can then be substituted into the Lennard-Jones equation:

$\varphi \left(r\right)=4ϵ(\frac{{\sigma }^{12}}{2^2{\sigma }^{12}}-\frac{{\sigma }^6}{2{\sigma }^6})$

$\varphi \left(r\right)=4ϵ(\frac{1}{4}-\frac{1}{2})$

$\varphi \left(r\right)=-ϵ$

The integral of the Lennard-Jones potential is:

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

If $\sigma =ϵ=1$

$=4ϵ(\frac{{\sigma }^{12}}{11r^{11}}-\frac{{\sigma }^6}{5r^5}){|}_{2.5\sigma }^{\mathrm{\infty }}=-0.00818$

$=4ϵ(\frac{{\sigma }^{12}}{11r^{11}}-\frac{{\sigma }^6}{5r^5}){|}_{2\sigma }^{\mathrm{\infty }}=-0.0248$

$=4ϵ(\frac{{\sigma }^{12}}{11r^{11}}-\frac{{\sigma }^6}{5r^5}){|}_{3\sigma }^{\mathrm{\infty }}=-0.00329$

Estimate the number of water molecules in 1ml of water under standard conditions. Estimate the volume of 10000 water molecules under standard conditions.

The density of water at 298K is 0.997 g/cm³, therefore the mass of 1 mL of water is 0.997 g. The mass of a water molecule is 2.99 x 10^-23.

$$\begin{gathered} \frac{0.997}{2.99\times 10^{-23}}=3.33\times 10^{23} \\ molecules \end{gathered}$$

10,000 water molecules would weight 2.99 x 10^-19 and would occupy 2.99 x 10^-19 cm³

Consider 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)$. In a single timestep, it moves along the vector $(0.7,0.6,0.2)$. At what point does it end up, after the periodic boundary conditions have been applied?.

The atom would be at position $(0.2,0.1,0.7)$.

The Lennard-Jones parameters for argon are $$\begin{gathered} \sigma =0.34\mathrm{n}\mathrm{m},ϵ \\ / \\ {k}_B=120\mathrm{K} \end{gathered}$$. If the LJ cutoff is $r^{*}=3.2$, what is it in real units? What is the well depth in $$\begin{gathered} {\mathrm{k}\mathrm{J} \\ \mathrm{m}\mathrm{o}\mathrm{l}}^{-1} \end{gathered}$$? What is the reduced temperature $T^{*}=1.5$ in real units?

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

Therefore,

$3.2\times 0.34=1.088nm$

The well depth is given by $ϵ=120K/cdot1.38\times 10^{-23}J{K}^{-1}=1.656\times 10^{-21}J$ NJ: And what's this in kJ mol^-1

And finally,

$T^{*}=\frac{k_{B}T}{ϵ}$

Therefore:

$\frac{1.5\times 1.656\times 10^{-21}}{1.38\times 10^{-23}}=180K$

NJ: Good, and nicely presented.

Equilibration

Why do you think giving atoms random starting coordinates causes problems in simulations? Hint: what happens if two atoms happen to be generated close together?

If atoms are initially randomly positioned, it is possible that configurations are created where atoms are placed very close to one another. Such configurations may be physically impossible in reality and give extremely high Lennard-Jones potentials, causing problems and inaccuracies in simulations.

NJ: Can you elaborate a bit more on this? What sort of problems?

Satisfy yourself that this lattice spacing corresponds to a number density of lattice points of 0.8. Consider instead a face-centred cubic lattice with a lattice point number density of 1.2. What is the side length of the cubic unit cell?

Lattice spacing in x,y,z = 1.07722 1.07722 1.07722

A simple cubic unit cell has one lattice point, thus the volume occupied by a lattice point is equal to:

$$\begin{gathered} V=1.07722^3=1.25000 \\ unit{s}^3 \end{gathered}$$

Therefore the number density of lattice points is given by:

$\rho =\frac{1}{1.25000}=0.8$

A face centered unit cell contains four lattice points, so the length of the cubic unit cell with lattice number density of 1.2 is equal to:

$1.2=\frac{4}{L^3}$

$L=\sqrt[3]{\frac{4}{1.2}}=1.494$

Consider again the face-centred cubic lattice from the previous task. How many atoms would be created by the create_atoms command if you had defined that lattice instead?

Created orthogonal box = (0 0 0) to (10.7722 10.7722 10.7722)

The face centred cubic cell has a unit length of 1.494, so the created box will contain:

$$\begin{gathered} (\frac{10.7722}{1.494}{)}^3=374.85 \\ unit \\ cells \end{gathered}$$

A face centred cubic cell has 4 lattice points, thus the number of atoms created would be equal to:

$$\begin{gathered} 374.85\times 4=1500 \\ atoms \end{gathered}$$

NJ: Unfortunately not. The input script you have always creates 1000 unit cells, no matter what type of lattice you ask for. The lattice type and density then specify how large the box is.

Using the LAMMPS manual, find the purpose of the following commands in the input script:

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

The first row defines only one type of atom and that they all have mass 1.0. The second row, pair_style, calculates the Lennard-Jones potential between pairs of atoms. The command has a cutoff argument, setting global cutoffs for all pairs of atoms, i.e. the maximum distance for which the Lennard-Jones potential can be calculated. In this case it is 3.0 but it can be smaller or larger than the dimensions of the simulation box. The final row, pair_coeff, overrides the global cutoff command for a specific pair of atoms.

Figure 6: Velocity Verlet algorithm.

NJ: Not in this case - it only specifies epsilon and sigma.

Given that we are specifying ${\mathbf{x}}_i\left(0\right)$ and ${\mathbf{v}}_i\left(0\right)$, which integration algorithm are we going to use?

The Velocity Verlet algorithm will be used as position and velocity are calculated at the same value of the time variable.

Look at the lines 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

The second line (starting "variable timestep...") tells LAMMPS that if it encounters the text ${timestep} on a subsequent line, it should replace it by the value given. In this case, the value ${timestep} is always replaced by 0.001. In light of this, what do you think the purpose of these lines is? Why not just write:

timestep 0.001
run 100000

This set of code is such that no matter what you define the timestep to be, the same overall amount of time is run. NJ: Good

Make plots of the energy, temperature, and pressure, against time for the 0.001 timestep experiment (attach a picture to your report). Does the simulation reach equilibrium? How long does this take? When you have done this, make a single plot which shows the energy versus time for all of the timesteps (again, attach a picture to your report). Choosing a timestep is a balancing act: the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality; short timesteps, however, mean that the same number of simulation steps cover a shorter amount of actual time, and this is very unhelpful if the process you want to study requires observation over a long time. Of the five timesteps that you used, which is the largest to give acceptable results? Which one of the five is a particularly bad choice? Why?

Figure 7: Plots of time against energy, temperature and pressure for the 0.001 timestep

If only the first second is plotted, it can be seen that the system reaches equilibrium in about 0.3 seconds: NJ: 0.3 reduced time units, be careful

Figure 8: One second plots of time against energy, temperature and pressure for the 0.001 timestep

A plot of time versus energy for timesteps from 0.001 to 0.015 helps us choose the most appropriate timestep:

Figure 9: Plot of time against energy for varying timesteps

The 0.015 timestep is clearly a bad choice because the system does not equilibrate and the energy gradually increases. The other timesteps all equilibrate and all look similar. A plot of the first second will differentiate between the remaining timesteps:

Figure 10: One second plot of time against energy for varying timesteps

0.01 is the largest timestep that equilibrates, however, 0.0025 appears to be a good choice because it matches the 0.001 timestep very closely, the smallest timestep but does not take as much time or oscillate as much. NJ: Good - a little more explanation would be nice (e.g. you know that theoretically 0.001 is more accurate, but this is not detectable at this precision.

Running simulations under specific conditions

Choose 5 temperatures (above the critical temperature $T^{*}=1.5$), and two pressures (you can get a good idea of what a reasonable pressure is in Lennard-Jones units by looking at the average pressure of your simulations from the last section). This gives ten phase points — five temperatures at each pressure. Create 10 copies of npt.in, and modify each to run a simulation at one of your chosen $(p,T)$ points. You should be able to use the results of the previous section to choose a timestep. Submit these ten jobs to the HPC portal. While you wait for them to finish, you should read the next section.

The chosen temperatures were T = 1.5, 2.0, 2.5, 3.0 and 6.0 at pressures p = 2.5 and 3.0. These simulations are to be performed with a 0.0025 timestep, owing to the previous conclusion.

We need to choose $\gamma$ so that the temperature is correct $T=\mathfrak{𝔗}$ if we multiply every velocity $\gamma$. We can write two equations:

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

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

Solve these to determine $\gamma$.

Equating the equations gives us:

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

So,

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

Therefore:

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

Use the manual page to find out the importance of the three numbers 100 1000 100000. How often will values of the temperature, etc., be sampled for the average? How many measurements contribute to the average? Looking to the following line, how much time will you simulate?

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

The penultimate line of code, starting fix aves, takes the defined inputs, such as density and pressure and averages them every set number of timesteps. How it does this is dictated by the numbers given in the line of code. These numbers represent Nevery, Nrepeat and Nfreq respectively. The Nevery number is after every however many timesteps an input value is taken, in this case 100. Nrepeat is how many times the input values are used for calculating averages, 1000, and Nfreq calculates an average every this many timesteps, 100000 for this code.

An average is taken every 100000 timesteps and an input value is taken every 100 timesteps, so 1000 measurements will contribute to the average. With a timestep of 0.0025 and 100000 timesteps to be run, 250 seconds will be simulated overall.

NJ: Not seconds! These units are reduced.

When your simulations have finished, download the log files as before. At the end of the log file, LAMMPS will output the values and errors for the pressure, temperature, and density $\left(\frac{N}{V}\right)$. Use software of your choice to plot the density as a function of temperature for both of the pressures that you simulated. Your graph(s) should include error bars in both the x and y directions. You should also include a line corresponding to the density predicted by the ideal gas law at that pressure. Is your simulated density lower or higher? Justify this. Does the discrepancy increase or decrease with pressure?

Figure 11: Density as a function of temperature plot

NJ: Nice plotThe plot for the simulated density as a function of temperature of the liquid at different pressures follow very similar curves. As expected, the higher pressure gives a higher density along the curve. Error bars are included on both axes but are small; to the order of around 10^-2 for density. The curves for the liquid at both pressures are lower than their respective ideal gas equivalent. This can be rationalised due to the fact that an ideal gas is considered to have no significant particle-particle interactions, while the liquid simulation takes into account Lennard-Jones interactions, so particles are more likely to be found further apart from one another. NJ: Actually, the ideal gas has no interactions between particles at all'. Can you say more about this? If you used a purely attractive potential between the particles, would you see the same thing? The discrepancy can be seen to increase with pressure, this fits with the conclusion. NJ: Why?

As temperature increases, the liquid and ideal gas curves converge, however. As the temperature increases, thermal kinetic energy overrides the intermolecular interaction effects and the liquid becomes more like an ideal gas. NJ: Good. Thermal energy and kinetic energy are essentially interchangeable terms.

Calculating heat capacities using statistical physics

As in the last section, you need to run simulations at ten phase points. In this section, we will be in density-temperature $({\rho }^{*},T^{*})$ phase space, rather than pressure-temperature phase space. The two densities required at $0.2$ and $0.8$, and the temperature range is $2.0,2.2,2.4,2.6,2.8$. Plot $C_V/V$ as a function of temperature, where $V$ is the volume of the simulation cell, for both of your densities (on the same graph). Is the trend the one you would expect? Attach an example of one of your input scripts to your report.

Many quantities in statistical thermodynamics can be calculated by considering how far the system is able to fluctuate from its average equilibrium state. For example, if the temperature is held "constant", then we know that the total energy must be fluctuating. The magnitude of these fluctuations actually tell us about the heat capacity of the system, according to the equation

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

$\mathrm{V}\mathrm{a}\mathrm{r}\left[E\right]$ is the variance in $E$, $N$ is the number of atoms, and it is a standard result from statistics that $\mathrm{V}\mathrm{a}\mathrm{r}\left[X\right]=⟨X^2⟩-{⟨X⟩}^2$. The $N^2$ is required because LAMMPS divides all energy output by the number of particles (so when you measure $E$, you are actually measuring $E/N$).

The script used to perform the simulation for density 0.8 and temperature 2.0 is shown below:

### DEFINE SIMULATION BOX GEOMETRY ###
variable d equal 0.8
lattice sc ${d}
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.0
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
fix nvt all nvt temp ${T} ${T} ${tdamp}
run 10000
reset_timestep 0

### MEASURE SYSTEM STATE ###
unfix nvt
fix nve all nve
thermo_style custom step etotal temp
variable e equal etotal
variable e2 equal etotal*etotal
variable temp equal temp
fix aves all ave/time 100 1000 100000 v_temp v_e v_e2
run 100000

variable avetemp equal f_aves[1]
variable avetemp2 equal f_aves[1]*f_aves[1]
variable avgesquare equal f_aves[3]
variable squareavge equal f_aves[2]*f_aves[2]
variable shola equal ${avgesquare}-${squareavge}
variable Cv equal atoms*atoms*${shola}/${avetemp2}
variable vol equal vol
variable CvV equal ${Cv}/${vol}

print "Temperature: ${avetemp}"
print "Cv: ${Cv}"
print "Cv/V: ${CvV}"

The data from the ten simulations could then be plotted on the same graph in Excel:

Figure 12: Cv/V as a function of temperature plot

NJ: Don't use the Excel "smoothed lines" feature for scientific plots. Do you think the maximum around T=2.4 is a real effect, or an artifact of the smoothing? When you only have 5 points, it's fine to display only the points.

The heat capacity of the system can be seen to decrease as the temperature increases for both densities. This is an observation that can be rationalised by quantum mechanics. As the temperature increases, higher energy levels start to populate and consequently reduce the ability of the system to absorb further energy. This reduces the heat capacity. NJ: This isn't a quantum system though, and there isn't a maximum energy level. You're right that you can use these ideas to rationalise the trend - the fluid has modes which are analogous to quantum energy levels. The decrease in heat capacity shows that they become more closely spaced in energy at higher temperatures. The amount of energy required to move from one mode to the next thus decreases as temperature increases, and so does the heat capacity.

As for the difference between the trends for the two densities, the higher heat capacity for the greater density is also to be expected. A greater number of particles in an equivalent volume would be able to absorb more energy than a system with a lower density. NJ: Can you say any more about this? Why is this the case?

Structural properties and the radial distribution function

Figure 13: Coexistence curve for the Lennard-Jones system (temperatures and densities in reduced units)

Perform simulations of the Lennard-Jones system in the three phases. When each is complete, download the trajectory and calculate

$g(r)$

and

$\int g(r)\mathrm{d}r$

. Plot the RDFs for the three systems on the same axes, and attach a copy to your report. Discuss qualitatively the differences between the three RDFs, and what this tells you about the structure of the system in each phase. In the solid case, illustrate which lattice sites the first three peaks correspond to. What is the lattice spacing? What is the coordination number for each of the first three peaks?

The density and temperature parameters for the three phases were chosen based on the phase diagram shown in Figure 13. The density of the gas, liquid and solid were 0.05, 0.8 and 1.2 respectively, while the temperature for all the phases was 1.2.

Figure 14: Radial Distribution Function for three phases (left) and solid (right)

Figure 15: Face centred cubic unit cell

The radial distribution function describes how the density in a system of particles varies as a function of distance from a reference particle, in this case a randomly chosen atom. Essentially, it is the probability of finding a particle at a certain distance away from the reference particle, relative to an ideal gas.

The densities of the solid and the liquid are similar, so the separation of molecules would be expected to be similar. NJ: The solid is 50% denser... This is the case but the fundamental difference between the two phases is the presence of long-range order in the solid state that is not seen in the liquid state. Three peaks can be seen for the liquid phase at low internuclear distances before the amplitude diminishes to 1. The short range order observed, despite the constant motion of particles in a liquid can be explained by coordination shells. The distances between the randomly chosen particle and the particles in the first coordination shell is much shorter than for the second coordination shell, which in turn is shorter than the third until there is essentially no coordination, hence the amplitude of 1.

For the face centred cubic solid state, however, the first three peaks are found at similar internuclear distances as that of the liquid but the amplitude does not diminish to the same extent and the peaks are more clearly defined than in the liquid. This is because in a solid, the atoms are in fixed positions and can be defined consistently regardless of the reference particle.

In the case of the gas there is only one peak due to very little 'local structure'. At distances shorter than the atomic diameter, g(r) is zero due to strong repulsive forces.

NJ: Very nice explanation. I would move the discussion of r < sigma to the very start, to make it clear that it applies to all three phases.

Due to the fixed nature of atoms in the face centred cubic lattice, the internuclear distance of the peaks, as shown in Figure 14, can be used to find the lattice spacing. Figure 15 shows the internuclear distances the peaks correspond to, in the context of a face centred cubic unit cell.

NJ: Earlier you works out that the spacing for this very lattice should be around 1.49. Why do you think this is larger? You should use all three peaks to get an "average" lattice spacing to minimise the error in reading off the graph. What about the coordination numbers?

Dynamical properties and the diffusion coefficient

In the D subfolder, there is a file liq.in that will run a simulation at specified density and temperature to calculate the mean squared displacement and velocity autocorrelation function of your system. Run one of these simulations for a vapour, liquid, and solid. You have also been given some simulated data from much larger systems (approximately one million atoms). You will need these files later.

As before, the density of the vapour is 0.05, the liquid is 0.8 and the solid 1.2.

Make a plot for each of your simulations (solid, liquid, and gas), showing the mean squared displacement (the "total" MSD) as a function of timestep. Are these as you would expect? Estimate $D$ in each case. Be careful with the units! Repeat this procedure for the MSD data that you were given from the one million atom simulations.

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

Figure 16: Mean square displacement for the vapour phase

Ultimately, all of the phases show a linear correlation between the mean square displacement and the number of timesteps. In the case of the vapour phase, this is after an initial 1500 timesteps of non-linearity. NJ: What do you think the reason for this is?This allows us to use the gradient of the slope in the equation shown above to estimate the diffusion coefficient, D. As the differential is with respect to time, t, we must divide the gradient by the length of the timestep, 0.002. With a gradient of 0.0374, we obtain for the diffusion coefficient:

$D=\frac{1}{6}\frac{0.0374}{0.002}=3.11$

For 1 million atoms, we obtain a very similar coefficient:

$D=\frac{1}{6}\frac{0.037}{0.002}=3.08$

Figure 17: Mean square displacement for the liquid phase

The liquid phase has a very close linear trend from the very first timestep. It also exhibits a much smaller overall mean square displacement at the final timestep compared to the vapour, thus also a smaller gradient. The value for the coefficient is calculated to be:

$D=\frac{1}{6}\frac{0.001}{0.002}=0.083$

The same diffusion coefficient is obtained for 1 million atoms as the gradient of the slope is the same.

Figure 18: Mean square displacement for the solid phase

Finally, the solid phase shows linearity after around 250 timesteps but in this case the gradient is very close to zero. A very small overall mean square displacement is to be expected for a face centred cubic lattice with fixed atoms. The diffusion coefficient is calculated to be:

$D=\frac{1}{6}\cdot \frac{1\times 10^{-8}}{0.002}=8.3\times 10^{-7}$

For 1 million atoms:

$D=\frac{1}{6}\cdot \frac{1\times 10^{-8}}{0.002}=8.3\times 10^{-6}$

It can also be seen that throughout the calculations of the three phases, simulating one million atoms provided very little benefit and is unnecessary. The only noticeable difference being the smoother gradient in the mean square displacement plot for the solid phase. NJ: Do you think the uncertainties in the measurements are the same for both numbers of atoms?

We can also calculate the diffusion coefficient using the velocity autocorrelation function, which we denote by

$C\left(\tau \right)=⟨\mathbf{v}\left(t\right)\cdot \mathbf{v}(t+\tau )⟩$

In essence, this function tells us how different the velocity of an atom will be at time $t+\tau$ to its velocity at time $t$. At long times, we expect $C\left(\mathrm{\infty }\right)=0$, because the velocities of atoms should become uncorrelated at very long times. By this, we mean that the velocity of an atom at a particular time should not depend on its velocity a long time in the past — the exact form of $C\left(\tau \right)$ tells us how long "a long time" is. Remarkably, we can calculate the diffusion coefficient by integrating this function:

$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 "optional output-3" file from your simulations contains measurements of the VACF. The format of the files is similar to the format of the MSD data files: the first column contains the timestep at which the VACF data was evaluated (this is $\tau$, in timestep units), the next three contain $C\left(\tau \right)$ for the Cartesian components of the velocity, and the final column contains the VACF for the whole velocity.

In the theoretical section at the beginning, the equation for the evolution of the position of a 1D harmonic oscillator as a function of time was given. Using this, evaluate the normalised velocity autocorrelation function for a 1D harmonic oscillator (it is analytic!):

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

Be sure to show your working in your writeup. On the same graph, with x range 0 to 500, plot $C\left(\tau \right)$ with $\omega =1/2\pi$ and the VACFs from your liquid and solid simulations. What do the minima in the VACFs for the liquid and solid system represent? Discuss the origin of the differences between the liquid and solid VACFs. The harmonic oscillator VACF is very different to the Lennard Jones solid and liquid. Why is this? Attach a copy of your plot to your writeup.

The positiion of a classical harmonic oscillator is given by:

$x(t)=A\cos (\omega t+\varphi )$

If

$t=t+\tau$ NJ: Be careful with your notation here. That equation states that $\tau =0$.

Then,

$x(t)=A\cos (\omega (t+\tau )+\varphi )$

This can be differentiated with respect to t, using the chain rule, to give us the velocity as:

$\frac{dx(t)}{dt}=v(t)=-A\omega \sin (\omega (t+\tau )+\varphi )$ NJ: Careful again. On the left you have a function of t only, and on the right one that involves $\tau$

This can be substituted into:

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

To give:

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

We can use the trigonometric identity:

$sin(u+v)=sin(u)cos(v)+cos(u)sin(v)$

To give:

$C\left(\tau \right)=\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 }^2(\omega t+\varphi )dt}$

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

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

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

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

The bottom integral does not converge but the top integral is a combination of an even and an odd function, therefore as long as the limits are equal and opposite, the integral will tend to zero.

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

Therefore:

$C\left(\tau \right)=\cos (\omega \tau )$ NJ: Good. A nicely present derivation, just be careful in future with those two lines at the start.

Figure 19: Velocity autocorrelation function for the solid and liquid phases and the harmonic oscillator

The minima on the plots represent collisions with other atoms. The solid plot shows a decaying nature after the initial collision on the solid plot due to its high degree of order, while the liquid exhibits only one minima. This represents a collision with the solvent cage. There are no collisions in the case of the harmonic oscillator, so the amplitude of the periodic wave remains constant.

NJ: Can you say a bit more about this? Why does a collision lead to a minimum? If the harmonic oscillator has no collisions, why does it too have minima?

Use the trapezium rule to approximate the integral under the velocity autocorrelation function for the solid, liquid, and gas, and use these values to estimate $D$ in each case. You should make a plot of the running integral in each case. Are they as you expect? Repeat this procedure for the VACF data that you were given from the one million atom simulations. What do you think is the largest source of error in your estimates of D from the VACF?

A running integral of the VACF can be approximated using the trapezium rule, allowing us the estimate the diffusion coefficient with the below equation:

$D=\frac{1}{6}\cdot \frac{1\times 10^{-8}}{0.002}=8.3\times 10^{-6}$

Figure 20: Running integral of the velocity autocorrelation function for the vapour phase

The diffusion coefficient works out to be 3.294 and for the 1 million atom system, 3.268. This is similar to the result obtained above from the mean square displacement method.

Figure 21: Running integral of the velocity autocorrelation function for the liquid phase

A value of 0.0977 is obtained for the liquid phase diffusion coefficient and 0.09 for the equivalent 1 million atom system, which is again similar to the previous method.

Figure 22: Running integral of the velocity autocorrelation function for the solid phase

The diffusion coefficient for the solid is again essentially zero which is to be expected. A value of 1.84 x 10^-4 is calculated for the diffusion coefficient and 4.53 x 10^-5 for the 1 million atom system. It can again be seen that there is little difference in the result from the 1 million atom simulation and it is therefore unnecessary to simulate such a large number of atoms. More accurate estimations of the diffusion coefficient could be obtained by using more sophisticated methods of approximating the integral.

NJ: It would be good to see a side-by-side comparison of the MSD and VACF results. Which do you think is more accurate? Do you think the integration method is really the largest source of error?