Compare, cluster, and deduplicate conformers¶

RMSD measures coordinate differences after establishing atom correspondence. matmalib aligns conformers with a proper Kabsch rotation and translation, clusters their pairwise RMSDs, and selects representative states. These operations preserve the input geometries and their order.

1. Separate rigid motion from structural change¶

import numpy as np
from matmalib import rmsd

xyz = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0],
                [0.2, 2.0, 0.0], [0.0, 0.3, 3.0]])
angle = np.radians(65)
rotation = np.array([[np.cos(angle), -np.sin(angle), 0],
                     [np.sin(angle), np.cos(angle), 0], [0, 0, 1]])
moved = xyz @ rotation + [20, -11, 3]
changed = moved.copy()
changed[3] += [0.0, 0.0, 0.7]
print("Moved, aligned:", rmsd(xyz, moved))
print("Moved, unaligned:", round(rmsd(xyz, moved, align=False), 3))
print("Changed, aligned:", round(rmsd(xyz, changed), 3))
print("Mirror image, aligned:", round(rmsd(xyz, xyz * [-1, 1, 1]), 3))

The translated/rotated geometry has aligned RMSD zero, while the changed geometry and mirror image remain different. align=False compares the coordinates directly, including absolute translation and orientation.

rmsd(first, second, atom_indices=None, heavy_only=False, align=True) accepts two State objects or two coordinate arrays of shape (n_atoms, 3). It returns one float in the input length unit, normally angstroms. Neither input is modified.

Atom correspondence is positional. The algorithm does not permute symmetric atoms and does not align a reflection. States must have the same element sequence; arrays have no elements to check. An equal element sequence alone cannot reveal a swap between two carbon atoms. Align and compare the same chemical atoms in the same order.

2. Inspect an ensemble's distance matrix¶

import matplotlib.pyplot as plt
from matmalib import Space

space = Space.from_xyz("data/glucose_puckers.xyz")
distances = space.rmsd_matrix()
print("Frames:", len(space.states))
print("Matrix shape:", distances.shape)
print("First two conformers, heavy-atom RMSD:", round(distances[0, 1], 4))

fig, ax = plt.subplots(figsize=(5, 4))
image = ax.imshow(distances, cmap="viridis", origin="lower")
fig.colorbar(image, ax=ax, label="Aligned heavy-atom RMSD (angstrom)")
ax.set(xlabel="Conformer index", ylabel="Conformer index", title="38 glucose conformers")
fig.tight_layout()

The matrix is symmetric, has a zero diagonal, and has shape (38, 38). Space.from_xyz reads every XYZ frame, assigns string IDs "0", "1", …, and retains the input path. glycans=True creates a Glycan for each frame; ordinary State objects are enough for RMSD and avoid repeated sugar detection on distorted geometries. XYZ comments are not parsed for energies.

rmsd_matrix(states, atom_indices=None, heavy_only=True, align=True) is also available as a standalone function accepting a Space or iterable. A dense N × N matrix requires quadratic memory and pairwise work; a large trajectory may need subsampling before this analysis.

Select the atoms that define similarity¶

The single-pair function defaults to all atoms. The matrix and clustering functions default to heavy atoms. Set the same selection explicitly when comparing their results.

from matmalib import Glycan, Space, rmsd, rmsd_matrix

space = Space.from_xyz("data/glucose_puckers.xyz")
reference = Glycan.from_file("data/glucose.xyz")
first, second = space.states[:2]
print("All atoms:", round(rmsd(first, second), 4))
print("Heavy atoms:", round(rmsd(first, second, heavy_only=True), 4))
print("Ring only:", round(rmsd(first, second, atom_indices=reference.rings[0]), 4))

# Raw arrays do not carry elements, so explicitly disable the heavy-atom filter.
array_distances = rmsd_matrix([first.xyz, second.xyz], heavy_only=False)
print("Array matrix:", array_distances.round(4))

atom_indices is a nonempty sequence of unique, in-range integer atom-list positions. heavy_only=True then removes hydrogen atoms from that selection; it requires states with element information. An empty final selection raises ValueError. Select a ring to compare ring shape or a whole glycan to include linkage and substituent motion; those are different definitions of similarity.

3. Cluster and inspect representative structures¶

from matmalib import Space, Plot

space = Space.from_xyz("data/glucose_puckers.xyz")
threshold = 0.5
clusters = space.cluster(threshold=threshold, method="average", energy=None)
print("Cluster count:", len(clusters.members))
print("Labels in original frame order:", clusters.labels.tolist())
print("Representative frame indices:", clusters.representative_indices)
for number, members in enumerate(clusters.members):
    print(f"Cluster {number}: {members}")

fig, ax = Plot().cluster_dendrogram(clusters, leaf_rotation=90, leaf_font_size=8)
fig.set_size_inches(10, 4)
ax.set_title(f"Average-linkage clusters at {threshold} angstrom")
fig.tight_layout()

At threshold=0.5, the fixture yields 13 average-linkage clusters. With energy=None, each representative is a medoid: the member with the smallest summed RMSD to other members of its cluster. Try 0.3 and 0.8; they produce 24 and 5 clusters respectively for this fixture. Thresholds are geometric choices, not universal chemical cutoffs.

The standalone cluster_conformers(states, ...) accepts a Space or iterable and returns the same ClusterResult as space.cluster(...).

Option	Default	Meaning
`threshold`	`1.0`	Finite nonnegative linkage-distance cutoff, in coordinate length units.
`method`	`"single"`	`"single"`, `"complete"`, `"average"`, or `"weighted"`.
`heavy_only`	`True`	Exclude hydrogen from the fit and RMSD.
`atom_indices`	`None`	Restrict fit and RMSD to selected positions.
`energy`	`"E"`	State attribute used for representative selection/filtering, or `None` to force medoids.
`energy_window`	`None`	Optional maximum energy above the minimum, always supplied in kcal/mol.
`energy_unit`	`"hartree"`	Units of the stored state energies: `"hartree"`, `"kcal/mol"`, or `"kJ/mol"`.

Clustering always uses aligned RMSD. Unlike the separate RMSD functions, it has no align=False option. A representative uses the lowest energy only when every input state has a finite scalar value for the requested energy attribute. Otherwise all clusters use medoids. A partially populated or nonfinite energy list does not silently become a partial energy-based ranking. Equal scores select the first member in the cluster's original-index order.

`ClusterResult` attribute	Meaning
`labels`	One zero-based cluster label per original input; energy-excluded inputs are `-1`.
`members`	Tuple of original-index tuples, one per cluster. Clusters are ordered by their smallest input index.
`representative_indices`	Original input index of each representative, in cluster order.
`representatives`	New `Space` sharing the selected `State` objects.
`distances`	Pairwise RMSD matrix for all inputs, including energy-excluded states.
`linkage`	SciPy linkage tree for included inputs; leaf positions refer to `included_indices`.
`included_indices`	Original indices retained after energy filtering.
`threshold`, `method`	Analysis settings recorded in the result.

The new representative Space shares its states with the original. It does not duplicate their atoms. Mutating a representative state therefore affects the same state in the input Space; make a deep copy if you intend independent edits.

Choose a linkage method¶

Method	Distance between clusters	Practical consequence
`single`	Closest pair of members	Can chain together a series whose endpoints are far apart.
`complete`	Farthest pair of members	A cut at the threshold limits every pair within a cluster.
`average`	Mean pairwise distance across both clusters	Groups by average cross-cluster similarity.
`weighted`	Equal-weight average of the two merged clusters' distances	Weights merged clusters equally, independent of their sizes.

Ward and centroid methods are not accepted for these RMSD distance inputs. Empty ensembles, invalid thresholds, incompatible atom counts/order, and nonfinite coordinates raise errors. A single conformer is valid and produces one cluster with an empty (0, 4) linkage tree.

4. Apply an energy window and select low-energy members¶

This deliberately small artificial example uses different two-atom bond lengths so the grouping is easy to audit. Its assigned energies are illustrative kcal/mol values, not calculation results.

from matmalib import Atom, State, Space

states = []
for length, energy in [(1.0, 0.4), (1.4, 0.0), (4.0, 3.0)]:
    states.append(State(atoms=[Atom([0, 0, 0], "C"), Atom([length, 0, 0], "C")],
                        E=energy))
space = Space(states=states)
all_states = space.cluster(threshold=0.25, energy_unit="kcal/mol")
filtered = space.cluster(threshold=0.25, energy_window=1.0, energy_unit="kcal/mol")
print("Before filtering:", all_states.members)
print("Lowest-energy representatives:", all_states.representative_indices)
print("After filtering labels:", filtered.labels.tolist())
print("Included original indices:", filtered.included_indices.tolist())
print("Filtered representative:", filtered.representative_indices)

Before filtering, states 0 and 1 form a cluster, represented by the lower-energy state 1. State 2 forms its own cluster. The 1 kcal/mol energy window excludes state 2, so the final labels are [0, 0, -1] and the representative is index 1.

An energy window requires complete finite energies on every input. It is always specified in kcal/mol even if stored E values are hartree or kJ/mol. The conversion uses 1 hartree = 627.509474 kcal/mol and 1 kcal = 4.184 kJ. Use energy="F" to select another state attribute. The algorithm calculates the full distance matrix before energy filtering, so filtering does not remove that initial memory cost.

Plot.cluster_dendrogram(result, labels=None, ax=None, **kwargs) labels leaves by original input index by default and draws a dashed cutoff line. Custom labels must cover all original inputs, including those later excluded. Excluded labels are omitted from the displayed tree. Other keywords pass to SciPy's dendrogram, for example orientation="left", leaf_rotation=90, or leaf_font_size=8. The method supports a single remaining state and returns (figure, axes), also stored on the Plot.

5. Remove near duplicates without chaining¶

space.remove_duplicates(threshold=0.1, **kwargs) uses complete linkage and returns only the representative Space. Every pair in each duplicate group is within the threshold. To inspect the groups too, call space.cluster(..., method="complete").

from matmalib import Atom, State, Space

space = Space(states=[
    State(atoms=[Atom([0, 0, 0], "C"), Atom([length, 0, 0], "C")])
    for length in (1.0, 2.0, 3.0)
])
threshold = 0.6
single = space.cluster(threshold, method="single", energy=None)
complete = space.cluster(threshold, method="complete", energy=None)
unique = space.remove_duplicates(threshold, energy=None)
print("Single-linkage groups:", single.members)
print("Complete-linkage groups:", complete.members)
print("Original conformers:", len(space.states))
print("Retained representatives:", len(unique.states))

Adjacent conformers have RMSD 0.5, but the two endpoints have RMSD 1.0. Single linkage joins all three through the middle member. Duplicate removal retains two representatives because it enforces the stricter pairwise limit. The original Space still contains all three states.

Duplicate-removal options match cluster, except method is fixed to complete linkage. Do not supply a second method. Selection depends on the same complete-energy or medoid policy; duplicate removal is not automatically “keep the first frame.”

6. Export representative coordinates locally¶

This writes a multi-frame XYZ file and keeps the original index in each comment for traceability.

from matmalib import Space

space = Space.from_xyz("conformers.xyz")
result = space.cluster(threshold=0.5, method="average")
with open("representatives.xyz", "w") as stream:
    for cluster, index in enumerate(result.representative_indices):
        state = space.states[index]
        stream.write(f"{len(state.atoms)}\nsource frame {index}; cluster {cluster}\n")
        for atom in state.atoms:
            x, y, z = atom.coord
            stream.write(f"{atom.type} {x:.8f} {y:.8f} {z:.8f}\n")

For each retained geometry, identify its sugar residues and draw SNFG, or classify its ring puckers. A cluster of similar full-system coordinates need not have a single exact pucker label, and a shared pucker label need not imply a low full-system RMSD.