View source for Talk:Modxp412ls

<span style="color:red">A very nice and well attempted report. Remember to be explicit about perhaps even minor details that could be important - why are there fluctuations in the g(r) for solid? Sometimes your explanation was very good but you forgot to make your conclusion - you'll use the 0.0025 timestep. Also be careful to check through for spelling and grammar. All in all though, really good effort! Mark: 80/100</span>


==Theory==
=== Velocity Verlet Algorithm ===
Newton’s second law, <math> F=ma</math>, allows for the calculation of acceleration of each atom in the system from known values of the force exerted on it and  atomic mass. Integration of the equations gives a trajectory that describes the position, velocity and acceleration of the atom with time. The velocity-Verlet algorithm is one of the ways that Newton’s second law may be integrated; it assumes the positions, velocities and accelerations can be approximated by applying a Taylor series expansion and outputs atomic positions, accelerations and velocities at the same instant of time t with very good precision. 

The position of an atom i, at time t is denoted by <math>x_i(t)</math> and the velocity of the atom at time t is <math>v_i(t)</math>. After time step <math>t + \delta t</math>, the position of the atom can be expressed by:

<center><math>\mathbf{x}_i\left(t + \delta t\right) = \mathbf{x}_i\left(t\right) + \frac{\mathrm{d}\mathbf{x}_i\left(t\right)}{\mathrm{d}t}\delta t + \frac{1}{2!}\frac{\mathrm{d}^2\mathbf{x}_i\left(t\right)}{\mathrm{d}t^2}\delta t^2 + \frac{1}{3!}\frac{\mathrm{d}^3\mathbf{x}_i\left(t\right)}{\mathrm{d}t^3} \delta t^3 + \mathcal{O}\left(\delta t^4\right) \ \ (1)</math></center>

A single timestep can be expressed as:

<math> x_(n+1) = x_n +v_n\Delta t+\frac{1}{2}a_n(\Delta t)^2  (2)</math>

<math>v_(n+1)=v_n+\frac{1}{2}(a_(n+1)+a_n)\Delta t  (3)</math> <span style="color:red">Very nice, be careful with your equation formatting (parentheses) </span>

[[File:Energyvvccx.png|550x550px|thumb|center|Fig 1:Velocity Verlet calcuated positions compared to classically calculated positions]]

Classically, the motion of a harmonic oscillator is given by <math> x(t) = A\cos(\omega t + \phi).</math> It can be seen from the position calculation comparisons(Fig 1) that the position values obtained through classical calculations and those obtained using the Velocity Verlet algorithm are in close agreement. This is because the algorithm is self-starting therefore minimising round-off errors. The maximum error of the Velocity Verlet algorithm in this simulation is <math>5.91\times10^-3</math>. Notably, the errors change follows a near oscillatory form has 5 maximas in the time simulated.

[[File:Errormaxfitx.png|551x551px|thumb|center|Fig 2:Error against time and error fit to a quadratic]]

In equation (1), the error of the Verlet algorithm is represented as <math>\Omicron(\delta t^4)</math>. The error per iteration of the algorithm can therefore be represented by:
 
<math>x(t_0+\delta t))= \Omicron(\delta t^4)</math>

After 2 timesteps of <math>\delta t</math>, the Taylor expansion becomes:

<math>x(t_0+2\delta t)= 2x(t_0+\delta t)-x(t_0)+\delta t^2x''(t_0+\delta t)+\Omicron(\delta t^4)</math>

In terms of error:

<math>error(x(t_0+2\delta t))=2error(x(t_0+\delta t)) + \Omicron(\delta t^4)=3\Omicron(\delta t^4)</math>

Through induction, the error can be generalised to:

<math> error(x(t_0 +n\delta t))= \frac{n(n+1)}{2}\Omicron(\delta t^4)</math> 

Considering the error in position between <math>x(t)</math> and <math>x(t+T)</math>, where <math>T=n\delta t</math>:

<math> error(x(t_0 +T))=(\frac{T^2}{2\delta t^2} + \frac{t}{2\delta t}) \mathcal{O}(\delta t^4) </math>

Therefore, the cumulative error over a constant interval is given by:

<math>error(x(t_0+T))=\Omicron(\Delta t^2)</math>

Therefore, the maxima of the error of the Velocity-Verlet algorithm increases quadratically with <math>\delta t</math> and the graph of the maxima of error against time can be fit to a quadratic equation as shown in Fig 2. <span style="color:red">IThis is really good - well done! </span>

The total energy of the oscillator can be calculated by: <math>E=\frac{p^2}{2m} + \frac{kx^2}{2}</math>, where <math>p</math> is the momentum. Since<math>m=1, k=1</math>, the equation can be simplified to <math>E=\frac{p^2+x^2}{2}</math>. To calculate the total energy in terms of velocity,<math>p=mv</math> can be substituted to give <math>E=\frac{(mv)^2+x^2}{2}</math>. The plot of the total energy of the simulated system 


{| class="wikitable" style="text-align: center
|[[File:Energy0.01xx.png|500x500px|thumb|left|Fig 3a: Energy vs Time at 0.01 timestep with error limits of 0.5% either side]]
||[[File:Energyerrorx.png|500x500px|thumb|center|Fig 3b: Energy vs Time at 0.02 timestep with error limits of 0.5% either side]]
|}

In order to ensure that the total energy does not change by more than 1% over the course of the "simulation", the timestep need so to be around 0.2. It is important to monitor the total energy when modelling a physical system numerically to ensure that energy has been conserved, as expected in real systems. <span style="color:red">Tick, correct. </span>

=== Atomic Forces ===
The [https://en.wikipedia.org/wiki/Lennard-Jones_potential Lennard Jones potential](LJ) approximates the potential energy of the interaction between a pair of uncharged atoms. It is mostly expressed in the 12/6 form:

<center><math>  \phi(r)=4\varepsilon\Bigg[\frac{\sigma^{12}}{r^{12}}\Bigg] - \Bigg[\frac{\sigma^6}{r^6}\Bigg]</math></center>

In this equation, <math>\varepsilon</math> is the potential well depth, <math>\sigma</math> is the distance at which the potential between the pair of particles is zero and <math>r</math> is the distance between the pair of particles. 

<b>For a single Lennard-Jones interaction, the potential energy is at zero when <math>r_0=\sigma</math>.</b> Since force is the negative derivative of potential energy, the relationship between force and the LJ potential can be expressed as:

<center><math>F_i=-\frac{d\phi(r^N)}{dr_i}</math></center>

Assuming <math>\varepsilon</math> and <math>\sigma</math> are constant, the derivative of the Lennerad-jones can be expressed as:

<center><math>\frac{d\phi (r^N)}{dr_i}=24\varepsilon\Bigg[-2\Bigg(\frac{\sigma^{12}}{r_i^{13}}\Bigg) + \Bigg(\frac{\sigma^6}{r_i^7}\Bigg)\Bigg]</math></center>

Therefore, when the potential energy is zero, <math>r_0</math> is substituted for <math>\sigma</math> in the equation.

<center><math>\frac{d\phi (r)}{dr_0}=24\varepsilon\Bigg[-2\Bigg(\frac{r_0^{12}}{r_0{13}}\Bigg) + \Bigg(\frac{r_0^6}{r_0^7}\Bigg)\Bigg]</math></center>


<center><math>\frac{d\phi (r)}{dr_0}=24\varepsilon\Bigg[\frac{1}{r_0}-\frac{2}{r_0}\Bigg]</math></center>


<center><math>\frac{d\phi (r)}{dr_0}=24\varepsilon\Bigg[-\frac{1}{r_0}\Bigg]</math></center> </center> <span style="color:red">Tick. </span>

Therefore, when the LJ potential is zero, the force can be expressed as:

<center><math>F = -\frac{24\varepsilon}{r_0}</math></center> </center> <span style="color:red">Tick. </span>

The equilibrium is reached when the attractive and repulsive forces balance, therefore the force becomes zero. Equating this to the derivative of LJ, the equation becomes

<center><math>24\varepsilon\Bigg[-2\Bigg(\frac{\sigma^{12}}{r^{13}}\Bigg) + \Bigg(\frac{\sigma^6}{r^7}\Bigg)\Bigg]=0</math></center>

Again, since <math>\varepsilon</math> is a constant, the equation can be divided by it on both sides to give

<center><math>-2\Bigg(\frac{\sigma^{12}}{r^{13}}\Bigg) + \Bigg(\frac{\sigma^6}{r^{7}}\Bigg)=0</math></center>

<center><math> -2\Bigg(\frac{\sigma^{6}}{r^{6}}\Bigg) +1 = 0 </math></center>

Therefore, the equilibrium separation is defined by:

<center><math>r_{eq}=\sqrt[6]{2}\sigma</math></center> <span style="color:red">Tick. </span>

The LJ potential at <math>r_{eq}</math> can be expressed as

<center><math>\phi(r_{eq})=4\varepsilon\Bigg(\frac{\sigma^{12}}{(\sqrt[6]{2}\sigma)^{12}}-\frac{\sigma^6}{(\sqrt[6]{2}\sigma)^6}\Bigg)</math></center></center> <span style="color:red">Tick. </span>

<center><math>\phi(r_{eq})=4\varepsilon\Bigg(\frac{1}{4}-\frac{1}{2}\Bigg)</math></center>

Finally, this can be rearranged for the well depeth, <math>\varepsilon</math>


<center><math>\varepsilon= -\phi(r_{eq})</math></center></center> <span style="color:red">Tick. </span>


Integrating the LJ potential gives
<center><math>\int\phi\left(r\right)\mathrm{d}r= - \frac{4\sigma^6\varepsilon(5\sigma-11r^6)}{55r^{11}}</math></center>.

When <math>\sigma=\varepsilon=1.0</math>, the equation becomes

<center><math> \int\phi\left(r\right)\mathrm{d}r= \frac{4(11r^6-5)}{55r^{11}}</math></center>

Since the LJ potential decays rapidly with the separation distance,<math> r</math>, the potential energy generated by atoms with <math>r</math> above a certain distance will be negligible. As this simulation will be modeled with an infinite number pair interactions, finding this cutoff distance can help avoid the need to calculate the potential energy for an infinite number of pair interactions. By trialing different lower limits for the LJ potential integral, a suitable cutoff distance can be found. 

<math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r= -\frac{669}{28160}= -2.48 x 10^-2</math>

<math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r= -8.18 x 10^-3</math>

<math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r= -\frac{32056}{9743085}= -3.29 x 10^-3</math> </center> <span style="color:red">Tick. </span>

=== Periodic Boundary Conditions ===
Since 1ml is equivalent to 1cm<sup>3</sup> and at standard conditions (298K, 1atm), the density of water is 1g/cm<sup>3</sup>, the number of molecules in 1ml of water can be found by diving the mass(m) and Avogdro's number(<math>N_A</math>) by the molecular weight of water (<math>Mw</math>).

<center><math>\frac{mN_A}{Mw}=\frac{1\times 6.02214\times 10^{23}}{18.01528}=3.34\times 10^{22}</math></center></center> <span style="color:red">Tick. </span>

Therefore, the volume of 10,000 molecules of water can be calculated by

<center><math>\frac{10000}{3.34\times 10^{22}}=2.99\times 10^{-19}cm^3</math></center></center> <span style="color:red">Tick. </span>

If an atom located at <math>(0.5,0.5,0.5)</math> is moved along the vector <math>(0.7,0.6,0.2)</math>, the resulting position is <math>(1.2,1.1,0.7)</math>. Applying the periodic boundary of <math>(0,0,0)</math> to <math>(1,1,1)</math>, the resulting position becomes <math>(0.2,0.1,0.7)</math>. </center> <span style="color:red">Tick. </span>

===Reduced Units===
The output form LAMMPS MD simulations are given in reduced units to generate dimensionless forms, allowing for better numerical analysis. The conversion of reduced units to real units is relatively simple; For example, considering the Lennard-Jones parameters for Argon, the Lennard-Jones cutoff is <math>r^*=3.2</math>, and in real units is <math>r=r^*\sigma=3.2\times3.4\times10\times10^-{10}=1.09nm</math>. The well depth is <math>\varepsilon=0.997kJmol^{-1}</math>.  
The reduced temperature of <math>T^*=1.5</math> corresponds to the real value of <math>T= 1.5 \times 120K=180K</math>. </center> <span style="color:red">Good, tick. </span>

==Equilibration==
===Creating the simulation box===

In order to model the system, a simulation box must be constructed in LAAMPS. The simulation box will contain the repeating lattice that will be used to study the system. The folling command defines the distance between each lattice point to be 1.07722(reduced units) in all directions of the 3D simulation. 

<pre>Lattice spacing in x,y,z = 1.07722 1.0722 1.0722</pre>

Since the system defined in the experiment is simple cubic, there will be 1 atom per unit cell. The density can therefore be calculated by:

<center><math>\textstyle \rho=\frac{1}{(1.07722)^3)}=0.799994\approx0.8</math></center>

The atoms coordinates were not randomly generated to eliminated the possibility of two atoms being generated too closely. If the positions of two atoms are too close, the two atoms may collide and the strong repulsive interactions bewteen the two centres will propel the atoms apart, leading to rapid increases in the potential energy of the system, causing the system to be energetically unstable, which is not an accurate representation of the systems under study. <span style="color:red">Good </span>

In a face-cetered cubic system, there are four lattice points per unit cell. Therefore, if the previous system was modelled on a FCC lattice with a lattice point number density of 1.2, the lattice spacing would be 1.49 in reduced units. 

<center>Lattice spacing<math> = \sqrt[3]{\frac{4}{1.2}}= 1.49</math></center>

If the 1000 unit cells were to be generated by the creat_atoms command, 4000 atoms would be generated for a FCC lattice. <span style="color:red">Tick </span>

 
<pre>Mass 1 1.0</pre> 
This indicates that the mass of the single type of atom in the simulation is 1.0. <span style="color:red">Tick </span>
<pre>Pair_style lj/cut 3.0</pre>
"Pair_style" sets the interactions to pairwise interactions. "lj.cut" is the specific pair style used that will compute the standrad 12/6 Lennard-Jones potential with no Colombic interaction."3.0" indicates that the cutoff for all atoms will be at 3.0.<span style="color:red">Tick </span>
<pre>Pair_coeff * *1.0 1.0</pre>
"Pair_coeff" specifies the pairwise force field coefficients, the functional form and parameters used to calculated the potential energy. the double asterisks indicate that the command will be applied to all atoms.<span style="color:red">Tick </span>

Since both <math>x_i</math> and <math>v_i</math> are bith <span style="color:red">*typo </span> being specified, the Velocity-Verlet algorithm will be used for the simulation. <span style="color:red">Tick </span>

To simplify the LAMMPS code, the timestep variable should be attributed a value, i.e. 0.01, that can be called upon by using "${timestep}". This can be excuted by writing:
<pre>variable timestep equals 0.01</pre>

The advantage of defining the variable is that if the timestep were to be changed, rather than having to individually write the numerical timestep over and over again in the script, which may result in human errors especially if the script is long and complex, the timestep would only need to be altered once, at where the variable was first defined. <span style="color:red">Tick </span>

{| class="wikitable" style="text-align: center
|[[File:PressurevsTimex.png|350x350px|thumb|left|Fig. 4a: Pressure vs Time at 0.01 timestep]]
||[[File:TemperaturevTimex.png|350x350px|thumb|center|Fig 4b: Temperature vs Time at 0.01 timestep]]||[[File:File.xpmodls0.01|400x400px|thumb|center|Fig 4c: Total energy vs Time at 0.01 timestep]]
||[[File:Alltimesteps.png|350px|thumb|left|Fig.5: Graph of densities obtained through molecular modelling and Ideal Gas Law]]

|}

It can be seen from the Pressure/Time and Temperature/Time graphs that the system reaches equilibrium at 0.01 timestep as the average pressure and temperature values become constant. Analysis of the output temperature and pressure data from the simulation shows that equilibrium is reached at around t=0.18. The approximate average pressure and temperature after equilibration is 1.26 and 2.61 respectively (both in reduced units).<span style="color:red">Good, agreed </span>

Since it is now known that the simulation reaches equilibrium at 0.01 timestep, longer timesteps can be trialed to find a more time efficient timestep. It is important to ensure that the timestep chosen for later experiments will also afford accurate results therefore simulations were ran at 4 further timesteps (0.0025, 0.0075, 0.01 and 0.15) and the total energy output from them were compared with that of the 0.01 timestep. It can be seen from Fig 5 that the total energy outputs from 0.0025 timestep are in very close agreement with those produced from 0.01 timestep. Systems simulated with 0.0075 and 0.01 both reach equilibrium but the output total energy values are slightly higher than those produced at the slowest timestep, 0.001, thus rendering them inaccurate. Therefore, the largest timestep that can be used to produce acceptable results is 0.025.The worst timestep to use would be 0.015 because the system does not reach equilibrium in the time simulated as indicated by the positive incline in total energy against time shown in Fig.5. <span style="color:red">Tick, what are you going to use and why? </span>

==Running simulations under specific conditions==
The velocity Verlet algorithm allows for the observation of systems that preserve total energy(E), volume(V) and the number of atoms(N), also known as the NVE ensemble. In order to investigate molecular dynamics with different experimental conditions, the systems can be set up with different parameters to mimic those conditions. In this section, systems will be subject to the constant-pressure, constant temperature(NpT) and the constant volume ensemble. 
===Barostat and Thermostat===
The energy of a molecule with N atoms, each with 3 degrees of freedom can be expressed as:

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


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


The temperature will fluctuate after each timestep therefore a velocity-scaling constant, <math>\gamma</math>, is introduced to keep the system isothermal.

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

To solve for <math>\gamma</math>, firstly, sum (5) and (6) to give

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

Since <math>\gamma</math> is a constant, it can be taken out of the summation and the equation can be rearranged to

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

Considering the relationship shown in eqn (5), the equation above can be rearranged to

<math>\frac{3}{2} N k_BT(1+\gamma^2)=  \frac{3}{2} N k_B (T + \mathfrak{T})</math>


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


<math>\gamma^2= \frac{\mathfrak{T}}{T}</math>


<math>\therefore \gamma= \sqrt{\frac{\mathfrak{T}}{T}}</math> <span style="color:red">Tick </span>

However, this is not the best method as rescaling at each step causes discontinuities in the momentum part of the phase space trajectory<ref><nowiki>http://home.gwu.edu/~yxzhao/ResearchNotes/ResearchNote007Thermostat.pdf</nowiki></ref>.

===Examining the Input Script===
<pre> fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_dens2 v_temp2 v_press_2</pre>

the command "fix aves" applies the average operation to the systems. "ave/time 100 1000 100000" computes the global time-averages of the variables inputted(density, pressure and temperature). "100 1000 100000" specify the timesteps the input values will be used in order to contribute to the average. Therefore for this simulation, the final average will be calculated using values at timestep 100, 200, 300....100000. In total, 1000 values will be used to calculate the global average.

<pre>run 100000</pre>

This command tells LAMMPS to run the simulation over 100,000 timesteps. Since a timestep of 0.0025 will be used, 250 time units simulated. <span style="color:red">Why do you use variables? </span>

===Plotting the Equations of State===
[[File:Density vs Temp graphxxx.png|450px|thumb|right|Fig.6:Density vs Time Graph at p=2.4 and p=2.6]]

The NpT ensemble (also known as the isothermal-isobaric ensemble) describes systems with a fixed number of particles, N, in contact with a thermostat at temperature <math>T</math> and a bariostat <span style="color:red">*typo </span> at pressure <math>p</math>. The system exchanges heat with the thermostat and pressure with the bariostat <span style="color:red">*barostat </span>. Therefore, the total energy,<math> E</math>, and volume, <math>V</math>, both fluctuate at thermal equilibrium. Simulations were conducted at two pressures (2.4, 2.6 in reduced units), and at 5 temperatures(2, 2.5, 3, 3.5, 5). In order to compare the results obtained through the Lennard Jones simulations to theoretical Idea Gas Law (<math>PV=Nk_BT</math>) results, a plot of density and temperature was constructed at both pressures(Fig.6). <span style="color:red">Good </span>

In plot of density against temperature (Fig.6), the error bars of the data generated do not overlap, indicating that the density outputs from the simulations are statistically significant, thus analysis of the graphs can be carried out with confidence.

The graph shows density decreases with increasing temperature. This is an expected phenomenon as the system will increase in volume to compensate for the increase in kinetic energy due to increase in temperature. The Ideal Gas Law can be rearranged to show the inverse relationship between temperature and density.

<center><math>\frac{N}{V}=\frac{P}{k_BT}</math>  (7)</center>


[[File:Npt vs Ideal Gas Lawx.png|450px|thumb|left|Fig.7: Graph of densities obtained through molecular modelling and Ideal Gas Law]]

Additionally, Fig.6 also shows that at constant temperature, density increases with increasing pressure. This is again expected because if increasing pressure was applied to a container with fixed number of particles, the volume would decrease, hence decreasing the density of the system. Again, this relationship can be justified by (7), where <math>\frac{N}{V}</math> is directly proportional to <math>P</math>. 


It can be seen from Fig.7 that the densities obtained through the Lennard-Jones molecular modelling is significantly lower than the densities obtained through Ideal Gas Law calculations. This is expected as the Ideal Gas Law assumes that molecules do not interact with each other, thus there will be no repulsion factor as there is in the LJ model, thus making the particles in the Ideal Gas system to be more compressible. This means that the molecules occupy negligible volume therefore for a given volume, the density will be higher. In the LJ model, the particles interact and repulsion occurs when <math>r</math> is small, hence the density at a given volume will be lower. <span style="color:red">Tick </span>

The Ideal Gas Law is most accurate in describing small gases at low pressure and low temperature

As seen from Fig.7, the largest discrepancies between the LJ simulation densities and the Ideal Gas densities occur at low temperature and high pressures(the gradient of Ideal Gas curve is higher at pressures whilst the gradient of the LJ simulation curve results do not differ as visibly). This is because the system behaves more like an ideal gas at high temperature and low pressure. The reason for this observation is that intermolecular interactions become less limiting when the thermal kinetic energy of the system is high i.e. high temperature; At lower pressure, assuming the volume is large, the intermolecular distance becomes larger therefore the density is less affected by the separation, <math> r</math>. <span style="color:red">Tick, very good </span>

In conclusion, discrepancies increase with an increase in pressure and a decrease in temperature.

==Heat Capacity Calculation==
The NVT(canonical) ensemble describes systems have a fixed number of particles, N, and volume, V, that are in contact with a thermostat at temperature T. In this section, the heat capacities under the NVT ensemble were calculated for two different pressures (p=0.2 and p=0.8) and at five different temperatures each (T=2, 2.2, 2.4, 2.6, 2.8) by using the equation

<center><math>C_V = \frac{\partial E}{\partial T} = \frac{\mathrm{Var}\left[E\right]}{k_B T^2} = N^2\frac{\left\langle E^2\right\rangle - \left\langle E\right\rangle^2}{k_B T^2}</math>(8)</center>

To run the simulation under NVT conditions, the following code was used.

<pre>### DEFINE SIMULATION BOX GEOMETRY ###
variable dens equal 0.2
variable T equal 2.0
variable timestep equal 0.0025

lattice sc ${dens}
region box block 0 15 0 15 0 15
create_box 1 box
create_atoms 1 box

### DEFINE PHYSICAL PROPERTIES OF ATOMS ###
mass 1 1.0
pair_style lj/cut/opt 3.0
pair_coeff 1 1 1.0 1.0
neighbor 2.0 bin

### ASSIGN ATOMIC VELOCITIES ###
velocity all create ${T} 12345 dist gaussian rot yes mom yes

### SPECIFY ENSEMBLE ###
timestep ${timestep}
fix nve all nve

### THERMODYNAMIC OUTPUT CONTROL ###
thermo_style custom time etotal temp press
thermo 10

### RECORD TRAJECTORY ###
dump traj all custom 1000 output-1 id x y z

### SPECIFY TIMESTEP ###

### RUN SIMULATION TO MELT CRYSTAL ###
run 10000
unfix nve
reset_timestep 0

### BRING SYSTEM TO REQUIRED STATE ###
variable tdamp equal ${timestep}*100
variable pdamp equal ${timestep}*1000
fix nvt all nvt temp ${T} ${T} ${tdamp} 
run 100000
reset_timestep 0

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp vol atoms
variable temp equal temp
variable temp2 equal temp*temps
variable energy equal etotal
variable energy2 equal etotal*etotal
variable atoms2 equal atoms*atoms
variable volume equal vol
fix aves all ave/time 100 1000 100000 v_temp v_temp2 v_energy v_energy2 v_volume
unfix nvt
fix nve all nve
run 100000

variable avetemp equal f_aves[1]
variable aveheatcap equal ${atoms2}*(f_aves[4]-(f_aves[3])(f_aves[3])f_aves[2]
variable avevolume equal f_aves[5]
variable aveheatcapv equal ${aveheatcap}/${avevolume}

print "Averages"
print "--------"
print "Temperature: ${avetemp}"
print "Total energy: ${etotal}"
print "Volume: ${avevolume}"
print "Heat Capacity: ${aveheatcap}"
print "Heat Capacity over V: ${aveheatcapv}"

</pre>

The simulation box was defined with parameters of <math>15\times 15\times 15</math> which resulted in a total of 3375 atoms. The volume for p=0.2 was and for p=0.8 was

In order to compare how heat capacities vary at different volumes, the heat capacity must be divided by the volume of each system thus <math>\frac{C_V}{V}</math> was plotted against temperature for each system (Fig.8).

[[File:Cvtgraphx.png|400px|thumb|left|Fig.8: C<sub>v</sub>/V against temperature for densities p=0.2 and p=0.8]]

The specific heat capacity defines how much heat is needed to raise the temperature of a unit mass by a unit degree of temperature. Fig.8 shows that the heat capacity  has an inversely proportional relationship with temperature (<span style="color:red">so why the linear plot? </span>), in agreement with the definition of equation (8). This can be explained by the band-gap (<span style="color:red">density of states </span>) between the lattice energy levels decreasing with increasing temperature. Therefore, less energy is required to heat the system; relating this back to equation (8), where it is shown that <math>C_v</math> is directly proportional to energy, increasing the temperature will decrease the heat capacity. 

The results also show that heat capacity decreases with decreasing density. This is expected as high density leads to lower volume for a fixed number of atoms so it its expected that the heat required to raise the temperature of the system will be lower as particles are closer together compared to a lower density system with the same number of atoms where the volume is larger.<span style="color:red">Tick </span>

==Structural properties and the radial distribution function==
The radial distribution function, g(r), describes the probability of finding a nearest neighbour from a reference particle. Thus it allows for the identification of the phase of the system as it provides information about the variation of density from a central reference particle. <span style="color:red">Tick</span>
[[File:RDFvsrxx.png ‎|500px|thumb|right|Fig 9 g(r) against distance for solid, liquid and gaseous phases ]]
The radial distribution function was plotted for solid, liquid and gas for a Lennard-Jones material(Fig.9).The parameters were chosen from the Lennard-Jones phase diagram reported by Hansen and Verlet<ref><nowiki>http://journals.aps.org/pr/abstract/10.1103/PhysRev.184.151</nowiki></ref>.. 
{| class="wikitable"
!Phase
!Density
!Temperature
|-
|Solid
|1.2
|1.0
|-
|Liquid
|0.8
|1.2
|-
|Vapour
|0.1
|1.6
|}

In all phases, a common factor was that g(r)=0 until r=0.9(reduced units). This is because at distances shorter than atomic diameter, strong repulsive forces causes the g(r) to be zero. <span style="color:red">Tick</span>

The qualitative appearance for the three systems have clear differences. The solid system has the most number of peaks followed by liquid then vapour. This appearance is due to the fact that the peaks indicate the density around each reference atom therefore, the phase with the most regular structure,solids, will have the most discrete peaks. The solid peaks show a general decrease in amplitude as the r increases with some fluctuations (<span style="color:red">Why does it fluctuate?</span>). The first three peaks in the solid RDF represent the first three coordination shells of the solid; the first peak represents the nearest neighbour of the reference molecule, the second peak is the second nearest neighbour and the third peak is the third nearest neighbour. There is zero probability of finding a particle in regions between the peaks as molecules in a solid cannot move. The distance between the zero probability minima gives the lattice spacing therefore the lattice spacing is defined by the second nearest neighbour; it is approximately 1.49 in reduced units. The appearance of the peaks in the solid system also indicate that the solid phase has long and short range order; the three tall peaks indicate short range order while the fluctuating smaller peaks indicate long range order.  

The liquid g(r)has three peaks which decrease in amplitude with increased interatomic distance. The peaks are less narrow that those in the solid phase thus suggesting the system to have more disorder. The amplitude of the peaks decrease with separation, showing order decreases with separation due to the random Brownian motion of of particles. The low number of peaks in the radial distribution function suggests that there is only short range order in the system, specifically between the first three neighbouring atoms.  

For gases, a single broad peak is observed at around r=1. The broad peak suggests that there is high disorder in the gaseous phase thus there is no short nor long range order. <span style="color:red">Tick</span>

[[File:Grintgraphxx.png|400px|thumb|left|Fig.10: Expansion of the integral of g(r) for solid phase]]

Plotting the integral of the RDF against distance will give information about the coordination numbers of the neighbours of the reference particle in a solid. The inflection points called out in the Fig.10 show the coordination number of the three nearest neighbours of the reference particle. Since the solid system was modeled on a FCC lattice, the reference particle is expected to have 12 nearest neighbours (this is confirmed by the point at r=1.225). The next inflection point at r=1.625 has a g(r) integral of <math>\approx18</math>. Since the plot is a running integral of g(r), the coordination number of the second nearest neighbours is <math>18-12=6</math>. Thus, the coordination number of the third nearest neighbours is <math>42-18=24</math>.<span style="color:red">Tick, nice</span>

==Dynamical properties and the diffusion coefficient==

===Mean Squared Displacement===

The [https://en.wikipedia.org/wiki/Mean_squared_displacement mean squared displacement] measures the average distance a particle travels in a given time from a reference position, due to random motion. It is defined:
<math>{\displaystyle {\rm {MSD}}\equiv \langle (x-x_{0})^{2}\rangle ={\frac {1}{N}}\sum _{n=1}^{N}(x_{n}(t)-x_{n}(0))^{2}}</math>

In this equation, <math>x_n(t)-x_n(0)</math> represents the vector distance traveled by a particle n over time interval of t. The square of the displacement is then averaged over N particles by summing n from 1 to N and then dividing by N.

{| class="wikitable"
| Smaller system simulation(3375 atoms) || Larger system simulation(1 million atoms)
|-
|[[File:Solmsdxx.png ‎|400px|thumb|Fig.11a:Solid Simulation at d=01.2, T=1.0]]  ||[[File: Solsimxx.png|400px|thumb|Fig.11b:Solid Simulation at d=0.8, T=1.2]] 
|-
|[[File:Liqsimmsdxx.png|400px|thumb|Fig.12a:Liquid Simulation at d=0.8, T=1.2]]  ||[[File:Largeliqxx.png|400px|thumb|Fig.12bLiquid Simulation at d=0.8, T=1.2]] 
|-
|[[File:Vapsimx.png|400px|thumb|Fig.13a:Vapour Simulation at d=0.1, T=1.6]]  ||[[File:Largervapsimx.png|400px|thumb|Fig.13b:Vapour Simulation at d=0.1, T=1.6]] 
|-

|}

The plots generated are shown above. In a solid system, the particles are bound to their position so there is not sufficient kinetic energy to reach a diffusive behaviour and this is shown in Fig.11a, as indicated by the nearly constant value of the MSD after around timestep 200.

The liquid plot(Fig.12) shows a clear linear relationship (<math>R^2=0.9999</math>) between MSD and time. This shows that the particles in the system move acccording to Brownian motion.

[[File:Vap2000simmsdxx.png|400px|thumb|Fig.14:Vapour Simulation at d=0.1, T=1.6 for smaller system at timesteps above 2000]] 

In the plot for the gaseous phase(Fig.13), relationship between MSD and time is initially curved, followed by a near linear relationship(shown by Fig.14 when only consdiering the relationship between MSD and time above timestep 2000). The inital curvation is due to the fact that atoms are placed randomly to simulate the system. Since it is a gaseous system, the interatomic separations are high. This results in a lower frequency of collisions thus the overall velocity the atoms travel at appears to be constant. Therefore, by definition, the distance travelled per unit time is also constant, hence MSD is proportional to <math>t^2</math>. As more time is simulated, more collisions occur and the atomic motion becomes more consistent with Brownian motion thus the MSD/Time graph tends to linearity. <span style="color:red">Takes a while for the system to act like a random walk</span>

The simulations performed at 3375 atoms and a million atoms both generated MSD/Time graphs that show the same general trends for each phase. In the solid simulations, there is less fluctuations in the MSD of the larger system, suggesting the system to be more stable for solids. The difference in the output values of simulation performed on 3375 atoms and the 1 million atoms are relatively small(around 0.02 difference in MSD). Since the phase parameters used for the larger simulation were not stated, the discrepancy could also be caused by the different simulation conditions used. However, this discrepancy could also be due to the difference in number of atoms used; if this were the case, it can concluded that there is no need to use such large systems as the difference in output data is not significantly difference thus the tradeoff for longer simulation times is not worthwhile. <span style="color:red">Tick</span>

Einstein showed that the the self-diffusion coefficient is related to the mean-squared displacement:

<math>D=\frac{1}{6}\frac{\delta\langle r^2\rangle}{\delta t}</math>

Substituting simulation values for the 3375 atoms system, the diffusion coefficients in reduced units are

Gas:<math>D=\frac{1}{6}\times \frac{0.0226}{0.002}=1.88</math>
Solid: <math>D=5.00\times10^{-5}</math>
Liquid:<math>D=8.33\times10^{-2}</math>

For the larger 1 million atoms system

Gas:<math>D=\frac{1}{6}\times \frac{0.0326}{0.002}=3.03</math>
Solid: <math>D=4.17\times10^{-7}</math>
Liquid:<math>D=8.33\times10^{-2}</math>

The diffusion coefficients are very similar for the smaller and larger model.<span style="color:red">Tick</span>

===Velocity Autocorrelation Function===
The motion of a 1D harmonic oscillator is given by

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

Velocity is the first derivative of positions and time:

<center><math>v(t)=\frac{dr}{dt}</math></center>

Therefore, the position of a 1D harmonic oscillator at a specific velocity can then be found by substituting

<center><math>v(t) = \frac{d(Acos(\omega t +\phi)}{dt}</math></center>

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

At timestep t+\tau , the velocity will therefore be:

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

Substituting this into the velocity correlation function, 

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


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

Considering the trigonometric identity <math>sin(x+y)=sin(x)cos(y)+cos(x)sin(y)</math>. The equation can therefore be written as:
<center><math>C(\tau) = \frac{\int_{-\infty}^{\infty}sin(\omega t+\phi)[sin(\omega t+ \phi)cos(\omega \tau)+cos(\omega t+\phi)sin(\omega \tau)]\mathrm{d}t}{\int_{-\infty}^{\infty} sin^2(\omega t+\phi)\mathrm{d}t}</math></center>

Since <math>\omega</math> and <math>\tau</math> are both constants, they can be taken out of the integral:

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

Since <math>sin(x)</math> is an odd function, evaluating its integral by limits that are symmetric across the origin,<math>-\infty</math> to <math>\infty</math>, will be zero. Thus the right hand side of the equation will equal zero and the equation becomes:

<center><math>C(\tau)=cos(\omega \tau)</math></center> <span style="color:red">Good</span>


[[File:VASCallgx.png|400px|thumb|left|Fig.15: Graph of VACF for solid, liquid and the Simple Harmonic Oscillator]]

A graph of VACF for the Lennard Jones liquid and solid simulations was plotted along with the VACF graph of the Simple Harmonic Oscillator. It can be observed in the graph that there is a maximum difference between the <math>\nu(t)</math> and <math>\nu(t+\tau)</math> at the minima of the solid and liquid plots. In the solid,, this indicates that atoms have collided and are changing directions of motion. In the liquid, this phenomena can be attributed to the atom colliding with their "solvent cage" and then changing direction. The minima of the solid is lower than the minima of the liquid as there are stronger interatomic forces in the solid. <span style="color:red">Tick, nice description</span>

The liquid VACF has one weak oscillation as the particles in the system only interact with their closest neighbours whereas in the solid, a series of oscillations are seen as the ordered lattice allows atoms to vibrate forwards and backwards. However, the amplitude of the oscillations become smaller.

For the Simple Harmonic Oscialltor, the appearance of the VACF is very different. This is due to the fact that the oscillator does not consider the interactions with other atoms in the system, and instead models a system with no energy loss, thus the velociy never decreases.

To find the running integral of under the VACF graphs, the trapezium rule was employed as the diffusion coefficient can also be defined 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>

The last point of data on the running integral of VACF graphs gives the total integral of the VACF graphs. 

For the smaller system the results are as follows
[[File:Gasvacfxx.png|400px|thumb|left|Fig.16a:Running intergral of VACF against time for Gas in small system]]

[[File:Liqvacfxx.png ‎|400px|thumb|left|Fig.16b:Running intergral of VACF against time for Liquid in small system]]

Solid:<math>D=4.63\times 10^{-4}</math>

Liquid:<math>D=7.93\times10^{-2}</math>

Gas:<math>D=2.35</math> <span style="color:red">Tick</span>

[[File:Largevapvacfxx.png  ‎|400px|thumb|right|Fig.18a:Running intergral of VACF against time for Gas in large system]]

[[File:Largeliqvacfxx.png ‎|400px|thumb|right|Fig.18b:Running intergral of VACF against time for Liquid in large system]]

For larger system, the diffusion coefficients are

Solid:<math>D=5.67\times 10^{-5}</math>

Liquid:<math>D=7.93\times 10^-2</math>

Gas:<math>D=2.35</math>

There is little discrepancy with the diffusion coefficients obtained from the larger and smaller system.

The coefficients calculated from the VACF method show the same trend as the Mean Displacement Squared approach. The diffusion coefficient was largest for gases, then liquids and smallest for solids. Neither methods were particularly affected by the sample size of the simulation box. The two methods produced <math>D</math> of similar values and of similar magnitudes for the liquid and gas simulations but the coefficients from the VACF method were larger significantly larger than those obatined thorugh the MSD method. Therefore, it must be noted that using the trapezium rule to integrate the VACF is likely to accumulate errors over time, especially for solids as solid VACF curves are non periodic. The accuracy of the method could be increased by using smaller timesteps but this would have implications on the simulation time. Another way to refne the VACF method would be to use a more sophisticated integration method. <span style="color:red">Good comparison</span>

==References==