Report

Abstract

Simulations were ran on the LAMMPS application for LJ solids, liquids and vapours. The results show that there is a large discrepency between simulated and ideal density because it doesn't take into account the interactions between atoms. Denser materials have a higher heat capacity because there is more atoms to store energy. The RDF show that solids have a fixed lattice structure, liquids have a more dynamic shell and vapours only have shells at very close interatomic distances. Through calculating diffusion coefficient using MSD and VACF, it shows that D for vapour>liquid>solid.

Introduction

Molecular dynamics simulation is "widely used to study a huge variety of chemical phenomena, from the behaviour of exotic materials under extreme conditions to protein folding and the properties of biological systems such as lipid membranes".[1] A recent discovery show that water diffusion between graphene oxide layers is slow because of the strong hydrogen bonding interactions. This causes most of the water molecules are bound to OH groups even at the highest hydration level. These graphene oxide layers will potentially revolutionize water filtration and desalination if selective transport of molecules can be controlled.[2] This experiment is going to investigate how number density and heat capacity vary with temperature. The structure of solid, liquid and vapour phase will be determined using RDF and the diffusion coefficient will be calculated using MSD and VACF.

Methods

All the simulations were ran through the LAMMPS application via the HPC portal.

Equation of State:

A simple cubic lattice was generated with 3375 atoms of mass 1.0, in a LJ forcefield with a cut-off distance of 3.0, under the NpT ensemble. The densities used were 2.8 and 1.8. The temperatures used were 1.8, 2.0, 2.2, 2.4 and 2.6. A timestep of 0.0025 was used because it gave accurate results without taking a long time to simulate.

Heat Capacity:

A simple cubic lattice was generated with 3375 atoms of mass 1.0, in a LJ forcefield with a cut-off distance of 3.0, under the NVT ensemble. The densities used were 0.2 and 0.8. The temperatures used were 2.0, 2.2, 2.4, 2.6 and 2.8. The input file could be found in TASKS. A timestep of 0.0025 was used.

Structural Properties and the Radial Distribution Function:

For the solid phase, a face-centered cubic lattice was generated with 8000 atoms of mass 1.0, in a LJ forcefield with a cut-off distance of 3.2. A density of 1.2 and temperature of 1.0 was used.

For the liquid and vapour phase, a simple cubic lattice was generated under the same parameters except for density, with a value of 0.8 and 0.01 respectively.

Dynamical Properties and the Diffusion Coefficient:

Same input file as the above was used.

Velocity Autocorrelation Function:

Same input file as the above was used.

Results and Discussion

Equation of State:

The equation of state is given by: $ρ = \frac{N}{V} = \frac{P}{k_{B} T}$

An ideal gas is considered to have no interatomic interactions.

The simulated density was lower than because it takes the LJ pair potential between atoms into account, therefore atoms can't be as close together as an ideal gas. The discrepancy increases with pressure because LJ pair potential increases when distance between atoms decreases therefore, there will be a higher repulsive force when the atoms are closer together due to higher pressure.

Heat Capacity:

Heat capacity is given by: $C_{V} = \frac{\partial E}{\partial T} = \frac{V a r [E]}{k_{B} T^{2}} = N^{2} \frac{⟨ E^{2} ⟩ - {⟨ E ⟩}^{2}}{k_{B} T^{2}}$

It is expected that heat capacity will increase with density because there is more molecules in the same volume so more energy can be stored. It is also expected that the heat capacity decreases with temperature because at higher temperature, higher energy levels starts getting occupied therefore the system won't be able to store as much energy.

Structural Properties and the Radial Distribution Function:

RDF provides information about the distances of the nearest neighbour, second nearest neighbour and etc relative to a certain origin in a solid lattice. The integral of the RDF gives the coordination number, the number of identical atoms with the same lattice spacing from the origin.

Solid phase's RDF fluctuates because it has an ordered structure so the atoms remain in roughly the same place relative to a certain origin and form distinct shells. The peaks refers to the location of lattice points in a face-centered cubic lattice. The RDF decreases with distance because the number of atoms in the shell increases at a slower rate than the increase in volume.

Liquid phase's RDF fluctuates less than solids because it has a less ordered structure and atoms are able to move around therefore the shells are less distinct and more spaced out.

In the vapour phase, atoms are not ordered and far apart therefore only one initial peak is shown because LJ pair potential only has a significant effect on small internuclear distances. At higher internuclear distances, RDF=1 because of the random movement of atoms.

The first peak corresponds to atom A, with a lattice spacing of 1.025 and a coordination number of 12. The second peak corresponds to atom B, with a lattice spacing of 1.475 and a coordination number of 18. The third peak corresponds to atom C, with a lattice spacing of 1.825 and a coordination number of 42.

Dynamical Properties and the Diffusion Coefficient:

The Mean Squared Displacement, MSD, is the average distance moved by an atom. It can be used to calculate the diffusion coefficient $D$ , using: $D = \frac{1}{6} \frac{\partial ⟨ r^{2} (t) ⟩}{\partial t}$

These are for 8000 atoms:

It is expected that the vapour phase will have the highest MSD because atoms are able to move freely. Solid phase has the lowest MSD because atoms are fixed in a lattice therefore there won't be much movement except for thermal vibrations. Liquid state has MSD values that are higher than but close to solids because it is able to move around but has similar density than solids.

The increase in MSD is parabolic initially for vapour phase because there is no collision due to its low density but it turns linear when collisions start to happen. Liquid phase's MSD increases linearly because it has similar density to solids so there is a lot of collisions as it travels.

These are for 1 million atoms:

The simulation with 1 million atoms has less fluctuation because the MSD are averaged over a lot more atoms, resulting in a more accurate representation of the diffusion coefficient, $D$ .

Number of Atoms	Solid	Liquid	Vapour
8000	9.57x10^-9	1.33x10^-4	8.07x10^-3
1 million	8.78x10^-9	1.67x10^-4	5.08x10^-3

Velocity Autocorrelation Function:

The position of a classical harmonic oscillator is given by:

$x (t) = A c o s (ω t + ϕ)$ .

By differentiating $x (t)$ with respect to time, it gives the velocity:

$v (t) = \frac{d x (t)}{d t} = - A ω s i n (ω t + ϕ)$ .

The velocity autocorrelation function is given by:

$C (τ) = ⟨ v (t) \cdot v (t + τ) ⟩$

The normalised function is given by:

$C (τ) = \frac{\int_{- \infty}^{\infty} v (t) v (t + τ) d t}{\int_{- \infty}^{\infty} v^{2} (t) d t}$

By substituting the velocity into the normalised function, it gives:

$C (τ) = \frac{\int_{- \infty}^{\infty} s i n (ω t + ϕ) s i n (ω (t + τ) + ϕ) d t}{\int_{- \infty}^{\infty} (s i n (ω t + ϕ))^{2} d t}$

Using the sine addition formula for $t + τ$ , it gives:

$C (τ) = \frac{\int_{- \infty}^{\infty} s i n (ω t + ϕ) (s i n (ω t + ϕ) c o s (ω τ) + s i n (ω τ) c o s (ω t + ϕ)) d t}{\int_{- \infty}^{\infty} (s i n (ω t + ϕ))^{2} d t} = c o s (ω τ) + \frac{\int_{- \infty}^{\infty} s i n (ω t + ϕ) s i n (ω τ) c o s (ω t + ϕ) d t}{\int_{- \infty}^{\infty} (s i n (ω t + ϕ))^{2} d t}$

$\frac{\int_{- \infty}^{\infty} s i n (ω t + ϕ) s i n (ω τ) c o s (ω t + ϕ) d t}{\int_{- \infty}^{\infty} (s i n (ω t + ϕ))^{2} d t} = 0$

Therefore:

$C (τ) = c o s (ω τ)$

The diffusion coefficient $D$ can be calculated by integrating the VACF:

$D = \frac{1}{3} \int_{0}^{\infty} d τ ⟨ v (0) \cdot v (τ) ⟩$

The minima in the VACF plot represents the collision with other atoms. In solid phase, the atoms can't move around and can only vibrate because there will be a restoring force, as shown by the fluctuations. However, in liquid phase, atoms can move around so there won't be much restoring force to move it back in place, as shown by the rapid damping of the velocity.

The harmonic oscillator VACF is different to the LJ solid and liquid because it doesn't take external forces from other atoms into account. It only considers an single atom that has a restoring force that is proportional to its displacement from the equilibrium position.

Integration of the VACF using the trapezium rule ( $A r e a = \frac{1}{2} \times w i d t h \times [($ 1st height $) + 2 ($ sum of other heights $)_{n - 1} + ($ last height $)_{n}]$ ) can be used to calculate the diffusion coefficient using: $D = \frac{1}{3} \int_{0}^{\infty} C (τ) d τ$

These are for 8000 atoms:

These are for 1 million atoms:

Number of Atoms	Solid	Liquid	Vapour
8000	1.57x10^-4	8.19x10^-2	6.82
1 million	4.79x10^-5	9.27x10^-2	3.49

These values are as expected because $D$ is greatest for vapour phase then liquid and solid for the same reasons as above.

The largest source of error in the estimate of $D$ for VACF will be the use of the trapezium rule, therefore some areas will not be accounted for, leading to a smaller than expected value. The accuracy can be increased by using smaller trapeziums.

Conclusion

The most useful timestep was found to be 0.0025 because it gives accurate results while not taking a long time to simulate. The results show that there is a large discrepency between simulated and ideal density because it doesn't take into account the interactions between atoms. Denser materials have a higher heat capacity because there is more atoms to store energy. The RDF show that solids have a fixed lattice structure, liquids have a more dynamic shell and vapours only have shells at very close interatomic distances. The integral of RDF in solid phase shows that there are 12 atoms that has a spacing of 1.025, 18 atoms with spacing of 1.475 and 42 atoms with a spacing of 1.825. Through calculating diffusion coefficient using MSD and VACF, it shows that D for vapour>liquid>solid. For future experiments, more than one atom type could be used to see the effect that it has on the RDF. Improvements could be made by using smaller trapeziums when integrating the VACF.

Reference

[1] https://wiki.ch.ic.ac.uk/wiki/index.php?title=Third_year_simulation_experiment

[2] https://www.nature.com/articles/srep29484?WT.feed_name=subjects_physical-sciences

[3] http://lammps.sandia.gov

TASKS

Theory

TASK: 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 (t) = A \cos (ω t + ϕ)$ (the values of $A$ , $ω$ , and $ϕ$ are worked out for you in the sheet).

TASK: 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.

TASK: 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?

A timestep of 0.0025 was used. It is important to monitor the total energy because of conservation of energy, therefore total energy should be constant.

TASK:For a single Lennard-Jones interaction, $ϕ (r) = 4 ϵ (\frac{σ^{12}}{r^{12}} - \frac{σ^{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_{e q}$ , and work out the well depth ( $ϕ (r_{e q})$ ). Evaluate the integrals $\int_{2 σ}^{\infty} ϕ (r) d r$ , $\int_{2.5 σ}^{\infty} ϕ (r) d r$ , and $\int_{3 σ}^{\infty} ϕ (r) d r$ when $σ = ϵ = 1.0$ .

When the potential energy is zero:

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

$\frac{σ^{12}}{r_{o}^{12}} - \frac{σ^{6}}{r_{o}^{6}} = 0$

$\frac{σ^{12}}{r_{o}^{12}} = \frac{σ^{6}}{r_{o}^{6}}$

$σ^{6} = r_{o}^{6}$

$σ = r_{o}$

To work out the force at this separation:

$F = - \frac{d ϕ (r)}{d r}$

$ϕ^{'} = 4 ϵ (- \frac{12 σ^{12}}{r^{13}} + \frac{6 σ^{6}}{r^{7}})$

Because $σ = r_{o}$ , therefore

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

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

At equilibrium, the net force will be zero therefore:

$F = - 4 ϵ (- \frac{12 σ^{12}}{r_{e q}^{13}} + \frac{6 σ^{6}}{r_{e q}^{7}}) = 0$

$\frac{12 σ^{12}}{r_{e q}^{13}} = \frac{6 σ^{6}}{r_{e q}^{7}}$

$r_{e q}^{6} = 2 σ^{6}$

$r_{e q} = 2^{\frac{1}{6}} σ$

To work out the well depth:

$ϕ (r_{e q}) = 4 ϵ (\frac{σ^{12}}{r_{e q}^{12}} - \frac{σ^{6}}{r_{e q}^{6}})$

$ϕ (r_{e q}) = 4 ϵ (\frac{σ^{12}}{4 σ^{12}} - \frac{σ^{6}}{2 σ^{6}})$

$ϕ (r_{e q}) = 4 ϵ (\frac{1}{4} - \frac{1}{2})$

$ϕ (r_{e q}) = - ϵ$

Evaluation of the integrals:

$\int_{2 σ}^{\infty} ϕ (r) d r = {[- \frac{4 r^{- 11}}{11} + \frac{4 r^{- 5}}{5}]}_{2 σ}^{\infty}$

$= [\frac{4}{11} \frac{1}{2^{11} σ^{11}} - \frac{4}{5} \frac{1}{2^{5} σ^{5}}] = - 0.0248$

$\int_{2.5 σ}^{\infty} ϕ (r) d r = {[- \frac{4 r^{- 11}}{11} + \frac{4 r^{- 5}}{5}]}_{2.5 σ}^{\infty}$

$= [\frac{4}{11} \frac{1}{2 . 5^{11} σ^{11}} - \frac{4}{5} \frac{1}{2 . 5^{5} σ^{5}}] = - 0.00817$

$\int_{3 σ}^{\infty} ϕ (r) d r = {[- \frac{4 r^{- 11}}{11} + \frac{4 r^{- 5}}{5}]}_{3 σ}^{\infty}$

$= [\frac{4}{11} \frac{1}{3^{11} σ^{11}} - \frac{4}{5} \frac{1}{3^{5} σ^{5}}] = - 0.00329$

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

Number of molecules in 1ml of water:

$n = \frac{m}{M_{r}} = \frac{1}{18} = 0.0555 m o l$

$N = n \times N_{A} = 0.0555 \times 6.022 \times 1 0^{23} = 3.34 \times 1 0^{22} m o l e c u l e s$

Volume of 10000 water molecules:

$n = \frac{N}{N_{A}} = \frac{10000}{6.022 \times 1 0^{23}} = 1.66 \times 1 0^{- 20} m o l$

$m = n \times M_{r} = 1.66 \times 1 0^{- 20} \times 18 = 2.99 \times 1 0^{- 19} g = 2.99 \times 1 0^{- 19} m l$ because density of water is 1g/ml

TASK: 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?.

After periodic boundary conditions have been applied, the new position is:

$(0.5, 0.5, 0.5) + (0.7, 0.6, 0.2) = (1.2, 1.1, 0.7) = (0.2, 0.1, 0.7)$

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

The LJ cutoff distance in real units:

$r = r^{*} \times σ = 3.2 \times 0.34 \times 1 0^{- 9} = 1.09 \times 1 0^{- 9} m$

The well depth in real units:

$ϵ = 120 \times k_{B} = \frac{1.66 \times 1 0^{- 21} \times 6.02 \times 1 0^{23}}{1000} = 0.997 k J m o l^{- 1}$

The temperature in real units:

$T = \frac{T^{*} ϵ}{k_{B}} = \frac{1.5 \times 1.66 \times 1 0^{- 21}}{1.38 \times 1 0^{- 23}} = 180 K$

Equilibration

TASK: 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?

When two atoms are generated close together, there will be a very high repulsive force because potential increases exponentially as the distance between the atoms approaches zero. This will lead to particles having a very highly velocity so unless there is a very small timestep, the simulation won't be able to model the behavior.

TASK: 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?

In a simple cubic lattice, there is only 1 lattice point per unit cell, therefore:

$D = \frac{N}{L^{3}} = \frac{1}{1.0772 2^{3}} = 0.800$

In a face-centred cubic lattice, there is 4 lattice point per unit cell, therefore:

$L = {(\frac{4}{1.2})}^{\frac{1}{3}} = 1.494$

TASK: 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?

The create_atom command will create 4000 atoms because there is 4x the lattice point per unit cell.

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

mass 1 1.0 = Atom type 1 with a mass of 1.0 in reduced units

pair_style lj/cut 3.0 = Lennard-Jones Potential between atoms with no coulombic interactions and a cut-off distance of 3.0

pair_coeff * * 1.0 1.0 = Pairwise force field coefficient between all types of atoms is 1.0

TASK: Given that we are specifying $x_{i} (0)$ and $v_{i} (0)$ , which integration algorithm are we going to use?

The Velocity-Verlet integration algorithm will be used because the initial postion and velocity is specified.

TASK: Look at the lines below.

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

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

By defining the variable "timestep", the number of simulation that it will run will change automatically according to the change in timestep, while keeping the total time constant.

TASK: 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?

The simulation reaches equilibrium in a time of 0.4.

The largest timestep to give acceptable results is 0.01. The 0.015 timestep was a bad choice because it doesn't converge and the total energy kept increasing.

Running simulations under specific conditions

TASK: 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.

TASK: We need to choose $γ$ so that the temperature is correct $T = 𝔗$ if we multiply every velocity $γ$ . 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} {(γ v_{i})}^{2} = \frac{3}{2} N k_{B} 𝔗$

Solve these to determine $γ$ .

$γ^{2} \frac{1}{2} \sum_{i} m_{i} v_{i}^{2} = \frac{3}{2} N k_{B} 𝔗$

$γ^{2} \frac{3}{2} N k_{B} T = \frac{3}{2} N k_{B} 𝔗$

$γ = {(\frac{𝔗}{T})}^{\frac{1}{2}}$

TASK: 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?

100 = Nevery = Use input values every 100 timesteps

1000 = Nrepeat = Number of times to use input value to calculate averages

100000 = Nfreq = Calculate averages every 10000 steps

Because each timestep is 0.0025, therefore the simulation will run for 250.

TASK: 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 $(\frac{N}{V})$ . 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?

The simulated density was lower than because it takes the LJ pair potential between atoms into account, therefore atoms can't be as close together as an ideal gas. The discrepancy increases with pressure because LJ pair potential increases when distance between atoms decreases therefore, there will be a higher repulsive force when the atoms are closer together due to higher pressure.

Calculating heat capacities using statistical physics

TASK: As in the last section, you need to run simulations at ten phase points. In this section, we will be in density-temperature $(ρ^{*}, 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.

### DEFINE SIMULATION BOX GEOMETRY ###
variable D equal 0.2

### DEFINE SIMULATION BOX GEOMETRY ###
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 p equal 2.5
variable timestep equal 0.0025

### ASSIGN ATOMIC VELOCITIES ###
 velocity all create ${T} 12345 dist gaussian rot yes mom yes
    
 ### SPECIFY ENSEMBLE ###
 timestep ${timestep}
 fix nve all nve
    
 ### THERMODYNAMIC OUTPUT CONTROL ###
 thermo_style custom time etotal temp press
 thermo 10
    
 ### RECORD TRAJECTORY ###
 dump traj all custom 1000 output-1 id x y z
    
 ### SPECIFY TIMESTEP ###
    
 ### RUN SIMULATION TO MELT CRYSTAL ###
 run 10000
 unfix nve
 reset_timestep 0
    
 ### BRING SYSTEM TO REQUIRED STATE ###
 variable tdamp equal ${timestep}*100
 variable pdamp equal ${timestep}*1000
 fix nvt all nvt temp ${T} ${T} ${tdamp} # system simulated in nvt with a specific density
 run 10000
 unfix nvt
 fix nve all nve
 reset_timestep 0
    
 ### MEASURE SYSTEM STATE ###
 thermo_style custom step etotal temp press density atoms
 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
 variable energy equal etotal
 variable energy2 equal etotal*etotal
 variable no_atoms2 equal atoms*atoms
 fix aves all ave/time 50 2000 100000 v_dens v_temp v_press v_dens2 v_temp2 v_press2 v_energy v_energy2
 run 100000
    
 variable avedens equal f_aves[1]
 variable avetemp equal f_aves[2]
 variable avepress equal f_aves[3]
 variable errdens equal sqrt(f_aves[4]-f_aves[1]*f_aves[1])
 variable errtemp equal sqrt(f_aves[5]-f_aves[2]*f_aves[2])
 variable errpress equal sqrt(f_aves[6]-f_aves[3]*f_aves[3])
 variable heatcap equal ${no_atoms2}*((f_aves[8]-(f_aves[7]*f_aves[7]))/f_aves[5])
    
 print "Averages"
 print "--------"
 print "Density: ${avedens}"
 print "Stderr: ${errdens}"
 print "Temperature: ${avetemp}"
 print "Stderr: ${errtemp}"
 print "Pressure: ${avepress}"
 print "Stderr: ${errpress}"
 print "heatcap: ${heatcap}"

It is expected that heat capacity will increase with density because there is more molecules in the same volume so more energy can be stored. It is also expected that the heat capacity decreases with temperature because at higher temperature, higher energy levels starts getting occupied therefore the system won't be able to store as much energy.

Structural properties and the radial distribution function

TASK: 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?

Solid phase's RDF fluctuates because it has an ordered structure so the atoms remain in roughly the same place relative to a certain origin and form distinct shells. The peaks refers to the location of lattice points in a face-centered cubic lattice. The RDF decreases with distance because the number of atoms in the shell increases at a slower rate than the increase in volume.

Liquid phase's RDF fluctuates less than solids because it has a less ordered structure and atoms are able to move around therefore the shells are less distinct and more spaced out.

In the vapour phase, atoms are not ordered and far apart therefore only one initial peak is shown because LJ pair potential only has a significant effect on small internuclear distances. At higher internuclear distances, RDF=1 because of the random movement of atoms.

The first peak corresponds to atom A, with a lattice spacing of 1.025 and a coordination number of 12. The second peak corresponds to atom B, with a lattice spacing of 1.475 and a coordination number of 18. The third peak corresponds to atom C, with a lattice spacing of 1.825 and a coordination number of 42.

Dynamical properties and the diffusion coefficient

TASK: 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.

TASK: 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.

These are for 8000 atoms:

MSD is the average distance moved by an atom therefore, it is expected that the vapour phase will have the highest MSD because atoms are able to move freely. Solid phase has the lowest MSD because atoms are fixed in a lattice therefore there won't be much movement except for thermal vibrations. Liquid state has MSD values that are higher than but close to solids because it is able to move around but has similar density than solids.

The increase in MSD is parabolic initially for vapour phase because there is no collision due to its low density but it turns linear when collisions start to happen. Liquid phase's MSD increases linearly because it has similar density to solids so there is a lot of collisions as it travels.

These are for 1 million atoms:

The simulation with 1 million atoms has less fluctuation because the MSD are averaged over a lot more atoms, resulting in a more accurate representation.

Number of Atoms	Solid	Liquid	Vapour
8000	9.57x10^-9	1.33x10^-4	8.07x10^-3
1 million	8.78x10^-9	1.67x10^-4	5.08x10^-3

Velocity Autocorrelation Function

TASK: 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 (τ) = \frac{\int_{- \infty}^{\infty} v (t) v (t + τ) d t}{\int_{- \infty}^{\infty} v^{2} (t) d t}$

Be sure to show your working in your writeup. On the same graph, with x range 0 to 500, plot $C (τ)$ with $ω = 1 / 2 π$ 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 position of a classical harmonic oscillator is given by:

$x (t) = A c o s (ω t + ϕ)$ .

By differentiating $x (t)$ with respect to time, it gives the velocity:

$v (t) = \frac{d x (t)}{d t} = - A ω s i n (ω t + ϕ)$ .

The velocity autocorrelation function is given by:

$C (τ) = ⟨ v (t) \cdot v (t + τ) ⟩$

The normalised function is given by:

$C (τ) = \frac{\int_{- \infty}^{\infty} v (t) v (t + τ) d t}{\int_{- \infty}^{\infty} v^{2} (t) d t}$

By substituting the velocity into the normalised function, it gives:

$C (τ) = \frac{\int_{- \infty}^{\infty} s i n (ω t + ϕ) s i n (ω (t + τ) + ϕ) d t}{\int_{- \infty}^{\infty} (s i n (ω t + ϕ))^{2} d t}$

Using the sine addition formula for $t + τ$ , it gives:

$C (τ) = \frac{\int_{- \infty}^{\infty} s i n (ω t + ϕ) (s i n (ω t + ϕ) c o s (ω τ) + s i n (ω τ) c o s (ω t + ϕ)) d t}{\int_{- \infty}^{\infty} (s i n (ω t + ϕ))^{2} d t} = c o s (ω τ) + \frac{\int_{- \infty}^{\infty} s i n (ω t + ϕ) s i n (ω τ) c o s (ω t + ϕ) d t}{\int_{- \infty}^{\infty} (s i n (ω t + ϕ))^{2} d t}$

$\frac{\int_{- \infty}^{\infty} s i n (ω t + ϕ) s i n (ω τ) c o s (ω t + ϕ) d t}{\int_{- \infty}^{\infty} (s i n (ω t + ϕ))^{2} d t} = 0$

Therefore:

$C (τ) = c o s (ω τ)$

The minima in the VACF plot represents the collision with other atoms. In solid phase, the atoms can't move around and can only vibrate because there will be a restoring force, as shown by the fluctuations. However, in liquid phase, atoms can move around so there won't be much restoring force to move it back in place, as shown by the rapid damping of the velocity.

The harmonic oscillator VACF is different to the LJ solid and liquid because it doesn't take external forces from other atoms into account. It only considers an single atom that has a restoring force that is proportional to its displacement from the equilibrium position.

TASK: 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?

These are for 8000 atoms:

These are for 1 million atoms:

Number of Atoms	Solid	Liquid	Vapour
8000	1.57x10^-4	8.19x10^-2	6.82
1 million	4.79x10^-5	9.27x10^-2	3.49

These values are as expected because $D$ is greatest for vapour phase then liquid and solid for the same reasons as above.

The largest source of error in the estimate of $D$ for VACF will be the use of the trapezium rule, therefore some areas will not be accounted for, leading to a smaller than expected value. The accuracy can be increased by using smaller trapeziums.