vibe-sketcherDocumentation
Source code for vibesketcher.render.symmetry_elements
"""Symmetry elements of a space group in a projection along a cell axis (#233).
Each operation (R, t) of the space group, combined with lattice translations,
is analysed geometrically: its type (rotation, screw, inversion, rotoinversion,
mirror, glide), the location of its fixed point set and its intrinsic
(screw or glide) translation, following International Tables for
Crystallography Vol. A, Section 1.2 (geometric interpretation of
symmetry operations). Elements parallel or perpendicular to the projection
axis are returned with their position in the projected cell, as in the
space-group diagrams of International Tables; oblique ones are counted.
Coordinates are fractional in the structure's cell; the projection axis is
one of its basis vectors ``k``, and the in-plane axes are the other two.
"""
from __future__ import annotations
import itertools
from dataclasses import dataclass
from fractions import Fraction
from typing import Any
import numpy as np
from vibesketcher.document import crystal
from vibesketcher.document.model import PeriodicStructure
_EPS = 1e-6
_ITA = "International Tables for Crystallography Vol. A, Sections 1.2 and 1.4"
[docs]
@dataclass(frozen=True)
class Element:
"""One symmetry element as it appears in the projection."""
kind: str # "axis", "inplane_axis", "plane", "horizontal_plane", "inversion", "rotoinversion"
symbol: str # 2, 2_1, 4_2, -4, m, a, b, c, n, d, e, -1
#: In-plane fractional coordinates: a point (axis, inversion), or two end points (lines).
points: tuple[tuple[float, float], ...]
#: Height along the projection axis (fraction of its period), where it is defined.
height: float | None = None
#: For in-plane axes and horizontal planes: the in-plane direction of a screw or glide.
detail: str = ""
[docs]
def operations_of(
structure: PeriodicStructure, symprec: float = 1e-3
) -> tuple[list[tuple[np.ndarray, np.ndarray]], str]:
"""The operations in the structure's basis and where they came from."""
record = crystal._builder(structure)
group = record.get("space_group")
setting = group.get("setting") if isinstance(group, dict) and "setting" in group else group
hall = setting.get("hall") if isinstance(setting, dict) else None
if hall and structure.cell_setting == "conventional" and record.get("kind") in (
"crystal", "cif_import", None): # fmt: skip
from vibesketcher.document.spacegroups import find
found = find(str(hall))
ops = [
(np.asarray(op.rotation, dtype=float), np.asarray([float(x) for x in op.translation]))
for op in found.operations()
]
return ops, f"{found.hm} (No. {found.number}), recorded"
from vibesketcher.document import symmetry
if not symmetry.available():
raise ValueError("this structure records no space group; install the symmetry extra "
"(spglib) to detect one") # fmt: skip
report = symmetry.detect(structure, symprec=symprec)
return (symmetry.operations(structure, symprec=symprec),
f"{report['international']} (No. {report['number']}), detected by spglib") # fmt: skip
def _order(r: np.ndarray) -> int:
m = np.eye(3)
for n in range(1, 7):
m = m @ r
if np.allclose(m, np.eye(3)):
return n
raise ValueError("not a crystallographic operation")
def _axis(r: np.ndarray) -> np.ndarray:
"""The direction kept by a proper rotation (or by -R for a rotoinversion), as integers."""
values, vectors = np.linalg.eig(r)
k = int(np.argmin(np.abs(values - 1)))
u = np.real(vectors[:, k])
u /= u[np.argmax(np.abs(u))]
return np.round(u * 6) / 6
def _solve(a: np.ndarray, b: np.ndarray) -> np.ndarray | None:
x, *_ = np.linalg.lstsq(a, b, rcond=None)
return x if np.allclose(a @ x, b, atol=1e-6) else None
def _frac(x: float) -> float:
return float(Fraction(x).limit_denominator(24))
[docs]
def analyse(r: np.ndarray, t: np.ndarray) -> dict[str, Any] | None:
"""The geometric element of x -> R x + t (fractional), or None for a pure translation."""
det = round(float(np.linalg.det(r)))
n = _order(r)
if n == 1:
return None
if det == 1:
w = sum(np.linalg.matrix_power(r, k) @ t for k in range(n)) / n
point = _solve(r - np.eye(3), -(t - w))
u = _axis(r)
# Rotation sense about u: R and R^-1 describe one element, so a screw is
# read for the positive sense (a -120 degree 3_1 is still 3_1).
v = np.eye(3)[int(np.argmin(np.abs(u)))]
sense = 1.0 if np.linalg.det(np.column_stack([u, v, r @ v])) > 0 else -1.0
return {"type": "rotation", "order": n, "axis": u, "screw": w * sense, "point": point}
if np.allclose(r, -np.eye(3)):
return {"type": "inversion", "point": t / 2}
if n == 2: # a mirror or glide plane: R^2 = 1, det -1
w = (t + r @ t) / 2
point = _solve(r - np.eye(3), -(t - w))
normal = (r - np.eye(3))[np.argmax(np.abs(r - np.eye(3)).sum(axis=1))]
return {"type": "plane", "normal": normal, "glide": w, "point": point}
# Rotoinversion -n = inversion after an n-fold rotation, so -R is that rotation;
# its fixed point is unique (ITA Vol. A, Section 1.2).
point = _solve(r - np.eye(3), -t)
return {"type": "rotoinversion", "order": _order(-r), "axis": _axis(-r), "point": point}
def _hexagonal(structure: PeriodicStructure) -> bool:
rows = crystal.lattice_rows(structure)
lengths = np.linalg.norm(rows, axis=1)
return any(
abs(float(rows[a] @ rows[b]) / (lengths[a] * lengths[b]) + 0.5) < 1e-4
for a, b in ((0, 1), (1, 2), (0, 2))
)
[docs]
def elements(
structure: PeriodicStructure, axis: int, *, symprec: float = 1e-3
) -> tuple[list[Element], dict[str, Any]]:
"""Symmetry elements in the projection along cell vector ``axis`` (0, 1 or 2)."""
if structure.pbc != (True, True, True):
raise ValueError("symmetry-element diagrams are drawn for three-dimensional crystals")
ops, source = operations_of(structure, symprec)
i, j = [a for a in range(3) if a != axis]
found: dict[tuple[Any, ...], Element] = {}
oblique: set[int] = set()
shifts = list(itertools.product((-1, 0, 1), repeat=3))
for index, ((r, t0), shift) in enumerate(itertools.product(ops, shifts)):
info = analyse(r, t0 + np.asarray(shift, dtype=float))
if info is None:
continue
made = _project(info, axis, i, j, _hexagonal(structure))
if made is None:
oblique.add(index // len(shifts))
continue
for element in made:
for di, dj in itertools.product((-1, 0, 1), repeat=2):
moved = _shifted(element, di, dj)
if moved is not None:
found.setdefault(_key(moved), moved)
summary = {"source": source, "projection_axis": "abc"[axis], "operations": len(ops),
"oblique_operations_not_drawn": len(oblique),
"reference": _ITA} # fmt: skip
return _highest(sorted(found.values(), key=_key)), summary
def _place(e: Element) -> tuple[Any, ...]:
points = tuple(sorted(e.points)) if e.kind == "plane" else e.points
return (e.kind, points, None if e.height is None else round(e.height % 1.0, 4))
def _rank(symbol: str) -> int:
return int(symbol.lstrip("-").split("_")[0])
def _highest(found: list[Element]) -> list[Element]:
"""At one point keep the highest axis (a 4 contains a 2), as International Tables draws.
A -4 or -6 on a plain 4 or 6 axis (4/m, 6/m) is dropped, since the axis, the
horizontal mirror and the inversion centres show it; a -n on a screw axis is
kept. A -3 replaces the 3 it contains. Inversion centres keep their own marks.
"""
axes: dict[tuple[float, float], Element] = {}
for e in found:
if e.kind == "axis":
point = e.points[0]
best = axes.get(point)
if best is None or _rank(e.symbol) > _rank(best.symbol):
axes[point] = e
# Coincident planes: a mirror stays m; two different glides make a double glide e.
planes: dict[tuple[Any, ...], set[str]] = {}
for e in found:
if e.kind in ("plane", "horizontal_plane"):
planes.setdefault(_place(e), set()).add(e.symbol)
shown: set[tuple[Any, ...]] = set()
out = []
for e in found:
if e.kind in ("plane", "horizontal_plane"):
place = _place(e)
if place in shown:
continue
shown.add(place)
letters = planes[place]
symbol = "m" if "m" in letters else "e" if len(letters) > 1 else e.symbol
out.append(Element(e.kind, symbol, e.points, e.height, e.detail))
continue
if e.kind == "axis" and axes[e.points[0]] is not e:
continue
if (
e.kind == "axis"
and e.symbol == "3"
and any(
r.kind == "rotoinversion" and r.symbol == "-3" and r.points == e.points
for r in found
)
):
continue # -3 contains the 3 and is what the tables draw
if e.kind == "rotoinversion" and e.symbol != "-3":
proper = axes.get(e.points[0])
if proper is not None and proper.symbol == e.symbol.lstrip("-"):
continue # 4 and -4 together are 4/m: the axis, the mirror and the centres show it
out.append(e)
return out
def _inside(p: tuple[float, float]) -> bool:
return all(-_EPS <= c <= 1 + _EPS for c in p)
def _shifted(e: Element, di: int, dj: int) -> Element | None:
if e.kind == "horizontal_plane":
return e if di == dj == 0 else None # one plane, drawn once at the cell corner
points = tuple((round(x + di, 6), round(y + dj, 6)) for x, y in e.points)
if e.kind in ("plane", "inplane_axis"):
clipped = _clip(points[0], points[1])
if clipped is None:
return None
points = clipped
elif not all(_inside(p) for p in points):
return None
return Element(e.kind, e.symbol, points, e.height, e.detail)
def _clip(a: tuple[float, float], b: tuple[float, float]) -> tuple[tuple[float, float], ...] | None:
"""The part of the infinite line through a, b inside the unit square."""
pa, d = np.asarray(a), np.asarray(b) - np.asarray(a)
lo, hi = -np.inf, np.inf
for k in range(2):
if abs(d[k]) < _EPS:
if pa[k] < -_EPS or pa[k] > 1 + _EPS:
return None
continue
s1, s2 = (0 - pa[k]) / d[k], (1 - pa[k]) / d[k]
lo, hi = max(lo, min(s1, s2)), min(hi, max(s1, s2))
if hi - lo < 1e-6:
return None
p, q = pa + lo * d, pa + hi * d
return (
(round(float(p[0]), 6), round(float(p[1]), 6)),
(round(float(q[0]), 6), round(float(q[1]), 6)),
)
def _key(e: Element) -> tuple[Any, ...]:
points = tuple(sorted(e.points)) if e.kind in ("plane", "inplane_axis") else e.points
height = None if e.height is None else round(e.height % 1.0, 4)
return (e.kind, e.symbol, points, height, e.detail)
def _screw_symbol(order: int, along: float) -> str:
m = round(order * (along % 1.0)) % order
return f"{order}_{m}" if m else str(order)
def _project(
info: dict[str, Any], axis: int, i: int, j: int, hexagonal: bool = False
) -> list[Element] | None:
kind = info["type"]
point = info.get("point")
if kind == "inversion":
p = info["point"]
return [Element("inversion", "-1", ((round(p[i], 6), round(p[j], 6)),), _frac(p[axis] % 1))]
if point is None:
return None
if kind in ("rotation", "rotoinversion"):
u = info["axis"]
if abs(u[i]) < _EPS and abs(u[j]) < _EPS: # parallel to the projection axis
# Move the point along the axis to height 0 (the axis is a line).
p = point - u * (point[axis] / u[axis])
if kind == "rotation":
symbol = _screw_symbol(info["order"], float(info["screw"][axis] / u[axis]))
return [Element("axis", symbol, ((round(p[i], 6), round(p[j], 6)),))]
where = ((round(point[i], 6), round(point[j], 6)),)
return [Element("rotoinversion", f"-{info['order']}", where, _frac(point[axis] % 1))]
if abs(u[axis]) < _EPS and kind == "rotation" and info["order"] == 2:
a = (point[i], point[j])
b = (point[i] + u[i], point[j] + u[j])
screw = info["screw"]
symbol = "2_1" if np.linalg.norm(screw) > _EPS else "2"
return [Element("inplane_axis", symbol, (a, b), _frac(point[axis] % 1))]
return None
# Planes: normal (a covector g with g . x = g . point).
g = info["normal"]
w = info["glide"]
letter = _glide_letter(w, hexagonal)
if abs(g[axis]) < _EPS: # the plane contains the projection axis: a line
if abs(g[i]) < _EPS and abs(g[j]) < _EPS:
return None
c = float(g @ point)
if abs(g[j]) > _EPS:
a = (0.0, (c - g[i] * 0.0) / g[j])
b = (1.0, (c - g[i] * 1.0) / g[j])
else:
a = (c / g[i], 0.0)
b = (c / g[i], 1.0)
return [Element("plane", letter, (a, b))]
if abs(g[i]) < _EPS and abs(g[j]) < _EPS: # perpendicular to the projection axis
height = float(g @ point / g[axis])
return [Element("horizontal_plane", letter, ((0.0, 0.0),), _frac(height % 1),
detail=_direction_text(w, i, j))] # fmt: skip
return None
def _glide_letter(w: np.ndarray, hexagonal: bool = False) -> str:
"""m, a, b, c, n, d or g from the intrinsic glide vector (ITA Vol. A, Table 1.2.1.1)."""
parts = [abs(x % 1.0) > _EPS and abs(x % 1.0 - 1) > _EPS for x in w]
if not any(parts):
return "m"
if sum(parts) == 1:
return "abc"[parts.index(True)]
if hexagonal:
return "g"
nonzero = [x % 1.0 for x, p in zip(w, parts, strict=True) if p]
if all(abs(x - 0.5) < 1e-6 for x in nonzero):
return "n"
if all(min(abs(x - 0.25), abs(x - 0.75)) < 1e-6 for x in nonzero):
return "d"
return "g"
def _direction_text(w: np.ndarray, i: int, j: int) -> str:
return ",".join(f"{_frac(w[k] % 1):g}" for k in (i, j))