   ## 4.3.1 Constraint dynamics with one bond length constraint

Molecular Dynamics for Linear Diatomics: Constraint method

The classical method to simulate diatomics is due to Singer (1977) and Fincham (1984), as also described in Allen-Tildesley, Computer Simulation of Liquids (Clarendon, Oxford 1987).

However, we prefer to use the method of Constraint Dynamics'' - which is particularly simple when applied to one single constraint.

Assume a diatomic with potential forces acting on the two atoms, respectively, and write etc. The only constraint is . The equations of motion are then   (4.6)   (4.7)

True to the idea of constraint dynamics, proceed as follows:
1. Given the positions at time , integrate the equations of motion for one time step only considering the potential forces. The resulting positions are denoted . These preliminary position vectors will not fulfill the constraint equation; rather, the value of will be non-zero.
2. Making the correction ansatz   (4.8)   (4.9)

with undetermined , requiring that the corrected positions fulfill the constraint equations we find (4.10)

with .
This equation for may be solved exactly: (4.11)

where the negative sign prevails.

Simple example: Free rotation of a 2D dumbbell

For this problem the constraint method should be exact, since neither the Verlet propagator nor the constraint reconstruction introduce errors.

At let , ; the masses are , the angular velocity is , and a time step of is assumed. Then the exact positions at time are (4.12)

and (4.13)

where we have written and .

Now let us test the constraint dynamics procedure against this exact result. For the Verlet integrator we need , for which we take the exact values and . After the free flight the new positions are (4.14)

with and . Using and we find for  (4.15)

Since we have . Insert this in the correction formula 4.8-4.9 to find (4.16)

which is identical to the exact values.

Thermostats / Nosé-Hoover method:
In a nonequilibrium process work must be invested which systematically heats up the system. In lab experiments a thermostat takes care of this; in simulations the problem of thermostating is non-trivial. The temperature of a molecular dynamics sample is not an input parameter to be manipulated at will; rather, it is a quantity to be measured'' in terms of an average of the kinetic energy of the particles, ( dimension). Suggestions as to how one can maintain a desired temperature in a dynamical simulation range from repeatedly rescaling all velocities (brute force thermostat'') to adding a suitable stochastic force acting on the molecules. Such recipes are unphysical and introduce an artificial trait of irreversibility and/or indeterminacy into the microscopic dynamics.

A interesting deterministic method of thermostating a simulation sample was introduced by Shuichi Nosé and William Hoover. Under very mild conditions the following augmented equations of motion will lead to a correct sampling of the canonical phase space at a given temperature :      The coupling parameter describes the inertia of the thermostat. The generalized viscosity is temporally varying and may assume negative values as well.

For systems of many degrees of freedom (particles) this is sufficient to produce a pseudo-canonical sequence of states. For low-dimensional systems it may be necessary to use two NH thermostats in tandem [MARTYNA 92].   F. J. Vesely / University of Vienna