Molecular Dynamics - N-body

Molecular Dynamics aims at describing the properties of an assembly of molecules in term of their structure and the microscopic interactions between them. It is an extremely powerful numerical approach to describe many-body systems.

Examples of Molecular Dynamics applications abound. For instance, bio-molecular simulations, complex fluids models, crystal models, descriptions of defects in crystals, models of surfaces ... all use variations of Molecular Dynamics. Each of these applications comes with its own set of issues.

The folding of proteins into their native geometrical configurations are often simulated by Molecular Dynamics. The model may assume a chain of beads, the amino acids sequence, surrounded by hundreds or thousand of water molecules. The interactions between the beads themselves and between the beads and the water molecules are adjusted phenomenologically. Newton equations, defining a system of ODEs, are solved numerically. The simulation can be done at constant energy, or at constant temperature using statistical physics concepts. The rotations of the particles (amino acids and water molecules) can be incorporated in the simulation with interactions that depend on the orientation of the molecules in the system. Boundary conditions can be reflexive (to model the containment of a vessel) or periodic to imitate an infinite system and reduce surface effects.

In Ab Initio Molecular Dynamics calculations the molecules are broken into individual atoms. In the Protein Folding problem, the amino acids can be described in terms of individual atoms. The same could be done for the water molecules. Incorporating the individual atoms and chemical bonds leads to more realism. Quantum Mechanics should eventually be used (at an extremely high computational cost) for ultimate reality.

A Simple Code

Learning about Molecular Dynamics is best achieved by looking at an example. The code Nbody_basic_verlet.c provides such an example. It solves the problem of N particles that interact via a two-body potential of finite range. Newton equations are solved. The masses are allowed variability. The 2-body interaction is of Lennard-Jones type, the 6-12 potential.

$$V(r) = 4 \epsilon \left( \frac{\sigma}{r}^{12} - \frac{\sigma}{r}^6 \right)$$

where $r$ is the inter-particle distance. This potential is described in Lennard-Jones potential. The Maple worksheet LJ_potential.mw analyzes this potential form. The constants $\sigma$ and $\epsilon$ define the geometry and strength of the interaction. They are both set to 1 in the code.

The code implements the following steps:

Step #1:

Set the model parameters, in particular the masses of the particles. The masses are randomly distributed with a Gaussian profile centered on a mass of 1 scaled unit.

Step #2:

The initial positions of the particles are set on a compact square lattice with cell spacing of 1.2 units, a distance which roughly corresponds to the minimum of the potential. This choice is reasonable given that the strong short distance repulsion of the potential should be avoided.

Step #3:

The initial velocities of the particles are set at random (uniform distribution) in the interval [-1,1] in scaled units.

Step #4:

The center of mass of the system is forced to be stationary by imposing its momentum to be zero.

Step #5:

Newton equations are based on the two-body interaction and the particle mass. The system of ODEs that follows is of dimension 6N ( $x$ , $y$ , $z$ , $v_x$ , $v_y$ and $v_z$ for each particle). The RHS of the velocity ODEs are the forces on each particle divided by its mass.

Step #6:

The Basic Verlet ODE integration schemes is used to solve the ODEs in order to maintain the total energy.

The two-body potential dictates the RHS of the ODEs. Recall that $f = - \frac{\partial U}{\partial r}$ . This yields

$$F_{x,i} = \sum_{j} 4 \left( \frac{12}{r_{ij}^{14}} - \frac{6}{r_{ij}^8} \right) x_{ij}$$

where $r_{ij}$ refer to the distance between particles $i$ and $j$ and $x_{ij}$ refers to its $x$ component. The same form applies to the $y$ ans $z$ components.

Note that the RHS of the ODEs are expensive to compute!

The Basic Verlet ODE integration scheme is used in the code to ensure constancy of the conserved quantities, i.e., energy and momentum. The interpretation of the results and the analysis of the solution rely on maintaining the conserved quantities.

Exercises

Exercise #1:

Run the code for $N=2$ . Plot the trajectories in $x$ , $y$ and $z$ . Plot the kinetic, potential, and total energy versus $t$ . Note the (slight lack of) constancy of the latter.

Exercise #2:

Run the code for $N=3$ . Plot the trajectories in $x$ , $y$ and $z$ . Plot the kinetic, potential, and total energy versus $t$ . Note the (slight lack of) constancy of the latter.

The trajectory definitely looks different than for $N=2$ . Poincare said that N-body systems with $N=3$ and larger are non-integrable (chaotic).

Exercise #3:

Repeat for $N = 4, 5, 6, ...$ Careful - the trajectory recording file can become extremely large!

Periodic Boundary Conditions

Some particles in the simulations for $N>2$ may take off. They will boil-off following a close encounter (collision) with other particles through a boomerang effect. For instance, Molecular Dynamics simulations of proteins will lead to an explosion of the molecules unless the relative amino acids positions are constrained externally or water molecules are included.

Boundary conditions can also constrain the particles in a finite volume. Reflexive surfaces model a confinement vessel. This approach should be used in gas models where the pressure needs to be calculated.

Periodic Boundary Conditions (PBC) model an infinite system via a finite volume. This is achieved by a cell tiling of the space while including one cell in the calculation. When a molecule escapes through a face of the basic cell, its mirror molecule reappears on the opposite face with the same velocity. This model is equivalent to a basic cell surrounded by image cells.

Fig.: Periodic boundary conditions. The central box is outlined by a thicker line. (Ref: PBC) )

Code Improvements

The code 'Nbody_basic_verlet.c is to be considered at most as a skeleton of a Molecular Dynamics code. Many aspects of it should be improved. In particular:

Exercises

Exercise #4:

Replace the ODE solver by a Velocity Verlet scheme.

Exercise #5:

Run the code for $N=2$ and $N=3$ . For the latter compare $x(t)$ , $y(t)$ , $z(t)$ and the energies $K(t)$ and $U(t)$ calculted by this more accurate code and the original one. Note the divergence of the curves. Using different solution schemes (leading to different accuracy) is equivalent to using different initial conditions in a chaotic regime.

Exercise #6:

Compute the averages of the kinetic energy $K(t)$ and potential energy $U(t)$ (beyond some transient times), say over t=[5:20]. Do that for the original code and the improved code. The expected values (time averages) of global variables should be mostly independent of the details of the results, thereby illustrating the validity of statistical physics.

Exercise #7:

Implement the Periodic Boundary Conditions in the improved code.

Exercise #8:

Run the code for a larger value of $N$ for which a particle was originally escaping. Plot trajectories and variables ($x(t)$ , ...) as a check.

Exercise #9:

Plot $U(t)$ , $K(t)$ and $e(t)$ for the case above. Note the gap(s) in $U(t)$ . Explain.

Exercise #10:

Run the code for large $N$ ( 9, 64, 125, 216, ... ) and plot the energies. Be careful to switch off the trajectory recording not to fill the disk! This option should be controlled by a line argument switch.

WEB references

Wikipedia: MD

MD primer

Wikipedia: Periodic Boundary Conditions