subgroup

Subgroup enumeration for crystallographic space groups.

This module provides tools to enumerate translationengleiche (t-) subgroups of a space group, identify their space group type, and expand the asymmetric unit when reducing symmetry.

Primary use case: expanding asymmetric units to contain complete molecules (e.g., urea with Z' = 0.5 -> Z' = 1).

StandardSettingResult dataclass

Result of identifying the standard ITA setting for a set of symops.

Attributes:

Name	Type	Description
sg_number	int	ITA space group number
sg_symbol	str	Short Hermann-Mauguin symbol
choice	str	Space group choice (e.g. "b", "2", "")
basis_transform	ndarray | None	P matrix (3x3), None if identity
origin_shift	ndarray | None	p vector (3,), None if zero shift
target_symops	tuple[int, ...]	Integer codes of the matched standard setting

Source code in chmpy/crystal/subgroup.py

@dataclass(frozen=True)
class StandardSettingResult:
    """Result of identifying the standard ITA setting for a set of symops.

    Attributes:
        sg_number: ITA space group number
        sg_symbol: Short Hermann-Mauguin symbol
        choice: Space group choice (e.g. "b", "2", "")
        basis_transform: P matrix (3x3), None if identity
        origin_shift: p vector (3,), None if zero shift
        target_symops: Integer codes of the matched standard setting
    """

    sg_number: int
    sg_symbol: str
    choice: str
    basis_transform: np.ndarray | None
    origin_shift: np.ndarray | None
    target_symops: tuple[int, ...]

SubgroupEnumerator

Enumerate and verify t-subgroups of a space group.

Uses BFS with closure checking to enumerate translationengleiche subgroups (same lattice, fewer point operations).

Attributes:

Name	Type	Description
sg_table		Precomputed SpaceGroupTable for the parent group
parent_sg		The parent SpaceGroup

Source code in chmpy/crystal/subgroup.py

class SubgroupEnumerator:
    """Enumerate and verify t-subgroups of a space group.

    Uses BFS with closure checking to enumerate translationengleiche
    subgroups (same lattice, fewer point operations).

    Attributes:
        sg_table: Precomputed SpaceGroupTable for the parent group
        parent_sg: The parent SpaceGroup
    """

    def __init__(self, sg_table: SpaceGroupTable, parent_sg: SpaceGroup):
        self.sg_table = sg_table
        self.parent_sg = parent_sg

    @classmethod
    def from_space_group(cls, sg: SpaceGroup) -> SubgroupEnumerator:
        """Create a SubgroupEnumerator from a SpaceGroup.

        Args:
            sg: The space group to enumerate subgroups of

        Returns:
            SubgroupEnumerator ready to enumerate subgroups
        """
        sg_table = SpaceGroupTable.from_space_group(sg)
        return cls(sg_table, sg)

    def enumerate_all(self, max_index: int = 8) -> list[SubgroupResult]:
        """Enumerate all t-subgroups up to a given index.

        Uses BFS: starts from the identity, iteratively adds generators,
        computes closure, and collects unique subgroups.

        Args:
            max_index: Maximum index [G:H] to consider. Subgroups with
                |G|/|H| > max_index are excluded.

        Returns:
            List of SubgroupResult for each unique subgroup found,
            sorted by index (ascending) then by number of symops.
        """
        n_ops = self.sg_table.n_ops
        identity_idx = self.sg_table.identity_idx()
        min_size = max(1, n_ops // max_index)

        # BFS: each state is (closure, last_generator_tried)
        # We only try generators with index > last_gen to avoid duplicates
        # visited tracks all closures we've seen (to avoid re-exploring)
        # found_subgroups collects closures that meet the min_size filter
        visited: set[frozenset[int]] = set()
        found_subgroups: set[frozenset[int]] = set()
        queue: list[tuple[frozenset[int], int]] = [
            (frozenset([identity_idx]), 0)
        ]

        while queue:
            current, last_gen = queue.pop(0)
            closure = compute_closure(current, self.sg_table)

            if closure in visited:
                continue
            visited.add(closure)

            if len(closure) >= min_size:
                found_subgroups.add(closure)

            # Try adding each element as new generator
            for g in range(last_gen, n_ops):
                if g not in closure:
                    queue.append((closure | {g}, g + 1))

        # Build results
        results = []
        full_group = frozenset(range(n_ops))
        for subgroup in found_subgroups:
            if subgroup == full_group:
                continue  # Skip the parent group itself

            index = n_ops // len(subgroup)
            if index > max_index:
                continue

            sg_number, sg_symbol, pg_symbol = _identify_space_group(
                tuple(subgroup), self.sg_table
            )
            results.append(
                SubgroupResult(
                    symop_indices=tuple(sorted(subgroup)),
                    index=index,
                    space_group_number=sg_number,
                    space_group_symbol=sg_symbol,
                    point_group_symbol=pg_symbol,
                    z_prime_factor=float(index),
                )
            )

        results.sort(key=lambda r: (r.index, len(r.symop_indices)))
        return results

    def find_by_index(self, index: int) -> list[SubgroupResult]:
        """Find all subgroups with a specific index.

        Args:
            index: The desired index [G:H]

        Returns:
            List of SubgroupResult with the specified index
        """
        all_subgroups = self.enumerate_all(max_index=index)
        return [sg for sg in all_subgroups if sg.index == index]

    def find_for_target_z_prime(
        self, current_z_prime: float, target: float
    ) -> list[SubgroupResult]:
        """Find subgroups that achieve a target Z'.

        Args:
            current_z_prime: Current Z' value
            target: Desired Z' value

        Returns:
            List of SubgroupResult that would produce the target Z'
        """
        if target <= current_z_prime:
            return []
        ratio = target / current_z_prime
        # ratio must be close to an integer (the index)
        index = round(ratio)
        if abs(ratio - index) > 0.01:
            return []
        return self.find_by_index(index)

    def verify_is_subgroup(self, symop_indices: Sequence[int]) -> bool:
        """Verify that a set of symop indices forms a valid subgroup.

        Checks closure, identity, and inverses using the Cayley table.

        Args:
            symop_indices: Indices to verify

        Returns:
            True if the indices form a valid subgroup
        """
        idx_set = set(symop_indices)
        n_ops = self.sg_table.n_ops

        # Check all indices are valid
        if not all(0 <= i < n_ops for i in idx_set):
            return False

        # Check identity
        if self.sg_table.identity_idx() not in idx_set:
            return False

        # Check closure and inverses
        for a in idx_set:
            # Inverse must be in set
            if int(self.sg_table.inverse_table[a]) not in idx_set:
                return False
            # All products must be in set
            for b in idx_set:
                if int(self.sg_table.mult_table[a, b]) not in idx_set:
                    return False

        # Check Lagrange's theorem: |H| divides |G|
        if n_ops % len(idx_set) != 0:
            return False

        return True

enumerate_all(max_index=8)

Enumerate all t-subgroups up to a given index.

Uses BFS: starts from the identity, iteratively adds generators, computes closure, and collects unique subgroups.

Parameters:

Name	Type	Description	Default
max_index	int	Maximum index [G:H] to consider. Subgroups with |G|/|H| > max_index are excluded.	8

Returns:

Type	Description
list[SubgroupResult]	List of SubgroupResult for each unique subgroup found,
list[SubgroupResult]	sorted by index (ascending) then by number of symops.

Source code in chmpy/crystal/subgroup.py

def enumerate_all(self, max_index: int = 8) -> list[SubgroupResult]:
    """Enumerate all t-subgroups up to a given index.

    Uses BFS: starts from the identity, iteratively adds generators,
    computes closure, and collects unique subgroups.

    Args:
        max_index: Maximum index [G:H] to consider. Subgroups with
            |G|/|H| > max_index are excluded.

    Returns:
        List of SubgroupResult for each unique subgroup found,
        sorted by index (ascending) then by number of symops.
    """
    n_ops = self.sg_table.n_ops
    identity_idx = self.sg_table.identity_idx()
    min_size = max(1, n_ops // max_index)

    # BFS: each state is (closure, last_generator_tried)
    # We only try generators with index > last_gen to avoid duplicates
    # visited tracks all closures we've seen (to avoid re-exploring)
    # found_subgroups collects closures that meet the min_size filter
    visited: set[frozenset[int]] = set()
    found_subgroups: set[frozenset[int]] = set()
    queue: list[tuple[frozenset[int], int]] = [
        (frozenset([identity_idx]), 0)
    ]

    while queue:
        current, last_gen = queue.pop(0)
        closure = compute_closure(current, self.sg_table)

        if closure in visited:
            continue
        visited.add(closure)

        if len(closure) >= min_size:
            found_subgroups.add(closure)

        # Try adding each element as new generator
        for g in range(last_gen, n_ops):
            if g not in closure:
                queue.append((closure | {g}, g + 1))

    # Build results
    results = []
    full_group = frozenset(range(n_ops))
    for subgroup in found_subgroups:
        if subgroup == full_group:
            continue  # Skip the parent group itself

        index = n_ops // len(subgroup)
        if index > max_index:
            continue

        sg_number, sg_symbol, pg_symbol = _identify_space_group(
            tuple(subgroup), self.sg_table
        )
        results.append(
            SubgroupResult(
                symop_indices=tuple(sorted(subgroup)),
                index=index,
                space_group_number=sg_number,
                space_group_symbol=sg_symbol,
                point_group_symbol=pg_symbol,
                z_prime_factor=float(index),
            )
        )

    results.sort(key=lambda r: (r.index, len(r.symop_indices)))
    return results

find_by_index(index)

Find all subgroups with a specific index.

Parameters:

Name	Type	Description	Default
index	int	The desired index [G:H]	required

Returns:

Type	Description
list[SubgroupResult]	List of SubgroupResult with the specified index

Source code in chmpy/crystal/subgroup.py

def find_by_index(self, index: int) -> list[SubgroupResult]:
    """Find all subgroups with a specific index.

    Args:
        index: The desired index [G:H]

    Returns:
        List of SubgroupResult with the specified index
    """
    all_subgroups = self.enumerate_all(max_index=index)
    return [sg for sg in all_subgroups if sg.index == index]

find_for_target_z_prime(current_z_prime, target)

Find subgroups that achieve a target Z'.

Parameters:

Name	Type	Description	Default
current_z_prime	float	Current Z' value	required
target	float	Desired Z' value	required

Returns:

Type	Description
list[SubgroupResult]	List of SubgroupResult that would produce the target Z'

Source code in chmpy/crystal/subgroup.py

def find_for_target_z_prime(
    self, current_z_prime: float, target: float
) -> list[SubgroupResult]:
    """Find subgroups that achieve a target Z'.

    Args:
        current_z_prime: Current Z' value
        target: Desired Z' value

    Returns:
        List of SubgroupResult that would produce the target Z'
    """
    if target <= current_z_prime:
        return []
    ratio = target / current_z_prime
    # ratio must be close to an integer (the index)
    index = round(ratio)
    if abs(ratio - index) > 0.01:
        return []
    return self.find_by_index(index)

from_space_group(sg) classmethod

Create a SubgroupEnumerator from a SpaceGroup.

Parameters:

Name	Type	Description	Default
sg	SpaceGroup	The space group to enumerate subgroups of	required

Returns:

Type	Description
SubgroupEnumerator	SubgroupEnumerator ready to enumerate subgroups

Source code in chmpy/crystal/subgroup.py

@classmethod
def from_space_group(cls, sg: SpaceGroup) -> SubgroupEnumerator:
    """Create a SubgroupEnumerator from a SpaceGroup.

    Args:
        sg: The space group to enumerate subgroups of

    Returns:
        SubgroupEnumerator ready to enumerate subgroups
    """
    sg_table = SpaceGroupTable.from_space_group(sg)
    return cls(sg_table, sg)

verify_is_subgroup(symop_indices)

Verify that a set of symop indices forms a valid subgroup.

Checks closure, identity, and inverses using the Cayley table.

Parameters:

Name	Type	Description	Default
symop_indices	Sequence[int]	Indices to verify	required

Returns:

Type	Description
bool	True if the indices form a valid subgroup

Source code in chmpy/crystal/subgroup.py

def verify_is_subgroup(self, symop_indices: Sequence[int]) -> bool:
    """Verify that a set of symop indices forms a valid subgroup.

    Checks closure, identity, and inverses using the Cayley table.

    Args:
        symop_indices: Indices to verify

    Returns:
        True if the indices form a valid subgroup
    """
    idx_set = set(symop_indices)
    n_ops = self.sg_table.n_ops

    # Check all indices are valid
    if not all(0 <= i < n_ops for i in idx_set):
        return False

    # Check identity
    if self.sg_table.identity_idx() not in idx_set:
        return False

    # Check closure and inverses
    for a in idx_set:
        # Inverse must be in set
        if int(self.sg_table.inverse_table[a]) not in idx_set:
            return False
        # All products must be in set
        for b in idx_set:
            if int(self.sg_table.mult_table[a, b]) not in idx_set:
                return False

    # Check Lagrange's theorem: |H| divides |G|
    if n_ops % len(idx_set) != 0:
        return False

    return True

SubgroupResult dataclass

Result of computing a t-subgroup.

Attributes:

Name	Type	Description
symop_indices	tuple[int, ...]	Indices into parent group's symop list
index	int	[G:H] = |G|/|H|, the index of the subgroup
space_group_number	int | None	ITA number if matches known group, else None
space_group_symbol	str | None	Symbol if identified, else None
point_group_symbol	str	Point group symbol (always identified)
z_prime_factor	float	Factor by which Z' increases (equals the index)

Source code in chmpy/crystal/subgroup.py

@dataclass(frozen=True)
class SubgroupResult:
    """Result of computing a t-subgroup.

    Attributes:
        symop_indices: Indices into parent group's symop list
        index: [G:H] = |G|/|H|, the index of the subgroup
        space_group_number: ITA number if matches known group, else None
        space_group_symbol: Symbol if identified, else None
        point_group_symbol: Point group symbol (always identified)
        z_prime_factor: Factor by which Z' increases (equals the index)
    """

    symop_indices: tuple[int, ...]
    index: int
    space_group_number: int | None
    space_group_symbol: str | None
    point_group_symbol: str
    z_prime_factor: float

compute_closure(generators, sg_table)

Compute group closure of generator indices using the Cayley table.

Uses the multiplication and inverse tables for O(1) composition, iteratively expanding the set until no new elements are produced.

Parameters:

Name	Type	Description	Default
generators	frozenset[int]	Frozenset of symop indices (into sg_table)	required
sg_table	SpaceGroupTable	Precomputed space group table with Cayley table	required

Returns:

Type	Description
frozenset[int]	Frozenset of all symop indices in the closure

Source code in chmpy/crystal/subgroup.py

def compute_closure(generators: frozenset[int], sg_table: SpaceGroupTable) -> frozenset[int]:
    """Compute group closure of generator indices using the Cayley table.

    Uses the multiplication and inverse tables for O(1) composition,
    iteratively expanding the set until no new elements are produced.

    Args:
        generators: Frozenset of symop indices (into sg_table)
        sg_table: Precomputed space group table with Cayley table

    Returns:
        Frozenset of all symop indices in the closure
    """
    closure = set(generators)
    changed = True
    while changed:
        changed = False
        new_elements = set()
        for a in closure:
            # Add inverse
            inv_a = int(sg_table.inverse_table[a])
            if inv_a not in closure:
                new_elements.add(inv_a)
            # Add all products with existing elements
            for b in closure:
                product = int(sg_table.mult_table[a, b])
                if product not in closure and product not in new_elements:
                    new_elements.add(product)
        if new_elements:
            closure.update(new_elements)
            changed = True
    return frozenset(closure)

expand_asymmetric_unit(asymmetric_unit, symops, subgroup_indices, sg_table, tolerance=0.0001)

Expand asymmetric unit when reducing to a subgroup.

First deduplicates the input asymmetric unit (removing atoms that are already related by parent group symmetry), then for each unique atom computes its stabilizer, finds coset representatives (symops in G but not accounted for by the subgroup), and generates new independent atoms.

Parameters:

Name	Type	Description	Default
asymmetric_unit	AsymmetricUnit	The original asymmetric unit (may contain redundancy)	required
symops	list[SymmetryOperation]	Parent group's symmetry operations	required
subgroup_indices	Sequence[int]	Indices of symops forming the subgroup	required
sg_table	SpaceGroupTable	Parent group's SpaceGroupTable	required
tolerance	float	Tolerance for position comparison	0.0001

Returns:

Type	Description
AsymmetricUnit	New AsymmetricUnit with expanded atom list

Source code in chmpy/crystal/subgroup.py

def expand_asymmetric_unit(
    asymmetric_unit: AsymmetricUnit,
    symops: list[SymmetryOperation],
    subgroup_indices: Sequence[int],
    sg_table: SpaceGroupTable,
    tolerance: float = 1e-4,
) -> AsymmetricUnit:
    """Expand asymmetric unit when reducing to a subgroup.

    First deduplicates the input asymmetric unit (removing atoms that are
    already related by parent group symmetry), then for each unique atom
    computes its stabilizer, finds coset representatives (symops in G but
    not accounted for by the subgroup), and generates new independent atoms.

    Args:
        asymmetric_unit: The original asymmetric unit (may contain redundancy)
        symops: Parent group's symmetry operations
        subgroup_indices: Indices of symops forming the subgroup
        sg_table: Parent group's SpaceGroupTable
        tolerance: Tolerance for position comparison

    Returns:
        New AsymmetricUnit with expanded atom list
    """
    from .asymmetric_unit import AsymmetricUnit
    from .site_symmetry import SiteSymmetry

    # Step 1: Deduplicate the input asymmetric unit
    unique_indices, _ = _deduplicate_asymmetric_unit(
        asymmetric_unit, symops, tolerance
    )
    if len(unique_indices) < len(asymmetric_unit.elements):
        LOG.info(
            "Deduplicated asymmetric unit: %d -> %d atoms",
            len(asymmetric_unit.elements),
            len(unique_indices),
        )

    subgroup_set = set(subgroup_indices)
    n_parent = len(symops)
    n_subgroup = len(subgroup_indices)

    new_elements = []
    new_positions = []
    new_labels = []

    for i in unique_indices:
        elem = asymmetric_unit.elements[i]
        pos = asymmetric_unit.positions[i]

        # Compute stabilizer in parent group
        site_sym = SiteSymmetry.from_position(pos, symops, tolerance)
        stabilizer_codes = set(site_sym.stabilizer_symop_codes)

        # Map stabilizer codes to indices
        stabilizer_indices = set()
        for code in stabilizer_codes:
            if code in sg_table.symop_to_idx:
                stabilizer_indices.add(sg_table.symop_to_idx[code])

        # Compute stabilizer in subgroup = Stab_G(x) ∩ H
        stab_in_subgroup = stabilizer_indices & subgroup_set

        # Multiplicity in parent = |G| / |Stab_G|
        mult_parent = n_parent // len(stabilizer_indices)
        # Multiplicity in subgroup = |H| / |Stab_H|
        mult_subgroup = n_subgroup // len(stab_in_subgroup) if stab_in_subgroup else n_subgroup

        # Number of new independent copies = mult_parent / mult_subgroup
        # This is how many orbits the single parent orbit splits into
        # under the subgroup. Equivalently, [Stab_G(x) : Stab_G(x) ∩ H].
        n_copies = mult_parent // mult_subgroup if mult_parent > mult_subgroup else 1

        label_base = (
            asymmetric_unit.labels[i]
            if i < len(asymmetric_unit.labels)
            else f"{elem}{i+1}"
        )

        if n_copies <= 1:
            # No expansion needed for this atom
            new_elements.append(elem)
            new_positions.append(pos.copy())
            new_labels.append(str(label_base))
            continue

        # Find coset representatives: we need n_copies distinct images.
        # Always start with the original position, then find additional
        # copies nearby (closest periodic image to original).
        seen_positions = [pos.copy()]
        new_elements.append(elem)
        new_positions.append(pos.copy())
        new_labels.append(str(label_base))

        for g_idx in range(n_parent):
            raw_pos = symops[g_idx].apply(pos.reshape(1, 3)).flatten()

            # Find nearest periodic image to original position
            new_pos = _nearest_image(raw_pos, pos)

            # Check if this position is equivalent to any seen position
            # under the subgroup operations
            is_new = True
            for seen_pos in seen_positions:
                for h_idx in subgroup_indices:
                    h_pos = symops[h_idx].apply(seen_pos.reshape(1, 3)).flatten()
                    diff = new_pos - h_pos
                    diff = diff - np.round(diff)
                    if np.allclose(diff, 0.0, atol=tolerance):
                        is_new = False
                        break
                if not is_new:
                    break

            if is_new:
                seen_positions.append(new_pos)
                new_elements.append(elem)
                new_positions.append(new_pos)
                suffix = chr(ord('a') + len(seen_positions) - 1)
                new_labels.append(f"{label_base}{suffix}")

                if len(seen_positions) >= n_copies:
                    break

    positions_array = np.array(new_positions)
    return AsymmetricUnit(
        new_elements, positions_array, labels=new_labels
    )

identify_standard_setting(symops)

Identify which standard ITA setting a set of symops corresponds to.

Uses a multi-tier strategy via _identify_match(): 1. Exact code match — already standard 2. Origin shift only (P=I) 3. Basis transform P + exact match 4. Basis transform P + origin shift p 5. Centering reduction + eigenvector transforms

Parameters:

Name	Type	Description	Default
symops	list[SymmetryOperation]	List of SymmetryOperation objects	required

Returns:

Type	Description
StandardSettingResult | None	StandardSettingResult if a match is found, None otherwise

Source code in chmpy/crystal/subgroup.py

def identify_standard_setting(
    symops: list[SymmetryOperation],
) -> StandardSettingResult | None:
    """Identify which standard ITA setting a set of symops corresponds to.

    Uses a multi-tier strategy via _identify_match():
    1. Exact code match — already standard
    2. Origin shift only (P=I)
    3. Basis transform P + exact match
    4. Basis transform P + origin shift p
    5. Centering reduction + eigenvector transforms

    Args:
        symops: List of SymmetryOperation objects

    Returns:
        StandardSettingResult if a match is found, None otherwise
    """
    sub_ops = [(s.rotation, s.translation) for s in symops]

    result = _identify_match(sub_ops)
    if result is None:
        return None

    sgdata, origin_shift, basis_transform = result
    has_shift = origin_shift is not None and not np.allclose(origin_shift, 0, atol=1e-6)
    return StandardSettingResult(
        sg_number=sgdata.number,
        sg_symbol=sgdata.short,
        choice=sgdata.choice,
        basis_transform=basis_transform,
        origin_shift=origin_shift if has_shift else None,
        target_symops=tuple(sgdata.symops),
    )