Ammm:Rigid body MD (Steven 12/2)

Introduction

Recall our primary motivations for coarse-graining:

1. Eliminate high frequency motions (if possible), AND,

2. Eliminate "uninteresting" degrees of freedom (both in translation and in rotation)

Rigid Body Equations of Motion

• Linear Momentum:
  

1. Definition

2. Rate of Change

 

which, because (for our purposes) mass is constant:

 

Translational motion can be considered as usual.

• Angular Momentum:
  

1a. Particle

      

1b. Distributed Masses at sites 'i'

$L=\sum _{i}^{}\mathbf {r} _{i}\times m_{i}\mathbf {V} _{i}=I\omega$

2. Rate of Change: Torque

$\tau ={\frac {\mathrm {d} \mathbf {L} }{\mathrm {d} t}}={\frac {\mathrm {d} I}{\mathrm {d} t}}\omega +I{\frac {\mathrm {d} \omega }{\mathrm {d} t}}$    

Thus, in its purest form the torque relates both a change in I AND a change in omega.

 

3. Moment of Inertia, I

            ;           

These Moments of Inertia (MOI) essentially describe the distribution of the mass of our rigid body (is the mass tightly packed or diffuse?). These MOI are functions of the position of the particles only.

                                                        

The closer the mass is to the rotation axis the smaller the I

The Moments of Inertia can be expressed together with their eigenvectors to define the 'principal axes' of the body:

It is extremely convenient to express rotations about these axes.

Frames of Reference

• Unfortunately these equations are only this simple in the Local Frame of reference (ie., a coordinate axis attached to the rigid body in question).
  

In the local frame our relation for torque:

$\tau ={\frac {\mathrm {d} \mathbf {L} }{\mathrm {d} t}}={\frac {\mathrm {d} I}{\mathrm {d} t}}\omega +I{\frac {\mathrm {d} \omega }{\mathrm {d} t}}$

reduces to

$\tau ={\frac {\mathrm {d} \mathbf {L} }{\mathrm {d} t}}=I{\frac {\mathrm {d} \omega }{\mathrm {d} t}}=I\alpha$

because the MOI is time-invariant in the local frame.  We know this because in the local frame we have a rigid body, and by definition the positions of particles do not change in something that is rigid.  Thus in this local frame the torque directly relates to the acceleration.

• However, when we have a rigid body rotating in a global (or laboratory) frame of reference we need a new relation for the change of angular momentum in these global coordinates:
  

The Euler Equations

Upon substituting   and simplifying gives the Euler Equations for rigid body motion:

Take note that the accelerations are functions of the velocities.  Recall the velocity verlet integrator for translation:

Where,

$a=-{\frac {1}{m}}{\frac {\partial U(x_{i})}{\partial x_{j}}}$

The dependence of the acceleration on the positions of the system alone permits the use of the Velocity Verlet integrator with translation because we can use the "new" postiions to find the new velocity.  In the case of the Euler Equations we're not so lucky -- the accelerations are a function of BOTH the positions (through the moments of inertia I) and the velocities (omega). As we'll see, this dual dependence of the acceleration will require iterative procedures.

Rotations

Frame Rotations v. Temporal Rotations

We have two distinct types of rotations to consider:

• Frame rotations
• Temporal rotations

Rotation Matrices

The rotation of a vector about some angle theta in two dimensions is:

                       

Pretty simple, eh?  Unfortunately we exist (and our simulations typically exist) in three dimensions...

Euler Angles

Euler angles are a traditional means for describing the rotation of a vector in three dimensions and follow directly from the formalism of rotation matrices.  For the three dimensional rotation of a vector 'p' we have:

In this 'zxZ' convention (there are 12 different permutations of rotations, but this is the most common) we have the following rotations:

-alpha: rotation about the original 'z' axis,

-beta: rotation about the new 'x' axis (which you may see referred to as the 'Line of Nodes'), and

-gamma: rotation about the new 'Z' axis

                                   File:Euler2.gif

[Left Image: z: Blue; x: Green; Z: Red]

Unfortunately at certain rotations we can run into a singular condition referred to as "gimbal lock" en.wikipedia.org/wiki/Gimbal_lock . This almost crashed Apollo 11 into the Moon.  Let's avoid this...

Quaternions

Quaternions are a special representation of rotations in which we lose the explicit representation of angles and instead focus on a four component "complex" vector:

$q=(w,x,y,z)=w+x\mathbf {i} +y\mathbf {j} +z\mathbf {k}$

'w' is a scalar and (x,y,z) form a "complex" vector [but unfortunately with strange multiplication rules en.wikipedia.org/wiki/Quaternions#Multiplication_of_basis_elements ]

A simple interpretation:

-the scalar 'w' is the magnitude of the rotation, and

-the complex vector (x,y,z) is the axis of that rotation.

This interpretation isn't too far off the truth:

To rotate a vector 'v' with a quaternion we use the following relation:

As an example, consider: en.wikipedia.org/wiki/Quaternions_and_spatial_rotation#Example

Integration Schemes

• As we'll see below the fact that we have two physically relevant frames of reference will force us into more complicated integrations than are necessary for translation (in which we don't have a "local" frame, meaning that our acceleration is not a function of position)
• The first two schemes presented use a "predictor - corrector" version of the familiar Velocity Verlet / Leapfrog integration alogrithm
• The first two schemes handle the frames of reference in different ways.

Quaternion: DL_Poly

DL_Poly is a British MD code: www.ccp5.ac.uk/

• The general scheme is as follows:

1) Translational motion (global reference frame)

2) Calculate torques in global reference frame

3) Transform torques to local reference frame via quaternions

4) Calculate angular velocities / positions in local reference frame

5) Update coordinates: find new quaternions.

As we see below, the ease of calculation in the local frame is paid for by extra computation in step 5. 

"Global": Tinker's RGDSTEP

The Rigid Body Dynamics integrator in Tinker follows a different path, notably never rotating into a local reference frame.

Recall our general relation for torque as the time derivative of angular momentum (tau and L above, T and J in the reference below):

$\tau ={\frac {\mathrm {d} \mathbf {L} }{\mathrm {d} t}}={\frac {\mathrm {d} I}{\mathrm {d} t}}\omega +I{\frac {\mathrm {d} \omega }{\mathrm {d} t}}$

The DL_Poly scheme above capitalizes on the fact that in the local reference frame ${\frac {\mathrm {d} I}{\mathrm {d} t}}=0$ .  This advantage comes at the price of having to (1) rotate into the local frame and (2) update the quaternions specifying that rotation through an iterative scheme.

RGDSTEP attempts to avoid these complications by remaining in the global (or "Laboratory") reference frame.

• General scheme:

1) Translational motion (global reference frame)

2) Calculate torques (global reference frame)

3) Directly calculate the new angular momentum: $\mathbf {L} (t+\Delta t)=\mathbf {L} (t)+\Delta t\mathbf {T}$

3b) Recognizing that in the global reference frame the moments of intertia change, we have: $\mathbf {L} (t+\Delta t)=\mathbf {I} (t+\Delta t)\mathbf {\omega } (t+\Delta t)$

4) Iterative scheme to simultaneously converge I and omega:

4a) Guess $\mathbf {\omega } (t+\Delta t)$ to estimate $\mathbf {I} (t+\Delta t)$

4b) With new $\mathbf {I} (t+\Delta t)$ estimate and KNOWN $\mathbf {L} (t+\Delta t)$ calculate new guess for $\mathbf {\omega } (t+\Delta t)$

4c) Repeat (4a) and (4b) until the $\mathbf {\omega } (t+\Delta t)$ estimate converges ("typically 3 to 4 steps")

Thus while we have the advantage of remaining in the global reference frame we're still forced into an iterative scheme due to the fundamental, inescapable fact that we have two distinct physical frames of reference (the global / laboratory frame and the principal frame of the body).

NOTE: in the below α,β are the x,y,z indices

Constraints: QSHAKE

• Before we move onto a more sophisticated method, let's discuss constraints
  
• Recall that (one of) our original motivation(s) for coarse-graining was to eliminate "fast" motions to enable a larger timestep.  Towards this goal in some applications it's desirable to consider the connections between rigid bodies to be "rigid" themselves.  In this case we won't automatically know the forces (and resulting torques) due to these rigid constraints because they're not part of our potential (from which we typically get our forces and torques).
• Therefore we need a scheme to calculate what are essentially "reaction" forces: equal and opposite in magnitude acting along a specific constraint direction.

SHAKE:

• The original SHAKE alogrithm has been popularly implemented to increase simulation timesteps by "freezing" the fast motions, typically the C-H bond stretching. In this context we constrain the length of the bond, avoiding instability from 'H' motions with our constraints.
• The general scheme:

1) Normal, unconstrained integration

2) Iterative calculation of the additional constraint force to be applied, using the temporary (unconstrained) position of the two atoms as an estimate of this force and waiting for convergence to the desired constraint length.

QSHAKE:

• QSHAKE is the extension of SHAKE to rigid bodies in which we recognize that the constraint forces introduce torques as well as translations on the constrained bodies.
• In essentially the same scheme as SHAKE we have:

Order-N: MBO(N)D / Lobatto III a-b Integrator

• In addition to the straight forward, classically-inspired MD simulation methods represented by those above there is a class of so-called "Order-N" "Efficient" multibody approaches.  The 'Order-N' reference does not refer to the error of the simulation (such as with a truncated series), but rather to the number of calculations required at each step.
• A hallmark of these approaches are 'internal' or 'generalized' coordinates, as opposed to Cartesian coordinates. Whereas a typical simulation might require 3N Cartesian coordinates to define the system an internal coordinate representation in which only torsion angles are considered (as in structure refinement with Xplor cns-online.org/v1.21/ ) could drop that requirement to N/3.

• Instead of worrying about referring the local frame of reference of each body to the global frame (ie., body(i) -> global) we reference each body to its neighbor (ie., body(i+1) -> body(i)).  

• Another hallmark of these schemes is the so-called Lobatto III a-b partitioned Runge-Kutta integration scheme.  Sounds fancy, but in reality it's pretty simple:

Unfortunately we're still stuck with the fact that our acceleration is a function of the velocity. In fact, since we reference each body to another, already-rotating body we're forced to consider the Coriolis acceleration ... which sounds fun because it makes toilets spin the other way in the Southern Hemisphere (supposedly), but the velocity dependence of the acceleration creates the same requirement for iterative / recursive approaches that we had when we stuck with the Euler equations .   en.wikipedia.org/wiki/Coriolis_effect

• The "efficiency" of these algorithms comes from the way in which the equations of motion of each body is solved:

Since each body is referenced to it's "base" neighbor, the schemes perform three "sweeps" :

Base -> Tip: propagates the velocities to compute spatial velocities

Tip -> Base: uses velocities to compute body properties (such as moment of inertia)

Base -> Tip: calculation of accelerations with updated body properties

• The details of these algorithms are heavy, but I refer you to the following references:

Anderson (2003) "Improved 'Order-N' Performance Algorithm for the Simulation of Constrained Multi-Rigid-Body Dynamic Systems"

Mukherjee (2008) "Substructured molecular dynamics using multibody dynamics algorithms"

Schwieters (2001) "Internal Coordinates for Molecular Dynamics and Minimization in Structure Determination and Refinement"

Chun (2000) "MBO(N)D: A Multibody Method for Long-Time Molecular Dynamics Simulations"

Multiscale OCTA http://octa.jp/OCTA/whatsOCTA.html
DL_CG DL_POLY http://www.ccp5.ac.uk/
Moldy http://www.ccp5.ac.uk/moldy/moldy.html
Mdynamix http://www.fos.su.se/~sasha/mdynamix/

Other useful info:

Xplor http://www.ocms.ox.ac.uk/mirrored/xplor/manual/htmlman/node221.html
http://kwon3d.com/theory/moi/triten.html