Hail Euler and Farewell: Rotational Motion in the Laboratory Frame.

W. Smith, Computational Science and Engineering Department, Daresbury Laboratory, Warrington WA4 4AD, United Kingdom (E-mail: w.smith@dl.ac.uk)

Abstract

In this paper we consider the treatment of rotational motion in molecular dynamics simulation and re-examine the Euler approach. It is shown that, for numerical integration, it is not necessary to couch the problem in terms of the principal reference frame in which the moment of inertia tensor is diagonal. From these considerations a new algorithm is developed which is solved entirely in the laboratory reference frame and avoids the pitfalls associated with applying Euler’s method in numerical simulations.

Introduction

In molecular dynamics the numerical integration schemes for integrating the rotational motion of a rigid body (in our case a rigid molecule), are conventionally discussed in terms of Euler’s rotational equations of motion [1]. In these equations it is assumed that, at a given instant, the rotation is described in a frame of reference in which the moment of inertia tensor is diagonal. This instantaneous frame of reference is known as the principal frame. In the principal frame the relationship between the applied torque and the rate of change of angular velocity is particularly simple and allows a simple numerical integration of the angular velocity. However, because the body rotates, it is necessary for the principal frame to rotate also, in order to ensure that the moment of inertia tensor remains diagonal.

Thus a numerical solution of Euler’s equations needs to be supplemented by a prescription for updating the principal frame as the simulation proceeds. While such a prescription can be easily obtained [2], it is well known that there are disadvantages to the scheme, due to the real possibility of the equations becoming indeterminate for certain angles of orientation. This was discussed in detail by Evans and Murad [3] some years ago and led to the introduction of quaternions, with singularity free equations of motion, into the methodology of molecular dynamics. There is however another approach, and that is to relinquish the concept of the principal frame of reference altogether and tackle the problem directly in the laboratory frame. It turns out that this approach has advantages for a numerical integration scheme. Firstly however we must re-examine the mathematical treatment of rotational motion as it is conventionally presented.

The Equations of Motion for Rigid Molecules

A molecule in the context of this treatment is a rigid assembly of point masses (atoms) of mass m_i, where i=1,…,n. From this simple definition a number of important molecular properties are obtained:

(1)

In which M is the molecular mass, R the centre-of-mass, d_i the location of an atom with respect to the centre-of-mass and I is the moment of inertia tensor.

The velocity and angular momentum associated with these variables are:

(2)

Where V is the molecular (translational) velocity, J is the angular momentum and ω is the angular velocity, while v_i is the velocity of an individual atom in the laboratory frame. (Note the use of the Einstein convention regarding tensor indices in equation (2c) here, and elsewhere in this paper.)

The dynamics of each molecule are driven by the molecular force F and torque T:

(3)

In which the atomic forces f_i are the vector sums of atom-atom forces acting between molecules.

It is important to note that in these definitions, moment of inertia, angular momentum, torque and angular velocity are defined with respect to the molecular centre-of-mass. It should also be noted that all have been defined with respect to the laboratory frame of reference, which means any inertial frame that happens to be convenient.

The equations of motion for each molecule (again in the laboratory frame) are well known:

(4)

The first of these is the standard Newtonian equation for translational motion, while the second embodies the rotational motion of the molecule. Following Euler this equation may be cast into the principal frame where it takes the familiar form know as Euler’s equations [1] but here we take an alternative approach. Equation (4b) can be expanded into:

(5)

Where the time derivative of the moment of inertia tensor can be obtained from equation (1d) as:

(6)

We now recognise that

(7)

in which the Levi-Civita tensor ε is 1 for a cyclic permutation of indices, -1 for a non-cyclic permutation and zero otherwise. Substitution of (7) into (6) provides the following result:

(8)

If equation (8) is substituted into (5) we obtain:

(9)

It is evident that (9) forms the basis for an integration of the angular velocity. We may therefore outline a possible numerical scheme for integrating the equations of motion for rigid molecules.

An integration scheme based on the velocity Verlet algorithm [2] is presented in Scheme 1. The scheme is remarkably simple, but experienced molecular dynamicists will notice a serious weakness at line g). This line is meant to compute the change in the orientation of the vector d_i as a result of the rotation generated by vector ω^n+1/2. However in this form it does not conserve the normalisation of d_i and consequently cannot maintain the structural integrity of the molecule during the simulation. Fortunately the geometric interpretation of this operation affords a reliable alternative prescription.

In Figure 1 (a) we provide a geometric interpretation of the action of an angular velocity vector ω on a typical displacement vector d represented the cross product ω x d. Under this operation the tip of vector d traces out the circle shown, the plane of the circle being perpendicular to the vector ω. The radius of the circle is represented by the vector u which is given by:

(10)

Figure 1 (b) shows the effect of the rotation on the vector u after a time interval of δt, which is a rotation of u to u’ though the angle ωδt. The change in the vector d is therefore correctly defined by the addition of the vector Δu=(u’-u) and not vector (ω x d)δt. By simple vector algebra we have:

(11)

Thus in place of step g) in Scheme 1 we have

These equations are guaranteed to conserve the norm of vector d_i.

Application

The proposed algorithm has been implemented in Fortran [4] and run on a set of model systems (all of which are defined in Lennard-Jones reduced units, and all pair interactions with ε=1 and σ=1 shifted to zero force at the cut off):

1. Case 1: 27 tetrahedral clusters of Lennard-Jones sites, with intramolecular site-site distance 0.7585, site masses {0.5,1.5,0.75,1.25}, cluster density ρ*=0.175, temperature T*=1.5 and time step Δt=0.001.

2. Case 2: pseudo-ammonia; 27 molecules, with H atoms at 3 vertices of a tetrahedron and the N atom at the tetrahedron centre. Site masses {1.0 (N),0.1 (H)}, N-H distance 0.7255, molecular density ρ*=0.2, temperature T*=1.0 and time step Δt=0.001.

3. Case 3: pseudo-water; 27 molecules, with H atoms at 2 vertices of a tetrahedron and the O atom at the tetrahedron centre, plus an additional L-J site bisecting O-H bonds at distance 0.1644 from oxygen, O-H bond length .7119, site masses {1.0 (O), 0.1 (H), 0.0 (additional site)}, molecular density ρ*=0.225, temperature T*=1.25 and time step Δt=0.001.

4. Case 4: As Case 3, with time step Δt=0.0005.

5. Case 5: As Case 1, with time step Δt=0.0005.

In each case a simulation of 1 million time steps was conducted, with the first 10,000 time steps reserved for equilibration. The average energy and drift over 990,000 time steps was obtained.

The results are presented in Table 1.

Table 1

The total drift presented in Table 1 was obtained from the slope of the Energy versus Time plot fitted by least squares and multiplied by 990,000Δt.

Cases 1, 2 and 3 represent molecular structures with decreasing moments of inertia, so it is anticipated that the algorithm will show greater drift progressing from 1 to 3. This is clearly seen in Table 1, The drift in Cases 1 and 2 is arguably negligible, but the result for Case 3 would not normally be acceptable. The particularly poor result obtained reflects, no doubt, the very low moment of inertia of this molecule. Nevertheless, reducing the time step by half, improves the result by at least two decades. The equivalent reduction in time step applied to Case 1 (i.e. Case 5) reduces the drift by a factor of ~4 (with a switch in sign also). For many applications the results obtained here would be quite acceptable.

Conclusion

The main conclusion to be drawn from this note is that casting the rotational equations of motion into a principal frame of reference is unnecessary if a numerical integration scheme is required. Solving the motion in the laboratory frame does not require recourse to quaternions as there are no difficulties with particular molecular orientations; the equations are singularity free. Finally, to this author at least, the method has an appealing directness that leans more on physical than mathematical intuition, which will be an advantage for others also.

References

[1] H. Goldstein, Classical Mechanics, Addison Wesley (1980).

[2] M.P. Allen and D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press (1980)

[3] D.J. Evans amd S. Murad, Molec. Phys. 34 (1977) 327.

[4] A copy of the test program revo.f is available from the author on request.