Third year simulation experiment/Equilibration

This is the third section of the third year simulation experiment. You can return to the previous page, Introduction to molecular dynamics simulation, or jump ahead to the next section, Running simulations under specific conditions.
We will be using the LAMMPS program to carry out our molecular dynamics simulations. In several places in this section, we will ask you to consult the LAMMPS manual to find out things about how the software works. You can find the manual here. We appreciate that the format of this document can make it a little hard to navigate, but it is the definitive resource on how different commands in LAMMPS work, and is therefore invaluable. The files you will need for this section can be found in the intro folder downloaded previously.

Creating the simulation box

In the previous section, it was pointed out that before we can start a simulation, we need to know the initial states of all of the atoms in the system. Exactly what information we need about each atom depends on which method of numerical integration we need, but at the very least we need to specify the starting position of each atom. If we wanted to simulate a crystal, this information would be quite easy to come by — we could just look up the crystal structure, and use that to generate coordinates for however many unit cells we wanted. For this purpose, LAMMPS includes a command which generates crystal lattice structures.

Generating coordinates for atoms in a liquid is more difficult. There is no long range order, so we can't use a single point of reference to work out the positions of every other atom like we can in a solid. We could generate a random position for each atom. This would certainly create a disordered structure, but causes larger problems when we try to run the simulation.

Instead, we are going to place the atoms on the lattice points of a simple cubic lattice. This, of course, is not a situation in which the system is likely to be found physically. It turns out, though, that if we simulate for enough time we will find that the atoms rearrange themselves into more realistic configurations. We will discuss towards the end of this section exactly what is meant by "enough time"!

Consider the line in the input file

lattice sc 0.8

This command (further info) creates a grid of points forming a simple cubic lattice (one lattice point per unit cell). The parameter

$0.8$

specifies the number density (number of lattice points per unit volume). In a corresponding output file, you will see the line

Lattice spacing in x,y,z = 1.07722 1.07722 1.07722

This indicates that the distance between the points of this lattice is

$1.07722$

(in reduced units, remember!).

The next lines in the input file are

region box block 0 5 0 5 0 5
create_box 1 box

The corresponding log file output is

Created orthogonal box = (0.0000000 0.0000000 0.0000000) to (5.3860867 5.3860867 5.3860867)

The region command (further info) simply defines a geometrical region in space, which we call "box". In this case, "box" is a cube extending ten lattice spacings from the origin in all three dimensions. The subsequent create_box command (further info) tells LAMMPS to use the geometrical region called "box" as a template for the simulation box. The number 1 between "create_box" and "box" indicates that our simulation will contain only one type (species) of atom.

So far we have defined a simulation box which is based around a virtual simple cubic lattice. Our box contains 125 (5x5x5) unit cells of this lattice, and so contains 125 lattice points. We now need to fill our simulation box with atoms. The input command is

create_atoms 1 box

(further info) while the log file simply contains an acknowledgement of this

Created 125 atoms

The create_atoms command has two arguments; the first tells LAMMPS that all of the atoms that we create will be of type 1. Every atom in the simulation has a type — because we will be simulating a pure fluid, containing only one chemical species, every atom will have the same type. The actual type that we assign to each atom is arbitrary — type 1 does not, for example, need to correspond to the element with atomic number 1 (hydrogen). If we wanted to simulate water, we might make the hydrogen atoms type 1 and the oxygen atoms type 2. We will specify the physical and chemical properties of each atom type later in the input script.

The remaining data in the log file isn't very instructive as it stands — it simply contains a list of the thermodynamic properties of the simulation at certain intervals. In a few sections time, we will plot this data, but for now you can close the log file. Keep the input script open.

Setting the properties of the atoms

In addition to their positions, we also need the physical properties of the atoms to be able to perform the simulation. We set these properties on a 'per-type' basis, so that every atom of the same type has the same mass and the same interactions.

So far we have created 125 atoms, and we know the starting ($t=0$) position for each of them. We have also set their masses, and told LAMMPS what sort of forces to calculate between them. The final thing we need to specify to completely specify the initial conditions is the velocity of each atom.

Choosing initial velocities for the atoms is a little easier than choosing initial positions. From the 1st year lectures, you should know that, at equilibrium, the velocities of atoms in any system must be distributed according to the Maxwell-Boltzmann (MB) distribution. If we know the masses of the atoms, and we know what temperature we want to simulate, then we can determine the relevant MB distribution function. LAMMPS is able to give every atom a random velocity whilst ensuring that overall the MB distribution is followed. This is the purpose of the line

velocity all create 1.5 12345 dist gaussian rot yes mom yes

You can see the manual page for this command here, but the key sections are:

• all: the group of atoms on which the command acts. all simply specifies that we want every atom to have a velocity assigned to it.
• 1.5: the temperature, $T$, needed to calculate the MB distribution(in reduced units, as always)

Monitoring thermodynamic properties

We need to be sure that our simulation is correctly modelling whatever physical system we want to study. It is relatively easy to set up simulations, but how can we be sure that the "results" we get make sense? One of the best ways is to calculate from the simulation things that we can measure in experiment, and see if they agree. For example, we might want to simulate our system at a particular temperature and pressure, and measure the resulting density. If we repeat this over a range of temperatures at the same pressure, we will be able to plot an equation of state, which we could compare to experimental measurements.

LAMMPS is able to calculate a great deal of thermodynamic information for us (you can see a full list of the properties it is able to calculate here), but in these first simulations we are only interested in those properties specified in these commands:

thermo_style custom time etotal temp press
thermo 10

The first controls which properties will be printed out in the log file. In this case, we print how much time we have simulated so far (which is not the same as how long it has taken us to simulate it!), the total energy of the atoms, their temperature, and their pressure. The second line tells LAMMPS to print this information on every 10th timestep.

Running the simulation

Look at the lines below.

### SPECIFY TIMESTEP ###
variable timestep equal 0.001
variable n_steps equal floor(50/${timestep})
timestep ${timestep}

### RUN SIMULATION ###
run ${n_steps}

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.

It is now time to run your first simulation, submit the input script with the data file in the intro folder of the files you have downloaded Try changing the timestep - what happens when you make the timestep larger?.

Visualising the trajectory

The trajectory files contain the positions of all the atoms in the simulation, recorded at a set interval (for all of these simulations, this was every ten timesteps — this is controlled by the dump command in the input scripts). We use a programme called VMD to view these trajectories, which you should find is already installed on both the desktop and laptop computers. You can run VMD from the start menu with Start -> All Programs -> University of Illinois -> VMD.

Loading a Trajectory

We'll start by looking at the output of the 0.02 timestep simulation. In the VMD Main window, select the menu option File -> New Molecule. Click the Browse button, then select the relevant trajectory file. In the Determine file type dropdown, select LAMMPS Trajectory. Then click Load.

You will see that the VMD 1.9.1 OpenGL Display window now shows a horrible mess. VMD's default behaviour is to draw lines between atoms which it thinks might be chemically bonded. Our system doesn't model chemical bonds, so we want to turn this off. In the VMD Main window, select the menu option Graphics -> Representations. This shows a list of "representations" of our atoms. You will see that at the moment, there is a single representation listed, and it is selected. It will have the Lines style, the Name colour, and the selection all. "Selection" simply tells VMD which atoms we want it to draw. We want to show every atom, so the current selection is fine. The name colouring method just makes VMD give atoms colours according to their specified type. The colour isn't important to us, so we can leave this be too. The "style" tells VMD what we want it to display for each atom. Change the Drawing Method from Lines to VDW. You will see that the mess of lines is replaced by a mess of low resolution, overlapping spheres. Change the Sphere Scale to 0.3, and the Sphere Resolution to 17. The result should look a little smoother. Close the Graphical Representations window. You will notice that in the bottom right of the VMD Main window, there is a small play button. Click this, and you will see the animated version of your simulation trajectory.

By clicking and dragging with the mouse, you can rotate the simulation box (though this may be sluggish). At any time, you can reset the view by pressing the equals key.

Tracking a Single Particle

To illustrate the periodic boundary conditions that we are using, we are going to draw almost all of the atoms as points, but we will pick a single atom at random to draw as a sphere. This will make it easy to see how a single atom moves through the box. Reset the display using the equals key, then use the VMD Main window controls to pause the trajectory and reset it to the first trajectory (play with the different buttons until you find the one that does this). You should see the perfect cubic lattice. Use the option Display -> Orthographic to change the drawing mode, then rotate the displayed crystal so that you are looking at one vertex (looking down the 111 direction, in crystallographic terms).

Open the Graphical Representations window again. Change the representation style from VDW to points, then click the Create Rep button. This creates a second representation, allowing a subset of the atoms to be drawn in a different way. The Selected Atoms box allows us to choose which atoms this representation applies to. We just want to pick two of them at random — VMD assigns every atom an index, from 0 to N-1. In our case, there are 125 atoms, so choose two numbers between 0 and 124. Changed the Selected Atoms field to

index i or index j

where i and j are your chosen numbers, press return, then change the Drawing Method to VDW. You should now see only two atoms represented by spheres, with the rest shown as small points. In the VMD Main window, click play. Try rotating the box, and changing the playback speed.

You will see that sometimes one of the spheres seems to change position across the box very rapidly — this occurs when it reaches one periodic boundary, and is reflected back across the other face. Try playing with some of the other representation types in VMD — it is a very powerful package, which is often used to render images of simulated proteins, so many of its options aren't relevant to our simple system!

Checking equilibration

When we first set up a simulation, it is very important to make sure that our system reaches an equilibrium state. We characterise equilibrium by the average values of thermodynamic quantities becoming constant (due to the approximations that we have made, there will always be fluctuations, but the average values will become constant).

In this section, we are going to plot the thermodynamic output of the simulation to see how long it takes to reach the equilibrium state (and indeed, whether this happens at all). Instructions are given below to import data from the LAMMPS log file into Microsoft Excel. Once you have the data in a spreadsheet, you can plot it. If you know how to use some of the other plotting software available on the chemistry computers (like Origin), you are welcome to use it.

1. Open a blank Excel workbook
2. Copy the data in the textfile into the first cell
3. With these data highlighted, click the Data tab and "Text to Columns"
4. Click "Delimited", continue and let it be space delimited
5. Click finish

You can then export this and data as a .csv and analyse in Python or Excel, as you wish.

Challenge: Can you write a python script or function that extracts this data for a given file automatically?

TASK 7:

Of the timesteps that you used, which timestep will you use for subsequent simulations and why? [6]

(Think about what is happening "physically" as you increase/decrease the timestep. Also, what features of each timeseries are indicative of the simulation's "health"?)

This is the third section of the third year simulation experiment. You can return to the previous page, Introduction to molecular dynamics simulation, or jump ahead to the next section, Running simulations under specific conditions.