vibe-sketcherDocumentation
Source code for vibesketcher.document.crystal
"""Bravais lattices, prototype crystals and cell choices for periodic structures.
Cells are right-handed row vectors in angstrom, as in `PeriodicStructure`.
Conventional cells use the standard settings of the International Tables for
Crystallography, Vol. A (2016); primitive cells use the transformations of
W. Setyawan and S. Curtarolo, Comput. Mater. Sci. 49, 299 (2010),
doi:10.1016/j.commatsci.2010.05.010. Wigner-Seitz faces are the
Voronoi-relevant lattice vectors of J. H. Conway and N. J. A. Sloane, Sphere
Packings, Lattices and Groups, ch. 21 (Springer, 1999): within each coset
T + 2L, only the shortest pair +-T bounds the cell.
Every builder records its recipe in ``provenance["builder"]``; the document
schema is unchanged. A Wigner-Seitz cell is an atom placement of the same
lattice, not a different lattice.
"""
from __future__ import annotations
import itertools
import math
import re
from dataclasses import dataclass, field
from fractions import Fraction
from typing import Any, Literal
import numpy as np
from vibesketcher.document.model import PeriodicSite, PeriodicStructure
CellChoice = Literal["primitive", "conventional"]
#: Pearson-style symbols. Three-dimensional lattices use the 14 Bravais types;
#: plane lattices use the 5 two-dimensional types; the line lattice is ``l``.
BRAVAIS_3D = ("aP", "mP", "mS", "oP", "oS", "oI", "oF", "tP", "tI", "hR", "hP", "cP", "cI", "cF")
BRAVAIS_2D = ("mp", "op", "oc", "tp", "hp")
BRAVAIS_1D = ("l",)
BRAVAIS = BRAVAIS_3D + BRAVAIS_2D + BRAVAIS_1D
_BRAVAIS_NAMES = {
"aP": "triclinic",
"mP": "primitive monoclinic",
"mS": "base-centred monoclinic (C, unique axis b)",
"oP": "primitive orthorhombic",
"oS": "base-centred orthorhombic (C)",
"oI": "body-centred orthorhombic",
"oF": "face-centred orthorhombic",
"tP": "primitive tetragonal",
"tI": "body-centred tetragonal",
"hR": "rhombohedral (hexagonal axes, obverse)",
"hP": "hexagonal",
"cP": "simple cubic",
"cI": "body-centred cubic",
"cF": "face-centred cubic",
"mp": "oblique plane",
"op": "rectangular plane",
"oc": "centred rectangular plane",
"tp": "square plane",
"hp": "hexagonal plane",
"l": "line",
}
_H = Fraction(1, 2)
_T = Fraction(1, 3)
#: Lattice points of the conventional cell (fractional) and the integer-free
#: primitive transform P with primitive rows = P @ conventional rows.
_CENTRING: dict[str, tuple[Any, Any]] = {
"F": (
((0, 0, 0), (0, _H, _H), (_H, 0, _H), (_H, _H, 0)),
((0, _H, _H), (_H, 0, _H), (_H, _H, 0)),
),
"I": (
((0, 0, 0), (_H, _H, _H)),
((-_H, _H, _H), (_H, -_H, _H), (_H, _H, -_H)),
),
"C": (
((0, 0, 0), (_H, _H, 0)),
((_H, -_H, 0), (_H, _H, 0), (0, 0, 1)),
),
"R": (
((0, 0, 0), (2 * _T, _T, _T), (_T, 2 * _T, 2 * _T)),
((2 * _T, _T, _T), (-_T, _T, _T), (-_T, -2 * _T, _T)),
),
"P": (((0, 0, 0),), ((1, 0, 0), (0, 1, 0), (0, 0, 1))),
}
def _centring_of(bravais: str) -> str:
if bravais in ("oF", "cF"):
return "F"
if bravais in ("oI", "tI", "cI"):
return "I"
if bravais in ("oS", "mS", "oc"):
return "C"
if bravais == "hR":
return "R"
return "P"
[docs]
def dimensionality(bravais: str) -> int:
if bravais in BRAVAIS_3D:
return 3
if bravais in BRAVAIS_2D:
return 2
if bravais in BRAVAIS_1D:
return 1
raise ValueError(f"unknown Bravais lattice {bravais!r}; choose one of {', '.join(BRAVAIS)}")
[docs]
@dataclass(frozen=True)
class LatticeParameters:
"""Conventional-cell parameters in angstrom and degrees.
``vacuum`` is the length of each nonperiodic vector of a plane or line
lattice; it is part of the recorded model, never invented later.
"""
a: float
b: float | None = None
c: float | None = None
alpha: float | None = None
beta: float | None = None
gamma: float | None = None
vacuum: float = 20.0
def _resolve(bravais: str, p: LatticeParameters) -> tuple[float, float, float, float, float, float]:
"""Fill the parameters a lattice type fixes and refuse those it forbids."""
family = bravais[0]
a = p.a
free: dict[str, tuple[str, ...]] = {
"a": ("a", "b", "c", "alpha", "beta", "gamma"),
"m": ("a", "b", "c", "beta"),
"o": ("a", "b", "c"),
"t": ("a", "c"),
"h": ("a", "c"),
"c": ("a",),
}
plane: dict[str, tuple[str, ...]] = {
"mp": ("a", "b", "gamma"),
"op": ("a", "b"),
"oc": ("a", "b"),
"tp": ("a",),
"hp": ("a",),
"l": ("a",),
}
allowed = plane[bravais] if bravais in plane else free[family]
given = {k for k in ("b", "c", "alpha", "beta", "gamma") if getattr(p, k) is not None}
extra = sorted(given - set(allowed))
if extra:
raise ValueError(f"{bravais} ({_BRAVAIS_NAMES[bravais]}) fixes {', '.join(extra)}")
missing = [k for k in allowed if k != "a" and getattr(p, k) is None]
if missing:
raise ValueError(f"{bravais} needs {', '.join(missing)}")
values = [a, p.b, p.c, p.alpha, p.beta, p.gamma, p.vacuum]
if not all(v is None or (math.isfinite(v) and v > 0) for v in values):
raise ValueError("lattice parameters must be finite and positive")
b = p.b if p.b is not None else a
c = p.c if p.c is not None else a
alpha = p.alpha if p.alpha is not None else 90.0
beta = p.beta if p.beta is not None else 90.0
gamma = p.gamma if p.gamma is not None else 90.0
if bravais in ("hP", "hR", "hp"):
gamma = 120.0
if bravais in ("mp", "op", "oc", "tp", "hp"):
c, alpha, beta = p.vacuum, 90.0, 90.0
if bravais in ("tp", "hp"):
b = a
if bravais == "l":
b = c = p.vacuum
alpha = beta = gamma = 90.0
if bravais in ("mP", "mS") and abs(beta - 90.0) < 1e-9:
raise ValueError("monoclinic cells need beta different from 90 degrees")
for angle in (alpha, beta, gamma):
if not 0 < angle < 180:
raise ValueError("cell angles must lie strictly between 0 and 180 degrees")
return a, b, c, alpha, beta, gamma
[docs]
def cell_from_parameters(
a: float, b: float, c: float, alpha: float, beta: float, gamma: float
) -> np.ndarray:
"""Rows a1 along x, a2 in the xy plane (the IUCr standard orientation)."""
ca, cb, cg = (math.cos(math.radians(x)) for x in (alpha, beta, gamma))
sg = math.sin(math.radians(gamma))
cy = (ca - cb * cg) / sg
cz2 = 1 - cb * cb - cy * cy
if cz2 <= 1e-12:
raise ValueError("cell angles do not describe a three-dimensional cell")
rows = np.array([[a, 0, 0], [b * cg, b * sg, 0], [c * cb, c * cy, c * math.sqrt(cz2)]])
rows[np.abs(rows) < 1e-12] = 0.0
return rows
[docs]
def conventional_lattice(bravais: str, parameters: LatticeParameters) -> np.ndarray:
dimensionality(bravais)
return cell_from_parameters(*_resolve(bravais, parameters))
[docs]
def primitive_transform(bravais: str) -> np.ndarray:
"""P with primitive rows = P @ conventional rows (Setyawan-Curtarolo 2010)."""
dimensionality(bravais)
return np.array(_CENTRING[_centring_of(bravais)][1], dtype=float)
[docs]
def centring_vectors(bravais: str) -> np.ndarray:
dimensionality(bravais)
return np.array(_CENTRING[_centring_of(bravais)][0], dtype=float)
[docs]
def pbc_of(bravais: str) -> tuple[bool, bool, bool]:
d = dimensionality(bravais)
return (True, d >= 2, d >= 3)
_TERM = re.compile(r"([+-]?)(\d+(?:/\d+)?|\d*\.\d+)?\*?([xyz])?")
[docs]
def parse_operation(text: str) -> tuple[np.ndarray, np.ndarray]:
"""Parse a coordinate triplet such as ``-x+y,-x,z+1/2`` (ITA notation)."""
parts = text.replace(" ", "").lower().split(",")
if len(parts) != 3 or not all(parts):
raise ValueError(f"symmetry operation {text!r} needs three components")
rotation = np.zeros((3, 3))
translation = np.zeros(3)
for row, part in enumerate(parts):
position = 0
while position < len(part):
match = _TERM.match(part, position)
if match is None or match.end() == position or not (match[2] or match[3]):
raise ValueError(f"cannot read symmetry operation {text!r}")
sign = -1.0 if match[1] == "-" else 1.0
number = float(Fraction(match[2])) if match[2] else 1.0
if match[3]:
rotation[row, "xyz".index(match[3])] += sign * number
else:
translation[row] += sign * number
position = match.end()
if abs(abs(np.linalg.det(rotation)) - 1) > 1e-9:
raise ValueError(f"symmetry operation {text!r} is not an isometry")
return rotation, translation
def _unique_fractional(
points: list[tuple[str, np.ndarray]],
tolerance: float,
pbc: tuple[bool, bool, bool] = (True, True, True),
) -> list[tuple[str, np.ndarray]]:
"""Wrap periodic components into [0, 1) and merge coincident sites."""
mask = np.asarray(pbc, dtype=bool)
kept: list[tuple[str, np.ndarray]] = []
for symbol, frac in points:
wrapped = np.asarray(frac, dtype=float).copy()
wrapped[mask] -= np.floor(wrapped[mask] + 1e-12)
wrapped[mask & (np.abs(wrapped - 1) < 1e-12)] = 0.0
duplicate = None
for other_symbol, other in kept:
delta = wrapped - other
delta[mask] -= np.round(delta[mask])
if np.max(np.abs(delta)) < tolerance:
duplicate = other_symbol
break
if duplicate is None:
kept.append((symbol, wrapped))
elif duplicate != symbol:
raise ValueError("two different elements occupy the same site")
return kept
[docs]
@dataclass(frozen=True)
class Prototype:
"""A crystal prototype: lattice type, default parameters and Wyckoff data."""
name: str
description: str
bravais: str
parameters: LatticeParameters
sites: tuple[tuple[str, tuple[float, float, float]], ...]
operations: tuple[str, ...] = ()
reference: str = ""
species: tuple[str, ...] = field(default=())
_WYCKOFF = (
"R. W. G. Wyckoff, Crystal Structures, 2nd ed., vol. 1 (Interscience, 1963); "
"default lattice constants are room-temperature experimental values"
)
_R3C = (
"x,y,z",
"-y,x-y,z",
"-x+y,-x,z",
"y,x,-z+1/2",
"x-y,-y,-z+1/2",
"-x,-x+y,-z+1/2",
"-x,-y,-z",
"y,-x+y,-z",
"x-y,x,-z",
"-y,-x,z+1/2",
"-x+y,y,z+1/2",
"x,x-y,z+1/2",
)
PROTOTYPES: dict[str, Prototype] = {
p.name: p
for p in (
Prototype(
"rocksalt",
"rocksalt AB (Fm-3m), e.g. MgO, NaCl",
"cF",
LatticeParameters(4.212),
(("Mg", (0, 0, 0)), ("O", (0.5, 0.5, 0.5))),
reference=_WYCKOFF,
species=("Mg", "O"),
),
Prototype(
"corundum",
"corundum A2B3 (R-3c, hexagonal axes), e.g. alpha-Al2O3",
"hR",
LatticeParameters(4.7602, c=12.9933),
(("Al", (0, 0, 0.35216)), ("O", (0.30624, 0, 0.25))),
operations=_R3C,
reference=(
"J. Lewis, D. Schwarzenbach and H. D. Flack, Acta Cryst. A38, 733 (1982), "
"doi:10.1107/S0567739482001478"
),
species=("Al", "O"),
),
Prototype(
"fcc",
"face-centred cubic element (Fm-3m), e.g. Cu, Al, Au, Pt",
"cF",
LatticeParameters(3.615),
(("Cu", (0, 0, 0)),),
reference=_WYCKOFF,
species=("Cu",),
),
Prototype(
"bcc",
"body-centred cubic element (Im-3m), e.g. Fe, W, Na",
"cI",
LatticeParameters(2.8665),
(("Fe", (0, 0, 0)),),
reference=_WYCKOFF,
species=("Fe",),
),
Prototype(
"hcp",
"hexagonal close-packed element (P6_3/mmc), e.g. Mg, Ti, Zn",
"hP",
LatticeParameters(3.209, c=5.211),
(("Mg", (1 / 3, 2 / 3, 0.25)), ("Mg", (2 / 3, 1 / 3, 0.75))),
reference=_WYCKOFF,
species=("Mg",),
),
Prototype(
"diamond",
"diamond (Fd-3m), e.g. C, Si, Ge",
"cF",
LatticeParameters(5.431),
(("Si", (0, 0, 0)), ("Si", (0.25, 0.25, 0.25))),
reference=_WYCKOFF,
species=("Si",),
),
Prototype(
"zincblende",
"zincblende AB (F-43m), e.g. GaAs, ZnS",
"cF",
LatticeParameters(5.653),
(("Ga", (0, 0, 0)), ("As", (0.25, 0.25, 0.25))),
reference=_WYCKOFF,
species=("Ga", "As"),
),
Prototype(
"wurtzite",
"wurtzite AB (P6_3mc), e.g. ZnO, GaN",
"hP",
LatticeParameters(3.250, c=5.207),
(
("Zn", (1 / 3, 2 / 3, 0)),
("Zn", (2 / 3, 1 / 3, 0.5)),
("O", (1 / 3, 2 / 3, 0.382)),
("O", (2 / 3, 1 / 3, 0.882)),
),
reference=_WYCKOFF,
species=("Zn", "O"),
),
Prototype(
"cscl",
"caesium chloride AB (Pm-3m)",
"cP",
LatticeParameters(4.123),
(("Cs", (0, 0, 0)), ("Cl", (0.5, 0.5, 0.5))),
reference=_WYCKOFF,
species=("Cs", "Cl"),
),
Prototype(
"fluorite",
"fluorite AB2 (Fm-3m), e.g. CaF2, CeO2",
"cF",
LatticeParameters(5.463),
(("Ca", (0, 0, 0)), ("F", (0.25, 0.25, 0.25)), ("F", (0.75, 0.75, 0.75))),
reference=_WYCKOFF,
species=("Ca", "F"),
),
Prototype(
"rutile",
"rutile AB2 (P4_2/mnm), e.g. TiO2",
"tP",
LatticeParameters(4.594, c=2.959),
(
("Ti", (0, 0, 0)),
("Ti", (0.5, 0.5, 0.5)),
("O", (0.3048, 0.3048, 0)),
("O", (0.6952, 0.6952, 0)),
("O", (0.8048, 0.1952, 0.5)),
("O", (0.1952, 0.8048, 0.5)),
),
reference=_WYCKOFF,
species=("Ti", "O"),
),
Prototype(
"perovskite",
"cubic perovskite ABX3 (Pm-3m), e.g. SrTiO3",
"cP",
LatticeParameters(3.905),
(
("Sr", (0, 0, 0)),
("Ti", (0.5, 0.5, 0.5)),
("O", (0.5, 0.5, 0)),
("O", (0.5, 0, 0.5)),
("O", (0, 0.5, 0.5)),
),
reference=_WYCKOFF,
species=("Sr", "Ti", "O"),
),
Prototype(
"graphene",
"graphene sheet (hexagonal plane)",
"hp",
LatticeParameters(2.46),
(("C", (1 / 3, 2 / 3, 0)), ("C", (2 / 3, 1 / 3, 0))),
reference=_WYCKOFF,
species=("C",),
),
Prototype(
"hbn",
"hexagonal boron nitride sheet (hexagonal plane)",
"hp",
LatticeParameters(2.504),
(("B", (1 / 3, 2 / 3, 0)), ("N", (2 / 3, 1 / 3, 0))),
reference=_WYCKOFF,
species=("B", "N"),
),
Prototype(
"chain",
"equally spaced atomic chain (line lattice), e.g. the H-chain model",
"l",
LatticeParameters(2.0),
(("H", (0, 0, 0)), ("H", (0.5, 0, 0))),
reference="model system; choose the spacing for the study",
species=("H",),
),
)
}
[docs]
def prototypes() -> list[dict[str, Any]]:
"""The prototype library as plain data (for the CLI, REST and MCP)."""
return [
{
"name": p.name,
"description": p.description,
"bravais": p.bravais,
"lattice": _BRAVAIS_NAMES[p.bravais],
"species": list(p.species),
"parameters": {
k: v for k, v in vars(p.parameters).items() if v is not None and k != "vacuum"
},
"reference": p.reference,
}
for p in PROTOTYPES.values()
]
def _low_dimensional_fraction(bravais: str, frac: tuple[float, float, float]) -> np.ndarray:
"""Plane/line bases give out-of-plane offsets as fractions of the vacuum vector.
The layer sits at the middle of each nonperiodic vector.
"""
out = np.asarray(frac, dtype=float)
d = dimensionality(bravais)
for axis in range(d, 3):
out[axis] += 0.5
return out
[docs]
def build_crystal(
bravais: str,
parameters: LatticeParameters,
sites: list[tuple[str, tuple[float, float, float]]],
*,
identifier: str,
operations: list[str] | tuple[str, ...] = (),
cell: CellChoice = "conventional",
tolerance: float = 1e-4,
species_map: dict[str, str] | None = None,
recipe: dict[str, Any] | None = None,
) -> PeriodicStructure:
"""Expand Wyckoff ``sites`` (conventional fractional) into a physical cell.
``operations`` are ITA coordinate triplets without the centring
translations, which the Bravais type adds. Without operations the sites are
taken as the full P1 motif of one lattice point. ``species_map`` renames
elements of a prototype (``{"Mg": "Ca"}`` turns MgO into CaO).
"""
if not sites:
raise ValueError("a crystal needs at least one site")
lattice = conventional_lattice(bravais, parameters)
ops = [parse_operation(op) for op in (operations or ("x,y,z",))]
centring = centring_vectors(bravais)
mapping = species_map or {}
points: list[tuple[str, np.ndarray]] = []
for symbol, frac in sites:
symbol = mapping.get(symbol, symbol)
position = np.asarray(frac, dtype=float)
if position.shape != (3,) or not np.isfinite(position).all():
raise ValueError("site coordinates must be three finite fractions")
for rotation, translation in ops:
image = rotation @ position + translation
for shift in centring:
points.append((symbol, image + shift))
d = dimensionality(bravais)
if d < 3 and operations:
raise ValueError("symmetry operations apply to three-dimensional lattices only")
unique = _unique_fractional(points, tolerance, pbc_of(bravais))
if d < 3:
for _, offset in unique:
if np.any(np.abs(offset[d:]) >= 0.5):
raise ValueError("out-of-plane offsets must stay within half the vacuum vector")
unique = [(s, _low_dimensional_fraction(bravais, tuple(f))) for s, f in unique]
if len(unique) > 25000:
raise ValueError("crystal exceeds 25000 sites")
record: dict[str, Any] = {
"kind": "crystal",
"bravais": bravais,
"lattice": _BRAVAIS_NAMES[bravais],
"parameters": {k: v for k, v in vars(parameters).items() if v is not None},
"cell": cell,
**(recipe or {}),
}
conventional = _structure(
identifier,
lattice,
[(s, f) for s, f in unique],
pbc_of(bravais),
"conventional",
record | {"cell": "conventional"},
)
if cell == "conventional":
return conventional
return primitive_cell(conventional, bravais, identifier=identifier, recipe=record)
def _site_ids(symbols: list[str]) -> list[str]:
counts: dict[str, int] = {}
ids = []
for symbol in symbols:
counts[symbol] = counts.get(symbol, 0) + 1
ids.append(f"{symbol}{counts[symbol]}")
return ids
def _structure(
identifier: str,
lattice: np.ndarray,
sites: list[tuple[str, np.ndarray]],
pbc: tuple[bool, bool, bool],
setting: Literal["primitive", "conventional", "unspecified"],
builder: dict[str, Any],
*,
charge: int | None = None,
multiplicity: int | None = None,
site_ids: list[str] | None = None,
) -> PeriodicStructure:
ids = site_ids or _site_ids([s for s, _ in sites])
return PeriodicStructure(
id=identifier,
lattice_angstrom=tuple(tuple(float(x) for x in row) for row in lattice),
pbc=pbc,
sites=[
PeriodicSite(id=i, symbol=s, fractional=tuple(_clean(x) for x in f))
for i, (s, f) in zip(ids, sites, strict=True)
],
cell_setting=setting,
charge=charge,
multiplicity=multiplicity,
provenance={"builder": _json(builder)},
)
def _clean(x: float) -> float:
value = float(x)
rounded = round(value, 12)
return 0.0 if rounded == 0 else rounded
def _json(value: Any) -> Any:
if isinstance(value, dict):
return {str(k): _json(v) for k, v in value.items()}
if isinstance(value, (list, tuple)):
return [_json(v) for v in value]
if isinstance(value, np.ndarray):
return _json(value.tolist())
if isinstance(value, (np.floating, np.integer)):
return value.item()
return value
def _builder(structure: PeriodicStructure) -> dict[str, Any]:
recorded = structure.provenance.get("builder")
return dict(recorded) if isinstance(recorded, dict) else {}
[docs]
def lattice_rows(structure: PeriodicStructure) -> np.ndarray:
return np.asarray(structure.lattice_angstrom, dtype=float)
[docs]
def cartesian(structure: PeriodicStructure) -> np.ndarray:
fractional = np.asarray([s.fractional for s in structure.sites], dtype=float)
return np.asarray(fractional @ lattice_rows(structure))
[docs]
def primitive_cell(
structure: PeriodicStructure,
bravais: str,
*,
identifier: str,
tolerance: float = 1e-4,
recipe: dict[str, Any] | None = None,
) -> PeriodicStructure:
"""Reduce a conventional cell of the given Bravais type to its primitive cell.
The site count must drop by exactly the number of lattice points of the
conventional cell, which proves that the motif carries the centring.
"""
transform = primitive_transform(bravais)
conventional = lattice_rows(structure)
lattice = transform @ conventional
inverse = np.linalg.inv(lattice)
points = [
(site.symbol, np.asarray(structure.cartesian(site)) @ inverse) for site in structure.sites
]
unique = _unique_fractional(points, tolerance, structure.pbc)
expected = len(structure.sites) / len(centring_vectors(bravais))
if len(unique) != expected:
raise ValueError(
f"the conventional cell does not carry the {bravais} centring "
f"({len(structure.sites)} sites, {len(unique)} after reduction)"
)
builder = _builder(structure)
builder.update(recipe or {})
builder.update(
{"cell": "primitive", "primitive_transform": transform, "conventional": conventional}
)
return _structure(
identifier,
lattice,
unique,
structure.pbc,
"primitive",
builder,
charge=structure.charge,
multiplicity=structure.multiplicity,
)
[docs]
def prototype(
name: str,
*,
identifier: str,
cell: CellChoice = "primitive",
a: float | None = None,
c: float | None = None,
vacuum: float | None = None,
species_map: dict[str, str] | None = None,
) -> PeriodicStructure:
"""Build a library prototype, optionally with other lattice constants or elements."""
try:
proto = PROTOTYPES[name.lower()]
except KeyError:
raise ValueError(
f"unknown prototype {name!r}; choose one of {', '.join(PROTOTYPES)}"
) from None
p = proto.parameters
parameters = LatticeParameters(
a if a is not None else p.a,
p.b,
c if c is not None else p.c,
p.alpha,
p.beta,
p.gamma,
vacuum if vacuum is not None else p.vacuum,
)
return build_crystal(
proto.bravais,
parameters,
list(proto.sites),
identifier=identifier,
operations=proto.operations,
cell=cell,
species_map=species_map,
recipe={"prototype": proto.name, "reference": proto.reference},
)
def _search_ranges(lattice: np.ndarray, pbc: tuple[bool, bool, bool], radius: float) -> list[range]:
"""Integer coefficient ranges that contain every translation shorter than ``radius``."""
reciprocal = np.linalg.inv(lattice).T # rows b_i with a_i . b_j = delta_ij
ranges = []
for axis in range(3):
if not pbc[axis]:
ranges.append(range(0, 1))
continue
n = int(math.ceil(radius * float(np.linalg.norm(reciprocal[axis])))) + 1
ranges.append(range(-n, n + 1))
if math.prod(len(r) for r in ranges) > 200000:
raise ValueError("the cell is too skewed for a bounded Wigner-Seitz search")
return ranges
[docs]
def wigner_seitz_cell(
structure: PeriodicStructure,
*,
identifier: str,
tolerance: float = 1e-6,
) -> PeriodicStructure:
"""Move each site to its lattice image nearest the origin lattice point.
The lattice, composition and site order are unchanged, so a PBC job sees
the same model; a cyclic cluster takes its shape from this placement.
Sites equidistant from two lattice points lie on a Wigner-Seitz face; the
image with the lexicographically smallest translation is kept and the site
is listed in ``provenance["builder"]["boundary_sites"]``.
"""
lattice = lattice_rows(structure)
periodic = [i for i in range(3) if structure.pbc[i]]
corner = np.abs(lattice[periodic]).sum(axis=0)
radius = float(np.linalg.norm(corner)) * 1.01 + 1.0
ranges = _search_ranges(lattice, structure.pbc, radius)
shifts = np.asarray(list(itertools.product(*ranges)), dtype=float)
boundary: list[str] = []
placed = []
for site in structure.sites:
frac = np.asarray(site.fractional, dtype=float)
base = frac.copy()
for axis in periodic:
base[axis] -= math.floor(base[axis] + 0.5)
candidates = (base + shifts) @ lattice
# Only periodic components measure the distance to the lattice point.
planar = candidates.copy()
if len(periodic) < 3:
projector = _periodic_projector(lattice, structure.pbc)
planar = candidates @ projector
distance = np.linalg.norm(planar, axis=1)
best = float(distance.min())
ties = np.flatnonzero(distance <= best + tolerance)
choice = min(ties, key=lambda i: tuple(-shifts[i]))
if len(ties) > 1:
boundary.append(site.id)
placed.append((site, base + shifts[choice]))
builder = _builder(structure)
builder.update({"placement": "wigner_seitz", "boundary_sites": boundary})
if builder.get("parent_cell") is None and structure.parent:
builder["parent_cell"] = structure.parent
data = structure.model_dump()
data["id"] = identifier
data["sites"] = [
{**s.model_dump(), "fractional": tuple(_clean(x) for x in f)} for s, f in placed
]
data["provenance"] = {**structure.provenance, "builder": _json(builder)}
data["source_cif"] = None
data["source_model_sha256"] = None
return PeriodicStructure.model_validate(data)
def _periodic_projector(lattice: np.ndarray, pbc: tuple[bool, bool, bool]) -> np.ndarray:
"""Orthogonal projector (row-vector convention) onto the span of periodic vectors."""
rows = lattice[[i for i in range(3) if pbc[i]]]
q, _ = np.linalg.qr(rows.T)
return np.asarray(q @ q.T)
[docs]
def wigner_seitz_polyhedron(structure: PeriodicStructure) -> dict[str, Any]:
"""Vertices (angstrom) and faces of the Wigner-Seitz cell of the periodic lattice.
3-D lattices give a polyhedron with faces listed as vertex cycles, plane
lattices a polygon in the periodic plane and the line lattice a segment.
Faces are the Voronoi-relevant vectors (Conway and Sloane, ch. 21).
"""
lattice = lattice_rows(structure)
periodic = [i for i in range(3) if structure.pbc[i]]
rows = lattice[periodic]
d = len(periodic)
if d == 1:
half = rows[0] / 2
return {"dimension": 1, "vertices": [(-half).tolist(), half.tolist()], "faces": []}
# Coefficients up to 2 on an LLL-free basis are enough for the coset minima
# of reasonably reduced cells; widen the search for skewed ones.
radius = 2.0 * float(np.linalg.norm(np.abs(rows).sum(axis=0)))
reciprocal = np.linalg.pinv(rows).T
spans = [int(math.ceil(radius * float(np.linalg.norm(reciprocal[k])))) + 1 for k in range(d)]
if math.prod(2 * s + 1 for s in spans) > 200000:
raise ValueError("the cell is too skewed for a bounded Wigner-Seitz construction")
best: dict[tuple[int, ...], list[tuple[float, np.ndarray]]] = {}
for n in itertools.product(*(range(-s, s + 1) for s in spans)):
if not any(n):
continue
parity = tuple(int(x) % 2 for x in n)
vector = np.asarray(n, dtype=float) @ rows
length = float(vector @ vector)
best.setdefault(parity, []).append((length, vector))
normals = []
for candidates in best.values():
candidates.sort(key=lambda item: item[0])
shortest = [v for length, v in candidates if length <= candidates[0][0] * (1 + 1e-9)]
if len(shortest) == 2: # +-T only: a facet
normals.extend(shortest)
planes = np.asarray(normals)
offsets = 0.5 * np.einsum("ij,ij->i", planes, planes)
vertices: list[np.ndarray] = []
basis = None
if d == 2:
q, _ = np.linalg.qr(rows.T)
basis = q[:, :2] # columns spanning the plane
groups = itertools.combinations(range(len(planes)), d)
for combo in groups:
idx = list(combo)
if d == 3:
matrix = planes[idx]
if abs(np.linalg.det(matrix)) < 1e-10:
continue
point = np.linalg.solve(matrix, offsets[idx])
else:
assert basis is not None
matrix = planes[idx] @ basis
if abs(np.linalg.det(matrix)) < 1e-10:
continue
point = basis @ np.linalg.solve(matrix, offsets[idx])
if np.all(planes @ point <= offsets + 1e-8) and not any(
np.linalg.norm(point - v) < 1e-7 for v in vertices
):
vertices.append(point)
faces = []
for k, normal in enumerate(planes if d == 3 else []):
on = [i for i, v in enumerate(vertices) if abs(normal @ v - offsets[k]) < 1e-7]
if len(on) < d:
continue
centre = np.mean([vertices[i] for i in on], axis=0)
u = vertices[on[0]] - centre
u /= np.linalg.norm(u)
w = np.cross(normal / np.linalg.norm(normal), u)
on.sort(
key=lambda i: math.atan2(
float((vertices[i] - centre) @ w), float((vertices[i] - centre) @ u)
)
)
faces.append(on)
if d == 2:
centre = np.zeros(3)
assert basis is not None
order = sorted(
range(len(vertices)),
key=lambda i: math.atan2(
float(vertices[i] @ basis[:, 1]), float(vertices[i] @ basis[:, 0])
),
)
vertices = [vertices[i] for i in order]
faces = [list(range(len(vertices)))]
volume = (
abs(float(np.linalg.det(rows)))
if d == 3
else float(np.linalg.norm(np.cross(rows[0], rows[1])))
)
return {
"dimension": d,
"vertices": [[_clean(x) for x in v] for v in vertices],
"faces": faces,
"neighbours": [[_clean(x) for x in v] for v in planes],
"measure": volume,
}
# Lattice constraints of each crystal system for space-group builds (#199).
_FIXED_ANGLES = {"orthorhombic", "tetragonal", "cubic"}
def _space_group_parameters(
setting: object, parameters: LatticeParameters
) -> tuple[float, float, float, float, float, float]:
"""Fill and check a, b, c, alpha, beta, gamma against the setting's crystal system."""
from vibesketcher.document.spacegroups import Setting
assert isinstance(setting, Setting)
system = setting.crystal_system
p = parameters
a = p.a
b = p.b if p.b is not None else (a if system in ("tetragonal", "trigonal", "hexagonal",
"cubic") else None) # fmt: skip
c = p.c if p.c is not None else (a if system == "cubic" or setting.extension == "R" else None)
if setting.extension == "R":
b = a if p.b is None else p.b
alpha = p.alpha
if alpha is None:
raise ValueError(f"{setting.symbol} on rhombohedral axes needs alpha")
angles = (alpha, p.beta or alpha, p.gamma or alpha)
if not (a == b == c and angles[0] == angles[1] == angles[2]):
raise ValueError("rhombohedral axes need a = b = c and alpha = beta = gamma")
return a, b, c, *angles
if b is None or c is None:
raise ValueError(f"{setting.symbol} ({system}) needs b and c" if b is None
else f"{setting.symbol} ({system}) needs c") # fmt: skip
alpha = p.alpha if p.alpha is not None else 90.0
beta = p.beta if p.beta is not None else 90.0
gamma = p.gamma if p.gamma is not None else (120.0 if system in ("trigonal", "hexagonal")
else 90.0) # fmt: skip
angles = (alpha, beta, gamma)
tolerance = 1e-6
if system in _FIXED_ANGLES and any(abs(x - 90) > tolerance for x in angles):
raise ValueError(f"{system} cells have alpha = beta = gamma = 90 degrees")
if system in ("tetragonal", "cubic", "trigonal", "hexagonal") and abs(a - b) > tolerance:
raise ValueError(f"{system} cells have a = b")
if system == "cubic" and abs(a - c) > tolerance:
raise ValueError("cubic cells have a = b = c")
if system in ("trigonal", "hexagonal") and (abs(alpha - 90) > tolerance
or abs(beta - 90) > tolerance
or abs(gamma - 120) > tolerance): # fmt: skip
raise ValueError(f"{system} cells on hexagonal axes have alpha = beta = 90, gamma = 120")
if system == "monoclinic":
unique = setting.qualifier[:1] or "b"
free = {"a": alpha, "b": beta, "c": gamma}[unique]
others = [x for name, x in zip("abc", angles, strict=True) if name != unique]
if any(abs(x - 90) > tolerance for x in others) or abs(free - 90) < tolerance:
raise ValueError(
f"{setting.symbol} has unique axis {unique}: only "
f"{ {'a': 'alpha', 'b': 'beta', 'c': 'gamma'}[unique] } differs from 90 degrees"
)
return a, b, c, alpha, beta, gamma
def _pearson(setting: object) -> str | None:
"""The Bravais symbol this module's primitive transforms know, if any."""
from vibesketcher.document.spacegroups import Setting
assert isinstance(setting, Setting)
system, centring = setting.crystal_system, setting.centring
if setting.extension == "R":
return None
table = {
("triclinic", "P"): "aP", ("monoclinic", "P"): "mP", ("orthorhombic", "P"): "oP",
("orthorhombic", "C"): "oS", ("orthorhombic", "I"): "oI", ("orthorhombic", "F"): "oF",
("tetragonal", "P"): "tP", ("tetragonal", "I"): "tI", ("trigonal", "P"): "hP",
("trigonal", "R"): "hR", ("hexagonal", "P"): "hP", ("cubic", "P"): "cP",
("cubic", "I"): "cI", ("cubic", "F"): "cF",
} # fmt: skip
if system == "monoclinic" and centring == "C" and (setting.qualifier[:1] or "b") == "b":
return "mS"
return table.get((system, centring))
[docs]
def from_space_group(
symbol: str | int,
parameters: LatticeParameters,
sites: list[tuple[str, tuple[float, float, float]]],
*,
identifier: str,
cell: CellChoice = "conventional",
tolerance: float = 1e-4,
near_special: float = 0.1,
species_map: dict[str, str] | None = None,
) -> PeriodicStructure:
"""Expand an asymmetric unit with the operations of a space-group setting (#199).
``sites`` are fractional coordinates of the setting's conventional cell.
Every generated site records its orbit (the asymmetric site it comes from)
and the orbit's multiplicity; a site that comes within ``near_special``
angstrom of a special position without lying on it is reported, never
moved (no snapping). Two elements on one position are refused.
"""
from vibesketcher.document.spacegroups import find
setting = find(symbol)
operations = setting.operations()
a, b, c, alpha, beta, gamma = _space_group_parameters(setting, parameters)
lattice = cell_from_parameters(a, b, c, alpha, beta, gamma)
mapping = species_map or {}
generated: list[tuple[str, np.ndarray, str]] = []
orbits: list[dict[str, object]] = []
warnings: list[str] = []
labels = [f"{mapping.get(s, s)}{i + 1}" for i, (s, _p) in enumerate(sites)]
for (symbol_in, frac), label in zip(sites, labels, strict=True):
element = mapping.get(symbol_in, symbol_in)
point = np.asarray(frac, dtype=float)
if point.shape != (3,) or not np.isfinite(point).all():
raise ValueError("site coordinates must be three finite fractions")
images = [np.asarray(op.apply(tuple(point))) for op in operations]
unique = _unique_fractional([(element, img) for img in images], tolerance)
# Distance of each image from the site itself (minimum image), in angstrom.
for img in images:
delta = img - point
delta -= np.round(delta)
distance = float(np.linalg.norm(delta @ lattice))
if tolerance < float(np.max(np.abs(delta))) and distance < near_special:
warnings.append(
f"{label} lies {distance:.3f} A from a special position; not snapped"
)
break
orbits.append({"site": label, "element": element, "fractional": list(point),
"multiplicity": len(unique),
"site_symmetry_order": len(operations) // len(unique)}) # fmt: skip
generated.extend((element, f, label) for _e, f in unique)
merged = _unique_fractional([(e, f) for e, f, _l in generated], tolerance)
if len(merged) != len(generated):
raise ValueError("two asymmetric sites generate the same position")
if len(generated) > 25000:
raise ValueError("the expanded crystal exceeds 25000 sites")
counts: dict[str, int] = {}
ids = []
for element, _f, _label in generated:
counts[element] = counts.get(element, 0) + 1
ids.append(f"{element}{counts[element]}")
orbit_of = {sid: label for sid, (_e, _f, label) in zip(ids, generated, strict=True)}
pearson = _pearson(setting)
builder: dict[str, object] = {
"kind": "crystal",
"space_group": {"number": setting.number, "hm": setting.hm,
"extension": setting.extension, "hall": setting.hall,
"operations": len(operations)},
"crystal_system": setting.crystal_system,
"parameters": {"a": a, "b": b, "c": c, "alpha": alpha, "beta": beta, "gamma": gamma},
"asymmetric_unit": orbits,
"orbit_of_site": orbit_of,
"tolerance": tolerance,
"warnings": warnings,
"cell": "conventional",
"reference": "S. R. Hall, Acta Cryst. A37, 517 (1981)",
} # fmt: skip
if pearson:
builder["bravais"] = pearson
structure = _structure(
identifier, lattice, [(e, f) for e, f, _l in generated], (True, True, True),
"conventional", builder, site_ids=ids,
) # fmt: skip
if cell == "conventional":
return structure
if pearson is None:
raise ValueError(f"no primitive transform is recorded for {setting.symbol}")
if len(centring_vectors(pearson)) == 1:
return structure.model_copy(update={"cell_setting": "primitive"})
return primitive_cell(structure, pearson, identifier=identifier)