Talk:Mod:Az1114

JC: General comments: All tasks attempted, but some details missing especially in the last few sections. Try to add a bit more to the explanations of your results and make sure that you understand the background theory behind the tasks and use this in your explanations.

Velocity-Verlot Algorithm versus the Harmonic Oscillator

Below in the two figures are the graphs for displacement and energy versus time of a classical harmonic oscillator compared to the values given by the velocity-Verlet algorithm. Both graphs are set at a timestep of 0.1, and it can be seen that for the displacement there is no significant deviation between the displacement of the two different methods.

The following graph is the error of the velocity-Verlot alorithm, or more specifically the absolute difference between the two methods versus time. Note that the value of error oscillates with time between 0 and a maximum, with the maximum generally increasing. The error maximum is found when displacement is zero, thus when at a velocity maximum.

A range of timesteps were used to find what is optimal, from 0.001 to 0.015. It was found that with higher timestep the error with displacement (and therefore energy) increased too. The total error of 1% was not crossed with the timesteps used, but it can be extrapolated to be around a timestep of 0.2.

JC: Good choice of timestep, why does the error oscillate and the maximum error increase over time?

Lennard-Jones Potential

The below equation calculates the Leonard-Jones potential:

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

By equating this to 0, we find that the $r_0$ - the separation at which the potential is 0 - is $r_0=\sigma$.

We can then find the force at this separation by differentiating the Lennard-Jones potential with respect to r, and substituting in the value $r_0$:

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

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

$r_0=\sigma$

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

The equilibrium separation ($r_{eq}$) is the separation of the particles at equilibrium, and represents the distance that represents the minimum in potential energy. It can therefore be found by equating the force to zero $\frac{\mathrm{d}\varphi }{\mathrm{d}r}$, (or $\mathbf{F}$) to 0.

Hence,

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

$r_{eq}=\sqrt[6]{2{\sigma }^6}$

Substituting $r_{eq}$ back into $\varphi (r)$ gives us the well depth:

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

---

To put this equation into perspective, let us calculate the L-J potential for some realistic cutoff values:

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

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

When $\sigma =ϵ=1.0:$

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

$=\underset{r\to \mathrm{\infty }}{\lim }4(\frac{1^{12}}{11r^{11}}-\frac{1^6}{5r^5})-4(\frac{1^{12}}{11\times 2^{11}}+\frac{1^6}{5\times 2^5})$

$=0-0.024822\cdots$

$\approx -0.0248$

Following the same calculations:

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

$=-0.0081767\cdots$

$\approx -0.00818$

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

$=-0.0032901\cdots$

$\approx -0.00329$

The calculations show that as the cutoff distance increases, the Lennard-Jones potential decreases, which agrees with the theory.

JC: All maths correct, except the force should be positive.

Periodic Boundary Conditions

It is very difficult to simulate realistic volumes of liquid and in the scope of this experiment, we have simulated liquids with N (number of atoms) between 1,000 and 10,000. To demonstrate why simulating realistic volumes of liquid is difficult, the number of molecules in 1 mL of water is given by:

(Assuming: Standard conditions; Density of water = 1.00 g/ml)

$1.00\mathrm{m}\mathrm{l}=1.00\mathrm{g}$

$18.02\mathrm{g}/\mathrm{m}\mathrm{o}\mathrm{l}=0.0555\text{ mols of water}$

$0.0555\times 6.022\times 10^{23}=3.34\times 10^{22}\text{ molecules of water.}$

Conversely, 10,000 water molecules is only $2.99\times 10^{-19}\text{ ml of water}$

For the simulation to be completed in a reasonable amount of time, the atoms in the simulation are placed in a box with boundary conditions. The box is repeated in all directions, similar to a unit cell in a lattice. If an atom moves outside of the box, another will enter the box from the opposite side, maintaining the same number of atoms in the box.

For example, if an atom at position $(0.5,0.5,0.5)$ in a cube of dimensions $(0,0,0)$ to $(1,1,1)$ moves along the vector $(0.7,0.6,0.2)$, its final position would be $(0.2,0.1,0.7)$, assuming the periodic boundary conditions have been upheld.

Reduced Units

The LAMMPS calculations run in this experiment are all done in reduced units

For example, the Lennard-Jones parameters for argon are $\sigma =0.34\mathrm{n}\mathrm{m}$; $ϵ/k_B=120K$.

If we set the cutoff to be $r*=3.2,$

$r=3.2\times 0.34=1.088\mathrm{n}\mathrm{m}$

Well depth: $ϵ=120K\times k_B\times 6.02\times 10^{23}$

$ϵ=0.998\mathrm{k}\mathrm{J}/\mathrm{m}\mathrm{o}\mathrm{l}$

and reduced temperature $T*=1.5$ would be $T=1.5\times 120=180\mathrm{K}.$

JC: All calculations correct and well laid out.

Equilibration

Creating the simulation box

When initiating the simulation, the atoms are given regular spacing intervals like that found in a square lattice. They are not given random positions, as if by chance two atoms occupy the same space, the program will encounter error.

For a simple cubic lattice (thus having 1 lattice point per unit cell) with number density 0.8:

$0.8\text{ points/volume}$

$\sqrt[3]{0.8}=0.9283\text{ points/length}$

$0.9283^{-1}=1.0772\text{ length/point}$

Therefore the spacing between atoms is 1.0772 unit lengths (in reduced units).

Giving the same treatment to a face-centred cubic lattice (4 lattice points/unit cell) with a number density of 1.2 will have a spacing of $0.5928$ as:

$\sqrt[3]{4.8^{-1}}=0.5928$

JC: This should be (4/1.2)^1/3 = 1.49.

In addition, the simple cubic lattice with a density of 0.8 would create 1,000 atoms in the given constraints of the box:

region box block 0 10 0 10 0 10

and a face centred cubic lattice with a density of 1.2 would create 4,000 atoms in the same box.

Setting the Properties of the Atoms

The input script contained a the following commands:

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

The mass command takes two arguments, setting the mass of type 1 atoms as $1.0$ (still in reduced units). The second line defines the type of interaction between the atoms (above set to 'lj' for Leonard-Jones); with a cutoff distance of <math?3.0</math> units. The third line defines the values for the Leonard-Jones interactions, the former $1.0$ setting the well depth $epsilon$ and the latter sets $sigma$. The two asterisks tell the program that this applies to all atoms in the system.

Once the properties of the particles are set, the initial positions and velocities are defined. As we have both the ${\mathbf{x}}_i(0)$ & ${\mathbf{v}}_i\left(0\right)$, the Velocity-Verlet integration half-step algorithm can be used.

JC: Good, why is a cutoff used?

Running the Simulation

In the following lines from the input script:

### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(100/${timestep})
timestep ${timestep}

### RUN SIMULATION ###
run ${n_steps}

We created and set a variable named 'timestep' to 0.001, and used that value to calculate how many runs (${n_steps}) to carry out. This is superior to writing

timestep 0.001
run 100000

This allows us only to change the value of the timestep variable to alter the number of runs, rather than having to calculate ${n_steps} each time the timestep is altered, especially in cases such as this computational experiment where the timestep can be changed multiple times.

JC: Correct.

Equilibration

The following 3 graphs check whether the system's total energy, temperature, and pressure (respectively), equilibrate over time during the 0.001 time-step simulation. y-axis values are in reduced units.

The following figure shows how total energy varies with time for many different timesteps. Note how at a large timestep of 0.15, the temperature actually diverges! Thus very high timesteps should not be used do to inaccuracies that can accrue.

JC: What about the other timesteps, which timestep did you choose and why? What is wrong with the 0.01 and 0.0075 - should the average total energy depend on the timestep?

Running Simulations Under Specific Conditions (NpT)

Temperature and Pressure Control

A timestep of 0.0025 was chosen as it is a good comprimise between accuracy of simulation and the length at which the simulation can run. The following Temperatures are pressures were chosen as were about the equilibration values: T: 1.5, 1.6, 1.7, 1.8, 1.9 P: 2.45, 2.65

In the simulation, the temperature was controlled deducing the (fluctuaing) instantaneous temperature $T,$ using the following equations:

$E_K=\frac{3}{2}N{k}_{B}T$

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

In order for $T$ to reach our target temperature $\mathfrak{𝔗},$ we introduce a constant $\gamma$ that controls the particles' velocities:

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

Thus we can solve for $\gamma .$:

From $(2)$:

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

Then substituting $\sum _i{m}_i{v}_i^2$ from $(1)$

${\gamma }^{2}3N{k}_{b}T=3N{k}_b\mathfrak{𝔗}$

$\gamma =\sqrt{\frac{ϵ}{T}}$

JC: Correct, but you have used epsilon rather than curly T in the final expression.

Explaining the Input Script

In order to calculate the average values for our simulation, the following parameters are set in the input script.

fix aves all ave/time Nevery Nrepeat Nfreq

This means that final averages are calculated only on timesteps that are multiples of Nevery, Nrepeat times; and an average is generated every Nfreq timesteps.

The following commands:

fix aves all ave/time 100 1000 100000
run 100000

Means that the results will be sampled every 100 timesteps for 1000 times (100, 200, 300, etc.). The averages are then calculated 100,000 timesteps, however as our simulation runs for the same amount of time, only a single average is produced with these parameters.

JC: Correct.

The timestep chosen for this simulation was 0.0025.

Equations of State

The above graph is a plot of the density versus temperature of the system for a liquid at two different pressures. There are plots for the values calculated from the simulation and from the ideal gas law.

There is a large difference between the simulated and ideal values, to an extent that they are much larger than the standard error of the simulated values, meaning the y-axis error bars are not visible. However they both follow a trend of the density decreasing with temperature, which follows the theory of thermal expansion.

A reason for the discrepancy will be due to the nature of the ideal gas law ignoring all interactions between particles, whereas in the simulation there is some Leonard Jones force, which when acting repulsively can push atoms further apart, resulting in a lower simulated density than in the ideal gas law.

Also note that the interactions in liquids are much greater than the interactions in gases, and so to use the ideal gas law for a liquid can quickly become a very crude approximation. It can also be seen that as the temperature increases, the discrepancy between the ideal and simulated densities becomes smaller.

JC: How does the discrepancy between the simulation and ideal gas results change with pressure? You could have chosen a larger difference in pressure to see the trend more clearly. Joining the ideal gas data points with straight lines is misleading as you know that the ideal gas law doesn't follow a straight line (density = 1/T).

Heat Capacity Calculations

The heat capacity to volume ratio versus temperature is plotted below. Theoretically, heat capacity will decrease with increasing temperature and increase with respect to increasing density, as there will be more molecules per unit volume that needs to increase in energy.

JC: Can you expand on this, why do you expect these trends theoretically?

Radial Distribution Functions

The radial distribution functions of a Leonard-Jones system in the solid, liquid and gas phase are distinctly different. This is due to the degrees of freedom that exist in each particular phase. In the solid phase, all the atoms are locked into a lattice structure, having only vibrational motion, with both long-range and short range order.

Liquids have more degrees of freedom, with no long-range order. They also exhibit rotational and to some extent translational motion. Then gases have the most freedom of all states, having a much larger translational motion.

From the RDF graph it can be seen that the solid function has a very distinct function with peaks and troughs. This is due to its atoms being fixed into a lattice with long-range order. For the liquid, there is a few peaks and troughs, representing the short range order of particles immediately surrounding the central particle, then averaging out to 1 - the system density. For the gaseous phase, the curve very quickly tends to the average value of 1, as there is no short- or long-range order in this phase.

JC: How can rotational motion alter the RDF considering that the simulations are of spherical particles?

From the curve for the solid the first 3 peaks at increasing distance can be attributed to the different lattice points in the fcc lattice about the original atom. For an atom on the corner of a cube, the 3 distances represent an atom at an adjacent vertex, an adjacent face, and an opposing face. These occur at a reduced distance of 1.06, 1.48 and 1.83 respectively. These agree with the distances noted on the graph. The coordination number can be deduced from the running integral of the curve and for the 3 peaks are given by 11, 6 & 24 respectively.

JC: Show the running integral graph and how you calculated at what distance the first 3 peaks should appear at. Can you calculate the lattice parameter from the RDF?

Dynamical Properties and the Diffusion Coefficient

The Diffusion Coefficients $D$ of a simulated gas, liquid and solid were calculated using several different methods:

1) From the Mean Squared Displacement (MSD) of 8,000 simulated atoms, according to the equation $D=\frac{1}{6}\frac{\partial ⟨r^2\left(t\right)⟩}{\partial t}$

2) From the MSD of a million atoms

3) From the Velocity Autocorrelation Function (VACF), by integration.

4) From the VACF of a million atoms

The results were as follows:

Timestep	Error in Displacement (absolute values)
0.01	0.000006
0.05	0.0005
0.075	0.002
0.1	0.005
0.15	0.08

Diffusion Coefficients (unit area/unit time)

State	From MSD (thousand atoms)	From MSD (million)	From VCAF (thousand atoms)	From VCAF (million)
Solid	$8.96\times 10^{-6}$	$4.5\times 10^{-6}$	$-1.66\times 10^{-4}$	$4.55\times 10^{-5}$
Liquid	$0.085$	$0.089$	$0.098$	$0.090$
Gas	$3.22$	$2.38$	$3.35$	$3.09$

We can see from the table that the diffusion coefficients calculated from the same number of atoms agree very well with each other, which suggests that the results are quite precise. However, we can assume that the million atoms simulations inevitably provide more accurate data.

The results also make sense logically, where the gas have the highest diffusion coefficient (they are able to diffuse much faster) than the liquid, and the diffusion coefficient of the solid is close to 0.

Mean Squared Displacement

From the graphs it can be concluded that the final 1000 timesteps follow linearity with which to apply the equations, and thus the diffusion coefficients can be calculated from the last timesteps.

JC: Show the lines of best fit on the graphs. Why does the gas MSD begin as a curve (ballistic motion).

Velocity Autocorrelation Function

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

Starting from a 1D Harmonic Oscillator:

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

Differentiating gives the velocity function:

$v\left(t\right)=-A\omega \sin (\omega t+\varphi )$

Substituting into the integrated form of the VAFC:

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

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

Then, using the trigonometric identity

$\sin (\alpha +\beta )=\sin \alpha \cos \beta +\sin \beta \cos \alpha ,$

$\sin (\omega t+\varphi +\omega \tau )=\sin (\omega t+\varphi )\cos \left(\omega \tau \right)+\sin \left(\omega \tau \right)\cos (\omega t+\varphi )$

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

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

$\tau$ is a constant, thus $\cos \left(\omega \tau \right)$ and $\sin \left(\omega \tau \right)$ can be taken out of the integral.

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

$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^2(\omega t+\varphi )$ can cancel out:

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

Since $\sin (x)$ is an odd function and $\cos (x)$ is an even function, $\sin (x)\cos (x)$ is odd.

$\therefore \cos (\omega t+\varphi )\sin (\omega t+\varphi )$ is an odd function (While the value of $\varphi$ can change whether or not the functions are even or odd, it shifts both functions by the same amount and the resulting $\cos (\omega t+\varphi )\sin (\omega t+\varphi )$ will still be odd)

$\therefore \underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\cos (\omega t+\varphi )\sin (\omega t+\varphi )=0$

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

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

JC: Good, clear derivation.

Comparison of Normalised VACF with VACF of Liquid and Solid

From the graph note that the liquid VACF tends quickly to zero, and the solid has small fluctuations. The difference arises as the in the liquid phase the particles have more translational freedom and can collide, exchanging energy between particles. Solid phase atoms retain a relatively static position around their equilibrium position. Also the solid rate will decay slower than in liquids.

For the harmonic oscillator, the VACF will not decay to zero as the there are atoms exchanging with their neighbours, which do not occur in harmonic oscillation.

File:Runningintegralaz.png

JC: The running integral plot has not uploaded properly. What does the minimum in the solid VACF represent physically?

The values for D from the VCAF model were greater in magnitude than those from the MSD calculations. The greatest source of error is most likely to be the use of the trapezium rule instead of an exact integral of each function.