Talk:Mod:APPLEBANANAKIWI

JC: General comments: Good report with detailed answers given to all tasks. Make sure that you understand the background theory behind the last few tasks so that you can fully explain what your results show.

Introduction

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\left(t\right)=A\cos (\omega t+\varphi )$ (the values of $A$, $\omega$, and $\varphi$ are worked out for you in the sheet).

Below in Figure 1 is the graph for the displacement $x\left(t\right)$ calculated using firstly the Verlet algorithm and secondly the analytical solution using the classical harmonic oscillator. Latter reduces down to $x\left(t\right)=\cos \left(t\right)$ as all the parameters are set to either zero or one. From the direct comparison on this relatively low resolution graph no significant deviation is observed. In Figure 2 the total energy versus time is shown at a time step of 0.1. The total energy was calculated as a sum of kinetic and potential energies and given in reduced form (as both the force constant k as well as the mass equals to one): $E_T=\frac{1}{2}(v^2+x^2)$. This expression originates from the kinetic energy expression $E_K=\frac{1}{2}m{v}^2$ and the potential energy expression $E_P=\frac{1}{2}k{x}^2$ which represents the potential energy gained as a function of displacement from the equilibrium distance.

Figure 1: Displacement

Figure 2: Total Energy

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.

The error between the Verlet algorithm and the classical harmonic oscillator seems to increase with time according to the equation in the diagram. It has to be noted that the error has a local maxima every time the displacement reaches zero and a local minima every time the displacement is at a turning point.

Figure 3: Error

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 wide range of time steps were used ranging from $t=0.01$ up to $t=0.5$. The table below represent the maximum error found in the simulations. As expected the error increased with higher time steps while it decreased with smaller time steps. The threshold of decreasing the error below 1% (0.01 in absolute values) of the displacement seemed to be just around $t=0.2$. It has to be noted that the percentage error of displacement will be proportionate to the percentage error in total energy as former is used to calculate the latter. Any calculation with a time step lower than 0.2 would therefore yield a sufficiently accurate result. However, it has to be noted that by decreasing the time step the length of the simulation decreases as well. In order to evaluate the same amount of steps one would have to run the calculation longer. It essentially comes down to a trade-off between computational duration and accuracy of results.

Timestep	Error in Displacement (absolute values)
0.01	0.000006
0.05	0.0005
0.075	0.002
0.1	0.005
0.2	0.09
0.5	0.7

JC: Good choice of timestep, but why is it important to check that total energy doesn't fluctuate too much?

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

Calculation of $r_0$	Comments
$\varphi \left(r\right)=4ϵ(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})=0$	For a single interaction the Lennard-Jones potential is shown. It consists of the combination of the repulsive part which dominates at small distances and the attractive term which dominates at larger distances. Here it is set to zero to find $r_0$.
$(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})=0$	In order for the function to equal zero either the inside of the brackets has to equal to zero (or the constant $ϵ=0$).
$\frac{{\sigma }^6}{r^6}(\frac{{\sigma }^6}{r^6}-1)=0$	Factorising and realising that the expression inside the brackets has to equal to zero (or the constant $\sigma =0$).
$\frac{{\sigma }^6}{r^6}=1$	rearranging and taking the sixth root leads to the next line
$r_0=\sigma$	Alternatively, when $r_0$ approaches infinity would also yield the potential to approach zero.

Calculation of ${\mathbf{F}}_i$	Comments
${\mathbf{F}}_i=-\frac{\mathrm{d}U\left({\mathbf{r}}^N\right)}{\mathrm{d}{\mathbf{r}}_i}$	The force at this separation is given by the derivative of the potential with respect to the distance.
${\mathbf{F}}_i=4ϵ(\frac{12{\sigma }^{12}}{r^{13}}-\frac{6{\sigma }^6}{r^7})$	This is the expression for that one interaction described in the previous example. Setting $r=\sigma$ and rearranging leads to the next line.
${\mathbf{F}}_i=\frac{24ϵ}{\sigma }$	The force thus depends on both the depth of the potential well as well as the distance $\sigma$ (which is essentially $r_0$).

Calculation of $r_{eq}$	Comments
$-{\mathbf{F}}_i=-4ϵ(\frac{12{\sigma }^{12}}{r^{13}}-\frac{6{\sigma }^6}{r^7})=0$	In order to find the equilibrium distance the derivative of the potential (negative force) has to be set to zero in order to find that minima. Using the same reasoning as in the above example, the expression inside the brackets should equal to zero.
$\frac{6{\sigma }^6}{r^7}(\frac{2{\sigma }^6}{r^6}-1)=0$	Using the expression of the inside of the brackets and setting it to zero again leads to the final line.
$r_{eq}=2^{\frac{1}{6}}\sigma$	The equilibrium distance is proportional to the parameter distance $\sigma$.

Calculation of the depth of the potential well	Comments
$U=4ϵ(\frac{{\sigma }^{12}}{4{\sigma }^{12}}-\frac{{\sigma }^6}{2{\sigma }^6})$	Substituting the value for $r_{eq}$ obtained above into the expression of the potential.
$U=4ϵ(\frac{1}{4}-\frac{1}{2})$	Simplifying the fractions leads to the final line.
$U=-ϵ$	The depth of the potential well is proportional to $ϵ$, which makes sense as this is the parameter which controls this property.

Evaluation of the total interaction integrals	Comments
$\underset{2\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=\underset{2\sigma }{\overset{\mathrm{\infty }}{\int }}4ϵ(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})\mathrm{d}r=\underset{\alpha }{\overset{\mathrm{\infty }}{\int }}4(\frac{1}{r^{12}}-\frac{1}{r^6})\mathrm{d}r$	Setting the parameters to one (as instructed) and changing the lower limit of the integral to $\alpha$ allows a general approach.
$\underset{\alpha }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=-4(\frac{1}{11r^{11}}-\frac{1}{5r^5}){|}_{\alpha }^{\mathrm{\infty }}$	This is a standard integration. For the upper bound the integral evaluates to zero, while the lower bound is shown in the next line.
$\underset{\alpha }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=4(\frac{1}{11{\alpha }^{11}}-\frac{1}{5{\alpha }^5})$	This expression can now be evaluated for different values of $\alpha$ as shown below.
$\alpha =2$ gives an integral of $-0.0248$, $\alpha =2.5$ gives $-0.00818$, $\alpha =3$ gives $-0.00329$	The integral will become smaller with increasing lower limit and eventually tend to zero. In the physical world this corresponds to the two particles losing all interaction at infinite distance.

JC: All maths correct and well explained.

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

Estimation of the water molecule properties	Comments
$N{o}_{Molecules}=\frac{6.023\times 10^{23}}{18}=3.35\times 10^{22}$	1 ml of water corresponds to 1 g at room temperature and pressure. Considering the molar mass of water is approximately 18.015 g/mol, means that the inverse of this number is the amount of matter (in mole) in 1 g. In order to find the number of molecules that inverse has to be multiplied by Avogado's constant.
$V_{10000Molecules}=\frac{m}{\rho }=\frac{10000\times 18}{6.023\times 10^{23}}=2.99\times 10^{-19}mL$	The volume is related to the mass through the density (equal to one at standard conditions). The mass of 10000 molecules obtained by dividing the molar mass by Avogado's constant to obtain the mass of one molecule and then multiplying it with 10000.

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

Calculation of a position Vector	Comments
$\overrightarrow{v}=(0.5,0.5,0.5)+(0.7,0.6,0.2)=(1.2,1.1,0.7)$	Simple vector addition before boundary conditions are applied.
$\overrightarrow{v}=(0.2,0.1,0.7)$	The boundary conditions impose a cutoff in all dimensions at one, which results in following position vector.

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

Calculation of LJ parameters in real units for argon	Comments
$r=\sigma r*=1.09nm$	LJ cutoff in real units
$ϵ=R\times 120K=0.998kJ.mo{l}^{-1}$	By multiplying with the gas constant one can immediately obtain the well depth.
$T=\frac{ϵT*}{k_B}=180K$	Temperature in real units.

JC: All calculations correct and working clearly shown.

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?

This could cause problems as the atoms could be initiated on top of each other, i.e. one atom would be occupying space which shouldn't be available. Likewise this would lead to a false representation of the interaction as it is unlikely that two atoms would ever get that close to each other (infinitely high potential).

JC: The high repulsive forces caused by overlapping particles would make the simulation unstable and could cause it to crash.

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?

Number density of lattice points	Comments
$\rho =\frac{n_{atoms}}{V}=\frac{1}{1.07722^3}=0.8$	The number density is given by the number of atoms per unit volume. For a simple cubic lattice with one atom and sides of equal length one can simplify the expression as shown. This leads to the number density of lattice points of 0.8.
$L={\left(\frac{n_{atoms}}{\rho }\right)}^{\frac{1}{3}}={\left(\frac{4}{1.2}\right)}^{\frac{1}{3}}=1.49$	Rearranging the equation above leads to a side length of 1.49. It has to be noted that the fcc-lattice contains 4 atoms instead of only 1.

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?

As 1000 unit cells are created with the command four times the amount of atoms are created. Hence 4000 atoms.

JC: Correct.

TASK: 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 command refers to the mass assignment of the atoms. The first parameter defines the type of atom (which in our case is always one) and the second command assigns that type a mass (in our case 1.0).

The second command defines the pairwise interaction for the atoms of the system (given by the lj). The cut 3.0 corresponds to the cutoff at a distance of 3.0 (reduced units) as the relevant part of the interaction happens in the region below the cutoff (as shown with the integral calculation in the previous section). This is done to save some computational time.

The third command defines the coefficients of the LJ potential $\sigma$ and $ϵ$ (sets them equal to 1 here). The two stars correspond to the atoms this constrain applies on (stars symbolise it applying to all atoms).

JC: Good explanations.

TASK: 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 given the specified starting parameters.

TASK: Look at the lines below. 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?

The reason the longer code is used is that it automatically calculates and changes the number of steps the simulation has to run when a time step is given as input. In the short code this would have to be done manually. Additionally the use of named variables makes the code more user friendly, easier to follow and easier for subsequent parameter changes.

JC: Good.

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 first three graphs are the plots for pressure, temperature and energy against time for the 0.001 time step in the time region of 0 to 5. This showed to be sufficient even though the simulation run further until 100. However the noise remained constant around the plateaued trend line and was cut off to identify the time it takes for the simulation to reach equilibrium. This was observed to be just below 0.5 and thus happens very fast. After that the values fluctuate around a constant average value for the duration of the simulation.

Figure 4: Pressure

Figure 5: Temperature

Figure 6: Total Energy

The time step that showed the biggest deviation is unsurprisingly the largest time step of 0.015. Both the two lowest time steps 0.001 and 0.0025 seem to have the same average total energy and seem to give similar data quality. Thus it would be sensible to chose the 0.0025 to still give acceptable results as it has a lower computational cost compared to a time step of 0.001. The reason for not choosing a larger time step is the possibility that the system diverges or doesn't give an accurate representation of the system.

JC: Good choice of timestep, the average total energy shouldn't depend on the timestep so 0.0075 and 0.01 are not suitable.

Figure 7: Total Energy of all time step systems

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.

A Temperature range of 1.5, 1.7, 1.9, 2.1 and 2.3 was chosen as well as the pressures 2.4 and 2.6.

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

$E_K=\frac{3}{2}N{k}_{B}T$

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

Solve these to determine $\gamma$.

Finding and expression for $\gamma$	Comments
$\frac{1}{2}\sum m_i{v}_i^2=\frac{3}{2}N{k}_{B}T$$\frac{1}{2}\sum m_i(\gamma v_i{)}^2=\frac{3}{2}N{k}_B\mathfrak{𝔗}$	The second equation originates from multiplying the velocity with $\gamma$ resulting in a new Temperature $\mathfrak{𝔗}$. Dividing the upper equation by the lower equation leads to the next line. (Note that gamma can be taken outside the summation as it is a constant.)
$\frac{1}{{\gamma }^2}=\frac{T}{\mathfrak{𝔗}}$	Everything cancels apart from the displayed equation.
$\gamma =\pm \sqrt{\frac{\mathfrak{𝔗}}{T}}$	Hence shown.

JC: Correct.

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?

The 100 indicates that every 100th input is sampled and used to calculate the average. The 1000 indicates that there will be 1000 data points used to calculate the average. The 100000 corresponds to the amount of steps after which an average is calculated. So using a time step of 0.0025 and a run of 100000 results in a total simulated time of 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 $\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?

Below are the graphs for the densities at pressures 2.4 and 2.6 respectively. The densities seem to decrease with increasing temperatures at both pressures. This is most likely due to the the expansion of the material with increasing temperature that leads to the decrease in density. The values of the two graphs only differ slightly from each other and both simulated densities lie lower than the data obtained through the ideal gas law. This is due to the ideal gas law assuming no interactions between atoms and no excluded volume which allows particles to be in very very close proximity of each other without being repelled. As those factors are present in the simulation it forces the particles at larger distance thus decreasing the overall density in comparison to the ideal gas law. The discrepancy between the two methods increases with pressure as the repulsion increases with higher pressure and thus the gap between ideal gas and simulation increases as well.

Figure 8: Density against Temperature at a Pressure of 2.4

Figure 9: Density against Temperature at a Pressure of 2.6

JC: Correct explanation, show all data on a single graph to make it clear that the simulated and ideal gas results are closer at lower 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 $({\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.

Firstly, the heat capacity per volume decreases with increasing temperature. This has to the do with the density of states of the system. At higher temperature the higher energy states become more populated and thus less energy is required to obtain the same final state so that less heat is absorbed. Secondly the heat capacity per volume increases with increasing density. This has to do with the increase in atoms in the system at constant volume. As there are more atoms they require more energy to raise the temperature. Thirdly, according to thermodynamic predictions the heat capacity was expected to be a constant independent of temperature (i.e. 1.5 R).

Figure 10: Heat capacity per volume vs Temperature

JC: Good suggestions, more analysis would be needed to investigate why the heat capacity has this trend with temperature.

The input script is given below. The main features that were changed from the script of the previous sections are in the #MEASURE SYSTEM STATE# section. New variables were introduced for additional data to be collected. Additionally new properties like the heat capacity and number of atoms were added. The heat capacity per volume was calculated directly through the fact that $\frac{C_V}{V}=N\rho \frac{⟨E^2⟩-{⟨E⟩}^2}{k_B{T}^2}$. This evaluates the expression directly without having to calculate the volume of the system first.

### DEFINE SIMULATION BOX GEOMETRY ###
variable sysdens equal 0.2
lattice sc ${sysdens}
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
unfix nvt
fix nve all nve
reset_timestep 0

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp atoms density
variable n equal atoms
variable n2 equal atoms*atoms
variable totenergy equal etotal
variable totenergy2 equal etotal*etotal
variable temp equal temp
variable temp2 equal temp*temp
variable dens equal density
fix aves all ave/time 100 1000 100000 v_totenergy v_totenergy2 v_temp v_temp2 v_dens
run 100000

variable heatcap equal ${n2}*(f_aves[2]-f_aves[1]*f_aves[1])/f_aves[4]
variable avetemp equal f_aves[3]
variable avedens equal f_aves[5]
variable heatcappervol equal ${n}*f_aves[5]*(f_aves[2]-f_aves[1]*f_aves[1])/f_aves[4]

print "Averages"
print "--------"
print "Heat Capacity: ${heatcap}"
print "Heat Capacity per Volume: ${heatcappervol}"
print "Density: ${avedens}"
print "Temperature: ${avetemp}"
print "Atoms: ${n}"

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?

The radial distribution function gives the electron density as a function of r. It is used to evaluate the order of the system. The RDF for the solid is the most distinct and complex and shows well defined peaks at various distances. This indicates that electron density is particularly high at those distances, probably due to neighbouring particles fixed in space through the solid structure (most ordered system). Both the liquid and the vapour phase have a plateaued line (at 1 - the system density) which indicates that their neighbours are at random distances and not fixed (i.e. constantly in motion). However, they have an initial peak which looks very similar to the inverted LJ potential. The maxima indicates that the electron density is highest in the potential well of the interaction with the nearest neighbour. When looking closely three additional small maxima can be discovered for the liquid phase which must correspond to some kind of ordered structure in the close environment of the particle. This indicates that the liquid phase is still slightly more ordered than the vapour phase. It has also to be noted that the RDF is zero below a distance of one due to an unrealistic overlap of particles.

The first three peaks of the solid state correspond to the three closest neighbours in the fcc crystal structure, which are at distance $$\begin{gathered} \frac{\sqrt{2}}{2}a, \\ a \end{gathered}$$and $\frac{\sqrt{6}}{2}a$. The lattice parameter a can be evaluate from the density of 1.2 using the equation of a previous section and yields 1.49. Substituting for a gives to those three distances $r_1=1.06$, $r_2=1.49$ and $r_3=1.83$, which vaguely correspond to the maxima of the RDF graph (~2-3% error). The coordination number of those three peaks can be found using the integral graph and relating each distance to a corresponding integral value. For $r_1$ this equals rounded to 12, for $r_2$ this is 6 and for $r_3$ this is 24.

Figure 11: RDF of the three phases

Figure 12: Integral of the RDF of the three phases

JC: The radial distribution function is not limited to electron density, there is no electron density in your simulations, just classical particles. The solid phase has long range order whereas the liquid phase has sort range order, but no long range order. Good idea to express all of the first 3 peaks in terms of the lattice parameter, but would have been good to have used these expression to calculate the lattice parameter from the first 3 peaks and then compare this with the initial value, rather than the opposite way round.

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.

The diffusion coefficient was calculated using $D=\frac{1}{6}\frac{\partial ⟨r^2\left(t\right)⟩}{\partial t}$. The calculated diffusion coefficients from the respective gradients of the graph trendiness are summarised below for the whole section.

Diffusion Coefficient (Area/Time)	MSD	MSD (Million atoms)	VCAF	VCAF (Million atoms)
Gas	3.19	2.38	3.42	3.16
Liquid	0.0700	0.0708	0.0699	0.0896
Solid	6.00E-6	4.43E-6	1.79E-4	4.83E-5

Increasing the number of atoms for the simulation only deviated the diffusion coefficient slightly but could lead to more accurate results due to a higher resolution of data. Furthermore, a clear trend in the size of the diffusion coefficient is visible as it increases drastically between the phases of liquid and gases. This is due to sudden increase in degrees of freedom a gas has in comparison to a liquid or solid (most ordered). The diffusion coefficients for liquid is also quite a bit larger than for a solid for the same reason.

The MSD graphs showed unsurprisingly that the MSD for vapour increased almost exponentially due to its high mobility and freedom. In the close up the MSD of the solid phase increases very quickly to a low equilibrium value and stays there quite constrained due to its high ordered nature. The MSD of the liquid phase seems to increase slowly and linearly over time. Running the simulation with more atoms decreased the rate at which the MSD for the vapour phase increased. This is due to the fact that the atoms have more collision partners available and thus an extra "constraint". By increasing the number of atoms infinitely the slope will approximate first that of a liquid and later of a solid (as the density increases and increases).

JC: Show the lines of best fit on the graphs, did you fit to the entire data range or just to the linear part? The vapour phase MSD begins as a quadratic curve which represents ballistic motion, before collisions occur and the diffusive regime is reached.

Figure 13: MSD for the small simulation

Figure 14: MSD for the million atom simulation

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!):

Evaluating the normalised velocity autocorrelation function	Comments
$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}$	Given $x\left(t\right)=A\cos (\omega t+\varphi )$ the derivative with respect to time is $v\left(t\right)=-\omega A\sin (\omega t+\varphi )$ and can be substituted in
$C\left(\tau \right)=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\omega A\sin (\omega t+\varphi )\times \omega A\sin (\omega (t+\tau )+\varphi )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\omega }^2{A}^2{\sin }^2(\omega t+\varphi )\mathrm{d}t}$	Factorising the constants and using the trigonometric identity $\sin (\omega t+\omega \tau +\varphi )=\sin (\omega t+\varphi )\cos \left(\omega \tau \right)+\sin \left(\omega \tau \right)\cos (\omega t+\varphi )$.
$C\left(\tau \right)=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi )\times [\sin (\omega t+\varphi )\cos \left(\omega \tau \right)+\sin \left(\omega \tau \right)\cos (\omega t+\varphi )]\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^2(\omega t+\varphi )\mathrm{d}t}$	Factorising the square brackets and using the fact that t is independent of $\tau$ to take some components out of the integral.
$C\left(\tau \right)=\frac{\cos \left(\omega \tau \right)\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^{}(\omega t+\varphi )\mathrm{d}t+\sin \left(\omega \tau \right)\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi )\cos (\omega t+\varphi )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^2(\omega t+\varphi )\mathrm{d}t}$	The integrals in the first term cancel out.
$C\left(\tau \right)=\cos \left(\omega \tau \right)+\frac{\sin \left(\omega \tau \right)\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (\omega t+\varphi )\cos (\omega t+\varphi )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^2(\omega t+\varphi )\mathrm{d}t}$	The second term equals to zero due to nature of the parity of the functions under the integral. As it consists of a sine function (odd) and a cosine function (even) the product of both is always an odd function which has the property of havinthere is more time between collisions and so it takes longer forg a zero integral from -infinity to infinity.
$C\left(\tau \right)=\cos \left(\omega \tau \right)$	Hence shown.

JC: Good.

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 minima in the VACF are related to the minima in the RDF from the previous section. They have to do with the interactions of particles in their respective phase as every interaction influences the velocity of the particles. For the liquid there is a relatively high rate of collision which causes a rapid decrease in the oscillation and so the VACF goes to zero very rapidly. Contrary, a solid which has spacial constrains imposed will slowly find its uncorrelated state. This is seen by the relative large fluctuations over time in comparison to the liquid. The harmonic oscillator does not include any collisions so that the VACF for it will remain constant forever. The reason it is defined by a cosine function rather than a straight line is that its displacement around the equilibrium distance has harmonic behaviour (and thus all the properties related to it do too, i.e. velocity, force, etc.).

JC: Can you be a bit more specific, what happens when the VACF become negative?

Figure 15: Autocorrelation Function

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?

The diffusion coefficients are listed in the table of the above section and are slightly larger than the ones obtained using MSD, which was unexpected. However, the largest error using the VACF method originates from using the trapezium rule. Rather than taking the exact integral it is approximated by a sum of rectangles with finite width. In this report they were calculated using the outside triangles thus making the integral seem bigger than it is. This would be one explanation for the deviation to the MSD method. The diffusion coefficient trends described earlier are still observed.

Figure 16: VACF for the small simulation

Figure 17: VACF for the million atom simulation

JC: Change the y axis scale to show the full running integral for the vapour phase, calculating the diffusion coefficient from the VACF relies on the running integral reaching a plateau in the simulation time.

Conclusion

This lab used LAMMPS to model and calculate structural and (thermo)dynamic properties of atom ensembles. It investigated the influence of outside parameters (temperature, pressure) on the system and showed the importance of the time step vs data accuracy trade off. Further investigations into the RDF and the diffusion coefficient of the ensembles revealed some interesting and well-known properties about how the interaction and order of a structure is reflected in those parameters. It also demonstrated the use of different methods to calculate the diffusion coefficient and showed the importance of understanding the underlying methods in order to evaluate which results are more suitable.