Wiki.ch.ic.ac.uk/wiki/index.php?title=Mod:PetrWiki

Optional experiment: Simulation of a simple liquid.

Introduction to molecular dynamics simulation.

The x(t) = Acos(ωt+φ) equation for a classical harmonic oscillator was used to fill the ANALYTICAL column.

Completed ERROR column by subtracting velocity-Verlet inter-atom distance solution from the calculated analytical inter-atom separation.

Completed ENERGY column by adding together the kinetic energy and potential energy of the oscillator, using the equations Kinetic energy = 1/2m(v(t))² and Potential energy = 1/2k(x(t))²

The graphs obtained from filling these columns are displayed below:

A timestep of 0.2 causes the energy of the system to oscillate between 0.495 and 0.5, this is a change in the total energy of 1%, therefore timesteps of 0.2 and above ensure that the total energy doesn't change by more than 1% over the course of the simulation.

It is important to monitor the total energy of a physical system when monitoring its behaviour numerically because this allows you to observe whether the system is at equilibrium or not, and therefore know whether the system is ready to have measurements taken from it that are representable of what the system is trying to model. Knowledge of the total energy also allows it (or related parameters) to be changed appropriately in a simulation, in order to constrain certain parameters within the simulation. For example, as temperature is related to total energy, energy data can be used to find temperature changes, and therefore be able to constrain temperature as required. In general, simulations are useful only when they are at equilibrium and represent a real body, which would tend to exist at equilibrium, and therefore it's important to know when a simulation is accurate by looking at energy. This is unless you are interested in observing a reacting system (or another unstable system), however even in a reacting system it is useful to know when the system stops reacting, or its rate of reaction speeds up/slows down by looking at total energy changes.

The peak values of the error within the system increase as time passes due to the analytically calculated position values for the harmonic oscillator becoming progressively more shifted away from those predicted by the velocity-Verlet algorithm, this is because the equations yield slightly different results, and the effects of these differences build up and become more and more noticeable with time. The error in the system oscillates together with the predictions for the inter-atom separation in the system, as the error calculation is just a subtraction, and therefore follows the osccilations of the two functions being substracted from each other in the calculation.

The equation for the potential energy of the Lennard-Jones interaction is:

To overcome the issue of being unable to directly simulate realistic quantities of liquid we model bulk liquid as an infinite series of boxes, each full of simulatable quantities of atoms, extending in all directions. When an atom hits the boundaries of a box, it gets teleported to the opposite side of the box whilst retaining its velocity from prior to its reappearing.

An atom in position (0.5,0.5,0.5), which moves along the vector (0.7,0.6,0.2) in a single timestep and resides within a cubic simulation box that runs from (0,0,0) to (1,1,1), will end up at the point (0.2,0.1,0.7) after one timestep and after periodic boundary conditions are applied.

Equillibration.

Giving atoms random starting coordinates causes problems in simulations because of the Lennard-Jones potential having very high potential energy values at low separations between atoms, meaning that if atoms are generated in random positions and happen to appear very close to each other, they will shoot off with very high kinetic energies away from each other in order to lose their high potential energy. This wouldn’t happen in real space, as atoms wouldn’t be able to get so close together so as to get such high potential energies in the first place.

There is one lattice point per unit cell in a simple cubic lattice. A lattice spacing of (1.07722, 1.07722, 1.07722) in (x,y,z) means that each unit cell has a volume of 1.25, and therefore a number density of lattice points of 1 lattice point per 1.25 volume = 0.8. There are four lattice points per unit cell in a face-centred cubic lattice. If the lattice point number density of such a lattice is 1.2 then:

If a box containing 10x10x10 of face centred cubic unit cells was created, and filled with atoms, it would have 4000 lattice points, as there are four lattice points per unit cell (and 1000 unit cells would be created), meaning that 4000 atoms would be created.

The command “mass 1 1.0” specifies the mass of a certain type of atom; in this case it specifies that the mass of atom type 1 (first argument) is equal to 1.0 LAMMPS units (second argument). The command “pair_style lj/cut 3.0” is used to compute interactions between atom pairs, in this case the lj/cut 3.0 term specifies that a Lennard-Jones potential is used to do this, with a maximum r value for calculating energies of interactions between atoms of 3.0 reduced distance units. The command pair coeff ** 1.0 1.0 specifies force field coefficients between specific atom types, in this case this is done for all atom type pairs possible in the simulation (denoted by double asterisk), with the force field coefficients between atom types always being 1.0.

Because we are specifying both x_i(0) and v_i(0), we are going to use the Velocity-Verlet integration algorithm, as it involves both terms in its use, and allows us to find velocities in the measured system.

We define a timestep variable in order that we can use the name of this variable instead of the numerical timestep value in code that involves the timestep in the input script, so that if we later want to change the timestep value used for the simulation, we can just change the value of the timestep in the variable definition, and all of the code using the timestep will automatically use the new updated timestep. Essentially, this procedure saves us time, and makes it easier to modify the timestep in existing scripts.

The system does reach equilibrium, this is almost instantaneous in the time periods covered by the simulation and takes about 0.38 reduced time units.

The 0.01 timestep is the largest to give acceptable results. The 0.015 timestep is a particularly bad choice because the system loses its equilibrium and starts increasing in energy due to lack of accuracy in the simulation (in regards to reality) due to overly high timestep, this happens as the simulation has to make predictions too far into the future at each calculation that it does.

Running Simulations Under Specific Conditions.

From the Equilibrium section simulations, it is evident that the timestep value of 0.001 is the most favourable, as it gives enough time for the system to equilibrate, and being the smallest timestep, yields results most similar to physical reality. It can be seen that the pressures of 2.4 and 2.6 in reduced units are reasonable for simulations (they are just above and just below the equilibrium pressure on the time against pressure graph for a timestep of 0.001). Therefore the combinations of (p,T) of (2.4,1.6), (2.4,1.7), (2.4,1.8), (2.4,1.9), (2.4,2.0), (2.6,1.6), (2.6,1.7), (2.6,1.8), (2.6,1.9), (2.6,2.0) are used for simulation, with a timestep of 0.001 and with all temperatures used being above the critical temperature of T=1.5 (all in reduced units).

In the command, “fix aves all ave/time 100 1000 100000 v_dens v_temp v_press”, the 100 argument signifies that the values for temperature etc. will be sampled every 100 timesteps. The 1000 argument shows that 1000 measurements will contribute to the average, and the 100000 argument shows that an average value will therefore be obtained after 100000 timesteps. The following “run 100000” term shows that 100000 timesteps will be simulated.

My simulated density is lower than the ideal gas law density, this is because there are no intermolecular forces acting between atoms in the ideal gas law model, meaning that atoms can approach each other as close as they like (and even pass through each other as they can’t collide), whereas in the Lennard-Jones model atoms won’t approach each other as closely due to repulsions between them at shorter separations, resulting in a less dense bulk species. The discrepancy between the models increases with pressure, as atoms are made more likely to approach each other closely, and the effect of repulsive forces at shorter separations becomes more evident.

Calculating heat capacities using statistical physics.

File:Heat Capacity Input - 0.2,2.0.in Petr

The output file gives the heat capacity for the entire modelled box, to convert this into heat capacity per unit volume, the heat capacity results for each simulation are divided by the volume of the box for the specific simulation.

The relationship between temperature and heat capacity is as expected, and implies that the simulations occur after the Schottky anomaly for the modelled system. Increasing the temperature of a system increases the average energy that particles in said system have, heat capacity reflects how easy it is to excite a system to a higher energy state; normally an increase in temperature means a higher proportion of particles exist at a higher energy and require less extra energy input to reach a higher energy state, and therefore heat capacity increases. However, after a certain point almost all particles have enough energy to freely move between energy levels, and therefore heat capacity begins to fall as it becomes harder to actually excite the system to a higher state, due to the separate energy states starting to act like a single energy state; this effect starts to happen after the Schottky anomaly and can be observed on the graph. The humps at temperature 2.4, are likely due to error. Increasing density causes there to be more atoms per unit volume in the simulated system, and therefore more things that can be excited to a higher energy state per unit volume. Heat capacity is calculated as being per unit volume, and therefore it increases together with density as shown on the graph.

Structural properties and the radial distribution function.

(The graphs should say "integrated" not "intergrated", I have noticed this too late).

The solid RDF shows a long range order with multiple peaks, implying that atoms are most likely to be found in fixed positions in regards to each other, this order continues up until the end of the graph. There is a uniform crystal lattice in this case. The Liquid RDF in comparison shows a shorter range order, with specific separations being more favoured, however not as perfectly as in the solid, and the distribution becoming uniform at higher separations, shown by the distribution becoming flat towards the end of the graph. The gas RDF has only one peak, which corresponds to the minimum on the Lennard-Jones curve for two atoms, the equilibrium distance on such a curve is given by r_eq=2^1/6σ. Because we are working in reduced units, this corresponds to r_eq=^21/6 for our simulation. 2^1/6=1.12 (3.s.f), which is very similar to the gas RDF peak separation of 1.175. The rest of the gas RDF is uniform, as the distribution levels off, this is due to the random free motion of atoms in a gas.

The integrated gas RDF goes up to ~213 (3.s.f), the integrated liquid RDF goes up to ~2850, and the solid goes up to ~4690, this shows that the solid has the most atoms closer to each other, the liquid has less, and the gas has the least. This matches theory as atoms in a solid are in a fixed lattice close to each other, and don’t move away but only vibrate. The atoms in a liquid can move away from each other, however tend to remain in close proximity to each other due to their movement away being significantly hindered by inter-atom attractions. Finally the atoms in a gas move almost freely relative to each other, and therefore are likely to be found far away from each other.

The first three peaks in the solid RDF have separations of 1.075, 1.525, and 1.875. These correspond to the smallest three inter-atom distances in the face centred cubic lattice which the solid is simulated as. These distances are depicted in the diagram below, with the smallest distance labelled a, the next b, and the longest c.

Therefore by trigonometry, the distance c corresponds to the RDF peak values (with slight deviation). The lattice spacing is equal to b, and is therefore 1.525.

To get the coordination number for each of the peaks, one needs to find the area underneath each peak, this can be done by looking at the integrated solid RDF distribution (easier to do on zoomed-in version), and finding the integral from 0 to the inter-atom separation value at the end of the peak in question. In order of increasing separation, the three peaks end at the points: (1.325, 12.07576227), (1.675, 18.0311431), and (2.025, 42.00453648) on the integrated solid RDF, the first peak has an area of ~12 under it, implying a coordination number of 12. The area under the second peak is 18.0311431 - 12.07576227 = 5.955 (4.s.f) implying a coordination of ~6. The third peak has an area underneath it of 42.00453648 - 18.0311431 = 23.97 (4.s.f) implying a coordination number of ~24. These coordination numbers match those that would be expected from a face centred cubic lattice, as mentally imagining one lattice point which is surrounded by eight unit cells (the lattice point is on the corner of each unit cell), yields 12 separations of distance a to other lattice points from the central lattice point, 6 separations by a distance of b, and 24 separations by a distance of c (3 such separations per surrounding unit cell).

Dynamical properties and the diffusion coefficient.

Generally, the mean squared displacement (MSD) plots agree with theory. The MSD of an atom is a quantification of the amount of the system that said atom has “explored”, this means the amount of space in the system that the atom has visited at some time during its motion during the period of time in question. Atoms in a gas move around freely, and therefore an increase in the MSD would be expected with increasing time, as the atoms explore more and more space. This is visible in the graphs below:

The rate of MSD increase increases with time on the graphs above, this may be explained by the level of disorder in the system progressively increasing, starting from the uniform lattice in which the atoms initially find themselves in. As the system becomes progressively more like true random free motion, and atoms being less likely to have their motion coupled to the motion of other atoms, the atoms start exploring the system more rapidly and the rate of MSD increase increases. This rate increase will level off over time as maximum disorder is reached, in the limited space of the modelled box.

In the liquid MSD graph, the MSD increases linearly with time. Motion in a liquid generally involves atoms moving around close to each other with a strong effect of the inter-atom forces on the distance covered by the motion of the coupled atoms. As the atoms constantly move around each other to cover more space, the rate of exploration stays constant for each atom, and therefore the MSD increases linearly with time.

In a solid, the atoms are in a fixed lattice, and therefore can’t explore regions outside of those that they cover by their vibration, therefore the solid initially has a rapid increase in its MSD, as the atoms explore the regions of space just around them (by vibrating), but then the MSD stops changing (on average and accounting for error) as the atoms can’t explore any more regions of space.

The MSD for a gas has the highest maximum value, followed by the liquid, and then the solid, with each being different to the prior by a factor of about 100. This is due to the motion of atoms being more free from each other with increasing disorder in the physical state that the atoms find themselves in, causing them to be less hindered and able to explore the system more rapidly. Also atoms in more disordered states are under higher temperatures, and therefore have more energy to explore the system more quickly (the same effect happens due to lower pressures).

With <r² (t)> =MSD, it is possible to estimate the values of D for each simulated system by measuring the gradient of the MSD against time.

For the gas, the gradient can either be taken as 15.216, from a linear line of best fit, giving a D value of 2.533 (3.s.f, no units because calculated from unit-less properties), or, by measuring the gradient from a line of best fit only on the second half of the graph where the trend becomes more linear, to give a gradient of 18.705, and a D value of 3.118 (4.s.f). This approach can be repeated with the million atom simulation output data, giving D values of 2.542, or 3.083. The million atom simulations tend to give similar results to the simulations with 8000 atoms (32000 atoms for solid as it starts from an fcc lattice arrangement), implying that using the smaller numbers of atoms in simulations is likely sufficient to get accurate results.

The gradient for the liquid simulation is 0.4804, giving a D value of 0.08007 (4.s.f), the million atoms simulation gives a D value of 0.08733. The solid simulation gives a gradient of ~0 for both simulations and D values of ~0. More exact D values can be calculated by taking a line of best fit on values past “Time = 2” for the solid (after the initial increase in MSD due to atom vibration), and are 1.5 x 10^-6 for the 34000 atom simulation and -8.33 x 10^-6 million atom simulation. The small non-zero nature of these D values is likely due to the solid lattice not being perfectly fixed, and some atoms changing position. Alternatively error could have caused this result, which is supported by the strange negative gradient value (and therefore D value) for the million atoms solid MSD, which implies that the atoms explored less and less overall space with passing time.

The minima in the VACF’s for the liquid and solid represent the atoms' velocities becoming negative as their movement direction oscillates due to the effect of inter-atom forces; this causes the scalar product of their initial and time τ velocities to be negative, and therefore C(τ) to be negative. This effect is less pronounced for the liquid, as a lot of the atomic motion is translational rather than vibrational due to the lack of a fixed lattice, resulting in fewer, less exaggerated minima. Overall the differences between the solid and liquid VACFs are due to the solid having a regular lattice of atoms, which can’t leave their lattice points, apart from by vibration around the centre of their lattice point, and the liquid having atoms which can move around one another and not just vibrate. This causes the liquid to equilibrate at C(τ) = 0 more rapidly as a lower proportion of its atoms’ motion is vibrational motion that can make the liquid atoms change direction and cause more oscillation in the VACF.

The Lennard-Jones potential is non-symmetric, with the left side of the curve having a high potential energy asymptote, and the right flattening out in energy. This effect causes atoms existing at the same potential energy to be more likely to be distributed on the right side of the curve at this potential energy than the left, after the system is allowed to equilibrate from the starting conditions, as atoms will tend to stay at separations away from the potential energy asymptote. This effect causes a constant change to be applied on the system, which will change the velocity of the involved atoms, as it changes their distance parameters. The velocity autocorrelation function represents how different the initial velocity of a species is to its velocity at time τ, if there is a constant pressure for change being applied to the velocity of a species, the autocorrelation function will decrease with time, this causes the damping effect observed in the VACF for the simulated solid and liquid. The harmonic oscillator on the other hand has a symmetrical potential energy curve, meaning that there is no directionality in the model, and no pressure for uncorrelation of its velocity at time τ from its initial velocity value.

The running integrals are as expected, and appear to represent the areas underneath their respective VACF curves.

By trapezium rule, the areas under the VACF curves for the gas, liquid and solid are: 9.883762, 0.261569, and 0.000367 respectably.

Both the MSD and VACF methods give similar results for each D value, implying that both work correctly. VACF D values are as expected, with more disordered physical states having greater particle motion, represented by higher values of D, and caused by the greater energy of particles allowing them to over come inter-atom forces more and travel around further, faster, and more freely (diffuse more rapidly). The solid D values are less similar between methods, however it is hard to predict values of such small magnitudes accurately due to error in the system.

The largest source of error in my estimates of D from the VACFs is due to the inaccuracy of the trapezium rule method, the VACF curves are not linear, and therefore can't be integrated by the trapezium rule with good accuracy unless much thinner trapezia are used.