2016-02-25T16:15:08Z

Wkl13: /* Task 1 */

=='''<big> Introduction </big>'''==
===Task 1===
'''Open the file HO.xls. In it, the velocity-Verlet algorithm is used to model the behaviour of a classical harmonic oscillator. Complete the three columns "ANALYTICAL", "ERROR", and "ENERGY": "ANALYTICAL" should contain the value of the classical solution for the position at time <math>t</math>, "ERROR" should contain the ''absolute'' difference between "ANALYTICAL" and the velocity-Verlet solution (i.e. ERROR should always be positive -- make sure you leave the half step rows blank!), and "ENERGY" should contain the total energy of the oscillator for the velocity-Verlet solution. Remember that the position of a classical harmonic oscillator is given by <math> x\left(t\right) = A\cos\left(\omega t + \phi\right)</math> (the values of <math>A</math>, <math>\omega</math>, and <math>\phi</math> are worked out for you in the sheet).'''

For a classical harmonic oscillator, each particles will experience a force. The analytical position x(t) can be described by the function

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

The error is the absolute of the differences between the Verlet algorithm and the analytical x. The graph below is error plotted against time.

The total energy is the sum of the kinetic and potential energy.
<ref name="Principles of Physics">D. Halliday, R. Resnick, J. Walker, ''Principles of Physics'', John Wiley & Sons, 9th edn., 2011, ISBN:0470561580</ref>
:<math>E(x) = \frac{1}{2} m v^2+ \frac{1}{2} k x^2</math>

For timestep 0.1,
<gallery widths=450px heights=150px>
File:X(t)againstt.png|x(t) against time graph
File:Eagainstt.png|Energy against time graph
File:Erroragainstt.png|Error against time
</gallery>

The total energy is almost constant. Analytically, it should be constant as energy is always conserved. Since the algorithm approach is used, there are sips in between, hence leading the energy fluctuates slightly.

===Task 2===
'''For the default timestep value, 0.1, estimate the positions of the maxima in the ERROR column as a function of time. Make a plot showing these values as a function of time, and fit an appropriate function to the data.'''
{| class="wikitable"
! local maxima
! 1st
! 2nd
! 3rd
! 4th
! 5th
|-
!positions
| 7.58E-04||2.01E-03||3.30E-03||4.60E-03||5.91E-03
|-
! time
| 2.0
| 4.9
| 8.0
| 11.1
| 14.2
|}
[[File:Errorplotmaxagainstt.png|none|300 × 200|left|Error plot max against time]]

A positive linear correlation is found. All the points satisfy into y =0.004x-0.00007 with <math> R^{2}= 1 </math>.

===Task 3===
'''Experiment with different values of the timestep. What sort of a timestep do you need to use to ensure that the total energy does not change by more than 1% over the course of your "simulation"? Why do you think it is important to monitor the total energy of a physical system when modelling its behaviour numerically?'''

Different values of timestep have been used to explore the percentage of total energy change. It is found x< 0.2 does not change the total energy by more than 1%. Analytically, by the conservation of energy, energy is always conserved. However, in the algorithm approach, the Taylor expansion only expands to the power of 3, in which the calculated value is not exact, hence there is a difference in energy. 1% is the benchmark.

===Task 4===
'''For a single Lennard-Jones interaction, <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>, find the separation, <math>r_0</math>, at which the potential energy is zero. What is the force at this separation? Find the equilibrium separation, <math>r_{eq}</math>, and work out the well depth (<math>\phi\left(r_{eq}\right)</math>). Evaluate the integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math> when <math>\sigma = \epsilon = 1.0</math>.'''

The potential energy is given by the Lennard-Jones equation<ref>ChemWiki, http://chemwiki.ucdavis.edu/Core/Physical_Chemistry/Physical_Properties_of_Matter/Atomic_and_Molecular_Properties/Intermolecular_Forces/Specific_Interactions/Lennard-Jones_Potential, (accessed February 2016)</ref>

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

when <math>\phi\left(r\right) = 0</math>,

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

<math>\frac{\sigma^{12}}{r_0^{12}} = \frac{\sigma^{6}}{r_0^{6}}</math>

<math>r_0 ^{6} = \sigma^{6} </math>

<math>r_0 = \sigma</math>

For a single interaction,

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

Consider the Lennard-Jones interaction is the only source of potential energy,

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

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

<math>= 48\epsilon \frac{\sigma^{12}}{r^{13}} - 24\epsilon\frac{\sigma^{6}}{r^{7}}</math>

for <math> r = \sigma </math>,

<math>\mathbf{F} = 48\epsilon \frac{\sigma^{12}}{\sigma^{13}} - 24\epsilon\frac{\sigma^{6}}{\sigma^{7}}</math>

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

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

At equilibrium <math>r_{eqm}</math>, the potential energy reaches the minimum, hence <math> \frac{\mathrm{d}U}{\mathrm{d}r}=0 </math>.

<math>\frac{\mathrm{d}U}{\mathrm{d}r} = - 4\epsilon \left( \frac{-12\sigma^{12}}{r_{eqm}^{13}} + 6\frac{\sigma^6}{r_{eqm}^7} \right) = 0</math>

<math>\frac{-12\sigma^{12}}{r_{eqm}^{13}} + 6\frac{\sigma^6}{r_{eqm}^7} = 0</math>

<math>\frac{12\sigma^{12}}{r_{eqm}^{13}} = \frac{6\sigma^6}{r_{eqm}^7}</math>

<math>2\sigma^{6} = {r_{eqm}^6}</math>

'''<math> r_{eqm} = ({2})^{1/6} \sigma </math>'''

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

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

<math> = 4\epsilon \left( \sigma^{12} \int_{2\sigma}^\infty \frac{1}{r^{12}} \mathrm{d}r - \sigma^{6} \int_{2\sigma}^\infty \frac{1}{r^{6}}\mathrm{d}r \right)</math>

===Task 5===
'''Estimate the number of water molecules in 1ml of water under standard conditions. Estimate the volume of <math>10000</math> water molecules under standard conditions.'''

Under standard conditions [101 kPa, 298K] ,

<math>\text {density of water} = 0.997 g ml^{-1}</math>

<math>Mr \left(water \right) = 18.02 g mol^{-1}</math>

<math> \text {no of moles} = \frac {mass}{Mr} = \frac {\text {no of molecules}}{\text {Avogadro's constant}} </math>

===Task 6===
'''Consider an atom at position <math>\left(0.5, 0.5, 0.5\right)</math> in a cubic simulation box which runs from <math>\left(0, 0, 0\right)</math> to <math>\left(1, 1, 1\right)</math>. In a single timestep, it moves along the vector <math>\left(0.7, 0.6, 0.2\right)</math>. At what point does it end up, ''after the periodic boundary conditions have been applied''?'''

<math>\left(0.5, 0.5, 0.5\right)+\left(0.7, 0.6, 0.2\right)</math>

<math>= \left(1.2, 1.1, 0.7\right)</math>

(Using periodic boundary condition*)

<math>= \left(0.2, 0.1, 0.7\right) </math>

*Periodic boundary condition: if particle moves out of the box <math>\left(1, 1, 1\right)</math>, the box will move back and particle will appear in the box again.

===Task 7===
'''The Lennard-Jones parameters for argon are <math>\sigma = 0.34\mathrm{nm}, \epsilon\ /\ k_B= 120 \mathrm{K}</math>. If the LJ cutoff is <math>r^* = 3.2</math>, what is it in real units? What is the well depth in <math>\mathrm{kJ\ mol}^{-1}</math>? What is the reduced temperature <math>T^* = 1.5</math> in real units?'''

<math> r = \frac{r^*}{\sigma} </math>

<math> = 3.2 \times 0.34 \mathrm{nm} </math>

<math> = 1.09 \mathrm{nm} </math>

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

<math>= 4 \times 120\ k_B\times \left( 3.2^{-12} {{-}} 3.2^{-6} \right) </math>

<math> = -6.17 \times 10^{-24} \mathrm{kJ mol^{-1}}</math>

<math> T^* = \frac{k_B T}{\epsilon} </math>

<math> T = \frac{T^* \epsilon} {k_B}</math>

<math> = 1.5 \times 120 \mathrm{K} </math>

<math> = 180 \mathrm{K}</math>

== Equilibration ==

===Task 1===

'''Why do you think giving atoms random starting coordinates causes problems in simulations? Hint: what happens if two atoms happen to be generated close together?'''

Atoms with random starting coordinates have a chance of superposing each other, which is impractical in reality. No two atoms can rest on the same space.

When two atoms are too close (the separation distance is small), according to the Lennard-Jones Potential, the potential energy is very high, as the atoms are repulsing each other. Therefore, the simulation does not fit into real case.

===Task 2===

'''Satisfy yourself that this lattice spacing corresponds to a number density of lattice points of <math>0.8</math>. Consider instead a face-centred cubic lattice with a lattice point number density of 1.2. What is the side length of the cubic unit cell?'''

<math> \text{number density} = \frac{\text{no of atoms}}{\text{volume}}</math>

<math> \text{volume} = \text{length}^{3} </math>

For a simple cubic lattice,

<math> \text{no of atoms} = 8 \times {\frac{1}{8}}= 1 </math>

'''<math> \text{number density} = \frac{1}{1.07722^{3}} = 0.800 \left(3 s.f.\right)</math>'''

For a face-centred cubic lattice,

<math> \text{no of atoms} = 8 \times {\frac{1}{8}}+ 6 \times {\frac{1}{2}}= 1 + 3 = 4</math>

<math> \text{volume} = \text{length}^{3} = \frac{\text{no of atoms}}{\text{number density}}= \frac{4}{1.2}</math>

<math> \text{length}^{3} = 3.\dot{3} </math>'''

'''<math> \text{length} = 1.49 \left(\text{in reduced units}\right)</math>'''

===Task 3===

'''Consider again the face-centred cubic lattice from the previous task. How many atoms would be created by the create_atoms command if you had defined that lattice instead?'''

<math> \text{no of atoms per FCC lattice} = 4 </math>

1 box creates 1000 unit cells

<math> \text{no of atoms created} = 4 \times 1000 = 4000</math>

===Task 4===

'''Using the [http://lammps.sandia.gov/doc/Section_commands.html#cmd_5 LAMMPS manual], find the purpose of the following commands in the input script:'''

<pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

{|-
|<pre>
mass 1 1.0
</pre>
|This is a mass command. '1' refers to the atom type while '1.0' refers to the mass value.
|-
|<pre>
pair_style lj/cut 3.0
</pre>
|This command implies only considering the Lennard-Jones potential up to 3.0 atomic units.
|-
|<pre>
pair_coeff * * 1.0 1.0
</pre>
|This is a pair coefficient command. The asterick (*) with no values refer all types of atom from 1 to N. The 1.0 1.0 refers to the atom type '1.0' interacting with '1.0'.
|}
<ref name="LAMMPS manual">LAMMPS manual, http://lammps.sandia.gov/doc/Section_commands.html#cmd_5, (accessed February, 2016)</ref>

===Task 5===

'''Given that we are specifying <math>\mathbf{x}_i\left(0\right)</math> and <math>\mathbf{v}_i\left(0\right)</math>, which integration algorithm are we going to use?'''

Velocity-Verlet algorithm <math>\mathbf{v}_i\left(0\right)</math>

For classical Verlet algorithm <math>\mathbf{x}_i\left(0\right)</math>,the initial conditions only require <math>\mathbf{x}_i\left(0\right)</math> and <math>\mathbf{x}_i\left(-\delta t\right)</math>;

while velocity Verlet algorithm <math>\mathbf{v}_i\left(0\right)</math> require <math>\mathbf{x}_i\left(0\right)</math> and <math>\mathbf{v}_i\left(0\right)</math> only.

===Task 6===

'''Look at the lines below.'''
<pre>
### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(100/${timestep})
variable n_steps equal floor(100/0.001)
timestep ${timestep}
timestep 0.001

### RUN SIMULATION ###
run ${n_steps}
run 100000
</pre>
'''The second line (starting "variable timestep...") tells LAMMPS that if it encounters the text ${timestep} on a subsequent line, it should replace it by the value given. In this case, the value ${timestep} is always replaced by 0.001. In light of this, what do you think the purpose of these lines is? Why not just write:'''
<pre>
timestep 0.001
run 100000
</pre>'''

The timestep command fixes the time interval for each run. The advantage of using the lines in the top box is flexibility. When the timestep is varied in the future experiment, the run will change automatically, and does not require manual input which could avoid human errors when running the LAMMPS simulation.

===Task 7===

'''make plots of the energy, temperature, and pressure, against time for the 0.001 timestep experiment (attach a picture to your report). Does the simulation reach equilibrium? How long does this take? When you have done this, make a single plot which shows the energy versus time for all of the timesteps (again, attach a picture to your report). Choosing a timestep is a balancing act: the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality; short timesteps, however, mean that the same number of simulation steps cover a shorter amount of actual time, and this is very unhelpful if the process you want to study requires observation over a long time. Of the five timesteps that you used, which is the largest to give acceptable results? Which one of the five is a ''particularly'' bad choice? Why?'''

For timestep = 0.001
<gallery widths = 450px heights = 350px>
File:Energty against time graph.tif|Energy against time graph
File:Temperature against time.tif|Temperature against time graph
File:Pressure agasint time.tif|Pressure against time graph
</gallery>

The simulation reached equilibrium, as the energy is almost constant. The fluctuation can be explained by the algorithm used, which only approximates up to the third term of Taylor's expansion. Equilibrium is reached after 0.3 time units.

A plot of energy against time with different timesteps is shown below:

[[File:Energy against time graph.tif|none|800x500 px|left|Energy against time for various timestep]]

Timestep of 0.015 does not reach equilibrium, as the energy is increasing.It is a "particularly" bad choice. Theoretically, by conservation of energy, the total energy should be a constant, which does not depend on the timestep. The timestep we choose for the experiment should give out the same energy. Among the other simulations, timestep of 0.0025 and 0.001 gives similar energies, which indicates the Verlet-algorithm gives a good approximation in these cases. In addition, in running the simulations,we always aim at longer observation time, which could go for longer time. Therefore, timestep of 0.0025 is the largest to give a result which looks like independent of the time.

== Running simulations under specific conditions ==

===Task 1===

'''Choose 5 temperatures (above the critical temperature <math>T^* = 1.5</math>), and two pressures (you can get a good idea of what a reasonable pressure is in Lennard-Jones units by looking at the average pressure of your simulations from the last section). This gives ten phase points — five temperatures at each pressure. Create 10 copies of npt.in, and modify each to run a simulation at one of your chosen <math>\left(p, T\right)</math> points. You should be able to use the results of the previous section to choose a timestep. Submit these ten jobs to the HPC portal. While you wait for them to finish, you should read the next section.'''

<math>T^* = 1.6, 1.8, 2.0, 2.2, 2.4</math> and <math> p^*=2.0, 3.0 </math> (fluctuate within this range as shown in previous part) were chosen.
As explained in previous section, timestep of <math> 0.0025 </math> was chosen.

===Task 2===

'''We need to choose <math>\gamma</math> so that the temperature is correct <math>T = \mathfrak{T}</math> if we multiply every velocity <math>\gamma</math>. We can write two equations:

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

'''Solve these to determine <math>\gamma</math>.'''

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

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

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

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

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

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

===Task 3===

<pre>
fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_dens2 v_temp2 v_press2
run 100000
</pre>
'''Use the Manual page<ref name="manual page">LAMMPS manual page, http://lammps.sandia.gov/doc/fix_ave_time.html, (accessed February, 2016)</ref>to find out the importance of the three numbers ''100 1000 100000''. How often will values of the temperature, etc., be sampled for the average? How many measurements contribute to the average? Looking to the following line, how much time will you simulate?'''

The lines indicate:
{|-
|<pre>
100
</pre>
|Nevery, refers to the use values every 100 timesteps.
|-
|<pre>
1000
</pre>
|Nrepeat, refers to the number of times to use the input values for calculating the averages.
|-
|<pre>
100000
</pre>
|Nfreq, calculates the averages every 100000 timesteps.
|}

The values of the temperature is recorded every 100 step.

The values are sampled every 100000 timesteps.

1000 measurements contribute to the average.

1000000 will be used to complete the simulation.

===Task 4===

'''When your simulations have finished, download the log files as before. At the end of the log file, LAMMPS will output the values and errors for the pressure, temperature, and density <math>\left(\frac{N}{V}\right)</math>. Use software of your choice to plot the density as a function of temperature for both of the pressures that you simulated. Your graph(s) should include error bars in both the x and y directions. You should also include a line corresponding to the density predicted by the ideal gas law at that pressure. Is your simulated density lower or higher? Justify this. Does the discrepancy increase or decrease with pressure?'''

The ideal gas model can also be used to approximate the gases, which assumes

(1) negligible intermolecular forces between the gas molecules

(2) perfectly elastic collisions between the molecule and with the container walls is the only interactions

(3) negligible volume occupied by the molecules compared to container

It can be described mathematically as shown below

<math> pV = nRT </math>

where p is pressure, V is the volume, n is number of moles, R is the gas constant and T is temperature.

Since <math> nR = Nk_B </math>, where N is the number of atoms and <math> k_B </math> is the Boltzmann's constant.

<math> k_B =1 (\text{in reduced units})</math>

<math> pV = NT </math>

<math> \frac{N}{V} = \frac{p}{T} </math>

<math>\text {density} \rho =\frac{p}{T} </math>

Apart from running simulation, density can also be calculated from the the ideal gas law by dividing pressure by temperature.

The error bars are based on standard deviations <math>\sigma </math> from the simulation, while <math> \sigma \left( \rho \right) </math> was calculated using the following equation <ref name= "Measurements and their Uncertainties A Practical Guide to Modern error Analysis" > I. HUGHES, T. HASE, ''Measurements and their Uncertainties A Practical Guide to Modern error Analysis'', Oxford University Press, OXFORD, 2010, ISBN:9780199566334</ref>

<math> \sigma \left( \rho\right) = \sqrt{\left( \frac{\sigma \left(p \right)}{p} \right) ^{2} +\left( \frac{\sigma \left(T \right)}{T} \right)^{2}} </math>

Assuming positive and negative error bar is the same in magnitude. The error bar stimulated density are too small (~<math> 10^{-2}</math>)to be shown on the graph.

Here is the plot of density as a function of temperature in two different (controlled) pressure, under simulation and ideal gas approximation.

[[File:Density against temp revised.tif|none|800x500 px|left|Density against temp revised]]

The stimulated density is lower than the calculated. This phenomenon can be explained by the assumption of gas molecules do not interact with each other in the ideal gas approximation model. In a fixed volume, more gas molecules can be found as the molecules do not repel with each other, hence density is higher.

The plot shows

(1)''Density is inversely proportional to temperature.''
''Negative gradients were found in the density against temperature graph.''

When temperature increases, the average kinetic energy increases, molecules are more energetic, and moving at a higher speed. In a fixed volume, fewer molecules are found (i.e. density decreases).

(2)''Discrepancy in density increases, as pressure increases.''
''Densities were observed to be higher at <math> p^{*} = 3.0 </math>than <math> p^{*} = 2.0 </math>(while other variables kept constant). The discrepancy between the simulation and the ideal gas model (i.e. less exact) increases when <math> p^{*} </math> increases from 2.0 to 3.0.''

When pressure increases, molecules are forced closer together and interact stronger. Consequently, volume at a given temperature is decreased, as the number of atoms N is fixed in the simulation (<math> pV = NT = const </math>). Density (<math>\rho =\frac{N}{V}</math>)at the given temperature increases, as the volume increases. There will be more interaction between molecules which has been neglected in the ideal gas assumption. Hence, the ideal gas model starts to break down, larger discrepancy will be observed.

== Calculating heat capacities using statistical physics ==

===Task 1===

'''As in the last section, you need to run simulations at ten phase points. In this section, we will be in density-temperature <math>\left(\rho^*, T^*\right)</math> phase space, rather than pressure-temperature phase space. The two densities required at <math>0.2</math> and <math>0.8</math>, and the temperature range is <math>2.0, 2.2, 2.4, 2.6, 2.8</math>. Plot <math>C_V/V</math> as a function of temperature, where <math>V</math> is the volume of the simulation cell, for both of your densities (on the same graph). Is the trend the one you would expect? Attach an example of one of your input scripts to your report.'''

Input script for <math> \rho^{*} = 0.2 </math> and <math> T^{*} =2.0 </math>,

<pre>
### DEFINE SIMULATION BOX GEOMETRY ###
lattice sc 0.2
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

### SPECIFY THE REQUIRED THERMODYNAMIC STATE ###
variable T equal 2.0
variable timestep equal 0.0025

### 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 10000
reset_timestep 0

### MEASURE SYSTEM STATE ###
unfix nvt
fix nve all nve
thermo_style custom step etotal temp press density atoms vol
variable dens equal density
variable dens2 equal density*density
variable temp equal temp
variable temp2 equal temp*temp
variable press equal press
variable press2 equal press*press
variable N equal atoms
variable N2 equal atoms*atoms
variable etotal equal etotal
variable etotal2 equal etotal*etotal
variable volume equal vol
fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_etotal v_dens2 v_temp2 v_press2 v_etotal2
run 100000

variable avedens equal f_aves[1]
variable avetemp equal f_aves[2]
variable avepress equal f_aves[3]
variable aveetotal equal f_aves[4]
variable aveetotal2 equal f_aves[8]
variable heatcapacity equal ${N2}*(${aveetotal2}-${aveetotal}*${aveetotal})/(${avetemp}*${avetemp})

print "Averages"
print "--------"
print "Density: ${avedens}"
print "Temperature: ${avetemp}"
print "Pressure: ${avepress}"
print "total energy: ${aveetotal}"
print "total energy^2: ${aveetotal2}"
print "heat capacity: ${heatcapacity}"
print "number of atoms: ${N2}"
print "Volume: ${volume}"
</pre>

[[File:Heatcappervolagainsttemp modified.tif|none|800x500 px|left|Heatcappervolagainsttemp_modified]]

The heat capacity at constant volume per unit volume against time plot is shown above.

The internal energy <math>U </math> varies with temperature and volume.

Heat capacity at constant volume <math>C_v </math> is defined as heat energy per unit temperature.

<math>C_V = \frac{\partial E}{\partial T} </math>

The volume is fixed, and we are interested in variation with the temperature. Since heat capacity is proportional to volume, we normalise the system to which is invariant with the size for comparison

''Higher <math>\frac{C_V}{V} </math> at specific temperature was observed in <math>\rho^{*} = 0.8 </math> than <math>\rho^{*} = 0.2 </math>.''

This can be explained molecules are closer at higher density. Higher potential energy is experienced as resulted from stronger interaction between molecules. Hence, the internal energy increases, consequently the thermal energy. Eventually, the <math>C_V </math> increases.

Another observation which out of expectation was also found.

[[File:Shottky Anomaly.png|thumb|250x150px|right|Schottky Anomaly]]

''Lower <math>\frac {C_V}{V} </math>, as temperature increases.''

According to the Schottky anomaly <ref name="Schottky anomaly"> The Schottky Anomaly, https://commons.wikimedia.org/wiki/File:Shottky_Anomaly.png#/media/File:Shottky_Anomaly.png, (accessed February, 2016)</ref>, which the ensemble in simulation is also NVT, at low temperature, there are more vacancies in higher energy states. Molecules are more likely to be populated from lower energy states to high energy states. Hence, heat capacity is expected to increase.

At high temperature (all simulations were run above
<math> T^{*)=2.0 = 2 \times 120 K =240 K </math>,
no quantum effect can be observed.

== Structural properties and the radial distribution function ==

===Task 1===

'''perform simulations of the Lennard-Jones system in the three phases. When each is complete, download the trajectory and 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, and attach a copy to your report. Discuss qualitatively the differences between the three RDFs, and what this tells you about the structure of the system in each phase. 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?'''

The simulations were run under these conditions,

solid <math> \rho = 1.2, T*= 0.8, \text {fcc}</math>;
liquid <math> \rho = 0.8, T*= 1.2 </math>; and
gas <math> \rho = 0.1 T*= 1.2</math>

Radial distribution function <math>g(r)</math> describes the variation of density of particles based on the distance from the reference point.

<gallery widths=500px heights=300px>
File:G(r) plot.tif
File:Int g(r).tif
File:Int g(r) zoom in.tif
</gallery>

In solid, atoms are tightly packed. Atoms are usually arranged in a regular manner. The peaks occur at the respective distance on the g(r) plot show the position of atoms. As we can see from the graph, more peaks are found at shorter distance. Solid has the highest density among three states. <math> g(r) </math> eventually fluctuates between 1.

In liquid, atoms are close together, but further apart than atoms in solids. Atoms are arranged irregularly. Less peaks are observed while the magnitude of the peak is lower than solid. Density of liquid is lower than solid. The peaks are located at distance shorter than 4 unit distance. <math> g(r) </math> eventually approaches 1.

In gas, atoms are very far apart. Only one peak at less than 2 unit distance on <math> g(r) </math> plot, is observed. <math> g(r) </math> approaches 1 after 2 unit distance. The <math> \int g(r) </math> is almost 0. This indicates there is no predictable arrangement of the system.

In the solid case, face-centred cubic structure was adopted. There are three major types of shortest distance in the unit cell. The respective distance can be found from the peaks on the g(r) plot.

<math> d_1 = 1.025 </math>

<math> d_2 = 1.475 </math>

<math> d_3 = 1.825 </math>

Let <math> y </math> be the length of the side of unit cell.

<math> a = b = c = y (\text{cubic}) </math>

<math> y^{2} + y^{2} = \left( 2 d_1 \right)^{2} \text{(Pythagoras's theorem)}</math>

<math> 2 y^{2} = 4 d_1^{2}</math>

<math> y^{2} = 2 \left( 1.025 \right)^{2}</math>

<math> y = \sqrt {2 \left( 1.025 \right)^{2}}</math>

<math> d_1 = 1.4496 \text{(5 s.f.)} </math>

<math> y = d_2 </math>

<math> = 1.475 </math>

<math> \left[ \left( \frac {y}{2} \right)^{2} + y^{2} \right] + \left( \frac {y}{2} \right)^{2} = d_3^{2} \text{ (Pythagoras's theorem)}</math>

<math> \frac {y^{2}}{4} + y^{2} + \frac {y^{2}}{4} = d_3 ^{2} </math>

<math> \frac {3 y^{2}}{2} + y^{2} = d_3 ^{2} </math>

<math> y = \sqrt {\frac {2 d_3 ^{2}}{3}} </math>

<math> y = \sqrt {\frac {2 1.825^{2}}{3}} </math>

<math> y = 1.4901 \text{(5 s.f.)} </math>

<math> \text{ average lattice spacing} = \frac {1.4496 + 1.475 + 1.4901}{3} </math>

<math> = 1.47 \text{(3 s.f.)} </math>

<math>\int g(r)\mathrm{d}r</math> gives the number of atoms per unit area (i.e. coordination number).

{| class="wikitable"
|-
! peak
! 1st
! 2nd
! 3rd
|-
! <math>\int g(r)\mathrm{d}r</math>
| 12.1 || 18.2||42.2|
|-
! coordination number
| 12
| 18-12= 6
| 42-18= 16
|}

== Dynamical properties and the diffusion coefficient ==

===Task 1===

''' make a plot for each of your simulations (solid, liquid, and gas), showing the mean squared displacement (the "total" MSD) as a function of timestep. Are these as you would expect? Estimate <math>D</math> 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.'''

The conditions for solid, liquid and gas simulations is the same as the previous part.

<gallery widths=450px heights=350px>
File:Solid modified.tif
File:Liquid modified.tif
File:Vapour modified.tif
</gallery>

The diffusion coefficient <math> D </math> can be defined from the mean squared displacement.

<math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>

<math>\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t} = 6 D</math>

Since the gradient of the "total" MSD against timestep gives <math> 6D </math>, a straight line graph is expected. Hence, <math> D </math> can be estimated from the slope of linear fitting. <math> D </math> is a measure of the rate of molecular motion<ref name = "Dynamic light scattering"> P.A.Atkins, ''Physical Chemistry'', Oxford University Press, OXFORD, 8th edn., 2006.</ref>. The more mobile the atoms in system, higher <math> D = \frac{1}{6} \text{grad.}</math> (larger slope in "total" MSD plot) should be observed.

{|class= "wikitable"
|-
!
! Gradient (simulation)
! D (simulation)
! Gradient (1 million atoms)
! D (1 million atoms)
|-
!solid
| 2.53E-8
| 4.22E-9
| 5.27E-8
| 8.78E-9
|-
!liquid
| 1.02E-3
| 1.70E-4
| 1.05E-3
| 1.75E-4
|-
!vapour
| 1.82E-2
| 3.03E-3
| 3.05E-2
| 5.08E-3
|}

The diffusion coefficient <math> D </math> increases as the the density of system increases (solid --> liquid --> gas). When the number of atoms in the simulation increases (simulation --> 10000 atoms vs 1 million atoms),<math> D </math> also increases. In addition, <math> D </math> increases more significantly as the density decreases (i.e. atoms tend to be in more random motion), esp from liquid to gas, <math> D </math> increases in the order of magnitude from <math> 10^{-9} </math> to <math> 10^{-4}</math>.

===Task 2===

''' In the theoretical section at the beginning, the equation for the evolution of the position of a 1D harmonic oscillator as a function of time was given. Using this, evaluate the normalised velocity autocorrelation function for a 1D harmonic oscillator (it is analytic!):

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

'''Be sure to show your working in your writeup. On the same graph, with x range 0 to 500, plot <math>C\left(\tau\right)</math> with <math>\omega = 1/2\pi</math> and the VACFs from your liquid and solid simulations. What do the minima in the VACFs for the liquid and solid system represent? Discuss the origin of the differences between the liquid and solid VACFs. The harmonic oscillator VACF is very different to the Lennard Jones solid and liquid. Why is this? Attach a copy of your plot to your writeup.'''

For 1D harmonic oscillator,

The displacement can be defined as

<math> x\left(t\right) = A\cos\left(\omega \mathit{t} \right)</math>

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

<math> = \frac {\mathrm{d}A\cos\left(\omega \mathit{t} \right)}{\mathrm{d}\mathit{t}}</math>

<math> = A \omega -\sin\left(\omega\mathit{t}\right)</math>

<math> = -A \omega \sin\left(\omega\mathit{t}\right)</math>

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

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

using <math> \sin(a)\sin(b) = \frac{\left( cos (a-b) - cos (a+b) \right)}{2} </math>

<math>= \frac{\int_{-\infty}^{\infty} \frac{1}{2} \times \left( \cos(\omega\mathit{t}-\omega\mathit{t}-\omega\tau)- \cos(\omega\mathit{t}+\omega\mathit{t}+\omega\tau) \right) \mathrm{d}t}{\int_{-\infty}^{\infty} \frac{1}{2} \times ( 1- \cos(2\omega\mathit{t})\mathrm{d}t}</math>

<math>= \frac{\int_{-\infty}^{\infty} \left( \cos(-\omega\tau)- \cos(2 \omega\mathit{t}+\omega\tau) \right) \mathrm{d}t}{\int_{-\infty}^{\infty} ( 1- \cos(2\omega\mathit{t})\mathrm{d}t}</math>

<math>= \lim_{R \to {\infty}} \frac{\int_{-\infty}^{\infty} \left( \cos(-\omega\tau)- \cos(2 \omega\mathit{t}+\omega\tau) \right) \mathrm{d}t}{\int_{-\infty}^{\infty} ( 1- \cos(2\omega\mathit{t})\mathrm{d}t}</math>

<math>= \lim_{R \to {\infty}} \frac { 2R\cos(-\omega\tau)- \frac {1}{2 \omega} \left( \sin (2 \omega R+\omega\tau)- \sin(-2 \omega R+ \omega\tau) \right)} {2R - \frac {1}{2 \omega} 2\sin (2\omega R)}</math>

<math>= \lim_{R \to {\infty}} \frac { 2R\cos(-\omega\tau)} {2R}</math>

<math>= cos(\omega \tau) </math>

[[File:VACF.tif|none|800x500px|left|Velocity autocorrelation function]]

The plot shows the velocity autocorrelation function (VACF) as a function of temperature in 1D harmonic oscillator, liquid and solid simulations.
The global minima occur at -1.30943 and -0.15104 for solid and liquid respectively. Since velocity is a vector quantity, the negative value indicates the change in direction.

The velocities (solid and liquid) generated from the simulation do not retain to that speed after reaching minima. After collision, the atom loses its initial velocity and moves round in random Brownian motion, which eventually will decay to zero. In a denser system (solids), the velocity correlation decays quicker to zero. System with fixed atomic positions (solid) diffuses faster after the initial simulation. This is analogous to a very damped oscillation. In liquid, atoms experience destructive intermolecular forces from the Lennard Jones potential, which takes longer time to decay. The harmonic oscillator, a single atomic system, never decays. It undergoes the simple harmonic motion with constant magnitude of amplitude and frequency. The change in VACF of the harmonic oscillator is the natural change in speed in the interchange between kinetic and potential energy.

===Task 3===

'''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 <math>D</math> in each case. You should make a plot of the running integral in each case. Are they as you expect? Repeat this procedure for the VACF data that you were given from the one million atom simulations. What do you think is the largest source of error in your estimates of D from the VACF?'''

The trapezium rule states the area under the curve can be estimated from the summing the trapezium stripes.

For a function <math> y = f(x) </math>,

<math> \int_{0}^{x_n} f \left(x\right) \, \mathrm{d}x \approx (b-a) \left[\frac{f(a) + f(b)}{2} \right]</math>

<gallery widths=450px heights=350px>
File:Simulation int VACF.tif
File:1m atoms int VACF.tif
</gallery>

Since diffusion coefficient can be described 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>

<math> = \frac{1}{3}\int_0^\infty C \left( \tau \right) \mathrm{d}\tau </math>

<math> \approx \frac{1}{3}\times \text {integral of VACF from trapezium rule} </math>

{| class="wikitable"
|-
! abs(int VACF)
! solid
! liquid
! gas
|-
! simulation
| 2.85E-03
| 2.56E-02
| 3.84E-02
|-
! 1m atoms
| 8.71E-05
| 2.81E-03
| 1.33E-01
|}

{| class="wikitable"
|-
! D
! solid
! liquid
! gas
|-
! simulation
| 9.50E-04
| 8.52E-03
| 1.28E-02
|-
! 1m atoms
| 2.90E-05
| 9.37E-04
| 4.44E-02
|}

==References==