I work in computational chemistry, and I want make a simple Python program that rotates a molecular structure by some angle $\alpha$ around some axis defined by the unit vector $\vec{v}$.
I define the axis of ration by computing the center of mass (COM) of the molecule, and then translating the coordinates such that the COM is at the origo. In doing so, the axis of rotation is the unit vector, and it passes through the COM. It is a bit arbitrary in which direction the axis of rotation points, as it will depend on the orientation of the initial molecular structure. But for my purposes, this does not matter.
Something is clearly incorrect with my derivation, as the atoms in the "rotated" structure actually all lie in the same plane, when I wanted to rotate by 90 degrees.
So, given the quaternion $q$
$$ \begin{align} q &= \cos\frac{\alpha}{2} + ix\sin\frac{\alpha}{2} + jy\sin\frac{\alpha}{2} + kz\sin\frac{\alpha}{2} \end{align} $$
where $x$, $y$, and $z$ are the elements of the vector defining the axis of rotation, and the set of atomic position vectors $\vec{p}_i$ defined by the origo and the initial atomic coordinates, we have to solve the two Hamilton products
$$ qpq^* = \mathcal{H(\mathcal{H}(q, p), q^*}) $$
The vectors $\vec{p}_i$ have elements $(p_x, p_y, p_z)$.
I find that the first product becomes
$$ \begin{align} \mathcal{H}(q, p) &= -sin\frac{\alpha}{2} (xp_x + yp_y + zp_z) \\ &+ i (p_x\cos\frac{\alpha}{2} + yp_z\sin\frac{\alpha}{2} - z_py\sin\frac{\alpha}{2}) \\ &+ j(p_y\cos\frac{\alpha}{2} - xp_z\sin\frac{\alpha}{2} + zp_x\sin\frac{\alpha}{2}) \\ &+ k(p_z\cos\frac{\alpha}{2} + xp_y\sin\frac{\alpha}{2} - yp_x\sin\frac{\alpha}{2}) \end{align} $$
Upon evaluating the full Hamilton product, I find that the real parts cancel (which to me indicates that I am on the right track). Note that I from here on have introduced the shorthands $\gamma = \cos\frac{\alpha}{2}$ and $\theta = \sin\frac{\alpha}{2}$.
I find that the quaternion part becomes
$$ \begin{align} qpq^* &= x^2p_x\theta^2i + xyp_y\theta^2i + xzp_z\theta^2i + xyp_x\theta^2j + y^2p_y\theta^2j + yzp_z\theta^2j + xzp_x\theta^2k + yzp_y\theta^2k + z^2p_z\theta^2k \\ &+ p_x\gamma^2i + yp_z\gamma\theta i- zp_y\gamma\theta i - yp_x\gamma\theta k - y^2p_z\theta^2k + yzp_y\theta^2k + zp_x\gamma\theta j + yzp_z\theta^2j - z^2p_y\theta^2j \\ &+ p_y\gamma^2j + zp_x\gamma\theta j - xp_z\gamma\theta j + xp_y\gamma\theta k + xzp_x\theta^2 k - x^2p_z\theta^2 k - zp_y\gamma\theta i - z^2p_x\theta^2 i + xzp_z\theta^2 i \\ &+ p_z\gamma^2 k + xp_y\gamma\theta k - yp_x\gamma\theta k - xp_z\gamma\theta j - x^2p_y\theta^2 j + xyp_x\theta^2j + yp_z\gamma\theta i + xyp_y\theta^2 i - y^2p_x\theta^2 i \end{align} $$
which when cleaned up a little becomes
$$ \begin{align} qpq^* &= i[p_x\theta^2(x^2-y^2-z^2) + 2x\theta^2(yp_y+zp_z) + 2\gamma\theta(yp_z-zp_y) \color{red}{+ p_x\gamma^2}] \\ &+ j[p_y\theta^2(y^2-z^2-x^2) + 2y\theta^2(xp_x+zp_z) + 2\gamma\theta(zp_x-xp_z) \color{red}{+ p_y\gamma^2}] \\ &+ k[p_z\theta^2(z^2-y^2-x^2) + 2z\theta^2(xp_x+yp_y) + 2\gamma\theta(xp_y-yp_x) \color{red}{+ p_z\gamma^2}] \end{align} $$
This is the result that I have implemented in Python, where I have taken the $i$ component to be the rotated x coordinate, $j$ component as y coordinates, and so on. Here is the relevant code
```python
def quaternion_rotation(angle=math.pi/2, axis=(1.0, 0.0, 0.0), coordinates=[]):
"""
Perform the quaternion rotation by "angle" radians around "axis"
"""
# Some shorthands
gamma = math.cos(angle/2)
theta = math.sin(angle/2)
x = axis[0] # Axis of rotation, x element
y = axis[1] # Axis of rotation, y element
z = axis[2] # Axis of rotation, z element
# Loop over atomic position vectors and apply rotation
coords_rot = []
for i, atom in enumerate(coordinates):
px = atom[0] # x element of vector to rotate
py = atom[1] # y element of vector to rotate
pz = atom[2] # z element of vector to rotate
# q p q*
px_rot = px*theta**2 * (x**2 - y**2 - z**2) + 2*x*theta**2 * (y*py + z*pz) + 2*gamma*theta * (y*pz -z*py) + px*gamma**2
py_rot = py*theta**2 * (y**2 - z**2 - x**2) + 2*y*theta**2 * (x*px + z*pz) + 2*gamma*theta * (z*px - x*pz) + py*gamma**2
pz_rot = pz*theta**2 * (z**2 - y**2 - x**2) + 2*z*theta**2 * (x*px + y*py) + 2*gamma*theta * (x*py - y*px) + pz*gamma**2
coords_rot.append([px_rot, py_rot, pz_rot])
return coords_rot
```
There may be trivial algebraic errors in the derivation, but I reckon it is my approach that somehow is wrong. Can anyone see why this is incorrect?
Edit:
In the second derivation, I got some factors missing in my initial final expression. These are indicated in red above. It seems they were in the "messy" expression, so I must just have forgotten about them. Here are two images of pre and post-rotated structures (90 degrees). The distance between the same two pairs of atoms are identical in both structures, so it seems a pure spatial rotation has taken place!
Your overall approach is fine. If you’re not getting the correct results, you’ve likely made an algebraic error along the way. The first product looks fine, so something went wrong when you calculated the second one. My guess would be a sign error or two as it’s very easy to get the factors in the wrong order or miss a double negative somewhere. What you should end up with at the end is Rodrigues’ rotation formula, so that’s a good sanity check. When I subtract your result from this formula, I don’t end up with zero. Examining the difference might give you a clue as to where you went wrong in your computation.
Most of the time I find that I make fewer computational errors if I delay expanding by coordinates until the very end. For this derivation, I’d work with a partial expansion of the quaternions into scalar and vector parts. The corresponding rule for multiplication is $$(r_1,\mathbf v_1)(r_2,\mathbf v_2) = (r_1r_2-\mathbf v_1\cdot\mathbf v_2,r_1\mathbf v_2+r_2\mathbf v_1+\mathbf v_1\times\mathbf v_2).$$ Taking the multiplications in the same order as yours, we have $$\left(\cos{\frac\alpha2},\mathbf k\sin{\frac\alpha2}\right)\left(0,\mathbf p\right) = \left(-\sin{\frac\alpha2}(\mathbf k\cdot\mathbf p),\mathbf p\cos{\frac\alpha2} + \sin{\frac\alpha2}(\mathbf k\times\mathbf p)\right),$$ which jibes with what you got. Continuing, we have for the scalar part of $qpq^{-1}$ $$-\cos{\frac\alpha2}\sin{\frac\alpha2}(\mathbf k\cdot\mathbf p) - \left(\mathbf p\cos{\frac\alpha2} + \sin{\frac\alpha2}(\mathbf k\times\mathbf p)\right)\cdot \left(-\mathbf k\sin{\frac\alpha2}\right) = 0$$ and for the vector part $$\left(\sin^2{\frac\alpha2}\right)(\mathbf k\cdot\mathbf p)\mathbf k + \left(\cos^2{\frac\alpha2}\right)\mathbf p + \left(\sin\alpha\right)\mathbf k\times\mathbf p - \left(\sin^2{\frac\alpha2}\right)(\mathbf k\times\mathbf p)\times\mathbf k.$$ This can be further simplified using $\mathbf k\times(\mathbf k\times\mathbf p) = \mathbf k(\mathbf k\cdot\mathbf p)-\mathbf p(\mathbf k\cdot\mathbf k) = \mathbf k(\mathbf k\cdot\mathbf p)-\mathbf p$ and $\cos\alpha = \cos^2(\alpha/2)-\sin^2(\alpha/2)$ to obtain Rodrigues’ rotation formula as mentioned earlier.