Jump to content

Talk:Mod:paxtonchia250893

From ChemWiki

JC: General comments: Good report with clearly written answers to all tasks. You could add a bit more detail to some of your discussions of your own results though and use the background theory to the experiment to support your explanations.

Introduction to molecular dynamics simulation

Velocity Verlet Algorithm

By using the data and formulas given in HO.xls, the three columns "ANALYTICAL", "ERROR", and "ENERGY" were completed. ANALYTICAL was filled in using the formula for the classical harmonic oscillator x(t)=Acos(ωt+ϕ). ERROR was determined by the absolute difference between the ANALYTICAL and the velocity-Verlet solution while ENERGY is calculated by adding up both potential and kinetic energy of the harmonic oscillator. The formula for ENERGY is given by 12kx2+12mv2. The graphs of the three variables with respect to time are shown below.


Graph of ANALYTICAL against Time
Graph of ENERGY against Time
Graph of ERROR against Time


For the default timestep, the positions of the maxima are estimated at t = 2.00, 4.90, 8.00, 11.10 and 14.20. These points are plotted against time and a linear function is fitted to the data as shown below.


Graph of ERROR MAXIMA against Time with fitted linear function


Different timesteps are experimented and the maximum change in energy is plotted as a percentage against time.


Graph of Maximum % Change in Energy against Timestep


As shown in the above graph, as the timestep increases, the change in energy also increases. Hence, a small timestep (≤0.200) is required to ensure that the total energy does not change by more than 1% over the course of the simulation. It is important to monitor the total energy of a physical system when modelling its behaviour numerically to ensure that the error in calculating the properties of the system are not large.

JC: Good thorough analysis of the timestep, why does the error oscillate? Energy should be conserved so we need to make sure this is approximately true in our simulations.

Single Lennard-Jones Interaction

Working with a single Lennard-Jones interaction, ϕ(r)=0 was substituted into the equation ϕ(r)=4ϵ(σ12r12σ6r6) to find the separation, r0, at which the potential energy is zero. The results obtained are as follows:


4ϵ(σ12r12σ6r6)=0


σ12r12σ6r6=0


σ12σ6r6=0


r6=σ6


r0=σ


The force at this separation, r0, can be calculated using the formula F=dϕ(r)dr.


dϕ(r)dr=4ϵ(12σ12r13+6σ6r7)


Substitute r=σ

F(r0) = dϕ(r)dr
= 4ϵ(12σ+6σ)
= 24ϵσ


The equilibrium separation, req, occurs when the force, F, is zero.


4ϵ(12σ12r13+6σ6r7)=0


12σ12r13+6σ6r7=0


12σ12+6σ6r6=0


r6=2σ6


req=26σ


The well-depth, ϕ(req), is calculated as follows:


ϕ(req) = 4ϵ(σ12(26σ)12σ6(26σ)6)
= 4ϵ(1412)
= ϵ


When σ=ϵ=1.0, ϕ(r)=4(1r121r6).


2σϕ(r)dr = 24(1r121r6)dr
= 4[111r11+15r5]2
= 4(111(2)1115(2)5)
= 0.0248

2.5σϕ(r)dr = 2.54(1r121r6)dr
= 4[111r11+15r5]2.5
= 4(111(2.5)1115(2.5)5)
= 0.00818

3σϕ(r)dr = 34(1r121r6)dr
= 4[111r11+15r5]3
= 4(111(3)1115(3)5)
= 0.00329

JC: All maths correct and well laid out.

Periodic Boundary Conditions

To estimate the number of water molecules given a volume under standard conditions, the following equation can be used:


Numberofwatermolecules = Volume×DensityMolarmass

Numberofwatermoleculesin1mLofwater = 1×118
= 0.0556mol
= 3.35×1022molecules


Similarly, the equation can be rearranged to find the volume of water given the number of molecules under standard conditions.


Volume = Numberofwatermolecules×MolarmassDensity

Volumeof10000watermolecules = 10000×181×6.022×1023
= 2.99×1019mL


Hence, realistic volumes of liquid cannot be simulated and periodic boundary conditions are used. This can be illustrated in the following example. 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). The atom's end point before applying the periodic boundary conditions can be found by adding the vector to the atom's original position, which gives us (1.2,1.1,0.7). Since the box has the dimensions (1,1,1), any value in the end position that is larger than 1 will have 1 subtracted from it after applying the periodic boundary conditions, resulting in (0.2,0.1,0.7) as the end point.


Reduced Units

By rearranging the equations given,

r = r*×σ
= 3.2×0.34
= 1.088nm

ϵkB = 120K
ϵ = 120×1.380×1023
= 1.656×1021J
= 0.997kJmol1

T = T*ϵkB
= 1.5×120
= 180K

JC: All calculations correct.

Equilibration

Simulation Box

Generating a random position for each atom to create a simulation box can cause problems in simulations if the atoms are generated close together. This results in a large repulsion between the atoms and may cause the atoms to move too far and out of range in a single timestep.

JC: This can make the simulation crash.

With a lattice spacing in x,y,z = 1.07722 1.07722 1.07722, the volume of the unit cell is 1.077223 = 1.25. Since there is one lattice point per unit cell in a simple cubic lattice, the number of lattice points per unit volume (number density) corresponds to 11.25=0.8.

For a face-centred cubic lattice, there are 4 lattice points per unit cell. With a lattice point density of 1.2, the volume of the unit cell is 41.2=3.33. Hence, the length of the unit cell is 3.333=1.49380.

If a face-centred cubic lattice was used instead, 4000 atoms will be created as there are 4 lattice points per unit cell and the box contains 1000 unit cells.

JC: Correct.

Properties of Atoms

The purpose of the following commands in the input script are found from the LAMMPS manual[1] as follows:

"mass 1 1.0" sets the mass for all atoms of type 1 to 1.0.

"pair_style lj/cut 3.0" tells LAMMPS to use the standard 12/6 Lennard-Jones potential between the atoms with a cutoff of 3.0 when calculating the forces. This means that the interaction between a pair of atoms is only calculated if their separation is less than 3.0.

"pair_coeff * * 1.0 1.0" sets both the parameters of the Lennard-Jones potential, ϵ and σ, to 1.0 for atom pairs of any type.

If we are specifying the initial position and velocity, we are going to use the velocity verlet algorithm.

JC: Good, why is a cutoff used for this potential?

Running the Simulation

The purpose of assigning a value to a text variable (e.g. timestep or n_steps) is for our own convenience. If we need to change the value of the variable, we only need to change the value assigned to the text variable. This only requires us to make a single edit instead of having to change the value every time it is used in the code.


Checking Equilibration

The thermodynamics output of the first simulations were analysed to see how long it takes to reach an equilibrium state. Graphs of energy, temperature and pressure against time for the 0.001 timestep experiment are shown below:



As shown above, all three thermodynamic output reached a constant average within the simulation time, which meant that an equilibrium state was reached. From the graph of energy against time, the system reached equilibrium at around t = 0.4.

A single graph of energy against time for all timesteps is plotted as shown:


Graph of Energy against Time for all timesteps


A timestep of 0.0025 is the largest to give acceptable results as larger timesteps give energies that are higher than that using a timestep of 0.001. Using a timestep of 0.015 is particularly bad because the system did not reach an equilibrium state.

JC: Good choice of timestep, the average total energy shouldn't depend on the choice of timestep so all timesteps above 0.0025 are no good. Choose the largest of the acceptable timesteps - 0.0025 - to allow the simulation to cover as much time as possible.

Running simulations under specific conditions

Thermostats and Barostats

Temperature is controlled by multiplying every velocity by a constant factor γ, which can be evaluated as follows, starting from the given equation:


12imi(γvi)2=32NkB𝔗


γ2(12imivi2)=32NkB𝔗


γ2(32NkBT)=32NkB𝔗


γ2=𝔗T


γ=𝔗T

JC: Correct.

Examining the Input Script

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

For the above line in the input script, "100" means that input values are used every 100 timesteps. "1000" means that 1000 input values are used to calculate the averages. "100000" means that averages are calculated every 100000 timesteps.[1] Values of each variable e.g. temperature will be sampled every 100 timesteps and 1000 measurements contribute to the average. The simulation will run for 100000 timesteps. With a timestep of 0.0025, 250 time units will be simulated.


Plotting the Equations of State

10 constant pressure/temperature simulations were ran with 5 different temperatures (T = 2, 3, 4, 5, 6) each at 2 different pressures (P = 2, 3). The value of the timestep chosen was 0.0025. Graphs of density against temperature were plotted separately for the 2 pressures. Error bars in both x and y directions and a line corresponding to the density predicted by the ideal gas law were included. As ϵ and σ are both 1.0, the ideal gas equation expressed in reduced units is (NV)=PT.


Graph of Density against Temperature at different Pressures
Pressure = 2
Pressure = 3


The simulated density is lower than the ideal gas density for both pressures. This is due to the approximation used in the ideal gas law, where there are no repulsive forces between the particles. However, in the simulation, the atoms are modeled using the Lennard-Jones potential, which accounts for the repulsion between atoms. This causes the atoms in the simulation to be further apart from one another when compared to the ideal gas under the same conditions. Hence, the density of the simulated liquid is lower than the ideal gas.

As pressure increases, the discrepancy increases. This is due to higher pressures causing the atoms to be closer to one another. Hence, the repulsion forces between the atoms in the simulation increases, resulting in larger deviations from the ideal gas.

JC: Good, plot results on a single graph to show the trend with pressure more clearly, what about the trend with temperature? The ideal gas is a good approximation to low density gases, where interparticle interactions are less significant.

Calculating heat capacities using statistical physics

10 constant volume/temperature simulations were ran with 5 different temperatures (T = 2.0, 2.2, 2.4, 2.6, 2.8) each at 2 different densities (0.2 and 0.8). A graph of CV/V was plotted against T.


Graph of CV/V against T


The expected trend is that heat capacity decreases with temperature, which is the trend shown in the graph above. Heat capacity measures the amount of energy required to increase the temperature of the system. As temperature increases, the vibrational energy levels of the system are more accessible and lesser energy will be required to increase the temperature of the system. Hence, heat capacity decreases as temperature increases.

Furthermore, a larger density means that there are more atoms per unit volume. Since heat capacity is an extensive property, CV/V is higher for the simulation with larger density. However, one would expect CV/V to be proportional to density but the heat capacity for ρ*=0.8 is not 4 times as large as that for ρ*=0.2.

JC: You are simulating a system of spherical particles so there are no vibrational energy levels, but the heat capacity can be related to the density of energy levels (density of states since the system is classical).

An example of my input script is included as follows:

### SPECIFY THE REQUIRED THERMODYNAMIC STATE ###
variable d equal 0.2
variable T equal 2.0
variable timestep equal 0.0025

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

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

###SWITCH OFF THERMOSTAT###
unfix nvt
fix nve all nve

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

variable avedens equal f_aves[1]
variable avetemp equal f_aves[2]
variable avepress equal f_aves[3]
variable aveetotal equal f_aves[4]
variable aveetotal2 equal f_aves[5]
variable heatcap equal ${atoms}*${atoms}*(f_aves[5]-f_aves[4]*f_aves[4])/f_aves[6]

print "Averages"
print "--------"
print "Number of atoms: ${atoms}"
print "Volume: ${vol}"
print "Energy: ${aveetotal}"
print "Energy2: ${aveetotal2}"
print "Heat Capacity: ${heatcap}"
print "Density: ${avedens}"
print "Temperature: ${avetemp}"
print "Pressure: ${avepress}"

Structural properties and the radial distribution function

The simulations of the Lennard-Jones system were ran for three phases using the values shown in the table below[2] and the RDFs for the three systems are plotted.

Phase Density Temperature
Vapour 0.1 1.2
Liquid 0.8 1.2
Solid 1.2 1.2
RDFs for vapour, liquid and solid phases


The solid phase RDF has the most number of peaks and the peaks are present even at long distances. The liquid phase RDF has a few peaks that are only present at short distances. The vapour phase RDF has only one broad peak and decays the fastest to a constant of one. This is due to the solid phase having long range order where atoms are arranged in a well-defined manner for long distances while the liquid phase only has short range order where nearby atoms are ordered. The vapour phase is the most disordered with no apparent regular arrangement of the atoms.


Face-centered cubic illustration of nearest neighbours


The first three peaks of the solid phase RDF corresponds to the position of the nearest neighbour, second nearest neighbour and third nearest neighbour respectively. The diagram above shows the positions of these three nearest neighbour with respect to the lattice site in yellow. The lattice site in red is the nearest neighbour, while the lattice site in blue is the second nearest neighbour and the lattice site in green is the third nearest neighbour.


Running Integral of RDF for solid phase


The lattice spacing is the distance of the second nearest neighbour, which is 1.475 as shown in the graph above. This also corresponds to the theoretical value calculated above (1.49380) for a FCC lattice with a lattice point density of 1.2. The coordination number for each of the three peaks can be obtained from the integrals of the peaks as shown above. The coordination number of the first peak is 12. As the graph shows a running integral, the coordination number of the second peak is 17.98 - 12.02 ≈ 6 while the coordination number of the third peak is 42.31 - 17.98 ≈ 24.

JC: Nice diagram to show which atoms are responsible for the first three peaks, do the coordination numbers that you've got from the integral of g(r) make sense based on the fcc structure? Could you have calculated the lattice parameter from each of the first three peaks separately, by expressing the distances to these atoms in terms of the lattice parameter using the geometry of an fcc lattice, and then calculated an average rather than just giving the value from the second peak?

Dynamical properties and the diffusion coefficient

Mean Squared Displacement

Mean squared displacement of solid, liquid and gas simulations as a function of timestep are plotted. The diffusion coefficient is estimated from the gradient of the total MSD against timestep graph.


Graph of Total MSD against Timestep for small systems Graph of Total MSD against Timestep for large systems
Solid Simulation. D0
Solid Simulation. D0
Liquid Simulation. D=16(0.0010190.002)=0.0849
Liquid Simulation. D=16(0.0010470.002)=0.0873
Gas Simulation. D=16(0.020060.002)=1.67
Gas Simulation. D=16(0.0369970.002)=3.08


The trend in the MSD graphs are as expected. The total MSD for the solid simulation stays approximately constant at 0.02 as the solid atoms are fixed in their lattice positions. This is also seen from its negligible diffusion coefficient. The total MSD for both liquid and gas increases linearly with time. However, the total MSD for gas increases much faster than that for liquid, which is also reflected in the larger diffusion coefficient. This is due to the gaseous atoms having weaker intermolecular forces between each other and being able to move more freely than liquid atoms.

For the large system, the trends in the MSD graphs are similar to that for a small system except for the gas phase. For solid and liquid, the MSD graphs and diffusion coefficients are approximately the same in both the large and small system. However, the total MSD for solid showed lesser fluctuations in the large system. On the other hand, the total MSD for gas increases much faster in the large system than the small system, which results in a larger diffusion coefficient as well.

JC: Why do you divide by 0.002 when you calculate the diffusion coefficient? Why is the gas MSD curved initially (ballistic motion), before eventually becoming linear?

Velocity Autocorrelation Function

The normalised velocity autocorrelation function for a 1D harmonic oscillator is evaluated as followsː


x(t) = Acos(ωt+ϕ)

v(t) = dx(t)dt=Aωsin(ωt+ϕ)

C(τ) = v(t)v(t+τ)dtv2(t)dt

= A2ω2sin(ωt+ϕ)sin(ωt+ωτ+ϕ)dtA2ω2sin2(ωt+ϕ)dt

= sin(ωt+ϕ)[sin(ωt+ϕ)cosωτ+cos(ωt+ϕ)sinωτ]dtsin2(ωt+ϕ)dt

= sin2(ωt+ϕ)cosωτ+sin(ωt+ϕ)cos(ωt+ϕ)sinωτdtsin2(ωt+ϕ)dt

= cosωτsin2(ωt+ϕ)dtsin2(ωt+ϕ)dt+sin(ωt+ϕ)cos(ωt+ϕ)sinωτdtsin2(ωt+ϕ)dt

= cosωτ+sinωτsin(ωt+ϕ)cos(ωt+ϕ)dtsin2(ωt+ϕ)dt

= cosωτ+sinωτ[sin2(ωt+ϕ)]2ωsin2(ωt+ϕ)dt

= cosωτ

JC: Correct answer, but the integration in the last step of your derivation is not correct. The function sin(x)*cos(x) is an odd function (odd*even=odd) and so the integral of this function with limits which are symmetric about zero will be zero.

The velocity autocorrelation function of the 1D harmonic oscillator, solid simulation and liquid simulation were plotted as shown below.


Graph of VACF against timestep


The minima in the VACFs for the liquid and solid system represent a reverse in the velocities of the atoms. In the solid system, the atoms vibrate about their fixed lattice positions and the interactions between neighbouring atoms disrupt this perfect oscillatory motion. Hence, the velocities of the solid system is not perfectly correlated to its initial velocity and the VACF for the solid system follows a damped oscillation. In the liquid system, the liquid atoms are free to move around so any oscillatory motion is rapidly disrupted by the diffusion of the liquid atoms. This results in the VACF for the liquid system to decay much faster to zero than the VACF for the solid system.[3] The harmonic oscillator VACF represents a perfect vibration where there are no interactions to disrupt the oscillatory motion. Its velocity reverses at the end of each oscillation and it does not decay after each oscillation. Therefore, the velocity is always correlated to its initial velocity and its VACF does not decay to zero, unlike the Lennard-Jones solid and liquid.

JC: Good, collisions between particles randomise particle velocities and cause the VACF to decay. The liquid does have a small minima before it becomes decorrelated.

The trapezium rule was used to approximate the integral under the VACF for solid, liquid and gas. The graph of running integral against time is plotted as shown and the diffusion coefficient is estimated using the values of the integral at t = 10.


Graph of Running Integral against Time for small systems Graph of Running Integral against Time for large systems
Solid Simulation. D0
Solid Simulation. D0
Liquid Simulation. D=13×0.293667=0.0979
Liquid Simulation. D=13×0.270274=0.0901
Gas Simulation. D=13×5.230557=1.74
Gas Simulation. D=13×9.805397=3.27


The diffusion coefficient of the solid system is approximately zero and increases from liquid to gas. Furthermore, using both VACF and MSD to calculate diffusion coefficients yield similar results, which is expected.

The largest source of error is probably not simulating a large enough amount of time for C(τ) to reach zero. The simulation is only ran until τ=10 and it can be seen from the running integral graphs that the integrals for the liquid and gas systems are still increasing. This means that the diffusion coefficients should be larger than my estimates.

JC: Good, using the trapezium rule for the integration also introduces some error.

References

  1. 1.0 1.1 LAMMPS Manual, http://lammps.sandia.gov/doc/Section_commands.html#cmd_5. Accessed on 17 Feb 2017.
  2. J. Hansen and L. Verlet, Phase Transitions of the Lennard-Jones System, Phys. Rev., 1969, 184, 151.DOI:10.1103/PhysRev.184.151
  3. The Velocity Autocorrelation Function, http://www.compsoc.man.ac.uk/~lucky/Democritus/Theory/vaf.html. Accessed 22 Feb 2017.