1549380323-Statistical Mechanics Theory and Molecular Simulation

Symplectic quaternions 119

3.12 Symplectic integration for quaternions

In Section 1.11, we showed that the rigid body equations of motion could be expressed
in terms of quaternions (cf. eqns. (1.11.44) to (1.11.47)). Eqn. (1.11.40) relates the
four-component quaternion vectorq= (q 1 ,q 2 ,q 3 ,q 4 ) to the Euler angles. We also saw
that when the rigid body equations of motion are expressed in termsof quaternions,
we must impose an additional constraint that the vectorqbe a unit vector:

∑^4

i=1

q^2 i= 1. (3.12.1)

As was noted in Section 3.10, the iterations associated with the imposition of holo-
nomic constraints affect the symplectic and time-reversibility properties of otherwise
symplectic solvers. However, the simplicity of the constraint in eqn.(3.12.1) allows the
quaternion equations of motion to be reformulated in such a way that the constraint is
satisfied automatically by the dynamics, thereby allowing symplectic integrators to be
developed using the Liouville operator. In this section, we will present such a scheme
following the formulation of Milleret al.(2002).
Recall that the angular vectorωwas defined byω= (0,ωx,ωy,ωz), which has
one trivial component that is always defined to be zero. Thus, there seems to be an
unnatural asymmetry between the quaternionqand the angularω. The idea of Miller
et al.is to restore the symmetry betweenqandωby introducing a fourth angular
velocity componentω 1 and redefining the angular velocity according to

ω(4)= 2ST(q)q ̇≡(ω 1 ,ωx,ωy,ωz). (3.12.2)

The idea of extending phase spaces is a common trick in molecular dynamics, and
we will examine numerous examples throughout the book. This new angular velocity
component can be incorporated into the Lagrangian for one rigid body according to

L=

1

2

[

I 11 ω^21 +Ixxωx^2 +Iyyω^2 y+Izzωz^2

]

−U(q). (3.12.3)

HereI 11 is an arbitrary moment of inertia associated with the new angular velocity
component. We will see that the choice ofI 11 has no influence on the dynamics.
From eqns. (3.12.2) and (3.12.3), the momentum conjugate to the quaternion and
corresponding Hamiltonian can be worked out to yield

p=

2

|q|^4

S(q)D−^1 ω(4)

H=

1

8

pTS(q)D−^1 ST(q)p+U(q), (3.12.4)

where the matrixDis given by

D=







I 11 0 0 0

0 Ixx 0 0

0 0 Iyy 0

0 0 0 Izz





. (3.12.5)