Quaternion velocity Verlet

Intro and notation

Rigid-body dynamics can be performed either by applied distance constraints to the bonds in the body with RATTLE or by using a “quaternion integrator” [AT87]. In the quaternion approach, the rigid body is completely determined by its center-of-mass position and velocity, by its quaternion orientation and by its angular velocity.

Here, we follow the velocity Verlet form of the algorithm proposed by Rozmanov and Kusalik [RK10].

For clarity, the notation for quaternions and the rigid body’s coordinates is given explicitly.

• \(r\) is a vector denoting the position of a point or particle in Euclidean space.

• \(r_i\) is as \(r\), but for the \(i^\mathrm{th}\) particle of the rigid body.

• \(r_{i,j}\) is the \(k^\mathrm{th}\) component of \(r_i\).

• \(v_i\) is the velocity of the \(i^\mathrm{th}\) particle of the rigid body.

• \(\omega\) is the angular velocity of the rigid body.

• \(I\) is the inertia tensor, with components \(I_{i,j}\).

• \(L\) is the angular momentum, a vector.

• \(T\) is the torque on the rigid body. Elsewhere in this documentation, \(T\) is the temperature. The context is sufficient to avoid ambiguity.

• \(t\) is the time.

• \(h\) is the integrator time step.

• \(\dot \bullet\) is the time derivative of \(\bullet\).

• \(\wedge\) is the cross product.

• \(\cdot\) is the dot product.

• \(\hat x\) is the unit vector parallel to \(x\).

• \(q\), \(q_1\), etc. are quaternions.

• \(q^\ast\) is the conjugate of \(q\).

• \(|q|\) is the norm of \(q\).

• The superscript \(^B\) is for quantities defined in the frame of reference of the rigid body. Other quantities are in the laboratory reference frame.

Quaternions are defined as the sum of a scalar part \(s\) and a 3-dimensional vector part \(w\) as

\[q = s + w\]

The multiplication of the elementary vectors \(w_i=(1, 0, 0)\), \(w_j=(0, 1, 0)\) and \(w_k=(0,0,1)\) is defined as

• \(w_i w_j = w_k\)

• \(w_j w_k = w_i\)

• \(w_k w_i = w_j\)

• \(w_j w_i = -w_k\)

• \(w_k w_j = -w_i\)

• \(w_i w_k = -w_j\)

• \(w_i^2 = w_j^2 = w_k^2 = -1\)

From these definitions, one can derive the product of quaternions:

\[q_1 q_2 = \left(s_1 s_2 - w_1\cdot w_2\right) + \left(q_1 w_2 + q_2 w_1 + w_1 \wedge w_2\right)\]

where the first parenthesized term is the scalar part and the second one the vector part. The addition of quaternions is done separately for the scalar and vector parts, following conventional algebra. The conjugate of a quaternion is \(q^\ast = s - w\), its norm is \(|q|=\sqrt{s^2 + w\cdot w}\) and its inverse is \(q^{-1} = q^\ast / |q|^2\).

The utility of quaternions arises from their ability to encode any solid rotation. The rotation operator about an axis \(n\) with an angle \(\theta\) is given by \(x' = q x q^\ast\) where \(q\) is a unit quaternion with scalar component \(\cos\theta/2\) and vector component \(\sin\theta/2 \ \hat n\).

All quaternion operations in RMPCDMD are performed with the Fortran module fortran_quaternion.

Definition of rigid-body quantities

The position in the laboratory frame of a member of the rigid body is

\[r_i(t) = q(t) r_i^B q^\ast(t)\]

The equation of motion for the angular momentum in the laboratory frame is

\[\dot L = T ,\]

and in the body frame

\[\dot L^B = T^B - \omega^B \wedge L^B .\]

The equation of motion for the quaternion is

\[\dot q = \frac{1}{2} \omega q\]

Algorithm

Here, we follow section IV.A of Ref. [RK10].

Angular momentum	Torque	Quaternion	Angular velocity
\(L(h/2) = L(0) + h T(0)/2\)
\(L^B(0) = q^\ast(0) L(0) q(0)\)	\(T^B(0) = q^\ast(0) T(0) q(0)\)
			\(\omega^B(0) = L^B(0) / I^B\)
\(\dot L^B(0) = T^B(0) - \omega^B(0)\wedge L^B(0)\)
\(L^B(h/2) = L^B(0) + h \dot L^B(0)/2\)
			\(\omega^B(h/2) = L^B(h/2)/I^B\)
		\(^0\dot q(h/2) = q(0) \omega^B(h/2)/2\)
		\(^0q(h/2) = q(0) + h ~^0\dot q(h/2)/2\)
Start of the iterative procedure
\(^{k+1}L^B(h/2) = ~^kq^\ast(h/2) L(h/2) ~^kq(h/2)\)
			\(^{k+1}\omega^B(h/2) = ~^{k+1}L^B(h/2) / I^B\)
		\(^{k+1}\dot q(h/2) = ~^k q(h/2) ~^{k+1}\omega^B(h/2)/2\)
		\(^{k+1}q(h/2) = q(0) + h~ ^{k+1}\dot q(h/2)/2\)
End of the iterative procedure
		\(q(h) = q(0) + h \dot q(h/2)\)
Update positions using \(q(h)\)
Compute the updated forces and torques in the lab frame
	\(T^B(h) = q^\ast(h)T(h)q(h)\)
\(L(h) = L(h/2) + h T(h)/2\)
Update the velocities using \(\omega^B(h) = L^B(h)/I^B\).

The steps until the update of \(q(h)\) and of the positions is the first step of the velocity Verlet algorithm. These steps are implemented in rigid_body_vv1.

The update of \(L(h)\) and of the velocities form the second part of the velocity Verlet algorithm and are implemented in rigid_body_vv2.