Skip to content

molclusters.cluster

cluster

Provides the groups of molecules that MolClusters tracks and analyses.

Classes:
- MolGroup: A group of molecules (residues) of a Universe, with properties such as
  mass, charge, dipole moment, radius of gyration, shape and density, computed on
  the group made whole across periodic boundaries. Nuclei are MolGroups.

- Cluster: A MolGroup that is a tracked molecular cluster: it has an id, a birth
  time and a read-only graph of the connections between its molecules. Clusters
  are created and updated by a `ClusterTracker`; everyone else only reads them.
Constants:
- EA2D: Conversion factor for dipole moments from e·Å (MDAnalysis's unit) to
  Debye.
Dependencies:
- MDAnalysis: For molecular dynamics trajectory and structure analysis.
- NetworkX: For the graph of a cluster's connections.
- NumPy: For numerical computations.

MolGroup

MolGroup(universe: Universe, residues: Iterable[int] | ResidueGroup)

A group of molecules (residues) of a Universe, and its properties.

Geometric properties are computed on the group made whole across periodic boundaries (see whole), and always for the Universe's current frame.

The group describes its residues, so it can't be changed through what it hands out: its arrays are read-only (to change one, work on a copy: array.copy()), and residues is a new group every time. Positions and topology attributes belong to the shared Universe, not the group, and changing them changes every group.

Attributes:

Name Type Description
_uni Universe

The MDAnalysis Universe object associated with the group.

_rg ResidueGroup

The group's residues, over indices no one else holds.

__cache_key tuple[int, bytes] | None

Frame and residue indices that __cache holds values for.

__cache dict[str, object]

Values computed for the current frame and residues (see _frame_cache).

Initialize the MolGroup with a universe and residues.

Parameters:

Name Type Description Default
universe Universe

The MDAnalysis Universe object associated with the group.

required
residues Union[Iterable[int], ResidueGroup]

The residues of the group, either as an iterable with residue IDs or a ResidueGroup. Residue IDs are taken as positions in the topology plus one, so they must be numbered 1 to N (see conntable.check_resids).

required
Source code in src/molclusters/cluster.py
def __init__(
    self, universe: mda.Universe, residues: Iterable[int] | core.groups.ResidueGroup
) -> None:
    """Initialize the MolGroup with a universe and residues.

    Parameters
    ----------
    universe : mda.Universe
        The MDAnalysis Universe object associated with the group.
    residues : Union[Iterable[int], core.groups.ResidueGroup]
        The residues of the group, either as an iterable with residue IDs or a
        ResidueGroup. Residue IDs are taken as positions in the topology plus
        one, so they must be numbered 1 to N (see `conntable.check_resids`).
    """
    self._uni = universe

    if isinstance(residues, core.groups.ResidueGroup):
        ix = residues.ix
    else:
        ix = np.array(residues) - 1
    self._rg = _own_residues(ix, universe)

    self.__cache_key: tuple[int, bytes] | None = None
    self.__cache: dict[str, object] = {}

universe property

universe: Universe

The Universe the group belongs to.

Returns:

Type Description
Universe

The group's Universe.

residues property

residues: ResidueGroup

The group's residues.

Returns:

Type Description
ResidueGroup

The residues of the group, as a new group: changing its indices (ix) in place leaves the group's own residues as they are.

atoms property

atoms: AtomGroup

The atoms of the group's residues.

Returns:

Type Description
AtomGroup

The atoms of the group, in the Universe's current positions (see whole for whole positions).

size property

size: int

The number of molecules in the group.

Returns:

Type Description
int

The number of residues in the group.

resnames property

resnames: ndarray

The residue names of the group.

Returns:

Type Description
ndarray

The residue names in the group (read-only).

resids property

resids: ndarray

The residue IDs of the group.

Returns:

Type Description
ndarray

The residue IDs in the group (read-only).

composition property

composition: Counter[str]

The number of molecules of each residue name in the group.

Returns:

Type Description
Counter[str]

Residue name -> number of molecules with that name.

mass property

mass: float

Calculate the total mass of the group.

Returns:

Type Description
float

The total mass of the group.

center_of_mass property

center_of_mass: ndarray

Calculate the center of mass of the group, made whole.

Unlike MDAnalysis's own center_of_mass, a group split across periodic boundaries gets the center of its whole shape, not a point between its pieces.

Returns:

Type Description
ndarray

The center of mass, wrapped into the primary unit cell (read-only).

sphericity property

sphericity: float

Calculate how spherical the group is, from its gyration tensor.

This is one minus the asphericity (Dima & Thirumalai), not Wadell's surface-area sphericity: 1 when the three principal moments of the gyration tensor are equal (a sphere, but also e.g. a cube), 0 for a rod (all mass on a line), and 0.75 for a flat disk or ring.

Returns:

Type Description
float

1 - asphericity, between 0 (rod) and 1 (spherical).

dipole_moment property

dipole_moment: float

Calculate the dipole moment of the group, about its center of mass.

For a neutral group the dipole is the same about any point. For a group with a net charge q (ions, charged nuclei or clusters) it isn't: moving the reference point by d changes the dipole by -q d, so its value is the dipole about the center of mass specifically, and only comparable between groups with that convention in mind.

Returns:

Type Description
float

The dipole moment of the group in Debye (D).

dipole property

dipole: ndarray

Calculate the dipole vector of the group, about its center of mass.

See dipole_moment for what that reference point means for charged groups.

Returns:

Type Description
ndarray

The dipole vector of the group in Debye (D) (read-only).

shape_parameter property

shape_parameter: float

Calculate the shape parameter of the group.

Returns:

Type Description
float

The shape parameter of the group, a measure of its geometric anisotropy.

bsphere property

bsphere: tuple[float, ndarray]

Calculate the bounding sphere of the group.

Returns:

Type Description
tuple[float, ndarray]

The radius and center of the bounding sphere, the center (read-only) in the group's whole positions (see whole).

radius_of_gyration property

radius_of_gyration: float

Calculate the radius of gyration of the group.

Returns:

Type Description
float

The radius of gyration of the group.

radius_buffer property

radius_buffer: float

Half the mean van der Waals radius of the group's atoms, by element.

The radius of gyration only sees the atoms' centers, so a sphere built from it ends at the outer atoms' centers and leaves their size out. Adding half, not all, of the atoms' van der Waals radius to it roughly accounts for that, since bonded atoms' van der Waals spheres overlap: on methanol-malic acid mixtures, the full radius puts a lone malic acid at 0.70 g/cm^3 and half at 1.21 (1.61 for the solid), and barely changes large clusters.

Returns:

Type Description
float

The buffer added to the equivalent sphere's radius, in angstroms.

Raises:

Type Description
ValueError

If the atoms have no elements, or an element has no van der Waals radius in MDAnalysis' table.

radius property

radius: float

The radius of the group's equivalent sphere.

That is the uniform sphere with the group's radius of gyration, which has a radius of sqrt(5/3) times it (a uniform sphere of radius R has a radius of gyration of sqrt(3/5) R), grown by radius_buffer to roughly account for the size of the atoms.

Returns:

Type Description
float

The group's equivalent sphere radius, in angstroms.

diameter property

diameter: float

Calculate the diameter of the group's equivalent sphere (see radius).

Returns:

Type Description
float

Twice the group's radius, in angstroms.

volume property

volume: float

Calculate the volume of the group's equivalent sphere (see radius).

Returns:

Type Description
float

Volume in cubic angstroms.

density property

density: float

Calculate the density of the group.

Returns:

Type Description
float

The density of the group in g/cm^3.

charge property

charge: float

Calculate the total charge of the group.

Returns:

Type Description
float

The total charge of the group.

__whole

__whole() -> tuple[ndarray, ndarray]

Get the group's whole positions and shift, computed once per frame.

Returns:

Type Description
tuple[ndarray, ndarray]

See __compute_whole_positions.

Source code in src/molclusters/cluster.py
def __whole(self) -> tuple[np.ndarray, np.ndarray]:
    """Get the group's whole positions and shift, computed once per frame.

    Returns
    -------
    tuple[np.ndarray, np.ndarray]
        See `__compute_whole_positions`.
    """
    cache = self._frame_cache()
    if "whole" not in cache:
        cache["whole"] = self.__compute_whole_positions()
    return cache["whole"]

__compute_whole_positions

__compute_whole_positions() -> tuple[ndarray, ndarray]

Compute the group's atom positions with the group made whole.

The group is translated so that its first residue sits at the box center, and each residue is made whole and wrapped into the box by its center of geometry, which gathers a group spanning less than half the box around the center. A larger group would be torn apart that way, so the residues are then placed one after another along a spanning tree of the group (see _placement_order), each at the periodic image nearest its tree neighbour. The Universe's positions are left exactly as they were.

Returns:

Type Description
tuple[ndarray, ndarray]

Whole positions, one row per atom of self._rg.atoms, and the translation applied to the group (the wrapping only adds whole box vectors on top of it).

Source code in src/molclusters/cluster.py
def __compute_whole_positions(self) -> tuple[np.ndarray, np.ndarray]:
    """Compute the group's atom positions with the group made whole.

    The group is translated so that its first residue sits at the box center,
    and each residue is made whole and wrapped into the box by its center of
    geometry, which gathers a group spanning less than half the box around the
    center. A larger group would be torn apart that way, so the residues are
    then placed one after another along a spanning tree of the group (see
    `_placement_order`), each at the periodic image nearest its tree neighbour.
    The Universe's positions are left exactly as they were.

    Returns
    -------
    tuple[np.ndarray, np.ndarray]
        Whole positions, one row per atom of ``self._rg.atoms``, and the
        translation applied to the group (the wrapping only adds whole box
        vectors on top of it).
    """
    atoms = self._rg.atoms
    original = atoms.positions  # a copy
    try:
        boxcenter = np.sum(self._uni.trajectory.ts.triclinic_dimensions, axis=0) / 2
        ref_mol_cm = self._rg[:1].center_of_mass(unwrap=True)
        shift = boxcenter - ref_mol_cm
        atoms.positions += shift
        atoms.unwrap(compound="residues", reference="cog", inplace=True)
        positions = atoms.positions
    finally:
        atoms.positions = original

    if len(self._rg) > 1:
        self.__place_residues(atoms, positions)
    return positions, shift

__place_residues

__place_residues(atoms: AtomGroup, positions: ndarray) -> None

Move each residue to the periodic image nearest its placed neighbour.

Changes positions in place, by whole box vectors per residue: none at all for a group spanning less than half the box.

Parameters:

Name Type Description Default
atoms AtomGroup

The group's atoms.

required
positions ndarray

Their positions, each residue whole and wrapped into the box.

required
Source code in src/molclusters/cluster.py
def __place_residues(
    self, atoms: core.groups.AtomGroup, positions: np.ndarray
) -> None:
    """Move each residue to the periodic image nearest its placed neighbour.

    Changes `positions` in place, by whole box vectors per residue: none at all
    for a group spanning less than half the box.

    Parameters
    ----------
    atoms : core.groups.AtomGroup
        The group's atoms.
    positions : np.ndarray
        Their positions, each residue whole and wrapped into the box.
    """
    box = self._uni.dimensions
    if box is None or not np.any(box[:3]):
        return
    # each atom's residue, as its row in self._rg
    by_ix = np.argsort(self._rg.ix)
    rows = by_ix[np.searchsorted(self._rg.ix[by_ix], atoms.resindices)]
    counts = np.bincount(rows, minlength=len(self._rg))
    centers = (
        np.stack(
            [np.bincount(rows, positions[:, k], len(self._rg)) for k in range(3)],
            axis=1,
        )
        / counts[:, None]
    )

    # in a rectangular box, residues spanning less than half of it along every
    # axis are all each other's nearest images already: nothing to place, and
    # nothing wrapping around the box (a connection torn by the wrapping into
    # the box would stretch the span past half of it)
    # (plain comparisons: np.allclose and np.ptp cost more than the rest here)
    rectangular = (np.abs(box[3:] - 90.0) < 1e-3).all()
    span = centers.max(axis=0) - centers.min(axis=0)
    if rectangular and (span < box[:3] / 2).all():
        return

    order, parent = self._placement_order(centers, box)
    steps = minimize_vectors(
        (centers[order[1:]] - centers[parent[order[1:]]]).astype(np.float32), box
    )
    placed = centers.copy()
    for row, step in zip(order[1:], steps, strict=True):
        placed[row] = placed[parent[row]] + step

    self._check_placement(placed, centers, box)
    positions += (placed - centers)[rows].astype(positions.dtype)

whole

whole() -> Iterator[AtomGroup]

Temporarily place the group's atoms in whole, centered positions.

Inside the with block the group's atoms are whole across periodic boundaries, so geometric properties (and writers) see the real shape. On exit, even on error, the atoms go back to their original positions.

Positions are never changed for good because they belong to the shared Universe: other groups over overlapping residues (e.g. a nucleus inside a cluster) would otherwise move atoms under each other and corrupt each other's geometry.

The whole positions are cached per frame and residue set (see _frame_cache), so the unwrap runs once per frame no matter how many properties are read.

Yields:

Type Description
AtomGroup

The group's atoms, in whole positions.

Source code in src/molclusters/cluster.py
@contextmanager
def whole(self) -> Iterator[core.groups.AtomGroup]:
    """Temporarily place the group's atoms in whole, centered positions.

    Inside the ``with`` block the group's atoms are whole across periodic
    boundaries, so geometric properties (and writers) see the real shape. On
    exit, even on error, the atoms go back to their original positions.

    Positions are never changed for good because they belong to the shared
    Universe: other groups over overlapping residues (e.g. a nucleus inside a
    cluster) would otherwise move atoms under each other and corrupt each
    other's geometry.

    The whole positions are cached per frame and residue set (see
    `_frame_cache`), so the unwrap runs once per frame no matter how many
    properties are read.

    Yields
    ------
    core.groups.AtomGroup
        The group's atoms, in whole positions.
    """
    positions, _ = self.__whole()
    atoms = self._rg.atoms
    original = atoms.positions  # a copy
    atoms.positions = positions
    try:
        yield atoms
    finally:
        atoms.positions = original

__add__

__add__(other: ResidueGroup | Self) -> MolGroup

Combine this group with a ResidueGroup or another MolGroup.

Parameters:

Name Type Description Default
other Union[ResidueGroup, MolGroup]

The ResidueGroup or MolGroup to combine with.

required

Returns:

Type Description
MolGroup

A new MolGroup with the residues of both.

Source code in src/molclusters/cluster.py
def __add__(self, other: core.groups.ResidueGroup | Self) -> "MolGroup":
    """Combine this group with a ResidueGroup or another MolGroup.

    Parameters
    ----------
    other : Union[core.groups.ResidueGroup, MolGroup]
        The ResidueGroup or MolGroup to combine with.

    Returns
    -------
    MolGroup
        A new MolGroup with the residues of both.
    """
    if isinstance(other, core.groups.ResidueGroup):
        return MolGroup(self._uni, self._rg + other)
    elif isinstance(other, MolGroup):
        return MolGroup(self._uni, self._rg + other.residues)

    return NotImplemented

__len__

__len__() -> int

Get the number of molecules in the group.

Returns:

Type Description
int

The number of residues in the group.

Source code in src/molclusters/cluster.py
def __len__(self) -> int:
    """Get the number of molecules in the group.

    Returns
    -------
    int
        The number of residues in the group.
    """
    return len(self._rg)

__iter__

__iter__() -> Iterator[int]

Iterate over the residue IDs of the group's molecules.

Returns:

Type Description
Iterator[int]

An iterator over the residue IDs of the group's molecules.

Source code in src/molclusters/cluster.py
def __iter__(self) -> Iterator[int]:
    """Iterate over the residue IDs of the group's molecules.

    Returns
    -------
    Iterator[int]
        An iterator over the residue IDs of the group's molecules.
    """
    return iter(self._rg.resids.tolist())

__contains__

__contains__(item: int) -> bool

Check if a molecule is in the group.

Parameters:

Name Type Description Default
item int

The residue ID of the molecule to check.

required

Returns:

Type Description
bool

True if the molecule is in the group, False otherwise.

Source code in src/molclusters/cluster.py
def __contains__(self, item: int) -> bool:
    """Check if a molecule is in the group.

    Parameters
    ----------
    item : int
        The residue ID of the molecule to check.

    Returns
    -------
    bool
        True if the molecule is in the group, False otherwise.
    """
    return bool(np.isin(item, self._rg.resids))

Cluster

Cluster(universe: Universe, subconntab: _SubConnTable, *, cluster_id: int)

Bases: MolGroup

A tracked molecular cluster.

A MolGroup with an id, a birth time and the graph of the connections between its molecules. Clusters are created and updated by a ClusterTracker and are otherwise read-only: the graph is frozen, attributes included, and a new one replaces it every frame; its arrays are read-only (see MolGroup).

A Cluster is the same object for as long as it lives, and its properties are those of the Universe's current frame, so it describes a given frame only while the Universe is at that frame.

Clusters compare and hash by identity, which is the tracked cluster: the tracker keeps one object per cluster for its whole life. To keep data per cluster, key it by id rather than by the Cluster, so a dead cluster (and its cached whole positions) isn't kept in memory.

Attributes:

Name Type Description
_graph Graph

The frozen graph of the cluster's molecules and their connections.

_id int

The cluster's id.

_birth_time float

The time of the frame the cluster was created at.

_wrap_reported bool

Whether the cluster was already reported wrapping around the box.

Initialize a Cluster from a connected group of the current frame.

Parameters:

Name Type Description Default
universe Universe

The MDAnalysis Universe object associated with the cluster.

required
subconntab _SubConnTable

The connected group of molecules the cluster is made of.

required
cluster_id int

The cluster's id, as handed out by whoever tracks it (e.g. a ClusterTracker).

required
Source code in src/molclusters/cluster.py
def __init__(
    self,
    universe: mda.Universe,
    subconntab: ConnectionTable._SubConnTable,
    *,
    cluster_id: int,
) -> None:
    """Initialize a Cluster from a connected group of the current frame.

    Parameters
    ----------
    universe : mda.Universe
        The MDAnalysis Universe object associated with the cluster.
    subconntab : ConnTable._SubConnTable
        The connected group of molecules the cluster is made of.
    cluster_id : int
        The cluster's id, as handed out by whoever tracks it (e.g. a
        `ClusterTracker`).
    """
    self._birth_time: float = universe.coord.time
    self._graph: nx.Graph = _frozen_graph(subconntab.graph)
    self._id = cluster_id
    self._wrap_reported = False

    super().__init__(universe, self._graph)

id property

id: int

The cluster's id, unique among the clusters of its tracker.

Returns:

Type Description
int

The id of the cluster.

graph property

graph: Graph

The molecules of the cluster and their connections, as a frozen graph.

Nodes are residue IDs, and each edge joins two connected molecules. A "cm" rule's edge has the distance between their centers of mass; an "hb" rule's has the donor-acceptor distance and D-H-A angle of the shortest of their H-bonds, and n_hbonds, how many they share.

Returns:

Type Description
Graph

The graph of the cluster, which can't be modified, nor can the attributes of its edges, nodes or itself (nx.Graph(graph) gives a modifiable copy).

birth_time property

birth_time: float

The time of the frame the cluster was created at.

Returns:

Type Description
float

The birth time of the cluster.

age property

age: float

The time since the cluster was created.

Returns:

Type Description
float

The age of the cluster at the Universe's current frame.

neighbors

neighbors(ref: int, *, level: int | None = None) -> set[int]

Get the neighbors of a molecule in the cluster.

Parameters:

Name Type Description Default
ref int

The reference molecule.

required
level int | None

How many connections away to look (1 when None).

None

Returns:

Type Description
set[int]

The molecules at most level connections away from ref.

Source code in src/molclusters/cluster.py
def neighbors(self, ref: int, *, level: int | None = None) -> set[int]:
    """Get the neighbors of a molecule in the cluster.

    Parameters
    ----------
    ref : int
        The reference molecule.
    level : int | None
        How many connections away to look (1 when None).

    Returns
    -------
    set[int]
        The molecules at most `level` connections away from `ref`.
    """
    if level:
        return set(
            nx.single_target_shortest_path(self._graph, ref, cutoff=level)
        ) - {ref}

    return set(self._graph.neighbors(ref))

distance

distance(moli: int, molj: int) -> float

Get the distance between two connected molecules of the cluster.

Parameters:

Name Type Description Default
moli int

The first molecule.

required
molj int

The second molecule.

required

Returns:

Type Description
float

The distance of the connection between the two molecules.

Raises:

Type Description
ValueError

If moli or molj is not in the cluster, or they aren't connected.

Source code in src/molclusters/cluster.py
def distance(self, moli: int, molj: int) -> float:
    """Get the distance between two connected molecules of the cluster.

    Parameters
    ----------
    moli : int
        The first molecule.
    molj : int
        The second molecule.

    Returns
    -------
    float
        The distance of the connection between the two molecules.

    Raises
    ------
    ValueError
        If ``moli`` or ``molj`` is not in the cluster, or they aren't connected.
    """
    for mol in (moli, molj):
        if mol not in self:
            raise ValueError(f"mol {mol} not in the cluster")
    if not self._graph.has_edge(moli, molj):
        raise ValueError(f"mols {moli} and {molj} are not connected")

    return self._graph[moli][molj]["distance"]

__contains__

__contains__(item: int) -> bool

Check if a molecule is in the cluster.

Parameters:

Name Type Description Default
item int

The residue ID of the molecule to check.

required

Returns:

Type Description
bool

True if the molecule is in the cluster, False otherwise.

Source code in src/molclusters/cluster.py
def __contains__(self, item: int) -> bool:
    """Check if a molecule is in the cluster.

    Parameters
    ----------
    item : int
        The residue ID of the molecule to check.

    Returns
    -------
    bool
        True if the molecule is in the cluster, False otherwise.
    """
    return item in self._graph

__iter__

__iter__() -> Iterator[int]

Iterate over the residue IDs of the molecules in the cluster.

Returns:

Type Description
Iterator[int]

An iterator over the residue IDs of the molecules in the cluster.

Source code in src/molclusters/cluster.py
def __iter__(self) -> Iterator[int]:
    """Iterate over the residue IDs of the molecules in the cluster.

    Returns
    -------
    Iterator[int]
        An iterator over the residue IDs of the molecules in the cluster.
    """
    return iter(self._graph)

__len__

__len__() -> int

Get the number of molecules in the cluster.

Returns:

Type Description
int

The number of molecules in the cluster.

Source code in src/molclusters/cluster.py
def __len__(self) -> int:
    """Get the number of molecules in the cluster.

    Returns
    -------
    int
        The number of molecules in the cluster.
    """
    return len(self._graph)

__str__

__str__() -> str

Get a string representation of the cluster.

Returns:

Type Description
str

The cluster's connections, with their attributes.

Source code in src/molclusters/cluster.py
def __str__(self) -> str:
    """Get a string representation of the cluster.

    Returns
    -------
    str
        The cluster's connections, with their attributes.
    """
    return self._graph.edges.data().__str__()