View source for Rep:Mod:Yr3LiquidSimJmsD

=Jack Dawson - Yr3 Comp: Liquid Sim.=

<span style="color:red">Not a bad attempt, although slightly uncomplete. You've attempted to give physical explanations for observed phenomena and generally I think you've understood the material to a good standard. Mark: 63/100</span>

== 1 Theory ==

==== 1.1 Task ====
'''''Compltete the Classical Harmonic Osciallator Spreadsheet'''''

[[File:JmsD_HO.xls|thumb|1250px|center|]]


==== 1.2 Task ====
'''''For δt=0.01, estimate the positions of maximum error and fit accordingly.'''''

[[File:JmsD Max Error Graph.png|thumb|800px|center|'''Figure 1.1:''' ''A plot of the error in VV algorithm data'']]

Figure 1.1 shows how the maximum error increases linearly. This can be rationalised by the fact that each cycle of simulation follows on from the last, so the error of the new cycle is combined to the uncertainty in a previous calculation for an older timestep and so on. <span style="color:red">Tick</span>

==== 1.3 Task ====
'''''Experiment to find the timestep range to ensure the total energy does not change by more than 1% over the course of the simulation. Why is this important?'''''

In simulations such as this, the physical model is considered an isolated system. We are modeling the changes in velocity and thus kinetic energy changes of induvidual particles/atoms over a total timeframe at constant timesteps. As such, the particles are moving around, interacting, transfering energy between each other. However, unless we are changing variables such as Temperature, Number of particles, Volume (density) etc. during a simulation, the overall energy should remain constant. Therefore it is important to monitor the total energy to ensure our model is fitting with our physical understanding.

From experimentation, it can be seen that to limit this total energy change over a simulation, the time step should really be smaller than 0.3. It can be seen in Figure 1.2 how the peaks and thus total energy begin to fluctuate significantly more for 0.3 and above.<span style="color:red">It's more like 0.2 but right idea</span>

[[File:JmsD_En_err_w_timestep.png|thumb|500px|center|'''Figure 1.2:''' ''Energy vs Time: a) δt=0.01, b) δt=0.03, c) δt=0.03'']]

==== 1.4 Task ====
'''''Maths based upon a single Lennard-Jones interaction'''''
 
[[File:JmsD_Maths_1.png|thumb|800px|center|'''Figure 1.3:''' ''Answers to the LJ based maths task'']] <span style="color:red">Tick</span>

==== 1.5 Task ====
'''''Estimate the number of water molecules in 1 mL under standard conditions.'''''
<center>Density of water = 1 g/ml</center>
<center>Molecular weight of wate = 18 g/mol</center>
<center>Moles in 1 ml, n = 1/18</center>
<center>Number of molecules = n x N<sub>A</sub> = 1/18 x 6.022x10<sup>23</sup></center>
<center>3.35x10<sup>22</sup></center><span style="color:red">Tick</span>


'''''Estimate the volume of 10000 water molecules under standard conditions.'''''
<center>10000/N<sub>A</sub> = moles of water, n</center>
<center>Mass of 10000 water molecules, m = n x 18</center>
<center>The volume, V = m / 1 g/ml </center>
<center>2.99x10<sup>-19</sup> </center><span style="color:red">Tick</span>

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

The particle should end up in position (0.2,0.1,0.7) after the application of periodic boundary conditions. 

This is because, if the number of particles in the box remains constant, if a particle excit the box on one face, it will return from the opposite. In this case, if the particle starts at (0.5,0.5,0.5) in the cubic box ((0,0,0):(1,1,1)) and moves along the vector (0.7,0.6,0.2), it will reappear from the other face, positioned now at (0.2,0.1,0.7).<span style="color:red">Tick</span>

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

<center>r = 1.088 x10<sup>-10</sup> m</center> <span style="color:red">Should be E-9 but again right idea</span>
<center>T = 180 K</center><span style="color:red">Tick</span>
<center><math>\epsilon</math> = 0.014 kJmol<sup>-1</sup></center><span style="color:red">not quite right, 0.997 kJ mol<sup>-1</sup></span>

== 2 Equilibriation ==
==== 2.1 Task ====
'''''Why do random starting coords cause problems in simulations'''''
If particles are placed randomly, there is chance for multiple particles to be placed in the same positions, or at least close too each other. It two particles are close or overlap, there would be very strong resulting interactions. These large forces would caus the model to deviate from its classical behaviour, and will affect properties from energy, to timestep. <span style="color:red">It doesn't deviate from classical behviour just the energies are so large that the corresponding timestep should be unreasonably small for a simulation without "blow-up"</span>

Therefore it is wise to start a liquid or gas simulation as a 'melt' from particles on lattice points (which can allow random positioning to some degree), e.g. Simple cubic. This ensures particles will begin a suitable distance away from each other.

==== 2.2 Task ====
''''' Prove that lattice spacing of 1.07722 gives a number density of 0.8'''''
<center> If lattice spacing is 1.07722, Volume of one cube = 1.07722<sup>3</sup></center>
<center> If number density is Number / Volume, and there is one particle per cube, number density = 1 / Volume</center>
<center> 1 / 1.07722<sup>3</sup> = 0.79999 </center><span style="color:red">Tick</span>

'''''If instead FCC, and number density of 1.2, what is the length of the cubic uni cell?'''''
<center> First calculate Volume of the cell: V = Number / number density </center>
<center> The Number of particles in a FCC unit cell are 4, so V = 4 / 1.2 </center>
<center> The cube root of the volume will give side length...</center>
<center> (4/1.2)<sup>1/3</sup> = 1.4938 </center><span style="color:red">Tick</span>

==== 2.3 Task ====
'''''For the FCC latice, how many atoms would be created by create_atoms'''''

Assuming no other parameters  (i.e. the box size [Below]) have changed...
<pre>
region box block 0 10 0 10 0 10
create_box 1 box
</pre>
There should be 4000 atoms created.

If there are 1000 atoms in a 10x10x10 box for a sc lattice that has one atom per unit cell, a lattice such as FCC that has 4 atoms per unit cell should yield 4000 atoms. <span style="color:red">Tick</span>

==== 2.4 Task ====
'''''Find and explain the purpose of:'''''
<pre>
mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0
</pre>

*''mass'' - Sets the mass for the types of atoms, argument(1) is the type, argument(2) is the mass.
*''pair_style'' - Calculates the pair interactions (pair potentials defined as interactions between atoms within a certain distance), argument(1) "lj/cut" means this line refers to the global cutoff for all Lennard-Jones interactions, argument(2) defines the distance.
*''pair_coeff'' - Provides the interaction coefficients for the interactions between atoms, arguments(1 & 2) refer to the attoms that interact (* means all types 1->N), and arguments (3 & 4) are those coefficients.<span style="color:red">Tick</span>

==== 2.5 Task ====
'''''Given that we are specifying x<sub>i</sub> and v<sub>i</sub>, which integration algorithm is to be used?'''''
<center> The Velocity-Verlet Algorithm </center>


==== 2.6 Task ====
<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 purpose of using ''${timestep}'', or for any other variable for that matter, is that it allows for that property in the simulation to be changed just once in the input file, as opposed to the file being manually altered on every line that variable is called.<span style="color:red">Tick</span>

==== 2.6 Task ====
'''''Plot Energy, Temperature and Pressure vs Time for timestep = 0.001'''''

[[File:JmsD_0.001all.png|thumb|600px|center|'''Figure 2.1:''' ''The data sets from the δt=0.001 initial simulation'']]

[[File:JmsD_IntroEnergy.png|thumb|600px|center|'''Figure 2.2:''' ''All initial timestep data, Total Energy vs Time'']]

[[File:JmsD_IntroEnergyZoomed.png|thumb|600px|center|'''Figure 2.3:''' ''Expansion of Figure 2.2 in region Time(0:35)'']]

'''''Plot all the inital simulations as Energy vs Time, discus the best timestep to progress with.'''''

The timestep chosen for ongoing calculations was 0.0025. The reason being that it was one of only two that provided data showing convergence within the full simulation time (Energy vs Time graph). Of these two timesteps - 0.0025 and 0.001 - the 0.0025 is the better choice to progres with due to it's slight advantage in computational efficiency. The larger the time step, the fewer data points collected.

[A note should be made about te 0.15 timestep - This is quite clearly a bad choice because instead of converging it diverged...] <span style="color:red">Tick but why do we see this trend of tot eners vs timesteps?</span>

== 3 Simmulations under specific conditions ==
==== 3.1 Task ====
'''''Choose 5 temperatures (above the critical temperature, T* = 1.5), and two reasonable pressures in Lennard-Jones units.'''''

'''''Run these simulations to plot the 'Equations of State' [Below].'''''

The 10 phase points used inn this task were determined by the 5 temperatures (T): 1.6, 1.8, 2.0, 2.2, 2.4 and the pressures (p) 2.4, 2.8. The timestep (δt) chosen was again 0.0025


Example input file: [[File:JmsD_Npt_1.6_2.4.in|thumb|1250px|center|]]

Example output file: [[File:JmsD_Npt_1.6_2.4.log|thumb|1250px|center|]]

==== 3.2 Task ====

The temperature ideally needs to be maintained as close to the target temperature, <math>\mathfrak{T}</math>, as possible. This is done by multplying every velocity, by a constant factor, <math>\gamma</math>.

The 'target temperature' is as stipulated in the input script

The instantaneous temperature, <math>T</math>, will fluctuate  due to approximations made in the algorithm. 
*If <math>T > \mathfrak{T}</math>, the E<sub>K</sub> is too high and thus <math>\gamma</math> must be less than 1.
*If, by contrast, <math>T < \mathfrak{T}</math>, <math>\gamma</math> must be greater than 1.


The 'Equipartition theorem' states that on average, every degree of freedom in a system at equilibrium will have K<sub>B</sub>T/2 of energy. This sytem has 3 degrees of freedom and as such:

<math>\frac{1}{2}\sum_i m_i \mathbf{v}_i^2 = \frac{3}{2} N k_B T</math>

Expressing in terms of the adjustment to target temperature as described above:

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


'''''Solve the equations for <math>\gamma</math>''''':


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

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

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

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

==== 3.3 Task ====
'''''Find out the importance of the three numbers 100 1000 100000.'''''

* 100 is used as the ''Nevery'' argument - This is the number of timesteps between inputing of values.
* 1000 is used as the ''Nrepeat'' argument - This is the number of times to use the input values to calculate averages.
* 100000 is used as the ''Nfreq'' argument - The number of timesteps between calculating averages.

<pre>
fix aves all ave/time 100 1000 100000 {...VARIABLES...}
run 100000
</pre>

'''''How much time will this simulate?'''''

* ''run 100000'' shows the number of timesteps that will be run in the simulation.

This shows that averages in this simulation are calculated only at the end of the simulations, as ''Nfreq''= ''N'' (run N ...).

==== 3.4 Task ====
'''''Using the output of the simulations in 3.1, plot density as a unction of temperature, for both pressures - (Incl. Error bars + Ideal Gas). Is the simulated density higher or lower than Ideal? Why? Finaly, does the discrepancy increase or decrease w. pressure?'''''

[[File:JmsD_NpT1.png|thumb|500px|center|'''Figure 3.1:''' ''A plot of Average Density vs Average Temperature for two different pressures'']]

[[File:JmsD_NpT_wideal.png|thumb|500px|center|'''Figure 3.2:''' ''A plot of Average Density vs Average Temperature for two different pressures and their relation to data with the same parameters calculated by the Ideal Gas Law'']]

Theoretical results , based on the ideal gas law, give a higher density than the simulated result. (It can also be added that the discrepancy between the ideal and simulated systems increases with pressure.) It can be reasoned that the results provided by the Ideal Gas Law are higher due to the fact that particle-particle interactions are not consided.This means the repulsion that would be experienced by the particles moving closer to each other are not expressed, as such, particles the distance between particles can decrease without energyetic effects, thus increasing the density. <span style="color:red">Also Charles' law!</span>

== 4 Calculation of heat capacities using statistical physics ==
==== 4.1 Task ====
'''''Run simulations and plot C<sub>v</sub>/V as a function of temperature. Is this as expected?'''''

Example input file: [[File:JmsD_Cv_2.0_0.2.in|thumb|1250px|center|]]

Example output file: [[File:JmsD_Cv_2.0_0.2.log|thumb|1250px|center|]]

For this section, 5 simulations were run for each of the two densities: 0.2, 0.8. These simulations corresponded to each of the respective temperatures: 2.0, 2.2, 2.4, 2.6, 2.8 (in reduced units). The desired results of these simulations were the calculated values of the Heat Capacity at constant volume, C<sub>v</sub>. The results can be seen in Figure 4.1 as a funtion of the average temperatures.

[[File:JmsD Cv graph1.png|thumb|500px|center|'''Figure 4.1:''' ''A plot comparing Heat Capacity/Volume vs Average Temperature for two different densities'']]

These graphs do indeed follow the trend expected, that is similar to the trend seen in a crystalline solid for example. The Heat Capacity is proportional to 1/T<sup>2</sup> [As seen in the Equations below] and as such a negative trend should be observed. <span style="color:red"> Does it vary with the reciprocal of the square? Importantly, the variance in energy varies with temperature and therefore, we cannot take this argument. Think more either internal energies or density of states..</span>

<math>C_V = \frac{\partial E}{\partial T} = N^2\frac{\left\langle E^2\right\rangle - \left\langle E\right\rangle^2}{k_B T^2}</math>

With reagrds to density, C<sub>v</sub>/V increases as the density increases - which is again as expected - because heat capacity is dependant on the amount of material per volume, and thus a denser material (with the same volume) has more mass or moles per unit voulume.

== 5 Radial Distribution Functions ==
==== 5.1 Task ====
''''' From simulations of the Lennard-Jones system in the three phases, Compute <math>g(r)</math> and <math>\int g(r)\mathrm{d}r</math> in VMD. Discuss the differences in the 3 RDFs and thus what this means about the structure in each phase.'''''

'''''For the solid phase, illustrate which lattics sites the first 3 peaks correspond to. What is the Lattice spacing? What is the Coordination number for each of the first three peaks?'''''

The initial parameters for the simulations, in order to specify the correct phases were:

<center> Temperature: (For all phases) = 1.2 </center>
<center> Densities: (Solid) = 1.50, (Liquid) = 0.85, (Vapour) = 0.01 </center>
<center> Finally, to provide the initial positions of atoms, the simple cubic model was used for Liquid and Vapour, but for Solid FCC was used.</center>

These parameters were calculated by examining the Lennar-Jones phase diagram found in literature. [http://journals.aps.org/pr/abstract/10.1103/PhysRev.184.151 <sup>[2] </sup>]


[[File:JmsD_g(r).png|thumb|500px|center|'''Figure 5.1:''' ''The Radial Distribution Functions of the three states'']]

[[File:JmsD_Int(g(r))dr.png|thumb|500px|center|'''Figure 5.2:''' ''The Integrals of the RDFs of the three states'']]

[[File:JmsD_Solid3(g(r)).png|thumb|500px|center|'''Figure 5.3:''' ''The RDF of the solid, paying particular attention to the first 3 peaks'']]


As it can be seen in Figure 5.1, the RDFs are very different for each phase. The area under the curves is the probability of finding another particle at that distance from the particle of reference, thus peaks are distances away from the reference that there is increased probability of encountering another.<span style="color:red">Tick for a normalised rdf</span>

For example, the Vapour phase [as shown in blue in Figure 5.1] is largely simple. Whereas, the RDF of the liquid [purple] has a series of fairly ordered peaks of decresing height. This can be explained by the ease of motion of these two phases. Vapour particles are very spread out in space and have a lot more motion than liquids which are held together and inplace by various intermolecular/atomic bonds. Because of this by moving away fro a reference particle in the vapour phase, there is a very low probability of encountaring another particle. Liquids are a lot more ordered than gases and coordination spheres are shown by the fairly smooth peaks as r increases. There is reasonable probability of finsing another particle in the first coordination sphere, but this decreases as you move away further.<span style="color:red">Tick</span>

Solids, on the other hand, have significantly greater order. The particles are generally held in a rather rigid lattice and their positions can be considered fairly fixed. In Figure 5.1 you can clearly see sharp peaks on the RDF as you move away from the reference, corresponding to neighbouring atoms. These come at regular distance intervals, which are the length of a unit cell away - This RDF in particular representing that of an FCC lattice. The area of these peaks does deacreas with distance somewhat as there is some movement, such as vibrations (especially if the particle's vibrations about the lattice point are coupled with vibrations of the lattice as a whole), that cause those further out to be proportionaly missalligned from their perfect lattice position.<span style="color:red">Tick, good explanation</span>

== 6 Dynamical properties and the diffusion coefficient ==
==== 6.1 Task ====
'''''Run simlations for the three phases, to calculate the mean squared displacement and velocity autocorrelation functions.'''''

Example input file: [[File:JmsD.in|thumb|1250px|center|]]

Example output file: [[File:JmsD.log|thumb|1250px|center|]]

==== 6.2 Task ====
'''''Measure D from it's connection to the MSD, for simulations and data for approx. one million atom simulations'''''

<span style="color:red">These diff coeffs seem rather small.. Could be use of timestep?</span>

The equation for the Diffusion Coefficient, <math>D = \frac{1}{6}\frac{\partial\left\langle r^2\left(t\right)\right\rangle}{\partial t}</math>.

However, <math>D</math> can be much more simply estimated by fitting the linear section of a Total MSD vs Timestep graph and dividing by 6.
This gives the following values for D (in reduced units) for the smaller simulations carried out by myself:
<center> D<sub>Solid</sub> = 1.31x10<sup>-10</sup> </center>
<center> D<sub>Liquid</sub> = 1.33<sup>-4</sup> </center>
<center> D<sub>Vapour</sub> = 1.55<sup>-2</sup> </center>
<center>[The graphs for these calculations can be seen in Figures 6.1-6.3]</center>

[[File:JmsD_my_sol_D.png|thumb|500px|center|'''Figure 6.1:''' ''The Total MSD vs timestep (with fitting) for my Solid simulation'']]
[[File:JmsD_my_liq_D.png|thumb|500px|center|'''Figure 6.2:''' ''The Total MSD vs timestep (with fitting) for my Liquid simulation'']]
[[File:JmsD_my_vap_D.png|thumb|500px|center|'''Figure 6.3:''' ''The Total MSD vs timestep (with fitting) for my Vapour simulation'']]

In comparison below are the values for D (in reduced units) for the larger (approx. 1 million atom) data supplied:
<center> D<sub>Solid</sub> = 1.71x10<sup>-8</sup> </center>
<center> D<sub>Liquid</sub> = 1.75<sup>-4</sup> </center>
<center> D<sub>Vapour</sub> = 6.28<sup>-3</sup> </center>
<center>[The graphs for these calculations can be seen in Figures 6.4-6.6]</center>

[[File:JmsD_mil_sol_D.png|thumb|500px|center|'''Figure 6.4:''' ''The Total MSD vs timestep (with fitting) for supplied Solid MSD data'']]
[[File:JmsD_mil_liq_D.png|thumb|500px|center|'''Figure 6.5:''' ''The Total MSD vs timestep (with fitting) for supplied Liquid MSD data'']]
[[File:JmsDil_vap_D.png|thumb|500px|center|'''Figure 6.6:''' ''The Total MSD vs timestep (with fitting) for supplied Vapour MSD data'']] <span style="color:red">Tick</span>

== 7 References ==
 1.  LAMMPS Manual, http://lammps.sandia.gov/doc/Section_commands.html#lammps-input-script , accessed November 2016
 2.  JP Hansen, L Verlet, Phys Rev, 184, 1969, http://journals.aps.org/pr/abstract/10.1103/PhysRev.184.151