PRACTICAL Implementation of Classical Molecular Dynamics Simulations using Python: (I) Initialization, Boundary conditions, Interaction Potentials, and Thermostat

Classical molecular dynamics is a computational simulation technique used to model structural and dynamical behaviors of particles over time. It is based on the principles of classical mechanics and the Newton's equations of motion to describe the time "evolution" of a system with interacting particles.

In this tutorial, we provide a practical implementation of classical molecular dynamics simulations using Python. This simulation involves the initialization of modelling systems, the selection of boundary conditions, the calculation of interaction potentials (energies and forces) between interacting particles, and their integrations over time so as to animate their dynamical evolution in the modelling system.

Contents of this documents and quicklinks：

• 1. Initialization
  • 1.1 Initialization of particle positions
  • 1.2 Animation of simulation trajectories
• 2. Boundary Conditions
  • 2.1 Motion of particles
  • 2.2 Periodic boundary condition
  • 2.3 Reflective boundary condition
• 3. Interaction Potentials and Thermostat
  • 3.1 The Lennard-Jones potential
  • 3.2 Implementation of Lennard-Jones potential
  • 3.3 Thermostat
• 4. Summary

1. Initialization

1.1 Initialization of particle positions

A classical molecular dynamics simulation starts with the initialization of positions for particles. For simplicity, we construct a two dimensional modelling system consisting of 16 particles in a simulation box with the box size of 10*10.

numAtoms = 16
boxLength = 10

Then we generate random positions for these particles using numpy,

import numpy as np
positions = np.random.rand(numAtoms, 2)*boxLength

in which the positions is a list containing numAtoms items with the X and Y coordinates of each particle being the first and the second one in each item. In order to ensure that everyone can get the same sequence of random numbers when running this script, we set a specific seed value through

np.random.seed(numAtoms)

, which is particularly useful for debugging and testing purposes, as it allows you to create predictable and repeatable results.

For molecules with complicated structures, it is feasible to read coordination info from an external file instead of generating random positions for molecules.

1.2 Animation of simulation trajectories

Data visualization and animation are powerful tools for exploring and presenting data in a visual and dynamic way. For classical molecular dynamics simulations, the animation of simulation trajectories provides a vivid representation of the movement of simulated particles within the modelling system as time evolves.

Matplotlib is a comprehensive and the most widely used library for creating static, animated, and interactive visualizations using Python. Celluloid is a Python library that simplifies the process of creating animations by automatically capturing each configuration in simulation trajectories as Matplotlib plots are generated, eliminating the need to manually manage configuration generation.

Herein, we use matplotlib and celluloid to animate the simulation trajectories of the modelling system with interacting particles.

The code example for the animation part is provided below:

from matplotlib import pyplot as plt
from celluloid import Camera

fig = plt.figure()
plt.xlim([0, boxLength])
plt.ylim([0, boxLength])
camera = Camera(fig)

colorsList = random.sample(range(1, numAtoms+1), k = len(positions) )
colors = colorsList

plt.scatter(positions[:, 0], positions[:, 1], s=256, c=colors, cmap='plasma')
camera.snap()

animation = camera.animate()
animation.save('simu_system.gif', writer = 'pillow', fps=64)

Combining all these code together into a file named ex1-initialization.py, we run this script and obtain the following figure.

---

2. Boundary Conditions

2.1 Motion of particles

In addition to positions, we will generate initial velocities for all particles. This is not necessary as in real modelling systems we can obtain these information via integrating particle positions over time. Herein for simplicity we generate a set of random velocities for all particles.

velocities = np.random.rand(numAtoms, 2) - 0.5

Here we assume that there is NO interactions between particles, meaning that the motion of all particles at all directions are independent.

We set all particles to run 100 steps with a timestep of 0.5, and we will calibrate these parameters when introducing appropriate interaction potentials.

iterations = 100
timestep = 0.5

A simple function is introduced to update positions of particles as time evolves.

def update_positions(positions, velocities, timeStep):
	for atom in range(numAtoms):
		for col in [0, 1]:
			positions[atom, col] += timeStep * velocities[atom, col]

Accordingly, a for loop will be provided in the animation part to capture instantaneous positions of all particles at each step.

for istep in range(numSteps):
	plt.scatter(positions[:, 0], positions[:, 1], s=256, c=colors, cmap='plasma')
	camera.snap()
	update_positions(positions, velocities, timeStep)

These code are assembled into a new script ex2-boundary-conditions-1-particle-motion.py to obtained the following animation.

From the above animation, it is shown that all particles gradually fly out of the simulation box. This is not what we expect as we want to track the motion of all particles within the simulation box. As such, we will introduce proper boundary conditions for the modelling system.

Boundary condition specifies the behavior of a modelling system at its boundaries or interfaces with the external environment. Boundary condition is crucial for accurately modeling the behavior of physical systems and ensuring that simulations produce meaningful results.

There are several types of boundary conditions in classical molecular dynamics simulations:

• Open Boundary Condition (OBC): An OBC allows particles to freely diffuse in or out of simulation box (as shown above) without any restriction at boundaries.
• Periodic Boundary Condition (PBC): A PBC describes an infinite system by replicating the simulation box in all directions. Particles leaving one side of the simulation box re-enter from the opposite side. PBC applies to simulations with a finite box size but can achieve an infinite system.
• Reflective Boundary Condition (RBC): In a RBC, particles are reflected back into the simulation box when they reach boundaries. This condition is commonly used in simulations of fluid particles near solid walls.

Besides these boundary conditions, there are some others that are widely used to simulate particular modelling systems. A comprehensive discussion of boundary conditions is available HERE.

2.2 Periodic boundary condition

PBC applies to simulations with a finite box size but can achieve an infinite system via replicating the system box in all directions. With PBC, particles leaving one side of the simulation box re-enter from the opposite side, which is achieved via the following code.

if positions[atom, col] > boxLength:
	positions[atom, col] -= boxLength
if positions[atom, col] < 0:
	positions[atom, col] += boxLength

This code is implemented into a new function boundary_conditions.
By passing periodic to a variable boundaryCondition, we run the script ex2-boundary-conditions-2-boundary.py and obtain the follow animation.

2.3 Reflective boundary condition

In a reflective boundary condition, particles are reflected back into simulation box when they reach boundaries, mimicking the behavior of fluid particles when they hit solid walls. Meanwhile, the velocities of particles at specific directions will be reversed. This condition can be achieved via the code below.

if positions[atom, col] > boxLength:
	positions[atom, col] -= 2.0*(positions[atom, col] - boxLength)
	velocities[atom, col] *= -1.0
if positions[atom, col] < 0:
	positions[atom, col] += 2.0*(0.0 - positions[atom, col])
	velocities[atom, col] *= -1.0

Passing reflective to the variable boundaryCondition, we run the script ex2-boundary-conditions-2-boundary.py and obtain the follow animation.

---

3. Interaction Potentials and Thermostat

In above code examples, there is no interactions between particles, indicating that all particles have independent diffusion at all directions. In addition, particles can also pass through others during diffusion, and can even have complete overlap with the other particles.

In real modelling simulations, we can include interaction potentials to simulate particles with particular characteristics.

3.1 The Lennard-Jones potential

Particle interaction potentials describe the energy and force associated with specific interactions between particles in simulation systems. These potentials are fundamental in determining structural and dynamical quantities of interacting particles in modelling systems.

In this tutorial, we implement the Lennard-Jones potential in the code example, which is defined by the following equation:

$$U(r)=4ϵ[(\frac{\sigma }{r}{)}^{12}-(\frac{\sigma }{r}{)}^6],$$
where

$r$ is the distance between interacting particles,

$ϵ$ is the depth of the potential well usually referred to as the strength of attractive interaction between particles, and

$\sigma$ is the finite distance at which the inter-particle potential

$U$ is zero, representing the distance at which the attractive and repulsive forces balance. The Lennard-Jones potential has its minimum at a distance of

$r=r_{min}=2^{1/6}\sigma$, where the potential energy has the value

$U=-ϵ$.

The code example for the Lennard-Jones energy and force (the derivative of the potential) with reduced units of

$ϵ=1.0$ and

$\sigma =1.0$ are available in the script ex3-potential-1-lj.py and representative plots are shown below.

3.2 Implementation of Lennard-Jones potential

Before the code implementation, we define parameters representing particle masses, interaction energies, forces, and accelerations acting on all particles.

masses = np.ones(numAtoms)
energies = np.zeros(numAtoms)
forces = np.zeros((numAtoms, 2))
accelerations = np.zeros((numAtoms, 2))

It should be noted that the Lennard-Jones potential has a singularity when

$r=0$, leading to an infinite force on particles. This singularity can cause numerical instabilities and simulation system crashes when particles get too close to each other. To avoid these issues, we generate particle positions on regular lattices instead of at random positions.

# 16 particles are generated on regular 4x4 lattices
positions = np.zeros((numAtoms, 2))
numLatticesX = round(math.sqrt(numAtoms))
numLatticesY = math.ceil(numAtoms/numLatticesX)
for iatom in range(numAtoms):
	positions[iatom, 0] = (iatom%numLatticesX+0.5)*(boxLength/numLatticesX) + np.random.rand()*0.5
	positions[iatom, 1] = (math.floor(iatom/numLatticesX)+0.5)*(boxLength/numLatticesY) + np.random.rand()*0.5

Additional note should be emphasized that when we use PBC in simulations, a particle moves beyond one side of the simulation box, it re-enters into the box from the opposite side as if the simulation box was replicated infinitely in all directions. This concept is often illustrated as "periodic images" or "image cells" to represent periodic replicas of the central simulation box. During the calculation of interaction energies and forces, we should determine the shortest distance between two interacting particles and their periodic images using the minimum image convention, which can be achieved via the code below.

for col in [0, 1]:
    dist_col = positions[iatom, col] - positions[jatom, col]
    if dist_col > boxLength/2.0:
        dist_col -= boxLength
    if dist_col < -boxLength/2.0:
        dist_col += boxLength

Combining the code mentioned above and the ex3-potential-1-lj.py together into a new function lj_potential in ex3-potential-2-potential-thermostat.py, we can determine the specific energies and forces acting on all interacting particles.

epsilon = 1.0  # energy in units of k_BT
sigma = 1.0    # distance in units of sigma
def lj_potential(energies, forces, positions, epsilon, sigma):
	for iatom in range(numAtoms-1):
		for jatom in range(iatom+1, numAtoms):
			dist_square = 0.0 # square of distance
			for col in [0, 1]:
				dist_col = positions[iatom, col] - positions[jatom, col]
				if dist_col > boxLength/2.0:
					dist_col -= boxLength
				if dist_col < -boxLength/2.0:
					dist_col += boxLength
				dist_square += dist_col**2
			dist = math.sqrt(dist_square)

			energy = 4.0*epsilon*((sigma/dist)**12 - (sigma/dist)**6)
			energies[iatom] += 0.5*energy # Accumulate energy
			energies[jatom] += 0.5*energy

			force_repulsion = 48.0*epsilon*(sigma/dist)**12
			force_attraction = -24.0*epsilon*(sigma/dist)**6
			force = force_repulsion + force_attraction

			for col in [0, 1]:
				dist_col = positions[iatom, col] - positions[jatom, col]
				if dist_col > boxLength/2.0:
					dist_col -= boxLength
				if dist_col < -boxLength/2.0:
					dist_col += boxLength
				forces[iatom, col] += (force*dist_col/dist_square)
				forces[jatom, col] -= (force*dist_col/dist_square)

3.3 Thermostat

Once we obtain the Lennard-Jones forces acting on all particles, we can calculate the accelerations on all particles, and therefore to update their velocities via a new function update_velocities.

def update_velocities(forces, masses, accelerations, velocities, timeStep):
	for iatom in range(numAtoms):
		for col in [0, 1]:
			accelerations[iatom, col] = forces[iatom, col]/masses[iatom]
			velocities[iatom, col] += accelerations[iatom, col]*timeStep

In classical molecular dynamics simulations, a thermostat is usually used to control the temperature of simulation system so as to reproduce statistical properties of certain thermodynamic ensembles. The velocity rescaling thermostat is a common method used to control the temperature of simulation systems via rescaling velocities of all particles at regular intervals to achieve a desired temperature. Below is the code to implement this method through a new function thermostat.

temperatureRef = 1.0
def thermostat(velocities, masses, temperatureRef):
	kineticEnergy = 0.0
	for iatom in range(numAtoms):
		ke_temp = 0.0
		for col in [0, 1]:
			ke_temp += velocities[iatom, col]**2
		ke_temp *= 0.5*masses[iatom]
		kineticEnergy += ke_temp
	temperatureTemp = kineticEnergy*2.0/(2.0*numAtoms)

	scaling_factor = math.sqrt(temperatureRef/temperatureTemp)
	for iatom in range(numAtoms):
		for col in [0, 1]:
			velocities[iatom, col] *= scaling_factor

Based on these implementations, we call these newly added functions from the iteration loop in the following code snippets.

for istep in range(numSteps):
	plt.scatter(positions[:, 0], positions[:, 1], s=256, c=colors, cmap='plasma')
	camera.snap()
	update_positions(positions, velocities, timeStep)
	boundary_conditions(positions, boxLength, boundaryCondition)
	lj_potential(energies, forces, positions, epsilon, sigma)
	update_velocities(forces, masses, accelerations, velocities, timeStep)
	thermostat(velocities, masses, temperatureRef)

In addition, a small timestep will be used in order to keep the simulation system relatively stable.

numSteps = 150
timeStep = 0.01

This is one of the first, of many, tradeoffs of accuracy vs. computational efficiency for classical molecular dynamics simulations. In fact using an extremely small time-step may lead to accurate simulation results, but at a cost of longer computation time or reduced total simulation time which can reduce statistical reliability. A larger time-step may help assist in larger or longer time-scales but at the cost of physical accuracy. It should be emphasized that such a tradeoff between accuracy and computational time is NOT unique to classical molecular dynamics, and can be found in practically every computational modeling technique.

Combining all code examples into the Python script ex3-potential-2-potential-thermostat.py, we run this Python script and get the following animation.

---

4. Summary

In this tutorial we have provided practical implementation of classical molecular dynamics simulations using Python focusing on

• initialization of modelling system
• application of periodic and reflective boundary conditions
• calculation of interaction energies and forces between particles using the Lennard-Jones potential
• implementation of the velocity rescaling thermostat to achieve certain thermodynamic ensembles

In the next tutorial, we will explore how to use additional simulation techniques to accelerate simulations of modelling systems containing thousands of particles.

---