Talk:Mod:NTH0510

JC: General comments: All tasks attempted and most completed, but some explanations were missing. In general try to make your explanations more detailed, complete and clearer. Also take a bit more care with writing your report to check the text and spellings.

Theory

Numerical Integration

TASK 1

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

• ANALYTICAL was calculated by the equation of $A\cos (\omega t+\varphi )$ (Figure 1)

File:ANALYTICAL 4.png

Figure.1 Analytical Positions VS Time

• ERROR was calculated by the equation of $|ANALYTICAL-x\left(t\right)|$ (Figure 2)

File:Error 4.png

Figure.2 Error VS Time

• ENERGY

Total energy= Potential Energy+ Kinetic energy (Figure 3)

For harmonic oscillation, potential energy obeys the Hooke’s law, and $PE=\int -Kxdx=1/2K{x}^2$

For kinetic energy, $KE=1/2m{v}^2$

File:EVST.png

Figure.3 Energy VS Time

TASK 2

TASK: 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.

The maxima error were extracted from time ranges of 0-3.2, 3.2-6.3, 6.3-9.4, 9.4-14.2 respectively. The plot of maxima error against time and the line of best fit are shown in Figure 4 below.

File:MAXIMA VS TIME.png

Figure.4 Maxima VS Time

Maxima	Time
12	0.000758457
4.9	0.00200771
8	0.00330076
11.1	0.004603854
14.2	0.005910508

JC: Is the linear fit really perfect (R^2 = 1)? Give the line of best fit to more than 1 s.f.. Why does the maximum error increase over time?

TASK 3

TASK: 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?

Nomally, the totally energy of the system is constant during harmonic oscillation. The fluctuation of the energy was due to the algorithm we used was not perfect. In order to minimise the change of energy is less than 1%, we need to change the timestep. As can be seen from Figure 5, when timestep equalled 0.06 or less, the energy difference was less than 1%.

File:EVSTS.png

Figure.5 Total Energy VS Timestep

JC: Change the x axis labels on figure 5 to show the size of the energy variation. What other timesteps did you try? You should find that a timestep around 0.2 produces energy fluctuations of ~1%.

Atomic Forces

TASK 4

• TASK: 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.

The separation at which the potential energy is zero,implying the at that point, the single Lennard-Jones interaction needs to be zero. Single Lennard-Jones interaction $\varphi \left(r\right)=4ϵ(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})$

Potential energy is zero, when $r$=$r_0$ = $\sigma$

• What is the force at this separation?

$F$ =$-\frac{\mathrm{d}\varphi \left(r\right)}{\mathrm{d}\mathbf{r}}$

$F$ = $-4ϵ(-12\frac{{\sigma }^{12}}{r^{13}}+6\frac{{\sigma }^6}{r^7})$

$F$ = $-24ϵ(2\frac{{\sigma }^{12}}{r^{13}}-\frac{{\sigma }^6}{r^7})$

Substitute $r$=$r_0$ = $\sigma$ into the equation

$F$ = $24\frac{ϵ}{\sigma }$

• Find the equilibrium separation.

When atoms are in equilibrium state, the potential energy is at minimum. In order to find the minimum potential energy, we need to differentiate the Lennard-Jones potential energy.

$\frac{\mathrm{d}\varphi \left(r\right)}{\mathrm{d}\mathbf{r}}=0$ ; $\varphi \left(r\right)=4ϵ(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})$

$\frac{\mathrm{d}\varphi \left(r\right)}{\mathrm{d}\mathbf{r}}=24ϵ(2\frac{{\sigma }^{12}}{r^{13}}-\frac{{\sigma }^6}{r^7})=0$

$$\begin{gathered} 2{\left(\frac{\sigma }{ \\ r}\right)}^6=1 \end{gathered}$$

$r=r_{eq}=2^{\frac{1}{6}}\sigma$

*Work out the well depth ($\varphi \left(r_{eq}\right)$).

$\varphi \left(r\right)=4ϵ(\frac{{\sigma }^{12}}{r_{eq}^{12}}-\frac{{\sigma }^6}{r_{eq}^6})$ ; $r=r_{eq}=2^{\frac{1}{6}}\sigma$

$\varphi \left(r\right)=4ϵ(\frac{{\sigma }^{12}}{{\left(2^{\frac{1}{6}}\sigma \right)}^{12}}-\frac{{\sigma }^6}{{\left(2^{\frac{1}{6}}\sigma \right)}^6})$

$\varphi \left(r\right)=4ϵ(\frac{1}{4}-\frac{1}{2})$

$\varphi \left(r\right)=ϵ-2ϵ$

$\varphi \left(r\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$.

$\underset{a\sigma }{\overset{b}{\int }}\varphi \left(r\right)\mathrm{d}r=\underset{a\sigma }{\overset{b}{\int }}4ϵ(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})\mathrm{d}r$, where a and b are boundary constant and $\sigma =ϵ=1.0$.

$\underset{a\sigma }{\overset{b}{\int }}\varphi \left(r\right)\mathrm{d}r=-\frac{4ϵ\sigma }{11}(\frac{1}{b^{11}}-\frac{1}{a^{11}})+\frac{4ϵ\sigma }{5}(\frac{1}{b^5}-\frac{1}{a^5})$

When $a=2$, $b=\mathrm{\infty }$, $\sigma =ϵ=1.0$

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

When $a=2.5$, $b=\mathrm{\infty }$, $\sigma =ϵ=1.0$

$\underset{2.5\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=-\frac{4ϵ\sigma }{11}(0-\frac{1}{2.5^{11}})+\frac{4ϵ\sigma }{5}(0-\frac{1}{2.5^5})=-0.008177$

When $a=3$, $b=\mathrm{\infty }$, $\sigma =ϵ=1.0$

$\underset{3\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=-\frac{4ϵ\sigma }{11}(0-\frac{1}{3^{11}})+\frac{4ϵ\sigma }{5}(0-\frac{1}{3^5})=-0.003290$

JC: All maths correct and clear working.

TASK 5_Periodic Boundary Conditions

TASK: Estimate the number of water molecules in 1ml of water under standard conditions.

• The density of water is 1g/ml, therefore, for 1ml of water, mass is 1g. According to $n=\frac{m}{M}$,$M=18g/mol$, the number of mole for 1 ml of water is 0.056 mol.

The number of water molecules equal to moles time advogado number ($6.02\times 10^{23}$), thus,the number of water molecules in 1ml of water is $3.371\times 10^{23}$.

JC: You are out by an order of magnitude, 0.056 x 6.02 x 10^23 = 3.37 x 10^22

• Estimate the volume of 10000 water molecules under standard conditions.

$\frac{3.34*10^{22}}{1mL}=\frac{10000}{V_{10000}}$

$V_{10000}=2.994*10^{-19}$

JC: Correct.

TASK 6_Periodic Boundary Conditions

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?

Originally, the atom was at position (0.5, 0.5, 0.5), then it moved along the vector(0.7, 0.6, 0.2), the subsequent new position was (1.2, 1.1, 0.7). The central box was repeated infinitely in all directions, and when an atom crosses the boundary of the box, one of its replicas enters the box through the opposite face. Therefore, the resultant position should be (0.2, 0.1, 0.7).

JC: Correct.

TASK 7_Reduced Units

TASK: 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? distance $r^{*}=\frac{r}{\sigma }$ energy $E^{*}=\frac{E}{ϵ}$ temperature $T^{*}=\frac{k_{B}T}{ϵ}$

• Distance

$r=r^{*}\sigma$ = $3.2\times 0.34\mathrm{n}\mathrm{m}$ = $1.088nm$

• The well depth

$ϵ=k_B\times 120K$ = $1.38\times 10^{-23}J/K\times 120K$= $1.656\times 10^{-21}J$ $3.371\times 10^{23}$

The well depth $\varphi \left(r\right)=4ϵ(\frac{{(0.34\mathrm{n}\mathrm{m})}^{12}}{{(1.088nm)}^{12}}-\frac{{(0.34\mathrm{n}\mathrm{m})}^6}{{(1.088nm)}^6})$ = $6.16\times 10^{-24}J$

• Temperature T = $$\begin{gathered} ϵT^{*} \\ / \\ {k}_B=180K \end{gathered}$$

JC: Values are correct, but the energy should have been given in kJ mol-1. Also don't give answers to more s.f. than used in the question (2 in this case).

Equilibration

Creating the simulation box

TASK 8

TASK: 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 atoms are given randomly at starting coordinates, two atoms got the chance to remain at the same position or approach each other which is unlikely to happen. There will be strong repulsion between two atoms if they are generated close together originally.

JC: Why is strong repulsion between atoms in a starting configuration a problem? A simulation may crash or require a much smaller timestep if starting from a configuration with high repulsive forces.

TASK 9

TASK: 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?

In a simple cubic lattice, there is 1 atom in the unit cell, while lattice spacing in x,y,z = 1.07722 1.07722 1.07722. The volume of unit cell=$1.07722^3=1.25000$ The number density of unit cell=$\frac{1}{1.25000}=0.8$

In a face-centered cubic lattice, there are 4 atoms in a unit cell, and the number density is 1.2.

Assuming lattice spacing is x

Therefore

$1.2=\frac{4}{x^3}=0.8$

$x=1.4938$

TASK 10

TASK: 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?

For a box containing 1000 unit cells, and 4 atoms in each unit cell, there will be 4000 atoms creat by the command.

JC: Correct.

Setting the properties of the atoms

TASK 11

TASK: Using the LAMMPS manual, find the purpose of the following commands in the input script:

mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0

• mass 1 1.0

mass 1 1.0 correspond to syntax of mass I value, where I represents the atom type, and value is the mass. In this case means 1 type of atom with mass equals to 1. The mass command aims to set the mass for all atoms of one or more atom types. The mass command can only be used if the atom style requires per-type atom mass to be set.

• pair_style lj/cut 3.0

Set the formula(s) LAMMPS uses to compute pairwise interactions. In LAMMPS, pair potentials are defined between pairs of atoms that are within a cutoff distance and the set of active interactions typically changes over time.

A command sets the 'pair style' to lj/'cut style, which is the standard 12/6 Lennard-Jones potential:

$E=4ϵ(\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})$, where $r<r_c$

$r_c$ is the cutoff for Lennard Jones interaction. In this case $r_c=3.0$. The cutoff means the the bond length where the Lennard-Jones potential is nearly zero.

JC: The potential is 0 at the cut off, why is a cut off used?

• pair_coeff * * 1.0 1.0

Specify the pairwise force field coefficients for one or more pairs of atom types. The number and meaning of the coefficients depends on the pair style.

For all of the lj/cut pair styles, pair_coeff command must be used to define coefficients for each pair of atoms types.

A wildcard asterisk (* *) is used to set the coefficients for all pairs of atom types, then followed by the values indicating the coefficients for the pairs of atom types. The coefficients indicate epsilon and sigma in the Lennard-Jones respectively, in this case are (1.0 1.0).

JC: Good remark that the pair coefficients in this case are epsilon and sigma.

TASK 12

TASK: 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?

The velocity Verlet algorithm.

JC: What would be needed for the standard Verlet algorithm to be used?

Running the simulation

TASK 13

TASK: 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

The previous version is more convenient as timestep is the only variable that we need to type; all other variables are dependent on the timestep. While for the later one, all various is independent on the timestep, thus we need to change every line to base on what early one does.

JC: Correct, changing every line would be time consuming and make errors more likely.

Checking equilibration

TASK 14

*Make plots of the energy, temperature, and pressure, against time for the 0.001 timestep experiment (attach a picture to your report). Does the simulation reach equilibrium? How long does this take?

Plots of the energy (Figure 6), temperature (Figure 7), and pressure (Figure 8), against time for the 0.001 timestep experiment. As can be seen from the graph below, the simulation reaches equilibrium after 0.3. The physical properties become stable after 0.3.

File:Total energy vs time.png

Figure.6 Total energy vs time (timesteps=0.001)

File:Temp vs time.png

Figure.7 Temperature vs time (timesteps=0.001)

File:Press vs time.png

Figure.8 Pressure vs time (timesteps=0.001)

*Make a single plot which shows the energy versus time for all of the timesteps.

File:AllE vs time.png

Figure.9 Energy VS Time under different timesteps

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

According to Fig.9, the largest acceptable timestep is 0.0025, as it has the same trend as 0.001, showing that it sufficient accuracy but not as short as 0.001, so we have longer time for simulation. 0.015 is a particularly bad choice for timestep because the total energy is not constant.

JC: Why are 0.0075 and 0.01 not suitable timesteps? Should total energy depend on the choice of timestep? Also the legend of your graph says 0.025, not 0.0025.

Running simulations under specific conditions

Temperature and Pressure Control

TASK 15_Thermostats and Barostats

TASK: 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$ (1)

$\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$.

$\frac{1}{2}\sum _i{m}_i{\gamma }^2{v}_i^2=\frac{3}{2}N{k}_B\mathfrak{𝔗}$

Since $\gamma$ is a constant:

${\gamma }^2\left(\frac{1}{2}\sum _i{m}_i{v}_i^2\right)=\frac{3}{2}N{k}_B\mathfrak{𝔗}$ (2)

Substituted (1) into (2): $\frac{1}{2}\sum _i{m}_i{v}_i^2=\frac{3}{2}N{k}_{B}T$

We get: ${\gamma }^2\left(\frac{3}{2}N{k}_{B}T\right)=\frac{3}{2}N{k}_B\mathfrak{𝔗}$

${\gamma }^2=\frac{\mathfrak{𝔗}}{T}$

$\gamma =\sqrt{\frac{\mathfrak{𝔗}}{T}}$

JC: Good, clear derivation.

TASK 16_Temperature and Pressure Control

TASK: 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 (p, T) points. You should be able to use the results of the previous section to choose a timestep.

The chosen conditions are shown below:

• Pressure: 2.6, 2.7 (selected values are base on previous section)

• Temperature: 1.6, 1.8, 2.0, 2.2, 2.4

• Timestep: 0.0025 (selected timestep are base on previous section)

TASK 17_Examining the Input Script

TASK: 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?

fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_dens2 v_temp2 v_press2
run 100000

• 100:Nevery (use input values every 100 timesteps)

• 1000:Nrepeat (use input values for calculating averages for 1000 times)

• 100000:Nfreq (calculate averages every 100000 timesteps)

In this case, Nevery=100, Nrepeat=1000, and Nfreq=100000. The values are taken every 100 timesteps. Fore examole, the values are taken on timesteps 100, 200,...100000, which will be used to compute the final average on every 100000 timestep. In other words, 1000 samples are obtained to calculate the average. The next line states 'run 100000', which means the simulation will run 100000 timesteps.

TASK 18_Plotting the Equations of State

TASK: 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. 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. Does the discrepancy increase or decrease with pressure?

• According to Figure.10, the higher the pressure, the higher the density, as higher pressure lead to a smaller volume. The is due to density equals to $\left(\frac{N}{V}\right)$, smaller the volume resulting in greater density.

File:Density vs Temperature 1.png

Fig.10 Density as a function of temperature for P=2.6 and P=2.7

JC: Why are there 2 colours for each pressure in the legend of figure 10?

• Figure. 11 indicates the comparison between simulated density and the predicted density from Ideal gas law. Figure. 10 also shows the trend of higher the pressure the higher the density and the simulated densities are smaller than the predicted values from ideal gas law. These observation can be explained as there is no interaction exist in ideal gas model; however, in real gas, interaction such as repulsion and attraction exits, meaning that the molecules cannot approach closely to each other. Thermodynamically, for ideal gas, pressure can be expressed in terms of entropy.:

$P=T{\left(\frac{\delta S}{\delta V}\right)}_{U,N}$

Ideal gas obeys the equaiton of $PV=nRT$ (3)

Since $n=\frac{N}{N_a}$ (4) ; $R=N_a{K}_b$ (5), where $K_b=1$ in reduced unit

Substituting (4) and (5) into (3), we get

$PV=N{k}_{B}T$ (6)

Rearrange (6)

$density=\left(\frac{N}{V}\right)=\frac{P}{K_{b}T}$

File:Density vs temperature ideal.png

Figure.11 The comparison between the densities of ideal gas and simulated gas as a function on temperature under P=2.6 and P=2.7.

JC: You should compare your simulation results with the values for the ideal gas. How do the differences between the results change with temperature and pressure?

Calculating Heat Capacities Using Statistical Physics

TASK 19_Heat Capacity Calculation

TASK: As in the last section, you need to run simulations at ten phase points. In this section, we will be in density-temperature $({\rho }^{*},T^{*})$ phase space, rather than pressure-temperature phase space. The two densities required at $0.2$ and $0.8$, and the temperature range is $2.0,2.2,2.4,2.6,2.8$. Plot $C_V/V$ as a function of temperature, where $V$ is the volume of the simulation cell, for both of your densities (on the same graph). Is the trend the one you would expect? Attach an example of one of your input scripts to your report.

• The script attached below is the my input file for ${\rho }^{*}=0.8T^{*}=2.2$ :

### DEFINE SIMULATION BOX GEOMETRY ###
lattice sc 0.8
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

### SPECIFY THE REQUIRED THERMODYNAMIC STATE ###
variable T equal 2.2
variable timestep equal 0.0025

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

### MEASURE SYSTEM STATE ###
unfix nvt
fix nve all nve
thermo_style custom step etotal temp press density vol atoms
variable temp equal temp
variable temp2 equal temp*temp
variable etotal equal etotal
variable etotal2 equal etotal*etotal
variable vol equal vol
fix aves all ave/time 100 1000 100000 v_temp v_temp2 v_etotal v_etotal2
run 100000

variable avetemp equal f_aves[1]
variable aveeng equal f_aves[3]
variable heatcap equal atoms*atoms*((f_aves[4]-(f_aves[3]*f_aves[3]))/f_aves[2])
variable heatcapvol equal ((atoms*atoms*((f_aves[4]-(f_aves[3]*f_aves[3]))/f_aves[2]))/vol)

print "Averages"
print "--------"
print "Temperature: ${avetemp}"
print "Etotal: ${aveeng}"
print "HeatCapacity: ${heatcap}"
print "HeatCapacityPerVol: ${heatcapvol}"

• Figure.12 shows as the temperature increases, the heat capacity per volume decrease.

Heat capaticy is a quantity which measures the amount of energy is required to raise the temperature of a system.

Heat can be calculated from equation $C_v$ is: $C_V=N^2\frac{⟨E^2⟩-{⟨E⟩}^2}{k_B{T}^2}$

where $k_B=1$ is in reduced unit.

Thus: $\frac{C_V}{V}=\frac{N^2}{V}\frac{⟨E^2⟩-{⟨E⟩}^2}{k_B{T}^2}=(density)N\frac{⟨E^2⟩-{⟨E⟩}^2}{k_B{T}^2}$

As expected from the equaiton above, heat capacity per volume is proportional to the density and inversely proportional to the square of temperature, and it is also the trend that is obtained from Figure.12.

JC: What about the temperature dependence of <E> and <E^2>? Can you explain how heat capacity changes with density?

File:Cv vs V.png

Figure.12 Heat capacity per volume as the function of temperature under density=0.2 and 0.8.

Structural properties and the radial distribution function

TASK 20_Simulations in this section

TASK: 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. 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?

Densities and temperatures for solid, liquid and vapor are indicated below:

Solid: density 1.25, temperature 1.35 fcc

Vapor: density 0.08, temperature 1.15 sc

Liquid: density 0.8, temperature 1.2 sc

File:G(r) VS distance.png

Fig.13 g(r) VS distance for solid, liquid and vapour.

In statistical mechanics, Radial distribution function (RDF) describes the behavior of density as a function of distance relative to reference particle. As can be seen form Figure. 13, the RDF is 0 when r is small (0.025-0.0825), as the atoms cannot approach to each other closely due to the present of repulsion.

• Solid (grey) phase has the highest density and the most ordered comparing to liquid and vapor . The fluctuation indicating atoms are closely packed together as the lattice structure is more ordered, therefore, shorter distance to find a neighbor atom.

• Liquid (blue) phase has lower density than solid but higher than vapor. only 3 obvious peaks can be found on Figure.13 for liquid, as the atom separations are greater than solid, and also less ordered, meaning that reference atom in liquid phase is more difficult to find a neighbor.

• Vapour (orange) phase has the lowest density as well as less ordered in contrast to other two phase. In Fig.13, only 1 peak can be observed, indicating the atom arrangement is highly disordered, and the probability that a reference atom to fin a neighbor is quite small.

As the distance increase further, the average RDF tend to be 1, as at that stage, RDF is dominated by the average density of that system.

JC: The RDF shows the extent of regular order in a system, the peaks in the solid RDF at large distance show that the solid has long range order.

File:Intedralgr.png

Fig.14 Integral of g(r) VS distance for solid, liquid and vapour.

• Fig. 14 shows the integral of g(r) as a function of distance, which is the area under g(r).

Solid has the highest probability to find a neighbour relative to the reference atom, resulting in the greater area under g(r) curve. The area under g(r) for liquid is less then solid then followed by vapour, due to the degree of disorder.

File:Cubic unit cell1.png

Figure.15 fcc cubic unic cell

In Figure. 15, atom a, b, and c are corresponding to the first 3 peak for solid state in Figure.13. Reference atom is indicating in red.

The distance from reference atom to atom a is defined as A.

The distance from reference atom to atom b is calculated to be $\frac{1}{\sqrt{2}}A$.

The distance from reference atom to atom b is calculated to be $\sqrt{\frac{3}{2}}A$.

Distance can be confirmed by determing A from the three peaks and compare with the theoretical A value.

• solid is possessed fcc lattice structure (Figure 15).

The coordination number for atom a is 6.

There are 3 equivalent atom b in a unit cell, and as reference atom is surranded by 8 unit cell---> $CN=\frac{3}{2}*8=12$

There are 3 equivalent atom c in a unit cell, unlike atom b, they do not share with unit cells ---> $CN=3*8=24$

File:Zoomin.png

Figure.16 Integral of g(r) VS distance for solid (r=0 to r=2)

The coordinatio number can be confirmed by the integral the area under g(r). In this case, as only focus in the first 3 peaks, the cure can be zoomed in into the region from r=0 to r=2 (see Figure. 16).

Atom	Peak Distance	Distance relatives to reference atom	Determined A	Coordination Number
a	1.025	A	1.025	6
b	0.976	$\frac{1}{\sqrt{2}}A$	1.381	12
c	1.675	$\sqrt{\frac{3}{2}}A$	1.490	24

File:Zoomin.png

Figure.16 Integral of g(r) VS distance for solid (r=0 to r=2)

JC: Good use of the inegral of g(r) to find the nearest neighbours. Can you comment on why you obtain 3 different values of A from the 3 g(r) peaks, which one is more accurate and closest to the theoretical value?

Dynamical properties and the diffusion coefficient

Mean Squared Displacement

TASK 21

TASK: 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.

The diffusion coefficient is given by the equation shows below:

$D=\frac{1}{6}\frac{\partial ⟨r^2\left(t\right)⟩}{\partial t}$

Thus

D can be calculated by $D=\frac{1}{6}\times$gradient of MSD as a function of time.

File:Msd sol.png

Figure.18 MSD VS Timesteps for Solid.

Figure.17 MSD VS Timesteps for Liquid.

File:Msd vap.png

Figure.19 MSD VS Timesteps for Vapour.

JC: The liquid movement is only diffusive when the MSD curve is linear, so to calculate the diffusion coefficient you should only fit a straight line to the linear part of the graph, not to the entire range.

	Gradient (timestep)	Gradient (time)	D
Solid (32000 atoms)	7.00E-08	3.00E-05	5.00E-06
Solid (1 million atoms)	5.00E-08	3.00E-05	5.00E-06
Liquid (8000 atoms)	1.00E-03	5.09E-01	8.49E-02
Liquid (1 million atoms)	1.00E-03	5.24E-01	8.73E-02
Vapour (8000 atoms)	2.70E-02	1.03E+01	1.72E+00
Vapour (1 million atoms)	3.05E-02	1.52E+01	2.54E+00

The plot of MSD vs Timestep shows trends as expected.

• For solid (Figure 17), MSD shows small fluctuations over the simulation, this is due to atoms in solid are well packed, and they only vibrate within small distance, therefore, the resultant diffusion coefficient is quite small (5.00E-06).

• For liquid (Figure 18), the MDS increase linearly over time, implying a greater freedom of freedom for liquid in contrast to solid.

• For gas(Figure 19), atoms in vapour have the greatest degree of freedom to move. It ca be seen from the plot that there is a rapid increase at the beginning of simulation. The high freedom of movement is resulting in greater diffusion coefficient.

Velocity Autocorrelation Function

TASK 22

TASK: 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}$

The position of a 1D harmonic oscillator is given by $x\left(t\right)=A\cos (\omega t+\varphi )$

Also, the velocity is $v\left(t\right)=\frac{x(t)}{dt}=-A\sin (\omega t+\varphi )$

$C(\tau )=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}(-Aw\sin (wt+\theta ))(-Aw\sin (w(t+\tau )+\theta )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}(-Aw\sin (wt+\theta ){)}^2\mathrm{d}t}$

$C(\tau )=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (wt+\theta )\sin (w(t+\tau )+\theta )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^2(wt+\theta )\mathrm{d}t}$

$C(\tau )=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (wt+\theta )\sin ((wt+w\tau )+\theta )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^2(wt+\theta )\mathrm{d}t}$

Since $\sin (a+b)=\sin (a)\cos (b)+\cos (a)\sin (b)$

$C(\tau )=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (wt+\theta )[\sin (wt+\theta )cos(w\tau )+\cos (wt+\theta )\sin (w\tau )]\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^2(wt+\theta )\mathrm{d}t}$

$C(\tau )=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^2(wt+\theta )\cos (w\tau )\mathrm{d}t+\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (wt+\theta )cos(wt+\theta )\sin (w\tau )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^2(wt+\theta )\mathrm{d}t}$

$C(\tau )=\frac{\cos (w\tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^{}(wt+\theta )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^2(wt+\theta )\mathrm{d}t}+\frac{\sin (w\tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (wt+\theta )\cos (wt+\theta )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{\sin }^2(wt+\theta )\mathrm{d}t}$

Since $\sin (wt+\theta )$ is an odd function and $\cos (wt+\theta )$ is an even function, $\sin (wt+\theta )\cos (wt+\theta )$ is consequently to be an odd function.

Hence $\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\sin (wt+\theta )\cos (wt+\theta )\mathrm{d}t=0$

$C(\tau )=\cos (w\tau )$

JC: Good derivation.

TASK: 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. 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.

File:V overall.png

Figure. 20 Velocity Autocorrelation Function for Solid and Liquid with C(t)

The minima in the VACFs plot of liquid and solid phase indicating the change of directions of movement.

In solid phase, the the change of direction is due to the vibrations of the atoms; however, in the liquid phase, the change is because to collisions of the atoms with each other. As can be seen from Figure.20, VACF decays as time goes by, this observation is due to the velocity decorrelates with time. In other words, the atom 'forget' its initial velocity.

In solid state, all atoms are well and closely pakced, therefore, the position of atoms in solid state is relatively fixed. Vibration is the only motion of solid atom. Atoms oscillates forwards and backwards and the velocities are reserved. This can be revealed from Figure. 20, showing that VACF oscillates strongly from positive to negative and so on. The decay of VACF is also associated with other factors.

Atoms in liquid state have similar behaviors as in solid state, except they can move freely and collide to each other. Therefore, VACF of atoms decay more rapidly in liquid state.

The harmonic oscillator VACF is very different to the Lennard Jones solid and liquid. The harmonic oscillator is assuming there is only one atom in the system and the atom moving forwards and backwards and has no interaction with others, meaning that,there is no energy transfer and energy is constant over time. Hence, the magnitudes of velocity remains constant and only signs change.

TASK 23

TASK: Use the trapezium rule to approximate the integral under the velocity autocorrelation function for the solid, liquid, and gas, and use these values to estimate D in each case. You should make a plot of the running integral in each case. Are they as you expect? Repeat this procedure for the VACF data that you were given from the one million atom simulations. What do you think is the largest source of error in your estimates of D from the VACF?

Since diffusion coefficient can be obtained by integrating this function

$D=\frac{1}{3}\underset{0}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}\tau ⟨\mathbf{v}\left(0\right)\cdot \mathbf{v}\left(\tau \right)⟩$

Thus

$D=\frac{1}{3}\times$running integral as a function of time.

Trapezium rule was employed to approximate the integral under the velocity autocorrelation function for the solid, liquid, and gas. The estimated integral and D are shown on the table below.

	Running integral (time)	D
Solid (32000 atoms)	7.73E-04	2.58E-04
Solid (1 million atoms)	-8.33E-04	-2.78E-04
Liquid (8000 atoms)	2.38E-01	7.94E-02
Liquid (1 million atoms)	2.47E-01	8.25E-02
Vapour (8000 atoms)	2.50E+00	8.34E-01
Vapour (1 million atoms)	2.93E+00	9.78E-01

The diffusion coefficient increase from solid to gas due various degree of freedom. The estimated D from VACF have the same trend as D which are obtained from MSD.

JC: You haven't said what you think is the largest source of error. You should include the plots of the VACF that you used to get the diffusion coefficients. Can you compare the diffusion coefficients calculated from the MSD and from the VACF, which method is more reliable?