ChemWiki - User contributions [en]

Rep:Mod:Az1114

2017-02-10T11:27:32Z

Az1114:

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

[[File:Displacementaz.png|thumb|600px|none]]

[[File:Energyintroaz.png|thumb|600px|none]]

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.

[[File:Erroraz.png|thumb|600px|none]]

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.

{| class="wikitable" border="1"
!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
|-

=== Lennard-Jones Potential ===
The below equation calculates the Leonard-Jones potential:

<math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>

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

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

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

<math>\frac{\mathrm{d}\phi}{\mathrm{d}r}=\mathbf{F}=-24\epsilon \left(\frac{2 \sigma^{12}}{r^{13}}-\frac{\sigma^6}{r^7}\right)</math>

<math>r_0=\sigma</math>

<math>\mathbf{F}=-24 \frac{\epsilon}{\sigma}</math>

The equilibrium separation (<math>r_{eq}</math>) 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 <math>\frac{\mathrm{d}\phi}{\mathrm{d}r}</math>, (or <math>\mathbf{F}</math>) to 0.

Hence,

<math>\mathbf{F}=-24\epsilon \left(\frac{2 \sigma^{12}}{r^{13}}-\frac{\sigma^6}{r^7}\right) 0</math>

<math>r_{eq}=\sqrt[6]{2\sigma^6}</math>

Substituting <math>r_{eq}</math> back into <math>\phi(r)</math> gives us the well depth:

<math>\phi(r_{eq})=-\epsilon</math>

----

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

<math>\int 4\epsilon \left(\frac{\sigma^{12}}{r^{12}}-\frac{\sigma^6}{r^6}\right)\mathrm{d}r</math>

<math>=4 \epsilon \left(\frac{\sigma^{12}}{11r^{11}}-\frac{\sigma^6}{5r^5}\right)\mathrm{d}r</math>

When <math>\sigma = \epsilon = 1.0:</math>

<math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=\lim_{r\rightarrow \infty} 4 \left(\frac{1^{12}}{11r^{11}}-\frac{1^6}{5r^5}\right) - 4\left(\frac{1^{12}}{11\times 2^{11}}+\frac{1^6}{5\times 2^5}\right)</math>

<math>=0-0.024822\cdots</math>

<math>\approx -0.0248</math>

Following the same calculations:

<math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=-0.0081767\cdots</math>

<math>\approx -0.00818</math>

<math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=-0.0032901\cdots</math>

<math>\approx -0.00329</math>

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

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

<math>1.00 \mathrm{ml} = 1.00 \mathrm{g}</math>

<math>18.02 \mathrm{g/mol} = 0.0555 \text{ mols of water}</math>

<math>0.0555 \times 6.022 \times 10^{23} = 3.34 \times 10^{22} \text{ molecules of water.}</math>

Conversely, 10,000 water molecules is only <math>2.99 \times 10^{-19} \text{ ml of water}</math>

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 <math>(0.5, 0.5, 0.5)</math> in a cube of dimensions <math>(0, 0, 0)</math> to <math>(1, 1, 1)</math> moves along the vector <math>(0.7, 0.6, 0.2)</math>, its final position would be <math>(0.2, 0.1, 0.7)</math>, 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 <math>\sigma=0.34 \mathrm{nm}</math>; <math>\epsilon/k_B=120K</math>.

If we set the cutoff to be <math>r*=3.2,</math>

<math>r=3.2\times0.34=1.088\mathrm{nm}</math>

Well depth: <math>\epsilon=120K\times k_B \times 6.02 \times 10^{23}</math>

<math>\epsilon=0.998 \mathrm{kJ/mol}</math>

and reduced temperature <math>T*=1.5</math> would be <math>T=1.5 \times 120=180\mathrm{K}.</math>

==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:

<math>0.8 \text{ points/volume}</math>

<math>\sqrt[3]{0.8}=0.9283\text{ points/length}</math>

<math>0.9283^{-1}=1.0772 \text{ length/point}</math>

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 <math>0.5928</math> as:

<math>\ \sqrt[3]{4.8^{-1}}=0.5928</math>

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

<pre>region box block 0 10 0 10 0 10</pre>

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:
<pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

The mass command takes two arguments, setting the mass of type 1 atoms as <math>1.0</math> (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 <math>1.0</math> setting the well depth <math>epsilon</math> and the latter sets <math>sigma</math>. 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 <math>\mathbf{x}_i(0)</math> & <math>\mathbf{v}_i\left(0\right)</math>, the Velocity-Verlet integration half-step algorithm can be used.

=== Running the Simulation ===
In the following lines from the input script:
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(100/${timestep})
timestep ${timestep}

### RUN SIMULATION ###
run ${n_steps}
</pre>

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
<pre>
timestep 0.001
run 100000
</pre>
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.

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

[[File:Total Energyaz.png|thumb|600px|none]] [[File:Temperatureaz.png|thumb|600px|none]] [[File:Pressureaz.png|thumb|600px|none]]

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.

[[File:Different Timesteps.png|thumb|600px|none]]

== 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 <math>T,</math> using the following equations:

<math>E_K = \frac{3}{2} N k_B T</math>

<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T \qquad\qquad\qquad(1)</math>

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

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B \mathfrak{T}\qquad\qquad(2)</math>

Thus we can solve for <math>\gamma.</math>:

From <math>(2)</math>:

<math>\gamma^2\sum_i m_i v_i^2 = 3N k_B \mathfrak{T}</math>

Then substituting <math>\sum_i m_i v_i^2</math> from <math>(1)</math>

<math>\gamma^2 3 N k_b T = 3N k_b \mathfrak{T}</math>

<math>\gamma = \sqrt{\frac{\epsilon}{T}}</math>

=== Explaining the Input Script ===

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

<pre>
fix aves all ave/time Nevery Nrepeat Nfreq
</pre>

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:

<pre>
fix aves all ave/time 100 1000 100000
run 100000
</pre>

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.

The timestep chosen for this simulation was 0.0025.

===Equations of State===

[[File:NpTaz.png|thumb|600px|none]]

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.

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

[[File:NVTaz.png|thumb|600px|none]]

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

[[File:RDFaz.jpg]]

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.

== Dynamical Properties and the Diffusion Coefficient ==

The Diffusion Coefficients <math>D</math> 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 <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

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:

{| class="wikitable" border="1"
|+ Diffusion Coefficients (unit area/unit time)
! State !! From MSD (thousand atoms) !! From MSD (million) !! From VCAF (thousand atoms) !! From VCAF (million)
|-
| Solid || <math>8.96 \times 10^{-6}</math> || <math>4.5 \times 10^{-6}</math> || <math>-1.66 \times 10^{-4}</math> || <math>4.55\times 10^{-5}</math>
|-
| Liquid || <math>0.085</math> || <math>0.089</math> ||<math>0.098</math> || <math>0.090</math>
|-
| Gas || <math>3.22</math> || <math>2.38</math> || <math>3.35</math> || <math>3.09</math>
|}

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 ===



[[File:MSD1az.jpg|thumb|600px|none]][[File:MSD2az.jpgthumb|600px|none]]

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.

=== Velocity Autocorrelation Function ===

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} v\left(t\right)v\left(t + \tau\right)\mathrm{d}t}{\int_{-\infty}^{\infty} v^2\left(t\right)\mathrm{d}t}</math>

Starting from a 1D Harmonic Oscillator:

<math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

Differentiating gives the velocity function:

<math> v\left(t\right) = -A\omega \sin\left(\omega t + \phi\right)</math>

Substituting into the integrated form of the VAFC:

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} A\omega \sin\left(\omega t + \phi \right) \cdot A \omega \sin \left(\omega \left(t+\tau\right)+\phi\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(A \omega \sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} \sin\left(\omega t + \phi \right) \cdot \sin \left(\omega \left(t+\tau\right)+\phi\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(\sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

Then, using the trigonometric identity

<math>\sin \left(\alpha + \beta \right)= \sin\alpha\cos\beta + \sin\beta\cos\alpha,</math>

<math>\sin\left(\omega t + \phi + \omega\tau\right) = \sin\left(\omega t+\phi\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\omega t + \phi\right)</math>

<math>\therefore C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} \sin\left(\omega t + \phi \right) \cdot \left(\sin\left(\omega t+\phi\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\omega t + \phi\right)\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(\sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

<math>C\left(\tau\right)=\frac{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)\cos\left(\omega \tau\right) + \int_{-\infty}^{\infty} \sin \left(\omega \tau\right) \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

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

<math>C\left(\tau\right)=\frac{\cos\left(\omega \tau\right)\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right) + \sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

<math>\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)</math> can cancel out:

<math>C\left(\tau\right)=\cos\left(\omega \tau\right)+\frac{\sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

Since <math>\sin(x)</math> is an odd function and <math>\cos(x)</math> is an even function, <math>\sin(x)\cos(x)</math> is odd.

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

<math>\therefore \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)=0</math>

<math>\therefore C\left(\tau\right)=\cos\left(\omega \tau\right)+\frac{\sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}=\cos\left(\omega \tau\right)</math>

<math>C\left(\tau\right)=\cos\left(\omega \tau\right)</math>

[[File:VACFaz.png|thumb|800px|Comparison of Normalised VACF with VACF of Liquid and Solid|left]]

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]]

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.

File:Different Timesteps.png

2017-02-10T11:27:03Z

Az1114: Az1114 uploaded a new version of File:Different Timesteps.png

Graph of Total energy versus time for different timesteps of simulation.

Rep:Mod:Az1114

2017-02-10T11:23:34Z

Az1114:

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

[[File:Displacementaz.png|thumb|600px|none]]

[[File:Energyintroaz.png|thumb|600px|none]]

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.

[[File:Erroraz.png|thumb|600px|none]]

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.

{| class="wikitable" border="1"
!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
|-

=== Lennard-Jones Potential ===
The below equation calculates the Leonard-Jones potential:

<math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>

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

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

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

<math>\frac{\mathrm{d}\phi}{\mathrm{d}r}=\mathbf{F}=-24\epsilon \left(\frac{2 \sigma^{12}}{r^{13}}-\frac{\sigma^6}{r^7}\right)</math>

<math>r_0=\sigma</math>

<math>\mathbf{F}=-24 \frac{\epsilon}{\sigma}</math>

The equilibrium separation (<math>r_{eq}</math>) 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 <math>\frac{\mathrm{d}\phi}{\mathrm{d}r}</math>, (or <math>\mathbf{F}</math>) to 0.

Hence,

<math>\mathbf{F}=-24\epsilon \left(\frac{2 \sigma^{12}}{r^{13}}-\frac{\sigma^6}{r^7}\right) 0</math>

<math>r_{eq}=\sqrt[6]{2\sigma^6}</math>

Substituting <math>r_{eq}</math> back into <math>\phi(r)</math> gives us the well depth:

<math>\phi(r_{eq})=-\epsilon</math>

----

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

<math>\int 4\epsilon \left(\frac{\sigma^{12}}{r^{12}}-\frac{\sigma^6}{r^6}\right)\mathrm{d}r</math>

<math>=4 \epsilon \left(\frac{\sigma^{12}}{11r^{11}}-\frac{\sigma^6}{5r^5}\right)\mathrm{d}r</math>

When <math>\sigma = \epsilon = 1.0:</math>

<math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=\lim_{r\rightarrow \infty} 4 \left(\frac{1^{12}}{11r^{11}}-\frac{1^6}{5r^5}\right) - 4\left(\frac{1^{12}}{11\times 2^{11}}+\frac{1^6}{5\times 2^5}\right)</math>

<math>=0-0.024822\cdots</math>

<math>\approx -0.0248</math>

Following the same calculations:

<math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=-0.0081767\cdots</math>

<math>\approx -0.00818</math>

<math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=-0.0032901\cdots</math>

<math>\approx -0.00329</math>

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

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

<math>1.00 \mathrm{ml} = 1.00 \mathrm{g}</math>

<math>18.02 \mathrm{g/mol} = 0.0555 \text{ mols of water}</math>

<math>0.0555 \times 6.022 \times 10^{23} = 3.34 \times 10^{22} \text{ molecules of water.}</math>

Conversely, 10,000 water molecules is only <math>2.99 \times 10^{-19} \text{ ml of water}</math>

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 <math>(0.5, 0.5, 0.5)</math> in a cube of dimensions <math>(0, 0, 0)</math> to <math>(1, 1, 1)</math> moves along the vector <math>(0.7, 0.6, 0.2)</math>, its final position would be <math>(0.2, 0.1, 0.7)</math>, 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 <math>\sigma=0.34 \mathrm{nm}</math>; <math>\epsilon/k_B=120K</math>.

If we set the cutoff to be <math>r*=3.2,</math>

<math>r=3.2\times0.34=1.088\mathrm{nm}</math>

Well depth: <math>\epsilon=120K\times k_B \times 6.02 \times 10^{23}</math>

<math>\epsilon=0.998 \mathrm{kJ/mol}</math>

and reduced temperature <math>T*=1.5</math> would be <math>T=1.5 \times 120=180\mathrm{K}.</math>

==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:

<math>0.8 \text{ points/volume}</math>

<math>\sqrt[3]{0.8}=0.9283\text{ points/length}</math>

<math>0.9283^{-1}=1.0772 \text{ length/point}</math>

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 <math>0.5928</math> as:

<math>\ \sqrt[3]{4.8^{-1}}=0.5928</math>

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

<pre>region box block 0 10 0 10 0 10</pre>

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:
<pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

The mass command takes two arguments, setting the mass of type 1 atoms as <math>1.0</math> (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 <math>1.0</math> setting the well depth <math>epsilon</math> and the latter sets <math>sigma</math>. 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 <math>\mathbf{x}_i(0)</math> & <math>\mathbf{v}_i\left(0\right)</math>, the Velocity-Verlet integration half-step algorithm can be used.

=== Running the Simulation ===
In the following lines from the input script:
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(100/${timestep})
timestep ${timestep}

### RUN SIMULATION ###
run ${n_steps}
</pre>

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
<pre>
timestep 0.001
run 100000
</pre>
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.

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

[[File:Total Energyaz.png|thumb|600px|none]]] [[File:Temperatureaz.png|thumb|600px|none]]] [[File:Pressureaz.png|thumb|600px|none]]]

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.

[[File:Different Timesteps.png|thumb|600px|none]]]

== 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 <math>T,</math> using the following equations:

<math>E_K = \frac{3}{2} N k_B T</math>

<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T \qquad\qquad\qquad(1)</math>

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

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B \mathfrak{T}\qquad\qquad(2)</math>

Thus we can solve for <math>\gamma.</math>:

From <math>(2)</math>:

<math>\gamma^2\sum_i m_i v_i^2 = 3N k_B \mathfrak{T}</math>

Then substituting <math>\sum_i m_i v_i^2</math> from <math>(1)</math>

<math>\gamma^2 3 N k_b T = 3N k_b \mathfrak{T}</math>

<math>\gamma = \sqrt{\frac{\epsilon}{T}}</math>

=== Explaining the Input Script ===

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

<pre>
fix aves all ave/time Nevery Nrepeat Nfreq
</pre>

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:

<pre>
fix aves all ave/time 100 1000 100000
run 100000
</pre>

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.

The timestep chosen for this simulation was 0.0025.

===Equations of State===

[[File:NpTaz.png|thumb|600px|none]]

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.

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

[[File:NVTaz.png|thumb|600px|none]]

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

[[File:RDFaz.jpg]]

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.

== Dynamical Properties and the Diffusion Coefficient ==

The Diffusion Coefficients <math>D</math> 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 <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

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:

{| class="wikitable" border="1"
|+ Diffusion Coefficients (unit area/unit time)
! State !! From MSD (thousand atoms) !! From MSD (million) !! From VCAF (thousand atoms) !! From VCAF (million)
|-
| Solid || <math>8.96 \times 10^{-6}</math> || <math>4.5 \times 10^{-6}</math> || <math>-1.66 \times 10^{-4}</math> || <math>4.55\times 10^{-5}</math>
|-
| Liquid || <math>0.085</math> || <math>0.089</math> ||<math>0.098</math> || <math>0.090</math>
|-
| Gas || <math>3.22</math> || <math>2.38</math> || <math>3.35</math> || <math>3.09</math>
|}

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 ===



[[File:MSD1az.jpg|thumb|600px|none]][[File:MSD2az.jpgthumb|600px|none]]

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.

=== Velocity Autocorrelation Function ===

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} v\left(t\right)v\left(t + \tau\right)\mathrm{d}t}{\int_{-\infty}^{\infty} v^2\left(t\right)\mathrm{d}t}</math>

Starting from a 1D Harmonic Oscillator:

<math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

Differentiating gives the velocity function:

<math> v\left(t\right) = -A\omega \sin\left(\omega t + \phi\right)</math>

Substituting into the integrated form of the VAFC:

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} A\omega \sin\left(\omega t + \phi \right) \cdot A \omega \sin \left(\omega \left(t+\tau\right)+\phi\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(A \omega \sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} \sin\left(\omega t + \phi \right) \cdot \sin \left(\omega \left(t+\tau\right)+\phi\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(\sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

Then, using the trigonometric identity

<math>\sin \left(\alpha + \beta \right)= \sin\alpha\cos\beta + \sin\beta\cos\alpha,</math>

<math>\sin\left(\omega t + \phi + \omega\tau\right) = \sin\left(\omega t+\phi\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\omega t + \phi\right)</math>

<math>\therefore C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} \sin\left(\omega t + \phi \right) \cdot \left(\sin\left(\omega t+\phi\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\omega t + \phi\right)\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(\sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

<math>C\left(\tau\right)=\frac{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)\cos\left(\omega \tau\right) + \int_{-\infty}^{\infty} \sin \left(\omega \tau\right) \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

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

<math>C\left(\tau\right)=\frac{\cos\left(\omega \tau\right)\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right) + \sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

<math>\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)</math> can cancel out:

<math>C\left(\tau\right)=\cos\left(\omega \tau\right)+\frac{\sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

Since <math>\sin(x)</math> is an odd function and <math>\cos(x)</math> is an even function, <math>\sin(x)\cos(x)</math> is odd.

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

<math>\therefore \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)=0</math>

<math>\therefore C\left(\tau\right)=\cos\left(\omega \tau\right)+\frac{\sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}=\cos\left(\omega \tau\right)</math>

<math>C\left(\tau\right)=\cos\left(\omega \tau\right)</math>

[[File:VACFaz.png|thumb|800px|Comparison of Normalised VACF with VACF of Liquid and Solid|left]]

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]]

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.

File:Pressureaz.png

2017-02-10T11:22:35Z

Az1114: Az1114 uploaded a new version of File:Pressureaz.png

File:Temperatureaz.png

2017-02-10T11:22:07Z

Az1114: Az1114 uploaded a new version of File:Temperatureaz.png

File:Total Energyaz.png

2017-02-10T11:21:20Z

Az1114: Az1114 uploaded a new version of File:Total Energyaz.png

File:MSD1az.jpg

2017-02-10T11:14:57Z

Az1114: Az1114 uploaded a new version of File:MSD1az.jpg

Rep:Mod:Az1114

2017-02-10T11:07:19Z

Az1114: /* Radial Distribution Functions */

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

[[File:Displacementaz.png|thumb|600px|none]]

[[File:Energyintroaz.png|thumb|600px|none]]

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.

[[File:Erroraz.png|thumb|600px|none]]

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.

{| class="wikitable" border="1"
!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
|-

=== Lennard-Jones Potential ===
The below equation calculates the Leonard-Jones potential:

<math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>

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

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

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

<math>\frac{\mathrm{d}\phi}{\mathrm{d}r}=\mathbf{F}=-24\epsilon \left(\frac{2 \sigma^{12}}{r^{13}}-\frac{\sigma^6}{r^7}\right)</math>

<math>r_0=\sigma</math>

<math>\mathbf{F}=-24 \frac{\epsilon}{\sigma}</math>

The equilibrium separation (<math>r_{eq}</math>) 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 <math>\frac{\mathrm{d}\phi}{\mathrm{d}r}</math>, (or <math>\mathbf{F}</math>) to 0.

Hence,

<math>\mathbf{F}=-24\epsilon \left(\frac{2 \sigma^{12}}{r^{13}}-\frac{\sigma^6}{r^7}\right) 0</math>

<math>r_{eq}=\sqrt[6]{2\sigma^6}</math>

Substituting <math>r_{eq}</math> back into <math>\phi(r)</math> gives us the well depth:

<math>\phi(r_{eq})=-\epsilon</math>

----

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

<math>\int 4\epsilon \left(\frac{\sigma^{12}}{r^{12}}-\frac{\sigma^6}{r^6}\right)\mathrm{d}r</math>

<math>=4 \epsilon \left(\frac{\sigma^{12}}{11r^{11}}-\frac{\sigma^6}{5r^5}\right)\mathrm{d}r</math>

When <math>\sigma = \epsilon = 1.0:</math>

<math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=\lim_{r\rightarrow \infty} 4 \left(\frac{1^{12}}{11r^{11}}-\frac{1^6}{5r^5}\right) - 4\left(\frac{1^{12}}{11\times 2^{11}}+\frac{1^6}{5\times 2^5}\right)</math>

<math>=0-0.024822\cdots</math>

<math>\approx -0.0248</math>

Following the same calculations:

<math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=-0.0081767\cdots</math>

<math>\approx -0.00818</math>

<math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=-0.0032901\cdots</math>

<math>\approx -0.00329</math>

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

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

<math>1.00 \mathrm{ml} = 1.00 \mathrm{g}</math>

<math>18.02 \mathrm{g/mol} = 0.0555 \text{ mols of water}</math>

<math>0.0555 \times 6.022 \times 10^{23} = 3.34 \times 10^{22} \text{ molecules of water.}</math>

Conversely, 10,000 water molecules is only <math>2.99 \times 10^{-19} \text{ ml of water}</math>

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 <math>(0.5, 0.5, 0.5)</math> in a cube of dimensions <math>(0, 0, 0)</math> to <math>(1, 1, 1)</math> moves along the vector <math>(0.7, 0.6, 0.2)</math>, its final position would be <math>(0.2, 0.1, 0.7)</math>, 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 <math>\sigma=0.34 \mathrm{nm}</math>; <math>\epsilon/k_B=120K</math>.

If we set the cutoff to be <math>r*=3.2,</math>

<math>r=3.2\times0.34=1.088\mathrm{nm}</math>

Well depth: <math>\epsilon=120K\times k_B \times 6.02 \times 10^{23}</math>

<math>\epsilon=0.998 \mathrm{kJ/mol}</math>

and reduced temperature <math>T*=1.5</math> would be <math>T=1.5 \times 120=180\mathrm{K}.</math>

==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:

<math>0.8 \text{ points/volume}</math>

<math>\sqrt[3]{0.8}=0.9283\text{ points/length}</math>

<math>0.9283^{-1}=1.0772 \text{ length/point}</math>

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 <math>0.5928</math> as:

<math>\ \sqrt[3]{4.8^{-1}}=0.5928</math>

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

<pre>region box block 0 10 0 10 0 10</pre>

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:
<pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

The mass command takes two arguments, setting the mass of type 1 atoms as <math>1.0</math> (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 <math>1.0</math> setting the well depth <math>epsilon</math> and the latter sets <math>sigma</math>. 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 <math>\mathbf{x}_i(0)</math> & <math>\mathbf{v}_i\left(0\right)</math>, the Velocity-Verlet integration half-step algorithm can be used.

=== Running the Simulation ===
In the following lines from the input script:
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(100/${timestep})
timestep ${timestep}

### RUN SIMULATION ###
run ${n_steps}
</pre>

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
<pre>
timestep 0.001
run 100000
</pre>
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.

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

[[File:Total Energyaz.png]] [[File:Temperatureaz.png]] [[File:Pressureaz.png]]

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.

[[File:Different Timesteps.png]]

== 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 <math>T,</math> using the following equations:

<math>E_K = \frac{3}{2} N k_B T</math>

<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T \qquad\qquad\qquad(1)</math>

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

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B \mathfrak{T}\qquad\qquad(2)</math>

Thus we can solve for <math>\gamma.</math>:

From <math>(2)</math>:

<math>\gamma^2\sum_i m_i v_i^2 = 3N k_B \mathfrak{T}</math>

Then substituting <math>\sum_i m_i v_i^2</math> from <math>(1)</math>

<math>\gamma^2 3 N k_b T = 3N k_b \mathfrak{T}</math>

<math>\gamma = \sqrt{\frac{\epsilon}{T}}</math>

=== Explaining the Input Script ===

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

<pre>
fix aves all ave/time Nevery Nrepeat Nfreq
</pre>

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:

<pre>
fix aves all ave/time 100 1000 100000
run 100000
</pre>

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.

The timestep chosen for this simulation was 0.0025.

===Equations of State===

[[File:NpTaz.png|thumb|600px|none]]

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.

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

[[File:NVTaz.png|thumb|600px|none]]

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

[[File:RDFaz.jpg]]

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.

== Dynamical Properties and the Diffusion Coefficient ==

The Diffusion Coefficients <math>D</math> 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 <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

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:

{| class="wikitable" border="1"
|+ Diffusion Coefficients (unit area/unit time)
! State !! From MSD (thousand atoms) !! From MSD (million) !! From VCAF (thousand atoms) !! From VCAF (million)
|-
| Solid || <math>8.96 \times 10^{-6}</math> || <math>4.5 \times 10^{-6}</math> || <math>-1.66 \times 10^{-4}</math> || <math>4.55\times 10^{-5}</math>
|-
| Liquid || <math>0.085</math> || <math>0.089</math> ||<math>0.098</math> || <math>0.090</math>
|-
| Gas || <math>3.22</math> || <math>2.38</math> || <math>3.35</math> || <math>3.09</math>
|}

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 ===



[[File:MSD1az.jpg|thumb|600px|none]][[File:MSD2az.jpgthumb|600px|none]]

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.

=== Velocity Autocorrelation Function ===

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} v\left(t\right)v\left(t + \tau\right)\mathrm{d}t}{\int_{-\infty}^{\infty} v^2\left(t\right)\mathrm{d}t}</math>

Starting from a 1D Harmonic Oscillator:

<math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

Differentiating gives the velocity function:

<math> v\left(t\right) = -A\omega \sin\left(\omega t + \phi\right)</math>

Substituting into the integrated form of the VAFC:

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} A\omega \sin\left(\omega t + \phi \right) \cdot A \omega \sin \left(\omega \left(t+\tau\right)+\phi\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(A \omega \sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} \sin\left(\omega t + \phi \right) \cdot \sin \left(\omega \left(t+\tau\right)+\phi\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(\sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

Then, using the trigonometric identity

<math>\sin \left(\alpha + \beta \right)= \sin\alpha\cos\beta + \sin\beta\cos\alpha,</math>

<math>\sin\left(\omega t + \phi + \omega\tau\right) = \sin\left(\omega t+\phi\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\omega t + \phi\right)</math>

<math>\therefore C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} \sin\left(\omega t + \phi \right) \cdot \left(\sin\left(\omega t+\phi\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\omega t + \phi\right)\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(\sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

<math>C\left(\tau\right)=\frac{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)\cos\left(\omega \tau\right) + \int_{-\infty}^{\infty} \sin \left(\omega \tau\right) \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

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

<math>C\left(\tau\right)=\frac{\cos\left(\omega \tau\right)\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right) + \sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

<math>\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)</math> can cancel out:

<math>C\left(\tau\right)=\cos\left(\omega \tau\right)+\frac{\sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

Since <math>\sin(x)</math> is an odd function and <math>\cos(x)</math> is an even function, <math>\sin(x)\cos(x)</math> is odd.

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

<math>\therefore \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)=0</math>

<math>\therefore C\left(\tau\right)=\cos\left(\omega \tau\right)+\frac{\sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}=\cos\left(\omega \tau\right)</math>

<math>C\left(\tau\right)=\cos\left(\omega \tau\right)</math>

[[File:VACFaz.png|thumb|800px|Comparison of Normalised VACF with VACF of Liquid and Solid|left]]

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]]

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.

Rep:Mod:Az1114

2017-02-10T10:58:49Z

Az1114: /* Radial Distribution Functions */

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

[[File:Displacementaz.png|thumb|600px|none]]

[[File:Energyintroaz.png|thumb|600px|none]]

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.

[[File:Erroraz.png|thumb|600px|none]]

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.

{| class="wikitable" border="1"
!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
|-

=== Lennard-Jones Potential ===
The below equation calculates the Leonard-Jones potential:

<math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>

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

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

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

<math>\frac{\mathrm{d}\phi}{\mathrm{d}r}=\mathbf{F}=-24\epsilon \left(\frac{2 \sigma^{12}}{r^{13}}-\frac{\sigma^6}{r^7}\right)</math>

<math>r_0=\sigma</math>

<math>\mathbf{F}=-24 \frac{\epsilon}{\sigma}</math>

The equilibrium separation (<math>r_{eq}</math>) 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 <math>\frac{\mathrm{d}\phi}{\mathrm{d}r}</math>, (or <math>\mathbf{F}</math>) to 0.

Hence,

<math>\mathbf{F}=-24\epsilon \left(\frac{2 \sigma^{12}}{r^{13}}-\frac{\sigma^6}{r^7}\right) 0</math>

<math>r_{eq}=\sqrt[6]{2\sigma^6}</math>

Substituting <math>r_{eq}</math> back into <math>\phi(r)</math> gives us the well depth:

<math>\phi(r_{eq})=-\epsilon</math>

----

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

<math>\int 4\epsilon \left(\frac{\sigma^{12}}{r^{12}}-\frac{\sigma^6}{r^6}\right)\mathrm{d}r</math>

<math>=4 \epsilon \left(\frac{\sigma^{12}}{11r^{11}}-\frac{\sigma^6}{5r^5}\right)\mathrm{d}r</math>

When <math>\sigma = \epsilon = 1.0:</math>

<math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=\lim_{r\rightarrow \infty} 4 \left(\frac{1^{12}}{11r^{11}}-\frac{1^6}{5r^5}\right) - 4\left(\frac{1^{12}}{11\times 2^{11}}+\frac{1^6}{5\times 2^5}\right)</math>

<math>=0-0.024822\cdots</math>

<math>\approx -0.0248</math>

Following the same calculations:

<math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=-0.0081767\cdots</math>

<math>\approx -0.00818</math>

<math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=-0.0032901\cdots</math>

<math>\approx -0.00329</math>

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

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

<math>1.00 \mathrm{ml} = 1.00 \mathrm{g}</math>

<math>18.02 \mathrm{g/mol} = 0.0555 \text{ mols of water}</math>

<math>0.0555 \times 6.022 \times 10^{23} = 3.34 \times 10^{22} \text{ molecules of water.}</math>

Conversely, 10,000 water molecules is only <math>2.99 \times 10^{-19} \text{ ml of water}</math>

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 <math>(0.5, 0.5, 0.5)</math> in a cube of dimensions <math>(0, 0, 0)</math> to <math>(1, 1, 1)</math> moves along the vector <math>(0.7, 0.6, 0.2)</math>, its final position would be <math>(0.2, 0.1, 0.7)</math>, 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 <math>\sigma=0.34 \mathrm{nm}</math>; <math>\epsilon/k_B=120K</math>.

If we set the cutoff to be <math>r*=3.2,</math>

<math>r=3.2\times0.34=1.088\mathrm{nm}</math>

Well depth: <math>\epsilon=120K\times k_B \times 6.02 \times 10^{23}</math>

<math>\epsilon=0.998 \mathrm{kJ/mol}</math>

and reduced temperature <math>T*=1.5</math> would be <math>T=1.5 \times 120=180\mathrm{K}.</math>

==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:

<math>0.8 \text{ points/volume}</math>

<math>\sqrt[3]{0.8}=0.9283\text{ points/length}</math>

<math>0.9283^{-1}=1.0772 \text{ length/point}</math>

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 <math>0.5928</math> as:

<math>\ \sqrt[3]{4.8^{-1}}=0.5928</math>

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

<pre>region box block 0 10 0 10 0 10</pre>

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:
<pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

The mass command takes two arguments, setting the mass of type 1 atoms as <math>1.0</math> (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 <math>1.0</math> setting the well depth <math>epsilon</math> and the latter sets <math>sigma</math>. 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 <math>\mathbf{x}_i(0)</math> & <math>\mathbf{v}_i\left(0\right)</math>, the Velocity-Verlet integration half-step algorithm can be used.

=== Running the Simulation ===
In the following lines from the input script:
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(100/${timestep})
timestep ${timestep}

### RUN SIMULATION ###
run ${n_steps}
</pre>

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
<pre>
timestep 0.001
run 100000
</pre>
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.

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

[[File:Total Energyaz.png]] [[File:Temperatureaz.png]] [[File:Pressureaz.png]]

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.

[[File:Different Timesteps.png]]

== 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 <math>T,</math> using the following equations:

<math>E_K = \frac{3}{2} N k_B T</math>

<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T \qquad\qquad\qquad(1)</math>

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

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B \mathfrak{T}\qquad\qquad(2)</math>

Thus we can solve for <math>\gamma.</math>:

From <math>(2)</math>:

<math>\gamma^2\sum_i m_i v_i^2 = 3N k_B \mathfrak{T}</math>

Then substituting <math>\sum_i m_i v_i^2</math> from <math>(1)</math>

<math>\gamma^2 3 N k_b T = 3N k_b \mathfrak{T}</math>

<math>\gamma = \sqrt{\frac{\epsilon}{T}}</math>

=== Explaining the Input Script ===

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

<pre>
fix aves all ave/time Nevery Nrepeat Nfreq
</pre>

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:

<pre>
fix aves all ave/time 100 1000 100000
run 100000
</pre>

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.

The timestep chosen for this simulation was 0.0025.

===Equations of State===

[[File:NpTaz.png|thumb|600px|none]]

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.

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

[[File:NVTaz.png|thumb|600px|none]]

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

[[File:RDFaz.jpg]]

== Dynamical Properties and the Diffusion Coefficient ==

The Diffusion Coefficients <math>D</math> 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 <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

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:

{| class="wikitable" border="1"
|+ Diffusion Coefficients (unit area/unit time)
! State !! From MSD (thousand atoms) !! From MSD (million) !! From VCAF (thousand atoms) !! From VCAF (million)
|-
| Solid || <math>8.96 \times 10^{-6}</math> || <math>4.5 \times 10^{-6}</math> || <math>-1.66 \times 10^{-4}</math> || <math>4.55\times 10^{-5}</math>
|-
| Liquid || <math>0.085</math> || <math>0.089</math> ||<math>0.098</math> || <math>0.090</math>
|-
| Gas || <math>3.22</math> || <math>2.38</math> || <math>3.35</math> || <math>3.09</math>
|}

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 ===



[[File:MSD1az.jpg|thumb|600px|none]][[File:MSD2az.jpgthumb|600px|none]]

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.

=== Velocity Autocorrelation Function ===

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} v\left(t\right)v\left(t + \tau\right)\mathrm{d}t}{\int_{-\infty}^{\infty} v^2\left(t\right)\mathrm{d}t}</math>

Starting from a 1D Harmonic Oscillator:

<math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

Differentiating gives the velocity function:

<math> v\left(t\right) = -A\omega \sin\left(\omega t + \phi\right)</math>

Substituting into the integrated form of the VAFC:

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} A\omega \sin\left(\omega t + \phi \right) \cdot A \omega \sin \left(\omega \left(t+\tau\right)+\phi\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(A \omega \sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} \sin\left(\omega t + \phi \right) \cdot \sin \left(\omega \left(t+\tau\right)+\phi\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(\sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

Then, using the trigonometric identity

<math>\sin \left(\alpha + \beta \right)= \sin\alpha\cos\beta + \sin\beta\cos\alpha,</math>

<math>\sin\left(\omega t + \phi + \omega\tau\right) = \sin\left(\omega t+\phi\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\omega t + \phi\right)</math>

<math>\therefore C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} \sin\left(\omega t + \phi \right) \cdot \left(\sin\left(\omega t+\phi\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\omega t + \phi\right)\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(\sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

<math>C\left(\tau\right)=\frac{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)\cos\left(\omega \tau\right) + \int_{-\infty}^{\infty} \sin \left(\omega \tau\right) \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

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

<math>C\left(\tau\right)=\frac{\cos\left(\omega \tau\right)\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right) + \sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

<math>\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)</math> can cancel out:

<math>C\left(\tau\right)=\cos\left(\omega \tau\right)+\frac{\sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

Since <math>\sin(x)</math> is an odd function and <math>\cos(x)</math> is an even function, <math>\sin(x)\cos(x)</math> is odd.

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

<math>\therefore \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)=0</math>

<math>\therefore C\left(\tau\right)=\cos\left(\omega \tau\right)+\frac{\sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}=\cos\left(\omega \tau\right)</math>

<math>C\left(\tau\right)=\cos\left(\omega \tau\right)</math>

[[File:VACFaz.png|thumb|800px|Comparison of Normalised VACF with VACF of Liquid and Solid|left]]

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]]

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.

Rep:Mod:Az1114

2017-02-10T10:53:13Z

2017-02-10T10:16:05Z

Az1114: /* Periodic Boundary Conditions */

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

[[File:Displacementaz.png|thumb|600px|none]]

[[File:Energyintroaz.png|thumb|600px|none]]

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.

[[File:Erroraz.png|thumb|600px|none]]

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.

{| class="wikitable" border="1"
!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
|-

=== Lennard-Jones Potential ===
The below equation calculates the Leonard-Jones potential:

<math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>

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

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

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

<math>\frac{\mathrm{d}\phi}{\mathrm{d}r}=\mathbf{F}=-24\epsilon \left(\frac{2 \sigma^{12}}{r^{13}}-\frac{\sigma^6}{r^7}\right)</math>

<math>r_0=\sigma</math>

<math>\mathbf{F}=-24 \frac{\epsilon}{\sigma}</math>

The equilibrium separation (<math>r_{eq}</math>) 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 <math>\frac{\mathrm{d}\phi}{\mathrm{d}r}</math>, (or <math>\mathbf{F}</math>) to 0.

Hence,

<math>\mathbf{F}=-24\epsilon \left(\frac{2 \sigma^{12}}{r^{13}}-\frac{\sigma^6}{r^7}\right) 0</math>

<math>r_{eq}=\sqrt[6]{2\sigma^6}</math>

Substituting <math>r_{eq}</math> back into <math>\phi(r)</math> gives us the well depth:

<math>\phi(r_{eq})=-\epsilon</math>

----

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

<math>\int 4\epsilon \left(\frac{\sigma^{12}}{r^{12}}-\frac{\sigma^6}{r^6}\right)\mathrm{d}r</math>

<math>=4 \epsilon \left(\frac{\sigma^{12}}{11r^{11}}-\frac{\sigma^6}{5r^5}\right)\mathrm{d}r</math>

When <math>\sigma = \epsilon = 1.0:</math>

<math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=\lim_{r\rightarrow \infty} 4 \left(\frac{1^{12}}{11r^{11}}-\frac{1^6}{5r^5}\right) - 4\left(\frac{1^{12}}{11\times 2^{11}}+\frac{1^6}{5\times 2^5}\right)</math>

<math>=0-0.024822\cdots</math>

<math>\approx -0.0248</math>

Following the same calculations:

<math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=-0.0081767\cdots</math>

<math>\approx -0.00818</math>

<math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=-0.0032901\cdots</math>

<math>\approx -0.00329</math>

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

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

<math>1.00 \mathrm{ml} = 1.00 \mathrm{g}</math>

<math>18.02 \mathrm{g/mol} = 0.0555 \text{ mols of water}</math>

<math>0.0555 \times 6.022 \times 10^{23} = 3.34 \times 10^{22} \text{ molecules of water.}</math>

Conversely, 10,000 water molecules is only <math>2.99 \times 10^{-19} \text{ ml of water}</math>

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 <math>(0.5, 0.5, 0.5)</math> in a cube of dimensions <math>(0, 0, 0)</math> to <math>(1, 1, 1)</math> moves along the vector <math>(0.7, 0.6, 0.2)</math>, its final position would be <math>(0.2, 0.1, 0.7)</math>, 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 <math>\sigma=0.34 \mathrm{nm}</math>; <math>\epsilon/k_B=120K</math>.

If we set the cutoff to be <math>r*=3.2,</math>

<math>r=3.2\times0.34=1.088\mathrm{nm}</math>

Well depth: <math>\epsilon=120K\times k_B \times 6.02 \times 10^{23}</math>

<math>\epsilon=0.998 \mathrm{kJ/mol}</math>

and reduced temperature <math>T*=1.5</math> would be <math>T=1.5 \times 120=180\mathrm{K}.</math>

==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:

<math>0.8 \text{ points/volume}</math>

<math>\sqrt[3]{0.8}=0.9283\text{ points/length}</math>

<math>0.9283^{-1}=1.0772 \text{ length/point}</math>

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 <math>0.5928</math> as:

<math>\ \sqrt[3]{4.8^{-1}}=0.5928</math>

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

<pre>region box block 0 10 0 10 0 10</pre>

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:
<pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

The mass command takes two arguments, setting the mass of type 1 atoms as <math>1.0</math> (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 <math>1.0</math> setting the well depth <math>epsilon</math> and the latter sets <math>sigma</math>. 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 <math>\mathbf{x}_i(0)</math> & <math>\mathbf{v}_i\left(0\right)</math>, the Velocity-Verlet integration half-step algorithm can be used.

=== Running the Simulation ===
In the following lines from the input script:
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(100/${timestep})
timestep ${timestep}

### RUN SIMULATION ###
run ${n_steps}
</pre>

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
<pre>
timestep 0.001
run 100000
</pre>
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.

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

[[File:Total Energyaz.png]] [[File:Temperatureaz.png]] [[File:Pressureaz.png]]

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.

[[File:Different Timesteps.png]]

== 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 <math>T,</math> using the following equations:

<math>E_K = \frac{3}{2} N k_B T</math>

<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T \qquad\qquad\qquad(1)</math>

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

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B \mathfrak{T}\qquad\qquad(2)</math>

Thus we can solve for <math>\gamma.</math>:

From <math>(2)</math>:

<math>\gamma^2\sum_i m_i v_i^2 = 3N k_B \mathfrak{T}</math>

Then substituting <math>\sum_i m_i v_i^2</math> from <math>(1)</math>

<math>\gamma^2 3 N k_b T = 3N k_b \mathfrak{T}</math>

<math>\gamma = \sqrt{\frac{\epsilon}{T}}</math>

=== Explaining the Input Script ===

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

<pre>
fix aves all ave/time Nevery Nrepeat Nfreq
</pre>

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:

<pre>
fix aves all ave/time 100 1000 100000
run 100000
</pre>

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.

The timestep chosen for this simulation was 0.0025. Thus, the simulation will run for <math>100000 \times 0.0025 = 250 \text{ unit time}</math>

===Plotting the equations of state===
[[File:Densityaz]]

The above graph displays the density (both ideal and from the simulation) versus temperature.

Note that there is a significant discrepancy between simulation and theory, the theoretical being greater in magnitude. The error bars cannot be seen in the graph as the difference between the ideal and simulated densities is so high that the error bars become too small to visualise.

Factors leading to this discrepancy may involve the model of the Leonard-Jones interaction in the simulation. The ideal gas law relies on there being no inter-particle interaction, and thus if there is a repulsive force due to the Leonard-Jones interactions, then the density may be lower than the ideal calculation.

It can also be seen that as the temperature increases, the discrepancy between the ideal and simulated densities becomes smaller.

==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 short-range order. They also exhibit rotational and to some extent translational motion.

[[File:RDFaz.jpg]]

== Dynamical Properties and the Diffusion Coefficient ==

The Diffusion Coefficients <math>D</math> 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 <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

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:

{| class="wikitable" border="1"
|+ Diffusion Coefficients (unit area/unit time)
! State !! From MSD (thousand atoms) !! From MSD (million) !! From VCAF (thousand atoms) !! From VCAF (million)
|-
| Solid || <math>8.96 \times 10^{-6}</math> || <math>4.5 \times 10^{-6}</math> || <math>-1.66 \times 10^{-4}</math> || <math>4.55\times 10^{-5}</math>
|-
| Liquid || <math>0.085</math> || <math>0.089</math> ||<math>0.098</math> || <math>0.090</math>
|-
| Gas || <math>3.22</math> || <math>2.38</math> || <math>3.35</math> || <math>3.09</math>
|}

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 ===



[[File:MSD1az.jpg|thumb|600px|none]][[File:MSD2az.jpgthumb|600px|none]]

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.

=== Velocity Autocorrelation Function ===

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} v\left(t\right)v\left(t + \tau\right)\mathrm{d}t}{\int_{-\infty}^{\infty} v^2\left(t\right)\mathrm{d}t}</math>

Starting from a 1D Harmonic Oscillator:

<math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

Differentiating gives the velocity function:

<math> v\left(t\right) = -A\omega \sin\left(\omega t + \phi\right)</math>

Substituting into the integrated form of the VAFC:

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} A\omega \sin\left(\omega t + \phi \right) \cdot A \omega \sin \left(\omega \left(t+\tau\right)+\phi\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(A \omega \sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} \sin\left(\omega t + \phi \right) \cdot \sin \left(\omega \left(t+\tau\right)+\phi\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(\sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

Then, using the trigonometric identity

<math>\sin \left(\alpha + \beta \right)= \sin\alpha\cos\beta + \sin\beta\cos\alpha,</math>

<math>\sin\left(\omega t + \phi + \omega\tau\right) = \sin\left(\omega t+\phi\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\omega t + \phi\right)</math>

<math>\therefore C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} \sin\left(\omega t + \phi \right) \cdot \left(\sin\left(\omega t+\phi\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\omega t + \phi\right)\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(\sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

<math>C\left(\tau\right)=\frac{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)\cos\left(\omega \tau\right) + \int_{-\infty}^{\infty} \sin \left(\omega \tau\right) \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

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

<math>C\left(\tau\right)=\frac{\cos\left(\omega \tau\right)\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right) + \sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

<math>\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)</math> can cancel out:

<math>C\left(\tau\right)=\cos\left(\omega \tau\right)+\frac{\sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

Since <math>\sin(x)</math> is an odd function and <math>\cos(x)</math> is an even function, <math>\sin(x)\cos(x)</math> is odd.

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

<math>\therefore \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)=0</math>

<math>\therefore C\left(\tau\right)=\cos\left(\omega \tau\right)+\frac{\sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}=\cos\left(\omega \tau\right)</math>

<math>C\left(\tau\right)=\cos\left(\omega \tau\right)</math>

[[File:VACFaz.png|thumb|800px|Comparison of Normalised VACF with VACF of Liquid and Solid|left]]

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.

Rep:Mod:Az1114

2017-02-10T10:14:33Z

Az1114: /* Dynamical Properties and the Diffusion Coefficient */

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

[[File:Displacementaz.png|thumb|600px|none]]

[[File:Energyintroaz.png|thumb|600px|none]]

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.

[[File:Erroraz.png|thumb|600px|none]]

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.

{| class="wikitable" border="1"
!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
|-

=== Lennard-Jones Potential ===
The below equation calculates the Leonard-Jones potential:

<math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>

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

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

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

<math>\frac{\mathrm{d}\phi}{\mathrm{d}r}=\mathbf{F}=-24\epsilon \left(\frac{2 \sigma^{12}}{r^{13}}-\frac{\sigma^6}{r^7}\right)</math>

<math>r_0=\sigma</math>

<math>\mathbf{F}=-24 \frac{\epsilon}{\sigma}</math>

The equilibrium separation (<math>r_{eq}</math>) 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 <math>\frac{\mathrm{d}\phi}{\mathrm{d}r}</math>, (or <math>\mathbf{F}</math>) to 0.

Hence,

<math>\mathbf{F}=-24\epsilon \left(\frac{2 \sigma^{12}}{r^{13}}-\frac{\sigma^6}{r^7}\right) 0</math>

<math>r_{eq}=\sqrt[6]{2\sigma^6}</math>

Substituting <math>r_{eq}</math> back into <math>\phi(r)</math> gives us the well depth:

<math>\phi(r_{eq})=-\epsilon</math>

----

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

<math>\int 4\epsilon \left(\frac{\sigma^{12}}{r^{12}}-\frac{\sigma^6}{r^6}\right)\mathrm{d}r</math>

<math>=4 \epsilon \left(\frac{\sigma^{12}}{11r^{11}}-\frac{\sigma^6}{5r^5}\right)\mathrm{d}r</math>

When <math>\sigma = \epsilon = 1.0:</math>

<math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=\lim_{r\rightarrow \infty} 4 \left(\frac{1^{12}}{11r^{11}}-\frac{1^6}{5r^5}\right) - 4\left(\frac{1^{12}}{11\times 2^{11}}+\frac{1^6}{5\times 2^5}\right)</math>

<math>=0-0.024822\cdots</math>

<math>\approx -0.0248</math>

Following the same calculations:

<math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=-0.0081767\cdots</math>

<math>\approx -0.00818</math>

<math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=-0.0032901\cdots</math>

<math>\approx -0.00329</math>

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

=== 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, let us first calculate the number of molecules in 1ml of water:

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

<math>1.00 \mathrm{ml} = 1.00 \mathrm{g}</math>

<math>18.02 \mathrm{g/mol} = 0.0555 \text{ mols of water}</math>

<math>0.0555 \times 6.022 \times 10^{23} = 3.34 \times 10^{22} \text{ molecules of water.}</math>

On the other hand, 10,000 water molecules is only <math>2.99 \times 10^{-19} \text{ ml of water}</math>

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 direction, 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 <math>(0.5, 0.5, 0.5)</math> in a cube of dimensions <math>(0, 0, 0)</math> to <math>(1, 1, 1)</math> moves along the vector <math>(0.7, 0.6, 0.2)</math>, its final position would be <math>(0.2, 0.1, 0.7)</math>, 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 <math>\sigma=0.34 \mathrm{nm}</math>; <math>\epsilon/k_B=120K</math>.

If we set the cutoff to be <math>r*=3.2,</math>

<math>r=3.2\times0.34=1.088\mathrm{nm}</math>

Well depth: <math>\epsilon=120K\times k_B \times 6.02 \times 10^{23}</math>

<math>\epsilon=0.998 \mathrm{kJ/mol}</math>

and reduced temperature <math>T*=1.5</math> would be <math>T=1.5 \times 120=180\mathrm{K}.</math>

==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:

<math>0.8 \text{ points/volume}</math>

<math>\sqrt[3]{0.8}=0.9283\text{ points/length}</math>

<math>0.9283^{-1}=1.0772 \text{ length/point}</math>

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 <math>0.5928</math> as:

<math>\ \sqrt[3]{4.8^{-1}}=0.5928</math>

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

<pre>region box block 0 10 0 10 0 10</pre>

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:
<pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

The mass command takes two arguments, setting the mass of type 1 atoms as <math>1.0</math> (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 <math>1.0</math> setting the well depth <math>epsilon</math> and the latter sets <math>sigma</math>. 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 <math>\mathbf{x}_i(0)</math> & <math>\mathbf{v}_i\left(0\right)</math>, the Velocity-Verlet integration half-step algorithm can be used.

=== Running the Simulation ===
In the following lines from the input script:
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(100/${timestep})
timestep ${timestep}

### RUN SIMULATION ###
run ${n_steps}
</pre>

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
<pre>
timestep 0.001
run 100000
</pre>
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.

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

[[File:Total Energyaz.png]] [[File:Temperatureaz.png]] [[File:Pressureaz.png]]

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.

[[File:Different Timesteps.png]]

== 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 <math>T,</math> using the following equations:

<math>E_K = \frac{3}{2} N k_B T</math>

<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T \qquad\qquad\qquad(1)</math>

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

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B \mathfrak{T}\qquad\qquad(2)</math>

Thus we can solve for <math>\gamma.</math>:

From <math>(2)</math>:

<math>\gamma^2\sum_i m_i v_i^2 = 3N k_B \mathfrak{T}</math>

Then substituting <math>\sum_i m_i v_i^2</math> from <math>(1)</math>

<math>\gamma^2 3 N k_b T = 3N k_b \mathfrak{T}</math>

<math>\gamma = \sqrt{\frac{\epsilon}{T}}</math>

=== Explaining the Input Script ===

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

<pre>
fix aves all ave/time Nevery Nrepeat Nfreq
</pre>

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:

<pre>
fix aves all ave/time 100 1000 100000
run 100000
</pre>

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.

The timestep chosen for this simulation was 0.0025. Thus, the simulation will run for <math>100000 \times 0.0025 = 250 \text{ unit time}</math>

===Plotting the equations of state===
[[File:Densityaz]]

The above graph displays the density (both ideal and from the simulation) versus temperature.

Note that there is a significant discrepancy between simulation and theory, the theoretical being greater in magnitude. The error bars cannot be seen in the graph as the difference between the ideal and simulated densities is so high that the error bars become too small to visualise.

Factors leading to this discrepancy may involve the model of the Leonard-Jones interaction in the simulation. The ideal gas law relies on there being no inter-particle interaction, and thus if there is a repulsive force due to the Leonard-Jones interactions, then the density may be lower than the ideal calculation.

It can also be seen that as the temperature increases, the discrepancy between the ideal and simulated densities becomes smaller.

==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 short-range order. They also exhibit rotational and to some extent translational motion.

[[File:RDFaz.jpg]]

== Dynamical Properties and the Diffusion Coefficient ==

The Diffusion Coefficients <math>D</math> 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 <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

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:

{| class="wikitable" border="1"
|+ Diffusion Coefficients (unit area/unit time)
! State !! From MSD (thousand atoms) !! From MSD (million) !! From VCAF (thousand atoms) !! From VCAF (million)
|-
| Solid || <math>8.96 \times 10^{-6}</math> || <math>4.5 \times 10^{-6}</math> || <math>-1.66 \times 10^{-4}</math> || <math>4.55\times 10^{-5}</math>
|-
| Liquid || <math>0.085</math> || <math>0.089</math> ||<math>0.098</math> || <math>0.090</math>
|-
| Gas || <math>3.22</math> || <math>2.38</math> || <math>3.35</math> || <math>3.09</math>
|}

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 ===



[[File:MSD1az.jpg|thumb|600px|none]][[File:MSD2az.jpgthumb|600px|none]]

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.

=== Velocity Autocorrelation Function ===

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} v\left(t\right)v\left(t + \tau\right)\mathrm{d}t}{\int_{-\infty}^{\infty} v^2\left(t\right)\mathrm{d}t}</math>

Starting from a 1D Harmonic Oscillator:

<math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

Differentiating gives the velocity function:

<math> v\left(t\right) = -A\omega \sin\left(\omega t + \phi\right)</math>

Substituting into the integrated form of the VAFC:

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} A\omega \sin\left(\omega t + \phi \right) \cdot A \omega \sin \left(\omega \left(t+\tau\right)+\phi\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(A \omega \sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} \sin\left(\omega t + \phi \right) \cdot \sin \left(\omega \left(t+\tau\right)+\phi\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(\sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

Then, using the trigonometric identity

<math>\sin \left(\alpha + \beta \right)= \sin\alpha\cos\beta + \sin\beta\cos\alpha,</math>

<math>\sin\left(\omega t + \phi + \omega\tau\right) = \sin\left(\omega t+\phi\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\omega t + \phi\right)</math>

<math>\therefore C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} \sin\left(\omega t + \phi \right) \cdot \left(\sin\left(\omega t+\phi\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\omega t + \phi\right)\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(\sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

<math>C\left(\tau\right)=\frac{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)\cos\left(\omega \tau\right) + \int_{-\infty}^{\infty} \sin \left(\omega \tau\right) \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

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

<math>C\left(\tau\right)=\frac{\cos\left(\omega \tau\right)\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right) + \sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

<math>\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)</math> can cancel out:

<math>C\left(\tau\right)=\cos\left(\omega \tau\right)+\frac{\sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

Since <math>\sin(x)</math> is an odd function and <math>\cos(x)</math> is an even function, <math>\sin(x)\cos(x)</math> is odd.

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

<math>\therefore \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)=0</math>

<math>\therefore C\left(\tau\right)=\cos\left(\omega \tau\right)+\frac{\sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}=\cos\left(\omega \tau\right)</math>

<math>C\left(\tau\right)=\cos\left(\omega \tau\right)</math>

[[File:VACFaz.png|thumb|800px|Comparison of Normalised VACF with VACF of Liquid and Solid|left]]

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.

2017-02-10T09:40:01Z

Az1114:

File:Energyintroaz.png

2017-02-10T09:39:17Z

Az1114:

2017-02-10T03:55:42Z

Az1114: /* Dynamical Properties and the Diffusion Coefficient */

'''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\left(t\right) = A\cos\left(\omega t + \phi\right) (the values of A, \omega, and \phi 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?
'''

'''TASK: For a single Lennard-Jones interaction, \phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right), 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 (\phi\left(r_{eq}\right)). Evaluate the integrals \int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r, \int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r, and \int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r when \sigma = \epsilon = 1.0.'''

=== Lennard-Jones Potential ===
The below equation calculates the Leonard-Jones potential:

<math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>

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

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

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

<math>\frac{\mathrm{d}\phi}{\mathrm{d}r}=\mathbf{F}=-24\epsilon \left(\frac{2 \sigma^{12}}{r^{13}}-\frac{\sigma^6}{r^7}\right)</math>

<math>r_0=\sigma</math>

<math>\mathbf{F}=-24 \frac{\epsilon}{\sigma}</math>

The equilibrium separation (<math>r_{eq}</math>) 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 <math>\frac{\mathrm{d}\phi}{\mathrm{d}r}</math>, (or <math>\mathbf{F}</math>) to 0.

Hence,

<math>\mathbf{F}=-24\epsilon \left(\frac{2 \sigma^{12}}{r^{13}}-\frac{\sigma^6}{r^7}\right) 0</math>

<math>r_{eq}=\sqrt[6]{2\sigma^6}</math>

Substituting <math>r_{eq}</math> back into <math>\phi(r)</math> gives us the well depth:

<math>\phi(r_{eq})=-\epsilon</math>

----

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

<math>\int 4\epsilon \left(\frac{\sigma^{12}}{r^{12}}-\frac{\sigma^6}{r^6}\right)\mathrm{d}r</math>

<math>=4 \epsilon \left(\frac{\sigma^{12}}{11r^{11}}-\frac{\sigma^6}{5r^5}\right)\mathrm{d}r</math>

When <math>\sigma = \epsilon = 1.0:</math>

<math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=\lim_{r\rightarrow \infty} 4 \left(\frac{1^{12}}{11r^{11}}-\frac{1^6}{5r^5}\right) - 4\left(\frac{1^{12}}{11\times 2^{11}}+\frac{1^6}{5\times 2^5}\right)</math>

<math>=0-0.024822\cdots</math>

<math>\approx -0.0248</math>

Following the same calculations:

<math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=-0.0081767\cdots</math>

<math>\approx -0.00818</math>

<math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>

<math>=-0.0032901\cdots</math>

<math>\approx -0.00329</math>

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

=== 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, let us first calculate the number of molecules in 1ml of water:

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

<math>1.00 \mathrm{ml} = 1.00 \mathrm{g}</math>

<math>18.02 \mathrm{g/mol} = 0.0555 \text{ mols of water}</math>

<math>0.0555 \times 6.022 \times 10^{23} = 3.34 \times 10^{22} \text{ molecules of water.}</math>

On the other hand, 10,000 water molecules is only <math>2.99 \times 10^{-19} \text{ ml of water}</math>

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 direction, 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 <math>(0.5, 0.5, 0.5)</math> in a cube of dimensions <math>(0, 0, 0)</math> to <math>(1, 1, 1)</math> moves along the vector <math>(0.7, 0.6, 0.2)</math>, its final position would be <math>(0.2, 0.1, 0.7)</math>, 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 <math>\sigma=0.34 \mathrm{nm}</math>; <math>\epsilon/k_B=120K</math>.

If we set the cutoff to be <math>r*=3.2,</math>

<math>r=3.2\times0.34=1.088\mathrm{nm}</math>

Well depth: <math>\epsilon=120K\times k_B \times 6.02 \times 10^{23}</math>

<math>\epsilon=0.998 \mathrm{kJ/mol}</math>

and reduced temperature <math>T*=1.5</math> would be <math>T=1.5 \times 120=180\mathrm{K}.</math>

==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:

<math>0.8 \text{ points/volume}</math>

<math>\sqrt[3]{0.8}=0.9283\text{ points/length}</math>

<math>0.9283^{-1}=1.0772 \text{ length/point}</math>

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 <math>0.5928</math> as:

<math>\ \sqrt[3]{4.8^{-1}}=0.5928</math>

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

<pre>region box block 0 10 0 10 0 10</pre>

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:
<pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

The mass command takes two arguments, setting the mass of type 1 atoms as <math>1.0</math> (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 <math>1.0</math> setting the well depth <math>epsilon</math> and the latter sets <math>sigma</math>. 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 <math>\mathbf{x}_i(0)</math> & <math>\mathbf{v}_i\left(0\right)</math>, the Velocity-Verlet integration half-step algorithm can be used.

=== Running the Simulation ===
In the following lines from the input script:
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(100/${timestep})
timestep ${timestep}

### RUN SIMULATION ###
run ${n_steps}
</pre>

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
<pre>
timestep 0.001
run 100000
</pre>
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.

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

[[File:Total Energyaz.png]] [[File:Temperatureaz.png]] [[File:Pressureaz.png]]

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.

[[File:Different Timesteps.png]]

== 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 <math>T,</math> using the following equations:

<math>E_K = \frac{3}{2} N k_B T</math>

<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T \qquad\qquad\qquad(1)</math>

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

<math>\frac{1}{2}\sum_i m_i \left(\gamma v_i\right)^2 = \frac{3}{2} N k_B \mathfrak{T}\qquad\qquad(2)</math>

Thus we can solve for <math>\gamma.</math>:

From <math>(2)</math>:

<math>\gamma^2\sum_i m_i v_i^2 = 3N k_B \mathfrak{T}</math>

Then substituting <math>\sum_i m_i v_i^2</math> from <math>(1)</math>

<math>\gamma^2 3 N k_b T = 3N k_b \mathfrak{T}</math>

<math>\gamma = \sqrt{\frac{\epsilon}{T}}</math>

=== Explaining the Input Script ===

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

<pre>
fix aves all ave/time Nevery Nrepeat Nfreq
</pre>

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:

<pre>
fix aves all ave/time 100 1000 100000
run 100000
</pre>

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.

The timestep chosen for this simulation was 0.0025. Thus, the simulation will run for <math>100000 \times 0.0025 = 250 \text{ unit time}</math>

===Plotting the equations of state===
[[File:Densityaz]]

The above graph displays the density (both ideal and from the simulation) versus temperature.

Note that there is a significant discrepancy between simulation and theory, the theoretical being greater in magnitude. The error bars cannot be seen in the graph as the difference between the ideal and simulated densities is so high that the error bars become too small to visualise.

Factors leading to this discrepancy may involve the model of the Leonard-Jones interaction in the simulation. The ideal gas law relies on there being no inter-particle interaction, and thus if there is a repulsive force due to the Leonard-Jones interactions, then the density may be lower than the ideal calculation.

It can also be seen that as the temperature increases, the discrepancy between the ideal and simulated densities becomes smaller.

==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 short-range order. They also exhibit rotational and to some extent translational motion.

[[File:RDFaz.jpg]]

== Dynamical Properties and the Diffusion Coefficient ==

The Diffusion Coefficients <math>D</math> 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 <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

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:

{| class="wikitable" border="1"
|+ Diffusion Coefficients (unit area/unit time)
! State !! From MSD (thousand atoms) !! From MSD (million) !! From VCAF (thousand atoms) !! From VCAF (million)
|-
| Solid || <math>0.000896</math> || <math>1.5 \times 10^{-5}</math> || <math>-1.66 \times 10^{-5}</math> || <math>4.55\times 10^{-5}</math>
|-
| Liquid || <math>0.0855</math> || <math>0.0892</math> ||<math>0.0979</math> || <math>0.09</math>
|-
| Gas || <math>6.05</math> || <math>3.22</math> || <math>6.32</math> || <math>3.27</math>
|}

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 ===



[[File:MSD1az.jpg|thumb|600px|none]][[File:MSD2az.jpgthumb|600px|none]]

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.

=== Velocity Autocorrelation Function ===

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} v\left(t\right)v\left(t + \tau\right)\mathrm{d}t}{\int_{-\infty}^{\infty} v^2\left(t\right)\mathrm{d}t}</math>

From the 1D Harmonic Oscillator:

<math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

Differentiating this gives the velocity function:

<math> v\left(t\right) = -A\omega \sin\left(\omega t + \phi\right)</math>

Substituting this into the integrated form of the VAFC:

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} A\omega \sin\left(\omega t + \phi \right) \cdot A \omega \sin \left(\omega \left(t+\tau\right)+\phi\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(A \omega \sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

<math>C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} \sin\left(\omega t + \phi \right) \cdot \sin \left(\omega \left(t+\tau\right)+\phi\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(\sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

Then, from the trigonometric identity

<math>\sin \left(\alpha + \beta \right)= \sin\alpha\cos\beta + \sin\beta\cos\alpha,</math>

<math>\sin\left(\omega t + \phi + \omega\tau\right) = \sin\left(\omega t+\phi\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\omega t + \phi\right)</math>

<math>\therefore C\left(\tau\right) = \frac{\int_{-\infty}^{\infty} \sin\left(\omega t + \phi \right) \cdot \left(\sin\left(\omega t+\phi\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\omega t + \phi\right)\right)\mathrm{d}t}{\int_{-\infty}^{\infty} \left(\sin\left(\omega t + \phi\right)\right)^2\mathrm{d}t}</math>

<math>C\left(\tau\right)=\frac{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)\cos\left(\omega \tau\right) + \int_{-\infty}^{\infty} \sin \left(\omega \tau\right) \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

<math>\tau</math> is a constant so <math>\cos\left(\omega\tau\right)</math> and <math>\sin\left(\omega\tau\right)</math> can be taken out of the integral.

<math>C\left(\tau\right)=\frac{\cos\left(\omega \tau\right)\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right) + \sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

<math>\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)</math> can cancel out:

<math>C\left(\tau\right)=\cos\left(\omega \tau\right)+\frac{\sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}</math>

Since <math>\sin(x)</math> is an odd function and <math>\cos(x)</math> is an even function, <math>\sin(x)\cos(x)</math> is an odd function.

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

<math>\therefore \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)=0</math>

<math>\therefore C\left(\tau\right)=\cos\left(\omega \tau\right)+\frac{\sin \left(\omega \tau\right) \int_{-\infty}^{\infty} \cos \left(\omega t + \phi\right)\sin\left(\omega t + \phi\right)}{\int_{-\infty}^{\infty} \sin^2\left(\omega t + \phi\right)}=\cos\left(\omega \tau\right)</math>

<math>C\left(\tau\right)=\cos\left(\omega \tau\right)</math>

[[File:VACF-SHL.png|thumb|800px|Comparison of Normalised VACF with VACF of Liquid and Solid|left]]

We can see from this graph that while the VACF of the liquid quickly decorrelates to 0, the solid still retains some minor fluctuations. This is due to the fact that the atoms in a liquid are more free to move around and collide with each other, "forgetting" what its initial velocity was from exchanging energy, whereas for atoms in the solid state, they are already in a relatively stable position (similar to <math>r_{eq}</math> in L-J interactions) and oscillate about 0. Eventually, atoms in solids, too, will decay, but at a much slower rate than in liquids.

The VACF for the harmonic oscillator does not decay at all since the decorrelation arises from the atoms exchanging energy with their neighbours, which does not occur in harmonic oscillators.

File:RDFaz.jpg

2017-02-10T03:52:29Z

Az1114:

File:MSD2az.jpgthumb

2017-02-10T03:46:08Z

Az1114:

2017-02-06T12:02:06Z

Az1114: /* Equilibration */