View source for Third year simulation experiment/Introduction to molecular dynamics simulation

'''<big><span style="color:blue; ">This is the second section of the third year simulation experiment. You can return to the previous section, [[Third year simulation experiment/Files to download|Downloading Files]], or jump ahead to the next section, [[Third year simulation experiment/Equilibration|Equilibration]].</span></big>'''

'''<big>This section contains background information about the theory of molecular dynamics simulations. It contains a number of relatively short exercises that you must complete as part of your lab write-up. These are labelled in bold and preceded by the word TASK in large print. It is recommended that you read the information on this page before carrying on with the rest of the experiment, but you are encouraged to save the TASKS for later; you can attempt them while you wait for long simulations to finish.</big>'''

In this section, we briefly discuss the theory behind molecular dynamics (MD) simulations. When we perform MD, we calculate how a particular set of atoms move over time. Using statistical physics, we can use the positions, velocities, and forces, of the atoms to calculate thermodynamic quantities like temperature and pressure.

==Theory==

===The Classical Particle Approximation===

As you may remember from your quantum chemistry lectures, it is very straightforward to write down the Schroedinger equation that describes the behaviour of any particular chemical system. For anything more complicated than a hydrogen atom, however, it is impossible to solve exactly. Even approximate solutions can be extremely computationally demanding. To be able to simulate a real system, we have to make some approximations.

It turns out that, to a very good approximation, we can assume that atoms behave as classical particles. Imagine a collection of <math>N</math> atoms. Each one of them will interact with all of the others, and so each atom will feel a force. Newton's second law tell us that that force causes the atom to accelerate.

Throughout this section, we are going to use the following notation:

* <math>\mathbf{F}_i</math> is the force acting on atom i.
* <math>m_i</math> is the mass of atom i.
* <math>\mathbf{a}_i</math> is the acceleration of atom i, the rate of change of its velocity.
* <math>\mathbf{v}_i</math> is the velocity of atom i.
* <math>\mathbf{x}_i</math> is the position of atom i.

<center><math>\mathbf{F}_i = m_i \mathbf{a}_i = m_i \frac{\mathrm{d}\mathbf{v}_i}{\mathrm{d}t} = m_i \frac{\mathrm{d}^2 \mathbf{x}_i}{\mathrm{d}t^2} \ \ (1)</math></center>

This is a second order differential equation for the positions of the atoms &mdash; if we know how the force, <math>\mathbf{F}_i</math>, behaves as a function of time, then we can determine the atomic positions and velocities at any time we like. Our system of <math>N</math> atoms has <math>N</math> of these equations, one for each atom. This is one of the reasons that computer simulations are needed &mdash; if we want to model the behaviour of a liquid, we can hardly solve the necessary number of equations by hand.

===Numerical Integration===

'''Numerical integration is a rather complex topic. In particular, the notation used below can be quite intimidating. Remember that you are encouraged to ask for help from the demonstrator if you want to discuss any part of this experiment!'''
<br>

There are a number of numerical algorithms to perform a molecular dynamics simulation, two are presented here- The Classical Verlet algorithm and The Velocity-Verlet algorithm.

====Verlet Algorithm====

To solve these equations numerically we have to ''discretise'' the problem: rather than treating the atomic positions, velocities, and forces as continuous functions of time, we break our simulation up into a sequence of '''timesteps''', each of length <math>\delta t</math>. This process is illustrated for a simple function in '''figure 1'''. The method that we are going to use to solve Newton's law for our atoms is usually called the Verlet algorithm (although it is an old method, and has been 'rediscovered' many times!). To understand its origin, we will begin with a brief derivation.

<br>
[[File:ThirdYearSimulationExpt-Intro-Discretisation.png|300px|thumb|center|'''Figure 1''': Discretisation of sin(x) between 0 and <math>2\pi</math>]]
<br>

We denote the position of an atom, <math>i</math>, at time <math>t</math> by <math>\mathbf{x}_i \left(t\right)</math>. Similarly, <math>\mathbf{v}_i \left(t\right)</math> is the velocity of that atom at the same time. What we want to know is the position of the atoms at the next timestep, <math>t + \delta t</math>. The basic Verlet algorithm is shown in '''figure 2''' - knowing a set of initial conditions the algorithm calculates forces and by Newton's second law, we can update the positions of a set of particles are a time <math>t + \delta t</math>.

<br>
[[File:ThirdYearSimulationExpt-Intro-Verlet-flowchart.svg|300px|thumb|center|'''Figure 2''': Steps to implement the classic Verlet algorithm.]]
<br>
'''<big> TASK 1: By taking Taylor expansions of  <math>x(t + \delta t)</math> and <math>x(t - \delta t)</math>, write general expressions for them up to the fourth order <math>\mathcal{O}\left(\delta t^4\right)</math> </big>'''

<br>
'''<big> Having written these expressions, derive the formula in figure 2 to update positions as used in the classical Verlet algorithm <math> x(t + \delta t) \approx 2x_{i}(t) - x_{i}(t-\delta t) + \frac{F_{i}(t)}{m} \delta t^{2}</math> by using Newton's second law to replace <math>\frac{\mathrm{d}^2\mathbf{x}_i\left(t\right)}{\mathrm{d}t^2}</math> [4 marks] </big>'''

Using this last equation, we can use a sequence of steps like those shown in '''figure 2''' to get the positions. At no point are the velocities calculated in this method!

====Velocity Verlet Algorithm====

If we assume that the acceleration of an atom depends only on its position and not its velocity, then we are able to come up with a new algorithm that lets us calculate atomic velocities explicitly as shown in '''figure 3'''.

<br>

[[File:ThirdYearSimulationExpt-Intro-VelocityVerlet-flowchart.svg|300px|thumb|center|'''Figure 3''': Steps to implement the velocity Verlet algorithm.]]

<br>

We start by noting that
<center><math>\mathbf{v}_i\left(t + \delta t\right) = \mathbf{v}_i\left(t\right) + \frac{\mathbf{a}_i\left(t\right) + \mathbf{a}_i\left(t + \delta t\right)}{2}\delta t \ \ (6)</math></center>

We can Taylor expand the velocity by half a step, instead of a full step.

<center><math>\mathbf{v}_i\left(t + \frac{1}{2}\delta t\right) = \mathbf{v}_i\left(t\right) + \frac{1}{2} \mathbf{a}_i\left(t\right)\delta t \ \ (7)</math></center>

We then substitute this into your expansion for <math> x_{i} (t + \delta t) </math> to obtain an accuracy up to <math> \delta t ^{2} </math>. Notice that terms up to <math> \delta t ^{2} </math> in your expansion can be written:

<center><math> \mathbf{x}_{i} (t + \delta t) = \mathbf{x}_{i} (t) + \mathbf{v}(t) \delta t + \frac{1}{2} \mathbf{a}(t) \delta t ^{2} </math></center>


<center><math>\mathbf{x}_i\left(t + \delta t\right) = \mathbf{x}_i\left(t\right) + \mathbf{v}_i\left(t + \frac{1}{2}\delta t\right)\delta t \ \ (8)</math></center>

When we know the updated atomic positions, we can calculate new forces, <math>\mathbf{a}_i\left(t + \delta t\right)</math>. Finally, we substitute equation (7) into equation (6) to get the new velocities <math>\mathbf{v}_i\left(t + \delta t\right)</math>.

<center><math>\mathbf{v}_i\left(t + \delta t\right) = \mathbf{v}_i\left(t + \frac{1}{2}\delta t\right) + \frac{1}{2}\mathbf{a}_i\left(t + \delta t\right)\delta t \ \ (9)</math></center>

Notice that for both numerical integration algorithms, the first step is "specify initial conditions". When using the Verlet algorithm, we need to know the starting positions of the atoms (<math>\mathbf{x}_i\left(0\right)</math>), and their positions one timestep in the past (<math>\mathbf{x}_i\left(-\delta t\right)</math>). If the velocity-Verlet algorithm is used, then we have to know the  the starting positions of the atoms (<math>\mathbf{x}_i\left(0\right)</math>) and their velocities at the same time (<math>\mathbf{v}_i\left(0\right)</math>). For this reason, we often start new simulations by using the output of older ones. If, however, you are performing your first simulations of a system (as we are now), then you must choose your initial conditions. The simulation software that we will use is able to do this for us, and this will be explained in the next section.

<br>
'''<big> TASK 2: What could be an advantage of the Velocity-Verlet algorithm over the classical Verlet algorithm? [2 marks]</big>'''
<br>

===Atomic Forces===

Since we can't reasonably solve the equations from quantum physics necessary to determine the forces acting on a given configuration of atoms, we have to make approximations. We know from classical physics that the force acting on an object is determined by the potential that it experiences:

<center><math>\mathbf{F}_i = - \frac{\mathrm{d}U\left(\mathbf{r}^N\right)}{\mathrm{d}\mathbf{r}_i}</math></center>

The shorthand notation <math>\mathbf{r}^N</math> stands for the position vectors of '''every''' atom in system. In principle, the force that a single atom feels is determined by the position of every other atom in the simulation. All we then need to do is to find a function <math>U</math> that captures all the key physics of the interatomic interactions in the system. For many simple liquids, it turns out that we can model the interactions between each pair of atoms extremely well using the Lennard-Jones potential. Overall, U takes the form:

<center><math>U\left(\mathbf{r}^N\right) = \sum_i^N \sum_{i \neq j}^{N} \left\{ 4\epsilon \left( \frac{\sigma^{12}}{r_{ij}^{12}} - \frac{\sigma^6}{r_{ij}^6} \right) \right\} \ \ (10)</math></center>

=== TASK 3: Consider the Lennard-Jones pair potential. What physical interaction(s) does it describe? What is the physical significance for the r^(-6) and r^(-12) terms? [3 marks] ===

'''<big>TASK 4: For a single Lennard-Jones interaction, <math>\phi\left(r\right) = 4\epsilon \left( \frac{\sigma^{12}}{r^{12}} - \frac{\sigma^6}{r^6} \right)</math>, find the separation, <math>r_0</math>, at which the potential energy is zero. What is the force at this separation? Find the equilibrium separation, <math>r_{eq}</math>, and work out the well depth (<math>\phi\left(r_{eq}\right)</math>). Evaluate the integrals <math>\int_{2\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, <math>\int_{2.5\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math>, and <math>\int_{3\sigma}^\infty \phi\left(r\right)\mathrm{d}r</math> when <math>\sigma = \epsilon = 1.0</math> [4 marks]</big>.'''

====Periodic Boundary Conditions====

[[File:ThirdYearSimulationExpt-Intro-Box.png|200px|thumb|right|'''Figure 4''': Diagram of a simulation box containing 2139 atoms. The blue lines indicate the boundaries of the box.]]
[[File:ThirdYearSimulationExpt-Intro-Periodic.svg|300px|thumb|left|'''Figure 5''': Periodic boundary conditions in two dimensions.]]

We cannot simulate realistic volumes of liquid. In fact, in our simulations, <math>N</math> will be between <math>1000</math> and <math>10000</math>. The following task should illustrate why this must be so.

In order for our simulations to approximate a bulk liquid, we have to use a computational trick. The atoms in the simulation are enclosed in a simulation box, of fixed dimensions ('''figure 4'''). This box is very often a cuboid, but parallelepipeds can also be used (and this can be very useful when simulating crystal structures). We pretend that we have repeated our box infinitely in all directions, so that the atoms at the very edges are not exposed to a vacuum. This is illustrated in two dimensions in '''figure 5'''. The darker coloured atoms in the central box are the "real" atoms. The faded atoms in the outer four boxes are the replicas. When an atom crosses the boundary of the box, one of its replicas enters the box through the opposite face. In this way, the number of atoms inside the box is always constant.

'''<big>TASK 5: Consider an atom at position <math>\left(0.5, 0.5, 0.5\right)</math> in a cubic simulation box which runs from <math>\left(0, 0, 0\right)</math> to <math>\left(1, 1, 1\right)</math>. In a single timestep, it moves along the vector <math>\left(0.7, 0.6, 0.2\right)</math>. At what point does it end up, ''after the periodic boundary conditions have been applied''? [1 marks]''' </big>.

====Truncation====

Periodic boundary conditions introduce their own problems. When we defined our potential function (equation 10), we specified that it depended on all possible pairs of atoms. If we have an infinite number of replicas of our system, how can we avoid calculating an infinite number of pair interactions?

Think about the three integrals you calculated for the Lennard-Jones potential task. They represent the area under the Lennard-Jones potential curve between some specified distance (<math>2\sigma</math>, <math>2.5\sigma</math>, and <math>3\sigma</math>), and infinite separation (where there is no interaction). You should find that this value becomes rather small as the near distance is increased! The attractive <math>\frac{1}{r^6}</math> part of the potential dominates here, and this decays rapidly with <math>r</math>. We assume that this means that there is a distance beyond which the interaction is so small that we can safely ignore it. In fact, in most simulations this is chosen to be something close to <math>2.5\sigma</math> or <math>3\sigma</math>. When the forces are calculated, we only calculate interactions between a pair of atoms if their separation is less than this cutoff.

====Reduced Units====

It is typical when using Lennard-Jones interactions to work in reduced units. By this, we mean that all quantities in our simulation are divided by scaling factors &mdash; for example, distances are divided by <math>\sigma</math>. The result of this is that the values become more manageable: all values that we might work out are typically around 1, rather than <math>1\times 10^{-10}</math> (in the case of distance), <math>300</math> (in the case of temperature), or <math>1\times 10^{-19}</math> (in the case of energy).

We denote these reduced quantities by a star, and they take the following conversion factors:

* distance <math>r^* = \frac{r}{\sigma}</math>
* energy <math>E^* = \frac{E}{\epsilon}</math>
* temperature <math>T^* = \frac{k_BT}{\epsilon}</math>

'''<big>TASK 6: The Lennard-Jones parameters for argon are <math>\sigma = 0.34\mathrm{nm}, \epsilon\ /\ k_B= 120 \mathrm{K}</math>. If the LJ cutoff is <math>r^* = 3.2</math>, what is it in real units? What is the well depth in <math>\mathrm{kJ\ mol}^{-1}</math>? What is the reduced temperature <math>T^* = 1.5</math> in real units? [1 marks]</big> '''

'''<big><span style="color:blue; ">This is the second section of the third year simulation experiment. You can return to the previous section, [[Third year simulation experiment/Running your first simulation|Files to Download]], or jump ahead to the next section, [[Third year simulation experiment/Equilibration|Equilibration]].</span></big>'''