Talk:Mod:xj1213

JC: General comments: Good report on the whole, you cover all of the tasks in the experiment well. Try to extend and improve your explanations of the underlying theory to demonstrate that you really understand the material and the point of each task.

THEORY

Numerical Integration

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

File:Xj1213positionvstime.png

Figure 1: Position Vs Time

All three column of 'ANALYTICAL', 'ERROR', and 'ENERGY' were calculated and plotted using the following equations:

File:Xj1213error.png

Figure 2: Errors vs Time

'ANALYTICAL' =

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

(Figure 1)

File:Xj1213energyvstime.png

Figure 3: Energy Vs Time

'ERROR' = |'ANALYTICAL' - $x\left(t\right)$| (Figure 2)

'ENERGY' = $\frac{1}{2}{m}_i{\mathbf{v}}_i^2$ (Figure 3)

JC: Energy should be total energy, not just kinetic energy.

$x\left(t\right)$ is the position for velocity-Verlet solution and the 'ANALYTICAL' is the classical solution position for the position at time $t$.

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.

File:Xj1213maxerror.png

Figure 4: Maximum Errors vs Time

JC: Why does the error graph have this oscillatory shape, but with the maximum error increasing over time?

The Error maxima was collected from 5 time ranges (0-3.20, 3.20-6.30, 6.30-9.40, 9.40-12.60) and the maxima vs Time graph was plotted.(Figure 4) The equation of the best fit line is shown on the graph.

JC: Use 3 s.f. for equation and R squared value in graph.

Time	Max
2.00	0.000758
4.90	0.002008
8.00	0.003301
11.10	0.004604
14.20	0.005911

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?

The timestep used to ensure the total energy which does not change by more than 1% is 1.98. (Figure 5)

Due to conservation of energy, the total energy of a system should be constant.

File:Xj12131perenergychange.png

Figure 5: Energy vs Time when timestep is changed

Atomic Force

Task 4

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

1. The separation, $r_0$

To find the separation, $r_0$, the Lennard-Jones potential energy is zero:

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

$r_0$ = $\sigma$

2. The Force at $r_0$

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

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

$\mathbf{F}=-24ϵ(2\frac{{\sigma }^{12}}{r^{13}}-\frac{{\sigma }^6}{r^7})$

Substitute $r_0$ = $\sigma$ into the equation

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

3. The equilibrium separation, $r_{eq}$

When the atom is in equilibrium, the potential energy of the atom is at the minimum. In order to do so, differentiate the Lennard-Jones potential energy, $\varphi \left(r\right)$

To find the minimum,

$\frac{\mathrm{d}\varphi \left(r\right)}{\mathrm{d}\mathbf{r}}=0$

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

$2={\left(\frac{r}{\sigma }\right)}^6$

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

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

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

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

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

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

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

5. 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$

$\underset{a\sigma }{\overset{b}{\int }}\varphi \left(r\right)\mathrm{d}r=\underset{a\sigma }{\overset{b}{\int }}4ϵ(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})\mathrm{d}r$, where a and b are constant

$\underset{a\sigma }{\overset{b}{\int }}\varphi \left(r\right)\mathrm{d}r=-\frac{4ϵ\sigma }{11}(\frac{1}{b^{11}}-\frac{1}{a^{11}})+\frac{4ϵ\sigma }{5}(\frac{1}{b^5}-\frac{1}{a^5})$

When a = 2, b = $\mathrm{\infty }$

$\underset{2\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=-\frac{4ϵ\sigma }{11}(0-\frac{1}{2^{11}})+\frac{4ϵ\sigma }{5}(0-\frac{1}{2^5})=-0.02482$

When a = 2.5, b = $\mathrm{\infty }$

$\underset{2.5\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=-\frac{4ϵ\sigma }{11}(0-\frac{1}{2.5^{11}})+\frac{4ϵ\sigma }{5}(0-\frac{1}{2.5^5})=-0.008177$

When a = 3, b = $\mathrm{\infty }$

$\underset{3\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=-\frac{4ϵ\sigma }{11}(0-\frac{1}{3^{11}})+\frac{4ϵ\sigma }{5}(0-\frac{1}{3^5})=-0.003290$

JC: All of the maths is correct, however there is a factor of $ϵ\sigma$ missing from the integral answers.

Periodic Boundary Conditions

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.

1. The number of water molecules in 1 ml of water under standard conditions, N

The density of liquid water under standard condition is approximately 1.0 g/cm. The molar mass of a water molecule is 18 g/mol.

Therefore, the mass of 1 mL of water is equal to 1.0 g/cm times 1 mL, which is 1.0 g.

Thus, number of moles = 1.0/18 = 0.056 mol

N = 0.056 x 6.02 x 10^23 = 3.34 x 10^22

2. The volume of $10000$ water molecules under standard conditions, $V_{10000}$

$\frac{3.34*10^{22}}{1mL}=\frac{10000}{V_{10000}}$

$V_{10000}=\frac{10000mL}{3.34*10^{22}}$

$V_{10000}=2.994*10^{-19}$

JC: Your answers are correct, but you need to be careful with units, density is measured in g/cm³ = g/ml since 1 ml = 1cm³. Also when you give your answer for the volume you need to give units.

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

For an atom positioned at $(0.5,0.5,0.5)$ initially then move alone the vector $(0.7,0.6,0.2)$, the new position will be $(1.2,1.1,0.7)$, which is outside of the boundary of $(0,0,0)$ to $(1,1,1)$. However, since the atom is in a box which is surrounded by repeated box infinitely in all directions, when an atom moves out of the boundary of the box, one of its replicas enters the box through the opposite face. In this way, the number of atoms inside the box stays constant. Therefore, the atom (which will be a replica in this case) will end up at position $(0.2,0.1,0.7)$.

Reduced Units

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?

• distance $r^{*}=\frac{r}{\sigma }$
• energy $E^{*}=\frac{E}{ϵ}$
• temperature $T^{*}=\frac{k_{B}T}{ϵ}$

Distance $r=r^{*}\sigma$ = $3.2$ x $0.34\mathrm{n}\mathrm{m}$ = $1.088nm$

$ϵ=k_B*120K$ = $1.38*10^{-23}J/K*120K$= $1.656*10^{-21}J$

The well depth $\varphi \left(r\right)=4ϵ(\frac{{(0.34\mathrm{n}\mathrm{m})}^{12}}{{(1.088nm)}^{12}}-\frac{{(0.34\mathrm{n}\mathrm{m})}^6}{{(1.088nm)}^6})$ = $6.16*10^{-24}J$

Temperature T = $$\begin{gathered} ϵT^{*} \\ / \\ {k}_B=180K \end{gathered}$$

JC: You have already shown that the well depth is $-ϵ$. Instead you have calculated the value of the potential at the cut off, a long way from the well minimum. The well depth should also have been given in kJ mol^-1. The values for r and T are correct.

EQUILIBRATION

Creating the Simulation Box

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 we give atoms random starting coordinates, there will be possibility that some of the atoms will be at the exact same position or will be very close together which is unlikely to happen. If two atoms happen to be generated close together there will be a very strong repulsion between the two atoms.

JC: Correct, the strong repulsion forces between these atoms mean that a smaller timestep is needed to run the simulation without it crashing.

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?

In a simple cubic lattice, there is only one lattice point in a unit cell. Lattice spacing in x,y,z = 1.07722 1.07722 1.07722, therefore the volume of the unit cell is $1.0772^3$.

The number density of lattice points = $\frac{1}{1.0772^3}=0.8$

In a face-centred cubic lattice there are 4 lattice points in a unit cell.

The side length of the cubic unit cell $x$ can be calculated from:

$1.2=\frac{4}{x^3}=0.8$

$x=1.4938$

JC: Your answer is correct, but try to be more careful with your presentation, what does $1.2=\frac{4}{x^3}=0.8$ mean?

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?

The face-centred cubic lattice has 4 lattice points, therefore for a box containing 1000 unit cells there will be 4000 atoms.

JC: Correct.

Setting The Properties of The Atoms

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

1. 'mass' is a command used to set the mass for all atoms of one or more atom types. In this case 1 means one atom type with the mass = 1.0

2. In this line it gives a command of setting the 'pair style' to lj/'cut style which is the standard 12/6 Lennard-Jones potential:

$E=4ϵ(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})$, where $r<r_c$

$r_c$ is the cutoff for Lennard Jones interaction. In this case $r_c=3.0$. The curoff means the the bond distance where the Lennard-Jones potential almost reaches zero. In order to save time the simulation stops after the cutoff.

JC: The simulation doesn't stop after the cut off, but the interaction energy is not calculated for atoms outside the cut off distance. This helps to speed up the calculation.

3. For all of the lj/cut pair styles, pair_coeff command must be used to define coefficients for each pair of atoms types.

A wildcard asterisk (* *) is used to set the coefficients for all pairs of atom types, which are followed by the values indicating the coefficients for the pairs of atom types. The coefficients used here is (1.0 1.0) which indicates the epsilon and sigma respectively in the Lennard-Jones potential.

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.

JC: Correct, but can you give a sentence to explain how you know this?

Running The Simulation

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

The former one is more convenient where the only variable that we need to type in is 'timestep'; all the other variables are dependent on the timestep. While the later one the number of steps is independent on the timestep and so we need to change every line chi achieve what the former one does.

Task 14

File:Xj1213totalengvstime0001.png

Figure 6: The Total Energy Vs Time when timestep is 0.001

File:Xj1213tempvstime0001.png

Figure 7: The Temperature Vs Time when timestep is 0.001

File:Xj1213pressvstime0001.png

Figure 8: The Pressure Vs Time when timestep is 0.001

Make plots of the energy, temperature, and pressure, against time for the 0.001 timestep experiment. Does the simulation reach equilibrium? How long does this take?

From the graphs (Figure 6, Figure 7 and Figure 8) we can see that the simulation reaches equilibrium when the time is approximately 0.3. After t = 0.3 all three properties reach a stable value with slight fluctuation.

When you have done this, make a single plot which shows the energy versus time for all of the timesteps.

File:Xj1213totengvstimedifftimestep.png

Figure 9: The Total Energy Vs Time with 5 different timesteps

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 largest acceptable timestep is 0.0025; it gives the same fluctuation as with 0.001 meaning that it is accurate enough but not as small as 0.001 so we have more time to observe. Timestep of 0.015 is a particularly bad choice because the the total energy is not constant.

JC: Good choice, but can you be a bit clearer in your justification, what does "it gives the same fluctuation" mean?

RUNNING SIMULATIONS UNDER SPECIFIC CONDITIONS

Temperature and Pressure Control

Thermostats and Barostats

Task 15

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{\gamma }^2{v}_i^2=\frac{3}{2}N{k}_B\mathfrak{𝔗}$

Because $\gamma$ is a constant: ${\gamma }^2\left(\frac{1}{2}\sum _i{m}_i{v}_i^2\right)=\frac{3}{2}N{k}_B\mathfrak{𝔗}$

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

We get: ${\gamma }^2\left(\frac{3}{2}N{k}_{B}T\right)=\frac{3}{2}N{k}_B\mathfrak{𝔗}$

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

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

JC: Correct.

Examining the Input Script

Task 16

Choose 5 temperatures (above the critical temperature $T^{*}=1.5$), and two pressures. 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.

The temperatures used are: 1.6, 1.8, 2.0, 2.2, 2.4

The pressures used are: 2.57 and 2.69 (Chosen based on the previous simulations)

The timestep used is: 0.0025

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?

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

100 -> Nevery = use input values every 100 timesteps

1000 -> Nrepeat = 1000 of times to use input values for calculating averages

100000 -> Nfreq = calculate averages every 100000 timesteps

In this case, Nevery=100, Nrepeat=1000, and Nfreq=100000, then values are taken every 100 timesteps (i.e. on timesteps 100, 200, 300, 400, 500...1100, 1200, 1300, 1400...100000), which will be used to compute the final average on timestep 100000. That means we will have 1000 samples to be used to calculate the average. The next line states 'run 100000', which means the simulation will run 100000 timesteps.

Plotting the Equations of State

Task 17

File:Xj1213density temp pre.png

Figure 10: The density as a function of temperature for both of the pressures that we simulated with error bars and best-fit lines

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)$. Plot the density as a function of temperature for both of the pressures that you simulated. Your graph(s) should include error bars in both the x and y directions. You should also include a line corresponding to the density predicted by the ideal gas law at that pressure. Is your simulated density lower or higher? Justify this. Does the discrepancy increase or decrease with pressure?

Figure 10 shows the density as a function of temperature for both of the pressures that we simulated. The higher the pressure, the higher the density the system has. This is because the higher the pressure, the smaller the volume. Density is inversely proportional to volume for a given mass, therefore an increase in pressure give rise to an increase in density.

File:Xj1213density temp pre with idealgas.png

Figure 11: The density as a function of temperature for both of the pressures that we simulated with with the comparison with thedensities of the ideal gas

Figure 11 shows the comparison of the simulated results with the density predicted by the ideal gas law at that pressure. It clearly shows that the simulated densities are lower than that of the ideal gas densities. This is because, for example, in real gas there are interactions such as attractions and repulsions. While in the ideal gas there is no interaction at all. That means the real gas cannot get as close as the ideal could because of repulsion force, which gives lower density of real gas then of ideal gas.

From Figure 11 it is also illustrated that the higher the pressure, the higher the density the system has (which is the same as for the simulated systems) resulting in the decrease in the discrepancy with the simulated results. Because ideal gas does not have interactions with each other, this phenomenon can be explained by entropy. Thermodynamically pressure can be expressed in terms of entropy:

$P=T{\left(\frac{\delta S}{\delta V}\right)}_{U,N}$

For a smaller volume, the way to arrange the particles decreases, giving rise to a decrease in entropy. Because the entropy of the system tends to increase therefore giving rise to a force on the wall of the box and thus the pressure.

From the ideal gas equation:

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

$P=\left(\frac{N}{V}\right)k_{B}T$, where $k_B=1$ when reduced.

$Density=\left(\frac{N}{V}\right)$

$P=\left(Density\right)k_{B}T$

Therefore, higher the pressure, higher the density.

JC: Good explanation of the origin of the pressure in an ideal gas. Do you expect simulations to agree better with the ideal gas model at lower or higher pressures?

CALCULATING HEAT CAPACITIES USING STATISTICAL PHYSICS

Heat Capacity Calculation

Task 18

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.

Below is an example of one of my input script with density 0.2 and temperature 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
variable pdamp equal ${timestep}*1000
fix nvt all nvt temp ${T} ${T} ${tdamp} 
run 10000
reset_timestep 0

### MEASURE SYSTEM STATE ###
unfix nvt
fix nve all nve
thermo_style custom step etotal temp press density vol atoms
variable temp equal temp
variable temp2 equal temp*temp
variable etotal equal etotal
variable etotal2 equal etotal*etotal
variable vol equal vol
fix aves all ave/time 100 1000 100000 v_temp v_temp2 v_etotal v_etotal2
run 100000

variable avetemp equal f_aves[1]
variable aveeng equal f_aves[3]
variable heatcap equal atoms*atoms*((f_aves[4]-(f_aves[3]*f_aves[3]))/f_aves[2])
variable heatcapvol equal ((atoms*atoms*((f_aves[4]-(f_aves[3]*f_aves[3]))/f_aves[2]))/vol)

print "Averages"
print "--------"
print "Temperature: ${avetemp}"
print "Etotal: ${aveeng}"
print "HeatCapacity: ${heatcap}"
print "HeatCapacityPerVol: ${heatcapvol}"

File:Xj1213hearcappervolwithdensnofit2.png

Figure 12: Heat Capacity per Volume as a Function of Temperature

JC: Script looks good and heat capacity is calculated correctly.

The equation we used to calculate the heat capacity $C_v$ is: $C_V=N^2\frac{⟨E^2⟩-{⟨E⟩}^2}{k_B{T}^2}$

where reduced $k_B=1$

Therefore: $\frac{C_V}{V}=\frac{N^2}{V}\frac{⟨E^2⟩-{⟨E⟩}^2}{k_B{T}^2}=(density)N\frac{⟨E^2⟩-{⟨E⟩}^2}{k_B{T}^2}$

As the above equations shown, $\frac{C_V}{V}$ is proportional to the densities and inversely proportional to $T^2$, therefore the trend shown in Figure 12 is exactly what we expect.

JC: Does the numerator also depend on temperature? What will happen to the energy fluctuations when temperature is increased.

However, in a more chemistry way to explain these trends, heat capacity is a measurement of the amount of energy required to raise the temperature of a system. With a higher density, i.e. more atoms in a unit volume, more energy is needed to heat up the entire system, thus higher the heat capacity the system has.

Moreover, the temperature dependency of the heat capacity is more complicated. From quantum chemistry we know that the higher the temperature a system has, the more atoms it is in the higher energy state. As the energy state goes up, the spacing between two energy states is smaller. Therefore, the energy needed to excite the system to a higher energy state is smaller due to the smaller spacing. This give rise to the decrease in heat capacity with increasing temperature.

JC: This is the right idea, but try to make your explanation a bit clearer.

STRUCTURAL PROPERTIES AND THE RADIAL DISTRIBUTION FUNCTION

Calculating $g(r)$ in VMD

Task 19

File:Xj1213grvsdistance2.png

Figure 13: A plot showing $g(r)$ varies with distance

File:Xj1213intgrvsdistance.png

Figure 14: A plot showing $\int g(r)\mathrm{d}r$ varies with distance

File:Xj1213intgrvsdistancezoom2.png

Figure 15: A plot showing $\int g(r)\mathrm{d}r$ varies with distance (zoomed in)

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 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 densities and temperatures for the three phases are chosen as following:

Vapour: density 0.1, temperature 1.25 sc

Liquid: density 0.8, temperature 1.2 sc

Solid: density 1.25, temperature 1.25 fcc

The Radial Distribution Function (RDF) $g(r)$ describes the variation of the density as a function of distance for a particular reference particle (r = interatomic separation). In other words, the RDF is an example of a pair correlation function, which describes, on average, how the particles in a systen are radially packed around each other.

As shown in Figure 13, for all three phases, the RDF is zero when r is small. This is due to the repulsion between atoms making them cannot go any closer than that (in this case r is around 0.775-0.825). This r indicates the effective atomic width.

In the solid phase where we have a higher density (blue line), there are many peaks shown at long range indicating that the atoms are packed continuously. This is because in the solid phase atoms are highly ordered and are closely packed next to each other.

In the liquid phase where we have a slightly low density (red line), there are much less peaks shown (only 3 obvious peaks) indicating that there are much less atoms at the neighbourhood due to continual movement of the atoms.

In the gas phase where we have a very low density (green line), there is only one peak shown due to the highly disordered atoms and the atoms are too far apart that the chance of meeting each other is very small.

At very long range every RDF tends to a value of 1 because at this range the RDF is governed by the average density of the system.

JC: Whether or not peaks are clear at large distances in the RDF depends really on the range over which the order in the system persists, rather than how many atoms there are at a particular distance or how close the atoms are together.

The $\int g(r)\mathrm{d}r$ shows the areas under each RDF. For solid which has a high density the chance of finding an atoms around the reference atoms is very high, therefore the area under the RDF is large. Similarly for liquid which has a slightly lower density, the end point of the curve is lower than that of the solid phase. For gas phase which has a very low density the probability of finding an atom around the reference is very low, therefore the curve is much lower than those for liquid and solid.

The RDF is very useful when identifying whether a system is in solid, liquid or gas phase.

File:Xj1213fccsolidstructurelab.png

Figure 16: Face Centred Cubic unit cell

JC: Great picture to illustrate which atoms are nearest neighbours.

The three atoms corresponding to the first three peaks on the solid RDF are illustrated in Figure 16, using the orange atom as the reference atom. The closest atoms is b which is $\frac{1}{\sqrt{2}}A$ away from the reference atom, where A is the length of the cubic. The second peak corresponds to atom a which has a distance of A between the reference atom. Atom c corresponds to the third peak which is $\sqrt{\frac{3}{2}}A$ away from the reference atom. These can be confirmed by calculating A from the three simulating peaks and compare with the theoretical A value.

Atom	Peak Distance	Relationship with A	A	Coordination Number
a	1.375	A	1.375	8
b	0.975	$\frac{1}{\sqrt{2}}A$	1.3788	12
c	1.675	$\sqrt{\frac{3}{2}}A$	1.3676	24

The average A is 1.37. Which is close to the theoretical $A^3=\frac{4}{1.25}$, $A=1.47$

From Figure 16 we can also find out the coordination number of the reference atom with respect to the three closest atoms respectively. In the same unit cell the closest atom b has 3 equivalent, each of which is shared between two unit cell. The reference atom is surrounded by 8 unit cells. Therefore the coordination number of the reference atom with b = $\frac{3}{2}*8=12$.

For the second closest atom a , it is clearly shown that the reference atom is attached to 6 different atoms with the same distance as a has. Therefore the coordination number of the reference atom with a is 6.

In the same unit cell the third closest atom c has 3 equivalent. The coordination number of the reference atom with c = $3*8=24$.

This can be confirmed by looking at the $\int g(r)\mathrm{d}r$ plot with distance in the region of 0-2 (Figure 15). The step-like graph indicates one peak for one step. The first step is 12. the second step is 18 which means coordination number is 18-12=6. The third step is 42 which shows that the coordination number is 42-12-6=24. This is consistent with the result stated above based on Figure 16.

DYNAMICAL PROPERTIES AND THE DIFFUSION COEFFICIENT

Mean Squared Displacement

Task 20

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.

File:Xj1213MSDsolid2.png

Figure 17: A plot for Mean Square Displacement vs Timesteps for 8000 atoms and 1 million atoms in solid phase

File:Xj1213MSDliquid.png

Figure 18: A plot for Mean Square Displacement vs Timesteps for 8000 atoms and 1 million atoms in liquid phase

File:Xj1213MSDvapour2.png

Figure 19: A plot for Mean Square Displacement vs Timesteps for 8000 atoms and 1 million atoms in gas phase

The diffusion coefficient is related to the mean square displacement:

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

Therefore $D=\frac{1}{6}*$ gradient of the MSD plot, with a unit in terms of timestep.

'	Gradient (timestep)	Gradient (time)	D
Gas (8000 atoms)	0.0204	1.02E+01	1.700E+00
Gas (1 million atoms)	0.0386	1.93E+01	3.217E+00
Liquid (8000 atoms)	0.001	5.00E-01	8.333E-02
Liquid (1 million atoms)	0.001	5.00E-01	8.333E-02
Solid (32000 atoms)	5.00E-08	2.50E-05	4.167E-06
Solid (1 million atoms)	5.00E-08	2.50E-05	4.167E-06

The graphs for Mean Square Displacement vs Time show an expected trend. For solid phase (Figure 17), the MSD's show small fluctuations over time as atoms in solid are closely packed, so they only tend to vibrate with a small distance. Therefore the diffusion coefficient for the MSD plot is almost zero. For liquid phase (Figure 18), the MDS increase almost linearly over time illustrating a greater freedom of movement for liquid. Atoms in gas phase have even greater freedom of movement compared to atoms in liquid as there is a more rapid increase in the MSD graph. The atoms in pas phase tends to stay as far apart from each other as possible. As the freedom of movement increases, the diffusion coefficient increases(Figure 19).

JC: Good reasoning, in solids atoms are localised to lattice points, but in liquids and gas phases atoms have freedom to move.

Velocity Autocorrelation Function

Task 21

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 a 1D harmonic oscillator is expressed as $x\left(t\right)=A\cos (\omega t+\varphi )$

Therefore the velocity is $v\left(t\right)=-A\sin (\omega t+\varphi )$

$C(\tau )=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}(-Aw\sin (wt+\theta ))(-Aw\sin (w(t+\tau )+\theta )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}(-Aw\sin (wt+\theta ){)}^2\mathrm{d}t}$

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

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

Because $\sin (a+b)=\sin (a)\cos (b)+\cos (a)\sin (b)$

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

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

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

Since $\sin (wt+\theta )$ is an odd function and $\cos (wt+\theta )$ is an even function, $\sin (wt+\theta )\cos (wt+\theta )$ is and odd function. Hence $\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (wt+\theta )\cos (wt+\theta )\mathrm{d}t=0$

Therefore $C(\tau )=\cos (w\tau )$

File:Xj1213VACFplot.png

Figure 20: The plot for Velocity Autocorrelation Function of solid and liquid with the plot of C(t)

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?

Each minima in the VACFs for the liquid and solid system represent the change of directions. In the solid system the change is due to the vibrations of the atoms, while in the liquid system the change is due to collisions of the atoms with each other. The VACF will decay over a long range of time. This is because the velocity decorrelates with time, which means the atom 'forgets' what its initial velocity was.

In solid structure, all the atoms are closely packed and their position is very stable. The only motion is oscillation; the atom vibrate backwards and forwards, reversing their velocity for each oscillation. Therefore the VACF oscillate strongly from positive to negative then positive again. However the oscillation will be disturbed by other factors and so over a period of time the function still decays.

Liquid behaves similarly to solid because the atoms are also very close to each other. The difference is in liquid the atoms can move freely and collide with each other, hence change of velocity. Therefore the VACF will decay very rapidly to zero giving maybe only one or two minima.

The harmonic oscillator VACF is very different to the Lennard Jones solid and liquid because the harmonic oscillator only considers an ideal situation where only one atom moving forwards and backwards without interactions with anything. This means that there is no transfer of energy and hence the velocity of the atom will remain the same with only the change of directions for each oscillation forever. On the other hand, the Lennard Jones system considers the interactions between the atoms; the atoms will collide with each other and hence change of energy and velocity. Therefore the velocity correlation will tend to zero over a long period of time.

JC: Good analysis, the key point is that the presence of collisions between atoms randomises particle velocities, causing them to decorrelate and for the VACF to tend to zero. There are no collisions in the harmonic oscillator VACF.

Task 22

File:Xj1213runningintegralgas80002.png

Figure 21: The Running Integral for Gas with 8000 atoms

File:Xj1213runningintegralgasm2.png

Figure 22: The Running Integral for Gas with 1 million atoms

File:Xj1213runningintegralliquid8000.png

Figure 23: The Running Integral for Liquid with 8000 atoms

File:Xj1213runningintegralliquidm.png

Figure 24: The Running Integral for Liquid with 1 million atoms

File:Xj1213runningintegralsolid320002.png

Figure 25: The Running Integral for Solid with 32000 atoms

File:Xj1213runningintegralsolidm2.png

Figure 26: The Running Integral for Solid with 1 million atoms

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

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

Therefore $D=\frac{1}{3}*$running integral

'	Running integral (timestep)	Running integral (time)	D
Gas (8000 atoms)	2605	5.21	1.737E+00
Gas (1 million atoms)	4902	9.804	3.268E+00
Liquid (8000 atoms)	120	0.24	8.000E-02
Liquid (1 million atoms)	123	0.246	8.200E-02
Solid (32000 atoms)	-0.044	-0.000088	-2.933E-05
Solid (1 million atoms)	0.068	0.000136	4.533E-05

The trend of D is exactly what we expect and is consistent with the D we obtained using MSD. The diffusion coefficient decreases from gas to liquid to solid because solid is fixed in one place while liquid is free to move and the gas has even more freedom of movement. The number of atoms also has influence on determining D because the more atoms the system has the more accurate the simulation. The main source of error with this method is that the simulation time may be not enough especially when simulating gas. Gas atoms have the highest degree of freedom when moving around so we expect it behaves most differently as time proceeds. Also the equation of calculating D include integrating the VACF from negative infinity to positive infinity while we only used 5000 timesteps.

JC: What about the error in the approximation of the VACF integral using the trapezium rule? How do the magnitudes of the diffusion coefficients compare with the values calculated using the MSD?