File:Maxima error vs time OH.png

2016-11-03T22:03:46Z

Sz3614:

File:When timestep=0.12 OH.png

2016-11-03T22:02:58Z

Sz3614:

File:Rdf and integral of rdf summary graph.png

2016-11-03T22:02:07Z

Sz3614:

Rep:SZ3614ls

2016-11-03T21:41:28Z

Sz3614: /* Task */

==Running Simulation==
* Using HPC system to perform input scripts which have different timesteps
* The timesteps are 0.001, 0.0025, 0.0075, 0.01 and 0.015

==Molecular dynamic simulation==
===Task===
*complete column ANALYTICAL, ERROR, and ENERGY
ANALYTICAL is calculated by the classical harmonic oscillator equation. The position of a classical harmonic oscillator is <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

ANALYTICAL vs time graph obtained[[File:Analytical.png]]

ERROR is the absolute value of the difference of ANALYTICAL and x(t). x(t) is already given in the file.

ERROR vs time graph is obtained [[File:Error.png]]

ENERGY is the total energy of the oscillator, <math>E_p + E_k = E_{total}</math>, which is the sum of kinetic energy <math>E_k = \frac{1}{2} mv^2</math> and the potential energy <math>E_p={1\over2}kx^2</math>

Energy vs time graph is obtained[[File:Energy.png]]

===Task===
*find a graph of maxima ERROR vs time

As shown in the graph below the appropriate function is y=0.0004x-0.00007
the graph is obtained [[File:Maxima error vs time.png]]

===Task===
*the timestep to ensure total energy changes not more than 1%

When the timestep is changed to 0.12, the total energy change is not more than 1%. Because the minimum and maximum on the energy graph is still within 1% from the middle value of the energy.

The energy vs time graph when timestep=0.12 is obtained[[File:When timestep=1.2.png]]
===Task===
*Find the separation, the force at this separation, the equilibrium separation, the well depth and the evaluate several integrals.
The separation <math>r_{0}</math> = σ when the potential energy is 0, via L-J potential equation<math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right) = 0 </math>

thus <math>\frac{\sigma^{6}}{r^{6}}=1</math>

thus <math>r_{0}</math> = σ

The force at this separation is 0 because the potential energy is 0 with the equation <math>\mathbf{F}_i = - \frac{\mathrm{d}U\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_i}</math>

The equilibrium separation is when the L-J potential is at the minimum.

Thus differentiation of L-J potential curve at this point should be equal to zero

Thus <math>\frac{\partial \phi \left ( r \right )}{\partial r}</math>=0

Then <math>4\varepsilon \left ( \sigma^{12}r^{-12}- \sigma^{6}r^{-6}\right )</math>=0

Therefore <math>r= \sqrt[6]{2}\sigma</math>

The well depth is the potential well, at the equilibrium, the well depth is -\varepsilon , calculated by substitute <math>r= \sqrt[6]{2}\sigma</math> back into <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>

The integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>=-0.02482, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>=<math>-8.177\times 10^{-3}</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math> =<math>-3.2901\times 10^{-3}</math> when <math>\sigma = \epsilon = 1.0</math>.

===Task===
*no. of water molecule in 1 ml water and volume of 10000 water molecules under under standard conditions.
The number of water molecule is <math>3.346\times 10^{22}</math>, via mass/volume equation and Avogadro's number with the calculation <math>1\div 18\times 6.022\times 10^{23}</math>= <math>3.346\times 10^{22}</math>

The volume is <math>2.99\times 10^{-19}</math> with the calculation <math>10000\div \left ( 6.022\times 10^{23}\right )\times 18</math>=<math>2.99\times 10^{-19}</math>

===Task===
*After boundary condition, the point where it ends up at
It first goes to (1.2, 1.1, 0.7). As the atom goes outside the box, there is another stom goes in the box.

The end up point is (0.2, 0.1, 0.7).

===Task===
*The real unit of distance and temperature. the well depth value.
the real r=1.088 m (<math>r=r^{*}\times \sigma= 3.2\times 0.34= 1.088</math>)and the real T= 180 K (<math>T=T^{*}\times \frac{\epsilon }{k_{B}}= 1.5\times 120= 180 K</math>), via the L-J parameters
Well depth Ԑ in <math>kJmol^{-1} </math>
<math>\epsilon\ /\ k_B= 120 \mathrm{K}</math> with Boltzmann's constant is approximately <math>1.3807\times 10^{-23}</math> joules per kelvin

thus Ԑ = <math>2.75\times 10^{-48}</math> <math>\mathrm{kJ\ mol}^{-1}</math>

==Equilibration==
===Task===
* Why random starting coordinates causes problem.
As two atoms get closer together due to the attraction, they lose potential energy. This energy is liberated as heat. [citation] The computational technique cannot handle the infinite potential therefore the atoms are going to blow up (going far away outside of the simulation box).
===Task===
*for ffc lattice, with lattice point number density 1.2, find the side length of the cubic unit cell
The ffc has four lattice point per one unit cell. Therefore the side length should be <math>\sqrt[3]{4\div 1.2}</math> = 1.494

===Task===
*The number of atoms created for ffc lattice via command
the command create_atom is to create 1 box which is 1000 unit cells. The ffc has four atoms per unit cell. Therefore, 4000 atoms would be created.
===Task===
*find the purpose of some commands
mass 1 1.0: 1 means the atome type is 1, 1.0 means the mass is 1.0

pair_style lj/cut 3.0: lj/cut is the style, 3.0 is the argument

pair_coeff * * 1.0 1.0: ** sets the coefficients for all I J pairs, 1.0 is the argument which is the coefficient for one or more pairs of atom types
===Task===
===Task===
*[[File:The question about the equilibrition.png]]
Use the variable instead just use the number, which allows us to change the variable efficiently. In a script, the variable will occur several times. If we want to change the variable value, we can simple change it at the second line instead of changing it several times in a scrpt.

===Task===
*plots of energy, temperature and pressure against time. equilibrium(how long to reach). plot of energy vs time. find a good timestep.
The plot of energy, temperature, pressure against time for 0.001 timestep is obtained.[[File:0.001 timestep nergy, temperature, and pressure, against time.png]]
The simulation reached equilibrium as the constant energy, temperature, and pressure have obtained with a little bit fluctuation.It takes 0.39s for energy, 0.95s for temperature, 2.24s for pressure.

The plot of energy vs time for all timestep is obtained.[[File:All timestep energy vs time.png]]

The largest timestep to give acceptable results is 0.01 because the simulation still reaches equilibrium even though very scattered graph.

The 0.015 timestep gives bad results. The slope of the energy vs time graph is too large and the simulation does not reach equilibrium

==Simulation with specific conditions==
===Task===
*choose five temperatures and two pressures
pressure: 2.65, 2.80 ( with the average pressure of simulations when the timestep is 0.001 as a basic idea)

Temperature: 1.5, 2.0, 2.5, 3.0, 3.5 (the critical temperature <math>T^* = 1.5</math>)

Timestep: 0.001(the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality)

===Task===
*find <math>\gamma</math> for <math>T = \mathfrak{T}</math>
<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T</math>

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

If we want <math>\gamma</math> for <math>T = \mathfrak{T}</math>, then two equations above should equal to each other.

After rearranging the equations, we get <math>\gamma =\sqrt{\frac{\mathfrak{T}}{T}}</math>

===Task===
*Find the importance of the three numbers 100 1000 100000. how often the values be sampled for average. How many measurements contribute to the average. How much time for simulation.
100 is the Nevery = use input values every this many timesteps
1000 is the Nrepeat = # of times to use input values for calculating averages
10000 is the Nfreq = calculate averages every this many timesteps

Every 100 steptime the values be sampled for average

1000 measurements contribute to the average

<math>100000\times 0.001= 100s</math> for simulation

===Task===
*plot of density vs temperature with error bars. a line corresponding to the density predicted by the ideal law.
The plot is obtained.

The density predicted by the ideal law can be calculated with <math>PV = nRT</math>, and <math>N/V</math> is density so <math>p=density*RT</math>

[[File:Correct density vs t and gas law.png]]

===Task===
*Is your simulated density lower or higher? Does the discrepancy increase or decrease with pressure?
Corresponding to the graph above, the simulated density is lower than the density calculated by ideal gas law.

As the gas law assumption says that there are no (or entirely negligible) intermolecular forces between the gas molecules. This means molecules in gas phases can stick together without considering the intermoleclar forces, which indicates that the density calculated by ideal gas law should be higher than the real density.

According to the graph above, the discrepancy increases with pressure.

==Heat Capacity==
===Task===
*Plot C_V/V as a function of temperature. Attach an input script.
The plot of Cv/V vs Temperature is obtained.[[File:CvV vs T.png]]

*Is the trend the one you would expect?

Yes. As shown in the graph, Cv/V decreases with the increasing temperature. (V remain the same when the density remain the same) As <math>C_V = \frac{\partial E}{\partial T}</math>, the graph is consistent with the theory.

*An input scripts when density=0.2 and temperature=2.0 (see below)

[[File:D=0.2 t=2.0 inputpart1.png]]
[[File:D=0.2 t=2.0 inputpart2.png]]

==Radical distribution Function==
===Task===
*calculate <math>g(r)</math> and <math>\int g(r)\mathrm{d}r</math>. Plot the RDFs for the three systems on the same axes.
The plot of RDFs for three systems is obtained [[File:Rdf for three systems.png]]

The RDF here is a function of the interatomic separation. Three RDFs all showing that at small r, the RDF is zero because of the atoms cannot approach any more closely. Also, significant peaks indicate that atoms pack around each other in 'shells' of neighbours. Finally, as r increases, RDFs tend to be 1 because RDF is the average density at this range.
*Qualitatively the differences between the three RDFs:

In solid, the atoms are positioned at lattice point, therefore, the peaks are very sharp.
As the distance between atoms become further, the liquid is behaving like the gas, as there is no peaks.
*The structure of the system in each phase.
In solid phase, the molecules are tightly bound to each other; in liquid phase, the molecules are more loosely bounded to each other which allows them to move; in gas phase, the molecules are free to move within a large distance.
*In the solid case, illustrate which lattice sites the first three peaks correspond to. What is the lattice spacing? What is the coordination number for each of the first three peaks?
Integral of g(r) obtained[[File:Integral of rdf in solid.png]]

Three plateaus in integral of g(r) graph indicated the number of atoms that can be found at certain distances.

The first three peaks are the three shortest distances, which are 1.325, 1.675 and 2.025, between two lattice point in a fcc unit. The picture below illustrated these three distance.
[[File:Fcc pic.png]]

The lattice spacing is 1.675 indicated by the second data point (1.325, 0.187694215) on the solid RDF graph.

==Dynamic==
===Task===
*a plot showing the "total" MSD vs timestep. Are these as you would expect? Estimate D in each case. Be careful with the units! Repeat this procedure for the MSD data that you were given from the one million atom simulations.
[[File:Msd vs timestep.png]]

The plot is not linear for a short time in the gas graph. This is because the path a molecule takes will only be an approximate straight line until it collides with its neighbour. In gas, the molecule ate far away from each other, therefore this is what we would expect.

*Estimate D in each case
As <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>, to find <math>D</math> is to plot a graph of total MSD vs actual time.

Then we get the slope from the graph which is <math>\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

Therefore we can calculate <math>D</math>. (the values are in the picture below)
The solid phase does not have a diffusion coefficient because the atoms are not diffusing.
[[File:Msd vs time and D value.png]]
*the same procedure with the one million atom simulations
Again, The solid phase does not have a diffusion coefficient because the atoms are not diffusing.
[[File:Msd vs time and D value for one million atom.png]]

===Task===
*evaluate C(τ) , plot C(τ) vs timestep and VACF vs timestep
The answer is <math>C\left ( \tau \right )= cos\left ( \omega \tau \right )</math>

The procedure of the evaluation is below

The position of a 1D harmonic oscillator as a function of time is <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

<math>V\left ( t \right )= \frac{\mathrm{d} x\left ( t \right )}{\mathrm{d} t}</math>

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

As given <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>

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

Fist, we work on the denominator which is <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt }</math>

As we know <math>sin^{2}\left ( \omega t+\phi \right )= \frac{1}{2}-\frac{1}{2}cos\left ( 2\omega t+2\phi \right )</math>

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

Therefore, <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt } = \left [ \frac{1}{2}t-\frac{1}{4} \omega \times sin\left ( 2\omega t+2\phi \right )\right ]_{-\infty}^{\infty}</math>

Secondly, we work on the numerator which is <math>\int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) sin\left ( \omega t+\omega \tau +\phi \right )dt</math>

Expand the sin equation, <math>\int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) sin\left ( \omega t+\omega \tau +\phi \right )dt = \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) \left ( sin\left ( \omega t+\phi \right )cos\left ( \omega \tau \right ) + cos\left ( \omega t+\tau \right )sin\left ( \omega \tau \right ))\right )dt</math>

As <math>cos\left ( \omega \tau \right )</math> and <math>sin\left ( \omega \tau \right )</math> are constant

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

We already know <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt } = \left [ \frac{1}{2}t-\frac{1}{4} \omega \times sin\left ( 2\omega t+2\phi \right )\right ]_{-\infty}^{\infty}</math>

Therefore, we are focusing on <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt</math>this part

Let <math>U= sin\left ( \omega t+\phi \right )</math>

Thus <math>\frac{\mathrm{d} U}{\mathrm{d} t}=\frac{1}{\omega }cos\left ( \omega t+\phi \right )</math>

Substitute <math>cos\left ( \omega t+\phi \right )</math> back in <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt</math>

Therefore we get <math>sin\left ( \omega \tau \right )\int_{-\infty}^{\infty}U\omega dU</math>

Thus <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt = \left [ \frac{1}{2}\omega sin\left ( \omega \tau \right )sin^{2}\left ( \omega t+\phi \right ) \right ]_{-\infty}^{\infty}</math>

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

As t tends to infinity, <math>C\left(\tau\right)=\frac{cos\left ( \omega \tau \right )\times \frac{1}{2}t}{\frac{1}{2}t}</math>

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

Plot obtained[[File:VACF and C(tau).png]]
The minima in the VACFs for the liquid and solid system, very damped oscillations, represent a collision between two atoms.

Even though liquid behave similarly to solids, atoms are close to each other, there are differences between the liquid and solid VACFs.

Because in solid the atoms are fixed at regularly positions, their motion is an oscillation. This is why we have a function that oscillates strongly from positive to negative value and back again.

In liquid, two atoms collide before they rebound from one another and diffuse away.

The harmonic oscillator VACF is cos(ωτ) vs time. Since <math>\omega = 1/2\pi</math>, τ the timestep is the variable here. Since cos function is a periodic function, the shape of harmonic oscillator VACF is very different to the Lennard Jones solid and liquid.

===Task===
*Use the trapezium rule to approximate the integral under the velocity autocorrelation function for the solid, liquid, and gas, and use these values to estimate D in each case. Are they as you expect? What do you think is the largest source of error in your estimates of D from the VACF?

With <math> \int_{a}^{b} f(x)\, dx \approx \frac{h}{2} \sum_{k=1}^{N} \left( f(x_{k+1}) + f(x_{k}) \right)</math> <math>{}= \frac{b-a}{2N}(f(x_1) + 2f(x_2) + 2f(x_3) + 2f(x_4) + \dotsb + 2f(x_N) + f(x_{N+1})).</math> equation, the integral under VACF is 0.803506332 in solid, 119.1372588 in liquid and 1277.155395 in gas.

Therefore, the D in solid phase is 0.267835444, in liquid phase is 39.71241959 and in gas phase is 425.718465.

As <math>D = \frac{1}{3}\int_0^\infty \mathrm{d}\tau \left\langle\mathbf{v}\left(0\right)\cdot\mathbf{v}\left(\tau\right)\right\rangle</math>, we can calculate D by dividing the integral under the velocity autocorrelation function by three.
*A plot of the running integral
[[File:Running integral summary.png]]

According to the graphs of VACF vs timestep, the running integral plots are reasonable.
*Repeat this procedure for the VACF data that you were given from the one million atom simulations.
With same approach as before, the integral under VACF is 0.068294227 in solid, 135.1372144 in liquid and 4902.698697in gas.

Therefore, the D in solid phase is 0.022764742, in liquid phase is 45.04573814 and in gas phase is 1634.232899.

[[File:Running integral summary for one million.png]]

==Conclusion&Summary==
*We run the simulations with different timesteps to compare which timestep is more accurate.
*With the velocity-Verlet algorithm and the Lennard-Jones potential, boundary conditions and reduce unit, we understand the theory for the simulation.
*Analyse the output of the simulation with thermodynamic properties to check the equilibrium.
*Under NpT conditions, plotting the equation of different states.
*Heat capacities are obtained in density-temperature phase space with the calculation of statistical physics.
*Radical distribution functions are obtained to understand how, on average, the atoms in a system are radially packed around each other.
*Mean Squared Displacement and Velocity Autocorrelation Function are obtained to understand how the square distance grows with time (since is squared, it will not be sums up to zero with both positive and negative values )

==References==
==Appendix==
All the graphs obtained

The specific calculation for all tasks

Rep:SZ3614ls

2016-11-03T21:31:34Z

Sz3614: /* Task */

==Running Simulation==
* Using HPC system to perform input scripts which have different timesteps
* The timesteps are 0.001, 0.0025, 0.0075, 0.01 and 0.015

==Molecular dynamic simulation==
===Task===
*complete column ANALYTICAL, ERROR, and ENERGY
ANALYTICAL is calculated by the classical harmonic oscillator equation. The position of a classical harmonic oscillator is <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

ANALYTICAL vs time graph obtained[[File:Analytical.png]]

ERROR is the absolute value of the difference of ANALYTICAL and x(t). x(t) is already given in the file.

ERROR vs time graph is obtained [[File:Error.png]]

ENERGY is the total energy of the oscillator, <math>E_p + E_k = E_{total}</math>, which is the sum of kinetic energy <math>E_k = \frac{1}{2} mv^2</math> and the potential energy <math>E_p={1\over2}kx^2</math>

Energy vs time graph is obtained[[File:Energy.png]]

===Task===
*find a graph of maxima ERROR vs time

As shown in the graph below the appropriate function is y=0.0004x-0.00007
the graph is obtained [[File:Maxima error vs time.png]]

===Task===
*the timestep to ensure total energy changes not more than 1%

When the timestep is changed to 0.12, the total energy change is not more than 1%. Because the minimum and maximum on the energy graph is still within 1% from the middle value of the energy.

The energy vs time graph when timestep=0.12 is obtained[[File:When timestep=1.2.png]]
===Task===
*Find the separation, the force at this separation, the equilibrium separation, the well depth and the evaluate several integrals.
The separation <math>r_{0}</math> = σ when the potential energy is 0, via L-J potential equation<math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right) = 0 </math>

thus <math>\frac{\sigma^{6}}{r^{6}}=1</math>

thus <math>r_{0}</math> = σ

The force at this separation is 0 because the potential energy is 0 with the equation <math>\mathbf{F}_i = - \frac{\mathrm{d}U\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_i}</math>

The equilibrium separation is when the L-J potential is at the minimum.

Thus differentiation of L-J potential curve at this point should be equal to zero

Thus <math>\frac{\partial \phi \left ( r \right )}{\partial r}</math>=0

Then <math>4\varepsilon \left ( \sigma^{12}r^{-12}- \sigma^{6}r^{-6}\right )</math>=0

Therefore <math>r= \sqrt[6]{2}\sigma</math>

The well depth is the potential well, at the equilibrium, the well depth is -\varepsilon , calculated by substitute <math>r= \sqrt[6]{2}\sigma</math> back into <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>

The integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>=-0.02482, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>=<math>-8.177\times 10^{-3}</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math> =<math>-3.2901\times 10^{-3}</math> when <math>\sigma = \epsilon = 1.0</math>.

===Task===
*no. of water molecule in 1 ml water and volume of 10000 water molecules under under standard conditions.
The number of water molecule is <math>3.346\times 10^{22}</math>, via mass/volume equation and Avogadro's number with the calculation <math>1\div 18\times 6.022\times 10^{23}</math>= <math>3.346\times 10^{22}</math>

The volume is <math>2.99\times 10^{-19}</math> with the calculation <math>10000\div \left ( 6.022\times 10^{23}\right )\times 18</math>=<math>2.99\times 10^{-19}</math>

===Task===
*After boundary condition, the point where it ends up at
It first goes to (1.2, 1.1, 0.7). As the atom goes outside the box, there is another stom goes in the box.

The end up point is (0.2, 0.1, 0.7).

===Task===
*The real unit of distance and temperature. the well depth value.
the real r=1.088 m (<math>r=r^{*}\times \sigma= 3.2\times 0.34= 1.088</math>)and the real T= 180 K (<math>T=T^{*}\times \frac{\epsilon }{k_{B}}= 1.5\times 120= 180 K</math>), via the L-J parameters
Well depth Ԑ in <math>kJmol^{-1} </math>
<math>\epsilon\ /\ k_B= 120 \mathrm{K}</math> with Boltzmann's constant is approximately <math>1.3807\times 10^{-23}</math> joules per kelvin

thus Ԑ = <math>2.75\times 10^{-48}</math> <math>\mathrm{kJ\ mol}^{-1}</math>

==Equilibration==
===Task===
* Why random starting coordinates causes problem.
As two atoms get closer together due to the attraction, they lose potential energy. This energy is liberated as heat. [citation] The computational technique cannot handle the infinite potential therefore the atoms are going to blow up (going far away outside of the simulation box).
===Task===
*for ffc lattice, with lattice point number density 1.2, find the side length of the cubic unit cell
The ffc has four lattice point per one unit cell. Therefore the side length should be <math>\sqrt[3]{4\div 1.2}</math> = 1.494

===Task===
*The number of atoms created for ffc lattice via command
the command create_atom is to create 1 box which is 1000 unit cells. The ffc has four atoms per unit cell. Therefore, 4000 atoms would be created.
===Task===
*find the purpose of some commands
mass 1 1.0: 1 means the atome type is 1, 1.0 means the mass is 1.0

pair_style lj/cut 3.0: lj/cut is the style, 3.0 is the argument

pair_coeff * * 1.0 1.0: ** sets the coefficients for all I J pairs, 1.0 is the argument which is the coefficient for one or more pairs of atom types
===Task===
===Task===
*[[File:The question about the equilibrition.png]]
Use the variable instead just use the number, which allows us to change the variable efficiently. In a script, the variable will occur several times. If we want to change the variable value, we can simple change it at the second line instead of changing it several times in a scrpt.

===Task===
*plots of energy, temperature and pressure against time. equilibrium(how long to reach). plot of energy vs time. find a good timestep.
The plot of energy, temperature, pressure against time for 0.001 timestep is obtained.[[File:0.001 timestep nergy, temperature, and pressure, against time.png]]
The simulation reached equilibrium as the constant energy, temperature, and pressure have obtained with a little bit fluctuation.It takes 0.39s for energy, 0.95s for temperature, 2.24s for pressure.

The plot of energy vs time for all timestep is obtained.[[File:All timestep energy vs time.png]]

The largest timestep to give acceptable results is 0.01 because the simulation still reaches equilibrium even though very scattered graph.

The 0.015 timestep gives bad results. The slope of the energy vs time graph is too large and the simulation does not reach equilibrium

==Simulation with specific conditions==
===Task===
*choose five temperatures and two pressures
pressure: 2.65, 2.80 ( with the average pressure of simulations when the timestep is 0.001 as a basic idea)

Temperature: 1.5, 2.0, 2.5, 3.0, 3.5 (the critical temperature <math>T^* = 1.5</math>)

Timestep: 0.001(the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality)

===Task===
*find <math>\gamma</math> for <math>T = \mathfrak{T}</math>
<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T</math>

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

If we want <math>\gamma</math> for <math>T = \mathfrak{T}</math>, then two equations above should equal to each other.

After rearranging the equations, we get <math>\gamma =\sqrt{\frac{\mathfrak{T}}{T}}</math>

===Task===
*Find the importance of the three numbers 100 1000 100000. how often the values be sampled for average. How many measurements contribute to the average. How much time for simulation.
100 is the Nevery = use input values every this many timesteps
1000 is the Nrepeat = # of times to use input values for calculating averages
10000 is the Nfreq = calculate averages every this many timesteps

Every 100 steptime the values be sampled for average

1000 measurements contribute to the average

<math>100000\times 0.001= 100s</math> for simulation

===Task===
*plot of density vs temperature with error bars. a line corresponding to the density predicted by the ideal law.
The plot is obtained.

The density predicted by the ideal law can be calculated with <math>PV = nRT</math>, and <math>N/V</math> is density so <math>p=density*RT</math>

[[File:Correct density vs t and gas law.png]]

===Task===
*Is your simulated density lower or higher? Does the discrepancy increase or decrease with pressure?
Corresponding to the graph above, the simulated density is lower than the density calculated by ideal gas law.

As the gas law assumption says that there are no (or entirely negligible) intermolecular forces between the gas molecules. This means molecules in gas phases can stick together without considering the intermoleclar forces, which indicates that the density calculated by ideal gas law should be higher than the real density.

According to the graph above, the discrepancy increases with pressure.

==Heat Capacity==
===Task===
*Plot C_V/V as a function of temperature. Attach an input script.
The plot of Cv/V vs Temperature is obtained.[[File:CvV vs T.png]]

*Is the trend the one you would expect?

Yes. As shown in the graph, Cv/V decreases with the increasing temperature. (V remain the same when the density remain the same) As <math>C_V = \frac{\partial E}{\partial T}</math>, the graph is consistent with the theory.

*An input scripts when density=0.2 and temperature=2.0 (see below)

[[File:D=0.2 t=2.0 inputpart1.png]]
[[File:D=0.2 t=2.0 inputpart2.png]]

==Radical distribution Function==
===Task===
*calculate <math>g(r)</math> and <math>\int g(r)\mathrm{d}r</math>. Plot the RDFs for the three systems on the same axes.
The plot of RDFs for three systems is obtained [[File:Rdf for three systems.png]]

The RDF here is a function of the interatomic separation. Three RDFs all showing that at small r, the RDF is zero because of the atoms cannot approach any more closely. Also, significant peaks indicate that atoms pack around each other in 'shells' of neighbours. Finally, as r increases, RDFs tend to be 1 because RDF is the average density at this range.
*Qualitatively the differences between the three RDFs:

In solid, the atoms are positioned at lattice point, therefore, the peaks are very sharp.
As the distance between atoms become further, the liquid is behaving like the gas, as there is no peaks.
*The structure of the system in each phase.
In solid phase, the molecules are tightly bound to each other; in liquid phase, the molecules are more loosely bounded to each other which allows them to move; in gas phase, the molecules are free to move within a large distance.
*In the solid case, illustrate which lattice sites the first three peaks correspond to. What is the lattice spacing? What is the coordination number for each of the first three peaks?
Integral of g(r) obtained[[File:Integral of rdf in solid.png]]

Three plateaus in integral of g(r) graph indicated the number of atoms that can be found at certain distances.

The first three peaks are the three shortest distances, which are 1.325, 1.675 and 2.025, between two lattice point in a fcc unit. The picture below illustrated these three distance.
[[File:Fcc pic.png]]

The lattice spacing is 1.675 indicated by the second data point (1.325, 0.187694215) on the solid RDF graph.

==Dynamic==
===Task===
*a plot showing the "total" MSD vs timestep. Are these as you would expect? Estimate D in each case. Be careful with the units! Repeat this procedure for the MSD data that you were given from the one million atom simulations.
[[File:Msd vs timestep.png]]

The plot is not linear for a short time in the gas graph. This is because the path a molecule takes will only be an approximate straight line until it collides with its neighbour. In gas, the molecule ate far away from each other, therefore this is what we would expect.

*Estimate D in each case
As <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>, to find <math>D</math> is to plot a graph of total MSD vs actual time.

Then we get the slope from the graph which is <math>\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

Therefore we can calculate <math>D</math>. (the values are in the picture below)
The solid phase does not have a diffusion coefficient because the atoms are not diffusing.
[[File:Msd vs time and D value.png]]
*the same procedure with the one million atom simulations
Again, The solid phase does not have a diffusion coefficient because the atoms are not diffusing.
[[File:Msd vs time and D value for one million atom.png]]

===Task===
*evaluate C(τ) , plot C(τ) vs timestep and VACF vs timestep
The answer is <math>C\left ( \tau \right )= cos\left ( \omega \tau \right )</math>

The procedure of the evaluation is below

The position of a 1D harmonic oscillator as a function of time is <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

<math>V\left ( t \right )= \frac{\mathrm{d} x\left ( t \right )}{\mathrm{d} t}</math>

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

As given <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>

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

Fist, we work on the denominator which is <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt }</math>

As we know <math>sin^{2}\left ( \omega t+\phi \right )= \frac{1}{2}-\frac{1}{2}cos\left ( 2\omega t+2\phi \right )</math>

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

Therefore, <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt } = \left [ \frac{1}{2}t-\frac{1}{4} \omega \times sin\left ( 2\omega t+2\phi \right )\right ]_{-\infty}^{\infty}</math>

Secondly, we work on the numerator which is <math>\int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) sin\left ( \omega t+\omega \tau +\phi \right )dt</math>

Expand the sin equation, <math>\int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) sin\left ( \omega t+\omega \tau +\phi \right )dt = \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) \left ( sin\left ( \omega t+\phi \right )cos\left ( \omega \tau \right ) + cos\left ( \omega t+\tau \right )sin\left ( \omega \tau \right ))\right )dt</math>

As <math>cos\left ( \omega \tau \right )</math> and <math>sin\left ( \omega \tau \right )</math> are constant

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

We already know <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt } = \left [ \frac{1}{2}t-\frac{1}{4} \omega \times sin\left ( 2\omega t+2\phi \right )\right ]_{-\infty}^{\infty}</math>

Therefore, we are focusing on <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt</math>this part

Let <math>U= sin\left ( \omega t+\phi \right )</math>

Thus <math>\frac{\mathrm{d} U}{\mathrm{d} t}=\frac{1}{\omega }cos\left ( \omega t+\phi \right )</math>

Substitute <math>cos\left ( \omega t+\phi \right )</math> back in <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt</math>

Therefore we get <math>sin\left ( \omega \tau \right )\int_{-\infty}^{\infty}U\omega dU</math>

Thus <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt = \left [ \frac{1}{2}\omega sin\left ( \omega \tau \right )sin^{2}\left ( \omega t+\phi \right ) \right ]_{-\infty}^{\infty}</math>

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

As t tends to infinity, <math>C\left(\tau\right)=\frac{cos\left ( \omega \tau \right )\times \frac{1}{2}t}{\frac{1}{2}t}</math>

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

Plot obtained[[File:VACF and C(tau).png]]
The minima in the VACFs for the liquid and solid system, very damped oscillations, represent a collision between two atoms.

Even though liquid behave similarly to solids, atoms are close to each other, there are differences between the liquid and solid VACFs.

Because in solid the atoms are fixed at regularly positions, their motion is an oscillation. This is why we have a function that oscillates strongly from positive to negative value and back again.

In liquid, two atoms collide before they rebound from one another and diffuse away.

The harmonic oscillator VACF is cos(ωτ) vs time. Since <math>\omega = 1/2\pi</math>, τ the timestep is the variable here. Since cos function is a periodic function, the shape of harmonic oscillator VACF is very different to the Lennard Jones solid and liquid.

===Task===
*Use the trapezium rule to approximate the integral under the velocity autocorrelation function for the solid, liquid, and gas, and use these values to estimate D in each case. Are they as you expect? What do you think is the largest source of error in your estimates of D from the VACF?

With <math> \int_{a}^{b} f(x)\, dx \approx \frac{h}{2} \sum_{k=1}^{N} \left( f(x_{k+1}) + f(x_{k}) \right)</math> <math>{}= \frac{b-a}{2N}(f(x_1) + 2f(x_2) + 2f(x_3) + 2f(x_4) + \dotsb + 2f(x_N) + f(x_{N+1})).</math> equation, the integral under VACF is 0.803506332 in solid, 119.1372588 in liquid and 1277.155395 in gas.

Therefore, the D in solid phase is 0.267835444, in liquid phase is 39.71241959 and in gas phase is 425.718465.

As <math>D = \frac{1}{3}\int_0^\infty \mathrm{d}\tau \left\langle\mathbf{v}\left(0\right)\cdot\mathbf{v}\left(\tau\right)\right\rangle</math>, we can calculate D by dividing the integral under the velocity autocorrelation function by three.
*A plot of the running integral
[[File:Running integral summary.png]]

According to the graphs of VACF vs timestep, the running integral plots are reasonable.
*Repeat this procedure for the VACF data that you were given from the one million atom simulations.
With same approach as before, the integral under VACF is in solid, in liquid and in gas.

Therefore, the D in solid phase is , in liquid phase is and in gas phase is .

[[File:Running integral summary for one million.png]]

==Conclusion&Summary==
*We run the simulations with different timesteps to compare which timestep is more accurate.
*With the velocity-Verlet algorithm and the Lennard-Jones potential, boundary conditions and reduce unit, we understand the theory for the simulation.
*Analyse the output of the simulation with thermodynamic properties to check the equilibrium.
*Under NpT conditions, plotting the equation of different states.
*Heat capacities are obtained in density-temperature phase space with the calculation of statistical physics.
*Radical distribution functions are obtained to understand how, on average, the atoms in a system are radially packed around each other.
*Mean Squared Displacement and Velocity Autocorrelation Function are obtained to understand how the square distance grows with time (since is squared, it will not be sums up to zero with both positive and negative values )

==References==
==Appendix==
All the graphs obtained

The specific calculation for all tasks

Rep:SZ3614ls

2016-11-03T21:27:22Z

Sz3614: /* Task */

==Running Simulation==
* Using HPC system to perform input scripts which have different timesteps
* The timesteps are 0.001, 0.0025, 0.0075, 0.01 and 0.015

==Molecular dynamic simulation==
===Task===
*complete column ANALYTICAL, ERROR, and ENERGY
ANALYTICAL is calculated by the classical harmonic oscillator equation. The position of a classical harmonic oscillator is <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

ANALYTICAL vs time graph obtained[[File:Analytical.png]]

ERROR is the absolute value of the difference of ANALYTICAL and x(t). x(t) is already given in the file.

ERROR vs time graph is obtained [[File:Error.png]]

ENERGY is the total energy of the oscillator, <math>E_p + E_k = E_{total}</math>, which is the sum of kinetic energy <math>E_k = \frac{1}{2} mv^2</math> and the potential energy <math>E_p={1\over2}kx^2</math>

Energy vs time graph is obtained[[File:Energy.png]]

===Task===
*find a graph of maxima ERROR vs time

As shown in the graph below the appropriate function is y=0.0004x-0.00007
the graph is obtained [[File:Maxima error vs time.png]]

===Task===
*the timestep to ensure total energy changes not more than 1%

When the timestep is changed to 0.12, the total energy change is not more than 1%. Because the minimum and maximum on the energy graph is still within 1% from the middle value of the energy.

The energy vs time graph when timestep=0.12 is obtained[[File:When timestep=1.2.png]]
===Task===
*Find the separation, the force at this separation, the equilibrium separation, the well depth and the evaluate several integrals.
The separation <math>r_{0}</math> = σ when the potential energy is 0, via L-J potential equation<math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right) = 0 </math>

thus <math>\frac{\sigma^{6}}{r^{6}}=1</math>

thus <math>r_{0}</math> = σ

The force at this separation is 0 because the potential energy is 0 with the equation <math>\mathbf{F}_i = - \frac{\mathrm{d}U\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_i}</math>

The equilibrium separation is when the L-J potential is at the minimum.

Thus differentiation of L-J potential curve at this point should be equal to zero

Thus <math>\frac{\partial \phi \left ( r \right )}{\partial r}</math>=0

Then <math>4\varepsilon \left ( \sigma^{12}r^{-12}- \sigma^{6}r^{-6}\right )</math>=0

Therefore <math>r= \sqrt[6]{2}\sigma</math>

The well depth is the potential well, at the equilibrium, the well depth is -\varepsilon , calculated by substitute <math>r= \sqrt[6]{2}\sigma</math> back into <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>

The integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>=-0.02482, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>=<math>-8.177\times 10^{-3}</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math> =<math>-3.2901\times 10^{-3}</math> when <math>\sigma = \epsilon = 1.0</math>.

===Task===
*no. of water molecule in 1 ml water and volume of 10000 water molecules under under standard conditions.
The number of water molecule is <math>3.346\times 10^{22}</math>, via mass/volume equation and Avogadro's number with the calculation <math>1\div 18\times 6.022\times 10^{23}</math>= <math>3.346\times 10^{22}</math>

The volume is <math>2.99\times 10^{-19}</math> with the calculation <math>10000\div \left ( 6.022\times 10^{23}\right )\times 18</math>=<math>2.99\times 10^{-19}</math>

===Task===
*After boundary condition, the point where it ends up at
It first goes to (1.2, 1.1, 0.7). As the atom goes outside the box, there is another stom goes in the box.

The end up point is (0.2, 0.1, 0.7).

===Task===
*The real unit of distance and temperature. the well depth value.
the real r=1.088 m (<math>r=r^{*}\times \sigma= 3.2\times 0.34= 1.088</math>)and the real T= 180 K (<math>T=T^{*}\times \frac{\epsilon }{k_{B}}= 1.5\times 120= 180 K</math>), via the L-J parameters
Well depth Ԑ in <math>kJmol^{-1} </math>
<math>\epsilon\ /\ k_B= 120 \mathrm{K}</math> with Boltzmann's constant is approximately <math>1.3807\times 10^{-23}</math> joules per kelvin

thus Ԑ = <math>2.75\times 10^{-48}</math> <math>\mathrm{kJ\ mol}^{-1}</math>

==Equilibration==
===Task===
* Why random starting coordinates causes problem.
As two atoms get closer together due to the attraction, they lose potential energy. This energy is liberated as heat. [citation] The computational technique cannot handle the infinite potential therefore the atoms are going to blow up (going far away outside of the simulation box).
===Task===
*for ffc lattice, with lattice point number density 1.2, find the side length of the cubic unit cell
The ffc has four lattice point per one unit cell. Therefore the side length should be <math>\sqrt[3]{4\div 1.2}</math> = 1.494

===Task===
*The number of atoms created for ffc lattice via command
the command create_atom is to create 1 box which is 1000 unit cells. The ffc has four atoms per unit cell. Therefore, 4000 atoms would be created.
===Task===
*find the purpose of some commands
mass 1 1.0: 1 means the atome type is 1, 1.0 means the mass is 1.0

pair_style lj/cut 3.0: lj/cut is the style, 3.0 is the argument

pair_coeff * * 1.0 1.0: ** sets the coefficients for all I J pairs, 1.0 is the argument which is the coefficient for one or more pairs of atom types
===Task===
===Task===
*[[File:The question about the equilibrition.png]]
Use the variable instead just use the number, which allows us to change the variable efficiently. In a script, the variable will occur several times. If we want to change the variable value, we can simple change it at the second line instead of changing it several times in a scrpt.

===Task===
*plots of energy, temperature and pressure against time. equilibrium(how long to reach). plot of energy vs time. find a good timestep.
The plot of energy, temperature, pressure against time for 0.001 timestep is obtained.[[File:0.001 timestep nergy, temperature, and pressure, against time.png]]
The simulation reached equilibrium as the constant energy, temperature, and pressure have obtained with a little bit fluctuation.It takes 0.39s for energy, 0.95s for temperature, 2.24s for pressure.

The plot of energy vs time for all timestep is obtained.[[File:All timestep energy vs time.png]]

The largest timestep to give acceptable results is 0.01 because the simulation still reaches equilibrium even though very scattered graph.

The 0.015 timestep gives bad results. The slope of the energy vs time graph is too large and the simulation does not reach equilibrium

==Simulation with specific conditions==
===Task===
*choose five temperatures and two pressures
pressure: 2.65, 2.80 ( with the average pressure of simulations when the timestep is 0.001 as a basic idea)

Temperature: 1.5, 2.0, 2.5, 3.0, 3.5 (the critical temperature <math>T^* = 1.5</math>)

Timestep: 0.001(the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality)

===Task===
*find <math>\gamma</math> for <math>T = \mathfrak{T}</math>
<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T</math>

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

If we want <math>\gamma</math> for <math>T = \mathfrak{T}</math>, then two equations above should equal to each other.

After rearranging the equations, we get <math>\gamma =\sqrt{\frac{\mathfrak{T}}{T}}</math>

===Task===
*Find the importance of the three numbers 100 1000 100000. how often the values be sampled for average. How many measurements contribute to the average. How much time for simulation.
100 is the Nevery = use input values every this many timesteps
1000 is the Nrepeat = # of times to use input values for calculating averages
10000 is the Nfreq = calculate averages every this many timesteps

Every 100 steptime the values be sampled for average

1000 measurements contribute to the average

<math>100000\times 0.001= 100s</math> for simulation

===Task===
*plot of density vs temperature with error bars. a line corresponding to the density predicted by the ideal law.
The plot is obtained.

The density predicted by the ideal law can be calculated with <math>PV = nRT</math>, and <math>N/V</math> is density so <math>p=density*RT</math>

[[File:Correct density vs t and gas law.png]]

===Task===
*Is your simulated density lower or higher? Does the discrepancy increase or decrease with pressure?
Corresponding to the graph above, the simulated density is lower than the density calculated by ideal gas law.

As the gas law assumption says that there are no (or entirely negligible) intermolecular forces between the gas molecules. This means molecules in gas phases can stick together without considering the intermoleclar forces, which indicates that the density calculated by ideal gas law should be higher than the real density.

According to the graph above, the discrepancy increases with pressure.

==Heat Capacity==
===Task===
*Plot C_V/V as a function of temperature. Attach an input script.
The plot of Cv/V vs Temperature is obtained.[[File:CvV vs T.png]]

*Is the trend the one you would expect?

Yes. As shown in the graph, Cv/V decreases with the increasing temperature. (V remain the same when the density remain the same) As <math>C_V = \frac{\partial E}{\partial T}</math>, the graph is consistent with the theory.

*An input scripts when density=0.2 and temperature=2.0 (see below)

[[File:D=0.2 t=2.0 inputpart1.png]]
[[File:D=0.2 t=2.0 inputpart2.png]]

==Radical distribution Function==
===Task===
*calculate <math>g(r)</math> and <math>\int g(r)\mathrm{d}r</math>. Plot the RDFs for the three systems on the same axes.
The plot of RDFs for three systems is obtained [[File:Rdf for three systems.png]]

The RDF here is a function of the interatomic separation. Three RDFs all showing that at small r, the RDF is zero because of the atoms cannot approach any more closely. Also, significant peaks indicate that atoms pack around each other in 'shells' of neighbours. Finally, as r increases, RDFs tend to be 1 because RDF is the average density at this range.
*Qualitatively the differences between the three RDFs:

In solid, the atoms are positioned at lattice point, therefore, the peaks are very sharp.
As the distance between atoms become further, the liquid is behaving like the gas, as there is no peaks.
*The structure of the system in each phase.
In solid phase, the molecules are tightly bound to each other; in liquid phase, the molecules are more loosely bounded to each other which allows them to move; in gas phase, the molecules are free to move within a large distance.
*In the solid case, illustrate which lattice sites the first three peaks correspond to. What is the lattice spacing? What is the coordination number for each of the first three peaks?
Integral of g(r) obtained[[File:Integral of rdf in solid.png]]

Three plateaus in integral of g(r) graph indicated the number of atoms that can be found at certain distances.

The first three peaks are the three shortest distances, which are 1.325, 1.675 and 2.025, between two lattice point in a fcc unit. The picture below illustrated these three distance.
[[File:Fcc pic.png]]

The lattice spacing is 1.675 indicated by the second data point (1.325, 0.187694215) on the solid RDF graph.

==Dynamic==
===Task===
*a plot showing the "total" MSD vs timestep. Are these as you would expect? Estimate D in each case. Be careful with the units! Repeat this procedure for the MSD data that you were given from the one million atom simulations.
[[File:Msd vs timestep.png]]

The plot is not linear for a short time in the gas graph. This is because the path a molecule takes will only be an approximate straight line until it collides with its neighbour. In gas, the molecule ate far away from each other, therefore this is what we would expect.

*Estimate D in each case
As <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>, to find <math>D</math> is to plot a graph of total MSD vs actual time.

Then we get the slope from the graph which is <math>\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

Therefore we can calculate <math>D</math>. (the values are in the picture below)
The solid phase does not have a diffusion coefficient because the atoms are not diffusing.
[[File:Msd vs time and D value.png]]
*the same procedure with the one million atom simulations
Again, The solid phase does not have a diffusion coefficient because the atoms are not diffusing.
[[File:Msd vs time and D value for one million atom.png]]

===Task===
*evaluate C(τ) , plot C(τ) vs timestep and VACF vs timestep
The answer is <math>C\left ( \tau \right )= cos\left ( \omega \tau \right )</math>

The procedure of the evaluation is below

The position of a 1D harmonic oscillator as a function of time is <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

<math>V\left ( t \right )= \frac{\mathrm{d} x\left ( t \right )}{\mathrm{d} t}</math>

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

As given <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>

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

Fist, we work on the denominator which is <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt }</math>

As we know <math>sin^{2}\left ( \omega t+\phi \right )= \frac{1}{2}-\frac{1}{2}cos\left ( 2\omega t+2\phi \right )</math>

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

Therefore, <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt } = \left [ \frac{1}{2}t-\frac{1}{4} \omega \times sin\left ( 2\omega t+2\phi \right )\right ]_{-\infty}^{\infty}</math>

Secondly, we work on the numerator which is <math>\int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) sin\left ( \omega t+\omega \tau +\phi \right )dt</math>

Expand the sin equation, <math>\int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) sin\left ( \omega t+\omega \tau +\phi \right )dt = \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) \left ( sin\left ( \omega t+\phi \right )cos\left ( \omega \tau \right ) + cos\left ( \omega t+\tau \right )sin\left ( \omega \tau \right ))\right )dt</math>

As <math>cos\left ( \omega \tau \right )</math> and <math>sin\left ( \omega \tau \right )</math> are constant

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

We already know <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt } = \left [ \frac{1}{2}t-\frac{1}{4} \omega \times sin\left ( 2\omega t+2\phi \right )\right ]_{-\infty}^{\infty}</math>

Therefore, we are focusing on <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt</math>this part

Let <math>U= sin\left ( \omega t+\phi \right )</math>

Thus <math>\frac{\mathrm{d} U}{\mathrm{d} t}=\frac{1}{\omega }cos\left ( \omega t+\phi \right )</math>

Substitute <math>cos\left ( \omega t+\phi \right )</math> back in <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt</math>

Therefore we get <math>sin\left ( \omega \tau \right )\int_{-\infty}^{\infty}U\omega dU</math>

Thus <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt = \left [ \frac{1}{2}\omega sin\left ( \omega \tau \right )sin^{2}\left ( \omega t+\phi \right ) \right ]_{-\infty}^{\infty}</math>

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

As t tends to infinity, <math>C\left(\tau\right)=\frac{cos\left ( \omega \tau \right )\times \frac{1}{2}t}{\frac{1}{2}t}</math>

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

Plot obtained[[File:VACF and C(tau).png]]
The minima in the VACFs for the liquid and solid system, very damped oscillations, represent a collision between two atoms.

Even though liquid behave similarly to solids, atoms are close to each other, there are differences between the liquid and solid VACFs.

Because in solid the atoms are fixed at regularly positions, their motion is an oscillation. This is why we have a function that oscillates strongly from positive to negative value and back again.

In liquid, two atoms collide before they rebound from one another and diffuse away.

The harmonic oscillator VACF is cos(ωτ) vs time. Since <math>\omega = 1/2\pi</math>, τ the timestep is the variable here. Since cos function is a periodic function, the shape of harmonic oscillator VACF is very different to the Lennard Jones solid and liquid.

===Task===
*Use the trapezium rule to approximate the integral under the velocity autocorrelation function for the solid, liquid, and gas, and use these values to estimate D in each case. Are they as you expect? What do you think is the largest source of error in your estimates of D from the VACF?

With <math> \int_{a}^{b} f(x)\, dx \approx \frac{h}{2} \sum_{k=1}^{N} \left( f(x_{k+1}) + f(x_{k}) \right)</math> <math>{}= \frac{b-a}{2N}(f(x_1) + 2f(x_2) + 2f(x_3) + 2f(x_4) + \dotsb + 2f(x_N) + f(x_{N+1})).</math> equation, the integral under VACF is 0.803506332 in solid, 119.1372588 in liquid and 1277.155395 in gas.

Therefore, the D in solid phase is 0.267835444, in liquid phase is 39.71241959 and in gas phase is 425.718465.

As <math>D = \frac{1}{3}\int_0^\infty \mathrm{d}\tau \left\langle\mathbf{v}\left(0\right)\cdot\mathbf{v}\left(\tau\right)\right\rangle</math>, we can calculate D by dividing the integral under the velocity autocorrelation function by three.
*A plot of the running integral
[[File:Running integral summary.png]]

According to the graphs of VACF vs timestep, the running integral plots are reasonable.
*Repeat this procedure for the VACF data that you were given from the one million atom simulations.
[[File:Running integral summary for one million.png]]

==Conclusion&Summary==
*We run the simulations with different timesteps to compare which timestep is more accurate.
*With the velocity-Verlet algorithm and the Lennard-Jones potential, boundary conditions and reduce unit, we understand the theory for the simulation.
*Analyse the output of the simulation with thermodynamic properties to check the equilibrium.
*Under NpT conditions, plotting the equation of different states.
*Heat capacities are obtained in density-temperature phase space with the calculation of statistical physics.
*Radical distribution functions are obtained to understand how, on average, the atoms in a system are radially packed around each other.
*Mean Squared Displacement and Velocity Autocorrelation Function are obtained to understand how the square distance grows with time (since is squared, it will not be sums up to zero with both positive and negative values )

==References==
==Appendix==
All the graphs obtained

The specific calculation for all tasks

2016-11-03T20:45:34Z

Sz3614: /* Task */

==Running Simulation==
* Using HPC system to perform input scripts which have different timesteps
* The timesteps are 0.001, 0.0025, 0.0075, 0.01 and 0.015

==Molecular dynamic simulation==
===Task===
*complete column ANALYTICAL, ERROR, and ENERGY
ANALYTICAL is calculated by the classical harmonic oscillator equation. The position of a classical harmonic oscillator is <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

ANALYTICAL vs time graph obtained[[File:Analytical.png]]

ERROR is the absolute value of the difference of ANALYTICAL and x(t). x(t) is already given in the file.

ERROR vs time graph is obtained [[File:Error.png]]

ENERGY is the total energy of the oscillator, <math>E_p + E_k = E_{total}</math>, which is the sum of kinetic energy <math>E_k = \frac{1}{2} mv^2</math> and the potential energy <math>E_p={1\over2}kx^2</math>

Energy vs time graph is obtained[[File:Energy.png]]

===Task===
*find a graph of maxima ERROR vs time

As shown in the graph below the appropriate function is y=0.0004x-0.00007
the graph is obtained [[File:Maxima error vs time.png]]

===Task===
*the timestep to ensure total energy changes not more than 1%

When the timestep is changed to 0.12, the total energy change is not more than 1%. Because the minimum and maximum on the energy graph is still within 1% from the middle value of the energy.

The energy vs time graph when timestep=0.12 is obtained[[File:When timestep=1.2.png]]
===Task===
*Find the separation, the force at this separation, the equilibrium separation, the well depth and the evaluate several integrals.
The separation <math>r_{0}</math> = σ when the potential energy is 0, via L-J potential equation<math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right) = 0 </math>

thus <math>\frac{\sigma^{6}}{r^{6}}=1</math>

thus <math>r_{0}</math> = σ

The force at this separation is 0 because the potential energy is 0 with the equation <math>\mathbf{F}_i = - \frac{\mathrm{d}U\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_i}</math>

The equilibrium separation is when the L-J potential is at the minimum.

Thus differentiation of L-J potential curve at this point should be equal to zero

Thus <math>\frac{\partial \phi \left ( r \right )}{\partial r}</math>=0

Then <math>4\varepsilon \left ( \sigma^{12}r^{-12}- \sigma^{6}r^{-6}\right )</math>=0

Therefore <math>r= \sqrt[6]{2}\sigma</math>

The well depth is the potential well, at the equilibrium, the well depth is -\varepsilon , calculated by substitute <math>r= \sqrt[6]{2}\sigma</math> back into <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>

The integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>=-0.02482, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>=<math>-8.177\times 10^{-3}</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math> =<math>-3.2901\times 10^{-3}</math> when <math>\sigma = \epsilon = 1.0</math>.

===Task===
*no. of water molecule in 1 ml water and volume of 10000 water molecules under under standard conditions.
The number of water molecule is <math>3.346\times 10^{22}</math>, via mass/volume equation and Avogadro's number with the calculation <math>1\div 18\times 6.022\times 10^{23}</math>= <math>3.346\times 10^{22}</math>

The volume is <math>2.99\times 10^{-19}</math> with the calculation <math>10000\div \left ( 6.022\times 10^{23}\right )\times 18</math>=<math>2.99\times 10^{-19}</math>

===Task===
*After boundary condition, the point where it ends up at
It first goes to (1.2, 1.1, 0.7). As the atom goes outside the box, there is another stom goes in the box.

The end up point is (0.2, 0.1, 0.7).

===Task===
*The real unit of distance and temperature. the well depth value.
the real r=1.088 m (<math>r=r^{*}\times \sigma= 3.2\times 0.34= 1.088</math>)and the real T= 180 K (<math>T=T^{*}\times \frac{\epsilon }{k_{B}}= 1.5\times 120= 180 K</math>), via the L-J parameters
Well depth Ԑ in <math>kJmol^{-1} </math>
<math>\epsilon\ /\ k_B= 120 \mathrm{K}</math> with Boltzmann's constant is approximately <math>1.3807\times 10^{-23}</math> joules per kelvin

thus Ԑ = <math>2.75\times 10^{-48}</math> <math>\mathrm{kJ\ mol}^{-1}</math>

==Equilibration==
===Task===
* Why random starting coordinates causes problem.
As two atoms get closer together due to the attraction, they lose potential energy. This energy is liberated as heat. [citation] The computational technique cannot handle the infinite potential therefore the atoms are going to blow up (going far away outside of the simulation box).
===Task===
*for ffc lattice, with lattice point number density 1.2, find the side length of the cubic unit cell
The ffc has four lattice point per one unit cell. Therefore the side length should be <math>\sqrt[3]{4\div 1.2}</math> = 1.494

===Task===
*The number of atoms created for ffc lattice via command
the command create_atom is to create 1 box which is 1000 unit cells. The ffc has four atoms per unit cell. Therefore, 4000 atoms would be created.
===Task===
*find the purpose of some commands
mass 1 1.0: 1 means the atome type is 1, 1.0 means the mass is 1.0

pair_style lj/cut 3.0: lj/cut is the style, 3.0 is the argument

pair_coeff * * 1.0 1.0: ** sets the coefficients for all I J pairs, 1.0 is the argument which is the coefficient for one or more pairs of atom types
===Task===
===Task===
*[[File:The question about the equilibrition.png]]
Use the variable instead just use the number, which allows us to change the variable efficiently. In a script, the variable will occur several times. If we want to change the variable value, we can simple change it at the second line instead of changing it several times in a scrpt.

===Task===
*plots of energy, temperature and pressure against time. equilibrium(how long to reach). plot of energy vs time. find a good timestep.
The plot of energy, temperature, pressure against time for 0.001 timestep is obtained.[[File:0.001 timestep nergy, temperature, and pressure, against time.png]]
The simulation reached equilibrium as the constant energy, temperature, and pressure have obtained with a little bit fluctuation.It takes 0.39s for energy, 0.95s for temperature, 2.24s for pressure.

The plot of energy vs time for all timestep is obtained.[[File:All timestep energy vs time.png]]

The largest timestep to give acceptable results is 0.01 because the simulation still reaches equilibrium even though very scattered graph.

The 0.015 timestep gives bad results. The slope of the energy vs time graph is too large and the simulation does not reach equilibrium

==Simulation with specific conditions==
===Task===
*choose five temperatures and two pressures
pressure: 2.65, 2.80 ( with the average pressure of simulations when the timestep is 0.001 as a basic idea)

Temperature: 1.5, 2.0, 2.5, 3.0, 3.5 (the critical temperature <math>T^* = 1.5</math>)

Timestep: 0.001(the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality)

===Task===
*find <math>\gamma</math> for <math>T = \mathfrak{T}</math>
<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T</math>

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

If we want <math>\gamma</math> for <math>T = \mathfrak{T}</math>, then two equations above should equal to each other.

After rearranging the equations, we get <math>\gamma =\sqrt{\frac{\mathfrak{T}}{T}}</math>

===Task===
*Find the importance of the three numbers 100 1000 100000. how often the values be sampled for average. How many measurements contribute to the average. How much time for simulation.
100 is the Nevery = use input values every this many timesteps
1000 is the Nrepeat = # of times to use input values for calculating averages
10000 is the Nfreq = calculate averages every this many timesteps

Every 100 steptime the values be sampled for average

1000 measurements contribute to the average

<math>100000\times 0.001= 100s</math> for simulation

===Task===
*plot of density vs temperature with error bars. a line corresponding to the density predicted by the ideal law.
The plot is obtained.

The density predicted by the ideal law can be calculated with <math>PV = nRT</math>, and <math>N/V</math> is density so <math>p=density*RT</math>

[[File:Correct density vs t and gas law.png]]

===Task===
*Is your simulated density lower or higher? Does the discrepancy increase or decrease with pressure?
Corresponding to the graph above, the simulated density is lower than the density calculated by ideal gas law.

As the gas law assumption says that there are no (or entirely negligible) intermolecular forces between the gas molecules. This means molecules in gas phases can stick together without considering the intermoleclar forces, which indicates that the density calculated by ideal gas law should be higher than the real density.

According to the graph above, the discrepancy increases with pressure.

==Heat Capacity==
===Task===
*Plot C_V/V as a function of temperature. Attach an input script.
The plot of Cv/V vs Temperature is obtained.[[File:CvV vs T.png]]

*Is the trend the one you would expect?

Yes. As shown in the graph, Cv/V decreases with the increasing temperature. (V remain the same when the density remain the same) As <math>C_V = \frac{\partial E}{\partial T}</math>, the graph is consistent with the theory.

*An input scripts when density=0.2 and temperature=2.0 (see below)

[[File:D=0.2 t=2.0 inputpart1.png]]
[[File:D=0.2 t=2.0 inputpart2.png]]

==Radical distribution Function==
===Task===
*calculate <math>g(r)</math> and <math>\int g(r)\mathrm{d}r</math>. Plot the RDFs for the three systems on the same axes.
The plot of RDFs for three systems is obtained [[File:Rdf for three systems.png]]

The RDF here is a function of the interatomic separation. Three RDFs all showing that at small r, the RDF is zero because of the atoms cannot approach any more closely. Also, significant peaks indicate that atoms pack around each other in 'shells' of neighbours. Finally, as r increases, RDFs tend to be 1 because RDF is the average density at this range.
*Qualitatively the differences between the three RDFs:

In solid, the atoms are positioned at lattice point, therefore, the peaks are very sharp.
As the distance between atoms become further, the liquid is behaving like the gas, as there is no peaks.
*The structure of the system in each phase.
In solid phase, the molecules are tightly bound to each other; in liquid phase, the molecules are more loosely bounded to each other which allows them to move; in gas phase, the molecules are free to move within a large distance.
*In the solid case, illustrate which lattice sites the first three peaks correspond to. What is the lattice spacing? What is the coordination number for each of the first three peaks?
Integral of g(r) obtained[[File:Integral of rdf in solid.png]]

Three plateaus in integral of g(r) graph indicated the number of atoms that can be found at certain distances.

The first three peaks are the three shortest distances between two lattice point in a fcc unit.

The lattice spacing is 1.325 indicated by the data point (1.325, 0.187694215) on the solid RDF graph.

==Dynamic==
===Task===
*a plot showing the "total" MSD vs timestep. Are these as you would expect? Estimate D in each case. Be careful with the units! Repeat this procedure for the MSD data that you were given from the one million atom simulations.
[[File:Msd vs timestep.png]]

The plot is not linear for a short time in the gas graph. This is because the path a molecule takes will only be an approximate straight line until it collides with its neighbour. In gas, the molecule ate far away from each other, therefore this is what we would expect.

*Estimate D in each case
As <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>, to find <math>D</math> is to plot a graph of total MSD vs actual time.

Then we get the slope from the graph which is <math>\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

Therefore we can calculate <math>D</math>. (the values are in the picture below)
The solid phase does not have a diffusion coefficient because the atoms are not diffusing.
[[File:Msd vs time and D value.png]]
*the same procedure with the one million atom simulations
Again, The solid phase does not have a diffusion coefficient because the atoms are not diffusing.
[[File:Msd vs time and D value for one million atom.png]]

===Task===
*evaluate C(τ) , plot C(τ) vs timestep and VACF vs timestep
The answer is <math>C\left ( \tau \right )= cos\left ( \omega \tau \right )</math>

The procedure of the evaluation is below

The position of a 1D harmonic oscillator as a function of time is <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

<math>V\left ( t \right )= \frac{\mathrm{d} x\left ( t \right )}{\mathrm{d} t}</math>

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

As given <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>

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

Fist, we work on the denominator which is <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt }</math>

As we know <math>sin^{2}\left ( \omega t+\phi \right )= \frac{1}{2}-\frac{1}{2}cos\left ( 2\omega t+2\phi \right )</math>

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

Therefore, <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt } = \left [ \frac{1}{2}t-\frac{1}{4} \omega \times sin\left ( 2\omega t+2\phi \right )\right ]_{-\infty}^{\infty}</math>

Secondly, we work on the numerator which is <math>\int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) sin\left ( \omega t+\omega \tau +\phi \right )dt</math>

Expand the sin equation, <math>\int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) sin\left ( \omega t+\omega \tau +\phi \right )dt = \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) \left ( sin\left ( \omega t+\phi \right )cos\left ( \omega \tau \right ) + cos\left ( \omega t+\tau \right )sin\left ( \omega \tau \right ))\right )dt</math>

As <math>cos\left ( \omega \tau \right )</math> and <math>sin\left ( \omega \tau \right )</math> are constant

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

We already know <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt } = \left [ \frac{1}{2}t-\frac{1}{4} \omega \times sin\left ( 2\omega t+2\phi \right )\right ]_{-\infty}^{\infty}</math>

Therefore, we are focusing on <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt</math>this part

Let <math>U= sin\left ( \omega t+\phi \right )</math>

Thus <math>\frac{\mathrm{d} U}{\mathrm{d} t}=\frac{1}{\omega }cos\left ( \omega t+\phi \right )</math>

Substitute <math>cos\left ( \omega t+\phi \right )</math> back in <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt</math>

Therefore we get <math>sin\left ( \omega \tau \right )\int_{-\infty}^{\infty}U\omega dU</math>

Thus <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt = \left [ \frac{1}{2}\omega sin\left ( \omega \tau \right )sin^{2}\left ( \omega t+\phi \right ) \right ]_{-\infty}^{\infty}</math>

C(τ) =cos(ωτ)

Plot obtained[[File:VACF and C(tau).png]]
The minima in the VACFs for the liquid and solid system, very damped oscillations, represent a collision between two atoms.

Even though liquid behave similarly to solids, atoms are close to each other, there are differences between the liquid and solid VACFs.

Because in solid the atoms are fixed at regularly positions, their motion is an oscillation. This is why we have a function that oscillates strongly from positive to negative value and back again.

In liquid, two atoms collide before they rebound from one another and diffuse away.

The harmonic oscillator VACF is cos(ωτ) vs time. Since <math>\omega = 1/2\pi</math>, τ the timestep is the variable here. Since cos function is a periodic function, the shape of harmonic oscillator VACF is very different to the Lennard Jones solid and liquid.

===Task===
*Use the trapezium rule to approximate the integral under the velocity autocorrelation function for the solid, liquid, and gas, and use these values to estimate D in each case. Are they as you expect? What do you think is the largest source of error in your estimates of D from the VACF?

With <math> \int_{a}^{b} f(x)\, dx \approx \frac{h}{2} \sum_{k=1}^{N} \left( f(x_{k+1}) + f(x_{k}) \right)</math> <math>{}= \frac{b-a}{2N}(f(x_1) + 2f(x_2) + 2f(x_3) + 2f(x_4) + \dotsb + 2f(x_N) + f(x_{N+1})).</math> equation, the integral under VACF is 0.803506332 in solid, 119.1372588 in liquid and 1277.155395 in gas.

Therefore, the D in solid phase is 0.267835444, in liquid phase is 39.71241959 and in gas phase is 425.718465.

As <math>D = \frac{1}{3}\int_0^\infty \mathrm{d}\tau \left\langle\mathbf{v}\left(0\right)\cdot\mathbf{v}\left(\tau\right)\right\rangle</math>, we can calculate D by dividing the integral under the velocity autocorrelation function by three.
*A plot of the running integral
[[File:Running integral summary.png]]

According to the graphs of VACF vs timestep, the running integral plots are reasonable.
*Repeat this procedure for the VACF data that you were given from the one million atom simulations.
[[File:Running integral summary for one million.png]]

==Conclusion&Summary==
*We run the simulations with different timesteps to compare which timestep is more accurate.
*With the velocity-Verlet algorithm and the Lennard-Jones potential, boundary conditions and reduce unit, we understand the theory for the simulation.
*Analyse the output of the simulation with thermodynamic properties to check the equilibrium.
*Under NpT conditions, plotting the equation of different states.
*Heat capacities are obtained in density-temperature phase space with the calculation of statistical physics.
*Radical distribution functions are obtained to understand how, on average, the atoms in a system are radially packed around each other.
*Mean Squared Displacement and Velocity Autocorrelation Function are obtained to understand how the square distance grows with time (since is squared, it will not be sums up to zero with both positive and negative values )

==References==
==Appendix==
The specific calculation for all tasks

Rep:SZ3614ls

2016-11-03T20:43:38Z

Sz3614: /* Task */

==Running Simulation==
* Using HPC system to perform input scripts which have different timesteps
* The timesteps are 0.001, 0.0025, 0.0075, 0.01 and 0.015

==Molecular dynamic simulation==
===Task===
*complete column ANALYTICAL, ERROR, and ENERGY
ANALYTICAL is calculated by the classical harmonic oscillator equation. The position of a classical harmonic oscillator is <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

ANALYTICAL vs time graph obtained[[File:Analytical.png]]

ERROR is the absolute value of the difference of ANALYTICAL and x(t). x(t) is already given in the file.

ERROR vs time graph is obtained [[File:Error.png]]

ENERGY is the total energy of the oscillator, <math>E_p + E_k = E_{total}</math>, which is the sum of kinetic energy <math>E_k = \frac{1}{2} mv^2</math> and the potential energy <math>E_p={1\over2}kx^2</math>

Energy vs time graph is obtained[[File:Energy.png]]

===Task===
*find a graph of maxima ERROR vs time

As shown in the graph below the appropriate function is y=0.0004x-0.00007
the graph is obtained [[File:Maxima error vs time.png]]

===Task===
*the timestep to ensure total energy changes not more than 1%

When the timestep is changed to 0.12, the total energy change is not more than 1%. Because the minimum and maximum on the energy graph is still within 1% from the middle value of the energy.

The energy vs time graph when timestep=0.12 is obtained[[File:When timestep=1.2.png]]
===Task===
*Find the separation, the force at this separation, the equilibrium separation, the well depth and the evaluate several integrals.
The separation <math>r_{0}</math> = σ when the potential energy is 0, via L-J potential equation<math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right) = 0 </math>

thus <math>\frac{\sigma^{6}}{r^{6}}=1</math>

thus <math>r_{0}</math> = σ

The force at this separation is 0 because the potential energy is 0 with the equation <math>\mathbf{F}_i = - \frac{\mathrm{d}U\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_i}</math>

The equilibrium separation is when the L-J potential is at the minimum.

Thus differentiation of L-J potential curve at this point should be equal to zero

Thus <math>\frac{\partial \phi \left ( r \right )}{\partial r}</math>=0

Then <math>4\varepsilon \left ( \sigma^{12}r^{-12}- \sigma^{6}r^{-6}\right )</math>=0

Therefore <math>r= \sqrt[6]{2}\sigma</math>

The well depth is the potential well, at the equilibrium, the well depth is -\varepsilon , calculated by substitute <math>r= \sqrt[6]{2}\sigma</math> back into <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>

The integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>=-0.02482, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>=<math>-8.177\times 10^{-3}</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math> =<math>-3.2901\times 10^{-3}</math> when <math>\sigma = \epsilon = 1.0</math>.

===Task===
*no. of water molecule in 1 ml water and volume of 10000 water molecules under under standard conditions.
The number of water molecule is <math>3.346\times 10^{22}</math>, via mass/volume equation and Avogadro's number with the calculation <math>1\div 18\times 6.022\times 10^{23}</math>= <math>3.346\times 10^{22}</math>

The volume is <math>2.99\times 10^{-19}</math> with the calculation <math>10000\div \left ( 6.022\times 10^{23}\right )\times 18</math>=<math>2.99\times 10^{-19}</math>

===Task===
*After boundary condition, the point where it ends up at
It first goes to (1.2, 1.1, 0.7). As the atom goes outside the box, there is another stom goes in the box.

The end up point is (0.2, 0.1, 0.7).

===Task===
*The real unit of distance and temperature. the well depth value.
the real r=1.088 m (<math>r=r^{*}\times \sigma= 3.2\times 0.34= 1.088</math>)and the real T= 180 K (<math>T=T^{*}\times \frac{\epsilon }{k_{B}}= 1.5\times 120= 180 K</math>), via the L-J parameters
Well depth Ԑ in <math>kJmol^{-1} </math>
<math>\epsilon\ /\ k_B= 120 \mathrm{K}</math> with Boltzmann's constant is approximately <math>1.3807\times 10^{-23}</math> joules per kelvin

thus Ԑ = <math>2.75\times 10^{-48}</math> <math>\mathrm{kJ\ mol}^{-1}</math>

==Equilibration==
===Task===
* Why random starting coordinates causes problem.
As two atoms get closer together due to the attraction, they lose potential energy. This energy is liberated as heat. [citation] The computational technique cannot handle the infinite potential therefore the atoms are going to blow up (going far away outside of the simulation box).
===Task===
*for ffc lattice, with lattice point number density 1.2, find the side length of the cubic unit cell
The ffc has four lattice point per one unit cell. Therefore the side length should be <math>\sqrt[3]{4\div 1.2}</math> = 1.494

===Task===
*The number of atoms created for ffc lattice via command
the command create_atom is to create 1 box which is 1000 unit cells. The ffc has four atoms per unit cell. Therefore, 4000 atoms would be created.
===Task===
*find the purpose of some commands
mass 1 1.0: 1 means the atome type is 1, 1.0 means the mass is 1.0

pair_style lj/cut 3.0: lj/cut is the style, 3.0 is the argument

pair_coeff * * 1.0 1.0: ** sets the coefficients for all I J pairs, 1.0 is the argument which is the coefficient for one or more pairs of atom types
===Task===
===Task===
*[[File:The question about the equilibrition.png]]
Use the variable instead just use the number, which allows us to change the variable efficiently. In a script, the variable will occur several times. If we want to change the variable value, we can simple change it at the second line instead of changing it several times in a scrpt.

===Task===
*plots of energy, temperature and pressure against time. equilibrium(how long to reach). plot of energy vs time. find a good timestep.
The plot of energy, temperature, pressure against time for 0.001 timestep is obtained.[[File:0.001 timestep nergy, temperature, and pressure, against time.png]]
The simulation reached equilibrium as the constant energy, temperature, and pressure have obtained with a little bit fluctuation.It takes 0.39s for energy, 0.95s for temperature, 2.24s for pressure.

The plot of energy vs time for all timestep is obtained.[[File:All timestep energy vs time.png]]

The largest timestep to give acceptable results is 0.01 because the simulation still reaches equilibrium even though very scattered graph.

The 0.015 timestep gives bad results. The slope of the energy vs time graph is too large and the simulation does not reach equilibrium

==Simulation with specific conditions==
===Task===
*choose five temperatures and two pressures
pressure: 2.65, 2.80 ( with the average pressure of simulations when the timestep is 0.001 as a basic idea)

Temperature: 1.5, 2.0, 2.5, 3.0, 3.5 (the critical temperature <math>T^* = 1.5</math>)

Timestep: 0.001(the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality)

===Task===
*find <math>\gamma</math> for <math>T = \mathfrak{T}</math>
<math>\frac{1}{2}\sum_i m_i v_i^2 = \frac{3}{2} N k_B T</math>

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

If we want <math>\gamma</math> for <math>T = \mathfrak{T}</math>, then two equations above should equal to each other.

After rearranging the equations, we get <math>\gamma =\sqrt{\frac{\mathfrak{T}}{T}}</math>

===Task===
*Find the importance of the three numbers 100 1000 100000. how often the values be sampled for average. How many measurements contribute to the average. How much time for simulation.
100 is the Nevery = use input values every this many timesteps
1000 is the Nrepeat = # of times to use input values for calculating averages
10000 is the Nfreq = calculate averages every this many timesteps

Every 100 steptime the values be sampled for average

1000 measurements contribute to the average

<math>100000\times 0.001= 100s</math> for simulation

===Task===
*plot of density vs temperature with error bars. a line corresponding to the density predicted by the ideal law.
The plot is obtained.

The density predicted by the ideal law can be calculated with <math>PV = nRT</math>, and <math>N/V</math> is density so <math>p=density*RT</math>

[[File:Correct density vs t and gas law.png]]

===Task===
*Is your simulated density lower or higher? Does the discrepancy increase or decrease with pressure?
Corresponding to the graph above, the simulated density is lower than the density calculated by ideal gas law.

As the gas law assumption says that there are no (or entirely negligible) intermolecular forces between the gas molecules. This means molecules in gas phases can stick together without considering the intermoleclar forces, which indicates that the density calculated by ideal gas law should be higher than the real density.

According to the graph above, the discrepancy increases with pressure.

==Heat Capacity==
===Task===
*Plot C_V/V as a function of temperature. Attach an input script.
The plot of Cv/V vs Temperature is obtained.[[File:CvV vs T.png]]

*Is the trend the one you would expect?

Yes. As shown in the graph, Cv/V decreases with the increasing temperature. (V remain the same when the density remain the same) As <math>C_V = \frac{\partial E}{\partial T}</math>, the graph is consistent with the theory.

*An input scripts when density=0.2 and temperature=2.0 (see below)

[[File:D=0.2 t=2.0 inputpart1.png]]
[[File:D=0.2 t=2.0 inputpart2.png]]

==Radical distribution Function==
===Task===
*calculate <math>g(r)</math> and <math>\int g(r)\mathrm{d}r</math>. Plot the RDFs for the three systems on the same axes.
The plot of RDFs for three systems is obtained [[File:Rdf for three systems.png]]

The RDF here is a function of the interatomic separation. Three RDFs all showing that at small r, the RDF is zero because of the atoms cannot approach any more closely. Also, significant peaks indicate that atoms pack around each other in 'shells' of neighbours. Finally, as r increases, RDFs tend to be 1 because RDF is the average density at this range.
*Qualitatively the differences between the three RDFs:

In solid, the atoms are positioned at lattice point, therefore, the peaks are very sharp.
As the distance between atoms become further, the liquid is behaving like the gas, as there is no peaks.
*The structure of the system in each phase.
In solid phase, the molecules are tightly bound to each other; in liquid phase, the molecules are more loosely bounded to each other which allows them to move; in gas phase, the molecules are free to move within a large distance.
*In the solid case, illustrate which lattice sites the first three peaks correspond to. What is the lattice spacing? What is the coordination number for each of the first three peaks?
Integral of g(r) obtained[[File:Integral of rdf in solid.png]]

Three plateaus in integral of g(r) graph indicated the number of atoms that can be found at certain distances.

The first three peaks are the three shortest distances between two lattice point in a fcc unit.

The lattice spacing is 1.325 indicated by the data point (1.325, 0.187694215) on the solid RDF graph.

==Dynamic==
===Task===
*a plot showing the "total" MSD vs timestep. Are these as you would expect? Estimate D in each case. Be careful with the units! Repeat this procedure for the MSD data that you were given from the one million atom simulations.
[[File:Msd vs timestep.png]]

The plot is not linear for a short time in the gas graph. This is because the path a molecule takes will only be an approximate straight line until it collides with its neighbour. In gas, the molecule ate far away from each other, therefore this is what we would expect.

*Estimate D in each case
As <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>, to find <math>D</math> is to plot a graph of total MSD vs actual time.

Then we get the slope from the graph which is <math>\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

Therefore we can calculate <math>D</math>. (the values are in the picture below)
The solid phase does not have a diffusion coefficient because the atoms are not diffusing.
[[File:Msd vs time and D value.png]]
*the same procedure with the one million atom simulations
Again, The solid phase does not have a diffusion coefficient because the atoms are not diffusing.
[[File:Msd vs time and D value for one million atom.png]]

===Task===
*evaluate C(τ) , plot C(τ) vs timestep and VACF vs timestep
The position of a 1D harmonic oscillator as a function of time is <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math>

<math>V\left ( t \right )= \frac{\mathrm{d} x\left ( t \right )}{\mathrm{d} t}</math>

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

As given <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>

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

Fist, we work on the denominator which is <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt }</math>

As we know <math>sin^{2}\left ( \omega t+\phi \right )= \frac{1}{2}-\frac{1}{2}cos\left ( 2\omega t+2\phi \right )</math>

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

Therefore, <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt } = \left [ \frac{1}{2}t-\frac{1}{4} \omega \times sin\left ( 2\omega t+2\phi \right )\right ]_{-\infty}^{\infty}</math>

Secondly, we work on the numerator which is <math>\int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) sin\left ( \omega t+\omega \tau +\phi \right )dt</math>

Expand the sin equation, <math>\int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) sin\left ( \omega t+\omega \tau +\phi \right )dt = \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right ) \left ( sin\left ( \omega t+\phi \right )cos\left ( \omega \tau \right ) + cos\left ( \omega t+\tau \right )sin\left ( \omega \tau \right ))\right )dt</math>

As <math>cos\left ( \omega \tau \right )</math> and <math>sin\left ( \omega \tau \right )</math> are constant

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

We already know <math>{\int_{-\infty}^{\infty}sin^{2}\left ( \omega t+\phi \right )dt } = \left [ \frac{1}{2}t-\frac{1}{4} \omega \times sin\left ( 2\omega t+2\phi \right )\right ]_{-\infty}^{\infty}</math>

Therefore, we are focusing on <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt</math>this part

Let <math>U= sin\left ( \omega t+\phi \right )</math>

Thus <math>\frac{\mathrm{d} U}{\mathrm{d} t}=\frac{1}{\omega }cos\left ( \omega t+\phi \right )</math>

Substitute <math>cos\left ( \omega t+\phi \right )</math> back in <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt</math>

Therefore we get <math>sin\left ( \omega \tau \right )\int_{-\infty}^{\infty}U\omega dU</math>

Thus <math>sin\left ( \omega \tau \right ) \int_{-\infty}^{\infty} sin\left ( \omega t+\phi \right )cos\left ( \omega t+\phi \right )dt = \left [ \frac{1}{2}\omega sin\left ( \omega \tau \right )sin^{2}\left ( \omega t+\phi \right ) \right ]_{-\infty}^{\infty}</math>

C(τ) =cos(ωτ)

Plot obtained[[File:VACF and C(tau).png]]
The minima in the VACFs for the liquid and solid system, very damped oscillations, represent a collision between two atoms.

Even though liquid behave similarly to solids, atoms are close to each other, there are differences between the liquid and solid VACFs.

Because in solid the atoms are fixed at regularly positions, their motion is an oscillation. This is why we have a function that oscillates strongly from positive to negative value and back again.

In liquid, two atoms collide before they rebound from one another and diffuse away.

The harmonic oscillator VACF is cos(ωτ) vs time. Since <math>\omega = 1/2\pi</math>, τ the timestep is the variable here. Since cos function is a periodic function, the shape of harmonic oscillator VACF is very different to the Lennard Jones solid and liquid.

===Task===
*Use the trapezium rule to approximate the integral under the velocity autocorrelation function for the solid, liquid, and gas, and use these values to estimate D in each case. Are they as you expect? What do you think is the largest source of error in your estimates of D from the VACF?

With <math> \int_{a}^{b} f(x)\, dx \approx \frac{h}{2} \sum_{k=1}^{N} \left( f(x_{k+1}) + f(x_{k}) \right)</math> <math>{}= \frac{b-a}{2N}(f(x_1) + 2f(x_2) + 2f(x_3) + 2f(x_4) + \dotsb + 2f(x_N) + f(x_{N+1})).</math> equation, the integral under VACF is 0.803506332 in solid, 119.1372588 in liquid and 1277.155395 in gas.

Therefore, the D in solid phase is 0.267835444, in liquid phase is 39.71241959 and in gas phase is 425.718465.

As <math>D = \frac{1}{3}\int_0^\infty \mathrm{d}\tau \left\langle\mathbf{v}\left(0\right)\cdot\mathbf{v}\left(\tau\right)\right\rangle</math>, we can calculate D by dividing the integral under the velocity autocorrelation function by three.
*A plot of the running integral
[[File:Running integral summary.png]]

According to the graphs of VACF vs timestep, the running integral plots are reasonable.
*Repeat this procedure for the VACF data that you were given from the one million atom simulations.
[[File:Running integral summary for one million.png]]

==Conclusion&Summary==
*We run the simulations with different timesteps to compare which timestep is more accurate.
*With the velocity-Verlet algorithm and the Lennard-Jones potential, boundary conditions and reduce unit, we understand the theory for the simulation.
*Analyse the output of the simulation with thermodynamic properties to check the equilibrium.
*Under NpT conditions, plotting the equation of different states.
*Heat capacities are obtained in density-temperature phase space with the calculation of statistical physics.
*Radical distribution functions are obtained to understand how, on average, the atoms in a system are radially packed around each other.
*Mean Squared Displacement and Velocity Autocorrelation Function are obtained to understand how the square distance grows with time (since is squared, it will not be sums up to zero with both positive and negative values )

==References==
==Appendix==
The specific calculation for all tasks