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.

Output of your First Simulations

In the first section, you submitted five simulations to the HPC queue. These should now be finished. Return to the HPC portal webpage and click Job List in the menu. A list of all the jobs that you have submitted will be displayed (figure XXX). The five jobs from the introduction section (which you named Intro0.001, Intro0.0025, ..., Intro0.015) should all display Finished in the Status column. Use the dropdown list and Download button in the Output files column to download the Log file for each job. By default, this will download a filename named logfile_XXXXXX.txt, where XXXXXX is a number unique to your job. You should also download the Trajectory file for each job - this contains the the positions of all the atoms in your simulation as a function of time, allowing you to visualise the behaviour of your system.

You should move these files somewhere safe, as you will need them later, and give them the same name as the job that they came from! For example, the log file from job Intro0.001 should be named Intro0.001.log. Give the trajectory files the same names, but with the file extension .lammpstrj — the trajectory file from job Intro0.001 should become Intro0.001.lammpstrj. This will make it easier to visualise the trajectory files later on.

NEED TO SPEAK TO MATT ABOUT GETTING THE ACTUAL LOG.LAMMPS FILE AS AN OUTPUT

Open the file "Intro0.001.log" in Notepad. This is a text file, just like the input files that we prepared. It contains the information generated by LAMMPS as it executed your simulation. Compare it to the input file, "0.001.in". LAMMPS reads the input file, "0.001.in", one line at a time. It prints that line out to the log file, and then runs whichever command that line tells it to. As it runs this command, it usually prints extra information to the log file. For example, where "0.001.in" contains the line

lattice sc 0.8

the log file contains the lines

lattice sc 0.8
Lattice spacing in x,y,z = 1.07722 1.07722 1.07722

In the section below, we will look at a few key points from the output file. If a line appears in the log file and is not mentioned below, it is because the purpose of that line is rather technical and not relevant to the key chemistry and physics of this experiment. The demonstrators will be happy to discuss the meaning of these lines with you if you are interested in them. Any line in the log file or input file which begins with a hash, #, is called a comment. These lines are only to make notes in the file about what different commands do - LAMMPS will read them, but ignore them completely.

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.

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.

TASK: Why do you think giving atoms random starting coordinates causes problems in simulations? Think about what happens to the Lennard-Jones interaction between two atoms as they get very close together.

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 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 the 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!).

TASK: Satisfy yourself that this lattice spacing corresponds to a number density of lattice points of $0.8$. 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?

The next lines in the log file are

region box block 0 10 0 10 0 10
create_box 1 box
Created orthogonal box = (0 0 0) to (10.7722 10.7722 10.7722)

The region command 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 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 1000 (10x10x10) unit cells of this lattice, and so contains 1000 lattice points. We now need to fill our simulation box with atoms.

create_atoms 1 box
Created 1000 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.

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?

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.

TASK: Using the LAMMPS manual, find the purpose of the following commands:

mass 1 1.0
pair_style lj/cut 3.0
pair_coeff * * 1.0 1.0

So far we have created 1000 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.

TASK: Given that we are specifying $\mathbf{x}_i\left(0\right)$ and $\mathbf{v}_i\left(0\right)$, which integration algorithm are we going to use?

Choosing initial velocities for the atoms is a little easier than choosing initial positions. From the statistical thermodynamics 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

TASK: Look at the lines below.

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

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.

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.

• loading a trajectory
• track a single particle or pair of particles to see how much it moves

Checking equilibration

When we first set up a simulation, it is very important to make sure that our system reaches an equilibrium state.

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.