Rep:Mod: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 $t$, "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 $x\left(t\right)=A\cos (\omega t+\varphi )$ (the values of $A$, $\omega$, and $\varphi$ 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.

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. If the total energy does not change by more than 1%, the minimum energy should be 0.5*(1-1%)=0.495. Therefor the timestep is required to change to 0.2 (or less than 0.2). Here is the new diagram when timestep equals to 0.2.

3. For a single Lennard-Jones interaction, $\varphi \left(r\right)=4ϵ(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})$, find the separation, $r_0$, at which the potential energy is zero. What is the force at this separation? Find the equilibrium separation, $r_{eq}$, and work out the well depth ($\varphi \left(r_{eq}\right)$). Evaluate the integrals $\underset{2\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r$, $\underset{2.5\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r$, and $\underset{3\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r$ when $\sigma =ϵ=1.0$.

1). when potential energy=0,

$(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})$=0. So,$r_0$= $\sigma$.

2). energy is zero, so net force is zero as well.

3). at $r_{eq}$, potential energy is minimum, $\frac{d\varphi }{dr}$ = 0,

$\frac{48ϵ{\sigma }^{12}}{{\sigma }^{13}}$ = $\frac{24ϵ{\sigma }^6}{{\sigma }^7}$

So, $r_{eq}$ = $\sqrt{2}\sigma$

4). $\underset{2\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r$= $\frac{-4r^{-11}}{11}+\frac{4r^{-5}}{5}=-0.0248$

$\underset{2.5\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r$=$-8.177*10^{-3}$

$\underset{3\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r$=$-3.29*10^{-3}$

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/c{m}^3$, no.(water)= $\frac{1}{18}*6.023*10^{23}$= $3.346*10^{22}$

2). $\frac{1}{3.346*10^{22}}$=$\frac{volume(H2O)}{10000}$, so volume of 10000 water molecules is $2.98*10^{-19}$

5. 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 new position $(0.5+0.7,0.5+0.6,0.5+0.2)=(1.2,1.1,0.7)$.

However it follows the periodic boundary condition, so it should be $(0.2,0.1,0.7)$.

6.The Lennard-Jones parameters for argon are $$\begin{gathered} \sigma =0.34\mathrm{n}\mathrm{m},ϵ \\ / \\ {k}_B=120\mathrm{K} \end{gathered}$$. If the LJ cutoff is $r^{*}=3.2$, what is it in real units? What is the well depth in $$\begin{gathered} {\mathrm{k}\mathrm{J} \\ \mathrm{m}\mathrm{o}\mathrm{l}}^{-1} \end{gathered}$$? What is the reduced temperature $T^{*}=1.5$ in real units?

1) r= 0.34* 3.2 nm = 1.088nm

2) $$\begin{gathered} ϵ \\ / \\ {k}_B=120\mathrm{K} \end{gathered}$$

$ϵ=120*1.38*10^{-23}*10^{-3}*6.023*10^{23}kJ/mol=0.9974kJ/mol$

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 in a face-centred cubic.

$\sqrt[3]{\frac{4}{1.2}}=1.4938$

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?

4*(10*10*10) = 4000 atoms will be created.

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.

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

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:

timestep 0.001
run 100000

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.

However, 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).

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

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.

Pressure	Temperature
2.7 , 2.8	1.8 , 1.9 , 2.0 , 2.1 , 2.2

2. We need to choose $\gamma$ so that the temperature is correct $T=\mathfrak{𝔗}$ if we multiply every velocity $\gamma$. We can write two equations:

$\frac{1}{2}\sum _i{m}_i{v}_i^2=\frac{3}{2}N{k}_{B}T$

$\frac{1}{2}\sum _i{m}_i{\left(\gamma v_i\right)}^2=\frac{3}{2}N{k}_B\mathfrak{𝔗}$. Solve these to determine $\gamma$.

So, $\sum _i{m}_i{\left(\gamma v_i\right)}^2=3N{k}_B\mathfrak{𝔗}$

$\gamma =\frac{1}{\sum _i{v}_i}\sqrt{\frac{3N{k}_B}{\sum _i{m}_i}}$

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?

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, $pV=N{k}_{B}T$

$density=\frac{N}{V}=\frac{p}{k_{B}T}=\frac{p}{T}$, $k_B=1$(reduced unit)

The simulated density is lower than the ideal density. Because there are interaction forces existing among particles. When the pressure is applied, the repulsion between particles becomes dominant, and the distance between particles will be larger, so the density will be lower than the ideal situation.

Calculating heat capacities using statistical physics

Task

1. input script

### DEFINE PARAMETERS ###
variable D equal 0.2
variable T equal 2.8
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
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}/(f_aves[5])
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}"

2. Plot C_V/V as a function of temperature

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 g(r) and $\int g(r)\mathrm{d}r$. Plot the RDFs for the three systems on the same axes, and attach a copy to your report.

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

g(r) shows several peaks in the liquid state, which means the probability at these radial distances will be much higher than the others. Secondly, there are much more peaks observed in solid phases. The peak differences are caused by the different interaction forces. In liquid phases, the attractive force is dominant, and in the solid phase, the structure is regular and periodic, the interaction is stronger than in liquid phase, and the repulsive force is dominant. For face centred cubic, regarding (0,0,0) as a center, the first three peaks correspond to (1/2,1/2,0), (1,0,0) and (1,1,0). The difference from 0 to the second peak is the lattice point, 1.475.

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

MSD(SOLID)	MSD(LIQUID)	MSD(GAS)

Because $MSD=⟨(x(t)-x_0{)}^2⟩=2Dt$,

$D=\frac{MSD}{2t}$

For solid, D = 0.0273/(2*0.016991) = 0.803

For liquid, D = 2.81/(2*1.67448) = 0.839

For gas, D = 12.9/(2*7.74)= 0.8333

the MSD data given

$D=\frac{1}{3}\underset{0}{\overset{\mathrm{\infty }}{\int }}C\left(\tau \right)$

2. 1). 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!):

$C\left(\tau \right)=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}v\left(t\right)v(t+\tau )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{v}^2\left(t\right)\mathrm{d}t}$

Be sure to show your working in your writeup.

For a 1D harmonic oscillator

$x\left(t\right)=A\cos (\omega t+\varphi )$

Therefore

$v\left(t\right)=-A\omega \sin (\omega t+\varphi )$

and

$v(t+\tau )=-A\omega \sin (\omega (t+\tau )+\varphi )$

Because $C\left(\tau \right)=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}v\left(t\right)v(t+\tau )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{v}^2\left(t\right)\mathrm{d}t}$

For the upper integration:

$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}v\left(t\right)v(t+\tau )dt$

$=A^2{\omega }^2\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}sin(\omega (t+\tau )+\varphi )sin(\omega t+\varphi )dt$

$=A^2{\omega }^2\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}sin(\omega t+\varphi )(sin(\omega t+\varphi )cos(\omega \tau )+sin(\omega \tau )cos(\omega t+\varphi ))dt$

$=A^2{\omega }^{2}cos(\omega \tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^2(\omega t+\varphi )dt+A^2{\omega }^{2}sin(\omega \tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}sin(\omega t+\varphi )cos(\omega t+\varphi )dt$

For the lower integration:

$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{v}^2\left(t\right)dt$

$=A^2{\omega }^2\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^2(\omega t+\varphi )dt$

So,

$C\left(\tau \right)=cos(\omega \tau )+sin(\omega \tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\frac{cos(\omega t+\varphi )}{sin(\omega t+\varphi )}$

$=cos(\omega \tau )+sin(\omega \tau )\frac{ln(|sin(\omega t+\varphi )|)}{\omega }$

from $-\mathrm{\infty }$ to $\mathrm{\infty }$,

$C\left(\tau \right)=cos(\omega \tau )$

2). On the same graph, with x range 0 to 500, plot $C\left(\tau \right)$ with $\omega =1/2\pi$ and the VACFs from your liquid and solid simulations.

$C\left(\tau \right)=cos(\frac{1}{2}\pi \tau )+2\frac{sin(\frac{1}{2}\pi \tau )}{\pi }$

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

The minima represents the starting points of self-diffusion. The origin of the differences is the interaction force difference. The interaction is larger , so the diffusion coefficient is smaller, which can be also observed by the area difference under the VACF(r) curve of two different phases.