Talk:Mod:mym113

JC: General comments: All tasks completed and report is well presented. Some of your explanations are quite confusing though so try to write more clearly and concisely and make sure that you understand the background to each task.

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

Following the classical harmonic oscillator equation above, the positions of atoms can be calculated (in the column "ANALYTICAL") by the equation below (Figure 1):

File:Position vs Time.PNG

Figure 1. Position vs Time

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

$A=1.00$

$\omega =1.00$

$\varphi =0.00$

"ERROR" is always positive as it is the absolute difference between "ANALYTICAL" and the velocity-Verlet solution.

"ENERGY" was calculated using the equation below (Figure 2):

File:Energy vs Time.PNG

Figure 2. Energy vs Time

$TotalEnergy=KineticEnergy+PotentialEnergy$

$$\begin{gathered} TotalEnergy=\frac{1}{2} \\ {m}_i{v_i}^2+\frac{1}{2} \\ {k}_i{x_i}^2 \end{gathered}$$

$m=1.00$

$k=1.00$

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 3 shows that error increases as time increases. From time to time, the error reaches a maximum point and fall back to zero. This repeats continuously. Maxima occurs when the simulated and velocity-Verlet solutions oscillate at the same time.

File:Error vs Time.PNG

Figure 3. Error vs Time

Maxima from the plot of "ERROR" against time (Figure 4) were found accurately using "Conditional Formatting" in Excel. The positions of the maxima are plotted as a function of time (Figure 4). The figure showed that the maximum of error increases as time increases with a linear relationship. The equation of the straight line was found to be $y=0.0004x-7\times 10^{-5}$.

JC: Good idea to use conditional formatting to find the maxima. Give your equation of best fit and R^2 value to more than 1 significant figure.

File:Maxima between classical and velocity-Verlet solution against time.PNG

Figure 4. Maxima between classical and velocity-Verlet solution against time

Table 1. Position of maxima and corresponding time

Time	Maxima
2.0	7.58E-004
4.9	2.01E-003
8.0	3.30E-003
11.1	4.60E-003
14.2	5.91E-003

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?

Different values of timestep was examined and the total energy was monitored. It was found that timestep less than or equal to 0.2 is needed to ensure the total energy does not change by more than 1% over the course of this simulation (Figure 5). As the timestep increases, the number of error increases and the magnitude of total energy changes significantly. Increasing the timestep increases the time interval between detection and therefore a larger difference in energy change was observed.

JC: Good choice of timestep. Increasing the timestep means that the simulation extrapolates further forward in time at each step and so the algorithm is less accurate, the time step doesn't just affect the interval between detections.

File:Energy vs time at different timestep.PNG

Figure 5. Energy vs time at different timestep

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.

At $r=r_0$, potential energy = 0

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

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

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

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

$r_0=\sigma$

TASK 4.2 What is the force at this separation?

${\mathbf{F}}_i=-\frac{\mathrm{d}U\left({\mathbf{r}}^N\right)}{\mathrm{d}{\mathbf{r}}_i}$

Differentiate the $\varphi \left(r_0\right)$ with respect to $r$:

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

To obtain the force, substitute $r_0=\sigma$:

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

$\mathbf{F}=-4ϵ(-\frac{12}{\sigma }+\frac{6}{\sigma })$

$\mathbf{F}=-4ϵ\left(\frac{-6}{\sigma }\right)$

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

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

$\mathbf{F}=-\frac{\mathrm{d}U}{\mathrm{d}\mathbf{r}}$

Equilibrium separation is when the force is at minimum.

$\mathbf{F}=-4ϵ(-\frac{12{\sigma }^{12}}{r_{eq}^{13}}+\frac{6{\sigma }^6}{r_{eq}^7})=0$

$\frac{-12{\sigma }^{12}}{r_{eq}^{13}}+\frac{6{\sigma }^6}{r_{eq}^7}=0$

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

$r_{eq}^6=2{\sigma }^6$

$r_{eq}=2^{1/6}\sigma$

The well depth ($\varphi \left(r_{eq}\right)$):

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

$\varphi \left(r_{eq}\right)=4ϵ(\frac{1}{4}-\frac{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$.

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

$\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.1 Estimate the number of water molecules in 1ml of water under standard conditions.

Density of water = $$\begin{gathered} 1 \\ g \\ c{m}^{-1} \end{gathered}$$

Mass of water in 1 mL of water = $$\begin{gathered} 1 \\ g \end{gathered}$$

No. of moles of water molecule in 1 mL of water = $$\begin{gathered} \frac{mass}{M_r}=\frac{1}{18}=0.0556 \\ mol \end{gathered}$$

No. of water molecules in 1 mL of water = $mol\times N_A=3.34\times 10^{22}$

TASK 5.2 Estimate the volume of $10000$ water molecules under standard conditions.

$V=\frac{10000}{3.34\times 10^{22}}=2.99\times 10^{-19}mL$

JC: All maths and calculations correct and nicely laid out.

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

The atom ends up at position $(1.2,1.1,0.7)$ which is outside of a 1 unit cubic box. After the periodic boundary conditions have been applied, the atom ends up at $(0.2,0.1,0.7)$ as one of its replicas enters the box through the opposite face.

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

$r^{*}$ in real units:

$$\begin{gathered} r=r^{*}\times \sigma =3.2\times 0.34=1.09 \\ nm \end{gathered}$$

Well depth in $$\begin{gathered} {\mathrm{k}\mathrm{J} \\ \mathrm{m}\mathrm{o}\mathrm{l}}^{-1} \end{gathered}$$ :

The well depth ($\varphi \left(r\right)$)$=-ϵ$

$ϵ=k_B\times 120=1.66\times 10^{-21}J$

$\varphi \left(r\right)=-1.66\times 10^{-21}J=-1.66\times 10^{-24}kJ$

$\varphi \left(r\right)=-1.66\times 10^{-24}\times N_A=-0.998kJmo{l}^{-1}$

$T^{*}=1.5$ in real units:

$T=\frac{T^{*}\times ϵ}{k_B}$

$T=1.5\times 120$

$$\begin{gathered} T=180 \\ K \end{gathered}$$

JC: Values are correct, but don't give the answers to more significant figures than are given in the question (2 in this case).

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?

If two atoms happened to be generated too close together, there will be repulsion between the electrons in the orbital. The two atoms will be pushed apart with a constant distance between. This situation is similar to atoms in solid state, where atoms are arranged in a regular way and has long range order. If two atoms happened to be generated far apart, the arrangement of atoms will be similar to a gas. In liquid state, atoms are close together and arranged in a random way, therefore the starting coordinates of the atoms need to be specified to make sure the simulation is simulating atoms in liquid state.

JC: The simulations are completely classical so electron electron repulsion is not simulated explicitly here. Your explanation is unclear. The problem with random coordinates is that if 2 atoms are placed close together the high repulsion can make the simulation unstable and will require a much smaller timestep.

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?

$$\begin{gathered} Number \\ density=\frac{Number \\ of \\ lattice \\ point}{Volume \\ of \\ unit \\ cell}=\frac{Number \\ of \\ lattice \\ point}{Lengt{h}^3} \end{gathered}$$

For simple cubic lattice, number of lattice point = $1$

$$\begin{gathered} Number \\ density=\frac{1}{1.0072^3}=0.8 \end{gathered}$$

For face-centred cubic lattice, number of lattice points = $4$, assume length = $a$

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

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

$$\begin{gathered} a=1.494 \\ in \\ reduced \\ units \end{gathered}$$

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?

In one unit of a face-centred cubic lattice, each atom at the corner contributes 1/8 to the unit cell and each atom on the face contributes 1/2 to the unit cell. This gives a total of 4 atoms per unit cell (8 x 1/8 + 6 x 1/2 = 4).

Since there are 1000 unit cells, there will be 4000 atoms (1000 x 4 = 4000) created by the create_atoms command.

JC: Correct.

TASK 11 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

mass 1 1.0

The mass command defines the mass of atoms. The number 1 between "mass" and "1.0" indicates that the simulation will contain only one type of atom. All atoms are identical and have a mass of 1.0 in reduced unit.

pair_style lj/cut 3.0

This command determines the potential between the particles. "lj" means that only Lennard Jones Potential is considered. "3.0" is the distance units between two particles. "cut 3.0" means that any interactions between particles that are more than 3 distance units apart will be negligible and the potential will equal to zero.

pair_coeff * * 1.0 1.0

This command determines the force field coefficient. The two numbers are the two parameters in the Lennard Jones Potential. The former "1.0" stands for the well depth ($ϵ$), the attraction between two atoms, and the latter "1.0" stands for the distance at which the potential between two atoms is zero ($\sigma$).

JC: What do the "*" mean?

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?

The velocity-Verlet algorithm as shown 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$

The velocity-Verlet algorithm is a modified version of Verlet's algorithm which takes velocity into account and can be used to calculate the kinetic energy of the system.

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

In both cases, the output values will be the same. However, the first method shown allows the timestep to be varied for different simulations in a much quicker way. By adding the second line, any variables that depend on the timestep value will change accordingly when the timestep value changes instead of changing the value manually line by line. The can reduce the probability of making mistakes in simulations.

JC: Good explanation.

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?

File:Total energy vs Time.PNG

Figure 6. Total energy vs Time

File:Temperature vs time.PNG

Figure 7. Temperature vs Time

File:Pressure vs Time.PNG

Figure 8. Pressure vs Time

After t=0.37, the energy remains fairly stable and mainly fluctuates between -3.183 and -3.185 (Figure 6). This happens when the simulation reaches equilibrium. Similar results were observed in Temperature vs Time (Figure 7) and Pressure vs Time (Figure 8). Both temperature and pressure fluctuate within a small range after t=0.37.

A plot showing energy versus time for 5 different timesteps was shown in Figure 9. Timestep = 0.0025 is the largest to give acceptable results. Timestep = 0.001 is the one with shortest timestep and therefore the most accurate results of the simulation. However, short timesteps mean that the same number of simulation steps cover a shorter amount of actual time and is very unhelpful when observation over a long time is needed. Similar results were obtained for timestep = 0.0025 which means that the data obtained using this timestep is still acceptable. Moreover, using timestep = 0.0025 require less processing power compared to timestep = 0.001. Timestep = 0.0075 and 0.001 also reaches equilibrium, however, the average energy value was higher and fluctuates more than the ones obtained with timestep = 0.001 and 0.0025.

Timestep = 0.015 is a particularly bad choice as it is the one with largest timestep and the only one that does not give similar result as the other four choices. The simulation with timestep = 0.015 does not reach equilibrium. Instead, the total energy increases linearly while fluctuating. The result shows that at such large interval of measurement, the total energy increases over time. At large timestep, significant errors in the simulation was observed. This was due to two neighbouring atoms being too close together when timestep = 0.015 and thus experience significant repulsive force leading to the unstabilised and deviated result. The graph of timestep = 0.015 also shows that error in the simulation accumulates over time, resulting in an increase in the deviation of total energy from the other four choices of timestep over time.

JC: Good analysis, timesteps of 0.0075 and 0.01 are not suitable as energy should not depend on timestep.

File:Total energy vs time at different timestep.PNG

Figure 9. Total energy vs time at different timestep

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.

Table 2. Conditions for simulations

Temperature	Pressure	Timestep
1.6	2.5	0.0025
1.9	2.5	0.0025
2.2	2.5	0.0025
2.5	2.5	0.0025
2.8	2.5	0.0025
1.6	2.7	0.0025
1.9	2.7	0.0025
2.2	2.7	0.0025
2.5	2.7	0.0025
2.8	2.7	0.0025

TASK16 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$ (1) $\frac{1}{2}\sum _i{m}_i{\left(\gamma v_i\right)}^2=\frac{3}{2}N{k}_B\mathfrak{𝔗}$ (2)

Solve these to determine $\gamma$.

Rearrange equation (2):

${\gamma }^2\frac{1}{2}\sum _i{m}_i{v}_i^2=\frac{3}{2}N{k}_B\mathfrak{𝔗}$

From equation (1):

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

Substitute equation (1) into equation (2) and rearrange to obtain $\gamma$:

${\gamma }^2\frac{3}{2}N{k}_{B}T=\frac{3}{2}N{k}_B\mathfrak{𝔗}$

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

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

JC: Good derivation.

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 values (100 1000 100000) specified are used to find the final average. The first value (100) represents Nevery and means use input values every this many timesteps. In this case, every 100 timesteps. The second value (1000) represents Nrepeat and means number of times to use input values for calculating averages. In this case, 1000. The third value (100000) stands for Nfreq which means calculate averages every this many timesteps.

The values of temperature, etc., will be sampled every 100000 timesteps for the average. Input values are used every 100 timesteps and therefore 1000 measurements contribute to the average.

The time stimulated will be dependent on the timestep chosen which will be 0.0025 in this simulation.

0.0025 x 100000 = 250 time in reduced units

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?

The ideal gas law equation was shown below:

$PV=N{k}_{B}T$

Since we are working in reduced units in the simulation, $k_B=1$

$\frac{P}{T}=\frac{N}{V}$

The density was found by pressure divided by temperature and substituting corresponding results obtained from simulations.

File:Density vs Temp at different Pressure.PNG

Figure 10. Density vs Temp at different Pressure

Both the simulated density and ideal gas law density are plotted against temperature at two different pressures (2.5, 2.7) in Figure 10. Both results showed that as pressure increases, the density increases. The simulated density is lower than the ideal gas law density. This is because ideal gas law density uses kinetic theory and an assumption of any interaction between particles is zero is made. This means that particles are tightly packed together and thus the density of ideal gas law is always higher. In the case of simulated density, the atoms are Lennard Jones particles, when the atoms are too close together, there will be strong repulsion between the atoms and cause the atoms to be further apart and thus lower density. Both simulated and ideal gas law results showed that as the temperature increases, the density decreases. This is because as temperature increases, atoms gain more kinetic energy and move away from the equilibrium position. At higher pressure, the density is always higher and this is caused by at high pressure, the atoms are forced to be closer together and thus result in a higher density. Error bars in both x and y direction obtained from simulation were included. However, the error bars in y direction are too small to be seen clearly. It is shown that as the temperature increases, the size of the error bar increases. This was suggested to be due to an increase in temperature causes the temperature to fluctuate more and therefore introduce a greater error during measurements. The discrepancy between simulated and ideal gas law densities increases with pressure. As pressure increases, the atoms are closer together and there will be more interaction between them and thus behave less like an ideal gas, which give rise to larger discrepancy.

A more sophisticated way to predict the density as temperature increases is to use the van der Waals equation shown below:

$(P+\frac{a{n}^2}{v^2})(V-nb)=nRT$

Unlike ideal gas law, the van der Waals equation takes into account attractive forces between atoms by including the correction term $\frac{a{n}^2}{v^2}$ and van der Waals coefficient $a$ and $b$ which depends on the type of atoms. The larger the attractive forces, the larger the coefficient $a$.

JC: Good explanation of the effects that temperature and pressure have on the discrepancy between the simulated and ideal gas results. Using the Van der Waals equation is great suggest to improve the agreement.

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.

Heat capacity can be calculated using the equation below:

$C_V=\frac{\partial E}{\partial T}=\frac{\mathrm{V}\mathrm{a}\mathrm{r}\left[E\right]}{k_B{T}^2}=N^2\frac{⟨E^2⟩-{⟨E⟩}^2}{k_B{T}^2}$

Where $\mathrm{V}\mathrm{a}\mathrm{r}\left[E\right]$ is the variance in $E$ and $N$ is the number of atoms.

A system at equilibrium by melting a crystal was prepared. 10 simulations including two different densities (0.2 and 0.8) and temperature range of 2.0 to 2.8 were run in the NVT ensemble.

Example input script for the simulation of ${\rho }^{*}$ = 0.2 and $T^{*}=2.0$ is shown below:

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

File:CvV vs Temperature.PNG

Figure 11. Cv/V vs Temperature

Heat capacity is a physical property and is related to how easy it is to excite the system to a higher energy state. The trend shown in Figure 11 is the one expected when $C_V/V$ plotted against temperature. As temperature increases, heat capacity decreases for both densities. This is due to an increase in temperature causes atoms to have more kinetic energy and thus a greater proportion of atoms are populated to a higher energy state. Energy levels at higher state are closer together and have a smaller separation; therefore less energy is required to promote another atom at higher energy level. Since the volume is constant at ${\rho }^{*}$ = 0.2 and 0.8 and heat capacity decreases, overall, $C_V/V$ decreases. The point at $T^{*}=2.4$ is an anomaly which might be because of Brownian motion of the liquid.

JC: Good explanation of decrease in heat capacity with temperature, but the point at 2.4 is not necessarily an anomaly. Can you explain the difference in heat capacity for different densities?/span>

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?

Table 3. Conditions for the Lennard-Jones system.

	Temperature	Density
Solid	1.4	1.4
Liquid	1.3	0.8
Gas	1.3	0.05

The radial distribution function shows how density changes as a function of interatomic separation, $r$. It gives the probability of finding an atom in the distance $r$ from another atom.

Looking at Figure 12, all three RDFs show three typical features:

File:Radial Distribution Function.PNG

Figure 12. Radial Distribution Function

1. When $r$ is small, the RDF is zero. This is due to strong repulsive force between the two positively charged nuclei and negatively charged electron clouds.
2. Obvious peaks were shown, indicating the chance of finding a pair of atom at a particular distance.
3. At long distances, the RDFs approach 1 as the overall density is equal to the local density.

The RDF of solid has the most number of peaks whereas the RDF of vapour only has one peak. The RDF of liquid is intermediate between solid and vapour, with a small number of peaks compared to the RDF of solid.

Atoms in solid are close together and are arranged in a regular way with regular spacing. Neighbouring atoms are found at a specific distance as indicated by the peaks in the RDF. Since atoms are close together in solid, there is a greater probability of finding a neighbouring atom compared to liquid and vapour, where atoms are further apart, and thus there are more peaks observed in RDF.

Atoms in liquid are close together, but arranged in a random way with irregular spacing. The RDF of liquid is intermediate between solid and vapour. It has much less peak than solid RDF, but more peaks than vapour RDF. The 4 distinct peaks showed that these are the 4 distances that are of favourable separation from another atom. Since atoms are less closer together compared to solid and have a lower probability of finding a neighbouring atom, the number of peaks observed decreases.

For vapour, only one peak was observed and the RDF is very smooth compared to the other two. This is because atoms in vapour are far apart and are arranged in a random way. The short range order disappears completely.

The face center cubic lattice structure has 4 lattice points are all in the same environment. The $r$ for the first three peaks of the solid RDF are 0.975, 1.425 and 1.725. These represents the three nearest atoms respectively. Distance from nearest atom: first layer < second layer < third layer. In the first layer, there are 12 neighbours with distance 0.975 (nearest neighbour) (Figure 13). In the second layer, there are 6 neighbours with distance 1.425 (second nearest neighbour) (Figure 14). In the third layer, there are 24 neighbours with distance 1.725 (Figure 15). The lattice spacing between the atom and its neighbour are calculated using Pythagoras theorem.

JC: The number of peaks in the RDF indicate how ordered a particular phase is, the position and areas of the peaks indicate how densely packed that phase is.

File:1st layer 12 neighbour.PNG

Figure 13. 1st layer 12 neighbour [1]

File:2nd layer 6 neighbour.PNG

Figure 14. 2nd layer 6 neighbour [1]

File:3rd layer 24 neighbour.PNG

Figure 15. 3rd layer 24 neighbour [1]

Table 4. Summary table of coordination number, distance and lattice spacing

Distance	Coordination Number	Distance (in reduced units)	Lattice spacing
Nearest neighbour	12	0.975	1.379
Second nearest neighbour	6	1.425	1.425
Neighbour	24	1.725	1.408

Average lattice spacing = (1.379 + 1.425 + 1.408)/3 = 1.404 in reduced units

JC: Good idea to take an average, show how you calculated the lattice spacing from each peak.

File:Int g(r) dr.PNG

Figure 16. Int g(r) dr

Figure 16 below shows $\int g(r)\mathrm{d}r$ of solid with $r$ equals 0 to 2. It can be used to verify the coordination number for each of the first 3 peaks.

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.

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.

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

For all three phases, solid, liquid and vapour, the total MSD is plotted against timestep (Figure 17-19) and time (timestep x 0.002) (Figure 20-22). The gradient of the linear behaviour is found by using Excel trendline option. The gradient of the linear line for Total MSD against time for each phases are then used to find the diffusion coefficient, $D$ using the equation below:

$D=\frac{1}{6}\times gradient$

File:MSD of solid vs Timestep.png

Figure 17. MSD of solid vs Timestep

File:MSD of liquid vs Timestep.png

Figure 18.MSD of liquid vs Timestep

File:MSD of vapour vs Timestep.png

Figure 19. MSD of vapour vs Timestep

File:MSD of solid vs Time.png

Figure 20. MSD of solid vs Time

File:MSD of liquid vs Time.png

Figure 21. MSD of liquid vs Time

File:MSD of vapour vs Time.png

Figure 22. MSD of vapour vs Time

JC: To get the diffusion coefficient from the MSD you should only fit a line to the linear part of the MSD graph as at shorter times the system has not reached the diffusive regime.

$D$ for vapour > liquid > solid is expected and this is the case observed in both simulated and one million atom results (Table 5). The results in the table above also shown agreement between simulated and much larger system results. MSD models the displacement of an atom in a finite volume. $D$ for solids are very small as atoms are arranged regularly and closely packed, atoms have minimal or no movement. The change in MSD is due to atoms vibrating on its lattice site and as atoms vibrate, neighbouring atoms will be affected via intermolecular interactions and spread vibrational energy through the lattice structure. This explains why the MSD fluctuates throughout the simulated time frame. In liquid, atoms are arranged randomly and there are more spaces in between compared to solid state. Atomic motions like Brownian motion allows continuous collision between atoms and give rise to diffusive behaviour in liquid. Liquid has greater displacement of atoms compared to solid and thus a larger $D$ is observed. In vapour, atoms are far apart with lots of spacing in between and there is not much intermolecular interactions. Atoms are free to travel in the simulation box and after a collision, there will be a much greater displacement and thus the largest $D$ is observed

Table 5. D of solid, liquid and gas from MSD

D from MSD	Solid	Liquid	Gas
Simulated result	$5*10^{-6}$	$0.0968$	$2.80$
1 Million Atom result	$5*10^{-6}$	$0.0873$	$2.54$

Units of the diffusion coefficient, $D$, is distance squared/ unit time.

TASK 23 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!). 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.

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

Position of an atom following the classical harmonic oscillator equation:

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

Differentiate to get the velocity of the atom.

$v=\frac{\mathrm{d}\mathbf{x}\left(t\right)}{\mathrm{d}t}=-A\omega \sin (\omega t+\varphi )$

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

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

$A^2$ and ${\omega }^2$ are taken out and cancelled

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

Using the following trigonometric 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}$

sin is an odd function and cosine is an even function. When an odd function is multiplied by an even function, the result is an odd function.

The integral of an odd function within symmetrical limits equal 0:

$\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+\theta )\mathrm{d}t}=0$

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

JC: Good, clear derivation.

File:VACF for solid, liquid and harmonic oscillator.png

Figure 23. VACF for solid, liquid and harmonic oscillator

Minima in the VACFs for the liquid and solid system represent collisions. (Figure 23) Velocity of atoms changes after collision and because the velocity of an atom before and after collisions is different, the VACF decorrelates. The differences between the liquid and solid VACFs are due to different ways the atoms are packed. Solids are closely packed and regularly arranged and atoms cannot move significantly away from their lattice site to collide with another atom, change in velocity is smaller. More oscillations were observed in the VACF of solid. There are both positive and negative correlations before the velocity drops to zero. The reverse in velocity after each collision is shown as a minimum. For liquid, atoms are randomly arranged and have spaces in between; atoms can diffuse more and collide with another atom. Therefore no oscillating behavior that is seen in solid is observed. The velocity decays rapidly with time and only one minimum is observed in the VACF for liquid. The harmonic oscillator VACF is very different to the Lennard Jones solid and liquid because the motion is periodic, as shown on the graph. There is a specific velocity at an instant time. There is no collision between atoms and velocity of the atom does not change until it hits the boundary, which reverses the velocity. Atoms in Lennard Jones solid and liquid collide and thus are different from harmonic oscillator VACF.

JC: This explanation is a bit vague, but you are right that collisions between particles cause the VACFs in the simulations to decorrelate, whereas this doesn't happen with the harmonic oscillator.

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?

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

$D=\frac{1}{3}\underset{0}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}\tau ⟨\mathbf{v}\left(0\right)\cdot \mathbf{v}\left(\tau \right)⟩$

Table 6. Summary of D obtained from VACF and MSD for smaller and larger system.

	D from MSD (simulated)	D from MSD (one million atoms)	D from VACF (simulated)'	D from VACF (one million atoms)
Solid	$5*10^{-6}$	$5*10^{-6}$	$-4.56*10^{-4}$	$-2.79*10^{-4}$
Liquid	$0.0968$	$0.0873$	$0.0794$	$0.0825$
Vapour	$2.80$	$2.54$	$0.813$	$0.978$

Using the trapezium rule, the running integral of solid, liquid and vapour are plotted against timestep (Figure 24-39) and the diffusion coefficients are found by multiplying the integral by $\frac{1}{3}$.

Expected result of $D$ for vapour > liquid > solid is observed in both simulated and one million atom results (Table). The results in the table below show agreement between simulated and much larger system results.

The largest source of error in estimates of $D$ from VACF is the function $C\left(\tau \right)=⟨\mathbf{v}\left(t\right)\cdot \mathbf{v}(t+\tau )⟩$ only takes into account of velocity and does not take into account the initial position of the atoms. Also, the VACF $D$ is estimated using trapezium rule. The trapezium rule is a way of estimating the area under a curve. It works by splitting the area under a curve into a number of trapeziums, which can be calculated then adding the area of the trapeziums to obtain an approximation of the integral. The accuracy of the area is dependent on the height of the trapezium, in this case, the timestep. To increase the accuracy of the result, smaller increments can be used for calculation.

JC: The initial positions of the atoms should have no effect on the diffusion coefficient as the diffusion coefficient should be measured when the system has reached equilibrium.

File:Running integral of vapour (solid).png

Figure 24. Running integral of solid (simulated)

File:Running integral of liquid(simulated).png

Figure 25. Running integral of liquid(simulated)

File:Running integral of vapour (simulated).png

Figure 26. Running integral of vapour (simulated)

File:Running integral of solid (one million atoms).png

Figure 27. Running integral of solid (one million atoms)

File:Running integral of liquid (one million atoms).png

Figure 28. Running integral of liquid (one million atoms)

File:Running integral of vapour (one million atoms).png

Figure 29. Running integral of vapour (one million atoms)

Reference

[1] FCC And Its Neighbours, http://www.ifmpan.poznan.pl/~urbaniak/fcc%20and%20its%20neigbours01.pdf, 24 Feb 2016.