Ammm:Md integrator

Integrator

For numerical integration techniques, the book Numerical Recipes in C is good reference.

As we already know MD simulations is solving Newton's law

${\frac {dU(x)}{dx}}=m*{\frac {d^{2}x}{dt^{2}}}$

using finite difference methods. The exact way to solve the ordinary differential equation depends on the integrator. A good integrator should allow

• large time step (However, if the timestep is too large for system, the simulation can "blow up")
• conserve energy (small energy)
• efficeint. Evalute fewer energy and force per step.

Popular integrators for MD

• Verlet leapfrog integrator

• Verlet velocity integrator

The Verlet velocity algorithm overcomes the out-of-synchrony
shortcoming of the Verlet leapfrog method 

• Runge–Kutta-4 integrator

Runge–Kutta-4 is a fourth-order Runge–Kutta method,
which is one of the oldest numerical methods for solving ordinary
differential equations. The method is self starting but requires four
energy evaluations per step.

• Beeman (http://en.wikipedia.org/wiki/Beeman%27s_algorithm)

Beeman is closely related to Verlet integration. It produces identical positions to Verlet, but is more accurate in velocities. The formula used to compute the positions at time $t+\Delta t$ is:

$x(t+\Delta t)=x(t)+v(t)\Delta t+{\frac {2}{3}}a(t)\Delta t^{2}-{\frac {1}{6}}a(t-\Delta t)\Delta t^{2}+O(\Delta t^{4})$

and this is the formula used to update the velocities:

$v(t+\Delta t)=v(t)+{\frac {1}{3}}a(t+\Delta t)\Delta t+{\frac {5}{6}}a(t)\Delta t-{\frac {1}{6}}a(t-\Delta t)\Delta t+O(\Delta t^{3})$

where

• $t$ is present time (ie: independent variable)
• $\Delta t$ is the time step size
• $x(t)$ is the position at time t
• $v(t)$ is the velocity at time t
• $a(t)$ is the acceleration at time t
• the last term is the error term, using the big O notation

which one is better?

• Beeman is my favorite for MD.
• try different time-steps and see whether the average energy change.
• for rigibd body there are different integrators
• for motorsport simulations: http://www.motorsport-sim.org/

---

Constraint

(http://en.wikipedia.org/wiki/Constraint_algorithm) To use even large time step, constraint is used to "fix" the fast vibrating bonds. In most Molecular dynamics simulation, constraints are enforced using the method of Lagrange multipliers. Given a set of n linear (holonomic) constraints at the time t,

$\sigma _{k}^{(t)}:=\|\mathbf {x} _{k\alpha }^{(t)}-\mathbf {x} _{k\beta }^{(t)}\|_{2}^{2}-d_{k}^{2}=0,\quad k=1\dots n$

where x ar ethe positions of the two atoms involved in the kth constraint (e.g. O and H in OH constraint) at the time t and "d" is the prescribed inter-particle distance.

These constraint equations, are added to the potential energy function in the equations of motion, resulting in, for each of the N particles in the system

${\frac {\partial ^{2}\mathbf {x} _{i}^{(t)}}{\partial t^{2}}}m_{i}={\frac {\partial }{\partial \mathbf {x} _{i}}}\left[V(\mathbf {x} _{i}^{(t)})+\sum _{k=1}^{n}\lambda _{k}\sigma _{k}^{(t)}\right],\quad i=1\dots N.$

Adding the constraint equations to the potential does not change it, since all $\sigma _{k}^{(t)}$ should, ideally, be zero.

Integrating both sides of the equations of motion twice in time yields the constrained particle positions at the time t + Δt

$\mathbf {x} _{i}^{(t+\Delta t)}={\hat {\mathbf {x} }}_{i}^{(t+\Delta t)}+\sum _{k=1}^{n}\lambda _{k}{\frac {\partial \sigma _{k}^{(t)}}{\partial \mathbf {x} _{i}}}\left(\Delta t\right)^{2}m_{i}^{-1},\quad i=1\dots N$

where ${\hat {\mathbf {x} }}_{i}^{(t+\Delta t)}$ is the unconstrained (or uncorrected) position of the ith particle after integrating the unconstrained equations of motion.

To satisfy the constraints $\sigma _{k}^{(t+\Delta t)}$ in the next timestep, the Lagrange multipliers must be chosen such that

$\sigma _{k}^{(t+\Delta t)}=\left\|\mathbf {x} _{k\alpha }^{(t+\Delta t)}-\mathbf {x} _{k\beta }^{(t+\Delta t)}\right\|_{2}^{2}-d_{k}^{2}=0.$

This implies solving a system of n non-linear equations

$\sigma _{j}^{(t+\Delta t)}=\left\|{\hat {\mathbf {x} }}_{j\alpha }^{(t+\Delta t)}-{\hat {\mathbf {x} }}_{j\beta }^{(t+\Delta t)}+\sum _{k=1}^{n}\lambda _{k}\left[{\frac {\partial \sigma _{k}^{(t)}}{\partial \mathbf {x} _{j\alpha }}}m_{j\alpha }^{-1}-{\frac {\partial \sigma _{k}^{(t)}}{\partial \mathbf {x} _{j\beta }}}m_{j\beta }^{-1}\right]\right\|_{2}^{2}-d_{j}^{2}=0,\quad j=1\dots n$

simultaneously for the n unknown Lagrange multipliers λ_k.

This system of n non-linear equations in n unknowns is best solved using Newton's method (oterwise matrix inversion)where the solution vector ${\underline {\lambda }}$ is updated iteratively using

${\underline {\lambda }}^{(l+1)}\leftarrow {\underline {\lambda }}^{(l)}-\mathbf {J} _{\sigma }^{-1}{\underline {\sigma }}$

where $\mathbf {J} _{\sigma }$ is the Jacobian of the equations σ_k:

$J=\left({\begin{matrix}{\frac {\partial {{\sigma }_{1}}}{\partial {{\lambda }_{1}}}}&{\frac {\partial {{\sigma }_{1}}}{\partial {{\lambda }_{2}}}}&\ldots &{\frac {\partial {{\sigma }_{1}}}{\partial {{\lambda }_{n}}}}\\{\frac {\partial {{\sigma }_{2}}}{\partial {{\lambda }_{1}}}}&{}&\ddots &{}\\\vdots &{}&{}&{}\\{\frac {\partial {{\sigma }_{n}}}{\partial {{\lambda }_{1}}}}&{}&{}&{\frac {\partial {{\sigma }_{n}}}{\partial {{\lambda }_{n}}}}\\\end{matrix}}\right)$

Since not all particles are involved in all constraints, $\mathbf {J} _{\sigma }$ is blockwise-diagonal and can be solved blockwise, i.e. molecule for molecule.

${\hat {\mathbf {x} }}_{i}^{(t+\Delta t)}\leftarrow {\hat {\mathbf {x} }}_{i}^{(t+\Delta t)}+\sum _{k=1}^{n}\lambda _{k}^{(l)}{\frac {\partial \sigma _{k}}{\partial \mathbf {x} _{i}}}$.

This is repeated until the constraint equations $\sigma _{k}^{(t+\Delta t)}$ are satisfied up to a prescribed tolerance \

or it will blow up after a number of itrations and two atoms are now far away from each other... 

Although there are a number of algorithms, SHAKE and RATTLE etc., to compute the Lagrange multipliers, they differ only in how they solve the system of equations, usually using Quasi-Newton methods.