Talk:Mod:Sh1994

NJ: Take care with your written English. Some sections of this report either don't make sense, or are hard to interpret.

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

Plot	Description
x200px200px	Plot illustrating analytical position against time (red) and velocity-Verlet position against time (blue).
x200px200px	Plot demonstrating how energy of the system changes with time. NJ: Isn't this simulation output from later, not from the harmonic oscillator?
x200px200px	Plot showing how error (absolute difference between analytical and velocity-Verlet solution) varies 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.

x200px200px

At a timestep of 0.1, there is a positive linear relationship between the maxima error in position and time. This relationship can be modelled by the function of y=0.004x-0.00007, where y=maximum error and x=time. NJ: WHy do you think the error accumulates like this?

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?

Total energy at time=0 is 0.5. Ensuring that the total energy does not change by more than 1% over the course of simulation means it must be kept above 0.495. It is found that the timestep not greater than the value of 0.2 (</= 0.2) fulfills such criterion. It is important to monitor the total energy of a physical system during simulation since in the real physical system total energy is conserved. Therefore, in order to model a physical system accurately, its total energy should not vary significantly - it should resemble the real physical system and have a roughly constant total energy.

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

To find the separation, $r_0$, at which the potential energy is zero, we set $\varphi \left(r\right)=0$, and hence we obtain:

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

Rearranging above gives:

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

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

This means $r_0=\sigma$

Force acting on an object is determined by the potential that it experiences:

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

It shows force is equivalent to negative differential of Lennard-Jones potential. In order to calculate force at $r_0$, Lennard-Jones potential is differentiated and the negative differential is calculated with $r$ = $r_0$ being substituted:

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

When $r$ = $r_0$ = $\sigma$,

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

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

${\mathbf{F}}_i=4ϵ(\frac{6}{\sigma })$

${\mathbf{F}}_i=(\frac{24ϵ}{\sigma })$

At equilibrium, net force acting on the system is equal to zero. Equilibrium separation, $r_{eq}$, is reached when ${\mathbf{F}}_i=0$:

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

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

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

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

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

Well depth, $\varphi \left(r_{eq}\right)$, is calculated by substituting $r$ with $r_{eq}$ into the Lennard-Jones potential equation:

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

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

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

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

$\varphi \left(r_{eq}\right)=-ϵ$

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

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

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

Given $\sigma =ϵ=1.0$

Substituting upper limit of $\mathrm{\infty }$ and lower limit of $2\sigma$:

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

Substituting upper limit of $\mathrm{\infty }$ and lower limit of $2.5\sigma$:

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

Substituting upper limit of $\mathrm{\infty }$ and lower limit of $3\sigma$:

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

Under standard conditions, water has a density of 1g/mL. Therefore, 1mL of water weighs 1g. Molecular mass of water, H2O, is 18:

Number of moles = Mass / Molecular Mass:

Number of moles of water molecules in 1mL = 1 / 18 = 0.0556 moles.

Number of molecules = Avogadro's constant x Number of moles. Given Avogadro's constant = 6.022 x 10^23:

Number of water molecules in 1mL = 6.022 x 10^23 x (0.0556) = 3.35 x 10^22

There are 3.35 x 10^22 water molecules in 1mL of water.

Volume of 10000 water molecules = 10000 / 3.35x10^22 = 2.99 x 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?

Applying periodic boundary conditions means that the total number of atoms in a box stays the same. When an atom moves to the right, leaving its reference box, periodic boundary conditions ensure another atom enters the reference box from the left. In this scenario, an atom at initial position (0.5, 0.5, 0.5) has moved to the position of (0.2, 0.1, 0.7).

x-coordinate = 0.5 + 0.7 - 1.0 = 0.2

y-coordinate = 0.5 + 0.6 - 1.0 = 0.1

z-coordinate = 0.5 + 0.2 = 0.7

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=r*/sigma$

r = 3.2 x 0.34

r = 1.088nm

The well depth, $ϵ$, is calculated by:

$ϵ$ = ${k}_B$ x 120

$ϵ$ = 1.38 x 10^-23 x 120

$ϵ$ = 1.66 x 10^-21 J

This is the energy for 1 molecule. To work out energy in kJ/mol, we multiply this molecular energy by the Avogadro's number and then convert it from J/mol to kJ/mol:

$ϵ$ = 1.66 x 10^-21 x (6.02 x 10^23) / 1000

$ϵ$ = 0.997 kJ mol^-1

T = T* ($ϵ$ / ${k}_B$)

T = 1.5 x (120)K

T = 180K

NJ: Well presented section, well done.

Equilibration

Task 1

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 atoms were given a random starting position, there is a chance of them superimposing on each other in the simulation. This is not possible in reality and energy calculated will not be accurate. Also in the real physical system, atoms are spread out randomly. Problem arises if two atoms are too close together. This is because such arrangement increases the potential energy. This 'excess' potential energy is later converted into a kinetic energy in the simulation. As the result, energy calculated using simulation will be inaccurate. Therefore, atoms are placed on the lattice points rather than being positioned randomly.

NJ: Good - the resulting large forces tend to make your simulation "explode" as atoms fly apart from each other, and a very small timestep becomes necessary to treat their motion.

Task 2

Satisfy yourself that this lattice spacing corresponds to a number density of lattice points of 0.8. Consider instead a face-centred cubic lattice with a lattice point number density of 1.2. What is the side length of the cubic unit cell?

Number Density = Number of particles per unit cell / volume

This is a simple cubic lattice, containing 1 atom per unit cell:

Number Density = 1 / (1.07722 ^ 3)

Number Density = 0.79999... = 0.8

Face-centred cubic lattice has 4 lattice points/ 4 atoms. Given lattice point density is 1.2:

1.2 = 4 / volume

volume = 3.3333...

(Side length of the cubic unit cell)^3 = 3.3333...

Side length of the cubic unit cell = 1.49

Task 3

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?

Face-centred cubic lattice contains 4 atoms and box contains 1000 unit cells of the lattice.

Total number of atoms = atoms per unit cell of lattice x number of unit cells

Total number of atoms = 4 x 1000 = 4000

Task 4

Using the LAMMPS manual, find the purpose of the given commands in the input script.

mass 1 1.0

This is known as the mass command, which sets the mass for all atoms of one or more atom types. '1' refers to the atom type and '1.0' refers to the mass value.

pair_style lj/cut 3.0

This indicates that the only Lennard-Jones interaction to be considered is those between two atoms with distance less or equal to 3 units apart. Interactions with neighboring atoms that position beyond distance of 3 units will be ignored and thus the potential energy will be set to zero.

pair_coeff * * 1.0 1.0

This is the pair_coeff command. The asterisks represent the atom types which can be anything from 1 to N. '1.0' and '1.0' is the coefficient value, corresponding to $\sigma$ and $ϵ$ involved in Lennard-Jones potential. $\sigma$ refers to the distance when the potential energy between atoms is zero while $ϵ$ represent the well depth calculated using Lennard-Jones interaction formula.

Task 5

Given that we are specifying Xi(0) and Vi(0), which integration algorithm are we going to use?

Velocity Verlet algorithm should be used. Verlet algorithm does not involve the calculation of the velocities.

Task 6

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 time used to run each simulation is kept constant. However, the number of steps used to run each simulation can be varied. The smaller the timestep is, the greater the number of steps used to run each simulation. The first set of codes set the timestep used to run each simulation. It allows timesteps to be modified instantaneously if a different value of timestep is being chosen. This is done by changing the numerical value stated after 'equal' in the first line of codes, e.g. 0.001 in this case, which would then automatically modify timestep in the later script. However, this cannot be done as easily using second set of codes. Second set of codes involves changing all lines related to timestep manually if timestep is to be changed. Hence, the first set of code is preferred. It allows modifications to be made more easily and reduces the risk of making any mistakes in the coding.

NJ: I think you're dancing around the answer a bit, but I think you've got the idea behind it - this section means that you always simulate the same time, not matter what timestep you choose.

Task 7

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?

Plot	Description
x200px200px	Plot illustrating variation in total energy against time, timestep 0.001.
x200px200px	Plot illustrating variation in temperature against time, timestep 0.001.
x200px200px	Plot illustrating variation in pressure against time, timestep 0.001.
x200px200px	Plot illustrating how total energy varies against time in all cases of timestep simulated.

Simulation run with timestep 0.001 quickly reaches equilibrium after about 0.3 time units, when the graph appears to oscillate about an average value. Total energy is conserved in the real physical system. Of the five timesteps used, timesteps of 0.0025 and 0.001 give energy with the smallest fluctuation over the simulation time. Hence, simulation run with 0.0025 and 0.001 resemble real physical system most closely in comparison to other timesteps. Since using short timesteps means that the same number of simulation steps would cover a shorter amount of simulation time, which makes simulation with long observation time undesirable. Therefore, timestep of 0.0025 is the largest one to give acceptable result. Timestep of 0.015 is a particular bad choice because equilibrium is never reached - total energy increases over the time of simulation and energy is not conserved.

NJ: Excellent

Running simulations under specific conditions

Task 1

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.

Temperature chosen are: 1.6, 1.8, 2.0, 2.2 and 2.4. Pressure chosen are 2.5 and 2.7. Timestep chosen is 0.0025 since it gives the smallest fluctuation in energy as discussed in the section above.

Task 2

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

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

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

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

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

Task 3

Use the manual page to find out the importance of the three numbers 100 1000 100000. How often will values of the temperature, etc., be sampled for the average? How many measurements contribute to the average? Looking to the following line, how much time will you simulate?

100 indicates that every 100th step during a simulation is being recorded. 100000 indicates that the average value is being calculated from the previous data recorded within the total steps of 100000. 1000 means 1000 previous values is used to take final average during the 100000 steps. 1000 x 100 = 100000, this means the simulation has been run on all data in this case.

Task 4

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?

Plot	Description
x200px200px	Plot illustrating how density varies with temperature at different pressure.

Density predicted by Ideal Gas Law can be worked out as follows:

$pV=N{K}_{B}T$

Since $K_B$ in reduced unit equals to 1,

$pV=NT$

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

Simulated density is lower than that predicted by Ideal Gas Law. Ideal Gas Law assumes there is no interaction between gas molecules. Hence, ideal gas molecules can approach each other with no repulsive force in between. As the result, gas molecules are closer to each other and more of them can be found in a given volume, giving greater a density. On the other hand, stimulation takes interaction between gas molecules into account. Gas molecules cannot approach as closely to each other in this case, due to the repulsion they experience as they get closer. This means fewer gas molecules can be found in a given volume and density is lower than that predicted using Ideal Gas Law.

The plot shows that as temperature increases, density decreases. This can be justified. As temperature increases, gas molecules gain a higher kinetic energy. As their kinetic energy increases, they move away from each other. Subsequently few molecules can be found in a given volume and density decreases.

At higher pressure, molecules experience a greater repulsive force, lowering the density. As discussed before, Ideal Gas Law assumes no interaction between gas molecules and so they can approach much closer to each other. Therefore, discrepancy between Ideal and stimulated cases increases with pressure.

NJ: Good. How does the discrepancy change with temperature?

Calculating heat capacities using statistical physics

Task

As in the last section, you need to run simulations at ten phase points. In this section, we will be in density-temperature $({\rho }^{*},T^{*})$ phase space, rather than pressure-temperature phase space. The two densities required at 0.2 and 0.8, and the temperature range is 2.0, 2.2, 2.4, 2.6, 2.8. Plot $C_V/V$ as a function of temperature, where V is the volume of the simulation cell, for both of your densities (on the same graph). Is the trend the one you would expect? Attach an example of one of your input scripts to your report.

Plot	Description
x200px200px	Plot illustrating how Cv/V (heat capacity/volume) varies with temperature at different densities.

Cv/V gives the heat energy required to increase the temperature of a substance by 1K per unit volume. Heat capacity is an extensive property and thus its magnitude depends on the size of the system. At a greater density, there are more molecules per unit volume and magnitude of Cv/V (heat capacity per unit volume) is expected to be higher. This is consistent with trend shown by the graph where heat capacity at a given temperature appears greater for the system with higher density. In addition, molecules are closer together at a higher density. This means the potential energy between molecules is higher. Thermal energy is the sum of potential energies and kinetic energies. Higher potential energy gives system a higher internal thermal energy, leading to a higher heat capacity of the system at a higher density.

NJ: Not necessarily, but it does mean that there are more degrees of freedom per unit volume.

The graph shows Cv/V decreases with increasing temperature, which is not expected. This is because as temperature increases, the degrees of freedom within a system increases, which results in a higher internal energy within the system. Higher internal energy makes it hard to raise temperature of substance, increasing heat capacity.

NJ: This is certainly true at very low temperatures - at higher temperatures, once you can access all modes in your system (translational, rotational, vibrational, ...), then you can observe this kind of behaviour. Can you rationalise it?

Script for density = 0.2 and T =2.2:

### 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 100000
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 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_dens2 v_temp2 v_press2 v_etotal 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[7]
variable aveetotal2 equal f_aves[8]
variable heatcapacity equal ${N2}*(${aveetotal2}-${aveetotal}*${aveetotal})/(${avetemp}*${avetemp})

print "Averages"
print "--------"
print "Density: ${avedens}"
print "Temperature: ${avetemp}"
print "pressure: ${avepress}"
print "total energy: ${aveetotal}"
print "total energy: ${aveetotal2}"
print "heatcapacity: ${heatcapacity}"
print "number of atoms: ${N2}"
print "volume: ${volume}"

Structural properties and the radial distribution function

Task

perform simulations of the Lennard-Jones system in the three phases. When each is complete, download the trajectory and calculate $g(r)$ and $\int g(r)\mathrm{d}r$. Plot the RDFs for the three systems on the same axes, and attach a copy to your report. Discuss qualitatively the differences between the three RDFs, and what this tells you about the structure of the system in each phase. In the solid case, illustrate which lattice sites the first three peaks correspond to. What is the lattice spacing? What is the coordination number for each of the first three peaks?

Plot	Description
x200px200px	Plot illustrating variation in g(r) with distance for solid, liquid and gas.
x200px200px	Plot illustrating variation in g(r) integral with distance for solid, liquid and gas.

Radial Distribution Function, g(r) describes how density of particles varies as a function of distance from the reference particle.

The plot of g(r) against distance shows many peaks for solid. Atoms can be found at the distances where graph peaks. Solid has a long range order where atoms within solid have an ordered structure, packed closely together. Hence atoms can be found more frequently within a given distance, resulting in many peaks. in comparison to liquid and gas.

The plot of g(r) against distance for liquid shows fewer peaks since atoms in liquid do not pack as closely together. For this reason, the densities, g(r), are lower than those found in solid. These peaks locate at distance shorter than 4 units, after which g(r) approaches 1. NJ: g(r) is not the density, it is the probability of observing a neighbouring atom.

Only one significant peak can be found within the distances being examined and it locates at distance shorter than 2 units, after which g(r) approaches 1. This shows atoms in gas are furthest apart from each other, giving the least ordered structure.

As distance from the reference particles increases, the density of atoms decreases - decreasing g(r) shows few atoms found. When g(r) approaches 1, it means there is no significant atoms found and there is no predictable arrangement of atoms at the given distance.

NJ: These two statements aren't the same - it's true that there is no order when g(r)=1, but that isn't the same as there being no atoms. You will observe a number of atoms consistent with the bulk density.

In solid, the three peaks indicate the position of atoms that are closest to the reference atom. Solid takes a face-centred cubic structure:

x200px200px	Diagram showing atoms nearest to the reference atom in solid that adopts face-centred cubic

Each side of face-centred cubic adopts a unit length of A. From the plot of g(r) against distance for solid, it is found that three nearest atoms locate at distance of 1.025, 1.525 and 1.825 respectively.

For the most nearest atom:

$1.025=A\sqrt{\frac{1}{2}}$

$A=1.450$

For the second nearest atom:

$1.525=A$

$A=1.525$

For the third nearest atom:

$1.825=A\sqrt{\frac{3}{2}}$

$A=1.490$

Taking average of A gives:

$\frac{1.450+1.525+1.490}{3}=1.49$

Hence, average lattice spacing equals to 1.49.

Plot of g(r) integral against distance gives coordination number at given distance.

Coordination number of first peak = 12

Coordination number of second peak = 18-12 = 6

Coordination number of third peak = 42-18 = 24

NJ: Excellent. And are these peaks consistent with what you would expect for an FCC lattice?

Dynamical properties and the diffusion coefficient

Task 1

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.

Plot	Description
x200px200px	Plot showing Mean Squared Displacement against Time for Solid
x200px200px	Plot showing Mean Squared Displacement against Time for Liquid
x200px200px	Plot showing mean Squared Displacement against Time for Gas

The trends illustrated by graphs are expected. The greater the gradient of the plot, the greater the diffusion coefficient. MSD against time plot for solid shows small fluctuations in MSD over time. This is because atoms in solid are closely packed - they are not able to move freely. On the other hand, the plot of MSD against time for liquid increases almost linearly with time, indicating the greater freedom of movement of atoms in liquid. Atoms in gas phase move most freely, moving as far apart from each other as possible over time, and this is indicated by the plot of MSD against time for gas that increases most rapidly.

NJ: Why is the gas curve not a straight line, do you think?

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

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

Phases	Gradient (Stimulated)	D (Stimulated)	Gradient (1 million atoms)	D (1 million atoms)
Solid	0.0002	3.33 x 10-5	0.00003	5 x 10-6
Liquid	1.0184	1.70 x 10-1	0.5236	8.73 x 10-2
Gas	11.389	1.90	15.249	2.54

Task 2

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.

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

Using $\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 this is an odd function, $\sin (w\tau )$ converges to zero and numerator of the second term becomes zero.

Hence,

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

Task 3

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.

x200px200px

Plot shows VACF against Timesteps for liquid and solid as well as $C(\tau )$ against Timesteps when $\omega =1/2\pi$

Minima and maxima in VACF plots for liquid and solid system represent the points at which the velocity changes the direction. In solid, this is due to the vibrations of atoms in a confined structure whereas in liquid, it is due to the collisions between the moving atoms.

C(τ) of solid system converges to zero over a longer period of time, in comparison to liquid system. This is because atoms within solid system have an ordered structure where atoms are tightly bonded to each other. It is harder for these tightly bonded atoms to become uncorrelated in velocity after the small vibrations they experience when they bump onto each other. However, atoms in liquid move relatively freely in a disordered system. It takes less time for liquid atoms to become uncorrelated in velocity since their velocities change by a greater amount when they collide onto each other due to the disordered nature of the system. Unlike solid, atoms in liquid are able to move away randomly after collision.

Harmonic oscillator is very different to the Lennard Jones Solid and Liquid. This is because system is modelled as a single-atom system using harmonic oscillator. Atoms in such system are assumed to move forward and backward only, without physically colliding with others and hence there is no transfer of energy. Harmonic oscillator is modelled using a Cosine Graph where magnitude of atomic velocity is fixed but the direction of the atom's motion reverses at a fixed time interval. Since there is no transfer of energy, atoms move in a forward and backward fashion for an infinite time steps. However, Lennard Jones considers interactions between two atoms. This means their velocity autocorrelation function goes to zero as time steps increases.

NJ: Good

Task 4

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?

Plots	Diagrams
Running Integral for Solid	x200px200px
Running Integral for Solid (1 million atoms)	x200px200px
Running Integral for Liquid	x200px200px
Running Integral for Liquid (1 million atoms)	x200px200px
Running Integral for Gas	x200px200px
Running Integral for Gas (1 million atoms)	x200px200px

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

$D=\frac{1}{3}$ x Integral of VACF

System	Integral (timesteps)	Integral (time)	D
Solid	-3.38	-0.00676	-0.00225
Liquid	477	0.954	0.318
Gas	3090	6.18	2.06
Solid (1 million atoms)	-2.65	-0.00531	-0.00177
Liquid (1 million atoms)	246	0.491	0.164
Gas (1 million atoms)	2930	5.86	1.95

Trends shown by results are expected. The value of diffusion coefficients increases as we move from a solid system to a gas system. This is because the systems become more disordered in sequence. Atoms in gas move quickly in such disordered system, resulting in high D value.

On the other hand, atoms in solid are tightly bonded together in a structured system. They are not free to move – they only vibrate. Hence, diffusion coefficient is the lowest for the solid system.

Values of D for each type of phases are very similar in both cases of simulation using MSD or VACF. They are expected to stay consistent since each phase system shows a distinct diffusion property. NJ: I see what you're getting at, but this sentence doesn't really make sense. I think you mean to say something like "They are expected to be consistent because the diffusion coefficient should depend only on the thermodynamic state, not on the method used to calculate it." Variation in the values is probably due to difference in simulation time period used as well as the number of atoms involved in the calculation.

One of the biggest sources of errors in VACF is due to the short simulation time. Simulation time used is particularly important for gas phase. In all examined cases of VACF, only 500 time steps have been considered in the calculation involving D. Atoms in each phase are expected to behave differently after some times. Atoms in gas system have the most freedom of movement and hence gas atoms are expected to behave most differently as time proceeds. This means inadequate amount of timesteps used will fail to examine the complete system (e.g. of gas) and will affect the accuracy of diffusion coefficient calculated. Also the number of atoms set for calculation can be a potential problem that lead to variation in magnitude of D values. Usually the greater the amount of atoms being modelled (VACF), the more accurate the results are. NJ: You're right about the time, but not about why. The more time you simulate, the more accurate the VACF will be. You know that it should decay to 0 at long times - in a finite simulation, it will never truly get there. The longer you simulate, the closer it will get.