Molecular dynamics simulation: A tool to view a system at the atomic level

Computer simulation mimics the behaviour of a real system. Simulation can verify a theoretical model which cannot be solved analytically. Simulations can also complement experiments to get a detailed view of a system. In some cases, experiments are impossible to perform, modeling and simulation is the only choice. The first and crucial step towards simulations is to develop a model for the system. The more accurate the model, the closer the simulation to a real system. It is essential to validate a model by comparing simulation results with the experimental results. The comparisons between experimental, simulation and theoretical results help build a better model. Once we have a good model of a system, we can explore a wide range of properties and characteristics of the system.

Molecular dynamics (MD) technique, which is based on the principle of classical mechanics is used to simulate a wide variety of systems, from a single crystal to a macromolecule such as DNA, RNA. For an ergodic system, the time average of a dynamical quantity is equal to its ensemble average. So, the MD technique may be used to calculate the macroscopic properties of an ergodic system.  

In all-atom molecular dynamics simulation, the basic units are atoms. Atoms bond together to make a molecule and the system has many such molecules. In coarse-grain molecular dynamics simulation, the basic units are beads. A bead consists of a group of atoms and a group of beads bond together to make a molecule.  We apply Newton's equation of motion to all the particles (atoms or beads) for some fixed interval of time, dt assuming the forces are constant during this interval. Then, we calculate the forces from the current coordinates of the particles. We apply Newton's equation of motion for the next time interval with the updated forces. We keep on repeating this procedure to find the trajectory of the system. This technique of simulation is known as molecular dynamics simulation. The molecular dynamics simulation technique is discussed briefly.

1. Initial configuration:

We need a realistic initial configuration to start with. If our initial configuration is far from the equilibrium configuration, the simulation may get stuck at some local minimum configuration. Apart from the initial coordinates of the particles, we need to provide a set of initial velocities. Identically, we can set an initial temperature and then, assign the initial velocities of the particles according to the Maxwell-Boltzmann distribution.

2. Interaction potential functions:

(a) Lennar-Jones potential:

If we have a system of mono-atomic uncharged molecules, the only relevant interaction potential is Lennard-Jones (LJ) Potential, which is a function of inter-particle distance. The LJ potential has a smooth and long-range attractive part (dipole-dipole, dipole-induced dipole and London dispersion interaction) and a steep and short-range repulsive part (Pauli repulsion). The attractive term varies inversely as the sixth power of distance whereas the repulsive term varies inversely as the twelfth power of distance. The LJ potential is truncated and shifted to save the computation time and satisfy the minimum image convention.

(b) Coulomb potential:

For a system of mono-atomic charged particles, we need Coulomb potential in addition to the LJ potential. The Onsager reaction field model is used to treat the long-range behaviour of Coulomb potential. In this model, the electrostatic interactions are calculated explicitly within a spherical cavity of finite radius and the medium is considered homogeneous and polarizable outside the cavity. This way the infinite Coulomb potential is replaced by a sum of a finite term and the reaction field. Both the LJ and the Coulomb potential are non-bonded pair-potential.

(c) Bond-potential:

For a system of diatomic charged molecules, we need bond-potential in addition to the LJ and Coulomb potential. The harmonic bond-potential is based on the fact there is an equilibrium bond length between two atoms, where the system energy is the minimum. The potential is proportional to the square of deviation of bond length from the equilibrium value. If the bond length increases or decreases, a restoring force arises to bring back the bond length to its equilibrium value. The harmonic bond-potential is not accurate when the bond deviates too much from the equilibrium bond length. The Morse potential, which includes the higher-order terms is often used to represent a more realistic bond in molecular dynamics simulation.

(d) Angle-potential:

For a system of triatomic molecules, we need angle-potential as well. The harmonic angle-potential is similar to the bond potential and based on the fact that there is an equilibrium angle between two bonds, where the system energy is minimum. The potential is proportional to the square of deviation of angle from its equilibrium value. The cosine harmonic angle potential and the Urey-Bradley potential are other two commonly used angle-potentials. 

(e) Torsion potential:

For molecules having more than three atoms, we may need to consider torsion potential in addition to the above-mentioned potentials. The torsion potentials are of two types - proper dihedral and improper dihedral potential. A proper dihedral angle depends on the four consecutive bonded atoms, whereas the improper dihedral angle is given by three atoms centered around a fourth atom. The proper dihedrals which are mostly used to constrain the rotation around a bond are of two types - periodic type and Ryckaert-Bellemans type. The improper dihedrals which keep planar groups planar are of two types - harmonic type and periodic type.

(g) Restraint and other special potentials:

Restraints are imposed on the motion, either to avoid disastrous deviations or to include knowledge from experimental data. There might be other special types of non-bonded and bonded potentials in some system depending on the type of interaction among the particles. Choosing the right potential functions is essential for an accurate representation of the system. 

3. Numerical integration to solve equations of motion:

The intermolecular and intramolecular interactions are functions of positions and orientations. So, once we have the interaction potentials and the system configuration, we can calculate forces on each and every atom. Once the forces for the configuration at the current time (t) are obtained, the co-ordinates and velocities for a later time, t+dt can be calculated by solving equations of motions. Various algorithms are available for solving equations of motion. Leap-frog integrator and Velocity Verlet integrator are two commonly used algorithms.  

4. Periodic Boundary Condition (PBC): 

The model system in MD simulation is very small compared to the real macroscopic system. The artefacts arise due to edge effects in this finite system and the model system cannot represent a bulk material. To avoid this problem, periodic boundary conditions are imposed. The non-bonded pair interaction in a periodic system is dealt with a cut-off distance and the minimum image convention method.

5. Treating long-range interactions:

The long-range interactions are computationally very expensive. To treat them efficiently, a twin-range method is employed. In this method, two cut-offs are specified. The long cut-off is used to prepare a non-bonded neighbor list whereas the short cut-off is used to calculate the interactions. The cut-off and minimum image methods are not always accurate. The electrostatic interactions in a periodic system are often treated by more accurate lattice sum methods such as Ewald summation, Particle-Mesh Ewald (PME) and Particle-Particle Particle Mesh (PPPM). They perform a full summation of the forces between the infinite periodic images of all atoms.

6. Temperature and Pressure Coupling:

We may need to control temperature and pressure of the system during a simulation. The Berendsen algorithm mimics the coupling of the system to an external heat bath to maintain a fixed temperature. In Nose and Hoover scheme, the thermal reservoir is considered to be an integral part of the system and it is represented by an additional degree of freedom. A similar Nose and Hoover scheme is available to maintain pressure during a simulation. The Parinello-Rahman algorithm is another commonly used pressure coupling method.

7. Calculation of statistical quantities:

Once we have the simulation trajectories of a system, we can extract its statistical information. The kinetic energy, potential energy, pressure tensor, coefficient of diffusivity, radial distribution function etc. can be calculated from a trajectory which keeps track of the coordinates and velocities of each and every particle with time.

The molecular dynamics simulation is computationally very expensive. It requires High Performance Computing (HPC) facility to run a molecular dynamics simulation. GROMACS, a popular molecular dynamics package by the University of Groningen, Netherlands supports CPU acceleration (SSE, AVX etc.), GPU acceleration, parallelization based on MPI, multithreading with thread-MPI and other performance enhancement techniques. Some other popular molecular dynamics packages, CHARMM, AMBER, LAMMPS, NAMD, OpenMM also use various acceleration and parallelization techniques.

Dr. U. K. Basak,

Professor

Techno India University.

Ex Senior Research Fellow, Saha Institute of Nuclear Physics.