Talk:Modxl9814

JC: General comments: All tasks answered, and results look good. Make sure that you understand the background theory behind these questions and use it to explain your result, and that your explanations don't contradict what your results show.

Theory

Velocity Verlet Algorithm

One way to solve Newton's Second law F=ma is the velocity-Verlet algorithm. By using a Taylor expansion,the atomic positions, velocities and accelerations can be approximated at time t with good precision. The position of atom i, at time t, is denoted by $x_i(t)$ and the velocity of the atom at time t is denoted by $v_i(t)$. Position at the next timestep $t+\delta t$ can be expressed by:

$x_i(t+\delta t)=x_i(t)+\frac{d{x}_i(t)}{dt}\delta t+\frac{1}{2!}\frac{d^2{x}_i(t)}{d{t}^2}\delta t^2+\frac{1}{3!}\frac{d^3{x}_i(t)}{d{t}^3}\delta t^3+O(\delta t^4)(1)$

A single timestep is expressed as:

$x(t+\delta t)=x_t+v_t\delta t+\frac{1}{2}{a}_t\delta t^2(2)$ $v(t+\delta t)=x_t+\frac{1}{2}(a_{t+\delta t}+a_t)\delta t(3)$

Fig 1: Classically calculated positions vs. velocity verlet calculated positions

The classical harmonic oscillator can be describe by $x(t)=Acos(\omega t+\varphi )$. The errors oscillate through 5 peaks in the simulated time. The plot of the total energy vs. time of the simulated system:

Fig 2:Error vs. time

The cumulative error over a constant interval of time is given by:

$error(x(t_0+n\delta t))=O(\delta t^2)$

Thus, it can be seen from this equation that the relation between the maxima of the error of the Velocity-Verlet algorithm and $\delta t$ is quadratically increasing. The graph of the maxima of error vs. time therefore can be fit into the quadratic equation in figure 2.

JC: Good additional research, but a linear fit would have been sufficient here.

The total energy of the oscillating system is the sum of the kinetic energy and the potential energy, with $E_k=\frac{1}{2}m{v}^2$ and $E_p=\frac{1}{2}k{x}^2$. In this case, m=1 and k=1, therefore the equation is $E=\frac{v^2+x^2}{2}$.

Fig 3a:Energy vs. Time at 0.1 timestep with error limites of 0.5% on either side

Fig 3b:Energy vs. Time at 0.2 timestep with error limites of 0.5% on either side

In order for the total energy not to change by more than 1% over the course of the simulation, the timestep needs to be 0.2. It is important to monitor the total energy of the system to ensure that energy conservation is obeyed, the same as the real system.

JC: Good choice of timestep.

Atomic Forces

The Lennard-Jones potential can tell the potential energy of the interaction between two uncharged atoms. It can be expressed in (12,6) form:

$\varphi (r)=4\epsilon (\frac{{\sigma }^{12}}{r^{12}}-\frac{{\sigma }^6}{r^6})$

In this equation,$\epsilon$ is the potential well depth,$\sigma$ is the distance where the potential between the pair of particles is zero and r is the distance between the pair of particles.

As force is the negative derivative of potential energy, the equation of force in terms of the Lennard-Jones potential is:

$F_i=-\frac{d\varphi (r^N)}{d{r}_i}$=$24\epsilon [2(\frac{{\sigma }^{12}}{r_i^{13}})-\frac{{\sigma }^6}{r_i^7}]$

When the potential energy is zero, $r_i=\sigma =r_0$, therefore by substitution we can get:

$F=24\epsilon [2(\frac{r_0^{12}}{r_0^{13}})-\frac{r_0^6}{r_0^7}]=24\epsilon [\frac{2}{r_0}-\frac{1}{r_0}]=\frac{24\epsilon }{r_0}$

The equilibrium is reached when the resultant force is zero, therefore:

$F=24\epsilon [2(\frac{{\sigma }^{12}}{r^{13}})-\frac{{\sigma }^6}{r^7}]=0$ $2(\frac{{\sigma }^{12}}{r^{13}})-\frac{{\sigma }^6}{r^7}=0$

Divide both sides by $\frac{{\sigma }^6}{r^7}$:

$2\frac{{\sigma }^6}{r^6}-1=0$

Therefore, the equilibrium separation is:

$r_{eq}=\sqrt[6]{2}\sigma$

The LJpotential at $r_{eq}$ is:

$\varphi (r_{eq})=4\epsilon (\frac{{\sigma }^{12}}{4{\sigma }^{12}}-\frac{{\sigma }^6}{2{\sigma }^6})=4\epsilon (-\frac{1}{4})=-\epsilon$ $\epsilon =-\varphi (r_{eq})$ $\int \varphi \left(r\right)\mathrm{d}r=-\frac{4}{11}\epsilon {\sigma }^{12}{r}^{-11}+\frac{4}{5}\epsilon {\sigma }^6{r}^{-5}$

$\sigma =\epsilon =1.0$, therefore:

$\int \varphi \left(r\right)\mathrm{d}r=-\frac{4}{11}{r}^{-11}+\frac{4}{5}{r}^{-5}$

$\underset{2\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=\frac{4}{11}\times \frac{1}{2^{11}}-\frac{4}{5}\times \frac{1}{2^5}=-2.48\times 10^{-2}$

$\underset{2.5\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=\frac{4}{11}\times \frac{1}{2.5^{11}}-\frac{4}{5}\times \frac{1}{2.5^5}=-8.18\times 10^{-3}$

$\underset{3\sigma }{\overset{\mathrm{\infty }}{\int }}\varphi \left(r\right)\mathrm{d}r=\frac{4}{11}\times \frac{1}{3^{11}}-\frac{4}{5}\times \frac{1}{3^5}=-3.29\times 10^{-3}$

JC: All maths correct.

Periodic Boundary Conditions

$1mL=1c{m}^3$. The density of water=$1g/c{m}^3$ under standard consitions (298K, 1atm). So the total mass of 1 mL water= 1g. The number of moles of water molecules=$\frac{1}{M_{H_{2}O}}=\frac{1g}{18g/mol}=0.056moles$

Therefore the total number of molecules in 1 mL of water=$n\times N_a=0.056\times 6.02\times 10^{23}=3.37\times 10^{22}$

The volume of 10,000 molecules of water=$\frac{10000}{3.37\times 10^{22}}=2.97\times 10^{-19}$

The initial position of atom is $(0.5,0.5,0.5)$. After it moves along the vector $(0.7,0.6,0.2)$, the position should be $(1.2,1.1,0.7)$. Applying the periodic boundary of $(0,0,0)$ to $(1,1,1)$, the position should be $(0.2,0.1,0.7)$.

Reduced Units

$r=r^{*}\sigma =0.34\times 10^{-9}\times 3.2=1.09nm$

The well depth=$\epsilon =120K\times K_B\times 10^{-3}\times 6.022\times 10^{23}=0.997KJ/mol$

$T^{*}=1.5$, therefore $T=T^{*}\times \frac{\epsilon }{K_B}=1.5\times 120K=180K$

JC: All calculations correct.

Equilibration

Creating the simulation box

Giving atoms random starting coordinates may make two atoms generated too close together. This will cause the two atoms to collide and arise the repulsion between the two atoms. The repulsive force between the atoms will drive them apart, leading to increase in the potential energy of the system and making it very unstable.

JC: High repulsive forces can cause the simulation to crash.

A face-centered cubic lattice has 4 lattice points per unit cell. The side length of the cubic unit cell=$\sqrt[3]{\frac{4}{1.2}}=1.49$

If 1000 unit cells were generated by the create_atoms command, 4000 atoms would be generated for a FCC lattice.

Setting the properties of the atoms

Mass 1 1.0

This means the mass of the single type of atom is 1.0.

Pair_style lj/cut 3.0

"Pair_style" indicates that the interaction is pairwise interaction. "lj.cut" describes the standard 12/6 Lennard-Jones potential without a Coulombic pairwise interaction. "3.0" indicates that the global cutoff for atoms is at 3.0.

Pair_coeff * *1.0 1.0

"pair_coeff" specifies the pairwise force field coefficients. The two asterisks indicate that the command will apply to all atoms.

If $x_i$ and $v_i$are specified,the Velocity-Verlet algorithm will be used for this simulation.

JC: Why do we use a cutoff? What are the forcefield coefficients in your simulations, considering that you are using a Lennard-Jones potential?

Running the simulation

The purpose of defining variable is that we don't need to manually change the numerical timestep each time the timestep needs to be changed. This reduces the human errors that may occur as the timestep only needs to be changed once to the value defined.

JC: Using variables means that all commands in the simulation that depend on a certain value, like the timestep in this case, are updated automatically when that value is changed.

Checking equilibration

Fig 4a:Total energy vs. time at 0.001 timestep

Fig 4b: Temperature vs Time at 0.001 timestep

Fig 4c: Pressure vs. Time at 0.001 timestep

The simulation reaches equilibrium at 0.001 timestep as pressure and temperature become constant. It can be seen from pressure and temperature data that the simulation reaches equilibrium at t=0.29.The average pressure value is about 2.61 and the average temperature value is about 1.26.

Fig.5: Graph of energies for all timesteps

It can be seen from Fig 5 that the total energy produced by 0.0025 timestep are very close to those produced by 0.001 timestep. Simulations at 0.0075 and 0.01 also reach equilibrium but the total energies are higher than those produced by 0.001 timestep, thus these timesteps are not very accurate. Therefore the largest timestep to get acceptable results is 0.0025 and the worst choice is 0.015 timestep as the simulation doesn't reach equilibrium.

JC: Good choice of timestep, average energy should not depend on the timestep, making timesteps above 0.0025 poor choices.

Running simulations under specific conditions

Barostat and Thermostat

In the system with N atoms, each with 3 degrees of freedom:

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

By multiplying every velocity by $\gamma$ and substituting T with $\mathfrak{𝔗}$ we can get the second equation:

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

By substituting (2) we can get:

$\frac{3}{2}N{k}_{B}T\times {\gamma }^2=\frac{3}{2}N{k}_B\mathfrak{𝔗}$ ${\gamma }^2=\frac{\mathfrak{𝔗}}{T}$ $\gamma =\sqrt{\frac{\mathfrak{𝔗}}{T}}$

JC: Correct.

Examining the Input Script

The numbers 100 1000 100000 indicate the timesteps the input values will be used to compute the averages of density, pressure and temperature. For this simulation, the average will be calculated using values produced by timestep 100,200,...100000. Therefore, 1000 values will be used to calculate the average. The following line tells LAMMPS to run the simulation for 100000 timesteps. 0.0025 timestep will be used. Therefore 250 time units are simulated.

JC: Be more specific about what each of the three numbers refers to.

Plotting the Equations of State

Fig.6:Density vs Temperature and Ideal Gas law at p=2.3 and p=2.6

Fig.7:Density calculated by Ideal Gas Law compared to LJ model at P=2.3 and P=2.6

Simulations were conducted at temperatures 2,2.5,3,3.5,5 and pressures 2.3 and 2.6.

Density can also be calculated by Ideal Gas Law $\frac{N}{V}=\frac{P}{k_{B}T}$ through the following:

${\rho }^{*}=\frac{N}{V^{*}}=N\frac{{\sigma }^3}{V}$

As $\frac{N}{V}=\frac{P}{k_{B}T}$,$P=P^{*}\frac{\epsilon }{{\sigma }^3}$ and $T=T^{*}\frac{\epsilon }{K_B}$, by substitution we can get:

${\rho }^{*}={\sigma }^3\frac{P}{k_{B}T}={\sigma }^3\frac{P^{*}\frac{\epsilon }{{\sigma }^3}}{k_B{T}^{*}\frac{\epsilon }{K_B}}=\frac{P^{*}}{T^{*}}$

Fig.7 shows that the simulated density is much lower than the density obtained by the Ideal Gas Law. This is because the Ideal Gas Law assumes that the molecules do not interact with each other and the repulsive force between the molecules is zero.This means that the particles in Ideal Gas system can be compressed to great extent, making the volume occupied very small for a given volume. Therefore the density is higher. In the Lennard-Jones model, the molecules will interact with each other and the repulsive force is greater when the distance between the molecules is smaller. Therefore for a given volume the molecules will rather stay far apart and the density is lower .

It can be seen from Fig.7 that the discrepancy increases with pressure. This is because at lower pressure, provided that the volume is large enough, the intermolecular distance is larger and the density will not change a lot by the distance between the particles.

JC: Good, how does the difference between the ideal gas and siulation results change with temperature?

Heat Capacity Calculation

In the NVT ensemble, pressures (0.2,0.8) and temperatures (2,2.2,2.4,2.6,2.8) were used to calculate the heat capacity by using the following equation:

$C_V=\frac{\partial E}{\partial T}=\frac{\mathrm{V}\mathrm{a}\mathrm{r}\left[E\right]}{k_B{T}^2}=N^2\frac{⟨E^2⟩-{⟨E⟩}^2}{k_B{T}^2}$

The code to run the simulation in the NVT ensemble is as following:

variable density equal 0.2

### DEFINE SIMULATION BOX GEOMETRY ###
lattice sc ${density}
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.0
variable p equal 4
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 ###
thermo_style custom step etotal temp press atoms density vol
variable dens equal density
variable dens2 equal density*density
variable temp equal temp
variable volume equal vol
variable temp2 equal temp*temp
variable press equal press
variable press2 equal press*press
variable N2 equal atoms*atoms
variable E2 equal etotal*etotal
variable E equal etotal
fix aves all ave/time 100 1000 100000 v_dens v_temp v_press v_dens2 v_temp2 v_press2 v_E v_E2
unfix nvt
fix nve all nve
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 heatcapacity equal ${N2}*(f_aves[8]-f_aves[7]*f_aves[7])/f_aves[5]
variable CV equal ${heatcapacity}/${volume}

print "Averages"
print "--------"
print "Density: ${avedens}"
print "Stderr: ${errdens}"
print "Temperature: ${avetemp}"
print "Stderr: ${errtemp}"
print "Pressure: ${avepress}"
print "Stderr: ${errpress}"
print "heatcapacity: ${heatcapacity}"
print "volume: ${volume}"
print "heatcapacity/volume: ${CV}"

Fig.8: Cv/V vs. temperature at densities 0.2 and 0.8

$\frac{C_V}{V}$ was plotted against temperature. The volume for $\rho =0.2$ is 16875. The volume for $\rho =0.8$ is 4218.75. The heat capacity is inversely proportional to temperature from Fig.8, the same as shown in the equation $C_V=\frac{\partial E}{\partial T}=\frac{\mathrm{V}\mathrm{a}\mathrm{r}\left[E\right]}{k_B{T}^2}=N^2\frac{⟨E^2⟩-{⟨E⟩}^2}{k_B{T}^2}$. This is because the lattice energy gap decreases with increasing temperature, so less energy will be required. This indicates that heat capacity is proportional to energy as shown in the equation. Also, it is shown that the lower the density, the lower the heat capacity.This is because high density means the particles will be closer together with lower volume, therefore less heat is required to heat the system. For the same number of particles, if the density is lower, that means the volume the particles take up is larger. Therefore the heat required is higher .

JC: The heat capacity is not necessarily inversely proportional to temperature because <E> and <E^2> are also dependent on temperature. You say that at higher density less heat is needed to raise the temperature which should mean the heat capacity is smaller, this is not what your results show.

Structural properties and the radial distribution function

Fig.9: g(r) vs. r for solid, liquid and gaseous phases

The radial distribution function was plotted for vapour, liquid and solid phases(Fig.9). The densities and temperatures were chosen from the phase diagram for the Lennard-Jones diagram.:

Phase	Density	Temperature
Vapour	0.1	1.2
Liquid	0.8	1.2
Solid	1.6	1.2

The radial distribution function indicates the probability of finding a nearest neighbor from a particle. It will reveal the phase of the system simulated. The RDFs for the three systems are very different. The solid has the largest number of peaks followed by liquid and then gas. The peaks represent the density around each atom and hence solid which has the highest density will have more peaks. The peaks in the solid phase has decreasing amplitude with increasing r. It can be seen that the probability of finding a particle between the first and second peak is zero. This is because particles in solid phase do not have brownian motion and can only vibrate in fixed positions. The solid phase has long and short range order and this can be indicated by the peaks. The short range order is shown by the first three tall peaks and the long range order is shown by the smaller peaks behind.

JC: More peaks don't indicate a higher density, the RDFs are normalised - they all tend to 1 - and so you cannot tell which has a higher density. The RDF just indicates order.

The RDF of the liquid phase has three peaks with decreasing amplitude as r increases. The wider peaks mean that the liquid phase is more disordered than the solid phase.The decreasing amplitude with increasing interatomic separation indicates that the Brownian motion of particles in the liquid phase makes the order decrease with increasing separation. The three peaks indicate that the liquid phase only has short range order.

The RDF of the gas phase only has one broad peak. This suggests the gas phase is highly disordered and there is no short nor long range order.

In the RDF of the solid phase, the first three peaks correspond to the nearest neighbor of the referenced particle, the second nearest particle and the third nearest particle respectively. The lattice spacing is the distance between the zero probability minima and is 1.275 in reduced units.

Fig.10: Integral of g(r) vs. interatomic distance for solid phase

The coordination number for the first three peaks can be calculated from the plot of the integral of g(r) against interatomic distance. The integral of g(r) at the inflection points represent the coordination number of the three nearest neighbors. As a FCC lattice is used in a solid system, there should be 12 neighboring particles around each particle (shown at r=1.275). So the coordination number of the first peak is 12. The next inflection number has a g(r) integral of 18. As it is a running integral, the coordination number of the second peak is $18-12=6$. The coordination number of the third peak is $42-18=24$.

JC: It would have been good to show on a diagram of an fcc lattice which atoms are responsible for the first 3 peaks. Did you calculate the lattice parameter?

Dynamical properties and the diffusion coefficient

Mean Squared Displacement

The mean squared displacement is given by:

$\mathrm{M}\mathrm{S}\mathrm{D}\equiv ⟨(x-x_0{)}^2⟩=\frac{1}{N}\sum _{n=1}^N(x_n(t)-x_n(0){)}^2$

For 3375 atoms:

Fig 11a: Liquid simulation at d=0.8, T=1.2

Fig 11b: Gas simulation at d=0.1, T=1.2

Fig 11c: Solid simulation at d=1.6, T=1.2

For 1 million atoms:

Fig 12a: Liquid simulation at d=0.8, T=1.2

Fig 12b: Gas simulation at d=0.1, T=1.2

Fig 12c: Solid simulation at d=1.6, T=1.2

All these plots are as expected. For the liquid phase, MSD is directly proportional to timestep, this is because the particles move in brownian motion. Foe gaseous phase, the first part (timestep 0-2000) is curved and the second part (above 2000) is linear. The curved part is because the particles move randomly in the system and the distance between them is very large. The frequency of collision between the particles is very low and thus the velocity of the atoms will be almost constant. The distance travelled per unit time is constant, thus MSD is proportional to $t^2$. As longer time is simulated, collisions will occur more frequently and the motion can be described by brownian motion and MSD changes linearly with timestep. For solid phase, the particles only vibrate in fixed positions and do not have enough kinetic energy to diffuse, thus MSD reaches at constant value at around timestep 200.

JC: The curved part of the gas plot is due to balistic motion, not because the particles move randomly. Show the lines of best fit on the graphs.

The diffusion coefficient is given by:

$D=\frac{1}{6}\frac{\delta ⟨r^2⟩}{\delta t}$

$\delta ⟨r^2⟩$ is the slope of the trendline of the mean squared displacement vs. timestep plot. The timestep $\delta t=0.002$.

For small system of 3375 atoms, the diffusion coefficient is: Liquid: $D=\frac{1}{6}\times \frac{0.001}{0.002}=0.083$; Gas: $D=\frac{1}{6}\times \frac{0.0245}{0.002}=2.042$; Solid: $D=\frac{1}{6}\times \frac{6\times 10^{-8}}{0.002}=5\times 10^{-6}$;

For large system of 1 million atoms, the diffusion coefficient is: Liquid: $D=\frac{1}{6}\times \frac{0.001}{0.002}=0.083$; Gas: $D=\frac{1}{6}\times \frac{0.0305}{0.002}=2.542$; Solid: $D=\frac{1}{6}\times \frac{6\times 10^{-8}}{0.002}=5\times 10^{-6}$.

All of the diffusion coefficients are in reduced units. The coefficients for the larger system were similar to the ones for the smaller system.

Velocity Autocorrelation Function

The equation of the position of a 1D harmonic oscillator is:

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

$v(t)=\frac{dx}{dt}$, thus:

$v(t)=\frac{d(Acos(\omega t+\varphi )}{dt}=-A\omega sin(\omega t+\varphi )$

and

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

Therefore, by substitution:

$C(\tau )=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}v(t)v(t+\tau )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{v}^2(t)\mathrm{d}t}=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}-A\omega sin(\omega t+\varphi )\times -A\omega sin(\omega (t+\tau )+\varphi )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}(-A\omega sin(\omega t+\varphi ){)}^2\mathrm{d}t}=\frac{(A\omega {)}^2\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}sin(\omega t+\varphi )sin(\omega (t+\tau )+\varphi )\mathrm{d}t}{(A\omega {)}^2\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^2(\omega t+\varphi )\mathrm{d}t}=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}sin(\omega t+\varphi )sin(\omega (t+\tau )+\varphi )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^2(\omega t+\varphi )\mathrm{d}t}$

Since $sin(x+y)=sin(x)cos(y)+cos(x)sin(y)$:

$C(\tau )=\frac{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}sin(\omega t+\varphi )[sin(\omega t+\varphi )cos(\omega \tau )+cos(\omega t+\varphi )sin(\omega \tau )]\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^2(\omega t+\varphi )\mathrm{d}t}=\frac{cos(\omega \tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^2(\omega t+\varphi )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^2(\omega t+\varphi )\mathrm{d}t}+\frac{sin(\omega \tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}sin(\omega t+\varphi )cos(\omega t+\varphi )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^2(\omega t+\varphi )\mathrm{d}t}=cos(\omega \tau )+\frac{sin(\omega \tau )\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}sin(\omega t+\varphi )cos(\omega t+\varphi )\mathrm{d}t}{\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}si{n}^2(\omega t+\varphi )\mathrm{d}t}$

Since $sin(2x)=2sin(x)cos(x)$:

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

As $sin(x)$ is an odd function, the area above the x-axis and below the x-axis cancel out from negative infinity to positive infinity. Thus, $\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}sin(2(\omega t+\varphi ))=0$. therefore:

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

JC: Nice clear derivation and well spotted that you can simplify the integrals. cos*sin is also an odd function (even x odd = odd) so you don't need to change this to sin(2x).

Fig.14: VACF for solid, liquid and 1D Harmonic Oscillator

The minima in the VACF for the liquid system represent the collisions between the atoms and the solvent molecules and change in direction. The minima in the VACF for the solid system represent the collisions between the atoms and change in direction. The minima for the solid system is lower than the minima for the liquid system because of the stronger interatomic forces. The VACF for the liquid system only has one weak oscillation, this is because the atoms only interact with their closest neighbor. In the VACF for the solid system, there are more oscillations as the atoms can vibrate in fixed positions. The harmonic oscillator VACF is very different to the Lennard Jones liquid and solid as there are no interactions between the atoms so the atoms will always vibrate with constant velocity without loss in energy. Therefore, the amplitude doe not change.

JC: There is not only one collision in the liquid, but the first collision randomises the particle velocities and the VACF decays to zero more quickly than in a solid.

By applying the trapezium rule, integral under VACF can be calculated and running integral can be plotted:

For 3375 atoms:

Fig.15a: running integral vs. time for liquid

Fig.15b: running integral vs. time for solid

Fig.15c: running integral vs. time for gas

For 1 million atoms:

Fig.16a: running integral vs. time for liquid

Fig.16b: running integral vs. time for solid

Fig.16c: running integral vs. time for gas

The diffusion coefficient is calculated by:

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

The last point of the running integral is $\underset{0}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}\tau ⟨\mathbf{v}\left(0\right)\cdot \mathbf{v}\left(\tau \right)⟩$.

For 3375 atoms: Liquid: $D=\frac{1}{3}\times 0.2937=9.79\times 10^{-2}$; Solid: $D=\frac{1}{3}\times 5.64\times 10^{-4}=1.88\times 10-4$; Gas: $D=\frac{1}{3}\times 7.054=2.351$;

For 1 million atoms: Liquid: $D=\frac{1}{3}\times 0.2703=9.01\times 10^{-2}$; Solid: $D=\frac{1}{3}\times 1.37\times 10^{-4}=4.57\times 10^{-5}$; Gas: $D=\frac{1}{3}\times 9.805=3.268$.

The diffusion coefficient calculated from this method was largest for gas, followed by liquid and then gas. The coefficients for the larger system were very similar to the ones for the smaller system. The coefficients calculated by MSD were similar to the ones calculated by VACF for liquid and gas, but the coefficient calculated by VACF was larger than the one calculated by MSD for solid. The largest source of error may be that the trapezium rule overestimates the area under the solid curve as the timestep is not small enough.

JC: The trapezium rule is a source of error, but does it definitely overestimate the integral?