# How to calculate quaternions from principal axes of an ellipsoid?

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,

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.

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.

• Hi @LonelyProf, I did try both yours and orion's methods to calculate the quaternions (and check it several times) but still haven't got a correct match. It is correct that I am doing MD simulations using Lammps. When I use a visualisation tool to displace the atomistic model and the ellipsoid model they don't agree with each other. I calculated the principal axes by taking the eigen-vectors of the moment of inertia tensor, and I'm pretty sure I've got them right because they agree with the molecule coordinates when I plotted them out on a 3D graph... I'm very confused... – Nguyen Thi Lan Huong Aug 22 '18 at 9:17
• I'm worried that the comments section may develop into a long question-and-answer thread which is not how maths stackexchange is intended to work. Also, I suspect that the real problem is no longer the one you had at the beginning: I believe that you have two satisfactory answers to that question. It looks more like your problem is computational: either "how to plot an ellipsoid" or "how to understand how LAMMPS works". I am not even sure that such questions would be welcome on scicomp.stackexchange.com because they will not help much with software packages. – user575517 Aug 22 '18 at 9:38
• Since you have insufficient reputation for chat, let me risk a question or two here. You plotted the inertia axes, along with the original atomistic coordinates, and they look right? But when you try and plot an ellipsoid with those same principal axes, it looks wrong? How are you trying to plot the ellipsoid? Does it rely on knowing the quaternions? – user575517 Aug 22 '18 at 9:46
• I plotted the atom coordinates and the principal axes using some python code: input coordinates, centre of mass (COM), axis vectors, draw a surface based on atom coordinates, draw lines from COM along the axis vectors, and the axes align right with the surface. For the ellipsoid model, I wrote a Lammps data file with calculated COM and quaternions, ran Lammps for 0 step to output a dump file and visualise it + the atomistic model using OVITO. Sorry about leading this to a long conversation about something else... I also think the problem is how they define ellipsoid orientation in Lammps. – Nguyen Thi Lan Huong Aug 22 '18 at 10:10
• I suspect that I can't help much more. If you have tried the transpose of the rotation matrix (or equivalently switching $q_1,q_2,q_3 \rightarrow -q_1,-q_2,-q_3$ leaving $q_0$ unchanged), then there's not much else that can be wrong. You can also try plotting the ellipsoid surface by hand using its definition here en.wikipedia.org/wiki/Ellipsoid in the section "As quadric" (defining a matrix $\mathbf{A}$ from your axes: no need for quaternions). – user575517 Aug 22 '18 at 10:29