vibe-sketcherDocumentation
Source code for vibesketcher.document.symmetry
"""Space-group detection, standardisation and cell reduction with spglib (#199).
Optional: ``pip install 'vibesketcher[symmetry]'`` installs spglib 2.8.0
(BSD-3-Clause; A. Togo, K. Shinohara and I. Tanaka, Sci. Technol. Adv. Mater.
Meth. 4, 2384822 (2024)). Building crystals *from* a space group does not need
it (:mod:`vibesketcher.document.spacegroups`).
Detection reports at a stated length tolerance and never changes the
structure. Standardisation returns a new structure in the standard
conventional or primitive setting; by default positions are *not* idealised,
and the displacement idealisation would apply is reported so it can be
reviewed first. Only three-dimensional, fully ordered cells are analysed: a
slab's layer group is not a space group of its vacuum-padded cell.
"""
from __future__ import annotations
import importlib
from typing import Any, Literal
import numpy as np
from vibesketcher.document import crystal
from vibesketcher.document.model import PeriodicStructure
from vibesketcher.elements import ATOMIC_NUMBER
HINT = "symmetry analysis needs the symmetry extra: pip install 'vibesketcher[symmetry]'"
REFERENCE = "A. Togo, K. Shinohara and I. Tanaka, Sci. Technol. Adv. Mater. Meth. 4, 2384822 (2024)"
def _spglib() -> Any:
try:
return importlib.import_module("spglib")
except ImportError:
raise SymmetryUnavailable(HINT) from None
def _cell(structure: PeriodicStructure) -> tuple[np.ndarray, np.ndarray, list[int]]:
if structure.pbc != (True, True, True):
raise ValueError(
"space groups describe three-dimensional crystals; a slab or wire has a layer or "
"rod group, which is not analysed here"
)
if any(site.occupancy != 1 for site in structure.sites):
raise ValueError("partially occupied sites have no single space group; resolve disorder")
lattice = crystal.lattice_rows(structure)
positions = np.asarray([s.fractional for s in structure.sites], dtype=float)
numbers = [ATOMIC_NUMBER[s.symbol] for s in structure.sites]
return lattice, positions, numbers
[docs]
def operations(
structure: PeriodicStructure, *, symprec: float = 1e-3
) -> list[tuple[np.ndarray, np.ndarray]]:
"""The space-group operations (R, t) in the structure's own fractional basis."""
if not 0 < symprec <= 1:
raise ValueError("symprec is a length tolerance between 0 and 1 angstrom")
dataset = _spglib().get_symmetry_dataset(_cell(structure), symprec=symprec)
if dataset is None:
raise ValueError("spglib found no space group at this tolerance")
return [(np.asarray(r, dtype=float), np.asarray(t, dtype=float) % 1.0)
for r, t in zip(dataset.rotations, dataset.translations, strict=True)] # fmt: skip
[docs]
def detect(structure: PeriodicStructure, *, symprec: float = 1e-3) -> dict[str, Any]:
"""The space group, setting, Wyckoff positions and equivalent atoms; read-only.
``symprec`` is a length tolerance in angstrom.
"""
if not 0 < symprec <= 1:
raise ValueError("symprec is a length tolerance between 0 and 1 angstrom")
spglib = _spglib()
dataset = spglib.get_symmetry_dataset(_cell(structure), symprec=symprec)
if dataset is None:
raise ValueError(f"spglib found no symmetry at symprec {symprec} A")
ids = [s.id for s in structure.sites]
sites = []
for i, sid in enumerate(ids):
sites.append(
{
"site": sid,
"symbol": structure.sites[i].symbol,
"wyckoff": str(dataset.wyckoffs[i]),
"site_symmetry": str(dataset.site_symmetry_symbols[i]).strip(),
"equivalent_to": ids[int(dataset.equivalent_atoms[i])],
}
)
orbits: dict[str, list[str]] = {}
for entry in sites:
orbits.setdefault(entry["equivalent_to"], []).append(entry["site"])
hall = spglib.get_spacegroup_type(int(dataset.hall_number))
return {
"number": int(dataset.number),
"international": str(dataset.international),
"hall": str(dataset.hall),
"hall_number": int(dataset.hall_number),
"choice": str(dataset.choice),
"point_group": str(dataset.pointgroup),
"crystal_system": _system(int(dataset.number)),
"operations": len(dataset.rotations),
"symprec_angstrom": symprec,
"sites": sites,
"independent_sites": len(orbits),
"standard_lattice": np.asarray(dataset.std_lattice).tolist(),
"transformation_matrix": np.asarray(dataset.transformation_matrix).tolist(),
"origin_shift": np.asarray(dataset.origin_shift).tolist(),
"arithmetic_crystal_class": str(hall.arithmetic_crystal_class_symbol),
"reference": REFERENCE,
}
def _system(number: int) -> str:
from vibesketcher.document.spacegroups import crystal_system
return crystal_system(number)
def _rebasis(
structure: PeriodicStructure, transform: np.ndarray, shift: np.ndarray, tolerance: float
) -> tuple[np.ndarray, list[tuple[str, np.ndarray]]]:
"""Measured atoms in the cell with rows ``transform @ L``, shifted by ``shift`` (fractional).
``transform`` may be rational (a larger or smaller cell); atoms are
enumerated in the new cell and the count must be N |det transform|.
"""
import itertools
lattice = crystal.lattice_rows(structure)
new = transform @ lattice
det = float(np.linalg.det(transform))
expected = len(structure.sites) * det
if det <= 0 or abs(expected - round(expected)) > 1e-6:
raise ValueError("the standard setting is not a lattice transform of this cell")
inverse = np.linalg.inv(new)
frac = np.asarray([s.fractional for s in structure.sites], dtype=float)
corners = np.asarray(list(itertools.product((0, 1), repeat=3)), dtype=float) @ transform
lo = np.floor(corners.min(axis=0)).astype(int) - 1
hi = np.ceil(corners.max(axis=0)).astype(int) + 1
shifts = np.asarray(list(itertools.product(*(range(lo[k], hi[k] + 1) for k in range(3)))))
eps = tolerance / float(np.linalg.norm(new, axis=1).min())
kept: list[tuple[str, np.ndarray]] = []
for site, f in zip(structure.sites, frac, strict=True):
for g in ((f + shifts) @ lattice) @ inverse + shift:
if np.any(g < -eps) or np.any(g >= 1 - eps):
continue
g = g - np.floor(g + eps)
if any(np.max(np.abs(((g - q) + 0.5) % 1 - 0.5)) < eps for _s, q in kept):
continue
kept.append((site.symbol, g))
if len(kept) != round(expected):
raise ValueError(
f"the standard cell holds {len(kept)} atoms, expected {round(expected)}"
) # pragma: no cover - guarded by spglib's own consistency
return new, kept
[docs]
def standardize(
structure: PeriodicStructure,
*,
identifier: str,
to: Literal["conventional", "primitive"] = "conventional",
idealize: bool = False,
symprec: float = 1e-3,
) -> PeriodicStructure:
"""The structure in the standard conventional or primitive setting of its space group.
Without ``idealize`` the measured atoms are only re-expressed in the
standard basis (spglib's transformation matrix and origin shift applied to
them, orientation kept); with it, spglib's symmetrised, rotated cell is
used. The record states how far idealisation moves atoms and cell vectors
either way, so the choice can be reviewed.
"""
spglib = _spglib()
cell = _cell(structure)
dataset = spglib.get_symmetry_dataset(cell, symprec=symprec)
ideal = spglib.standardize_cell(cell, to_primitive=False, no_idealize=False, symprec=symprec)
if dataset is None or ideal is None:
raise ValueError(f"spglib found no space group at symprec {symprec} A")
# spglib: (a_s b_s c_s) = (a b c) P^-1 and x_s = P x + p, so rows_s = P^-T rows.
p_matrix = np.asarray(dataset.transformation_matrix, dtype=float)
rows_transform = np.linalg.inv(p_matrix).T
lattice, sites = _rebasis(structure, rows_transform, np.asarray(dataset.origin_shift), symprec)
measured = (lattice, np.asarray([f for _s, f in sites]), [ATOMIC_NUMBER[s] for s, _f in sites])
displacement = _largest_displacement(measured, ideal)
if idealize:
lattice = np.asarray(ideal[0], dtype=float)
sites = [
(_symbol(int(z)), np.asarray(f, dtype=float))
for f, z in zip(ideal[1], ideal[2], strict=True)
]
number = int(dataset.number)
builder = crystal._builder(structure)
for key in ("primitive_transform", "conventional", "bravais"):
builder.pop(key, None)
builder.update(
{
"cell": f"standard_{to}",
"standardized_from": structure.id,
"space_group": {
"number": number,
"international": str(dataset.international),
"hall": str(dataset.hall),
},
"idealized": idealize,
"idealization_displacement_angstrom": displacement,
"symprec_angstrom": symprec,
"reference": REFERENCE,
}
)
pearson = _pearson(number, str(dataset.international))
if pearson:
builder["bravais"] = pearson
conventional = crystal._structure(
identifier, lattice, sites, (True, True, True), "conventional", builder
)
if to == "conventional":
result = conventional
elif pearson is None:
raise ValueError(f"no primitive transform is recorded for space group {number}")
elif len(crystal.centring_vectors(pearson)) == 1:
result = conventional.model_copy(update={"cell_setting": "primitive"})
else:
result = crystal.primitive_cell(conventional, pearson, identifier=identifier)
return result.model_copy(
update={
"charge": _scaled_charge(structure, len(result.sites)),
"multiplicity": structure.multiplicity
if len(result.sites) == len(structure.sites)
else None,
}
)
def _scaled_charge(structure: PeriodicStructure, count: int) -> int | None:
if structure.charge is None:
return None
ratio = count / len(structure.sites)
charge = structure.charge * ratio
return int(round(charge)) if abs(charge - round(charge)) < 1e-9 else None
def _symbol(number: int) -> str:
for symbol, z in ATOMIC_NUMBER.items():
if z == number:
return symbol
raise ValueError(f"unknown atomic number {number}") # pragma: no cover
def _pearson(number: int, international: str) -> str | None:
from vibesketcher.document.crystal import _pearson as pearson
from vibesketcher.document.spacegroups import find
try:
return pearson(find(number))
except ValueError: # pragma: no cover
return None
def _largest_displacement(kept: Any, ideal: Any) -> dict[str, float]:
"""How far idealisation moves atoms and cell vectors (angstrom).
Both cells are in the same standard basis, so fractional coordinates are
comparable; each measured atom is matched to the nearest idealised atom of
the same element (minimum image), since the atom order may differ, after
the best common translation (the origin may move to an equivalent atom).
"""
k_lattice, k_pos, k_num = kept
i_lattice, i_pos, i_num = ideal
k_pos, i_pos = np.asarray(k_pos), np.asarray(i_pos)
k_num, i_num = np.asarray(k_num), np.asarray(i_num)
lattice = np.asarray(i_lattice)
def worst(shift: np.ndarray) -> float:
moves = []
for f, z in zip(k_pos + shift, k_num, strict=True):
delta = i_pos[i_num == z] - f
delta -= np.round(delta)
moves.append(delta[int(np.argmin(np.linalg.norm(delta @ lattice, axis=1)))])
moved = np.asarray(moves)
moved -= moved.mean(axis=0) # a rigid translation is not a displacement
return float(np.max(np.linalg.norm(moved @ lattice, axis=1)))
# The idealised cell may sit on another origin: allow one common translation.
shifts = i_pos[i_num == k_num[0]] - k_pos[0]
largest = min(worst(t) for t in shifts)
vectors = float(
np.max(np.abs(np.linalg.norm(i_lattice, axis=1) - np.linalg.norm(k_lattice, axis=1)))
)
return {"atoms": largest, "cell_vectors": vectors}
[docs]
def reduce(
structure: PeriodicStructure, *, identifier: str, method: Literal["niggli", "delaunay"]
) -> PeriodicStructure:
"""The same crystal in its Niggli- or Delaunay-reduced cell (det P = +1, same atoms)."""
from vibesketcher.document.periodic import supercell
spglib = _spglib()
if sum(structure.pbc) != 3:
raise ValueError("cell reduction here needs a three-dimensional cell")
lattice = crystal.lattice_rows(structure)
reduced = (spglib.niggli_reduce if method == "niggli" else spglib.delaunay_reduce)(lattice)
if reduced is None:
raise ValueError(f"spglib could not {method}-reduce the cell") # pragma: no cover
p = np.rint(np.asarray(reduced) @ np.linalg.inv(lattice)).astype(int)
if round(np.linalg.det(p)) == -1:
p = -p # the same lattice, right-handed
if round(np.linalg.det(p)) != 1:
raise ValueError("the reduced cell is not a unimodular transform") # pragma: no cover
result = supercell(structure, p.tolist(), identifier=identifier)
builder = crystal._builder(structure)
builder.update(
{
"cell": f"{method}_reduced",
"reduction": {"method": method, "matrix": p.tolist()},
"reference": REFERENCE,
}
)
return result.model_copy(
update={"provenance": {**result.provenance, "builder": crystal._json(builder)}}
)