Talk:Mod:sh5214 liquid simulation

JC: General comments: All tasks attempted, but a few mistakes. Add a bit more detail to your written answers to show that you have a good overall understanding of the experiment and show your working for calculations.

TASK: Open the file HO.xls. In it, the velocity-Verlet algorithm is used to model the behaviour of a classical harmonic oscillator. Complete the three columns "ANALYTICAL", "ERROR", and "ENERGY": "ANALYTICAL" should contain the value of the classical solution for the position at time $t$ , "ERROR" should contain the absolute difference between "ANALYTICAL" and the velocity-Verlet solution (i.e. ERROR should always be positive -- make sure you leave the half step rows blank!), and "ENERGY" should contain the total energy of the oscillator for the velocity-Verlet solution. Remember that the position of a classical harmonic oscillator is given by $x (t) = A \cos (ω t + ϕ)$ (the values of $A$ , $ω$ , and $ϕ$ are worked out for you in the sheet).

TASK: For the default timestep value, 0.1, estimate the positions of the maxima in the ERROR column as a function of time. Make a plot showing these values as a function of time, and fit an appropriate function to the data.

TASK: Experiment with different values of the timestep. What sort of a timestep do you need to use to ensure that the total energy does not change by more than 1% over the course of your "simulation"? Why do you think it is important to monitor the total energy of a physical system when modelling its behaviour numerically?

The maximum timestep that doesn't result in more than an 1% change in energy is around 0.45. The reason why the energy is monitored is to make sure the energy change between each timestep is not too large. If there is a large energy difference, it would mean that the overall change in velocity of the particles in the system is too large, which would lead to not obtaining the full amount of information from the simulation.

JC: Can you be more specific, energy should be conserved so we need to check that energy is approximately constant to ensure that the simulation results are physical.

TASK: For a single Lennard-Jones interaction, $ϕ (r) = 4 ϵ (\frac{σ^{12}}{r^{12}} - \frac{σ^{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_{e q}$ , and work out the well depth ( $ϕ (r_{e q})$ ). Evaluate the integrals $\int_{2 σ}^{\infty} ϕ (r) d r$ , $\int_{2.5 σ}^{\infty} ϕ (r) d r$ , and $\int_{3 σ}^{\infty} ϕ (r) d r$ when $σ = ϵ = 1.0$ .

The value of $r_{0}$ , for which potential energy is 0, is $r_{0} = σ$ .

The force at this seperation is $F = \frac{24 ϵ}{σ}$ .

The equilibrium separation is $r_{e q} = \sqrt[6]{2} σ$ .

The well depth is $ϕ (r_{e q}) = - 2 ϵ$

$\int_{2 σ}^{\infty} ϕ (r) d r = - 0.0248 (3 s f)$

$\int_{2.5 σ}^{\infty} ϕ (r) d r = - 0.00818 (3 s f)$

$\int_{3 σ}^{\infty} ϕ (r) d r = - 0.00329 (3 s f)$

JC: The well depth is epsilon, not 2*epsilon. Otherwise the maths is correct, but show your working.

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

Number of water molecules in 1mL water = 0.0556 moles (3sf) = 3.35e22 molecules (3sf)

Volume of 10000 water molecules = 2.99e-19 mL (3sf) = cube with side length of ~ 6.69 nm (3sf)

JC: Correct, but show some working.

TASK: Consider an atom at position $(0.5, 0.5, 0.5)$ in a cubic simulation box which runs from $(0, 0, 0)$ to $(1, 1, 1)$ . In a single timestep, it moves along the vector $(0.7, 0.6, 0.2)$ . At what point does it end up, after the periodic boundary conditions have been applied?.

After PBC has been applied the atom will be at (0.2,0.1,0.7)

JC: Correct.

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

$r = r^{*} \times σ = 3.2 \times 0.34 = 1.088 n m$

$T = T^{*} \frac{ϵ}{k_{B}} = 1.5 \times 120 = 180 K$

$ϕ (r) = - 6.16 \times 1 0^{- 24} J m o l^{- 1}$

JC: r and T are correct, but phi is not - again show some working.

TASK: Why do you think giving atoms random starting coordinates causes problems in simulations? Hint: what happens if two atoms happen to be generated close together?

If atoms are created too close, then we will have a very large potential, meaning that if the timestep is not reduced, then we will have a large change in energy between timesteps which could mean that we do not obtain all the information from the simulation system. However if we decrease the timestep so that the change in energy is reduced to a suitable size, then it would require much more calculations to run the simulation for the same of time.

JC: Good.

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

Simple cubic: 0.8 density = 1.25 cell volume = $\sqrt[3]{1.25}$ = 1.07722 cell length.

Face centre cubic: 1.2 density = 0.833 cell volume = $\sqrt[3]{0.833}$ = 0.941 cell length.

JC: An FCC lattice has 4 atoms per unit cell, so the side length should be (4/1.2)^(1/3) = 1.49.

TASK: Consider again the face-centred cubic lattice from the previous task. How many atoms would be created by the create_atoms command if you had defined that lattice instead?

FCC contains 4 atoms in each lattice cell, SC contains 1 atom in each lattice cell, thus 4000 atoms would be created if FCC was used.

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

mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0

mass 1 1.0: mass of atom type 1 is 1.0

pair_style lj/cut 3.0: cut-off for lennard jones interaction is 3 unit length

pair_coeff * * 1.0 1.0: pair_coeff refers to the force field coeffecients for pairs of atoms, * * means between any two atom types, 1.0 is the coefficient

JC: What are the force field coefficients for a Lennard Jones simulation?

TASK: Given that we are specifying $x_{i} (0)$ and $v_{i} (0)$ , which integration algorithm are we going to use?

Velocity velvet integration

TASK: 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 reason why a variable is used to hold the value of timestep is so that it is more convenient to change the value of timestep, as you would only have to change the line that defines the timestep instead of going through the whole code and changing all instances where timestep is used.

JC: Good explanation.

TASK: make plots of the energy, temperature, and pressure, against time for the 0.001 timestep experiment (attach a picture to your report). Does the simulation reach equilibrium? How long does this take? When you have done this, make a single plot which shows the energy versus time for all of the timesteps (again, attach a picture to your report). Choosing a timestep is a balancing act: the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality; short timesteps, however, mean that the same number of simulation steps cover a shorter amount of actual time, and this is very unhelpful if the process you want to study requires observation over a long time. Of the five timesteps that you used, which is the largest to give acceptable results? Which one of the five is a particularly bad choice? Why?

From the plot of energy vs time for timestep = 0.001 we can see that equilibrium is reached, and from the last plot it can be seen that it takes about 0.4 units of time to reach equilibrium. The largest timestep that results in the equilibration of energy is timestep = 0.01. Timestep = 0.015 is a bad choice because the system does not equilibrate and instead diverges.

JC: Which timestep did you choose? Average total energy shouldn't depend on choice of timestep so 0.01 and 0.0075 are not suitable. 0.0025 and 0.001 have the same average energy so choose the larger timestep so that simulations will cover more time.

TASK: We need to choose $γ$ so that the temperature is correct $T = 𝔗$ if we multiply every velocity $γ$ . 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} {(γ v_{i})}^{2} = \frac{3}{2} N k_{B} 𝔗$

Solve these to determine $γ$ .

$\frac{1}{2} \sum_{i} m_{i} {(γ v_{i})}^{2} = \frac{3}{2} N k_{B}$

$⟹ \frac{1}{2} \sum_{i} m_{i} v_{i}^{2} = \frac{3 N k_{B} 𝔗}{2 γ^{2}} = \frac{3}{2} N k_{B} T$

$⟹ γ^{2} = \frac{𝔗}{T} ⟹ γ = \sqrt{\frac{𝔗}{T}}$

JC: Correct.

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

100: number of timesteps between each sample

1000: number of samples for each average

100000: number of timesteps required to obtain one average

TASK: When your simulations have finished, download the log files as before. At the end of the log file, LAMMPS will output the values and errors for the pressure, temperature, and density $(\frac{N}{V})$ . 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?

My simulated density is lower than the density predicted by the idea gas law, this discrepancy increases with pressure. The reason why the deal gas law gives a higher density is due to the fact that it doesn't take into account the interactions between atoms, hence repulsions between atoms are not accounted for in the ideal gas law, this means that the atoms would be closer together in the idea gas law and so a higher density.

JC: Good, how does the discrepancy change with temperature?

TASK: As in the last section, you need to run simulations at ten phase points. In this section, we will be in density-temperature $(ρ^{*}, 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.

From the graph it can be seen that specific heat capacity decreases with temperature, this is expected because there is a limit to how much energy a molecule can hold, as the temperature of the molecule increases, the amount of available rotational/vibrational energy level that the molecule can be excited into decreases, meaning that the molecule will be able to store less energy as its temperature increases.

We can also see that the specific heat capacity for the system with 0.8 density is higher than the system with 0.2 density, this is expected because a higher density means that there are more molecules per unit volume, and therefore more energy could be stored per unit volume.

JC: You are just simulating spherical particles, not molecules, so there is no rotational or vibrational energy. Also remember that your simulations are classical.

### DEFINE SIMULATION BOX GEOMETRY ###
variable D equal 0.2
lattice sc ${D}
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 vdamp equal ${timestep}*1000
fix nvt all nvt temp ${T} ${T} ${tdamp}
run 10000
reset_timestep 0

## SWITCH OFF THERMOSTAT ##	
unfix nvt
fix nve all nve

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp press density atoms vol
variable vol equal vol
variable N2 equal atoms*atoms
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 etotal equal etotal
variable etotal2 equal etotal*etotal
fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_etotal v_dens2 v_temp2 v_press2 v_etotal2 v_vol 
run 100000

variable avetemp equal f_aves[2]
variable errtemp equal sqrt(f_aves[6]-f_aves[2]*f_aves[2])
variable cv equal (${N2}*(f_aves[8]-f_aves[4]*f_aves[4])/(${T}*${T}))
variable V equal f_aves[9]

print "Averages"
print "--------"
print "Temperature: ${avetemp}"
print "Stderr: ${errtemp}"
print "cv: ${cv}"
print "volume: ${V}"

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

Firstly the range of peaks for solids is longer than liquids which is longer than gas, this means that in the system of solids, the atoms are more ordered than the atoms in the system of liquids, which are more ordered than the atoms in the system of gases. Another observation that could be made is that the magnitude of the peaks decreases as r increases, this is because g(r) is calculated by number of atoms in the shell divided by volume of shell and the density of the system, as r increases, the volume of the shell will increase at a faster rate than the number of atoms in the shell, and thus g(r) decreases. Also it can be seen that magnitude of the peaks of the RDF of the solid is greater than both liquid and gas, this is because solids are more closely packed and more ordered than both gas and liquid, hence at certain r values there will be a high number of atoms and at other r values there will be a very low number of atoms, and that's why solids have higher peaks and lower troughs.

The first peak for the solid RDF is at r=1.025 and the second peak is at r=1.425, thus the lattice spacing is 0.4 times the interatomic separation. There are three different planes in a fcc crystal, they have miller indices of (1,1,1) (1,1,0) and (1,0,0). The (1,1,1) plane corresponds to the first peak and has a coordination number of 6 atoms, the (1,1,0) plane corresponds to the second peak and has a coordination number of 2, the (1,0,0) plane corresponds to the third peak and has a coordination number of 4. The coordination number of lattice sites will affect the magnitude of the peak, and that's why magnitude of first peak > third peak > second peak

JC: Did you plot the integral of g(r)? This will tell you the area under each peak which is the coordination number for that peak, the coordination numbers should be 12 for the first peak, 6 for the second and 24 for the third. Drawing a diagram of an fcc lattice to show which atoms are responsible for each peak might help. The lattice spacing is the width of the unit cell. You can calculate this from each of the first three peaks separately and then calculate the average.

TASK: make a plot for each of your simulations (solid, liquid, and gas), showing the mean squared displacement (the "total" MSD) as a function of timestep. Are these as you would expect? Estimate $D$ in each case. Be careful with the units! Repeat this procedure for the MSD data that you were given from the one million atom simulations.

Gas MSD (my data):

t=2000, MSD=29.2
t=5000, MSD=104.7
gradient = (104.7-29.2)/3000 = 0.0252 (3sf)
D = gradient/6 = 0.00419 (3sf)

Gas MSD (1e6 atoms):

t=2000, MSD=36.4
t=5000, MSD=144.4
gradient = (144.4-36.4)/3000 = 0.036
D = gradient/6 = 0.006 (3sf)

Liquid MSD (my data):

t=0, MSD=0
t=5000, MSD=11
gradient = 11/5000 = 0.0022
D = gradient/6 = 0.000367 (3sf)

Liquid MSD (1e6 atoms):

t=0, MSD=0
t=5000, MSD=5.19
gradient = 5.19/5000 = 0.00104 (3sf)
D = gradient/6 = 0.000173

Solid MSD (my data): D=0

Solid MSD (1e6 atoms): D=0

MSD is the average distance traveled by a particle in a system, it can be seen that the MSD of the solid plateaus very early on and at a very low value this is due to the fact that atoms in a solid are rigid and fixed in place therefore average distance traveled by a atom in a solid is very small. We see that the MSD plot of the gas starts off as a parabola before turning linear, the parabola part corresponds to when gas atom doesn't collide with anything else meaning that during this time the velocity of the atom is constant and so distance increases proportional to time squared, when more and more collisions occur the curve goes from parabolic to linear. In the liquid MSD plot, we can see that at the start of the curve for a short amount of time the curve is not linear, this is because little to no collisions occur right at the start of the simulation, for most of the plot the line is linear due to the fact that there are collisions between liquid atoms.

The magnitude of MSD for gas > liquid > solid, this is due to the fact that density of gas < liquid < solid, a lower density means that less collisions occur and thus the average translational velocity of atoms is higher, resulting in greater distance traveled.

JC: Good explanation of the curve in the gas MSD. Show the lines of best fit that you used to calculate D on your graphs.

TASK: In the theoretical section at the beginning, the equation for the evolution of the position of a 1D harmonic oscillator as a function of time was given. Using this, evaluate the normalised velocity autocorrelation function for a 1D harmonic oscillator (it is analytic!):

$C (τ) = \frac{\int_{- \infty}^{\infty} v (t) v (t + τ) d t}{\int_{- \infty}^{\infty} v^{2} (t) d t}$

Be sure to show your working in your writeup. On the same graph, with x range 0 to 500, plot $C (τ)$ with $ω = 1 / 2 π$ 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.

$x (t) = A c o s (ω t + ϕ)$

$v (t) = \frac{d x}{d t} = - A ω s i n (ω t + ϕ)$

$C (τ) = \frac{\int_{- \infty}^{\infty} - A ω s i n (ω t + ϕ) . - A ω s i n (ω t + ω τ + ϕ) d t}{\int_{- \infty}^{\infty} A^{2} ω^{2} s i n^{2} (ω t + ϕ) d t}$

$= \frac{\int_{- \infty}^{\infty} 0.5 A^{2} ω^{2} [c o s (ω τ) - c o s (2 ω t + ω τ + 2 ϕ)] d t}{\int_{- \infty}^{\infty} 0.5 A^{2} ω^{2} - 0.5 A^{2} ω^{2} c o s (2 ω t + 2 ϕ)}$

$= \frac{[0.5 A^{2} ω^{2} c o s (ω τ) t + 0.25 A^{2} ω s i n (2 ω t + ω τ + 2 ϕ)]_{- \infty}^{\infty}}{[0.5 A^{2} ω^{2} t + 0.25 A^{2} ω c o s (2 ω t + 2 ϕ)]_{- \infty}^{\infty}}$

We know that sine and cosine has a range of -1 to 1, hense the above equation can be simplified to:

$\approx \frac{[0.5 A^{2} ω^{2} c o s (ω τ) t]_{- \infty}^{\infty}}{[0.5 A^{2} ω^{2} t]_{- \infty}^{\infty}}$

$\approx c o s (ω τ)$

JC: Correct final result, but you need to check the derivation again. You don't have to calculate any of the integrals if you use a trigonometric identity for sin(A+B) to split up one of the sin terms and then realise that the integral of an odd function over limits which are symmetric about zero is zero.

The minima in the VACFs for the liquid and solid system represents the time at which the velocity of the atoms in the system is at a maximum. We can see that the VACF for solids oscillates about zero, with the amplitude of the oscillation will decrease over time, this is due to the fact that the atoms in the solids are locked up in a lattice, these atoms in the lattice will want to be at the optimal position where the repulsive forces and attractive forces balance out, they will therefore oscillate about this position which will also result in an oscillation in velocity, the reason why the magnitude of these oscillation decreases over time is because of other forces that perturb the oscillation. We can see that there is also a very slight minima on the liquid VACF curve, which shows that liquids also oscillate but the oscillation disappears very quickly, the reason why atoms in liquids don't oscillate for as long as solids is because liquid atoms are not locked in a lattice like solids and so are free to move, this means that any oscillation is very quickly dampened out by diffusion.

The reason why the harmonic oscillator VACF is so different is due to the fact that it doesn't take into account dampening of the oscillation, hence the magnitude of the oscillation never decreases.

JC: The harmonic oscillator doesn't decay because there are no collisions. In the solid and liquid VAFCs, successive collisions randomise the particle velocities and so, since the VACF is calculated as an average over all particles, it decays to zero.

TASK: Use the trapezium rule to approximate the integral under the velocity autocorrelation function for the solid, liquid, and gas, and use these values to estimate $D$ in each case. You should make a plot of the running integral in each case. Are they as you expect? Repeat this procedure for the VACF data that you were given from the one million atom simulations. What do you think is the largest source of error in your estimates of D from the VACF?

Gas VACF (my data):

t=5000, Running Integral=6.54 (3sf)
D = Running Integral/3 = 2.18 (3sf)

Gas VACF (1e6 atoms):

t=5000, Running Integral=9.81 (3sf)
D = Running Integral/3 = 3.27 (3sf)

Liquid VACF (my data):

t=5000, Running Integral=0.569 (3sf)
D = Running Integral/3 = 0.190 (3sf)

Liquid VACF (1e6 atoms):

t=5000, Running Integral=0.270 (3sf)
D = Running Integral/3 = 0.0901 (3sf)

Solid VACF (my data): D = 0

Solid VACF (1e6 atoms): D = 0

JC: How does this method compare to the MSD and what was the main source of error?