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,
Best Answer
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.