diff --git a/docs/api/api_symmetry.md b/docs/api/api_symmetry.md index 5592ce8..73ba880 100644 --- a/docs/api/api_symmetry.md +++ b/docs/api/api_symmetry.md @@ -46,3 +46,11 @@ :members: :undoc-members: ``` + +## Subgroup enumeration + +```{eval-rst} + .. automodule:: spgrep.symmetry.subgroup + :members: + :undoc-members: +``` diff --git a/pyproject.toml b/pyproject.toml index add0864..0597906 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -30,6 +30,7 @@ dynamic = [ "version", ] dependencies = [ + "hsnf>=0.3.16", "numpy>=1.20.1", "spglib>=2.7", "typing-extensions", diff --git a/src/spgrep/symmetry/subgroup.py b/src/spgrep/symmetry/subgroup.py new file mode 100644 index 0000000..d25ea1c --- /dev/null +++ b/src/spgrep/symmetry/subgroup.py @@ -0,0 +1,125 @@ +"""Isotropy subgroup of space group.""" + +from __future__ import annotations + +from queue import Queue + +from spgrep.rep.group import get_identity_index, get_inverse_index +from spgrep.utils import ( + NDArrayInt, +) + + +def enumerate_point_subgroup( + table: NDArrayInt, preserve_sublattice: list[bool], return_conjugacy_class: bool = True +) -> list[list[int]]: + """Enumerate conjugacy subgroups of point group. + + Parameters + ---------- + table: array[int], (order, order) + Multiplication table of group + preserve_sublattice: list[bool] + Specify ``preserve_sublattice[i] = True`` if the ``i``-th operation preserves translational subgroup of isotropy subgroup + return_conjugacy_class: bool, default=True + If true, return representatives of conjugacy classes. + + Returns + ------- + subgroups: list[list[int]] + """ + order = len(table) + identity = get_identity_index(table) + # Represent choice of elements by bit array + st = {1 << identity} + for i in range(order): + if (not preserve_sublattice[i]) or (i == identity): + continue + if (1 << i) in st: + # Already visited + continue + + next_st = set() + for bits in st: + elements = _decode_bits(bits, order) + assert _is_subgroup(elements, table) + generated = _traverse(elements + [i], identity, table) + next_st.add(sum(1 << idx for idx in set(generated))) + + st = st.union(next_st) + + if not return_conjugacy_class: + subgroups = [] + for bits in sorted(st): + subgroups.append(_decode_bits(bits, order)) + return subgroups + + # Group by conjugacy classes + found = set() + ret = [] + for bits in sorted(st): + if bits in found: + continue + elements = _decode_bits(bits, order) + ret.append(elements) + for i in range(order): + if not preserve_sublattice[i]: + continue + inv = get_inverse_index(table, i) + conj = [int(table[inv, table[idx, i]]) for idx in elements] + found.add(sum(1 << idx for idx in set(conj))) + + assert found == st + return ret + + +def enumerate_point_subgroup_naive(table, preserve_sublattice: list[bool]): + """Enumerate conjugacy subgroups of point group in brute force.""" + order = len(table) + ret = [] + for bits in range(1, 1 << order): + elements = _decode_bits(bits, order) + if not all([preserve_sublattice[idx] for idx in elements]): + continue + + if _is_subgroup(elements, table): + ret.append(bits) + + return ret + + +def _decode_bits(bits: int, order: int) -> list[int]: + elements = [idx for idx in range(order) if (bits >> idx) & 1 == 1] + return elements + + +def _is_subgroup(elements: list[int], table: NDArrayInt) -> bool: + subtable = table[elements][:, elements] + for i in range(len(subtable)): + if (set(subtable[i]) != set(elements)) or (set(subtable[:, i]) != set(elements)): + return False + return True + + +def _traverse( + generators: list[int], + identity: int, + table: NDArrayInt, +) -> list[int]: + """Traverse group elements from generators.""" + visited = [False for _ in range(len(table))] + que = Queue() # type: ignore + que.put(identity) + + while not que.empty(): + g = que.get() + if visited[g]: + continue + visited[g] = True + + for h in generators: + gh = int(table[g, h]) # cast np.int64 to int + if not visited[gh]: + que.put(gh) + + return sorted([i for i, v in enumerate(visited) if v]) diff --git a/tests/symmetry/test_symmetry_isotropy.py b/tests/symmetry/test_symmetry_isotropy.py new file mode 100644 index 0000000..6671a34 --- /dev/null +++ b/tests/symmetry/test_symmetry_isotropy.py @@ -0,0 +1,16 @@ +from __future__ import annotations + +import numpy as np + +from spgrep.symmetry.group import get_cayley_table +from spgrep.symmetry.pointgroup import pg_dataset +from spgrep.symmetry.subgroup import enumerate_point_subgroup, enumerate_point_subgroup_naive + + +def test_enumerate_point_subgroup(): + pointgroup = pg_dataset["4/mmm"][0] + table = get_cayley_table(np.array(pointgroup)) + flags = [True for _ in range(len(pointgroup))] + subgroups_actual = enumerate_point_subgroup(table, flags, return_conjugacy_class=False) + subgroups_expect = enumerate_point_subgroup_naive(table, flags) + assert len(subgroups_actual) == len(subgroups_expect) diff --git a/uv.lock b/uv.lock index 21c6cd4..2c7ceb5 100644 --- a/uv.lock +++ b/uv.lock @@ -887,6 +887,25 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/d3/b7/4a806f85d62c20157e62e58e03b27513dc9c55499768530acc4f4c5ce4be/h5py-3.15.1-cp314-cp314-win_arm64.whl", hash = "sha256:a6d8c5a05a76aca9a494b4c53ce8a9c29023b7f64f625c6ce1841e92a362ccdf", size = 2465544, upload-time = "2025-10-16T10:35:25.695Z" }, ] +[[package]] +name = "hsnf" +version = "0.3.16" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, + { name = "numpy", version = "2.4.1", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, + { name = "scipy", version = "1.15.3", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, + { name = "scipy", version = "1.17.0", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, + { name = "setuptools" }, + { name = "setuptools-scm" }, + { name = "typing-extensions" }, + { name = "wheel" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/54/88/7461240d61cfaa7301163e6fc6b78e8ac5ced3dbda2b708476cec2b6978b/hsnf-0.3.16.tar.gz", hash = "sha256:ae4f99c89076f734969d6aedaa0d3390f667929a2aa4862c7e03c067cd063207", size = 21153, upload-time = "2023-02-24T14:14:49.914Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/b5/90/f6f56e32a71a091ff0afa377f92a97a2a3804dc8ab8210134a94144e9bf8/hsnf-0.3.16-py3-none-any.whl", hash = "sha256:906d0ad4af691e1c45e51fd6c99d8b3ed0be1179cbcc21cbd4dda11cc97dc1e8", size = 11053, upload-time = "2023-02-24T14:14:47.633Z" }, +] + [[package]] name = "httpcore" version = "1.0.9" @@ -2960,6 +2979,20 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/a3/dc/17031897dae0efacfea57dfd3a82fdd2a2aeb58e0ff71b77b87e44edc772/setuptools-80.9.0-py3-none-any.whl", hash = "sha256:062d34222ad13e0cc312a4c02d73f059e86a4acbfbdea8f8f76b28c99f306922", size = 1201486, upload-time = "2025-05-27T00:56:49.664Z" }, ] +[[package]] +name = "setuptools-scm" +version = "9.2.2" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "packaging" }, + { name = "setuptools" }, + { name = "tomli", marker = "python_full_version < '3.11'" }, +] +sdist = { url = "https://files.pythonhosted.org/packages/7b/b1/19587742aad604f1988a8a362e660e8c3ac03adccdb71c96d86526e5eb62/setuptools_scm-9.2.2.tar.gz", hash = "sha256:1c674ab4665686a0887d7e24c03ab25f24201c213e82ea689d2f3e169ef7ef57", size = 203385, upload-time = "2025-10-19T22:08:05.608Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/3d/ea/ac2bf868899d0d2e82ef72d350d97a846110c709bacf2d968431576ca915/setuptools_scm-9.2.2-py3-none-any.whl", hash = "sha256:30e8f84d2ab1ba7cb0e653429b179395d0c33775d54807fc5f1dd6671801aef7", size = 62975, upload-time = "2025-10-19T22:08:04.007Z" }, +] + [[package]] name = "six" version = "1.17.0" @@ -3044,6 +3077,7 @@ wheels = [ name = "spgrep" source = { editable = "." } dependencies = [ + { name = "hsnf" }, { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, { name = "numpy", version = "2.4.1", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, { name = "spglib" }, @@ -3093,6 +3127,7 @@ testing = [ [package.metadata] requires-dist = [ { name = "docutils", marker = "extra == 'docs'" }, + { name = "hsnf", specifier = ">=0.3.16" }, { name = "ipykernel", marker = "extra == 'dev'" }, { name = "ipython", marker = "extra == 'dev'" }, { name = "linkify-it-py", marker = "extra == 'docs'" }, @@ -3733,3 +3768,12 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/9a/3f/f70e03f40ffc9a30d817eef7da1be72ee4956ba8d7255c399a01b135902a/websockets-16.0-pp311-pypy311_pp73-win_amd64.whl", hash = "sha256:a653aea902e0324b52f1613332ddf50b00c06fdaf7e92624fbf8c77c78fa5767", size = 178735, upload-time = "2026-01-10T09:23:42.259Z" }, { url = "https://files.pythonhosted.org/packages/6f/28/258ebab549c2bf3e64d2b0217b973467394a9cea8c42f70418ca2c5d0d2e/websockets-16.0-py3-none-any.whl", hash = "sha256:1637db62fad1dc833276dded54215f2c7fa46912301a24bd94d45d46a011ceec", size = 171598, upload-time = "2026-01-10T09:23:45.395Z" }, ] + +[[package]] +name = "wheel" +version = "0.45.1" +source = { registry = "https://pypi.org/simple" } +sdist = { url = "https://files.pythonhosted.org/packages/8a/98/2d9906746cdc6a6ef809ae6338005b3f21bb568bea3165cfc6a243fdc25c/wheel-0.45.1.tar.gz", hash = "sha256:661e1abd9198507b1409a20c02106d9670b2576e916d58f520316666abca6729", size = 107545, upload-time = "2024-11-23T00:18:23.513Z" } +wheels = [ + { url = "https://files.pythonhosted.org/packages/0b/2c/87f3254fd8ffd29e4c02732eee68a83a1d3c346ae39bc6822dcbcb697f2b/wheel-0.45.1-py3-none-any.whl", hash = "sha256:708e7481cc80179af0e556bbf0cc00b8444c7321e2700b8d8580231d13017248", size = 72494, upload-time = "2024-11-23T00:18:21.207Z" }, +]