Rep:Mod:sfs114

2016-12-15T17:29:42Z

Sfs114: /* Structural Properties and the Radial Distribution Function */

== Theory ==

=== Numerical Integration ===

The classical solution for the position at time <math>t</math> compares well with the velocity-Verlet solution:

[[File:Classicalsolution1.jpg|thumb|centre|800px|A plot of x(t) against t comparing analytical and velocity-Verlet solutions]]

The total energy for the oscillator varies as shown:

[[File:Energy111.jpg|thumb|centre|800px|A plot of energy against time]]

An approximate linear fit has been performed on the maxima of the error of the calculations; the absolute difference between classical and velocity-Verlet solutions. Iterations of using previous results causes error to propagate and increase.

[[File:Errorfit.jpg|thumb|centre|800px|A plot of absolute error, with a linear function fitted to the error maxima]]

The smaller the timestep, the smaller fluctuations in total energy. Calculations over larger timesteps causes a greater error, as particles could end up too close together and face extremely large forces, for example. It is important to monitor the total energy of a physical system to ensure energy is conserved, however infinitesimal timesteps greatly increase time needed to run simulations. Larger timesteps allow a longer length of time to be simulated. A timestep of '''0.028s''' allows energy fluctuations to be as low as 1% and this is a good trade off with simulation run time.

[[File:energytimestep.jpg|thumb|centre|600px|Energy oscillation with a timestep of 0.028s]]
[[File:energytimestep2.jpg|thumb|centre|600px|Energy oscillation with a timestep of 0.5s]]

=== Atomic Forces ===

For a single Lennard-Jones interaction, <math>\phi (r) = 4\epsilon (\frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6})</math>.

:* When potential energy is 0, φ(r) = 0 and r=r0
::If <math>\phi=0, (\frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6}) = 0</math>,
:::<math>\frac{\sigma^{12}}{r^{12}} = \frac{\sigma^{6}}{r^{6}}</math>
:::<math>r_0 = \sigma</math>

:*The force is given by <math>F= \frac{d\phi}{dr}=4\epsilon(-12\sigma^{12}r^{-13}+6\sigma^{6}r^{-7})</math> and at a potential energy of 0, <math>F=\frac{-24\epsilon}{\sigma}</math>.

:*At equilibrium separation <math>r_{eq}</math>, <math>\frac{d\phi}{dr}=0</math>.
::<math>0=4\epsilon(\frac{-12\sigma^{12}}{r^{13}}+\frac{6\sigma^6}{r^7})</math>
::<math>(\frac{12\sigma^{12}}{r^{13}} = \frac{6\sigma^6}{r^7})</math>
::<math>r_{eq}=2^{\frac{1}{6}}\sigma</math>

:*At equilibrium separation, the well depth:
::<math>\phi(r_{eq})=\phi(2^{\frac{1}{6}}\sigma)=-\epsilon</math>

:*<math>\int\phi(r)\, dr=4\epsilon[\frac{-\sigma^{12}}{11r^{11}}+\frac{\sigma^6}{5r^5}]+c</math> and given that <math>\sigma = \epsilon = 1.0</math> so <math>\int\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]+c</math>
::*<math>\int\limits_{2\sigma}^{\infty}\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]_2^{\infty} = -4(\frac{1}{5(2)^5}-\frac{1}{11(2)^{11}}) = -0.02482</math>
::*<math>\int\limits_{2.5\sigma}^{\infty}\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]_{2.5}^{\infty} = -4(\frac{1}{5(2.5)^5}-\frac{1}{11(2.5)^{11}}) = -0.008177</math>
::*<math>\int\limits_{3\sigma}^{\infty}\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]_{3}^{\infty} = -4(\frac{1}{5(3)^5}-\frac{1}{11(3)^{11}}) = -0.00329</math>

=== Periodic Boundary Conditions ===

As <math> pV=Nk_BT </math>, the number of water molecules in 1 mL of water is approximately <math>2.46x10^{19}</math> and 10000 molecules takes up an approximate volume of <math>4.06x10^{-22} m^3</math>.

In a simulation box which runs from (0,0,0) to (1,1,1), an atom that starts at (0.5,0.5,0.5) and moves along vector (0.7,0.6,0.2), will end up at '''(0.2,0.1,0.7)''' once periodic boundary conditions have been applied.

=== Reduced Units ===

The LJ parameters for Argon are: <math>\sigma=0.34nm, \frac{\epsilon}{k_B}=120K</math>.
*<math>r=\sigma \cdot r*= 0.34 \cdot 3.2 = 1.088 nm</math>.
*<math> Well Depth = \phi(r_{eq}) = - \epsilon = - k_B \cdot 120 = -1.656^{-21}J</math>
:<math>-1.656^{-21} \cdot N_A = -997.4 = -0.997 kJ mol^{-1}</math>.
*<math>T=\frac{\epsilon \cdot T*}{k_B}=120 \cdot 1.5 = 180 K</math>.

== Equilibriation ==

=== Creating the Simulation Box ===

If two atoms are generated too close together, the LJ potential shows that the potential between the two would be infinitely large, making force simulations between these two atoms too large to realistically simulate. The LJ cutoff also ensures that LJ potentials are only calculated for atoms that are near enough, and not every other atom in the infinitely repeating lattice, which would greatly increase simulation run time.

A lattice spacing of 1.07722 corresponds to a lattice number density of <math>\frac{1}{1.07722^3}=0.8</math> for a simple cubic lattice. A face centred cubic lattice has 4 lattice points per cell, and thus would require a lattice spacing of <math> (\frac{4}{1.2})^{\frac{1}{3}}= 1.4938</math>. A 10x10x10 box would contain 1000 unit cells, and 4000 lattice points, so the create_atoms command for such a lattice would create 4000 atoms.

=== Setting the Properties of the Atoms ===

The command '''mass 1 1.0''' assigns all atoms of type 1 a mass of 1.0.
The command '''pair_style lj/cut 3.0''' defines the cutoff distance between atoms that have a potential between them to be 3.0 (ie. the simulation does not run for atoms farther apart or closer than this distance).
The command '''pair_coeff ** 1.0 1.0''' specifically defines the pairwise force field coefficients for multiple pairs atoms.

The velocity-Verlet algorithm is the numerical integration method that will be used if <math>x_i(0)</math> and <math>v_i(0)</math> are defined.

=== Running the Simulation ===

Calling upon variables, instead of assigning numbers, makes it much easier to change these variables for every simulation that is run.

=== Checking Equilibriation ===

The simulation takes about 0.3 seconds to equilibriate energy, temperature, and pressure, as shown below:

[[File:energyeq.jpg|thumb|400px|left|Energy equilbriation]]
[[File:energyeqzoom.jpg|thumb|400px|right|Closer look at energy equilibriation]]
[[File:tempeq.jpg|thumb|400px|left|Temperature equilbriation]]
[[File:tempeqzoom.jpg|thumb|400px|right|Closer look at temperature equilibriation]]
[[File:presseq.jpg|thumb|400px|left|Pressure equilbriation]]
[[File:presseqzoom.jpg|thumb|400px|right|Closer look at pressure equilibriation]]

[[File:Equilibriation3.jpg|thumb|800px|centre|Energy equilibriation for multiple timesteps]]
Of the five timsteps used, 0.0025 is the largest acceptable timestep to use as a smaller timestep of 0.01 results in a very similar equilibriation. 0.015 does not equilbriate at all as the time steps are too large for the numerical integration to accurately find an average for the ensemble, and energy drifts; diverging instead of converging to an average value.

== Running Simulations Under Specific Conditions ==

=== Thermostats & Barostats ===

<math> E_K = \frac{3}{2} Nk_B\mathfrak{T} = \frac{1}{2}\sum_{i} m_i (\gamma v_i)^2</math>

:<math> \gamma^2 = \frac{\frac{3}{2}Nk_B\mathfrak{T}}{\frac{1}{2}\sum_{i}m_i v_i^2}</math>

:<math> \gamma^2 = \frac{\frac{3}{2}Nk_B\mathfrak{T}}{\frac{3}{2}Nk_BT}</math>

:<math> \gamma = (\frac{\mathfrak{T}}{T})^\frac{1}{2} </math>

=== Examining the Input Script ===

The command '''fix aves all ave/time 100 1000 100000''' means that values will be sampled every 100 timesteps; in total 1000 readings will be taken to compute a final average on the 100000th timestep.

'''run 100000''' indicates that 100000 timesteps will be simulated.

=== Plotting the Equations of State ===
[[File:download1.jpg|thumb|center|800px|A plot of density against temperature compared to ideal gas law]]

Higher pressures lead to higher densities, both in theory and in these simulations. Our simulated density is higher than that given by the ideal gas law because the simulation takes particle interactions into account. The error increases at higher pressures, when more collisions are likely to occur.

== Calculating Heat Capacities Using Statistical Physics ==
[[File:HEATCAP.jpg|thumb|centre|800px|A plot of heat capacity/volume against temperature]]

Higher pressure results in higher heat capacity as the increased number of molecules per unit volume that can absorb energy to their vibrational excited states. As the simulation is in a lattice, rotational degrees of freedom are not available to the atoms, and so heat capacity decreases as temperature increases, despite expectations.

An example of the input scripts is below:

[[Media:inputnpt.in]]

== Structural Properties and the Radial Distribution Function ==
[[File:RDF1.jpg|thumb|centre|800px|A plot of RDFs]]

The RDF shows the probability of finding a particle at a distance r from a reference particle, relative to an ideal gas. In a gas, there is little order and minimal structure to particles and so the graph has minimal features.

Liquids are slightly more ordered and the decreasing heights of peaks of the RDF correlate to coordination spheres. There is a high probability of finding another particle in a primary coordination sphere but this probability decreases as you go farther away from the reference particle.

The solid FCC lattice has a much higher order, and the RDF peak separation and heights define the lattice structure. The first, second, and third sharp peaks refer to different sets of nearest neighbours, while their heights show how many of those nearest neighbours there are. The lattice spacing is the same as the distance to the second nearest neighbour, 1.475.

[[File:nearestn.jpg|thumb|centre|400px|A diagram of first (left) and second (right) nearest neighbours (shown in blue) with respect to a reference particle (red) in an FCC lattice (other atoms shown in black)]]
[[File:nearestn3.jpg|thumb|centre|250px|A diagram of third nearest neighbours (blue) with respect to a reference particle (red) in an FCC lattice (other atoms shown in black)]]

[[File:Int1.jpg|thumb|centre|600px|A plot of the running integral of the RDF of a solid]]

The coordination numbers are 12 (Int(g(1.205)=12, 12 neighbours), 6 (Int(g(1.475)=18, 6 neighbours), and 24 (Int(g(1.775)=42, 24 neighbours) respectively.

== Dynamical Properties and the Diffusion Coefficient ==

=== Mean Squared Displacement ===
*<math> D_{solid} = \frac{1}{6} \cdot 0.00382752458475 = 6.38x10^{-4}</math>
*<math> D_{liquid} = \frac{1}{6} \cdot 0.509774821123 = 0.085</math>
*<math> D_{gas} = \frac{1}{6} \cdot 7.99193006423 = 1.33</math>

The diffusion coefficient increases as entropy of the phase increases, which matches expectations as gas particles are much more likely to diffuse than a rigid lattice of solid molecules.

[[File:solid.jpg|thumb|centre|Total MSD of simulated solid]]
[[File:liquid1.jpg|thumb|centre|Total MSD of simulated liquid]]
[[File:gas11.jpg|thumb|centre|Total MSD of simulated gas]]

==== 1000000 Atoms ====

[[File:MSD11.jpg|thumb|800px|Total MSD as a function of time for 1000000 atoms]]

*<math> D_{solid} = \frac{1}{6} \cdot 2.79195534196x10^5 = 4.65x10^{-6} </math>
*<math> D_{liquid} = \frac{1}{6} \cdot 0.531614512766 = 0.0886</math>
*<math> D_{gas} = \frac{1}{6} \cdot 18.0968139669 = 3.02</math>

The MSD graph for a gas is curved at first, indicating ballistic motion proportional to <math>T^{2}</math>. After enough collisions have occurred, diffusion is linear, as it is for a liquid which constantly has the same collisions. The diffusion coefficient is close to 0 for solids which is as expected.

=== Velocity Autocorrelation Function ===

::*<math> x(t)=Acos(\omega t + \phi) </math> and <math> v(t) = \frac{dx(t)}{dt} = -A\omega sin (\omega t+ \phi) </math>

<math> C(\tau) = \frac{\int\limits_{-\infty}^{\infty}v(t)-v(t+\tau)dt}{\int\limits_{-\infty}^{\infty}v^2(t)dt} = \frac{\int\limits_{-\infty}^{\infty}sin(\omega t+\phi)sin(\omega(t+\tau)+\phi)}{\int\limits_{-\infty}^{\infty}sin^2(\omega t + \phi)dt} </math>

::*<math> sin(A+B) = sinAcosB + cosAsinB </math>

<math> C(\tau) = \frac{[cos (\omega \tau) \int\limits_{-\infty}^{\infty} sin^2(\omega t + \phi) dt] + [sin (\omega \tau) \int\limits_{-\infty}^{\infty} sin(\omega t + \phi) cos(\omega t + \phi) dt]}{\int\limits_{-\infty}^{\infty}sin^2(\omega t + \phi)dt} </math>

::*<math> sin^2(x) = \frac{1}{2} (1-cos(2x)) </math>

<math> C(\tau) = \frac{[cos (\omega \tau) \int\limits_{-\infty}^{\infty} \frac{1}{2} (1-cos(2(\omega t + \phi)] + [sin (\omega \tau) \int\limits_{-\infty}^{\infty} sin(\omega t + \phi) cos(\omega t + \phi) dt]}{\int\limits_{-\infty}^{\infty} \frac{1}{2} (1-cos(2(\omega t + \phi) dt} </math>

::*<math> \int \frac{1}{2} (1-cos(2(\omega t + \phi)) dt = \frac{t}{2} + \frac{sin(2(\omega t + \phi)}{4\omega} + c </math>
::*<math> \int sin(\omega t + \phi) cos(\omega t + \phi) dt = \frac{sin^2(\omega t + \phi)}{2\omega} + c</math>

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

::*sin(x) is an odd function and integrating between <math>-\infty</math> and <math>\infty</math> will result in 0

<math> C(\tau) = \frac{\frac{tcos(\omega \tau)}{2}}{\frac{t}{2}} = cos (\omega \tau)</math>

VACF minima refer to collisions of particles where velocity is instantaneously 0, negative as they are in the opposite direction to <math>\tau</math>. As VACF is averaged over all molecules, they cancel out once they are out of phase, which happens faster for liquids than it does for solids. In comparison to the harmonic oscillator, which only models one particle without any collisions, no convergence to 0 occurs.

[[File:VACF1.jpg|thumb|centre|800px|A plot of total VACF against <math>\tau</math>]]

Diffusion coefficient estimations, using the trapezium rule:
*<math> D_{solid} = \frac{1}{3} \cdot 0.43184744815700105 = 0.144 </math>
*<math> D_{liquid} = \frac{1}{3} \cdot 146.83331703729999 = 48.9 </math>
*<math> D_{gas} = \frac{1}{3} \cdot 1451.848385 =484 </math>
[[File:runint4.jpg|thumb|centre|Plot of running integral of VACF of a solid]]
[[File:runint5.jpg|thumb|centre|Plot of running integral of VACF of a liquid]]
[[File:runint6.jpg|thumb|centre|Plot of running integral of VACF of a gas]]

==== 1000000 Atoms ====
[[File:VACF1111.jpg|thumb|centre|800px|A plot of total VACF against <math>\tau</math> for 1000000 atoms]]

*<math> D_{solid} = \frac{1}{3} \cdot 0.41626945400000459 = 0.139</math>
*<math> D_{liquid} = \frac{1}{3} \cdot 123.7270701106= 41.2 </math>
*<math> D_{gas} = \frac{1}{3} \cdot 1466.443215 = 489 </math>
[[File:runint7.jpg|thumb|centre|Plot of running integral of VACF of a solid]]
[[File:runint8.jpg|thumb|centre|Plot of running integral of VACF of a solid]]
[[File:runint9.jpg|thumb|centre|Plot of running integral of VACF of a solid]]

The estimated diffusion coefficients for the two simulations follow the same trend and agree quite well, however the values obtained for the larger 1000000 atom simulations are generally larger. The largest sources of error include the trapezium rules used to calculate the integral and the simulation assumption that velocities do not change upon collisions.

Rep:Mod:sfs114

2016-12-15T17:22:23Z

Sfs114: /* Periodic Boundary Conditions */

== Theory ==

=== Numerical Integration ===

The classical solution for the position at time <math>t</math> compares well with the velocity-Verlet solution:

[[File:Classicalsolution1.jpg|thumb|centre|800px|A plot of x(t) against t comparing analytical and velocity-Verlet solutions]]

The total energy for the oscillator varies as shown:

[[File:Energy111.jpg|thumb|centre|800px|A plot of energy against time]]

An approximate linear fit has been performed on the maxima of the error of the calculations; the absolute difference between classical and velocity-Verlet solutions. Iterations of using previous results causes error to propagate and increase.

[[File:Errorfit.jpg|thumb|centre|800px|A plot of absolute error, with a linear function fitted to the error maxima]]

The smaller the timestep, the smaller fluctuations in total energy. Calculations over larger timesteps causes a greater error, as particles could end up too close together and face extremely large forces, for example. It is important to monitor the total energy of a physical system to ensure energy is conserved, however infinitesimal timesteps greatly increase time needed to run simulations. Larger timesteps allow a longer length of time to be simulated. A timestep of '''0.028s''' allows energy fluctuations to be as low as 1% and this is a good trade off with simulation run time.

[[File:energytimestep.jpg|thumb|centre|600px|Energy oscillation with a timestep of 0.028s]]
[[File:energytimestep2.jpg|thumb|centre|600px|Energy oscillation with a timestep of 0.5s]]

=== Atomic Forces ===

For a single Lennard-Jones interaction, <math>\phi (r) = 4\epsilon (\frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6})</math>.

:* When potential energy is 0, φ(r) = 0 and r=r0
::If <math>\phi=0, (\frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6}) = 0</math>,
:::<math>\frac{\sigma^{12}}{r^{12}} = \frac{\sigma^{6}}{r^{6}}</math>
:::<math>r_0 = \sigma</math>

:*The force is given by <math>F= \frac{d\phi}{dr}=4\epsilon(-12\sigma^{12}r^{-13}+6\sigma^{6}r^{-7})</math> and at a potential energy of 0, <math>F=\frac{-24\epsilon}{\sigma}</math>.

:*At equilibrium separation <math>r_{eq}</math>, <math>\frac{d\phi}{dr}=0</math>.
::<math>0=4\epsilon(\frac{-12\sigma^{12}}{r^{13}}+\frac{6\sigma^6}{r^7})</math>
::<math>(\frac{12\sigma^{12}}{r^{13}} = \frac{6\sigma^6}{r^7})</math>
::<math>r_{eq}=2^{\frac{1}{6}}\sigma</math>

:*At equilibrium separation, the well depth:
::<math>\phi(r_{eq})=\phi(2^{\frac{1}{6}}\sigma)=-\epsilon</math>

:*<math>\int\phi(r)\, dr=4\epsilon[\frac{-\sigma^{12}}{11r^{11}}+\frac{\sigma^6}{5r^5}]+c</math> and given that <math>\sigma = \epsilon = 1.0</math> so <math>\int\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]+c</math>
::*<math>\int\limits_{2\sigma}^{\infty}\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]_2^{\infty} = -4(\frac{1}{5(2)^5}-\frac{1}{11(2)^{11}}) = -0.02482</math>
::*<math>\int\limits_{2.5\sigma}^{\infty}\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]_{2.5}^{\infty} = -4(\frac{1}{5(2.5)^5}-\frac{1}{11(2.5)^{11}}) = -0.008177</math>
::*<math>\int\limits_{3\sigma}^{\infty}\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]_{3}^{\infty} = -4(\frac{1}{5(3)^5}-\frac{1}{11(3)^{11}}) = -0.00329</math>

=== Periodic Boundary Conditions ===

As <math> pV=Nk_BT </math>, the number of water molecules in 1 mL of water is approximately <math>2.46x10^{19}</math> and 10000 molecules takes up an approximate volume of <math>4.06x10^{-22} m^3</math>.

In a simulation box which runs from (0,0,0) to (1,1,1), an atom that starts at (0.5,0.5,0.5) and moves along vector (0.7,0.6,0.2), will end up at '''(0.2,0.1,0.7)''' once periodic boundary conditions have been applied.

=== Reduced Units ===

The LJ parameters for Argon are: <math>\sigma=0.34nm, \frac{\epsilon}{k_B}=120K</math>.
*<math>r=\sigma \cdot r*= 0.34 \cdot 3.2 = 1.088 nm</math>.
*<math> Well Depth = \phi(r_{eq}) = - \epsilon = - k_B \cdot 120 = -1.656^{-21}J</math>
:<math>-1.656^{-21} \cdot N_A = -997.4 = -0.997 kJ mol^{-1}</math>.
*<math>T=\frac{\epsilon \cdot T*}{k_B}=120 \cdot 1.5 = 180 K</math>.

== Equilibriation ==

=== Creating the Simulation Box ===

If two atoms are generated too close together, the LJ potential shows that the potential between the two would be infinitely large, making force simulations between these two atoms too large to realistically simulate. The LJ cutoff also ensures that LJ potentials are only calculated for atoms that are near enough, and not every other atom in the infinitely repeating lattice, which would greatly increase simulation run time.

A lattice spacing of 1.07722 corresponds to a lattice number density of <math>\frac{1}{1.07722^3}=0.8</math> for a simple cubic lattice. A face centred cubic lattice has 4 lattice points per cell, and thus would require a lattice spacing of <math> (\frac{4}{1.2})^{\frac{1}{3}}= 1.4938</math>. A 10x10x10 box would contain 1000 unit cells, and 4000 lattice points, so the create_atoms command for such a lattice would create 4000 atoms.

=== Setting the Properties of the Atoms ===

The command '''mass 1 1.0''' assigns all atoms of type 1 a mass of 1.0.
The command '''pair_style lj/cut 3.0''' defines the cutoff distance between atoms that have a potential between them to be 3.0 (ie. the simulation does not run for atoms farther apart or closer than this distance).
The command '''pair_coeff ** 1.0 1.0''' specifically defines the pairwise force field coefficients for multiple pairs atoms.

The velocity-Verlet algorithm is the numerical integration method that will be used if <math>x_i(0)</math> and <math>v_i(0)</math> are defined.

=== Running the Simulation ===

Calling upon variables, instead of assigning numbers, makes it much easier to change these variables for every simulation that is run.

=== Checking Equilibriation ===

The simulation takes about 0.3 seconds to equilibriate energy, temperature, and pressure, as shown below:

[[File:energyeq.jpg|thumb|400px|left|Energy equilbriation]]
[[File:energyeqzoom.jpg|thumb|400px|right|Closer look at energy equilibriation]]
[[File:tempeq.jpg|thumb|400px|left|Temperature equilbriation]]
[[File:tempeqzoom.jpg|thumb|400px|right|Closer look at temperature equilibriation]]
[[File:presseq.jpg|thumb|400px|left|Pressure equilbriation]]
[[File:presseqzoom.jpg|thumb|400px|right|Closer look at pressure equilibriation]]

[[File:Equilibriation3.jpg|thumb|800px|centre|Energy equilibriation for multiple timesteps]]
Of the five timsteps used, 0.0025 is the largest acceptable timestep to use as a smaller timestep of 0.01 results in a very similar equilibriation. 0.015 does not equilbriate at all as the time steps are too large for the numerical integration to accurately find an average for the ensemble, and energy drifts; diverging instead of converging to an average value.

== Running Simulations Under Specific Conditions ==

=== Thermostats & Barostats ===

<math> E_K = \frac{3}{2} Nk_B\mathfrak{T} = \frac{1}{2}\sum_{i} m_i (\gamma v_i)^2</math>

:<math> \gamma^2 = \frac{\frac{3}{2}Nk_B\mathfrak{T}}{\frac{1}{2}\sum_{i}m_i v_i^2}</math>

:<math> \gamma^2 = \frac{\frac{3}{2}Nk_B\mathfrak{T}}{\frac{3}{2}Nk_BT}</math>

:<math> \gamma = (\frac{\mathfrak{T}}{T})^\frac{1}{2} </math>

=== Examining the Input Script ===

The command '''fix aves all ave/time 100 1000 100000''' means that values will be sampled every 100 timesteps; in total 1000 readings will be taken to compute a final average on the 100000th timestep.

'''run 100000''' indicates that 100000 timesteps will be simulated.

=== Plotting the Equations of State ===
[[File:download1.jpg|thumb|center|800px|A plot of density against temperature compared to ideal gas law]]

Higher pressures lead to higher densities, both in theory and in these simulations. Our simulated density is higher than that given by the ideal gas law because the simulation takes particle interactions into account. The error increases at higher pressures, when more collisions are likely to occur.

== Calculating Heat Capacities Using Statistical Physics ==
[[File:HEATCAP.jpg|thumb|centre|800px|A plot of heat capacity/volume against temperature]]

Higher pressure results in higher heat capacity as the increased number of molecules per unit volume that can absorb energy to their vibrational excited states. As the simulation is in a lattice, rotational degrees of freedom are not available to the atoms, and so heat capacity decreases as temperature increases, despite expectations.

An example of the input scripts is below:

[[Media:inputnpt.in]]

== Structural Properties and the Radial Distribution Function ==
[[File:RDF1.jpg|thumb|centre|800px|A plot of RDFs]]

The RDF shows the probability of finding a particle at a distance r from a reference particle, relative to an ideal gas. In a gas, there is little order and minimal structure to particles and so the graph has minimal features.

Liquids are slightly more ordered and the decreasing heights of peaks of the RDF correlate to coordination spheres. There is a high probability of finding another particle in a primary coordination sphere but this probability decreases as you go farther away from the reference particle.

The solid FCC lattice has a much higher order, and the RDF peak separation and heights define the lattice structure. The first, second, and third sharp peaks refer to different sets of nearest neighbours, while their heights show how many of those nearest neighbours there are. The lattice spacing is the same as the distance to the second nearest neighbour, 1.475.

[[File:nearestn.jpg|thumb|centre|400px|A diagram of first (left) and second (right) nearest neighbours (shown in blue) with respect to a reference particle (red) in an FCC lattice (other atoms shown in black)]]
[[File:nearestn2.jpg|thumb|centre|250px|A diagram of third nearest neighbours (blue) with respect to a reference particle (red) in an FCC lattice (other atoms shown in black)]]

[[File:Int1.jpg|thumb|centre|600px|A plot of the running integral of the RDF of a solid]]

The coordination numbers are 12 (Int(g(1.205)=12, 12 neighbours), 6 (Int(g(1.475)=18, 6 neighbours), and 24 (Int(g(1.775)=42, 24 neighbours) respectively.

== Dynamical Properties and the Diffusion Coefficient ==

=== Mean Squared Displacement ===
*<math> D_{solid} = \frac{1}{6} \cdot 0.00382752458475 = 6.38x10^{-4}</math>
*<math> D_{liquid} = \frac{1}{6} \cdot 0.509774821123 = 0.085</math>
*<math> D_{gas} = \frac{1}{6} \cdot 7.99193006423 = 1.33</math>

The diffusion coefficient increases as entropy of the phase increases, which matches expectations as gas particles are much more likely to diffuse than a rigid lattice of solid molecules.

[[File:solid.jpg|thumb|centre|Total MSD of simulated solid]]
[[File:liquid1.jpg|thumb|centre|Total MSD of simulated liquid]]
[[File:gas11.jpg|thumb|centre|Total MSD of simulated gas]]

==== 1000000 Atoms ====

[[File:MSD11.jpg|thumb|800px|Total MSD as a function of time for 1000000 atoms]]

*<math> D_{solid} = \frac{1}{6} \cdot 2.79195534196x10^5 = 4.65x10^{-6} </math>
*<math> D_{liquid} = \frac{1}{6} \cdot 0.531614512766 = 0.0886</math>
*<math> D_{gas} = \frac{1}{6} \cdot 18.0968139669 = 3.02</math>

The MSD graph for a gas is curved at first, indicating ballistic motion proportional to <math>T^{2}</math>. After enough collisions have occurred, diffusion is linear, as it is for a liquid which constantly has the same collisions. The diffusion coefficient is close to 0 for solids which is as expected.

=== Velocity Autocorrelation Function ===

::*<math> x(t)=Acos(\omega t + \phi) </math> and <math> v(t) = \frac{dx(t)}{dt} = -A\omega sin (\omega t+ \phi) </math>

<math> C(\tau) = \frac{\int\limits_{-\infty}^{\infty}v(t)-v(t+\tau)dt}{\int\limits_{-\infty}^{\infty}v^2(t)dt} = \frac{\int\limits_{-\infty}^{\infty}sin(\omega t+\phi)sin(\omega(t+\tau)+\phi)}{\int\limits_{-\infty}^{\infty}sin^2(\omega t + \phi)dt} </math>

::*<math> sin(A+B) = sinAcosB + cosAsinB </math>

<math> C(\tau) = \frac{[cos (\omega \tau) \int\limits_{-\infty}^{\infty} sin^2(\omega t + \phi) dt] + [sin (\omega \tau) \int\limits_{-\infty}^{\infty} sin(\omega t + \phi) cos(\omega t + \phi) dt]}{\int\limits_{-\infty}^{\infty}sin^2(\omega t + \phi)dt} </math>

::*<math> sin^2(x) = \frac{1}{2} (1-cos(2x)) </math>

<math> C(\tau) = \frac{[cos (\omega \tau) \int\limits_{-\infty}^{\infty} \frac{1}{2} (1-cos(2(\omega t + \phi)] + [sin (\omega \tau) \int\limits_{-\infty}^{\infty} sin(\omega t + \phi) cos(\omega t + \phi) dt]}{\int\limits_{-\infty}^{\infty} \frac{1}{2} (1-cos(2(\omega t + \phi) dt} </math>

::*<math> \int \frac{1}{2} (1-cos(2(\omega t + \phi)) dt = \frac{t}{2} + \frac{sin(2(\omega t + \phi)}{4\omega} + c </math>
::*<math> \int sin(\omega t + \phi) cos(\omega t + \phi) dt = \frac{sin^2(\omega t + \phi)}{2\omega} + c</math>

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

::*sin(x) is an odd function and integrating between <math>-\infty</math> and <math>\infty</math> will result in 0

<math> C(\tau) = \frac{\frac{tcos(\omega \tau)}{2}}{\frac{t}{2}} = cos (\omega \tau)</math>

VACF minima refer to collisions of particles where velocity is instantaneously 0, negative as they are in the opposite direction to <math>\tau</math>. As VACF is averaged over all molecules, they cancel out once they are out of phase, which happens faster for liquids than it does for solids. In comparison to the harmonic oscillator, which only models one particle without any collisions, no convergence to 0 occurs.

[[File:VACF1.jpg|thumb|centre|800px|A plot of total VACF against <math>\tau</math>]]

Diffusion coefficient estimations, using the trapezium rule:
*<math> D_{solid} = \frac{1}{3} \cdot 0.43184744815700105 = 0.144 </math>
*<math> D_{liquid} = \frac{1}{3} \cdot 146.83331703729999 = 48.9 </math>
*<math> D_{gas} = \frac{1}{3} \cdot 1451.848385 =484 </math>
[[File:runint4.jpg|thumb|centre|Plot of running integral of VACF of a solid]]
[[File:runint5.jpg|thumb|centre|Plot of running integral of VACF of a liquid]]
[[File:runint6.jpg|thumb|centre|Plot of running integral of VACF of a gas]]

==== 1000000 Atoms ====
[[File:VACF1111.jpg|thumb|centre|800px|A plot of total VACF against <math>\tau</math> for 1000000 atoms]]

*<math> D_{solid} = \frac{1}{3} \cdot 0.41626945400000459 = 0.139</math>
*<math> D_{liquid} = \frac{1}{3} \cdot 123.7270701106= 41.2 </math>
*<math> D_{gas} = \frac{1}{3} \cdot 1466.443215 = 489 </math>
[[File:runint7.jpg|thumb|centre|Plot of running integral of VACF of a solid]]
[[File:runint8.jpg|thumb|centre|Plot of running integral of VACF of a solid]]
[[File:runint9.jpg|thumb|centre|Plot of running integral of VACF of a solid]]

The estimated diffusion coefficients for the two simulations follow the same trend and agree quite well, however the values obtained for the larger 1000000 atom simulations are generally larger. The largest sources of error include the trapezium rules used to calculate the integral and the simulation assumption that velocities do not change upon collisions.

Rep:Mod:sfs114

2016-12-15T17:19:07Z

Sfs114: /* 1000000 Atoms */

== Theory ==

=== Numerical Integration ===

The classical solution for the position at time <math>t</math> compares well with the velocity-Verlet solution:

[[File:Classicalsolution1.jpg|thumb|centre|800px|A plot of x(t) against t comparing analytical and velocity-Verlet solutions]]

The total energy for the oscillator varies as shown:

[[File:Energy111.jpg|thumb|centre|800px|A plot of energy against time]]

An approximate linear fit has been performed on the maxima of the error of the calculations; the absolute difference between classical and velocity-Verlet solutions. Iterations of using previous results causes error to propagate and increase.

[[File:Errorfit.jpg|thumb|centre|800px|A plot of absolute error, with a linear function fitted to the error maxima]]

The smaller the timestep, the smaller fluctuations in total energy. Calculations over larger timesteps causes a greater error, as particles could end up too close together and face extremely large forces, for example. It is important to monitor the total energy of a physical system to ensure energy is conserved, however infinitesimal timesteps greatly increase time needed to run simulations. Larger timesteps allow a longer length of time to be simulated. A timestep of '''0.028s''' allows energy fluctuations to be as low as 1% and this is a good trade off with simulation run time.

[[File:energytimestep.jpg|thumb|centre|600px|Energy oscillation with a timestep of 0.028s]]
[[File:energytimestep2.jpg|thumb|centre|600px|Energy oscillation with a timestep of 0.5s]]

=== Atomic Forces ===

For a single Lennard-Jones interaction, <math>\phi (r) = 4\epsilon (\frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6})</math>.

:* When potential energy is 0, φ(r) = 0 and r=r0
::If <math>\phi=0, (\frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6}) = 0</math>,
:::<math>\frac{\sigma^{12}}{r^{12}} = \frac{\sigma^{6}}{r^{6}}</math>
:::<math>r_0 = \sigma</math>

:*The force is given by <math>F= \frac{d\phi}{dr}=4\epsilon(-12\sigma^{12}r^{-13}+6\sigma^{6}r^{-7})</math> and at a potential energy of 0, <math>F=\frac{-24\epsilon}{\sigma}</math>.

:*At equilibrium separation <math>r_{eq}</math>, <math>\frac{d\phi}{dr}=0</math>.
::<math>0=4\epsilon(\frac{-12\sigma^{12}}{r^{13}}+\frac{6\sigma^6}{r^7})</math>
::<math>(\frac{12\sigma^{12}}{r^{13}} = \frac{6\sigma^6}{r^7})</math>
::<math>r_{eq}=2^{\frac{1}{6}}\sigma</math>

:*At equilibrium separation, the well depth:
::<math>\phi(r_{eq})=\phi(2^{\frac{1}{6}}\sigma)=-\epsilon</math>

:*<math>\int\phi(r)\, dr=4\epsilon[\frac{-\sigma^{12}}{11r^{11}}+\frac{\sigma^6}{5r^5}]+c</math> and given that <math>\sigma = \epsilon = 1.0</math> so <math>\int\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]+c</math>
::*<math>\int\limits_{2\sigma}^{\infty}\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]_2^{\infty} = -4(\frac{1}{5(2)^5}-\frac{1}{11(2)^{11}}) = -0.02482</math>
::*<math>\int\limits_{2.5\sigma}^{\infty}\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]_{2.5}^{\infty} = -4(\frac{1}{5(2.5)^5}-\frac{1}{11(2.5)^{11}}) = -0.008177</math>
::*<math>\int\limits_{3\sigma}^{\infty}\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]_{3}^{\infty} = -4(\frac{1}{5(3)^5}-\frac{1}{11(3)^{11}}) = -0.00329</math>

=== Periodic Boundary Conditions ===

As <math> pV=Nk_BT </math>, the number of water molecules in 1 mL of water is approximately <math>2.46x10^{19}</math> and 10000 molecules takes up an approximate volume of <math>4.06x10^{-22} m^3</math>.

In a simulation box which runs from (0,0,0) to (1,1,1), an atom that starts at (0.5,0.5,0.5) and moves along vector (0.7,0.6,0.2), will end up at (0.2,0.1,0.7) once periodic boundary conditions have been applied.

=== Reduced Units ===

The LJ parameters for Argon are: <math>\sigma=0.34nm, \frac{\epsilon}{k_B}=120K</math>.
*<math>r=\sigma \cdot r*= 0.34 \cdot 3.2 = 1.088 nm</math>.
*<math> Well Depth = \phi(r_{eq}) = - \epsilon = - k_B \cdot 120 = -1.656^{-21}J</math>
:<math>-1.656^{-21} \cdot N_A = -997.4 = -0.997 kJ mol^{-1}</math>.
*<math>T=\frac{\epsilon \cdot T*}{k_B}=120 \cdot 1.5 = 180 K</math>.

== Equilibriation ==

=== Creating the Simulation Box ===

If two atoms are generated too close together, the LJ potential shows that the potential between the two would be infinitely large, making force simulations between these two atoms too large to realistically simulate. The LJ cutoff also ensures that LJ potentials are only calculated for atoms that are near enough, and not every other atom in the infinitely repeating lattice, which would greatly increase simulation run time.

A lattice spacing of 1.07722 corresponds to a lattice number density of <math>\frac{1}{1.07722^3}=0.8</math> for a simple cubic lattice. A face centred cubic lattice has 4 lattice points per cell, and thus would require a lattice spacing of <math> (\frac{4}{1.2})^{\frac{1}{3}}= 1.4938</math>. A 10x10x10 box would contain 1000 unit cells, and 4000 lattice points, so the create_atoms command for such a lattice would create 4000 atoms.

=== Setting the Properties of the Atoms ===

The command '''mass 1 1.0''' assigns all atoms of type 1 a mass of 1.0.
The command '''pair_style lj/cut 3.0''' defines the cutoff distance between atoms that have a potential between them to be 3.0 (ie. the simulation does not run for atoms farther apart or closer than this distance).
The command '''pair_coeff ** 1.0 1.0''' specifically defines the pairwise force field coefficients for multiple pairs atoms.

The velocity-Verlet algorithm is the numerical integration method that will be used if <math>x_i(0)</math> and <math>v_i(0)</math> are defined.

=== Running the Simulation ===

Calling upon variables, instead of assigning numbers, makes it much easier to change these variables for every simulation that is run.

=== Checking Equilibriation ===

The simulation takes about 0.3 seconds to equilibriate energy, temperature, and pressure, as shown below:

[[File:energyeq.jpg|thumb|400px|left|Energy equilbriation]]
[[File:energyeqzoom.jpg|thumb|400px|right|Closer look at energy equilibriation]]
[[File:tempeq.jpg|thumb|400px|left|Temperature equilbriation]]
[[File:tempeqzoom.jpg|thumb|400px|right|Closer look at temperature equilibriation]]
[[File:presseq.jpg|thumb|400px|left|Pressure equilbriation]]
[[File:presseqzoom.jpg|thumb|400px|right|Closer look at pressure equilibriation]]

[[File:Equilibriation3.jpg|thumb|800px|centre|Energy equilibriation for multiple timesteps]]
Of the five timsteps used, 0.0025 is the largest acceptable timestep to use as a smaller timestep of 0.01 results in a very similar equilibriation. 0.015 does not equilbriate at all as the time steps are too large for the numerical integration to accurately find an average for the ensemble, and energy drifts; diverging instead of converging to an average value.

== Running Simulations Under Specific Conditions ==

=== Thermostats & Barostats ===

<math> E_K = \frac{3}{2} Nk_B\mathfrak{T} = \frac{1}{2}\sum_{i} m_i (\gamma v_i)^2</math>

:<math> \gamma^2 = \frac{\frac{3}{2}Nk_B\mathfrak{T}}{\frac{1}{2}\sum_{i}m_i v_i^2}</math>

:<math> \gamma^2 = \frac{\frac{3}{2}Nk_B\mathfrak{T}}{\frac{3}{2}Nk_BT}</math>

:<math> \gamma = (\frac{\mathfrak{T}}{T})^\frac{1}{2} </math>

=== Examining the Input Script ===

The command '''fix aves all ave/time 100 1000 100000''' means that values will be sampled every 100 timesteps; in total 1000 readings will be taken to compute a final average on the 100000th timestep.

'''run 100000''' indicates that 100000 timesteps will be simulated.

=== Plotting the Equations of State ===
[[File:download1.jpg|thumb|center|800px|A plot of density against temperature compared to ideal gas law]]

Higher pressures lead to higher densities, both in theory and in these simulations. Our simulated density is higher than that given by the ideal gas law because the simulation takes particle interactions into account. The error increases at higher pressures, when more collisions are likely to occur.

== Calculating Heat Capacities Using Statistical Physics ==
[[File:HEATCAP.jpg|thumb|centre|800px|A plot of heat capacity/volume against temperature]]

Higher pressure results in higher heat capacity as the increased number of molecules per unit volume that can absorb energy to their vibrational excited states. As the simulation is in a lattice, rotational degrees of freedom are not available to the atoms, and so heat capacity decreases as temperature increases, despite expectations.

An example of the input scripts is below:

[[Media:inputnpt.in]]

== Structural Properties and the Radial Distribution Function ==
[[File:RDF1.jpg|thumb|centre|800px|A plot of RDFs]]

The RDF shows the probability of finding a particle at a distance r from a reference particle, relative to an ideal gas. In a gas, there is little order and minimal structure to particles and so the graph has minimal features.

Liquids are slightly more ordered and the decreasing heights of peaks of the RDF correlate to coordination spheres. There is a high probability of finding another particle in a primary coordination sphere but this probability decreases as you go farther away from the reference particle.

The solid FCC lattice has a much higher order, and the RDF peak separation and heights define the lattice structure. The first, second, and third sharp peaks refer to different sets of nearest neighbours, while their heights show how many of those nearest neighbours there are. The lattice spacing is the same as the distance to the second nearest neighbour, 1.475.

[[File:nearestn.jpg|thumb|centre|400px|A diagram of first (left) and second (right) nearest neighbours (shown in blue) with respect to a reference particle (red) in an FCC lattice (other atoms shown in black)]]
[[File:nearestn2.jpg|thumb|centre|250px|A diagram of third nearest neighbours (blue) with respect to a reference particle (red) in an FCC lattice (other atoms shown in black)]]

[[File:Int1.jpg|thumb|centre|600px|A plot of the running integral of the RDF of a solid]]

The coordination numbers are 12 (Int(g(1.205)=12, 12 neighbours), 6 (Int(g(1.475)=18, 6 neighbours), and 24 (Int(g(1.775)=42, 24 neighbours) respectively.

== Dynamical Properties and the Diffusion Coefficient ==

=== Mean Squared Displacement ===
*<math> D_{solid} = \frac{1}{6} \cdot 0.00382752458475 = 6.38x10^{-4}</math>
*<math> D_{liquid} = \frac{1}{6} \cdot 0.509774821123 = 0.085</math>
*<math> D_{gas} = \frac{1}{6} \cdot 7.99193006423 = 1.33</math>

The diffusion coefficient increases as entropy of the phase increases, which matches expectations as gas particles are much more likely to diffuse than a rigid lattice of solid molecules.

[[File:solid.jpg|thumb|centre|Total MSD of simulated solid]]
[[File:liquid1.jpg|thumb|centre|Total MSD of simulated liquid]]
[[File:gas11.jpg|thumb|centre|Total MSD of simulated gas]]

==== 1000000 Atoms ====

[[File:MSD11.jpg|thumb|800px|Total MSD as a function of time for 1000000 atoms]]

*<math> D_{solid} = \frac{1}{6} \cdot 2.79195534196x10^5 = 4.65x10^{-6} </math>
*<math> D_{liquid} = \frac{1}{6} \cdot 0.531614512766 = 0.0886</math>
*<math> D_{gas} = \frac{1}{6} \cdot 18.0968139669 = 3.02</math>

The MSD graph for a gas is curved at first, indicating ballistic motion proportional to <math>T^{2}</math>. After enough collisions have occurred, diffusion is linear, as it is for a liquid which constantly has the same collisions. The diffusion coefficient is close to 0 for solids which is as expected.

=== Velocity Autocorrelation Function ===

::*<math> x(t)=Acos(\omega t + \phi) </math> and <math> v(t) = \frac{dx(t)}{dt} = -A\omega sin (\omega t+ \phi) </math>

<math> C(\tau) = \frac{\int\limits_{-\infty}^{\infty}v(t)-v(t+\tau)dt}{\int\limits_{-\infty}^{\infty}v^2(t)dt} = \frac{\int\limits_{-\infty}^{\infty}sin(\omega t+\phi)sin(\omega(t+\tau)+\phi)}{\int\limits_{-\infty}^{\infty}sin^2(\omega t + \phi)dt} </math>

::*<math> sin(A+B) = sinAcosB + cosAsinB </math>

<math> C(\tau) = \frac{[cos (\omega \tau) \int\limits_{-\infty}^{\infty} sin^2(\omega t + \phi) dt] + [sin (\omega \tau) \int\limits_{-\infty}^{\infty} sin(\omega t + \phi) cos(\omega t + \phi) dt]}{\int\limits_{-\infty}^{\infty}sin^2(\omega t + \phi)dt} </math>

::*<math> sin^2(x) = \frac{1}{2} (1-cos(2x)) </math>

<math> C(\tau) = \frac{[cos (\omega \tau) \int\limits_{-\infty}^{\infty} \frac{1}{2} (1-cos(2(\omega t + \phi)] + [sin (\omega \tau) \int\limits_{-\infty}^{\infty} sin(\omega t + \phi) cos(\omega t + \phi) dt]}{\int\limits_{-\infty}^{\infty} \frac{1}{2} (1-cos(2(\omega t + \phi) dt} </math>

::*<math> \int \frac{1}{2} (1-cos(2(\omega t + \phi)) dt = \frac{t}{2} + \frac{sin(2(\omega t + \phi)}{4\omega} + c </math>
::*<math> \int sin(\omega t + \phi) cos(\omega t + \phi) dt = \frac{sin^2(\omega t + \phi)}{2\omega} + c</math>

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

::*sin(x) is an odd function and integrating between <math>-\infty</math> and <math>\infty</math> will result in 0

<math> C(\tau) = \frac{\frac{tcos(\omega \tau)}{2}}{\frac{t}{2}} = cos (\omega \tau)</math>

VACF minima refer to collisions of particles where velocity is instantaneously 0, negative as they are in the opposite direction to <math>\tau</math>. As VACF is averaged over all molecules, they cancel out once they are out of phase, which happens faster for liquids than it does for solids. In comparison to the harmonic oscillator, which only models one particle without any collisions, no convergence to 0 occurs.

[[File:VACF1.jpg|thumb|centre|800px|A plot of total VACF against <math>\tau</math>]]

Diffusion coefficient estimations, using the trapezium rule:
*<math> D_{solid} = \frac{1}{3} \cdot 0.43184744815700105 = 0.144 </math>
*<math> D_{liquid} = \frac{1}{3} \cdot 146.83331703729999 = 48.9 </math>
*<math> D_{gas} = \frac{1}{3} \cdot 1451.848385 =484 </math>
[[File:runint4.jpg|thumb|centre|Plot of running integral of VACF of a solid]]
[[File:runint5.jpg|thumb|centre|Plot of running integral of VACF of a liquid]]
[[File:runint6.jpg|thumb|centre|Plot of running integral of VACF of a gas]]

==== 1000000 Atoms ====
[[File:VACF1111.jpg|thumb|centre|800px|A plot of total VACF against <math>\tau</math> for 1000000 atoms]]

*<math> D_{solid} = \frac{1}{3} \cdot 0.41626945400000459 = 0.139</math>
*<math> D_{liquid} = \frac{1}{3} \cdot 123.7270701106= 41.2 </math>
*<math> D_{gas} = \frac{1}{3} \cdot 1466.443215 = 489 </math>
[[File:runint7.jpg|thumb|centre|Plot of running integral of VACF of a solid]]
[[File:runint8.jpg|thumb|centre|Plot of running integral of VACF of a solid]]
[[File:runint9.jpg|thumb|centre|Plot of running integral of VACF of a solid]]

The estimated diffusion coefficients for the two simulations follow the same trend and agree quite well, however the values obtained for the larger 1000000 atom simulations are generally larger. The largest sources of error include the trapezium rules used to calculate the integral and the simulation assumption that velocities do not change upon collisions.

File:Runint9.jpg

2016-12-15T17:18:06Z

Sfs114:

File:Runint8.jpg

2016-12-15T17:17:42Z

Sfs114:

2016-12-15T15:24:37Z

Sfs114:

2016-12-14T19:49:41Z

Sfs114: Sfs114 uploaded a new version of File:Runint3.jpg

File:Runint2.jpg

2016-12-14T19:49:06Z

Sfs114: Sfs114 uploaded a new version of File:Runint2.jpg

File:Runint4.jpg

2016-12-14T19:48:39Z

Sfs114: Sfs114 uploaded a new version of File:Runint4.jpg

Rep:Mod:sfs114

2016-12-14T19:48:12Z

Sfs114: /* 1000000 Atoms */

== Theory ==

=== Numerical Integration ===

The classical solution for the position at time <math>t</math> compares well with the velocity-Verlet solution:

[[File:Classicalsolution1.jpg]]

The total energy for the oscillator varies as shown:

[[File:Energy111.jpg]]

An approximate linear fit has been performed on the maxima of the error of the calculations; the absolute difference between classical and velocity-Verlet solutions. Iterations of using previous results causes error to propagate and increase.

[[File:Errorfit.jpg]]

The smaller the timestep, the smaller fluctuations in total energy. It is important to monitor the total energy of a physical system to ensure energy is conserved, however infinitesimally small timesteps greatly increase time needed to run simulations. A time step of '''0.028s''' allows energy fluctuations to be as low as 1% and this is a good trade off with simulation run time.

=== Atomic Forces ===

For a single Lennard-Jones interaction, <math>\phi (r) = 4\epsilon (\frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6})</math>.

:* When potential energy is 0, φ(r) = 0 and r=r0
::If <math>\phi=0, (\frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6}) = 0</math>,
:::<math>\frac{\sigma^{12}}{r^{12}} = \frac{\sigma^{6}}{r^{6}}</math>
:::<math>r_0 = \sigma</math>

:*The force is given by <math>F= \frac{d\phi}{dr}=4\epsilon(-12\sigma^{12}r^{-13}+6\sigma^{6}r^{-7})</math> and at a potential energy of 0, <math>F=\frac{-24\epsilon}{\sigma}</math>.

:*At equilibrium separation <math>r_{eq}</math>, <math>\frac{d\phi}{dr}=0</math>.
::<math>0=4\epsilon(\frac{-12\sigma^{12}}{r^{13}}+\frac{6\sigma^6}{r^7})</math>
::<math>(\frac{12\sigma^{12}}{r^{13}} = \frac{6\sigma^6}{r^7})</math>
::<math>r_{eq}=2^{\frac{1}{6}}\sigma</math>

:*At equilibrium separation, the well depth:
::<math>\phi(r_{eq})=\phi(2^{\frac{1}{6}}\sigma)=-\epsilon</math>

:*<math>\int\phi(r)\, dr=4\epsilon[\frac{-\sigma^{12}}{11r^{11}}+\frac{\sigma^6}{5r^5}]+c</math> and given that <math>\sigma = \epsilon = 1.0</math> so <math>\int\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]+c</math>
::*<math>\int\limits_{2\sigma}^{\infty}\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]_2^{\infty} = -4(\frac{1}{5(2)^5}-\frac{1}{11(2)^{11}}) = -0.02482</math>
::*<math>\int\limits_{2.5\sigma}^{\infty}\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]_{2.5}^{\infty} = -4(\frac{1}{5(2.5)^5}-\frac{1}{11(2.5)^{11}}) = -0.008177</math>
::*<math>\int\limits_{3\sigma}^{\infty}\phi(r)\, dr=4\left[\frac{-1}{11r^{11}}+\frac{1}{5r^5}\right]_{3}^{\infty} = -4(\frac{1}{5(3)^5}-\frac{1}{11(3)^{11}}) = -0.00329</math>

=== Periodic Boundary Conditions ===

As <math> pV=Nk_BT </math>, the number of water molecules in 1 mL of water is approximately <math>2.46x10^{19}</math> and 10000 molecules takes up an approximate volume of <math>4.06x10^{-22} m^3</math>.

In a simulation box which runs from (0,0,0) to (1,1,1), an atom that starts at (0.5,0.5,0.5) and moves along vector (0.7,0.6,0.2), will end up at (0.2,0.1,0.7) once periodic boundary conditions have been applied.

=== Reduced Units ===

The LJ parameters for Argon are: <math>\sigma=0.34nm, \frac{\epsilon}{k_B}=120K</math>.
*<math>r=\sigma \cdot r*= 0.34 \cdot 3.2 = 1.088 nm</math>.
*<math> Well Depth = \phi(r_{eq}) = - \epsilon = - k_B \cdot 120 = -1.656^{-21}J</math>
:<math>-1.656^{-21} \cdot N_A = -997.4 = -0.997 kJ mol^{-1}</math>.
*<math>T=\frac{\epsilon \cdot T*}{k_B}=120 \cdot 1.5 = 180 K</math>.

== Equilibriation ==

=== Creating the Simulation Box ===

If two atoms are generated too close together, the LJ potential shows that the potential between the two would be infinitely large, making force simulations between these two atoms too large to realistically simulate. The LJ cutoff also ensures that LJ potentials are only calculated for atoms that are near enough, and not every other atom in the infinitely repeating lattice, which would greatly increase simulation run time.

A lattice spacing of 1.07722 corresponds to a lattice number density of <math>\frac{1}{1.07722^3}=0.8</math> for a simple cubic lattice. A face centred cubic lattice has 4 lattice points per cell, and thus would require a lattice spacing of <math> (\frac{4}{1.2})^{\frac{1}{3}}= 1.4938</math>. A 10x10x10 box would contain 1000 unit cells, and 4000 lattice points, so the create_atoms command for such a lattice would create 4000 atoms.

=== Setting the Properties of the Atoms ===

The command '''mass 1 1.0''' assigns all atoms of type 1 a mass of 1.0.
The command '''pair_style lj/cut 3.0''' defines the cutoff distance between atoms that have a potential between them to be 3.0 (ie. the simulation does not run for atoms farther apart or closer than this distance).
The command '''pair_coeff ** 1.0 1.0''' specifically defines the pairwise force field coefficients for multiple pairs atoms.

The velocity-Verlet algorithm is the numerical integration method that will be used if <math>x_i(0)</math> and <math>v_i(0)</math> are defined.

=== Running the Simulation ===

Calling upon variables, instead of assigning numbers, makes it much easier to change these variables for every simulation that is run.

=== Checking Equilibriation ===

The simulation takes about 0.3 seconds to equilibriate energy, temperature, and pressure, as shown below:

[[File:energyeq.jpg]]
[[File:energyeqzoom.jpg]]
[[File:tempeq.jpg]]
[[File:tempeqzoom.jpg]]
[[File:presseq.jpg]]
[[File:presseqzoom.jpg]]

Of the five timsteps used, 0.0025 is the largest acceptable timestep to use as a smaller timestep of 0.01 results in a very similar equilibriation. 0.015 does not equilbriate at all as the time steps are too large for the numerical integration to accurately find an average for the ensemble, and energy accumulates.

[[File:Equilibriation3.jpg]]

== Running Simulations Under Specific Conditions ==

=== Thermostats & Barostats ===

<math> E_K = \frac{3}{2} Nk_B\mathfrak{T} = \frac{1}{2}\sum_{i} m_i (\gamma v_i)^2</math>

:<math> \gamma^2 = \frac{\frac{3}{2}Nk_B\mathfrak{T}}{\frac{1}{2}\sum_{i}m_i v_i^2}</math>

:<math> \gamma^2 = \frac{\frac{3}{2}Nk_B\mathfrak{T}}{\frac{3}{2}Nk_BT}</math>

:<math> \gamma = (\frac{\mathfrak{T}}{T})^\frac{1}{2} </math>

=== Examining the Input Script ===

The command '''fix aves all ave/time 100 1000 100000''' means that values will be sampled every 100 timesteps; in total 1000 readings will be taken to compute a final average on the 100000th timestep.

'''run 100000''' indicates that 100000 timesteps will be simulated.

=== Plotting the Equations of State ===
[[File:download1.jpg]]

Higher pressures lead to higher densities, both in theory and in these simulations. Our simulated density is higher than that given by the ideal gas law because the simulation takes particle interactions into account. The error increases at higher pressures, when more collisions are likely to occur.

== Calculating Heat Capacities Using Statistical Physics ==
[[File:HEATCAP.jpg]]

Higher pressure results in higher heat capacity as the increased number of molecules per unit volume that can absorb energy to their vibrational excited states. As the simulation is in a lattice, rotational degrees of freedom are not available to the atoms, and so heat capacity decreases as temperature increases, despite expectations.

An example of the input scripts is below:

[[Media:inputnpt.in]]

== Structural Properties and the Radial Distribution Function ==
[[File:RDF1.jpg]]

The RDF shows the probability of finding a particle at a distance r from a reference particle, relative to an ideal gas. In a gas, there is little order and minimal structure to particles and so the graph has minimal features.

Liquids are slightly more ordered and the decreasing heights of peaks of the RDF correlate to coordination spheres. There is a high probability of finding another particle in a primary coordination sphere but this probability decreases as you go farther away from the reference particle.

The solid FCC lattice has a much higher order, and the RDF peak separation and heights define the lattice structure. The first, second, and third sharp peaks refer to different sets of nearest neighbours, while their heights show how many of those nearest neighbours there are. The lattice spacing is the same as the distance to the second nearest neighbour, 1.475.

The following images show a reference particle in red and first, second, and third nearest neighbouring particles respectively. The rest of the lattice is shown in black.
[[File:nearestn.jpg]]
[[File:nearestn1.jpg]]

The coordination numbers are 12 (Int(g(1.205)=12, 12 neighbours), 6 (Int(g(1.475)=18, 6 neighbours), and 24 (Int(g(1.775)=42, 24 neighbours) respectively.

[[File:Int1.jpg]]

== Dynamical Properties and the Diffusion Coefficient ==

=== Mean Squared Displacement ===
*<math> D_{solid} = \frac{1}{6} \cdot 0.00382752458475 = 6.38x10^{-4}</math>
[[File:solid.jpg]]

*<math> D_{liquid} = \frac{1}{6} \cdot 0.509774821123 = 0.085</math>
[[File:liquid1.jpg]]

*<math> D_{gas} = \frac{1}{6} \cdot 7.99193006423 = 1.33</math>
[[File:gas11.jpg]]

The diffusion coefficient increases as entropy of the phase increases, which matches expectations as gas particles are much more likely to diffuse than a rigid lattice of solid molecules.

==== 1000000 Atoms ====

[[File:MSD11.jpg]]

*<math> D_{solid} = \frac{1}{6} \cdot 2.79195534196x10^5 = 4.65x10^{-6} </math>
*<math> D_{liquid} = \frac{1}{6} \cdot 0.531614512766 = 0.0886</math>
*<math> D_{gas} = \frac{1}{6} \cdot 18.0968139669 = 3.02</math>

The MSD graph for a gas is curved at first, indicating ballistic motion proportional to <math>T^{2}</math>. After enough collisions have occurred, diffusion is linear, as it is for a liquid which constantly has the same collisions. The diffusion coefficient is close to 0 for solids which is as expected.

=== Velocity Autocorrelation Function ===

::*<math> x(t)=Acos(\omega t + \phi) </math> and <math> v(t) = \frac{dx(t)}{dt} = -A\omega sin (\omega t+ \phi) </math>

<math> C(\tau) = \frac{\int\limits_{-\infty}^{\infty}v(t)-v(t+\tau)dt}{\int\limits_{-\infty}^{\infty}v^2(t)dt} = \frac{\int\limits_{-\infty}^{\infty}sin(\omega t+\phi)sin(\omega(t+\tau)+\phi)}{\int\limits_{-\infty}^{\infty}sin^2(\omega t + \phi)dt} </math>

::*<math> sin(A+B) = sinAcosB + cosAsinB </math>

<math> C(\tau) = \frac{[cos (\omega \tau) \int\limits_{-\infty}^{\infty} sin^2(\omega t + \phi) dt] + [sin (\omega \tau) \int\limits_{-\infty}^{\infty} sin(\omega t + \phi) cos(\omega t + \phi) dt]}{\int\limits_{-\infty}^{\infty}sin^2(\omega t + \phi)dt} </math>

::*<math> sin^2(x) = \frac{1}{2} (1-cos(2x)) </math>

<math> C(\tau) = \frac{[cos (\omega \tau) \int\limits_{-\infty}^{\infty} \frac{1}{2} (1-cos(2(\omega t + \phi)] + [sin (\omega \tau) \int\limits_{-\infty}^{\infty} sin(\omega t + \phi) cos(\omega t + \phi) dt]}{\int\limits_{-\infty}^{\infty} \frac{1}{2} (1-cos(2(\omega t + \phi) dt} </math>

::*<math> \int \frac{1}{2} (1-cos(2(\omega t + \phi)) dt = \frac{t}{2} + \frac{sin(2(\omega t + \phi)}{4\omega} + c </math>
::*<math> \int sin(\omega t + \phi) cos(\omega t + \phi) dt = \frac{sin^2(\omega t + \phi)}{2\omega} + c</math>

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

::*sin(x) is an odd function and integrating between <math>-\infty</math> and <math>\infty</math> will result in 0

<math> C(\tau) = \frac{\frac{tcos(\omega \tau)}{2}}{\frac{t}{2}} = cos (\omega \tau)</math>

VACF minima refer to collisions of particles where velocity is instantaneously 0, negative as they are in the opposite direction to <math>\tau</math>. As VACF is averaged over all molecules, they cancel out once they are out of phase, which happens faster for liquids than it does for solids. In comparison to the harmonic oscillator, which only models one particle without any collisions, no convergence to 0 occurs.

[[File:VACF1.jpg]]

Further diffusion coefficient estimations, using the trapezium rule, and running integrals are shown below:

*<math> D_{solid} = \frac{1}{3} \cdot 0.43184744815700105 = 0.144 </math>
[[File:runint4.jpg]]
*<math> D_{liquid} = \frac{1}{3} \cdot 146.83331703729999 = 48.9 </math>
[[File:runint5.jpg]]
*<math> D_{gas} = \frac{1}{3} \cdot 1451.848385 =484 </math>
[[File:runint6.jpg]]

==== 1000000 Atoms ====
[[File:VACF1111.jpg]]

*<math> D_{solid} = \frac{1}{3} \cdot -0.41626945400000459 = -0.139</math>
[[File:runint1.jpg]]
*<math> D_{liquid} = \frac{1}{3} \cdot 123.7270701106= 41.2 </math>
[[File:runint2.jpg]]
*<math> D_{gas} = \frac{1}{3} \cdot 1466.443215 = 489 </math>
[[File:runint3.jpg]]

File:VACF1111.jpg

2016-12-14T19:44:54Z

Sfs114: