Molecular Dynamics: Lab 1

In part based on Fortran code from Furio Ercolessi.

Basics of Molecular Dynamics

The basic idea of molecular dynamics is to follow the location of the particles in time using Newton's laws. Each particle has a location $\vec{X}_i$ which obeys

$$ \begin{equation} m_i \frac{d^2}{d t^2} \vec{X}_i = m_i \vec{A}_i = -\nabla V \left( \vec{X}_1, \dots, \vec{X}_N \right). \end{equation} $$

Most of the time we will work in computational coordinates, rescaled by a reference length $L$, so that the particles obey

$$ \begin{equation} \frac{d^2}{d t^2} \vec{x}_i = \vec{a}_i = -\nabla V \left( \vec{x}_1, \dots, \vec{x}_N \right). \end{equation} $$

This gives us a computational domain, once translated to the centre-of-mass, of $\vec{x} \in [-0.5, 0.5]^3$.

The simplest possibility is that the interaction potential $V$ is the sum of pairwise interactions,

$$ \begin{equation} V \left( \vec{x}_1, \dots, \vec{x}_N \right) = \sum_i \sum_{j>i} \phi \left( L \left| \vec{x}_i - \vec{x}_j \right| \right). \end{equation} $$

The most commonly used pairwise interaction potential is the Lennard-Jones potential

$$ \begin{equation} \phi(r) = 4 \varepsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^{6} \right]. \end{equation} $$

To stop this potential extending for infinite distance a cutoff is used; $\phi(r) \to \phi(r) - \phi(R_c)$, where $R_c$ is a constant (here $2.5$ when $\sigma = 1 = \varepsilon = M_i$).

As, in this case,

$$ \begin{equation} \frac{d}{dr} \phi(r) = 24 \left[ 2 \left( \frac{1}{r} \right)^{14} - \left( \frac{1}{r} \right)^{8} \right], \end{equation} $$

we have

$$ \begin{equation} \vec{a}_i = -\left. \nabla V \right|_{i} = \sum_{j>i} (\vec{x}_i - \vec{x}_j) \, \frac{d}{dr} \phi(r) - \sum_{j<i} (\vec{x}_i - \vec{x}_j) \, \frac{d}{dr} \phi(r). \end{equation} $$

Note that Newton's third law means that the force on particle $i$ from particle $j$ must be exactly opposite to the force on particle $j$ from particle $i$. So the algorithm can be implemented by setting $\vec{a}_i = 0$ and, for each particle $i$, and for all $j > i$, setting both

1. $\vec{a}_i = \vec{a}_i + (\vec{x}_i - \vec{x}_j) \, \frac{d}{dr} \phi(r)$,
2. $\vec{a}_j = \vec{a}_j - (\vec{x}_i - \vec{x}_j) \, \frac{d}{dr} \phi(r)$.

Implement an algorithm to calculate the acceleration given the positions.

Time evolution

The problem with time integrating is that most numerical algorithms do not conserve the total energy of the system, which can be be crucial. Certain algorithms, such as the velocity Verlet algorithm, do conserve total energy.

Velocity Verlet

Given initial positions $\vec{x}$ and velocities $\vec{v}$ of each particle, the algorithm can be written as

1. $x^{n+1} = \vec{x}^n + \Delta t \, \vec{v}^n + \tfrac{1}{2} \Delta t^2 \, \vec{a}^n$
2. $\vec{v}^{*} = \vec{v}^n + \tfrac{1}{2} \Delta t \, \vec{a}^n$;
3. Derive $\vec{a}^{n+1}$ from the interaction potential using $\vec{x}^{n+1}$;
4. $\vec{v}^{n+1} = \vec{v}^{*} + \tfrac{1}{2} \Delta t \vec{a}^{n+1}$.

Implement this algorithm.

Periodic boundaries

If we just evolved the particles, they could fly off to infinity. Instead it is typical to work inside a box with periodic boundaries, so that there is an "infinite lattice" of particles.

Implement a function that resets the position of particles that leave the given domain; the velocities and particles should remain unchanged.

What to compute

The fundamental quantity is the location of the particles. However, the thing to analyze will be average properties of the particles; their mean temperature, pressure and energies.

Implement a function, given the particle locations $\vec{x}$ to compute the temperature, where

$$ \begin{align} % E_{\text{potential},\ i} &= \sum_{j /ne i} \phi \left( L \left| \vec{x}_i - \vec{x}_j \right| \right), \\ E_{\text{kinetic},\ i} &= \tfrac{1}{2} L^2 \left| \vec{v}_i \right|^2, \\ T & = \frac{2}{3 N} \sum_{i=1}^N E_{\text{kinetic},\ i}. \end{align} $$

Run in a box

Fix the domain to be $\vec{x} \in [0, 6.1984]^3$. Fix the timestep to be $0.032$ and take 100 timesteps, given the initial data prescribed in input.dat (see data/MD/input.dat on github). Remember to rescale the initial coordinates so that

1. The computational domain is $[-0.5, 0.5]^3$;
2. The average particle location is the origin.

That is, given the data $\vec{X}$ from input.dat do

1. $\vec{x} = \vec{X} / L$ where $L = 6.1984$;
2. $\vec{x} \to \vec{x} - \text{average}(\vec{x})$.

Then plot some basic quantities, such as

1. The temperature as a function of time;
2. The trajectory of the first two particles as a function of time (using plot3D);
3. The trajectories of all the particles.