The Lennard-Jones potential

Molecular dynamics (MD) is extremely widely used in computational chemistry, physics, biology, material science etc. It was also one of the first scientific uses of computers going back to work by Alder and Wainwright who, in 1957, simulated the motions of about 150 argon atoms. (See Alder, B. J. and Wainwright, T. E. J. Chem. Phys. 27, 1208 (1957) and J. Chem. Phys. 31, 459 (1959))

In MD each atom is represented as a point in space. The atoms interact with each other via pairwise potentials, in a similar manner to gravity interacting between planets, sun, moons and stars. In practise this is an approximation since for atoms quantum effects are potentially important - but we'll ignore this and confine our discussions to ``classical'' MD.

Gravitational bodies interact via a potential that goes as $\frac{m_1m_2}{r_{12}}$ where $m_1$ and $m_2$ are the masses of the bodies involved and $r_{12}$ is the distance between their centres of mass. The potential that we will use is called a Lennard-Jones potential. This is in the form:

$v_{12} = \frac{1}{r_{12}^{12}} - \frac{2}{r_{12}^6}$

To evaluate the interaction for many bodies we must simply perform a sum over all unique two body interactions. Thus for 3 particles we will have interactions 1-2, 1-3 and 2-3.

Given a pairwise potential there will be a pairwise force. This corresponds to the derivative of the potential with respect to $r_{12}$. Assuming the location of each atom is given in Cartesian coordinates the Cartesian forces are given by:

$v_{12} = (\frac{12}{r_{12}^{14}} - \frac{12}{r_{12}^8})(x_1 - x_2)$

The total force on a given atom will be the sum of all the pairwise forces.

Implementation of Lennard-Jones

To find the total potential of the system, we use the following algorithm:

<Evaluate Potential>=
def evalPotential(atoms):
    # Evaluate the total potential energy
    potential = 0
    for i in range(len(atoms)):
        for j in range(i):
            dist2 = ((atoms[i][0] - atoms[j][0]) ** 2 +
                     (atoms[i][1] - atoms[j][1]) ** 2 +
                     (atoms[i][2] - atoms[j][2]) ** 2)
            dist6 = dist2 ** 3
            dist12 = dist6 ** 2
            potential = potential + ((1.0 / dist12) - (2.0 / dist6))
    return potential

Used below.

To find the force on each atom, we use the following algorithm:

<Evaluate Force>=
def evalForce(atoms):
    # Evaluate the force on each atom
    force = []
    for i in range(len(atoms)):
        force.append((0.0, 0.0, 0.0))
    for i in range(len(atoms)):
        for j in range(i):
            dist2 = ((atoms[i][0] - atoms[j][0]) ** 2 +
                     (atoms[i][1] - atoms[j][1]) ** 2 +
                     (atoms[i][2] - atoms[j][2]) ** 2)
            dist6 = dist2 ** 3
            dist8 = dist2 * dist6
            dist14 = dist8 * dist6
            a = 12.0 / dist14 - 12.0 / dist8
            b0 = (a * (atoms[i][0] - atoms[j][0]))
            b1 = (a * (atoms[i][1] - atoms[j][1]))
            b2 = (a * (atoms[i][2] - atoms[j][2]))
            force[i] = (force[i][0] + b0, force[i][1] + b1, force[i][2] + b2)
            force[j] = (force[j][0] - b0, force[j][1] - b1, force[j][2] - b2)
    return force

Used below (1), below (2).