View source for Rep:Tfy286

==Introduction to molecular dynamics simulation==
===Task===
1. 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).

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.

[[file:A.e.e.PNG]]

2. 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?
  
''''The maximum energy is 0.5, so the minimum energy should be 0.5*(1-1%)=0.495. Therefor the timestep is required to change to 0.2 , and here is the new diagram.''''

[[file:Energy.change.PNG|x400px]]

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

1). when potential energy=0, <math>\left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>=0, so <math>r_0</math>= <math>\sigma</math>.

2). energy is zero, so net force is zero as well.
  
3). at <math>r_{eq}</math>, potential energy is minimum, so <math> \frac{d\phi}{dr}\ </math> = 0, so,<math>\dfrac{48\epsilon\sigma^{12}}{\sigma^{13}} </math> = <math>\dfrac{24\epsilon\sigma^{6}}{\sigma^{7}} </math>, so <math>r_{eq}</math> = <math>\sqrt{2}\sigma</math>

4). <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>= <math>\dfrac{-4r^{-11}}{11}+\dfrac{4r^{-5}}{5}</math>= -0.0248

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

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

4. Estimate the number of water molecules in 1ml of water under standard conditions. Estimate the volume of 10000 water molecules under standard conditions.

1). density of water is 1g/cm^3, no.(water)= <math>\frac{1}{18}* 6.023* {10^{23}}</math>= <math>3.346*10^{22}</math>
2). <math>\dfrac{1}{3.346*10^{22}}</math>=<math>\dfrac{volume(H2O)}{10000}</math>, so volume of 10000 water molecules is  <math>2.98*10^{-19}</math>

5. 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''?

the new position <math>\left(0.5+0.7, 0.5+0.6, 0.5+0.2\right)=\left(1.2, 1.1,0.7\right)</math>, 
because it follows the periodic boundary condition, so it should be <math>\left(0.2, 0.1,0.7\right)</math>

6.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?


1) r= 0.34* 3.2 nm = 1.088nm
2) 
3) T= 120* 1.5 K = 180K



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

If two atoms are too close together, there will be a huge repulsion existing between them. In order to calculate the huge energy change, the timestep should be very small. Therefore the computation process will be much more difficult. 

2. Satisfy yourself that this lattice spacing corresponds to a number density of lattice points of 0.8. 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?
there are four lattice points per volume per unit cell.

<math>\sqrt[3]{\tfrac{4}{1.2}} = 1.4938</math>

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?

create 4*1000 = 4000 atoms. 

4. Using the LAMMPS manual, find the purpose of the following commands in the input script:
   mass 1 1.0
1 indicates the atom type, here is type 1.
1.0 indicates the atom mass, here is 1.0
   pair_style lj/cut 3.0
pair_style lj/cut means cutoff Lennard-Jones potential with no Coulomb
3.0 is the arguments used by a particular style 
   pair_coeff * * 1.0 1.0
asterisks here are used in place of or in conjunction with the I,J arguments to set the coefficients for multiple pairs of atom types. 
1.0 and 1.0 here are the coefficients for the atom types

5. Given that we are specifying \mathbf{x}_i\left(0\right) and \mathbf{v}_i\left(0\right), which integration algorithm are we going to use?
velocity-Verlet algorithm 

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>

if the timestep is fixed at 0.001, once the timestep is required to change, all the timesteps will change into other values by hand.
if the timestep is set as a variable, only the first line is required to change, and the rest of them will change automatically. 

7.1) make plots of the energy, temperature, and pressure, against time for the 0.001 timestep experiment (attach a picture to your report).

[[File:E.t.p]]

2) Does the simulation reach equilibrium? How long does this take? 

Yes, the simulation reaches equilibrium when time= 0.43

3) make a single plot which shows the energy versus time for all of the timesteps (again, attach a picture to your report).

[[File:Energy.timestep.PNG|x600px]]

4)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?

0.01 is the largest to give acceptable result. And 0.015 is a particular bad choice because there is a big deviation from the others.

==Running simulations under specific conditions==
===Task===
1. Choose 5 temperatures (above the critical temperature T^* = 1.5), 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 \left(p, T\right) 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.
{| class="wikitable"
! style="background: #0D4F8B; color: white;" | Pressure
! style="background: #0D4F8B; color: white;" | Temperature
|-
| 2.7 , 2.8
| 1.8 , 1.9 , 2.0 , 2.1 , 2.2
|}

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

<math>\gamma = \frac{1}{\sum_i v_i} \sqrt{\dfrac{3 N k_B}{\sum_i m_i}}</math>

3. Use the manual page 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?

100 stands for the number of input values every this many timesteps

1000 stands for the number of times to use input values for calculating averages

100000 stands for calculate averages every this many timesteps

Therefore the value of temperature will sampled for 1000 times. And there are 1000 measurements contributing to the average.

4. 1). 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 \left(\frac{N}{V}\right). 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. Does the discrepancy increase or decrease with pressure?

[[File:Density. temperature.PNG]]

The higher the pressure the discrepancy increases. 

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

according to the ideal gas law, <math> p V = N k_B T </math>

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



==Calculating heat capacities using statistical physics==
===Task===
1. input script
 <pre>### DEFINE PARAMETERS ###
variable D equal 0.8
variable T equal 2.0
variable timestep equal 0.0075

### DEFINE SIMULATION BOX GEOMETRY ###
lattice sc ${D} 
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 10000
reset_timestep 0
unfix nvt
fix nve all nve

### MEASURE SYSTEM STATE ###
thermo_style custom step etotal temp press density atoms
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 atomss equal atoms
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 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 energyvariance equal (f_aves[8]-f_aves[7]*f_aves[7])
variable heatcap equal ${atomss}* ${avedens}* ${energyvariance}/(${avetemp}*${avetemp})
print "Averages"
print "--------"
print "Density: ${avedens}"
print "Stderr: ${errdens}"
print "Temperature: ${avetemp}"
print "Stderr: ${errtemp}"
print "Pressure: ${avepress}"
print "Stderr: ${errpress}"
print "heatcap per volume: ${heatcap}"
</pre>

2. Plot C_V/V as a function of temperature

[[File:Heat capacity.PNG]]

==Structural properties and the radial distribution function==
===Task===
For a 1D harmonic oscillator

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

Therefore 

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

and

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

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

For the upper integration:

<math> \int_{-\infty}^{\infty} v\left(t\right)v\left(t + \tau\right) dt </math>

<math> = A^2 \omega^2 \int_{-\infty}^{\infty} sin\left(\omega (t+\tau) + \phi\right) sin\left(\omega t + \phi\right) dt </math>

<math> = A^2 \omega^2 \int_{-\infty}^{\infty} sin\left(\omega t + \phi\right) (sin\left(\omega t + \phi\right) cos(\omega\tau)+ sin\left(\omega\tau)cos(\omega t + \phi\right)) dt</math>

<math> = A^2 \omega^2 cos(\omega\tau)\int_{-\infty}^{\infty} sin^2 \left(\omega t + \phi\right) dt + A^2 \omega^2 sin(\omega\tau)\int_{-\infty}^{\infty}sin\left(\omega t + \phi\right) cos\left(\omega t + \phi\right) dt</math>


For the lower integration:

<math> \int_{-\infty}^{\infty} v^2\left(t\right) dt </math>

<math> = A^2 \omega^2 \int_{-\infty}^{\infty} sin^2 \left(\omega t + \phi\right) dt</math>


So, 

<math>C\left(\tau\right)= cos(\omega\tau) + sin(\omega\tau) \int_{-\infty}^{\infty}\dfrac{cos(\omega t + \phi)}{sin\left(\omega t + \phi\right)}</math>

<math> =cos(\omega\tau) + sin(\omega\tau)\dfrac{ln(|sin\left(\omega t + \phi\right)|)}{\omega}</math>