molclusters.conntable

conntable

Defines the ConnectionTable class, which manages molecular connectivity tables for analyzing molecular clusters. It uses NetworkX for graph-based operations and MDAnalysis for molecular dynamics trajectory analysis.

Classes:

- ConnectionTable: Represents a connectivity table for molecular clusters, allowing operations
  such as retrieving connections, subgraphs, and constructing connectivity graphs.

Functions:

- check_resids: Checks that the residues are numbered 1 to N, as MolClusters needs.

Dependencies:

- MDAnalysis: For molecular dynamics trajectory and structure analysis.
- NetworkX: For graph-based operations on molecular connectivity.
- NumPy: For numerical computations.
- HydrogenBondAnalysis: For analyzing hydrogen bonds between molecules.

ConnectionTable

ConnectionTable(universe: Universe, cluster_args: SymmetricDict[str, Rule], selections: dict[str, AtomGroup], backend: DistanceBackend = 'serial')

Represents a connectivity table for molecular clusters.

This class manages molecular connectivity data using NetworkX graphs and provides methods for retrieving connections, constructing subgraphs, and updating connectivity tables.

Attributes:

Name	Type	Description
uni	Universe	The MDAnalysis Universe object associated with the molecular system.
clst_args	SymmetricDict[str, Rule]	The clustering arguments specifying connectivity rules.
sels	dict[str, AtomGroup]	The atom groups for each residue type.
conntab	MultiGraph	The connectivity graph representing molecular connections.
cms	dict[str, ndarray]	The center of mass for each residue type.
hbs	SymmetricDict[str, HydrogenBondAnalysis]	The hydrogen bond analysis objects for residue pairs.
backend	DistanceBackend	The MDAnalysis acceleration backend used for "cm" rule distance calculations. MDAnalysis silently ignores this for its "nsgrid" method, which it auto-selects whenever the cutoff is much smaller than the box (the common case for realistic systems), so switching this doesn't always change timing. Even when it is used, "OpenMP" is a no-op on builds where MDAnalysis's c_distances_openmp extension was compiled without OpenMP actually enabled (see _warn_if_openmp_unavailable). "hb" rules are unaffected, since HydrogenBondAnalysis doesn't expose a backend option.
rule_connections	Counter[tuple[str, str]]	The connections each rule (keyed as in clst_args) has found, summed over every frame the table was built for; for "hb" rules, the H-bonds (a pair of molecules can share several).

Initialize the ConnectionTable.

Parameters:

Name	Type	Description	Default
universe	Universe	The MDAnalysis Universe object associated with the molecular system.	required
cluster_args	dict[str, dict[str, tuple[str, float]]]	The clustering arguments specifying connectivity rules.	required
selections	dict[str, AtomGroup]	The atom groups for each residue type.	required
backend	DistanceBackend	The MDAnalysis acceleration backend for "cm" rule distance calculations (see the class docstring's backend attribute for its caveats).	'serial'

Source code in src/molclusters/conntable.py

def __init__(
    self,
    universe: mda.Universe,
    cluster_args: SymmetricDict[str, Rule],
    selections: dict[str, core.groups.AtomGroup],
    backend: DistanceBackend = "serial",
) -> None:
    """Initialize the ConnectionTable.

    Parameters
    ----------
    universe : mda.Universe
        The MDAnalysis Universe object associated with the molecular system.
    cluster_args : dict[str, dict[str, tuple[str, float]]]
        The clustering arguments specifying connectivity rules.
    selections : dict[str, core.groups.AtomGroup]
        The atom groups for each residue type.
    backend : DistanceBackend
        The MDAnalysis acceleration backend for "cm" rule distance calculations
        (see the class docstring's `backend` attribute for its caveats).
    """
    check_resids(universe)
    self.uni = universe
    self.clst_args = cluster_args
    self.sels = selections
    self.backend = backend
    _warn_if_openmp_unavailable(backend)

    self.cms: dict[str, np.ndarray] = {}
    self.rule_connections: Counter[tuple[str, str]] = Counter()
    self.hbs: SymmetricDict[str, HydrogenBondAnalysis] = SymmetricDict()
    self.__start_hbonds()
    self.update()

__get_mass_centers

__get_mass_centers() -> None

Calculate the whole center of mass of each molecule, per residue type.

Source code in src/molclusters/conntable.py

def __get_mass_centers(self) -> None:
    """Calculate the whole center of mass of each molecule, per residue type."""
    for res in self.clst_args.all_keys():
        self.cms[res] = _whole_centers_of_mass(self.sels[res])

__start_hbonds

__start_hbonds() -> None

Initialize hydrogen bond analysis for residue pairs.

Source code in src/molclusters/conntable.py

def __start_hbonds(self) -> None:
    """Initialize hydrogen bond analysis for residue pairs."""
    for resi, resj in self.clst_args:
        rule = self.clst_args[resi, resj]
        if not isinstance(rule, HBRule):
            continue

        if (resi, resj) not in self.hbs:
            # Left unset, MDAnalysis would guess hydrogens and acceptors over
            # the whole system and search all of them every frame, only to
            # discard the pairs outside `between` afterwards.
            pair_sel = f"resname {resi} {resj}"
            hydrogens_sel, acceptors_sel = _hb_selections(rule)
            # TODO: activate supported backend
            hb = HydrogenBondAnalysis(
                self.uni,
                between=[f"resname {resi}", f"resname {resj}"],
                hydrogens_sel=f"({pair_sel}) and {hydrogens_sel}",
                acceptors_sel=f"({pair_sel}) and {acceptors_sel}",
                d_a_cutoff=rule.dist,
                d_h_a_angle_cutoff=rule.ang,
                update_selections=False,
            )
            _check_hb_private_api(hb)

            hb._prepare()

            self.hbs[resi, resj] = hb

__cm_connections

__cm_connections(resi: str, resj: str, box: ndarray) -> tuple[ndarray, dict[str, ndarray]]

Find the molecules connected by a center-of-mass ("cm") rule.

Parameters:

Name	Type	Description	Default
resi	str	The rule's first residue name.	required
resj	str	The rule's second residue name.	required
box	ndarray	The simulation box dimensions, used for periodic boundary conditions.	required

Returns:

Type	Description
tuple[ndarray, dict[str, ndarray]]	An (n, 2) array of connected resids, and each connection's "distance".

Source code in src/molclusters/conntable.py

def __cm_connections(
    self, resi: str, resj: str, box: np.ndarray
) -> tuple[np.ndarray, dict[str, np.ndarray]]:
    """Find the molecules connected by a center-of-mass ("cm") rule.

    Parameters
    ----------
    resi : str
        The rule's first residue name.
    resj : str
        The rule's second residue name.
    box : np.ndarray
        The simulation box dimensions, used for periodic boundary conditions.

    Returns
    -------
    tuple[np.ndarray, dict[str, np.ndarray]]
        An ``(n, 2)`` array of connected resids, and each connection's
        ``"distance"``.
    """
    cutoff = self.clst_args[resi, resj].dist
    if resi != resj:
        pairs, distances = mda.lib.distances.capped_distance(
            self.cms[resi], self.cms[resj], cutoff, box=box, backend=self.backend
        )
    else:
        pairs, distances = mda.lib.distances.self_capped_distance(
            self.cms[resi], cutoff, box=box, backend=self.backend
        )

    # pairs index each selection's residues, not the resids they carry
    connections = np.column_stack(
        (
            self.sels[resi].residues.resids[pairs[:, 0]],
            self.sels[resj].residues.resids[pairs[:, 1]],
        )
    )
    return connections, {"distance": distances}

__hb_connections

__hb_connections(resi: str, resj: str) -> tuple[ndarray, dict[str, ndarray]]

Find the molecules connected by a hydrogen bond ("hb") rule.

Parameters:

Name	Type	Description	Default
resi	str	The rule's first residue name.	required
resj	str	The rule's second residue name.	required

Returns:

Type

Description

tuple[ndarray, dict[str, ndarray]]

An (n, 2) array of connected molecule pairs (resids, as the first of their H-bonds joins them: hydrogen's molecule first), one row per pair however many H-bonds join it, and for each pair: "distance" and "angle", the donor-acceptor distance and D-H-A angle of its shortest (strongest) H-bond, and "n_hbonds", the number of H-bonds between the two, either way. Only H-bonds between two molecules count.

Source code in src/molclusters/conntable.py

def __hb_connections(
    self, resi: str, resj: str
) -> tuple[np.ndarray, dict[str, np.ndarray]]:
    """Find the molecules connected by a hydrogen bond ("hb") rule.

    Parameters
    ----------
    resi : str
        The rule's first residue name.
    resj : str
        The rule's second residue name.

    Returns
    -------
    tuple[np.ndarray, dict[str, np.ndarray]]
        An ``(n, 2)`` array of connected molecule pairs (resids, as the first
        of their H-bonds joins them: hydrogen's molecule first), one row per
        pair however many H-bonds join it, and for each
        pair: ``"distance"`` and ``"angle"``, the donor-acceptor distance and
        D-H-A angle of its shortest (strongest) H-bond, and ``"n_hbonds"``, the
        number of H-bonds between the two, either way. Only H-bonds between
        two molecules count.
    """
    # TODO: implement own HB analysis as HydrogenBondAnalysis from mda repeats
    # distance and angle calculations. Until then, this depends on the private
    # API guarded by _check_hb_private_api (see its docstring for why).
    hb = self.hbs[resi, resj]
    hb._ts = self.uni.trajectory.ts

    # a frame without hydrogen bonds is normal here, not worth a warning
    with warnings.catch_warnings():
        warnings.filterwarnings("ignore", message="No hydrogen bonds were found")
        hb._single_frame()

    _, _, hydrogens, acceptors, distances, angles = hb.results.hbonds
    hb._prepare()

    atoms = self.uni.atoms
    connections = np.column_stack(
        (
            atoms[np.asarray(hydrogens, dtype=np.intp)].resids,
            atoms[np.asarray(acceptors, dtype=np.intp)].resids,
        )
    )
    # HydrogenBondAnalysis also finds H-bonds within a molecule, which connect
    # it to nothing else
    between = connections[:, 0] != connections[:, 1]
    connections = connections[between]
    distances = np.asarray(distances)[between]
    angles = np.asarray(angles)[between]
    if len(connections) == 0:
        return connections, {
            "distance": distances,
            "angle": angles,
            "n_hbonds": np.zeros(0, dtype=int),
        }

    # a pair of molecules can share several H-bonds (e.g. a carboxylic acid
    # dimer), but the graph holds one connection per pair: keep the shortest
    # H-bond's geometry and count them all. Each pair keeps the place and
    # orientation of its first H-bond, as the graph got them when it kept every
    # H-bond: the order molecules enter the graph breaks ties between clusters
    # born together (see ClusterTracker)
    _, first, pair = np.unique(
        np.sort(connections, axis=1), axis=0, return_index=True, return_inverse=True
    )
    pair = pair.ravel()
    n_hbonds = np.bincount(pair)
    by_distance = np.lexsort((distances, pair))  # by pair, the shortest first
    shortest = by_distance[np.cumsum(n_hbonds) - n_hbonds]
    in_order = np.argsort(first)  # the pairs, as first found
    return connections[first[in_order]], {
        "distance": distances[shortest[in_order]],
        "angle": angles[shortest[in_order]],
        "n_hbonds": n_hbonds[in_order],
    }

__construct_table

__construct_table() -> None

Construct the connectivity table as a graph.

Source code in src/molclusters/conntable.py

def __construct_table(self) -> None:
    """Construct the connectivity table as a graph."""
    self.conntab = nx.Graph()

    for resi, resj in self.clst_args:
        if self.clst_args[resi, resj].type == "cm":
            pairs, attribs = self.__cm_connections(resi, resj, self.uni.dimensions)
        else:
            pairs, attribs = self.__hb_connections(resi, resj)
        # for "hb" rules, every H-bond rather than every connected pair
        self.rule_connections[resi, resj] += int(
            attribs["n_hbonds"].sum() if "n_hbonds" in attribs else len(pairs)
        )

        # a single bulk insert of plain Python values is much faster than an
        # add_edge() call per connection
        columns = (*pairs.T.tolist(), *(a.tolist() for a in attribs.values()))
        self.conntab.add_edges_from(
            (ri, rj, dict(zip(attribs, vals, strict=True)))
            for ri, rj, *vals in zip(*columns, strict=True)
        )

update

update() -> None

Update the connectivity table.

Source code in src/molclusters/conntable.py

def update(self) -> None:
    """Update the connectivity table."""
    self.__get_mass_centers()
    self.__construct_table()

__getitem__

__getitem__(key: tuple[int, int]) -> int

__getitem__(key: int) -> list[int]

__getitem__(key: tuple[int, int] | int) -> float | list[int]

Get the attributes of a connection or molecule in the connectivity table.

Parameters:

Name	Type	Description	Default
key	(tuple[int, int], int)	The molecule or connection to retrieve.	required

Returns:

Type	Description
(float, list[int])	The attributes of the connection or molecule.

Raises:

Type	Description
KeyError	If the key is invalid or not found in the connectivity table.

Source code in src/molclusters/conntable.py

def __getitem__(self, key: tuple[int, int] | int) -> float | list[int]:
    """Get the attributes of a connection or molecule in the connectivity table.

    Parameters
    ----------
    key : tuple[int, int], int
        The molecule or connection to retrieve.

    Returns
    -------
    float, list[int]
        The attributes of the connection or molecule.

    Raises
    ------
    KeyError
        If the key is invalid or not found in the connectivity table.
    """
    if isinstance(key, tuple):
        if len(key) > 2:
            raise KeyError(
                "ConnTable accepts only one or two parameters to access data!"
            )

        if key[0] not in self:
            raise KeyError(f"{key[0]} not in ConnTable")

        if key[1] not in self:
            raise KeyError(f"{key[1]} not in ConnTable")

        return self.conntab[key[0]][key[1]]["distance"]

    if key not in self:
        raise KeyError(f"{key} not in ConnTable")

    return self.mols_connected_to(key)

__contains__

__contains__(item: int) -> bool

Check if a molecule is in the connectivity table.

Parameters:

Name	Type	Description	Default
item	int	The molecule to check.	required

Returns:

Type	Description
bool	True if the molecule is in the connectivity table, False otherwise.

Source code in src/molclusters/conntable.py

def __contains__(self, item: int) -> bool:
    """Check if a molecule is in the connectivity table.

    Parameters
    ----------
    item : int
        The molecule to check.

    Returns
    -------
    bool
        True if the molecule is in the connectivity table, False otherwise.
    """
    return item in self.conntab

__iter__

__iter__() -> Iterator[int]

Iterate over the molecules in the connectivity table.

Returns:

Type	Description
Iterator[int]	An iterator over the molecules in the connectivity table.

Source code in src/molclusters/conntable.py

def __iter__(self) -> Iterator[int]:
    """Iterate over the molecules in the connectivity table.

    Returns
    -------
    Iterator[int]
        An iterator over the molecules in the connectivity table.
    """
    return iter(self.conntab)

connections_from

connections_from(mol: int) -> list[tuple[int, int]]

Get the connections from a molecule.

Parameters:

Name	Type	Description	Default
mol	int	The molecule to retrieve connections from.	required

Returns:

Type	Description
list[tuple[int, int]]	A list of connections from the molecule.

Source code in src/molclusters/conntable.py

def connections_from(self, mol: int) -> list[tuple[int, int]]:
    """Get the connections from a molecule.

    Parameters
    ----------
    mol : int
        The molecule to retrieve connections from.

    Returns
    -------
    list[tuple[int, int]]
        A list of connections from the molecule.
    """
    return list(self.conntab.edges(mol))

connection_tree_from

connection_tree_from(mol: int) -> list[tuple[int, int]]

Get the connection tree from a molecule.

Parameters:

Name	Type	Description	Default
mol	int	The molecule to retrieve the connection tree from.	required

Returns:

Type	Description
list[tuple[int, int]]	A list of connections in the tree.

Source code in src/molclusters/conntable.py

def connection_tree_from(self, mol: int) -> list[tuple[int, int]]:
    """Get the connection tree from a molecule.

    Parameters
    ----------
    mol : int
        The molecule to retrieve the connection tree from.

    Returns
    -------
    list[tuple[int, int]]
        A list of connections in the tree.
    """
    g = self.__graph_from_mol(mol)
    return list(g.edges)

mols_connected_to

mols_connected_to(mol: int) -> list[int]

Get the molecules connected to a given molecule.

Parameters:

Name	Type	Description	Default
mol	int	The molecule to retrieve connected molecules for.	required

Returns:

Type	Description
list[int]	A list of connected molecules.

Source code in src/molclusters/conntable.py

def mols_connected_to(self, mol: int) -> list[int]:
    """Get the molecules connected to a given molecule.

    Parameters
    ----------
    mol : int
        The molecule to retrieve connected molecules for.

    Returns
    -------
    list[int]
        A list of connected molecules.
    """
    return list(self.conntab[mol])

mols_connected_tree_to

mols_connected_tree_to(mol: int) -> list[int]

Get the molecules in the connection tree of a given molecule.

Parameters:

Name	Type	Description	Default
mol	int	The molecule to retrieve the connection tree for.	required

Returns:

Type	Description
list[int]	A list of molecules in the connection tree.

Source code in src/molclusters/conntable.py

def mols_connected_tree_to(self, mol: int) -> list[int]:
    """Get the molecules in the connection tree of a given molecule.

    Parameters
    ----------
    mol : int
        The molecule to retrieve the connection tree for.

    Returns
    -------
    list[int]
        A list of molecules in the connection tree.
    """
    g = self.__graph_from_mol(mol)
    return list(g.nodes)

__subgraphs

__subgraphs() -> Generator[Graph, None, None]

Generate subgraphs of the connectivity table.

Yields:

Type	Description
Graph	A subgraph of the connectivity table.

Source code in src/molclusters/conntable.py

def __subgraphs(self) -> Generator[nx.Graph, None, None]:
    """Generate subgraphs of the connectivity table.

    Yields
    ------
    nx.Graph
        A subgraph of the connectivity table.
    """
    for c in nx.connected_components(self.conntab):
        yield self.conntab.subgraph(c)

subconntables

subconntables() -> Generator[_SubConnTable, None, None]

Generate sub-connectivity tables from the connectivity table.

Yields:

Type	Description
_SubConnTable	A sub-connectivity table.

Source code in src/molclusters/conntable.py

def subconntables(self) -> Generator["_SubConnTable", None, None]:
    """Generate sub-connectivity tables from the connectivity table.

    Yields
    ------
    _SubConnTable
        A sub-connectivity table.
    """
    # largest first: clusters born in the same frame are numbered by size
    # (largest gets the lowest id), and exact ties in MolClusters go to the first
    for s in sorted(self.__subgraphs(), key=lambda x: len(x), reverse=True):
        if len(s) == 1:
            # TODO: find a way to not change how HB calculations are done to avoid this check
            continue
        yield self._SubConnTable(s, self)

check_resids

check_resids(universe: Universe) -> None

Check that the residues are numbered 1, 2, ..., N in topology order.

MolClusters identifies molecules by resid (in the connection graphs, the clusters and the outputs) and finds their residues back by position, as resid - 1. Any other numbering (an offset, gaps, repeats) would silently put molecules in the wrong clusters and compute their properties on the wrong atoms.

Parameters:

Name	Type	Description	Default
universe	Universe	The Universe to check.	required

Raises:

Type	Description
ValueError	If a residue's resid isn't its position in the topology plus one.

Source code in src/molclusters/conntable.py

def check_resids(universe: mda.Universe) -> None:
    """Check that the residues are numbered 1, 2, ..., N in topology order.

    MolClusters identifies molecules by resid (in the connection graphs, the
    clusters and the outputs) and finds their residues back by position, as
    ``resid - 1``. Any other numbering (an offset, gaps, repeats) would silently
    put molecules in the wrong clusters and compute their properties on the
    wrong atoms.

    Parameters
    ----------
    universe : mda.Universe
        The Universe to check.

    Raises
    ------
    ValueError
        If a residue's resid isn't its position in the topology plus one.
    """
    resids = universe.residues.resids
    wrong = np.flatnonzero(resids != np.arange(1, len(resids) + 1))
    if wrong.size:
        first = int(wrong[0])
        raise ValueError(
            "MolClusters identifies molecules by residue id, and needs the residues "
            f"numbered 1 to {len(resids)} in topology order, but residue {first + 1} "
            f"has resid {resids[first]} ({wrong.size} residue(s) are numbered "
            "otherwise). Renumber the residues in the topology, or from Python with "
            "`universe.residues.resids = numpy.arange(1, len(universe.residues) + 1)`."
        )