Rotations with scipy.spatial.transform: superpose two conformers

The problem: two turns that do not commute

Lay a book flat on the desk, turn it 90° about the x axis, then 90° about the y axis. Take a second copy and do the same two turns in the other order. The two books end up facing different ways, and the turn that takes one onto the other is itself a rotation: 120° about the axis \((1, -1, 1)/\sqrt{3}\). Rotations in three dimensions do not commute, and most bugs in code that turns things come from forgetting it. SciPy's answer is one object, scipy.spatial.transform.Rotation, that you build from Euler angles, a matrix, an axis and angle, or a quaternion, and that gives back any of them. Each representation is explained in the step that first needs it.

The practical end is a chemist's daily chore: two conformers of ibuprofen, 33 atoms with hydrogens, one as RDKit generated it, the other turned into a random orientation and shifted by (6, −3, 2) Å, as if it came from another program. How different are their shapes once position and orientation are taken out? The usual measure is the root-mean-square deviation,

\[\mathrm{RMSD} = \sqrt{\frac{1}{N}\sum_{i=1}^{N} \lVert \mathbf a_i - \mathbf b_i \rVert^2},\]

with \(\mathbf a_i\) and \(\mathbf b_i\) the positions of atom \(i\) in the two conformers and \(N = 33\). The sum pairs atom \(i\) with atom \(i\), so it only means something when both lists name the same atoms in the same order.

Two ibuprofen conformers projected on the first two principal axes of the reference, in Å. Left: as placed, B off to the side and turned, RMSD 8.53 Å. Right: superposed, the skeletons follow each other but the ring is tilted and the acid group at the right end turned apart, RMSD 2.23 Å.

This is where we end up: before the superposition the RMSD is 8.53 Å, after it 2.23 Å, and what remains is a real difference of shape. Rotation.align_vectors found the turn. The rest of this tutorial builds that figure.

Setup

RDKit installs with pip install rdkit. Three of its calls make the molecule: Chem.MolFromSmiles reads the structure from a SMILES string, Chem.AddHs adds the hydrogens that SMILES leaves implicit, and AllChem.EmbedMultipleConfs generates 3D coordinates for two conformers, reproducibly with a seed. RDKit from the ground up does the rest of what RDKit can do. The rmsd function is the formula above.

import warnings
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.transform import Rotation
from rdkit import Chem
from rdkit.Chem import AllChem, rdMolAlign

plt.rcParams.update({
    "figure.figsize": (7, 3.6), "figure.dpi": 110,
    "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25,
    "font.size": 11, "lines.linewidth": 1.8,
})
INK, ACCENT, SECOND, MUTED = "#1f2a44", "#c8553d", "#2a7f9e", "#8a8f98"
SEED = 42

mol = Chem.AddHs(Chem.MolFromSmiles("CC(C)Cc1ccc(cc1)C(C)C(=O)O"))   # ibuprofen
AllChem.EmbedMultipleConfs(mol, numConfs=2, randomSeed=SEED)
A = mol.GetConformer(0).GetPositions()      # the reference, in Å
B0 = mol.GetConformer(1).GetPositions()     # the second conformer, before we move it

def rmsd(a, b):
    return np.sqrt(np.mean(np.sum((a - b) ** 2, axis=1)))

print(f"{mol.GetNumAtoms()} atoms, A {A.shape}, B0 {B0.shape}, RMSD(A, B0) = {rmsd(A, B0):.2f} Å")

33 atoms, A (33, 3), B0 (33, 3), RMSD(A, B0) = 2.60 Å

The two conformers differ by 2.60 Å before anything is turned, so there is a difference of shape for the superposition to find.

Step 1: Build a rotation four ways and convert it

A quaternion is four numbers that encode a turn by the angle \(\theta\) about the unit axis \(\mathbf n\):

\[q = (x, y, z, w) = \left(n_x \sin\tfrac{\theta}{2},\; n_y \sin\tfrac{\theta}{2},\; n_z \sin\tfrac{\theta}{2},\; \cos\tfrac{\theta}{2}\right).\]

The first three, \(x\), \(y\), \(z\), carry the axis, the fourth, \(w\), the amount of turning, and their squares sum to 1, so a unit quaternion is a point on a sphere in four dimensions. SciPy stores every rotation this way because composing two quaternions is cheap and has no special angles, which Euler angles do have (Step 3). A rotation vector is the axis times the angle in radians. Here is the same 90° turn about x, built four ways:

Rx = Rotation.from_euler("x", 90, degrees=True)
same = [Rotation.from_rotvec([np.pi / 2, 0, 0]),
        Rotation.from_quat([np.sin(np.pi / 4), 0, 0, np.cos(np.pi / 4)]),
        Rotation.from_matrix([[1, 0, 0], [0, 0, -1], [0, 1, 0]])]
print("all equal:", all(r.approx_equal(Rx) for r in same))
print("quaternion     ", Rx.as_quat().round(4))
print("rotation vector", Rx.as_rotvec().round(4))
print("matrix\n", Rx.as_matrix().round(4))
print("y axis goes to ", Rx.apply([0, 1, 0]).round(4))

all equal: True
quaternion      [0.7071 0.     0.     0.7071]
rotation vector [1.5708 0.     0.    ]
matrix
 [[ 1.  0.  0.]
 [ 0.  0. -1.]
 [ 0.  1.  0.]]
y axis goes to  [0. 0. 1.]

The quaternion is \((\sin 45°, 0, 0, \cos 45°)\), as the formula says, with \(w\), the one component without an axis in it, last. apply turns vectors: the y axis lands on z, which is what a right-handed 90° turn about x does.

Step 2: Compose two turns and find the difference

The product Ry * Rx is the rotation that applies Rx first, read from right to left like matrices. Let the book's length, width, and thickness lie along x, y, and z, and follow the three edges through both orders:

Ry = Rotation.from_euler("y", 90, degrees=True)
edges = np.eye(3)                            # length, width, thickness of the book
for label, r in [("x, then y", Ry * Rx), ("y, then x", Rx * Ry)]:
    moved = r.apply(edges).round(3) + 0.0    # + 0.0 turns -0 into 0
    print(f"{label}:  length -> {moved[0]}  width -> {moved[1]}  thickness -> {moved[2]}")

D = (Rx * Ry) * (Ry * Rx).inv()              # undo the first order, do the second
print(f"difference: {np.degrees(D.magnitude()):.1f}° about {(D.as_rotvec() / D.magnitude()).round(3)}")

x, then y:  length -> [ 0.  0. -1.]  width -> [1. 0. 0.]  thickness -> [ 0. -1.  0.]
y, then x:  length -> [0. 1. 0.]  width -> [0. 0. 1.]  thickness -> [1. 0. 0.]
difference: 120.0° about [ 0.577 -0.577  0.577]

Not one edge lands in the same place. The difference rotation D undoes the first order with .inv() and applies the second; magnitude is its angle and the normalized rotation vector its axis, 120.0° about \((0.577, -0.577, 0.577)\). Drawn, the book turned x first stands on its short edge with the cover toward −y, the other lies on its long edge with the cover toward +x:

Show code

from mpl_toolkits.mplot3d.art3d import Poly3DCollection

def draw_book(ax, r, label):
    half = np.array([1.5, 1.0, 0.2])                   # a 3 x 2 x 0.4 book centered at the origin
    corners = np.array([[i, j, k] for i in (-1, 1) for j in (-1, 1) for k in (-1, 1)]) * half
    c = r.apply(corners)
    faces = {"cover": [1, 3, 7, 5], "back": [0, 2, 6, 4], "spine": [0, 1, 5, 4],
             "edge": [2, 3, 7, 6], "top": [0, 1, 3, 2], "bottom": [4, 5, 7, 6]}
    polys = [c[idx] for idx in faces.values()]
    fill = [SECOND if name == "cover" else "white" for name in faces]
    ax.add_collection3d(Poly3DCollection(polys, facecolors=fill, edgecolors=INK, linewidths=0.8))
    for k, name in enumerate("xyz"):
        tip = [2.2, 2.6, 1.9][k] * np.eye(3)[k]       # lengths that keep each letter clear of the book
        ax.plot(*zip(np.zeros(3), tip), color=MUTED, lw=1)
        ax.text(*(tip * 1.25), name, color=MUTED, ha="center", va="center")
    ax.set(xlim=(-2, 2), ylim=(-2, 2), zlim=(-1.6, 1.6))
    ax.set_box_aspect((1, 1, 0.8), zoom=1.1)          # equal scale on all three axes
    ax.set_axis_off()
    ax.text2D(0.04, 0.97, label, transform=ax.transAxes, color=INK, va="top")

fig, axes = plt.subplots(1, 3, figsize=(7.5, 2.4), subplot_kw=dict(projection="3d"))
fig.subplots_adjust(left=0, right=1, bottom=0, top=1, wspace=0)
for ax, r, label in zip(axes, [Rotation.identity(), Ry * Rx, Rx * Ry],
                        ["start", "x, then y", "y, then x"]):
    draw_book(ax, r, label)
plt.show()

A flat book in three 3D panels with fixed x, y, z axes: at the start, after 90° about x then 90° about y, and after the same turns in the other order. The blue cover faces up at the start, toward −y after x then y, and toward +x after y then x.

Euler angles describe a rotation as three successive turns about coordinate axes, and the axis string names those axes in order. Lowercase letters mean axes fixed in space (extrinsic), uppercase letters axes that turn with the body (intrinsic):

print('"xyz" [90, 90, 0] equals Ry * Rx:', Rotation.from_euler("xyz", [90, 90, 0], degrees=True).approx_equal(Ry * Rx))
print('"XYZ" [90, 90, 0] equals Rx * Ry:', Rotation.from_euler("XYZ", [90, 90, 0], degrees=True).approx_equal(Rx * Ry))

"xyz" [90, 90, 0] equals Ry * Rx: True
"XYZ" [90, 90, 0] equals Rx * Ry: True

The same three numbers give both books, depending on the case of three letters. Turning about body axes in the order X, Y is the same as turning about fixed axes in the reverse order, y then x, which is why "XYZ" gives Rx * Ry.

Step 3: Euler angles at 90° pitch, and q against −q

For the string "ZYX" the three angles are yaw (the heading), pitch (nose up or down), and roll (bank), as for an aircraft. Pitch is the middle one. Pitch by 90°, so that the nose swings onto the z axis, and ask for the angles back:

r = Rotation.from_euler("ZYX", [30, 90, 10], degrees=True)
with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    print("angles back:", r.as_euler("ZYX", degrees=True).round(3))
print("warning:    ", caught[0].message)
print("[40, 90, 20] is the same rotation:", Rotation.from_euler("ZYX", [40, 90, 20], degrees=True).approx_equal(r))

angles back: [20. 90.  0.]
warning:     Gimbal lock detected. Setting third angle to zero since it is not possible to uniquely determine all angles.
[40, 90, 20] is the same rotation: True

You put in 30° yaw and 10° roll and get back 20° and 0°. At 90° pitch the roll axis, the body's x axis, points straight down along −z, opposite to the yaw axis. Yaw and roll then turn about one line in opposite senses, so only their difference counts: 30° − 10°, 20° − 0°, and 40° − 20° are the same 20°. One of the three angles is lost: gimbal lock. SciPy warns and picks one answer. The fix is to keep Rotation objects or quaternions throughout and convert to Euler angles only for display.

Quaternions have their own ambiguity, a harmless one:

q = Rx.as_quat()
flipped = Rotation.from_quat(-q)
print("-q is the same rotation:", flipped.approx_equal(Rx))
print("as_quat()              ", flipped.as_quat().round(4) + 0.0)
print("as_quat(canonical=True)", flipped.as_quat(canonical=True).round(4))

-q is the same rotation: True
as_quat()               [-0.7071  0.      0.     -0.7071]
as_quat(canonical=True) [0.7071 0.     0.     0.7071]

Adding 360° to \(\theta\) is the same turn, but it adds 180° to \(\theta/2\), where sine and cosine both change sign: every component of \(q\) flips. Two quaternions per rotation, so compare rotations with approx_equal, never quaternions element by element. canonical=True picks the one with \(w \ge 0\).

Step 4: Throw a conformer into a random orientation

Rotation.random draws an orientation uniformly from all possible ones. Turn conformer B with it and shift it, as if it came from another program:

turn = Rotation.random(rng=np.random.default_rng(1))
shift = np.array([6.0, -3.0, 2.0])           # Å
B = turn.apply(B0) + shift
print(f"random turn {np.degrees(turn.magnitude()):.1f}°,  RMSD(A, B) as placed = {rmsd(A, B):.2f} Å")

Ac = A - A.mean(axis=0)                       # each set about its own centroid
Bc = B - B.mean(axis=0)
print(f"RMSD after centering = {rmsd(Ac, Bc):.2f} Å")

random turn 72.2°,  RMSD(A, B) as placed = 8.53 Å
RMSD after centering = 4.77 Å

The turn is 72.2°, and the RMSD as placed is 8.53 Å, more than three times the 2.60 Å of shape difference. Subtracting each set's centroid removes the shift and brings it to 4.77 Å. What is left is orientation and shape.

Step 5: Superpose with align_vectors

Finding the best turn is the Kabsch problem (Kabsch 1976): for centered points, find the rotation \(R\) that minimizes

\[\sum_{i=1}^{N} \lVert \mathbf a_i - R\,\mathbf b_i \rVert^2 .\]

Rotation.align_vectors(a, b) solves it and returns the rotation that turns b onto a, here B onto A, so the reference goes first. It pairs row \(i\) of a with row \(i\) of b: both arrays must list the same atoms in the same order. Two conformers of one RDKit molecule do. Two files from different programs often do not and must be reordered first; RDKit's GetBestRMS also tries symmetry-equivalent orders.

rot, rssd = Rotation.align_vectors(Ac, Bc)
N = len(Ac)
print(f"RMSD from rssd   = {rssd / np.sqrt(N):.4f} Å")
print(f"RMSD by hand     = {rmsd(Ac, rot.apply(Bc)):.4f} Å")

RMSD from rssd   = 2.2269 Å
RMSD by hand     = 2.2269 Å

The second return value, rssd, is the root of the sum of squared distances after the fit, so dividing by \(\sqrt{N}\) gives the RMSD: 2.2269 Å, the same as the formula applied to the turned coordinates. To put the superposed B back into A's original frame, add A's centroid: rot.apply(Bc) + A.mean(axis=0).

RDKit does the same superposition in one call. Write B into the second conformer of a copy and align it to the first:

probe = Chem.Mol(mol)
probe.GetConformer(1).SetPositions(B)
print(f"rdMolAlign.AlignMol: {rdMolAlign.AlignMol(probe, mol, prbCid=1, refCid=0):.4f} Å")

rdMolAlign.AlignMol: 2.2269 Å

The same 2.2269 Å. It is not noise: the two conformers really differ, and no rotation removes a difference of shape.

Step 6: Project both conformers on their principal axes

A molecule is three-dimensional and the page is not. The principal axes are the directions along which the reference conformer is most extended, the longest first; np.linalg.svd of the centered coordinates returns them as the rows of Vt (the same axes as in PCA with scikit-learn). Project both conformers onto the first two, before and after, with the bonds drawn from mol.GetBonds():

_, _, Vt = np.linalg.svd(Ac)
axes2 = Vt[:2].T                              # columns: the two longest principal axes
heavy = np.array([atom.GetAtomicNum() > 1 for atom in mol.GetAtoms()])
bonds = [(b.GetBeginAtomIdx(), b.GetEndAtomIdx()) for b in mol.GetBonds()]

def draw(ax, xyz, color):
    p = xyz @ axes2
    for i, j in bonds:
        ax.plot(p[[i, j], 0], p[[i, j], 1], color=color, lw=1.2)
    ax.plot(*p[heavy].T, "o", color=color, ms=6)
    ax.plot(*p[~heavy].T, "o", color=color, ms=3)

fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(7.5, 3.6), sharex=True, sharey=True)
for ax, a, b in [(ax0, A, B), (ax1, Ac, rot.apply(Bc))]:
    draw(ax, a, INK)
    draw(ax, b, ACCENT)
    ax.text(0.03, 0.96, f"RMSD {rmsd(a, b):.2f} Å", color=INK, transform=ax.transAxes, va="top")
    ax.text(0.90, 0.96, "A", color=INK, transform=ax.transAxes, va="top", ha="right")
    ax.text(0.97, 0.96, "B", color=ACCENT, transform=ax.transAxes, va="top", ha="right")
    ax.set(xlabel="PC 1 / Å", aspect="equal")
ax0.set(ylabel="PC 2 / Å", ylim=(-5.5, 5.5))   # room above the molecules for the labels
plt.show()

Two ibuprofen conformers projected on the first two principal axes of the reference, PC 1 and PC 2 in Å. Left: as placed, B off to the side and turned, RMSD 8.53 Å. Right: superposed, the skeletons follow each other but the ring is tilted and the acid group at the right end turned apart, RMSD 2.23 Å.

Before, B sits off to the side and turned. After, the two skeletons follow each other from end to end, though not atom on atom, and they part most at the right end. Where the remaining 2.23 Å sit, atom by atom (RDKit numbers the atoms in SMILES order, hydrogens last):

d = np.linalg.norm(Ac - rot.apply(Bc), axis=1)          # per-atom distance after the fit
groups = {"isobutyl C": range(0, 4), "ring C": range(4, 10), "alpha C, methyl": range(10, 12),
          "carboxyl C, O": range(12, 15), "hydrogens": range(15, 33)}
for name, idx in groups.items():
    k = max(idx, key=lambda i: d[i])
    print(f"{name:15s} largest distance {d[k]:.2f} Å  (atom {k:2d}, {mol.GetAtomWithIdx(k).GetSymbol()})")

rot_h, rssd_h = Rotation.align_vectors(Ac[heavy] - Ac[heavy].mean(axis=0), Bc[heavy] - Bc[heavy].mean(axis=0))
print(f"length along PC 1 {np.ptp(Ac @ Vt[0]):.1f} Å,  RMSD of the {heavy.sum()} heavy atoms alone {rssd_h / np.sqrt(heavy.sum()):.2f} Å")

isobutyl C      largest distance 1.22 Å  (atom  3, C)
ring C          largest distance 1.95 Å  (atom  5, C)
alpha C, methyl largest distance 1.83 Å  (atom 11, C)
carboxyl C, O   largest distance 3.76 Å  (atom 13, O)
hydrogens       largest distance 4.04 Å  (atom 32, H)
length along PC 1 10.2 Å,  RMSD of the 15 heavy atoms alone 1.58 Å

Every isobutyl carbon lies within 1.22 Å of its partner, the ring is tilted against its partner by up to 1.95 Å, and the carboxyl oxygen atom 13 is 3.76 Å from its own: the acid group is turned differently in the two conformers. The hydrogens move most of all. In a molecule 10.2 Å long, 2.23 Å is a fifth of the length, a real change of shape. Fitted and measured on the 15 heavy atoms alone, the RMSD is 1.58 Å, so always say whether the hydrogens are in when you quote one.

Pitfalls

Lowercase and uppercase axis strings. Angles copied from a paper or from another program reproduce the wrong orientation, with no error. The source uses the other convention, and Step 2 showed that "xyz" and "XYZ" with the same numbers are different rotations. Find out which convention the source uses, test it once against an orientation you know from a picture, and put it in the variable name: euler_XYZ_deg.

Forgetting to center. align_vectors finds only a rotation about the origin. Fed the raw coordinates, it still returns a rotation and no error:

rot_raw, _ = Rotation.align_vectors(A, B)
print(f"without centering: RMSD = {rmsd(A, rot_raw.apply(B)):.2f} Å")

without centering: RMSD = 7.41 Å

7.41 Å instead of 2.23 Å, because the best rotation about the origin is mostly busy with the shift. Subtract the centroids of both sets first, and add the reference centroid back afterwards.

The order of the quaternion components. SciPy writes \((x, y, z, w)\), scalar last; many papers and other libraries write \((w, x, y, z)\). A quaternion copied across turns into a different rotation, again with no error:

q_wxyz = Rx.as_quat(scalar_first=True)
wrong = Rotation.from_quat(q_wxyz)            # read as if scalar last
print(f"{q_wxyz.round(4)} read scalar-last: {np.degrees(wrong.magnitude()):.0f}° about {(wrong.as_rotvec() / wrong.magnitude()).round(3)}")

[0.7071 0.7071 0.     0.    ] read scalar-last: 180° about [0.707 0.707 0.   ]

The 90° turn about x became a 180° turn about \((1, 1, 0)/\sqrt{2}\). Pass scalar_first=True to from_quat and as_quat whenever the other side writes the scalar first.

Variations

Weight the atoms by mass. Center on the center of mass and pass weights=masses to align_vectors; the heavy atoms then dominate the fit. Report the plain RMSD anyway when you compare with the literature or with AlignMol, because that is what most papers quote.
Average a cluster of orientations. Repeated measurements of one crystal's orientation, or the frames of a trajectory: stack them in one Rotation with Rotation.concatenate and call .mean(). It is meaningful for a cluster, not for orientations spread over all directions.
Interpolate between two orientations. Slerp(times, rotations) from the same module gives the turns in between along the shortest path, for an animation or a camera path.
The attitude of a satellite from two star directions. This is Wahba's problem: the same align_vectors, with two unit vectors measured in the spacecraft's frame and the same two from the star catalog, weighted by their accuracy.

Cheat sheet

r = Rotation.from_euler("xyz", [a, b, c], degrees=True)  # lowercase: fixed axes; "XYZ": body axes
r = Rotation.from_quat(q, scalar_first=False)            # SciPy default (x, y, z, w)
p * r                                                    # r first, then p
r.inv(), r.apply(v), r.magnitude(), r.as_rotvec()        # inverse, turn vectors, angle (rad), axis * angle
r.as_quat(canonical=True)                                # q and -q are the same turn; this picks w >= 0
r.approx_equal(other)                                    # compare rotations, never quaternions
Rotation.random(rng=np.random.default_rng(1))            # uniform random orientation
rot, rssd = Rotation.align_vectors(a - a.mean(0), b - b.mean(0))  # turns b onto a; row i pairs atom i
rmsd = rssd / np.sqrt(len(a))                            # superposed b in a's frame: rot.apply(b - b.mean(0)) + a.mean(0)