View source for Talk:Mod:fpp1994

==Introduction to Molecular Dynamics==


[[File:IEIEIEIEIIEIEIEIE.PNG]]
[[File:FrankTASK1.xls]]



[[File:Lelelelele.PNG]][[File:Momomom.PNG]] 
[[File:FppTASK2.xls]]


[[File:Hadfdfhaha.PNG]][[File:Eroroor.PNG]]

[[File:FppTASK3.xls]] Here, it can be seen that when the time-step is 0.2 the total energy fluctuates by 1%. It is important to model the energy of a system you are modelling because in some cases it can be used to calculate thermodynamic quantities as energy is what will dictate how things behave.
'''<span style="color:red">NJ: It's more important to monitor the conservation of energy, this is an indicator of how well the numerical algorithm is modelling the system.</span>'''
'''<span style="color:red">NJ: Why do you think the error accumulates in the way that it does?</span>'''
----


'''Task 4'''

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

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

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

Thus <math> \sigma = r_0</math> ie <math>\sigma</math> is the distance when the potential is zero  

*At this seperation, the force is <math> \frac{-24\epsilon}{r_0}</math>

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

<math> F = 24\epsilon \left(\frac{-2}{r_0} + \frac{1}{r_0}\right) = -\frac{24\epsilon}{r_0}</math>

'''<span style="color:red">NJ: Be careful with the notation, <math>F = -\frac{dU}{dr}</math>, so this is out by a factor of -1.</span>'''

*The equilibrium seperation, <math>r = r_{eq}</math> occurs at the bottom of the well, so force is at a minimum.
'''<span style="color:red">NJ: The potential is at a minimum, so the force is zero.</span>'''

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

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

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

<math>r_{eq} = r_0\sqrt[6]{2} = 1.122 r_0</math>

*<math>\phi\left(r_{eq}\right) = 4\epsilon\left(\frac{r_0^{12}}{r_{eq}^{12}} - \frac{r_0^6}{r_{eq}^6}\right) </math> and from just above <math>r_{eq} = r_0\sqrt[6]{2}</math> 
'''<span style="color:red">NJ: Good</span>'''

So <math>\phi\left(r_{eq}\right) = 4\epsilon\left(\frac{r_0^{12}}{2^6r_{0}^{12}} - \frac{r_0^6}{2r_0^6}\right) = 4\epsilon\left(\frac{1}{2^6} - \frac{1}{2}\right) = 1.936\epsilon</math>
'''<span style="color:red">NJ: You should get a value of <math>\epsilon</math> here </span>'''


----

<math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r = 4\epsilon \left[- \frac{\sigma^{12}}{11r^{11}} + \frac{\sigma^6}{5r^5} \right]_{2\sigma}^\infty</math>

The expansion of the above bracket is <math> 4\epsilon \left( \left( - \frac{\sigma^{12}}{11 \left( \infty \right)^{11}} + \frac{\sigma^6}{5 \left( \infty \right)^{5}} \right) - \left( - \frac{\sigma^{12}}{11 \left( 2\sigma \right)^{11}} + \frac{\sigma^6}{5 \left( 2\sigma \right)^{5}} \right) \right) </math> but since <math> \frac{1}{\infty^{n}} = 0 </math> and <math> \epsilon = \sigma = 1 </math> this becomes <math> 4 \left( \frac{1}{11 \left( 2 \right)^{11}} - \frac{1}{5 \left( 2 \right)^{5}} \right) = -0.0248 </math>

Using the same expansion, but <math>2.5\sigma</math> and <math>3\sigma</math> as the lower limits gives <math>-0.00818</math> and <math>-0.00329</math> respectively.

'''<span style="color:red">NJ: Good</span>'''
----


'''Task 5:''' 

The density of water is <math>1000 \mathrm{ kg\ m}^{-3}</math> and so <math>1 \mathrm{mL}</math> would weigh <math>1 \mathrm{g}</math>. The molecular mass of water is <math> 18.0 \mathrm{g mol}^{-1}</math> and so the amount of molecules in <math>1 \mathrm{mL}</math> of water is <math>1/18 \times N_{A} = 3.35 \times 10^{22}</math>. Therefore <math>10000</math> molecules would occupy <math>\frac{10000}{3.35 \times 10^{22}} = 2.99 \times 10^{-18}\mathrm{mL} </math>.

'''<span style="color:red">NJ: Good</span>'''
----


'''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?. It would be at <math>\left(0.2, 0.1, 0.7\right)</math>


----

'''Task 7'''

*<math>r^* = \frac{r}{\sigma}</math> therefore <math> r^*\sigma = r = 0.34 \times 10^{-9} \times 3.2 = 1.08\times 10^{-9} \mathrm{m}</math>
'''<span style="color:red">NJ: 1.088nm</span>'''

*If <math>\frac{\epsilon}{k_B} = 120\mathrm{K}</math> then since <math>k_B = 1.38 \times 10^{-23}</math> the well depth in <math>\mathrm{kJmol^{-1}}</math> is <math>\epsilon = 120 \times 1.38 \times 10^{-23} \times 1000 = 1.66 \times 10^{-18} \mathrm{kJmol^{-1}}
'''<span style="color:red">NJ: This is just in kJ, not kJ per mole.</span>'''
</math>
* temperature <math>T^* = \frac{k_BT}{\epsilon}</math> therefore <math> T = \frac{T^{*}\epsilon}{k_B} = 1.5 \times 120 = 180 \mathrm{K} </math>

==Equilibration==

'''Task 1''' Giving atoms random starting positions in simulations can cause problems because if there are two atoms that are placed too close together (ie if they are within <math> r_0 </math> of each other) then their potential energies - and thus initial accelerations and velocities -will be extremely high and they will move through the sample with an unrealistically high speed and this will also disrupt other atoms.

----

'''Task 2''' 

If each the spacing lattice point is 1.07722 and the lattice is three dimensional then the number density of lattice point is given by <math> \frac{1}{spacing} = \frac{1}{1.07722^3} = 0.8 </math>

If a face-centred cubic lattice has a number density of 1.2 then because it has four atoms per cell its volume can be given by <math> V = \frac{number}{density} = \frac{4}{1.2} </math> and so the side length is just the cube root of this and so <math> l = \sqrt[3]{\frac{4}{1.2}} = 1.49 </math>

'''Task 3''' 

In the same lattice if the command "create_atoms" was used, 4000 atoms would be created as there are 4 atoms per unit cell.


----

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

The first command sets the relative mass of Atom 1 as 1.

In the second command "pair_style" is used to describe interactions between two particles and the "lj/cut" part tells LAMMPS that this interaction follows the Lennard-Jones potential and nothing more and that it cuts of when <math> r = 3 </math>.

"pair_coeff" tells LAMMPS what coefficients to use in the style defined above. In this case it is the Lennard-Jones potential, and so it is telling LAMMPS to set <math> \epsilon </math> and <math>\sigma</math> as 1.

----

'''Task 5''' Specifying velocity and position as starting conditions means the Velocity Verlet Algorithm should be used.

----

'''Task 6'''

Specifying variables rather than numerical values useful is because you may need to use that value at various points throughout the script. Using numerical values all over the script would become particularly annoying if you were running lots of simulations with the same variable multiple times as you would have to look through the whole script to find each value but with the dollar command you can just change it once. 

'''Task 7'''

[[File:ALL3GRAPH1.PNG]]

The system reaches equilibrium after a time of about 0.5
[[File:Total energy.PNG]]

5 simulations were run with timestops of 0.001, 0.01, 0.0025, 0.0075, 0.015. The worst choice is the one with the highest timestop, 0.015 because ideally a balance between resolution and the amount of time the system can be monitored is preferred but all of the timestops show the eventual equilibrium and so there is little benefit to a large value as instead it does not show the fact that equilibrium had to be reached.

==Running simulations under specific conditions==

'''Task 1''' 

10 phase points (two pressures and 5 temperatures) were selected, <math>T^* = 1.6,  1.8,  2.0,  2.2,  2.4 </math> and <math> p^* = 2.4,  2.8</math> and the timestep was 0.001

'''Task 2''' 

<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 = \gamma^2 \frac{1}{2}\sum_i m_i v_i^2\ = \frac{3}{2} N k_B \mathfrak{T}</math>



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

Cancel <math>\frac{1}{2}\sum_i m_i v_i^2</math> to give <math>\gamma^2 = \frac{\mathfrak{T}}{T}
</math>

'''Task 3'''

From the script

<pre> fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_dens2 v_temp2 v_press2 </pre>

From the LAMMPS manual

<pre> fix ID group-ID ave/time Nevery Nrepeat Nfreq value1 value2 ... keyword args ...</pre>

The numbers 100, 1000 and 100000 correspond to Nevery, Nrepeat and Nfreq respectively. These commands tell LAMMPS on which timesteps to calculate the averages of the properties that follow. Nrepeat tells LAMMPS how many averages to calculate, starting from the value given by Nfreq and working back in multiples of Nevery. In this case 1000 averages are taken in total, once every 100 timesteps up to 100000 which is as far as the script will run.

[[File:GRAPHS NPT 1.PNG]]

'''Task 4'''

The bottom two lines are the density as a function of temperature as calculated by LAMMPS whereas the top lines are calculated using the ideal gas law using the equation derived below where pressure was 2.4 and 2.8 and the temperature was the LAMMPS average. The simulated density is much lower than the one predicted by the gas law because the gas is not ideal as there are repulsive and attractive terms in the Lennard-Jones potential and ideal gases have no interaction between the particles. And as the pressure increases the gap between the simulated and calculated densities increases too because at higher pressures there are more interaction between the particles and so the ideal gas approximation gets worse.

<math>PV = Nk_BT</math>
 
<math>P = \frac{N}{V}k_BT </math>   and   <math>   \rho = \frac{N}{V}</math>


so <math>P = \rho k_BT      </math>     but on LAMMPS in lj style, <math> \epsilon, \sigma </math> and <math> k_B </math> are all equal to 1

so <math> \frac{P}{T} = \rho</math>

==Calculating heat capacities using statistical physics==

The script used to calculate the heat capacities is below

<pre>

### DEFINE SIMULATION BOX GEOMETRY ###
variable density equal  0.2
lattice sc ${density}
region box block 0 15 0 15 0 15
create_box 1 box
create_atoms 1 box
variable atoms equal 15*15*15*${density}

### 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.001

### 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
unfix nvt
fix nve all nve

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp press density
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 etotal equal etotal
variable etotal2 equal etotal*etotal
fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_dens2 v_temp2 v_press2 v_etotal 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[7]
variable aveetotal2 equal f_aves[8]
variable errdens equal sqrt(f_aves[4]-f_aves[1]*f_aves[1])
variable errtemp equal sqrt(f_aves[5]-f_aves[2]*f_aves[2])
variable errpress equal sqrt(f_aves[6]-f_aves[3]*f_aves[3])
variable heatcap equal ${atoms}*${atoms}*((f_aves[8]-f_aves[7]*f_aves[7])/(f_aves[2]*f_aves[2]))

print "Averages"
print "--------"
print "Density: ${avedens}"
print "Stderr: ${errdens}"
print "Temperature: ${avetemp}"
print "Stderr: ${errtemp}"
print "Pressure: ${avepress}"
print "Stderr: ${errpress}"
print "Heat capacity: ${heatcap}"

</pre>

[[File:Heatcap2.PNG]]

==Structural properties and the radial distribution function==

[[File:RDFs.PNG]]

The solid is face-centred cubic and so when a reference atom is defined, there are 3 different distances that the other atoms in the unit cell can be away from each other. If the reference atom is in the middle of the top face of the cube, out of the first three peaks, the smallest (and second) peak will be from the atom in the middle of the opposite face of the cell because there are 4 atoms that fit this description, also the RDF has also been normalised and this is a relationship between two of the same lattice points and so it makes sense that the density at this point would be 1. The first peak will be the atoms on the corners of the same face as the reference atom and the atoms that are in the middle of the cell’s remaining faces as they are all the same distance from the reference atom; this is because the distance as you go to an edge and go down is the same as if you go to the edge and go across since the cell is a cube, also there are 12 such atoms and this is the tallest peak. The third peak arises from the atoms that are on the corners of the opposite face of which there are 8. The numbers quoted all arise from considering the actual lattice rather than the cell and also compare well with the relative heights of the peaks. 


The RDF for liquids and gases look quite similar to each other but different to that of the solid as they are much less ordered. In both cases, the first peak is quite large before the curve quickly levels out at 1. The RDF measures the amount of atoms in a shell of radius, r, around the reference atom and so when the shell is small the concentration of nearby atoms is high and they resemble clusters. However, as the shell gets bigger it starts to become more diffuse and eventually reaches a point where the concentration of atoms in the shell cannot be distinguished from the density of the system as a whole. 


By eye, the liquid and gas RDFs would seem to suggest that the gas is ever so slightly denser than the liquid, however, the functions are normalised and the actual density of gases is much lower and so there are far fewer atoms than in the liquid, as can be seen in the plots of their integrations.

[[File:Integrations.PNG]]

From the log file, the lattice spacing of the solid is 

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


==Dynamical properties and the diffusion coefficient==

[[File:MSDgas.PNG]]
<math>
D = \frac{1}{6} \times 0.0025 = 4.17 \times 10^{-4} m^2 timestep{^-1}
</math>

[[File:MSDliq.PNG]]
<math>
D = \frac{1}{6} \times 0.001 = 1.66 \times 10^{-4} m^2 timestep^{-1}
</math>


[[File:MSDsolid.PNG]]
<math>
D = \frac{1}{6} \times 0.0025 = 1.33 \times 10^{-6} m^2 timestep^{-1}
</math>

[[File:Millison msd gas.PNG]]
<math>
D = \frac{1}{6} \times 0.0305 = 0.00508 m^2 timestep^{-1}
</math>

[[File:Million msd liquid.PNG]]
<math>
D = \frac{1}{6} \times 0.001 = 1.66 \times 10^{-4} m^2 timestep^{-1}
</math>

[[File:Million msd solid.PNG]]
<math>
D = \frac{1}{6} \times 6 \times 10^{-8} = 1 \times 10^{-8} m^2 timestep^{-1}
</math>



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

<math> x'\left(t\right) = v(t) = -A\omega\sin\left(\omega t + \phi\right)</math>

<math> v^2(t) = A^2\omega\sin^2\left(\omega t + \phi\right)</math>          

<math> v(t+\tau) = -A\omega\sin\left(\omega t + \omega\tau + \phi\right)</math>  define: <math>\alpha = \omega t + \phi</math>

<math> v(t+\tau) = -A\omega\sin\left(\alpha + \omega\tau \right) = - A\omega\sin\left(\alpha\right)\cos\left(\omega\tau\right)\sin\left(\omega\tau\right)\cos\left(\alpha\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} = \frac{\int_{-\infty}^{\infty} A^2\omega^2\sin^2\left(\alpha\right)\cos\left(\omega\tau\right) + \sin\left(\omega\tau\right)\cos\left(\alpha\right)\mathrm{d}\alpha}{\int_{-\infty}^{\infty}A^2\omega\sin^2\left(\alpha\right)\mathrm{d}\alpha}</math>