How to calculate quaternions from principal axes of an ellipsoid?

Question

I am working a project that involves calculating quaternions of ellipsoids providing that I already know the unit vectors of their principal axes. The reason for this is the program I am using to work with these ellipsoids require quaternions. The description from the manual of that program of quaternions is: "The particles must define a quaternion for their orientation. Note that particles defined as ellipsoid have 3 shape parameters. The 3 values must be non-zero for each particle set by this command. They are used to specify the aspect ratios of an ellipsoidal particle, which is oriented by default with its x-axis along the simulation box’s x-axis, and similarly for y and z. If this body is rotated (via the right-hand rule) by an angle theta around a unit rotation vector (a,b,c), then the quaternion that represents its new orientation is given by $(cos(\theta/2), a*sin(\theta/2), b*sin(\theta/2), c*sin(\theta/2))$"

From some source in the literature (which I am not sure if is correct), I calculate quaternions like this:

• Construct a orientation matrix whose columns are the principal axes
• Solve for e-values and e-vectors of that matrix. The e-values should be (1, $cos(\theta) -i sin(\theta)$,$cos(\theta) +i sin(\theta)$) where theta is the rotational angle. The unit e-vector $e = [e_1 e_2 e_3]^T$ that correspond to e-value 1 is an invariant principal axis of rotation.
• The quaternions are: $q_w = cos(\theta/2)$, $q_i = e_1sin(\theta/2)$, $q_j = e_2sin(\theta/2)$, $q_k = e_3sin(\theta/2)$

However, when I try to check if the ellipsoids using the quaternions I calculated agree with the principal axes (using a graphic tool), they don't agree with each other.

I am not an expert in Maths, and only do this because of my project. I really appreciate if someone can help me to point out what I did wrong, or give me some good references, or even just explain to me what the manual description means.

Cheers,

Answer 1

Knowing the orientation vectors is essentially the same as knowing the transformation matrix, which is constructed by putting these vectors (they have to be normalized and orthogonal) into matrix columns.

First possible issue: arrange the vectors so that they form a positively oriented triad (the determinant of the matrix should be 1, not -1). Otherwise, your matrix is not a rotation but a rotation+reflection. Fixing it is easy, as it's an ellipsoid and you can flip the sign of any main axis vector without changing anything. Second possible issue is that the order of axes is different and you are putting the long axis where the short axis should be (it depends what is the "home" position of the ellipsoid before rotation).

Then, the easiest way is to use the invariant properties of the matrix to construct the axis and angle of rotation.

The trace of the matrix (sum of diagonal) always equals $2\cos\theta+1$, so this is the cheapest way of getting the rotation angle. Note that you still don't know if it's clockwise or anticlockwise (cosine destroys the sign), but you didn't choose the axis orientation either, so it's fine.

The axis is easiest to get from the antisymmetric part of the matrix. If your matrix is $R$, then calculate

$$ R-R^T= \begin{bmatrix} 0&c&-b\\ -c & 0 &a \\ b & -a & 0 \end{bmatrix} $$ So, taking the difference of the matrix to its transpose, you can read the axis $(a,b,c)$ from the off-diagonal elements. Essentially, you have $$a=R_{23}-R_{32},\quad b=R_{31}-R_{13},\quad c=R_{12}-R_{21}$$

The length of this vector is not unit, but $2\sin\theta$, so you must normalize it (you can also verify this angle if it matches the previous one calculated from the trace).

Now that you have the normalized $(a,b,c)$ axis (unit vector) and the angle $\theta$ you can construct the quaternion.

There's always a possibility that there is another minus that we lost when defining the angle and axis direction - try both and see if it helps.

Answer 2

The phrases you quote seem to be from the LAMMPS manual, so it is possible that you will be repeating this calculation for a large number of ellipsoids in a molecular dynamics simulation. If that's so, I feel that it may be helpful to provide an extended version of the excellent answer of @orion. This is more numerically robust against the danger of dividing by a number close to zero. It goes back to S. W. Shepperd, "Quaternion from rotation matrix", J. Guidance and Control, 1, 223-224 (1978). The paper itself is tricky to find online, but if you search for "shepperd quaternion" you can turn up several places where the method is described. But anyway, here it is.

There is no need to introduce any trigonometric functions (unless you really need to know $\theta$, which I doubt). Instead of your $(q_w,q_x,q_y,q_z)$ I am going to use the notation $(q_0,q_1,q_2,q_3)$. The rotation matrix takes the form (which you can find in books such as Classical Mechanics by H Goldstein) $$ \begin{pmatrix} R_{11} & R_{12} & R_{13} \\ R_{21} & R_{22} & R_{23} \\ R_{31} & R_{32} & R_{33} \end{pmatrix} = \begin{pmatrix} q_0^2+q_1^2-q_2^2-q_3^2 & 2(q_1q_2+q_0q_3) & 2(q_1q_3-q_0q_2) \\ 2(q_1q_2-q_0q_3) & q_0^2-q_1^2+q_2^2-q_3^2 & 2(q_2q_3+q_0q_1) \\ 2(q_1q_3+q_0q_2) & 2(q_2q_3-q_0q_1) & q_0^2-q_1^2-q_2^2+q_3^2 \end{pmatrix} $$ (It is possible that LAMMPS will have used the transpose of this matrix: if so, the following algebra is easily modified).

Start by calculating the rotation matrix from the unit vectors of the ellipsoid, exactly as advised in the answer of @orion.

Let the trace of the matrix be $T = R_{11} + R_{22} + R_{33}$. Remember that the quaternions are normalized: $q_0^2+q_1^2+q_2^2+q_3^2 = 1$. Then, looking at the diagonal elements of $\mathbf{R}$, \begin{align*} 4q_0^2 & = 1 + 2T - T \qquad = 1+T \\ 4q_1^2 & = 1 + 2R_{11} - T \\ 4q_2^2 & = 1 + 2R_{22} - T \\ 4q_3^2 & = 1 + 2R_{33} - T \end{align*} The first equation may look a bit odd, but it is written this way to make the following choice clear: if you choose the largest (most positive) of $\{ T,R_{11},R_{22},R_{33} \}$, then this clearly corresponds to the largest of $\{ q_0^2, q_1^2, q_2^2, q_3^2 \}$. Make that choice: calculate that component $q_i$ from the corresponding equation. You will need to take a square root: you can define this component $q_i$ to be positive without loss of generality.

Now look at the off-diagonal elements of $\mathbf{R}$: \begin{align*} 4q_0q_1 &= R_{23} - R_{32} & 4q_2q_3 &= R_{23} + R_{32} \\ 4q_0q_2 &= R_{31} - R_{13} & 4q_3q_1 &= R_{31} + R_{13} \\ 4q_0q_3 &= R_{12} - R_{21} & 4q_1q_2 &= R_{12} + R_{21} . \end{align*} Exactly three of these six equations will involve our chosen component $q_i$ (the one whose magnitude we already know is largest). Use these three equations to calculate the other components: this will involve dividing by $4q_i$.

If it turns out that the trace $T$ is larger than any of the diagonal elements of $\mathbf{R}$, then $q_0$ will be the largest (in magnitude) of the quaternion parameters, and can be calculated as $\sqrt{(1+T)/4}$; the components $q_1$, $q_2$, $q_3$ follow (from the three equations on the left of the above set of six) and the procedure reduces to the answer of @orion (once the functions involving $\theta$ have been eliminated). However, Shepperd's algorithm guarantees that we avoid numerical troubles if $q_0$ turns out to be small.