Jump to content

Rep:Mod:IOPJKL

From ChemWiki

Liquid Simulations Computational Experiment

Theory

Task 1: Calculations of a Harmonic Oscillator using both Classical and Velocity-Verlet solutions

Figure 1.a) Positions calculated via the two solutions vs time b) Difference between the positions calculated via different methods against time c) Total energy of the harmonic oscillator calculated using Velocity-Verlet solutions vs time

Position calculations from the two methods: The position x at time t was determined using x(t)=Acos(ωt+ϕ) where; A=1,ω=1 and ϕ=0. These were plotted against time(Figure 1.a).


Error analysis of the two positions: The difference between the two positions was calculated and plotted against time(Figure 1.b). The error increases with time and position, this loss of precision could come from two sources. As the position gets further from the original position the a(t)dt2 part of the equation has less importance also the errors could build up since the previous position is used in the preceding position calculation.


Total Energy of the Harmonic-Oscillator vs time: Values for velocity v and position x were obtained from the Velocity-Verlet calculation, values for the force-constant k and mass m were both given as 1. Using E=12mv2+12kx2 the total energies were calculated and plotted against time(Figure 1.c).


Task 2: Plot of maxima error positions vs time

The coordinates of the 5 maximas in figure 2 were estimated and plotted on a new graph. This yielded a linear function where

x(t)=0.0004x7*105

.

Figure 2. Linear plot of the maxima error positions(figure 1.b) against time


Task 3: Timesteps to minimise variation in total energy

Multiple timesteps were used to calculate the total energy of the harmonic oscillator. The percentage difference was taken between the initial energy assumed to be correct and the energy of the first minima, these results were then plotted(figure 3). In order to keep the change in total energy below 1% a timestep equal to or below 0.2 should be used. The change in total energy of the simulation should be kept to a minimum as the total energy should not vary over time in the harmonic oscillator.

Figure 3. Plot of percentage energy difference and timestep


Task 4: Lennard-Jones interaction

1. Calculation of distance at zero-potential, r = r0: ϕ(r0)=4ϵ(σ12r012σ6r06)

where φ = potential energy r = the distance ε = well depth and σ = distance from zero potential distance
                                ϕ(r0)=4ϵ(σ12σ6r06)
=4ϵ(1r06σ6)=0
=1r06σ6=0
σ6=r06
σ=r0

2. Force at r0:

F0=dϕ(rN)dr0=ddr04ϵ(σ12r012σ6r06)=4ϵ(12σ12r013+6σ6r07)

Since σ=r0

=4ϵ(12σ12σ13+6σ6σ7)=4ϵ6σ=24ϵσ


3. Finding the equilibrium separation req:

Feq=dϕ(rN)dreq=ddreq4ϵ(σ12req12σ6req6)=4ϵ(12σ12req13+6σ6req7)=48ϵ(σ12req13)+24ϵ(σ6req7)=0
48ϵ(σ12req13)=24ϵ(σ6req7)
2ϵ(σ12req13)=ϵ(σ6req7)
2ϵ(σ12)=ϵ(σ6req6)
2ϵ(σ6)=ϵ(req6)
2(σ6)=(req6)
26σ=req

4. Finding the well-depth at equilibrium:

ϕ(req)=4ϵ(σ12req12σ6req6)=4ϵ(σ12(26σ)12σ6(26σ)6)=4ϵ(σ1222σ12σ62σ6)=4ϵ(1412)=4ϵ(14)=ϵ

5. Evaluation of integrals:

a) 2σϕ(r)dr

24ϵ(σ12r12σ6r6)dr=24ϵ(112r121r6)dr=24ϵ(1r121r6)dr
=[4ϵ(111r1115r5)]2
=4ϵ(111()1115()5)4ϵ(111(2)1115(2)5)
=0+4ϵ(111(2)1115(2)5)=24.82KJ/mol

b) 2.5σϕ(r)dr

2.5ϕ(r)dr=2.54ϵ(σ12r12σ6r6)dr=2.54ϵ(112r121r6)dr=2.54ϵ(1r121r6)dr
=[4ϵ(111r1115r5)]2.5
=4ϵ(111()1115()5)4ϵ(111(2.5)1115(2.5)5)
=0+4ϵ(111(2.5)1115(2.5)5)=8.18KJ/mol

c) 3σϕ(r)dr

3ϕ(r)dr=34ϵ(σ12r12σ6r6)dr=34ϵ(112r121r6)dr=34ϵ(1r121r6)dr
[4ϵ(111r1115r5)]3
4ϵ(111()1115()5)4ϵ(111(3)1115(3)5)
=0+4ϵ(111(3)1115(3)5)=3.29KJ/mol


Task 5: Estimating the number of molecules in volumes of H2O

1. Number of molecules in 1mL:

Mr(H2O)=18.02g/mol,ρ=1g/cm3

n=mMr=118.02=5.55x102 Numberofmolecules=5.55x102x6.022x1023=3.34192x1022

2. Volume of 1000 molecules:

n=10006.022x1023=1.66x1021mol m=1.66x1021x18.02=2.99x1020g V=2.99x10201=2.99x1020cm3


Task 6: Periodic boundary conditions

When the atom at (0.5,0.5,0.5) in the cubic simulation box (0,0,0,) to (1,1,1) moves along the vector (0.7,0.6,0.2) it ends up at (0.2,0.1,0.7) (when the periodic boundry condition is taken into consideration).


Task 7: Converting reduced values to real values:

1. Reduced distance

r*=3.2=rσ
r=3.2x0.34x109=1.09x109=1.09nm

2. Well-depth

ϵKB=120
ϵ=120KB=8.69x10246.022x1023=997.71Jmol1=9.98x101KJmol1

3. Temperature

T*=1.5=KBTϵ
T=T*ϵKB=1.5x1.657x1021KB=180K

Equilibration

Task 8: Problems with random starting coordinates in non-solid simulations

If atoms were given random starting coordinates then they could end up within a contact proximity, they could bump into each other and then they would not only experience the potential of each but also the kinetic energy.

Task 9: The number density of lattice points calculations:

The number density n is given by the following equation

n=NV where N is the number of atoms/lattice points per unit cell and V is volumes of the unit cell.

1. The number density of a simple cubic lattice with inter lattice point distance of 1.07722

For a simple cubic lattice there is one atom/lattice point, therefore;

n=NV=Nr*3=11.077223=0.8

2. Side length of a face-centred cubic lattice with a number density of 1.2

A face-centred cubic lattice has 4 atoms/lattice points, therfore;

r*=V3=Nn3=41.23=1.49380

Task 10: How many atoms are created by create_atoms 1 box command for a face-centred cubic lattice?

Since there are 4 atoms per lattice and the box is ten lattice spacings larger than the single lattice spacing then there would be 4000 atoms.

Task 11: LAMMPS commands for setting the properties of atoms:

1. Mass command

mass I value = mass 1 1.0  

This command is used for setting the mass of the atom, in the general command I is the atom type and value is the mass(in reduced units). In the example the atom is type 1 and the mass is 1.0.

2. Pair_style command

pair_style style args = pair_style lj/cut 3.0

The pair_style command is used to describe paired interactions, style is the interaction which is being looked at and args is the argument used for that style. In the example above the style is defining the cutoff distance, that is the distance at which interactions become negligible. The argument for the cutoff distance is 3.0, after which the simulation will not compute.

3. pair_coeff command

pair_coeff * * args = pair_coeff * * 1.0 1.0 This command is used to define the force field coefficients between pairs of atoms, the asterix terms are to define the coefficients for all atoms, the second part of this command, args, defines the coeffecients.

Task 12: Specifying the velocity of each atom

The Velocity-Verlet algorithm is used to specify velocity and position vi(t+δt)=vi(t+12δt)+12ai(t+δt)δt(9)

Task 13: Running the simulation:

The command $(100/timestep) command allows for the immediate evaluation of the the formula with the pre-defined variable i.e. the timestep. The following line variable n_steps equal floor(100/$({timestep}) will be broken down into individual contributions. The variable component is used to allows for the return to re-run a command, this allows creates a loop. The n_steps part describes the number of timesteps you wish to calculate which is defined in the 'Run Simulation' command. Equal is used to create a mathematical argument. Floor(11/$({timestep}) is useful in setting the maximum timestep bigger than the original one. In the run simulation part first we define the number of timesteps to be considered and then report the floor value i.e. the total time it takes. These lines are essential in creating a loop where each new timestep may be calculated subsequently. Make it easier to change the timestep, all values dependedent on the timestep will be changed with it and so it reduces possible human errors when changing all variables and also makes it easier and faster.

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

Task 14: Checking equilibration

1. Pressure, Temperature and energy against time for simulation ran at 0.001 timestep

Figure 5. Temperature, pressure and total energy against time for the simulation using a timestep of 0.001

2. Equilibration

The simulation reaches equilibration at approximately 0.28, this is the timestep at which the simulation becomes relatively constant.

3. Maximising results accuracy whilst keeping the timestep within a reasonable timescale

The best timestep for maximum accuracy and sufficient information is 0.0025. It overlaps with 0.001 and so although it is more than twice the time as 0.001 it maintains a satisfactory level of accuracy. The others give equilibrium energies which are too far off the most accurate value, with the worst being for 0.01 as it diverges quite dramatically.

Figure 5. The energy vs time plots of simulations with varying timesteps


Running simulations under specific conditions

Task 15: Finding the expression for factor γ

γ is a factor by which the velocity can be corrected.

12mivi2=32NKBT (1) 12mivi2=32NKBT (2)

Where T' is the target temperature and T is the instantaneous temperature

12mivi2=32NKBT 12mivi2=12NKBT mivi2=NKBT vi2=NKBTmi vi=NKBTmi2

Insert v_i into equation x

12miγ2NKBTmi=32NKBT 12miγ2NKBTmi=12NKBT γ2NKBT=NKBT γ2=NKBTNKBT γ2=TT γ=TT2


Task 16: The importance of number, repetition and frequency of average calculations with timesteps

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

100 sets the number of timesteps after which another thermodynamic value is used to calculate the average, 1000 is the number of thermodynamic values which are used for the average and 1000000 sets the frequency at which averages are calculated.

100000 timesteps will be run.


Task 17: Equation of state plots

5 temperatures above the critical temeperature were chosen(1.6, 1.8, 2.0, 3.0 and 5.0) these were all simulated with a pressure of 2.6(n-1) and also 3.3(n-2). The timestep chose was 0.01 as in previous experiments it was the largest timesteps yielding the most accuracy. The pressures were plotted against the temperatures along with their error bars.

The ideal gas densities were calculated with the following equation: ρ=NV=PKBT=PT where K_B = 1 in reduced terms

Figure 6. Densities n-1 and n-2 calculated via simulation and Ideal Gas methods plotted against temperature

The simulated density is lower than that of the ideal gas method, this can be explained by an assumption to satisy an ideal gas. The atoms must not have any forces upon them and also they mustn't take up any space. In the simulation the atoms do take up space and therefore they have a lower density. This discrepancy increases with pressure. The general trend of the two plots is a decrease in heat capacity as temperature is increased.


Heat Capacity Calculation

Task 18: The temperature dependence of the isochoric heat capacity

Figure 7. Cv/V as a function of T* at P= 0.2 and 0.8

A LAAMPS script which could calculate the isochoric heat capacity from calculated thermodynamic values was used figure 7. The Cv/V was calculated with the following equation; Cv=N2<E2><E>2KBT2

Figure 8. Cv/V as a function of T* at P= 0.2 and 0.8

The effect on heat capacity at different temperatures and pressures

The decrease in heat capacity at higher temperatures is a result of the non-classical behaviour of the system at low temperatures. The energy levels behave as if they are discrete, the lower the temperature the more discrete they are. As illustrated here, the energy levels are closer together at higher temperatures e.g.T2 and further apart at lower temperatures e.g.T1. When a heat increment, dq, is introduced into the two systems we see that it has a different effect in both. At T2, dq will give a greater temperature increase than T1, this causes the denominator in the equation below to become more important thus decreasing the isochoric heat capacity.

CV=dqdT


The difference between the two pressured systems in figure 7 is what would be expected, the higher pressure plot has a higher isochoric heat capacity. At a higher pressure the atoms are closer together, they can thus transfer heat more effectively with a given temperature change than a system with a lower pressure, since the atoms are further apart.

Structural properties and the radial distribution function

Task 18: The Radial Distribution Functions of a solid, liquid and gas

1. Comparisons of the solid, liquid and gas radial distributions functions as well as the integrated function

figure 9
figure 10

The distribution plot(figure 9) of the gas has one short, broad peak which is close to the reference atom; as you move further away the distribution becomes constant. This characteristic shows the small probability of finding an atom close to the reference atom; also when it becomes constant this shows how there is negligible probability of finding another atom at the distance studied. When looking at the integrated function(figure 10), it emphasizes how few atoms you will find in the distance, this is characteristic of a gas' dispersed atoms.

The liquid's distribution has a taller first peak than the gas' and then multiple peaks tapering off as you move further away; it is evident that there is a higher probability of finding an atom throughout the space than in a gas. Additionally the integrated distribution(figure 10) shows a larger number of atoms as you move through the space. This demonstrates the higher proximity of the atoms in the liquid state.

Finally the solid's distribution has a very sharp and tall peak close to the reference atom but then as you move further away the peaks become smaller but continuously fluctuate. The sharpness of the peaks and their continuous distribution arise from the highly ordered structure of the solid. The peaks are sharp unlike the gas and liquid because all of the atoms are in fixed positions and so at a certain position there will be a number of atoms in an equivalent position around the sphere of the distribution function. Furthermore the solid's first peak is considerably larger than the other phases, the probability is much greater due to the closeness of the atoms. The integrated function shows a very large number of atoms in the space, this again illustrates the closeness of the atoms in the solid state.

2. Using the solid RDF to extract information on the lattice structure

figure 11

The three first peaks of the solid's radial distribution function are distances 1.0, 1.5 and 1.8. The lattice spacing is 1.5 and the calculated distances ra,rb and rc(where r is the distance from the reference atom to the atoms a, b and c) were found to be 1.0, 1.5 and 1.8(respectively). Therefore the first three peaks 1, 2 and 3 correspond to the lattice points a, b and c(respectively).

The coordination of atom a, b and c is 12, 6 and 24(respectively). These values were interpolated from the magnification of the integrated solid distribution function.

Dynamical properties and the diffusion coefficient

Task 19: The mean squared displacement(MSD) of large and small simulations

figure 12

MSD(t)=<((t0+t)(t0))>2 where; t0= the initial timestep and t=the new timestep.

The mean squared displacement(equation above) is a measure of how much random motion there is within the system. In this simulation the ability to move randomly is determined by the amount of free space and therefore the density. In the gaseous phase there is a lot of free space, the atoms move freely and encounter only a small number of atoms to collide with thus impeding it's motion; the linear increase of the plot reflects this behavior. The liquid has a lower degree of random motion than the gas, due closer proximity of the atoms, forces such as drag will now effect the atom's motion as it moves past other atoms. The two first plots show that the atom is able to move in drift, however the liquid experiences more drag and has a smaller gradient. The solid displays the least amount of translational freedom, the close packed structure prevents any random translational motion; it does however show some random movement, arising from vibrations of the constrained atoms. Therefore the plots in figure 12 are a reasonable representation of the random movement of the three phases.

Finding the diffusion coefficient D

D=16δ<r2(t)>δt> where D= the diffusion coefficient in ms2,<r2(t)= MSD and t= timestep The gradient of the slopes give the diffusion coefficient, however this uses a dimensionless timestep denominator, therefore we must convert it to time units by dividing by the timestep(0.002).

The gradient was found by fitting a linear trendline through the function.

Solid

Conversion of the gradient to give time units gave D=4.17x106ms2 for the million atom simulation and D=5.83x107ms2 for the 3200 atom simulation.

Liquid

Conversion of the gradient to give time units gave D=8.33x102ms2 for the million atom simulation and D=8.33x102ms2 for the 3200 atom simulation.

Gas

A trendline was fitted to the function, since the simulation took some time to produce the linear behavior it was fitted to be linear with the latter part of the function. Conversion of the gradient to give time units gave D=3.09ms2 for the million atom simulation and D=2.43ms2 for the 3200 atom simulation.

Task 20: The velocity auto-correlation function(VACF)

Evaluating the 1D harmonic oscillator VACF using the expression for the position as a function of time

The normalized VACF is given by;

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

We must simplify this expression in order for it to be useful, to do this we must first find expressions for v(t) and v(t+τ):

The expression for position as a function of time for the HO is defined as; x(t)=Acos(ωt+ϕ)

This can be substituted into the expression of velocity as a function of time;

v(t)=dx(t)dt=dAcos(ωt+ϕ)dt

This is then differentiated using the chain rule;

dx(t)dt=dx(t)duxdudt

Where u=ωt+ϕ

v(t)=ωAsin(ωt+ϕ)

We also need our velocity as a function of time after τ, this is acheived by substituting in t=t+τ;

v(t+τ)=ωAsin(ω(t+τ)+ϕ) =ωAsin(ωt+ωτ+ϕ)

The two velocity expressions are substituted into the VACF;

C(τ)=ωAsin(ωt+ϕ)ωAsin(ωt+ωτ+ϕ)dtωAsin(ωt+ϕ)2dt

=sin(ωt+ϕ)sin(ωt+ωτ+ϕ)dtsin(ωt+ϕ)2dt

Using the trigonometric identity sin(α+β)=sin(α)cos(β)+cosαsinβ we can re-write our equation;

where ωt+ϕ=α and ωτ=β

C(τ)=sin(ωt+ϕ)sin(ωt+ϕ)cos(ωτ)+cos(ωt+ϕ)sin(ωτ)dtsin(ωt+ϕ)2dt

C(τ)=sin(ωt+ϕ)2cos(ωτ)+cos(ωt+ϕ)sin(ωτ)dtsin(ωt+ϕ)2dt

The fraction can be expanded;

C(τ)=sin(ωt+ϕ)2cos(ωτ)dtsin(ωt+ϕ)2dt+cos(ωt+ϕ)sin(ωτ)dtsin(ωt+ϕ)2dt

Then since cos(ωτ) does not contain a t term it can be moved outside of the integral;

=cos(ωτ)sin(ωt+ϕ)2dtsin(ωt+ϕ)2dt+cos(ωt+ϕ)sin(ωτ)dtsin(ωt+ϕ)2dt

The alike terms cancel to give;

=cos(ωτ)+cos(ωt+ϕ)sin(ωτ)dtsin(ωt+ϕ)2dt

Finally cos(ωt+ϕ)sin(ωτ)dtsin(ωt+ϕ)2dt=0, this is because sin and cossin are odd functions. When odd functions are integrated between symmetrical boundaries the product is zero.

Therefore;

C(τ)=cos(ωτ)


Analysis of the Lennard-Jones(LJ) solid, liquid and harmonic osillator(HO) VACFs against position

figure 13 VACF plots for LJ solid and liquid simulations and the HO

The minima in the VACFs represent a change in the direction of the velocity, this caused by collisions with other atoms.

The liquid and solid VACFs are similar but have some important differences. The solid VACF reaches its minima much quicker than the liquid; the atoms in the solid are much closer together and so have more surrounding atoms, it therefore takes less time for atoms to collide. The decorrelated solid function shows some oscillation about C(τ)=0 this arises from the jostling about with other atoms in the fixed lattice and so we see a sort of harmonic oscillator motion. The liquid does not display this behaviour because they are free to move and so will be able to "rapidly destroy any oscillatory motion."[1] by the "diffusive motion".

The LJ VACFs become decorrelated after colliding with other atoms, after collision the subsequent velocities no longer depend on the initial velocity. The Harmonic Oscillator’s VACF is fully correlated with the initial velocity. It behaves periodically, always returning through the same points after a certain amount of time. The differences seen between these plots are rooted in the essence of the representations. The HO models two sphere’s connected by a spring, the spring may expand and contract with perfect elasticity. In addition to this there are no external forces which could cause the velocities to decorelate. The LJ representions are effected by surrounding atoms, the velocities loose their correlation after a collision.


Task 19: Integrating using the trapezium rule to get D

The plots of C(τ)vstime of the three phases in a small and large system were plotted along with their running integral.

D=130dτ<v(0)xv(τ)>=Vn+Vn+12x(Tn+1Tn)+rn13

The first part of this equation is the relationship between D, the diffusion coefficient and the VACF. The second part is the equation used in excel to calculate the running integral. Where V is the velocity, T is the timestep and r is the last integral calculated. The last value of the running integral was divided by 3 to give D.


Solid(figures 14-15)

figure 14 VACF and running integral vs time for 3200 atom solid simulation
D=1.8x104ms2 for the small solid simulation and 4.55x105ms2 for the million atom simulation. 
figure 15 VACF and running integral vs time for 1000000 atom solid simulation


Liquid(figures 16-17)

[[File:liquid8000.png|200px|thumb|left|figure 16 VACF and running integral vs time for 8000 atom liquid simulation]D=9.8x102ms2 for the small liquid simulation and 9.0x102sm2 for the million atom simulation.

figure 17 VACF and running integral vs time for 1000000 atom liquid simulation


Gas(figures 18-19)

figure 18 VACF and running integral vs time for 8000 atom gas simulation
D=2.34ms2

for the small gas simulation and

3.27ms2

for the million atom simulation.

figure 19 VACF and running integral vs time for 1000000 atom gas simulation


The diffusion coefficients are as would be expected for the respective states. An atom in a tight lattice is not free to diffuse and is locked in position, this is reflected in its very small value. The liquid has a more important diffusion coefficient since the atoms are free to move, however they do experience more impeding forces such as drag while passing other moving atoms. Finally the gas is also a reasonable estimation, it is much greater than the solid and liquid since the atoms are a lot more disperse and so move very quickly and encountering relatively few restrictive forces.


Sources of error in the estimates

A larger number of atoms included in the simulation gives a more accurate result, therefore the million atom simulation and the smaller one do have a small difference in values which can be accounted for by this. Another important source of error comes from the running integral calculation, each new area will have the previous area added to it, this results in a build-up of errors. The two solid diffusion coefficients calculated from the running integral are an order larger than the coefficient calculated using the mean squared displacement method.






References: