Rep:Mod:Yr3LiquidSimJmsD

Rep:Mod:Yr3LiquidSimJmsD

2017-02-24T09:26:46Z

Jmd114: /* 2 Equilibriation */

=Jack Dawson - Yr3 Comp: Liquid Sim.=

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

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

[[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'']]

==== 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 NA = 1/18 x 6.022x1023</center>
<center>3.35x1022</center>

'''''Estimate the volume of 10000 water molecules under standard conditions.'''''
<center>10000/NA = 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-19 </center>

==== 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).

==== 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-1? What is the reduced temperature T*=1.5 in real units? '''''

<center>r = 1.088 x10-10 m</center>
<center>T = 180 K</center>
<center><math>\epsilon</math> = 0.014 kJmol-1</center>

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

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.077223</center>
<center> If number density is Number / Volume, and there is one particle per cube, number density = 1 / Volume</center>
<center> 1 / 1.077223 = 0.79999 </center>

'''''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)1/3 = 1.4938 </center>

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

==== 2.3 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.

== 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 EK 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 KBT/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>

==== 3.3 Task ====
'''''Find out the importance of the three numbers 100 1000 100000.'''''
'''''How often are values sampled for the average?'''''
'''''How many measurments contribute to the average?'''''
<pre>
fix aves all ave/time 100 1000 100000 {...VARIABLES...}
run 100000
</pre>
'''''How much time will this simulate?'''''

==== 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'']]

== 4 Calculation of heat capacities using statistical physics ==
==== 4.1 Task ====
'''''Run simulations and plot Cv/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, Cv. 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/T2 [As seen in the Equations below] and as such a negative trend should be observed.

<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, Cv/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 [2] ]

[[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'']]

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

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> DSolid = 1.31x10-10 </center>
<center> DLiquid = 1.33-4 </center>
<center> DVapour = 1.55-2 </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> DSolid = 1.71x10-8 </center>
<center> DLiquid = 1.75-4 </center>
<center> DVapour = 6.28-3 </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'']]

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

Rep:Mod:Yr3LiquidSimJmsD