Mod:MH2015LS

Molecular Dynamics - Liquid Simulation

A third year computational chemistry exercise by Michael Howlett.

Abstract

wordswordswordswords wordswordswordswordswordswordswordswordswordswordswordswordswordswordswordswords wordswordswordswordswordswordswordswordswordswordswordswordswordswordswordswords wordswordswordswordswordswordswordswords wordswordswordswordswordswordswordswordswordswordswordswords wordswordswordswordswordswordswordswords wordswordswordswords

Classical Harmonic Oscillator

TASK: Complete the three columns "ANALYTICAL", "ERROR" and "ENERGY" in the excel file provided.

A comparison of the classical harmonic oscillator with the velocity-Verlet method.

'ANALYTICAL' values for the position of a harmonic oscillator given by: $x(t)=Acos(\omega t+\phi)$

'ERROR' values are given by the modulus of the difference between the harmonic oscillator and velocity-Verlet results for position:

'ENERGY' values represent the total energy of the oscillator as calculated in the velocity-Verlet system by the sum of the potential and kinetic components:

Graphical comparison of these results to the Velocity-Verlet solution gives:

[[File:Energy graph.jpg|thumb|Energy vs time elapsed]]

[[File:Image name.jpg|thumb|right|Caption for the image]]

TASK: For the default timestep value, 0.1, estimate the positions of the maxima in the ERROR column as a function of time. 

Plot of the error maxima as a function of time:

time	2.00	4.90	8.00	11.10	14.20
Error value /x10-3	0.758	2.00	3.30	4.60	5.91

TASK: Experiment with different values of the timestep.

The limiting timestep to prevent >1% change in the total energy during the simulation is = 0.200 as shown by the inserted graph.

Why do you think it is important to monitor the total energy of a physical system when modelling its behaviour numerically?

Conservation of energy: For an isolated simulated system, the total internal energy should remain constant to reflect that seen in experiment? Any fluctuations seen are likely due to the approximations made in the computation and should be minimised in order to increase the accuracy of the simulation. ?

Lennard-Jones Potential

TASK: For a single Lennard-Jones interaction, , find the separation, , at which the potential energy is zero. What is the force at this separation? Find the equilibrium separation, , and work out the well depth ( ). Evaluate the integrals , , and when .

TASK: The Lennard-Jones parameters for argon are . If the LJ cutoff is , what is it in real units? What is the well depth in  ? What is the reduced temperature in real units? r = r*σ = 3.2 * 0.34 nm = 1.088 nm = 1.088 x 10^-9 m From the derivative of V(r) we know that for rmin, the well depth is = ε Therefore for the argon well depth, find ε. From the argon L-J parameters ε = kb * 120 = 8.314E-3 *120 = 9.98E-1 kJ/mol T = Tε/Kb = 1.5*120 = 180 K We denote these reduced quantities by a star, and they take the following conversion factors: distance energy temperature

Simulation theory - additional exercises and questions

TASK: Estimate the number of water molecules in 1mL of water under standard conditions. Estimate the volume of water molecules under standard conditions. Density of pure water at 0°C (STP) = 0.999 g/mL, at 25°C (SATP) = 0.997 g/mL STP: Mr = 18.015 g/mol Therefore for 0.999/ 18.015 = 0.05545 moles 0.05545 x 6.022 x 10^23 = 3.3394 x 10^22 molecules of H2O. SATP: Mr = 18.015 g/mol Therefore for 0.997/ 18.015 = 0.05534 moles 0.05534 x 6.022 x 10^23 = 3.3327 x 10^22 molecules of H2O.

10000/6.022 x 10^23 = 1.6606 x 10^-21 moles x 18.015 g/mol = 2.9915 x 10^-20 g / 0.999 = 2.9945 x10^-20 mL / 0.997 = 3.0005 x10^-20 mL Absolutely tiny fraction of 1 mL, i.e. an impossibly small volume for real world use. TASK: Consider an atom at position in a cubic simulation box which runs from to . In a single timestep, it moves along the vector . At what point does it end up, after the periodic boundary conditions have been applied?. Without periodic boundary conditions, a real atom would end up at position vector [1.2, 1.1, 0.7] with respect to the origin in the starting box. However with the periodic boundary conditions applied in a simulation, the particle will 'reappear' at the opposite side of the box as it crosses the boundary into the adjacent simulated box. Therefore, traveling [0.7, 0.6, 0.2] from the middle of the original box, [0.5,0.5,0.5], the atom will end the timestep in position [0.2, 0.1, 0.7].

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? Were atoms to be too close to each other in their initial state then there is the possibility that this could lead to strong intermolecular interactions between these (due to the large dependency of the L-J potential on distance r) and atom pairs could be formed at the equilibrium bond length. Taking this idea further, clusters of atoms could become strongly attracted to oneanother and this setup wouldn’t accurately simulate liquid behavior.

TASK: Satisfy yourself that this lattice spacing corresponds to a number density of lattice points of . 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? lattice point density= (number of lattice points per lattice ×total number of lattices)/〖(number of lattices on side of cube ×lattice side length)〗^3 lattice point density= (1 ×〖10〗^3)/〖(10 ×1.07722)〗^3 =0.8 10 x 10 x 10 cubic lattices of side length 1.07722 each. Therefore the simulation comprises of a volume of 10.77223 = 1250.009 volume units. For the 1.2 fcc density case: lattice side length = ∛(((lattice point density)/(number of lattice points per lattice ×total number of lattices))^(-1) )/(number of lattices on a side of cube)

lattice side length = ∛((1.2/(4 ×〖10〗^3 ))^(-1) )/10=1.49380 , confirmed by simulation of fcc. As shown in diagram, with an atom at each lattice point, the inter-atom spacing is half that of the lattice unit cell side length. Lattice spacing = 1.07722 Lattice spacing = 1.49380 / 2 = 0.7469

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 10 x 10 x 10 cube of fcc lattices. Atom at every lattice point and 4 lattice points per fcc lattice. Therefore for the pre-set simulation box of 10 x 10 x 10 for a 0.8 lattice point density with a fcc system the number of atoms produced = 4000. 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 Specify the pairwise force field coefficients for one or more pairs of atom types.

Simple cubic (left) and face-centred cubic (right) lattice types. [image]

TASK: Look at the lines below.

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

1. 1. 1. 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 Ask the demonstrator if you need help.

Because it enables you to change a single line that specifies the value of the variable ‘timestep’ and the change will be propagated through the rest of the script without extra work.

Equilibrium simulations

TASK: Given that we are specifying and , which integration algorithm are we going to use? Velocity-Verlet Method

TASK: 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? The system does reach an equilibrium in which the temperature and pressure fluctuate about an average value. The total energy is, as expected, constant throughout. As shown by the plot the system reaches this average equilibrium state within 0.2 - 0.4 seconds. How long does this take? When you have done this, make a single plot which shows the energy versus time for all of the timesteps (again, attach a picture to your report). Choosing a timestep is a balancing act: the shorter the timestep, the more accurately the results of your simulation will reflect the physical reality; short timesteps, however, mean that the same number of simulation steps cover a shorter amount of actual time, and this is very unhelpful if the process you want to study requires observation over a long time. Of the five timesteps that you used, which is the largest to give acceptable results? Which one of the five is a particularly bad choice? Why? 0.0025 is a good one to choose. 2.5x the actual time covered vs the 0.001 but with a very small change in the average energy value according to the graphs.

Running simulations under specific conditions

TASK: Choose 5 temperatures (above the critical temperature ), 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 points. You should be able to use the results of the previous section to choose a timestep. Submit these ten jobs to the HPC portal. While you wait for them to finish, you should read the next section.

TASK: We need to choose so that the temperature is correct if we multiply every velocity . We can write two equations:

Solve these to determine . γ^2 1/(2 ) ∑_i▒〖m_i v_i^2= 3/2〗 Nk_B T Given

γ^2 3/2 Nk_B T= 3/2 Nk_B T Therefore γ= √(T/T) 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? Nevery = use input values every this many timesteps Nrepeat = # of times to use input values for calculating averages Nfreq = calculate averages every this many timesteps 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 . 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?

Calculating heat capacities using statistical physics

TASK: As in the last section, you need to run simulations at ten phase points. In this section, we will be in density-temperature phase space, rather than pressure-temperature phase space. The two densities required at and , and the temperature range is . Plot as a function of temperature, where 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.

Structural properties and the radial distribution function

TASK: perform simulations of the Lennard-Jones system in the three phases. When each is complete, download the trajectory and calculate and . 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?

Dynamical properties and the diffusion coefficient

TASK: In the D subfolder, there is a file liq.in that will run a simulation at specified density and temperature to calculate the mean squared displacement and velocity autocorrelation function of your system. Run one of these simulations for a vapour, liquid, and solid. You have also been given some simulated data from much larger systems (approximately one million atoms). You will need these files later. 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 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.

〈r^2 (t)〉=2dDt+C 1/2d ∂〈r^2 (t)〉/∂t=D D = diffusion coefficient Units: distance2/unit time, d = number of dimensions, C = constant (representing in the solid the mean squared displacement of a simulated fcc lattice at equilibrium) D = gradient /2d Own simulations: Ds = 3.82x10-8/ (2 X 3) = 6.367 x10-9 Dl = 5.09x10-1/ (2 X 3) = 8.483 x10-2 Dv = 9.41/ (2 X 3) = 1.568 Million atoms: Ds = 3.18x10-8/ (2 X 3) = 5.300 x10-9 Dl = 5.24x10-1/ (2 X 3) = 8.733 x10-2 Dv = 1.88x101/ (2 X 3) = 3.133 Unitless but not dimensionless (m^2/s) i.e length^2/time

"uncorrelated” http://www.ccp5.ac.uk/DL_POLY/Democritus/Theory/msd.html Write in own words “At very short times however, the plot is not linear. This is because the the path a molecule takes will be an approximate straight line until it collides with its neighbour. Only when it starts the collision process will its path start to resemble a random walk. Until it makes that first collision, we may say it moves with approximately constant velocity, which means the distance it travels is proportional to time, and its MSD is therefore proportional to the time squared. Thus at very short time, the MSD resembles a parabola. This is of course a simplification - the collision between molecules is not like the collision between two pebbles, it is not instantaneous in space or time, but is `spread out' a little in both. This means that the behaviour of the MSD at short time is sometimes more complicated than this MSD plot shows.” Democritus

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!):

Be sure to show your working in your writeup. On the same graph, with x range 0 to 500, plot with 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.

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

Bold text