Molecular Dynamics and Statistical Mechanics

This is a simple MD simulation using 49 particle interacting through a repulsive 1/r^{12} potential. The equations of motion are integrated using the symplectic Stormer-Verlet method with periodic boundary conditions. Here are a few snapshots displaying particle locations as they evolve in time.

one snapshot per 0.001 time units,
one snapshot per 0.01 time units.

A hallmark of symplectic integration is their excellent conservation of energy. This can be explained via backward error analysis. See the work by Benettin et al, Hairer and Lubich, and myself. Here is an example of that behavior already noted by Verlet in 1967 when he did the first MD simulations of Argon:

computed energy.

Still, it is quite obvious that all accuracy is pretty quickly lost in a MD simulation due to its inherent chaotic nature. But accurate trajectory simulations never have been and will never be the goal of MD. Instead one wishes to extract information on the thermodynamic behavior of the molecular system. This gets us into statistical mechanics and the computation of time-averages. One thermodynamic variable of interest is temperature which can be defined at any moment in time as the average kinetic energy over all particles. Unless an infinite number of particles is taken, this instanteneous temperature will fluctuate. To demonstrate this, we take two sets of initial conditions with exactly the same energy and plot the instanteneous energy along the numerically computed trajectories:

instanteneous temperature I,
instanteneous temperature II.

It appears that the instanteneous temperature fluctuates about a mean value. Indeed, under the assumption of ergodicity, this mean value exists and is independent of the initial data with identical energy. Can we reproduce this fact numerically?

time averaged instanteneous temperature I,
time averaged instanteneous temperature II.

But can we trust all these numbers? Well, under the assumption that the system is hyperbolic, we indeed can. The numerically computed time averages are to high probability within a certain distance to the exact time-average. For a requested (fixed) error tolerance, the probability increases exponentially in the simulation time. This can be shown combining backward error analysis, shadowing, and a large deviation theorem. See my 1999 paper in SIAM J. Numer. Anal.

The simulations were performed on a HP laptop using MATLAB.