Rep:Mod:nhl13

Introduction to Molecular Dynamics Simulation

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

Fig.1 Analytical position vs. Time.

Fig.2 Total Energy from velocity-Verlet solution vs. Time.

Fig.3 Error vs. Time.

Analytical position of a classical harmonic oscillator was plotted in figure 1, where the analytical position was calculated by equation:

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

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

The total energy of the classical harmonic oscillator for the velocity-Verlet solution was calculated by equation:

$TotalEnergy=KineticEnergy+PotentialEnergy$,

For a classical harmonic oscillator, $$\begin{gathered} KE=\frac{1}{2} \\ m{v}^2 \end{gathered}$$ and $$\begin{gathered} PE=\frac{1}{2} \\ k{x}^2 \end{gathered}$$,

and x(t) from velocity-Verlet solution was used instead of analytical x(t).

The total energy of the oscillator was plotted in figure 2. The energy fluctuated between 0.500 and 0.49875, which was relatively small when compared to the maximum total energy, i.e. 0.25% of the total energy.

The error in position between classical solution and velocity-Verlet solution was plotted in figure 3. It showed that the error gradually increased with time.

TASK 2: 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 showed the maxima of error in position between the classical and velocity-Verlet solution respectively was plotted against time. It gave a straight line with a positive slope that indicated the error was directly proportional to time. The corresponding function was y = 0.0004x, when the line was set to intercept at origin.

Fig 4.Position of error maxima vs. Time.

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

Fig.5 Fluctuation in energy with different timestep.

The maximum total energy of the harmonic oscillator was 0.500. From figure 5, the degree of fluctuation increased when time step increased. The total energy had a fluctuation of 5x10^-3, i.e. 1% of the maximum energy, when the timestep reached 0.20 s. Hence, the timestep had to be equal or smaller than 0.20 s in order to obtain an accurate simulation.

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

When potential energy $\varphi =0$, $r=r_0$

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

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

$(\frac{{\sigma }^{12}}{r_0^{12}}=\frac{{\sigma }^6}{r_0^6})$

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

$r_0=\sigma$

TASK 4.2: What is the force at this separation?

F (r) = - dU/dr

F (r) = - d [4 ε (σ¹² / r¹²) + 4 ε (σ⁶ / r₀⁶)] / dr

F (r) = d [ - 4 ε (σ¹² / r¹²) + 4 ε (σ⁶ / r₀⁶)] / dr

F (r) = d [ - 4 ε (σ¹² r^-12) + 4 ε (σ⁶ r^-6)] / dr

F (r) = 48 ε σ¹² r^-13 - 24 ε σ⁶ r^-7

When r = r₀, r₀ = σ

F (r₀) = 48 ε σ¹² σ^-13 - 24 ε σ⁶ σ^-7

F (r₀) = 48 ε σ^-1 - 24 ε σ^-1

F (r₀) = 24 ε / σ

TASK 4.3: Find the equilibrium separation, $r_{eq}$, and work out the well depth ($\varphi \left(r_{eq}\right)$).

When F = 0, r = r_eq

F (r) = 48 ε σ¹² r^-13 - 24 ε σ⁶ r^-7

0 = 48 ε σ¹² r_eq^-13 - 24 ε σ⁶ r_eq^-7

48 ε σ¹² r_eq^-13 = 24 ε σ⁶ r_eq^-7

2 σ¹² / σ⁶ = r_eq¹³ / r_eq⁷

2 σ⁶ = r_eq⁶

r_eq = 2^1/6 σ

When r_eq = 2^1/6 σ,

Φ = 4 ε [(σ¹² /r¹² ) - (σ⁶ / r⁶)]

Φ = 4 ε {[σ¹² /(2^1/6 σ)¹²] - [σ⁶ / (2^1/6 σ)⁶]}

Φ = 4 ε [(1/4) - (1/2)]

Φ = - ε

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

When $\sigma =ϵ=1.0$

= ʃ 4 ε [(σ¹² / r¹² ) - (σ⁶ / r⁶)] dr

= [-4 / 11 r¹¹ + 4 / 5 r⁵]

Therefore,

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

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

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

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

Density of water = 1 g cm^-3

Mass of water in 1 ml of water = 1 g

No. of moles in 1 ml of water = 1 g ÷ 18.0 g mol^-1 = 0.0556 mol

No. of water molecules in 1 ml of water = 0.0556 mol × 6.022×10²³ = 3.34×10²² molecules

No. of moles of 10000 water molecules = 10000 ÷ 6.022×10²² = 1.66×10^-20 mol

Mass of 10000 water molecules = 1.66×10^-20 mol × 18.0 g mol^-1 = 2.99×10^-19 g

volume of 10000 water molecules = 2.99×10^-19 g ÷ 1 g cm^-1 = 2.99×10^-19 ml

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

Periodic boundary condition states that boxes are repeated infinitely in all direction, hence, if one molecule exits the boundary of the original box, one molecule from the replicated box enters the original box through the opposite face to retain the number of molecules. As a result, we could work out the final position of the atom simply by adding the movement vector to its initial position vector, i.e. (0.5, 0.5, 0.5) + (0.7, 0.6, 0.2) = (1.2, 1.1, 0.7). As the box had a dimension of (1.0, 1.0, 1.0), so the resulting position vector was (0.2, 0.1, 0.7).

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

r* = r / σ

r = r* σ

r = (3.2)(0.34×10^-9)

r = 1.088 nm

TASK 7.2: What is the well depth in $$\begin{gathered} {\mathrm{k}\mathrm{J} \\ \mathrm{m}\mathrm{o}\mathrm{l}}^{-1} \end{gathered}$$?

Φ = - ε, as calculated above, and ε / k_B = 120K

Φ = - 120 × k_B

Φ = - 1.66×10^-21 J ÷ 1000 × N_A

Φ = 0.998 kJmol^-1

TASK 7.3: What is the reduced temperature $T^{*}=1.5$ in real units?

T* = k_BT ÷ ε

T = (T*)(ε) ÷ k_B

T = (1.5)(120)(k_B) ÷ k_B

T = (1.5)(120)

T = 180 K

Equilibration

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

In liquid, molecules are held together by intermolecular bonds. They are not fixed and can flow around and have a certain average intermolecular distance. If the molecules are too far apart, the attractive interaction between them is not as strong as it used to and behaves as gas instead of liquid. In contrast, if molecules are too close together, they experience a repulsion force between them. The liquid simulation would be affected in both extremes.

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

$numberdensity=\frac{numberofatom}{volume}$ can be used to rationalise the relationship between the number density, and the lattice spacing, a.

Case 1: In a simple cubic lattice, there is only 1 atom inside a lattice unit cell.

$numberdensity=\frac{numberofatom}{volume}$

$0.8=\frac{0.8}{a^3}$

$a=1.0772m$

This result agreed with the value obtained from LAMMPS simulation, i.e. lattice spacing in x,y,z = 1.07722 1.07722 1.07722

Case 2: In a face-centred cubic lattice, there are 1 atom at each corner and 1 atom on each face. Atom at the corner and on the face contribute to $\frac{1}{4}$ and $\frac{1}{2}$ atom to the lattice. Hence, number of atom in a fcc lattice = $\frac{1}{8}\times 8+\frac{1}{2}\times 6=4$

$numberdensity=\frac{numberofatom}{volume}$

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

$a=1.494m$

Therefore, a fcc lattice with a number density of 1.2 will have a side length of 1.494 m.

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

Each fcc unit cell consists of 4 atoms. With the same 10×10×10 box, i.e. 1000 lattice unit cells, hence, fcc lattice will consist of 4 × 1000 = 4000 atoms.

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

mass 1 1.0

Syntax of the command

mass I value

Mass command defines the mass of atoms, where I is the number of atom type and value is the mass of atom. In our case, there is only 1 type of atom and its mass is 1.0 in reduced unit.

pair_style lj/cut 3.0

Syntax of the command

pair_style style args

pair_style command compute different styles of pairwise interaction. The second ''style'' defines a particular type of pairwise interaction, and the ''args'' is the arguments used by that particular style. In our case, lj/cut is used to compute the cutoff Lennard-Jones potential with no Coulomb, and the cutoff potential is set at 3.0 reduced unit.

pair_coeff * * 1.0 1.0

Syntax of the command

pair_coeff command specifies the pairwise force field coefficients for one or more pairs of atom types. ''I,J'' is the atom types and ''args'' is the coefficients for one or more pairs of atom types. In out case, I and J equal * which means there is N types of atom and means all types from 1 to N, and both coefficients are set at 1.0.

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

We are using Velocity Verlet Algorithm, see below.

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

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

Those extra lines would allow the variable to change automatically when timestep in the second line is changed.

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

Fig6. Equilibration of total energy.

Fig.7 Equilibration of Temperature.

Fig.8 Equilibration of Pressure.

Fig.9 Total energy at different timestep.

All three properties reached equilibrium after short period of time when timestep was 0.001 s. The total energy (figure 6), temperature (figure 7) and pressure (figure 8) reached the equilibrium after approximately 0.42, 0.32 and 0.25 s respectively. Figure 9 showed the total energy fluctuation against time simulated at different timestep. When timestep was as large as 0.015 s, the total energy did not reach equilibrium. This timestep was not suitable for simulation. When timestep was equal to or lower than 0.01 s, the energy did reach equilibrium. However, the fluctuation of the energy level was quite large when timestep equal 0.01 and 0.075 s. On the other hand, timestep cannot be too short, as it lowers the efficiency of the simulation. As a result, the most acceptable timestep for our simulation to compensate both disadvantages was 0.0025 s.

Running simulations under specific conditions

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

In our simulation, timestep of 0.0025 s was used, as it allowed the oscillator to reach equilibrium and had a small fluctuation in energy level. The two pressures used were 2.55 and 2.65 atm, as there were within the range of the pressure equilibrium calculated in the last simulation. Finally, the temperature used were 1.7, 1.9, 2.1, 2.3 and 2.5, where all of them were above the critical temperature.

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

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

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

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

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

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

$\gamma =(\frac{\mathfrak{𝔗}}{T}{)}^{1/2}$

TASK 17: 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 three number in the command corresponded to Nevery, Nrepeat and Nfreq respectively. Nevery tells how often are the input values to be used for calculation, Nrepeat tells the number of times of each set of input values is used to calculate the average, and Nfreq tells the number of timestep for calculating an average. In our case, i.e. 100 1000 100000, the input values are used every 100 timesteps, the average calculated from input values are repeated for 1000 times, and the averages are calculated every 100000 timesteps. As timestep of 0.0025 s was used, our simulation time would be 100000 x 0.0025 = 250 s

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

Fig.10 Density vs. Temp at different pressure.

Figure 10 showed that the simulated density decreases with temperature, while increases with pressure. These phenomenons could be rationalised easily. As temperature increases, molecules have higher kinetic energy to overcome the attractive interaction between them, and move apart easier, thus resulting a lower density. On the other hand, when pressure increases, molecules experience a greater compression force and are packed closer together, thus resulting a higher density. Figure 10 also showed that the error in density was much smaller than the error in temperature. On the other hand, the simulated density was much lower than the density calculated from the ideal gas law. It was because the ideal gas law assume there were no repulsive forces between molecules, while the simulation by LAMMPS was able to mimic the repulsive forces between molecules. As consequence, the simulated density was lower as repulsive force between molecules pushed each other away. Moreover, the discrepancy between simulated and calculated density increased when pressure increased. It was because gas would behave more as a real gas rather than ideal gas under higher pressure, hence, the discrepancy increased as ideal gas was no longer relevant to the situation.

Calculating heat capacities using statistical physics

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

Fig.11 heat capacity against temperature at different density

Cv/V was plotted against temperature in figure 11. For both density, the simulation indicated the volume of the liquid remain unchanged, i.e. 16875 and 4218.75. The graph showed that the Cv/V decrease with temperature. The extent of decrease was greater when temperature was small and started to decrease slower as temperature increase gradually. However, there was an anomalous when T=2.4, D=0.8, which possibly arose due Brownian motion of the liquid.

Example of the script is shown below, where density is 0.2 and temperature is 2.0.

### DEFINE SIMULATION BOX GEOMETRY ###
lattice sc 0.2
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
unfix nvt
fix nve all nve

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp press density atoms vol
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 N equal atoms
variable N2 equal atoms*atoms
variable etotal equal etotal
variable etotal2 equal etotal*etotal
variable volume equal vol

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

variable avedens equal f_aves[1]
variable avetemp equal f_aves[2]
variable avepress equal f_aves[3]
variable aveetotal equal f_aves[4]
variable aveetotal2 equal f_aves[8]

variable Cv equal ${N2}*(${aveetotal2}-${aveetotal}*${aveetotal})/(${avetemp}*${avetemp})

print "Averages"
print "--------"
print "Density: ${avedens}"
print "Temperature: ${avetemp}"
print "Pressure: ${avepress}" 
print "Heat Capacity: ${Cv}"
print "Total Energy: ${aveetotal}"
print "Total Energy2: ${aveetotal2}"
print "Volume: ${volume}"
print "No of atoms square: ${N2}"

Structural properties and the radial distribution function

TASK 20: 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 following table showed the conditions under which the simulations of the Lennard-Jones system were performed.

	Temperature	Density
Solid	1.25	1.25
Liquid	1.1	0.7
Gas	1.1	0.05

Fig.12 RDF for solid, liquid and gas.

Fig.13 Integral of RDF for solid, liquid and gas.

RDF provides the density of atoms at particular radial distance to a reference atom. Each peak in RDF indicates a particularly favoured separation distance of a neighbouring atom to a reference atom. Thus, RDF reveals the atomic structure of the system being simulated. Figure 12 showed that RDF of solid had lots of sharp peaks, while there were only a few broad peaks for liquid, and just one large broad peak for gas. It revealed atoms in solid were highly ordered and stayed at particular distance to the reference atom. Sharp peaks also indicated the atoms were not spread out randomly but were concentrated at several favoured particular distances. Moreover, the coordination number at favoured particular distances could be determined from the integral of RDF. In contrast, a few peaks in liquid RDF revealed its system was less ordered and atoms were able to move around freely and rarely stay at a particular distance to the reference atom. However, peaks still showed up illustrated that the atoms were still held together by some sort of attraction forces. On the other hand, broad peaks were observed in RDF of liquid and gas, while peak for gas was broader than those for liquid. This could be rationalised by the Brownian motion of the atoms in the system. The atoms in gas had higher kinetic energy, therefore, they could move faster and more randomly and resulted in a broad distribution at particular distance, where the only peak observed indicated the gas atoms were still holding with each other at a maximum interacting distance.

Fig.14 1st and 2nd coordinating neighbours of fcc.

Fig.15 3rd coordinating neighbours of fcc.

The first three peaks at the RDF of soild appeared at 1.025, 1.475 and 1.825. These corresponded to the first three sets of nearest coordinating neighbours to a reference atom. The simulation shown the lattice spacing, a, was 1.47361, and by visualising the fcc lattice model, we knew the nearest neighbour had a separation of $\frac{\sqrt{2}}{2}a$, the second nearest neighbour had a separation of $a$, and the third nearest neighbour had a separation of $\sqrt{\frac{3}{2}}a$. When we combined all informations, the lattice spacing between the reference atom and the nearest neighbours could then be calculated by Pythagoras theorem. The results were summarised in the table below.It was found that the peaks in RDF were similar and comparable to the calculated results. On the other hand, the integral of RDF, which was accumulative, showed the coordination number of the neighbouring atoms. This also agreed with the theoretical fcc model. See figure 14 and 15. To conclude, both peak and integral of RDF gave evidence that RDF could be used to determine lattice structure.

no of peak	RDF Peak, Lattice spacing	Calculated lattice spacing	Coordination Number	Integral of RDF
First	1.025	1.042	12	11.9$\approx$12
Second	1.475	1.474	6	17.8-11.9$\approx$6
Third	1.825	1. 805	24	41.8-17.8$\approx$24

Dynamical properties and the diffusion coefficient

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

The following table shows the conditions under which the simulations were performed.

	Temperature	Density
Solid	1.25	1.25
Liquid	1.1	0.7
Gas	1.1	0.05

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

Fig.16 MSD of gas against timestep.

Fig.17 MSD of liquid against timestep.

Fig.18 MSD of solid against timestep.

From equation $D=\frac{1}{6}\frac{\partial ⟨r^2\left(t\right)⟩}{\partial t}$, we can calculate D in reduced unit, D*, from the MSD-timestep graph, i.e. $D^{*}=\frac{1}{6}\times gradient$. The timestep of the simulation was 0.002 s, therefore, $D=\frac{D^{*}}{0.002}$. The results are illustrated in the table below.

D from MSD	Solid	Liquid	Gas
Simulated Data	$4.28*10^{-6}$	$0.121$	$2.92$
1 Million Atom Data	$4.39*10^{-6}$	$0.0873$	$3.00$

The D values obtained were expected. Gas had the highest D among the three, as gas molecules had the highest kinetic energy, and can move around freely with only a little restriction. Liquid had a smaller D value, while liquid molecules can move around freely, but have lower kinetic energy, pack closer together and have less space to move around. In contrast, solid had a very low D value of $10^{-6}$. It is because solid molecules have the least kinetic energy. Also, they are ordered in a rigid shape, and can only oscillate about its equilibrium position. On the other hand, Simulated data was similar to the 1 million atom data given, which showed consistency of the the simulation.

TASK 23.1: 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}$

The position of the oscillator is defined as $x\left(t\right)=A\cos (\omega t+\varphi )$

The velocity is the first derivative of the position equation: $v=\frac{\mathrm{d}\mathbf{x}\left(t\right)}{\mathrm{d}t}=-A\omega \sin (\omega t+\varphi )$

Then substitute $v=-A\omega \sin (\omega t+\varphi )$ into the $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}$

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

And$-A^2{\omega }^2$ can be cancelled out.

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

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

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

Using trigonometry identity: $\sin (A+B)=\sin (A)\cos (B)+\cos (A)\sin (B)$

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

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

$C(\tau )=\frac{\cos (\omega \tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^{}(\omega t+\varphi )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^2(\omega t+\varphi )\mathrm{d}t}+\frac{\sin (\omega \tau )\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}$

$C(\tau )=\cos (\omega \tau )+\frac{\sin (\omega \tau )\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}$

We know that Sine is an odd function and cosine is an even function. We also know an odd function multiplies an even function gives an odd function, and the integral of an odd function with a symmetrical limit becomes 0.

Therefore, $\frac{\sin (\omega \tau )\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}=0$

Therefore, $C(\tau )=\cos (\omega \tau )$.

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

Figure 19 compared the VACF of a solid and a liquid simulation to the normalised classic harmonic oscillator. The solid simulation showed a minima at timestep = 47, while liquid simulation did not give a obvious minima, its function fluctuated about 0 after timestep = 109. Minima was observed because there were atom collisions in the solid and liquid simulations, and caused a change in velocity direction. Moreover, collisions in Lennard-Jones system led to loss of velocity and loss of energy, hence, the VACF was no longer correlated to time at very long time and remained at 0. In contrast, the VACF of harmonic oscillator gave a periodic function because there was no atom collision within the system. Therefore, no lost of energy and the oscillator would oscillate forever.

Fig.19 VACF of solid, liquid and harmonic oscillator against timestep.

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

Diffusion coefficient could be calculated from equation: $D=\frac{1}{3}\underset{0}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}\tau ⟨\mathbf{v}\left(0\right)\cdot \mathbf{v}\left(\tau \right)⟩$

Hence, $D*=\frac{1}{3}\times integral$

$D=D*\times 0.002$

Fig.20 Integral of VACF of gas against timestep.

The results were compared with the D calculated from MSD and were showed in the table below.

D from VACF	Solid	Liquid	Gas	D from MSD	Solid	Liquid	Gas
Simulated Data	$5.51*10^{-4}$	$0.113$	$0.885$	Simulated Data	$4.28*10^{-4}$	$0.121$	$2.92$
1 Million Atom Data	$2.78*10^{-4}$	$0.0825$	$0.978$	1 Million Atom Data	$4.39*10^{-4}$	$0.0873$	$3.00$

Fig.21 Integral of VACF of liquid against timestep.

The diffusion coefficient, D, obtained from MSD and VACF were not quite the same, however, they did have comparable value in term of the magnitude. D of solid was almost zero for both calculations. D calculated from VACF was 100 times larger than that from MSD. In contrast, two calculated D values for liquid were similar and were comparable, even for both 1 atom simulation and 1 million atom simulation. The D value obtained for 1 million atom simulation were similar to those for 1 atom simulation. The largest source of error for estimate D from VACF was that it was only an approximation of the integral.

Fig.22 Integral of VACF of solid against timestep.

Fig.23 Integral of VACF of 1 million gas molecules against timestep.

Fig.24 Integral of VACF of 1 million liquid molecules against timestep.

Fig.25 Integral of VACF of 1 million solid against timestep.