Quaternion rotation in 3D

The purpose of this notebook was to explore the behaviour of quaternion rotation in 3 dimensional (3D) space.

Start with a 3D unit vector u = (u_x, u_y, u_z) - note a unit vector has length 1. This unit vector represents the axis of rotation. In the chemical context that I will use the unit vector, this unit vector lies along the bc bond, pointing from atom b to atom c. The purpose of learning rotations is to apply rotations to atoms in 3D space to modify dihedral angles. A dihedral angle is described by the angle between 2 planes, one plane formed by atoms a, b, and c, and the other plane formed by atoms b, c, and d.

a - b - c - d

Atoms connected to c (directly or indirectly), except those connected through the b-c bond, will need to rotate to establish a new dihedral angle.

Note that if the bc bond is part of a ring, then this method of rotation should not be used (i.e. b is connected to c in another route other than bc bond).

Here are the basic steps that need to be followed:

1. Determine if bc bond is contained within a ring. If the bond is inside a ring, cancel calculation and use another algorithm.
2. Calculate the existing dihedral angle.
3. Calculate how much rotation is required to get to the new dihedral angle, which is calculated by subtraction. If the amount of rotation needed is zero, quit the algorithm.
4. Determine which atoms need to be rotated as a result of the dihedral change. These atoms should be connected to atom c (directly or indirectly), except those connected through the b-c bond.
5. Determine the axis of rotation. This axis is the unit vector in the direction of the b-c bond, pointing from b to c, and has length 1.
6. For each atom that needs to be rotated, it must be translated relative to atom b. Atom b in this case is treated as the "origin". For example, if atom b is at point (2,4,1) and atom x is atom point (4,-1,2), the vector BX is (4-2,-1-4,2-1) = (2,-5,1).
7. Using quaternion rotation, apply the formula q p q-1, where p-1 is the inverse of q.
8. The resulting vector from step 7 will point in the correct direction and have the correct magnitude; however, it will still be relative to the origin (0,0,0) and not atom b. Therefore, the last step is to add the two vectors together to get the final point. For example, if the result of step 7 for atom x is (1, 0.5, 0) and atom b is located at (5, 6, 2), the new location for atom x after rotation is (1+5, 0.5+6, 0+2).
9. Steps 6, 7, and 8 are repeated for each atom identified for rotation in step 4.

Notes on contents of notebook

The purpose of the notebook was to identify if the translation steps (6 and 8) were required. It was determined by visual inspection (using an online 3D plotting tool - https://www.geogebra.org/3d?lang=en) that the results were incorrect unless translation was applied for both stages. More explicitly, the first plot_atoms function (see below) placed atom b at the origin and atom c was along the x axis. The calculations were performed to determine the new position for 2 points, which were not translated. The result was correct. However, when this same procedure was used in the second version of plot_atoms (called plot_atoms_new), the positions of the rotated atoms were incorrect. The difference between plot_atoms and plot_atoms_new was that for the new version the entire "molecule" was translated to a random position. Therefore, atom b was no longer at the origin.

Next test required

Perform the rotation on a new set of atomic coordinates that have been translated and do not use a primary axis (i.e. the x, y, or z axis) for the b-c bond. After visualising the results, they look okay as well.

More on quaternions, aka step 7

A quaternion is a special type of number that has 1 real component and 3 imaginary components.

quaternion = [a, bi, cj, dk]

Where a,b,c,d are real numbers, and i^2 = j^2 = k^2 = -1.

To calculate the result of rotation, we use the following formula:

result = q p q-1

Where q-1 is q inverse. q and p-1 are unit quaternions, represented by the following equations:

q = [cos(0.5#), (u_x, u_y, u_z) * sin(0.5#)]

q-1 = [cos(0.5#), -(u_x, u_y, u_z) * sin(0.5#)]

||q|| = 1

p is a vector (p_x, p_y, p_z); in this representation as a quaternion it has a=0.

q p q-1 gives a quaternion with a zero real component and the i,j,k component being equal to the result after rotating by angle # about the axis described by vector u.

Quaternion multiplication follows a special formula, as follows:

[a,u] [b,v] = [a b - u v, (a v + b * u) + u x b]

Note that a and b are real numbers, whereas u and v are vectors with 3 components. The x in this formula means performing a cross product. In the example code below, the cross product is calculated by using the formula for a 3 by 3 matrix determinant. u * v indicates dot product.

Helpful lecture here: https://www.youtube.com/watch?v=rJqII1FyrdM

from math import pi, cos, sin, sqrt
from numpy import allclose
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def vector_magnitude(vector):
    """
    The magnitude of a vector (or length) is denoted:
    ||v|| = sqrt (x_0 ^ 2 + x_1 ^ 2 + ... + x_n-1 ^ 2)
    """
    result = 0
    for element in vector:
        result += element ** 2
    return sqrt(result)


def unit_vector(vector):
    """
    unit vector formula:
    v / ||v||
    
    a unit vector has a magnitude of 1
    """
    magnitude = vector_magnitude(vector)
    result = []
    for element in vector:
        result.append(element / magnitude)
    assert allclose(vector_magnitude(result), 1.0)
    return result


def cross_product(v1, v2):
    """calculate cross product b/w 2 3-element vectors"""
    a1 = v1[0]
    a2 = v1[1]
    a3 = v1[2]
    
    b1 = v2[0]
    b2 = v2[1]
    b3 = v2[2]
    
    return [
        a2 * b3 - a3 * b2,
        -1 * (a1 * b3 - a3 * b1),
        a1 * b2 - a2 * b1
    ]

def plot_atoms():
    atoms = [
        [0, 0, 0],
        [1, 0, 0],
        [1, 0.5, 0],
        [1, 1, 0],
        [-1, 0.5, 0]
    ]
    x = [atom[0] for atom in atoms]
    y = [atom[1] for atom in atoms]
    z = [atom[2] for atom in atoms]
    fig = plt.figure()
    ax = fig.gca(projection="3d")
    ax.plot(x, y, z, "o")
    #plt.show()
    
    atom_a = [0, 0, 0]
    atom_b = [1, 0, 0]
    
    bond_ab = [b - a for a, b in zip(atom_a, atom_b)]
    unit_ab = unit_vector(bond_ab)
    
    angle = 3*pi/2 # radians
    
    # quaternion q
    q = [
        cos(0.5 * angle),
        unit_ab[0] * sin(0.5 * angle),
        unit_ab[1] * sin(0.5 * angle),
        unit_ab[2] * sin(0.5 * angle)
    ]
    
    # quaternion q inverse
    qi = [
        cos(0.5 * angle),
        -1 * unit_ab[0] * sin(0.5 * angle),
        -1 * unit_ab[1] * sin(0.5 * angle),
        -1 * unit_ab[2] * sin(0.5 * angle)
    ]
    
    # real part of these quaternions is 0
    p_quaternions = [[0, 1, 0.5, 0], [0, 1, 1, 0]]
    
    # quaternion multiplication
    # q1 * q2 = [a, u] * [b, v]
    #         = [a*b - u*v, (a*v + b*u) + u x v]
    # where x is cross product, a and b are scalars
    # and u and v are 3 component vectors
    
    for p in p_quaternions:
        # q * p * q-1
        # q * p first
        qp_real = q[0] * p[0] - sum([u*v for u,v in zip(q[1:], p[1:])])
        qp_img_1 = [q[0] * element for element in p[1:]]
        qp_img_2 = [p[0] * element for element in q[1:]]
        qp_img_3 = [a + b for a,b in zip(qp_img_1, qp_img_2)]
        qp_img_4 = cross_product(q[1:], p[1:])
        qp_img = [a + b for a,b in zip(qp_img_3, qp_img_4)]
        
        qp = [qp_real] + qp_img
        
        # qp * p inverse
        qpqi_real = qp[0] * qi[0] - sum([u*v for u,v in zip(qp[1:], qi[1:])])
    
        qpqi_img_1 = [qp[0] * element for element in qi[1:]]
        qpqi_img_2 = [qi[0] * element for element in qp[1:]]
        qpqi_img_3 = [a + b for a,b in zip(qpqi_img_1, qpqi_img_2)]
        qpqi_img_4 = cross_product(qp[1:], qi[1:])
        qpqi_img = [a + b for a,b in zip(qpqi_img_3, qpqi_img_4)]
        
        qpqi = [qpqi_real] + qpqi_img
        assert allclose(qpqi[0], 0.0)
        ax.plot([qpqi[1]], [qpqi[2]], [qpqi[3]], "or")
        print(qpqi)
    
    

plot_atoms()

[0.0, 1.0, -1.1102230246251565e-16, -0.5]
[0.0, 1.0, -2.220446049250313e-16, -1.0]

def plot_atoms_new():
    atoms = [
        [5, 2, -2],
        [6, 2, -2],
        [6, 2.5, -2],
        [6, 3, -2],
        [4, 2.5, -2]
    ]
    x = [atom[0] for atom in atoms]
    y = [atom[1] for atom in atoms]
    z = [atom[2] for atom in atoms]
    fig = plt.figure()
    ax = fig.gca(projection="3d")
    ax.plot(x, y, z, "o")
    #plt.show()
    
    atom_a = [5, 2, -2]
    atom_b = [6, 2, -2]
    
    bond_ab = [b - a for a, b in zip(atom_a, atom_b)]
    unit_ab = unit_vector(bond_ab)
    
    angle = 3*pi/2 # radians
    
    # quaternion q
    q = [
        cos(0.5 * angle),
        unit_ab[0] * sin(0.5 * angle),
        unit_ab[1] * sin(0.5 * angle),
        unit_ab[2] * sin(0.5 * angle)
    ]
    
    # quaternion q inverse
    qi = [
        cos(0.5 * angle),
        -1 * unit_ab[0] * sin(0.5 * angle),
        -1 * unit_ab[1] * sin(0.5 * angle),
        -1 * unit_ab[2] * sin(0.5 * angle)
    ]
    
    # real part of these quaternions is 0
    # need to calculate them as vectors using atom_a
    # as "origin"
    # p_quaternions = [[0, 6, 2.5, -2], [0, 6, 3, -2]]
    
    p_quaternions = [
        [0, 6 - atom_a[0], 2.5 - atom_a[1], -2 - atom_a[2]],
        [0, 6 - atom_a[0], 3 - atom_a[1], -2 - atom_a[2]]
    ]
    
    # quaternion multiplication
    # q1 * q2 = [a, u] * [b, v]
    #         = [a*b - u*v, (a*v + b*u) + u x v]
    # where x is cross product, a and b are scalars
    # and u and v are 3 component vectors
    
    for p in p_quaternions:
        # q * p * q-1
        # q * p first
        qp_real = q[0] * p[0] - sum([u*v for u,v in zip(q[1:], p[1:])])
        qp_img_1 = [q[0] * element for element in p[1:]]
        qp_img_2 = [p[0] * element for element in q[1:]]
        qp_img_3 = [a + b for a,b in zip(qp_img_1, qp_img_2)]
        qp_img_4 = cross_product(q[1:], p[1:])
        qp_img = [a + b for a,b in zip(qp_img_3, qp_img_4)]
        
        qp = [qp_real] + qp_img
        
        # qp * p inverse
        qpqi_real = qp[0] * qi[0] - sum([u*v for u,v in zip(qp[1:], qi[1:])])
    
        qpqi_img_1 = [qp[0] * element for element in qi[1:]]
        qpqi_img_2 = [qi[0] * element for element in qp[1:]]
        qpqi_img_3 = [a + b for a,b in zip(qpqi_img_1, qpqi_img_2)]
        qpqi_img_4 = cross_product(qp[1:], qi[1:])
        qpqi_img = [a + b for a,b in zip(qpqi_img_3, qpqi_img_4)]
        
        qpqi = [qpqi_real] + qpqi_img
        assert allclose(qpqi[0], 0.0)
        # ax.plot([qpqi[1]], [qpqi[2]], [qpqi[3]], "or")
        # need to add change existing atom_a
        qpqi[1] += atom_a[0]
        qpqi[2] += atom_a[1]
        qpqi[3] += atom_a[2]
        ax.plot([qpqi[1]], [qpqi[2]], [qpqi[3]], "or")
        print(qpqi)

plot_atoms_new()

[0.0, 6.0, 2.0, -2.5]
[0.0, 6.0, 1.9999999999999998, -3.0]

def plot_atoms_test():
    atoms = [
        [5, 2, -2], # b
        [7.5, 6, 5],# c
        [6, 2.5, -1],# connected to c
        [6, 3, -5],  # connected to c
        [4, 4, -3] # connected to b
    ]
    x = [atom[0] for atom in atoms]
    y = [atom[1] for atom in atoms]
    z = [atom[2] for atom in atoms]
    fig = plt.figure()
    ax = fig.gca(projection="3d")
    ax.plot(x, y, z, "o")
    
    # these 2 atoms make up the bc bond.
    atom_a = atoms[0]
    atom_b = atoms[1]
    
    bond_ab = [b - a for a, b in zip(atom_a, atom_b)]
    unit_ab = unit_vector(bond_ab)
    
    angle = 3*pi/2 # radians
    
    # quaternion q
    q = [
        cos(0.5 * angle),
        unit_ab[0] * sin(0.5 * angle),
        unit_ab[1] * sin(0.5 * angle),
        unit_ab[2] * sin(0.5 * angle)
    ]
    
    # quaternion q inverse
    qi = [
        cos(0.5 * angle),
        -1 * unit_ab[0] * sin(0.5 * angle),
        -1 * unit_ab[1] * sin(0.5 * angle),
        -1 * unit_ab[2] * sin(0.5 * angle)
    ]
    
    # real part of these quaternions is 0
    # need to calculate them as vectors using atom_a
    # as "origin"
    # p_quaternions = [[0, 6, 2.5, -2], [0, 6, 3, -2]]
    
    p_quaternions = [
        [0, atoms[2][0] - atom_a[0], atoms[2][1] - atom_a[1], atoms[2][2] - atom_a[2]],
        [0, atoms[3][0] - atom_a[0], atoms[3][1] - atom_a[1], atoms[3][2] - atom_a[2]]
    ]
    
    # quaternion multiplication
    # q1 * q2 = [a, u] * [b, v]
    #         = [a*b - u*v, (a*v + b*u) + u x v]
    # where x is cross product, a and b are scalars
    # and u and v are 3 component vectors
    
    for p in p_quaternions:
        # q * p * q-1
        # q * p first
        qp_real = q[0] * p[0] - sum([u*v for u,v in zip(q[1:], p[1:])])
        qp_img_1 = [q[0] * element for element in p[1:]]
        qp_img_2 = [p[0] * element for element in q[1:]]
        qp_img_3 = [a + b for a,b in zip(qp_img_1, qp_img_2)]
        qp_img_4 = cross_product(q[1:], p[1:])
        qp_img = [a + b for a,b in zip(qp_img_3, qp_img_4)]
        
        qp = [qp_real] + qp_img
        
        # qp * p inverse
        qpqi_real = qp[0] * qi[0] - sum([u*v for u,v in zip(qp[1:], qi[1:])])
    
        qpqi_img_1 = [qp[0] * element for element in qi[1:]]
        qpqi_img_2 = [qi[0] * element for element in qp[1:]]
        qpqi_img_3 = [a + b for a,b in zip(qpqi_img_1, qpqi_img_2)]
        qpqi_img_4 = cross_product(qp[1:], qi[1:])
        qpqi_img = [a + b for a,b in zip(qpqi_img_3, qpqi_img_4)]
        
        qpqi = [qpqi_real] + qpqi_img
        assert allclose(qpqi[0], 0.0)
        # ax.plot([qpqi[1]], [qpqi[2]], [qpqi[3]], "or")
        # need to add change existing atom_a
        qpqi[1] += atom_a[0]
        qpqi[2] += atom_a[1]
        qpqi[3] += atom_a[2]
        ax.plot([qpqi[1]], [qpqi[2]], [qpqi[3]], "or")
        print(qpqi)

plot_atoms_test()

[0.0, 5.344273884153915, 2.1125000451045364, -0.5443835558289907]
[0.0, 6.7421538056599895, -0.5318468332206661, -3.2468567401810438]