3.3.0 Molecular Dynamics Overview

Starting from a (more or less arbitrary) point in phase space (p, q), the time evolution of a trajectory of a system of particles is followed, usually by applying the classical mechanical laws of motion.

For example, given N identical sphere interacting through pair potentials which depend only on the distance r_ij separating particles i and j:

Pair force:

Newton's equation of motion:

6N coupled first order differential equations.

Method: the 6N coupled equations are supplemented with 6N initial conditions (for example, 3N coordinates at one time and 3N coordinates at a later time, or 3N coordinates and 3N velocities at a common time) and solved numerically. Common algorithms are the Verlet algorithms or Predictor-Corrector algorithms (e.g. Euler-type predictor with trapezoid corrector).

Properties are then computed from samples along the trajectory:

for simple thermodynamic averages, like <U>, <V> etc. The total internal energy is the sum of the potential and kinetic contributions:

<tau> is a "configuration temperature".

Pressure (if volume is allowed to vary during the simulation):

Pressure (if volume is held constant during the simulation):

Note that for small molecules it is conventional to associate pressure with intermolecular forces, where f_ij is the force between molecules i and j. For polymers, this distinction becomes vague, and it is important to distinguish between the "molecular virial" and the "atomic virial".

Combinatorial properties are implicit in the MD ensemble, but may be obtained relative to a known thermodynamic state by the use of thermodynamic integration:

Other quantities may be found from fluctuations in quantities like energy, volume, etc. For more detail, refer to Allen and Tildesley, Chapter 2.

Molecular Dynamics, by nature, samples phase space using equations of motion which are energy-conserving, that is, constant NVE. Methods have been proposed to "rescale" the equations of motion so that other ensembles are samples (e.g. for NVT, see S. Nose, J. Chem. Phys., 81, 571, 1984; for NPT see M. Parinello and A. Rahman, J. Appl. Phys., 52, 7182, 1981.)

However, one of the real advantages of MD is that it allows estimation of dynamic quantities (i.e. transport coefficients) in NVE.

3.3.1 Transport Coefficients

Einstein Relation (ca. 1900)

Green Kubo Relation (ca. 1955)

In this manner, any molecular property F can be related to the corresponding transport coefficient b:

coefficient of self diffusion

coefficient of thermal conductivity

static electrical conductivity

3.3.2 The Correlation and Autocorrelation Functions

The quantity <F(t).F(0)> is a special case of the general formula for the Correlation Function <B(t)A(0)>. If B=A, this is called the Autocorrelation Function. In the limit of very long time, B and A become uncorrelated, so that the limiting value of the correlation function is the product of the averages:

It is common to subtract the limiting value from the correlation function and then to normalize this function to get a function which decays from 1 to 0. A related function is the Fluctuation Correlation:

3.3.3 Issues in MD

In systems with simple, non-continuous potentials ("hard body" systems), acceleration only occurs at collisions, so solution of the 6N differential equations is simplified to determining the time of the next collision. impact occurs when solve quadratic in t to get time of collision:

Upon impact, collision is assumed elastic, so:
Systems with continuous potentials are integrated numerically. Accuracy of integration depends on the time step delta_t. Usually, delta_t is desired as large as possible, to reduce computation time, but too large delta_t leads to numerical inaccuracy in solving the set of differential equations. Thus, delta_t must be samll enough to ensure energy is conserved, typically O(10^-15 second) for polymers, since the CH vibration frequency is O(10^+14 Hz). This is one of the greatest limitations of MD!
Equilibration - there is no guarantee that an arbitrarily chosen initial state belongs to the desired ensemble with correct volume and energy. An equilibration period is required at the start of every MD run during which the system is "adjusted" as necessary (e.g. velocity rescaling to add or subtract energy) while it settles to equilibrium, so that the equilibrium corresponds to the desired ensemble. To check this equilibration, one can monitor energy or pressure as a function of time, or the velocity distribution. One function often used in the past to monitor equilibration is the Boltzmann H-function, and measure of "entropy density":
Phase Lock - This is a confinement or"quasi-ergodic" problem. It is easily possible, especially in dense systems, that a trajectory never leaves a small region of configuration space, due to the presence of high energy barriers separating it from other low energy regions of phase space:

MIT Home page | UP one level | Ask the Professor

Last modified 6/20/97 - GCR