Measure and classify ring puckering

Cremer–Pople coordinates summarize the out-of-plane shape of an ordered five- or six-membered ring. The amplitude describes how far the ring departs from planarity; angles describe its position in conformational space. matmalib provides both continuous coordinates and nearest-reference labels.

1. Analyze a glucose ring

from matmalib import Glycan

glucose = Glycan.from_file("data/glucose.xyz")
coordinates = glucose.cremer_pople()[0]
classified = glucose.identify_ring_puckers()[0]
for name in ("Q", "q2", "q3", "theta", "phi"):
    print(f"{name:>5}: {float(coordinates[name]):.6f}")
print("Label:", classified["label"])
print("Family:", classified["family"])
print("Reference distance (degrees):", round(classified["distance"], 3))
print("CarPpy angle vector:", classified["carppy_theta"])

Expect a 4C1 chair with theta near zero. A real geometry generally has a small nonzero q2, so it can have a finite phi even near a chair. An exact mathematical chair has undefined phi.

Glycan.cremer_pople(degrees=True) returns one coordinate dictionary per monomer. Glycan.identify_ring_puckers(**kwargs) returns one classification dictionary per monomer and accepts the same options as classify_ring_pucker.

Key	Five-membered ring	Six-membered ring	Units
Q	Equal to q2	sqrt(q2**2 + q3**2)	Coordinate length unit
q2	Nonnegative amplitude of the second mode	Nonnegative amplitude of the second mode	Coordinate length unit
q3	Absent	Signed alternating-height amplitude	Coordinate length unit
phi	Pseudorotation phase	Second-mode phase	Degrees by default, in [0, 360)
theta	Absent	Polar angle, atan2(q2, q3)	Degrees by default, in [0, 180]

Set degrees=False for radians: phi lies in [0, 2π) and theta in [0, π]. Length units are never converted. Atom order is part of the definition; moving the starting atom or reversing ring traversal changes the angular convention. See Conventions before comparing results across programs.

2. Build a ring with known coordinates

This geometric experiment inverts the Cremer–Pople height expansion over a regular clockwise polygon. It is a controllable shape model, not an energy-minimized molecule. Change q2, q3, and phase_degrees to see the recovered coordinates.

import numpy as np
import matplotlib.pyplot as plt
from matmalib import cremer_pople

q2 = 0.30
q3 = 0.40
phase_degrees = 60.0

j = np.arange(6)
angle = 2 * np.pi * j / 6
height = np.sqrt(2 / 6) * q2 * np.cos(2 * angle + np.radians(phase_degrees))
height += q3 / np.sqrt(6) * (-1.0) ** j
xyz = np.column_stack((1.4 * np.cos(angle), -1.4 * np.sin(angle), height))
result = cremer_pople(xyz)
print({key: round(float(value), 6) for key, value in result.items()})

fig = plt.figure(figsize=(6, 4))
ax = fig.add_subplot(projection="3d")
closed = np.vstack((xyz, xyz[0]))
ax.plot(*closed.T, "o-", color="#087f8c")
for index, point in enumerate(xyz):
    ax.text(*point, str(index))
ax.set(xlabel="x", ylabel="y", zlabel="height",
       title=f"Q={float(result['Q']):.2f}; theta={float(result['theta']):.1f}°")
ax.set_box_aspect((1, 1, 0.5))
fig.tight_layout()

The starting values give Q=0.5, theta approximately 36.8699°, and phi 60°. Try the following:

• Set q2=0, q3=0.6: a pure chair, theta 0°, phi NaN.
• Set q2=0, q3=-0.6: the opposite pure chair, theta 180°, phi NaN.
• Set q3=0, q2=0.4: theta 90°; changing phase moves around the equator.
• Set both amplitudes to zero: a planar ring, with theta and phi NaN.

A planar, noncollapsed ring is valid. Collinear or coincident atoms can make the mean plane degenerate; those coordinates raise ValueError instead of returning a misleading orientation.

3. Use five-membered rings and vectorized batches

The low-level cremer_pople(coordinates, degrees=True) accepts shape (5, 3) or (6, 3), or any leading batch dimensions: (..., 5, 3) or (..., 6, 3). All rings in one batch must have the same size. Results have the corresponding leading dimensions.

import numpy as np
from matmalib import cremer_pople, classify_ring_pucker

phases = np.array([0.0, 18.0, 90.0, 180.0])
j = np.arange(5)
angle = 2 * np.pi * j / 5
rings = []
for phase in phases:
    height = np.sqrt(2 / 5) * 0.4 * np.cos(2 * angle + np.radians(phase))
    rings.append(np.column_stack((1.4 * np.cos(angle), -1.4 * np.sin(angle), height)))

batch = cremer_pople(np.array(rings))
print("Amplitudes:", np.round(batch["Q"], 6))
print("Recovered phases:", np.round(batch["phi"], 6))
print("Available keys:", list(batch))
for phase, xyz in zip(phases, rings):
    result = classify_ring_pucker(xyz)
    print(f"Phase {phase:5.1f}: {result['label']:>3} ({result['family']})")

There is no theta or q3 for five-membered rings. At phases 0°, 18°, and 180°, expect 0E, 0T1, and E0. The indices 0 through 4 are input ring positions. 0E raises atom 0 relative to the CP normal; 0T1 raises atom 0 and lowers atom 1.

classify_ring_pucker analyzes one ring at a time. For a batch, loop over rings as above or use the trajectory interface.

4. Understand canonical classification

Classification combines CP coordinates with the following metadata:

Key	Meaning
label	Accepted nearest label, "unassigned" when distance rejection applies, or "planar".
nearest_label	Closest candidate even if rejected; None for planar rings.
family	Chair, boat, skew, half-chair, envelope, twist, or planar, as applicable. A rejected match retains its nearest candidate's family.
distance	Reference distance in degrees, even when CP angles are requested in radians. NaN for planar rings.
method	"carppy-three-angle" for nonplanar six-membered rings; "cremer-pople-phase" for nonplanar five-membered rings. Absent for planar rings.
carppy_theta	Three transformed diagonal dihedrals in degrees for nonplanar six-membered rings. This is distinct from the CP polar theta.

Six-membered rings are compared with 38 canonical CarPpy references: two chairs and the reference boat, skew, half-chair, and envelope conformations. Distance is the Euclidean distance between three-angle vectors, not a distance between CP polar angles. Names use positional numbers 1–5 for the five carbons after ring oxygen, and 6 for oxygen. Use oxygen-first cyclic order for this classification. For ketoses, these positions are not chemical carbon numbers.

Five-membered rings are compared with 20 ideal envelope/twist phases. Distance is the shortest wrapped phase difference, in degrees. This distance and the six-membered three-angle distance are different metrics and should not share an assumed universal cutoff.

Parameter	Default	Use
degrees	True	Units for returned CP angles only.
planar_tolerance	1e-6	Rings with Q <= planar_tolerance are labeled planar. Uses coordinate length units; must be finite and nonnegative.
max_distance	None	Accept every nearest nonplanar reference by default. Supply a finite nonnegative cutoff to reject distant matches.

from matmalib import Glycan, classify_ring_pucker, ring_dihedrals

glycan = Glycan.from_file("data/glucose.xyz")
xyz = glycan.xyz[list(glycan.rings[0])]
for cutoff in (None, 5.0, 0.0):
    result = classify_ring_pucker(xyz, max_distance=cutoff, degrees=False)
    print("Cutoff:", cutoff, "label:", result["label"],
          "nearest:", result["nearest_label"],
          "distance:", round(result["distance"], 4), "degrees")
print("CP theta in radians:", float(result["theta"]))
print("Three-angle vector in degrees:", ring_dihedrals(xyz))

This glucose geometry is about 10.57° from its nearest three-angle reference, so both the 5.0 and 0.0 cutoffs reject it. The zero-distance cutoff rejects a real geometry unless it exactly matches a canonical reference. nearest_label lets you inspect the candidate without treating every distorted ring as an exact canonical conformation.

ring_dihedrals(coordinates) is also public. It accepts one six-membered ring or a batch (..., 6, 3) and always returns degrees. Its transformed diagonal definitions use ring-order positions (5, 3, 1, 2), (1, 5, 3, 4), and (3, 1, 5, 0); each signed dihedral is transformed to sign(angle) * (180 - abs(angle)).

Validation and limitations

CP calculations reject invalid shapes, nonfinite coordinates, and degenerate mean planes. They do not detect chemistry, infer atom ordering, or correct periodic boundaries. Whole, unwrapped rings are required. Rigid translations and proper rotations preserve the results; reflecting a geometry is a change of handedness, not an alignment operation.

The bundled glucose_puckers.xyz contains 38 reference conformers. Repository tests check their canonical labels, motion invariance, five-membered phase conventions, planar behavior, and trajectory consistency. These references validate the implementation's convention; a nearest-reference classification does not establish stability, population, or a reaction pathway.

Next, analyze many frames and display their populations, or cluster complete conformers by RMSD.