vibe-sketcherDocumentation
Source code for vibesketcher.document.bilbao
"""Bilbao Crystallographic Server structure format and setting transformations (#230).
The BCS structure format, as shown on the server's input pages (TRANSTRU,
VISUALIZE): the space-group ITA number, the lattice parameters (angstrom,
degrees), the number of independent atoms, and one line per atom of the
asymmetric unit::
# Space Group ITA number
221
# Lattice parameters (lengths in angstroems and angles in degrees)
5.0 5.0 5.0 90 90 90
# Number of independent atoms in the asymmetric unit
3
# [atom type] [number] [WP] [x] [y] [z]
Ba 1 1a 0.0 0.0 0
Ti 2 1b 0.5 0.5 0.5
O 3 3c 0.5 0.0 0.5
Lines starting with ``#`` are labels. The server takes the default setting
of the space group, which is also what :func:`spacegroups.find` chooses for a
bare number (origin choice 2, hexagonal axes). The second atom column is kept
as a label only; a Wyckoff field such as ``3c`` (or ``-``) is checked against
the multiplicity of the generated orbit. Nothing is retrieved from the
server: users paste or open files they saved themselves.
Setting transformations ``(P, p)`` follow International Tables for
Crystallography Vol. A, Section 1.5: the new basis is (a', b', c') =
(a, b, c) P, the origin shift p is given in the old coordinates, and
x' = P^-1 (x - p). Written as ``a-b,a+b,c;0,0,1/4``.
References: M. I. Aroyo et al., Z. Kristallogr. 221, 15 (2006);
M. I. Aroyo et al., Bulg. Chem. Commun. 43, 183 (2011).
"""
from __future__ import annotations
import re
from collections import Counter
from fractions import Fraction
from typing import Any
import numpy as np
from vibesketcher.document import crystal
from vibesketcher.document.model import PeriodicStructure
from vibesketcher.document.spacegroups import find
from vibesketcher.elements import SYMBOLS
REFERENCE = (
"M. I. Aroyo, J. M. Perez-Mato, C. Capillas, E. Kroumova, S. Ivantchev, G. Madariaga, "
"A. Kirov, H. Wondratschek, Z. Kristallogr. 221, 15 (2006)"
)
_ITA = "International Tables for Crystallography Vol. A, Section 1.5"
def _element(label: str, line: int) -> str:
match = re.match(r"([A-Za-z]{1,2})", label)
if match:
for candidate in (match[1][:2], match[1][:1]):
symbol = candidate[0].upper() + candidate[1:].lower()
if symbol in SYMBOLS:
return symbol
raise ValueError(f"BCS line {line}: {label!r} is not an element")
def _fraction(text: str, line: int) -> float:
try:
return float(Fraction(text))
except (ValueError, ZeroDivisionError):
raise ValueError(f"BCS line {line}: {text!r} is not a number") from None
[docs]
def read_bcs(text: str, *, identifier: str) -> PeriodicStructure:
"""A crystal from BCS structure-format text, in the default setting."""
from vibesketcher.limits import check_size
check_size(text, "BCS structure")
lines = [(n + 1, raw.split("#", 1)[0].split())
for n, raw in enumerate(text.replace("\r", "").split("\n"))] # fmt: skip
lines = [(n, words) for n, words in lines if words]
if len(lines) < 4:
raise ValueError("BCS structure needs the group, the cell, the atom count and atoms")
(n_group, group), (n_cell, cell), (n_count, count) = lines[0], lines[1], lines[2]
if len(group) != 1 or not group[0].isdigit() or not 1 <= int(group[0]) <= 230:
raise ValueError(f"BCS line {n_group}: the space-group ITA number must be 1 to 230")
if len(cell) != 6:
raise ValueError(f"BCS line {n_cell}: give a b c alpha beta gamma")
a, b, c, alpha, beta, gamma = (_fraction(x, n_cell) for x in cell)
if len(count) != 1 or not count[0].isdigit():
raise ValueError(f"BCS line {n_count}: the number of independent atoms is missing")
atoms = lines[3:]
if len(atoms) != int(count[0]):
raise ValueError(f"BCS structure announces {count[0]} atoms but lists {len(atoms)}")
setting = find(int(group[0]))
sites: list[tuple[str, tuple[float, float, float]]] = []
wyckoff: list[tuple[int, str, str]] = []
labels = []
for n, words in atoms:
if len(words) != 6:
raise ValueError(f"BCS line {n}: give atom type, number, Wyckoff position, x y z")
label, index, position, x, y, z = words
sites.append((_element(label, n), (_fraction(x, n), _fraction(y, n), _fraction(z, n))))
wyckoff.append((n, position, label))
labels.append({"type": label, "number": index, "wyckoff": position})
structure = crystal.from_space_group(
setting.number, crystal.LatticeParameters(a, b, c, alpha, beta, gamma), sites,
identifier=identifier, cell="conventional",
) # fmt: skip
builder = dict(crystal._builder(structure))
orbits = builder.get("asymmetric_unit", [])
assert isinstance(orbits, list)
for (n, position, label), orbit in zip(wyckoff, orbits, strict=True):
match = re.fullmatch(r"(\d+)([a-zA-Z])", position)
if position != "-" and match is None:
raise ValueError(f"BCS line {n}: Wyckoff position {position!r} is not like 4a")
if match and int(match[1]) != orbit["multiplicity"]:
raise ValueError(
f"BCS line {n}: {label} is given as {position} but its orbit has "
f"{orbit['multiplicity']} sites in {setting.symbol}"
)
builder["bcs"] = {"atoms": labels, "setting": "default", "reference": REFERENCE}
return structure.model_copy(
update={"provenance": {**structure.provenance, "format": "bcs", "builder": builder}}
)
[docs]
def write_bcs(structure: PeriodicStructure, *, symprec: float = 1e-3) -> str:
"""BCS structure-format text of a 3D crystal.
A structure built from a space group in its default setting (a template,
the space-group builder or a BCS file) is written from its recorded
asymmetric unit. Anything else is standardised and idealised with spglib
(the ``symmetry`` extra) first.
"""
if sum(structure.pbc) != 3:
raise ValueError("the BCS structure format describes 3D crystals")
record = crystal._builder(structure)
group = record.get("space_group")
orbits = record.get("asymmetric_unit")
usable = (
isinstance(group, dict) and isinstance(orbits, list) and structure.cell_setting ==
"conventional" and record.get("kind") in ("crystal", None) and
find(int(group["number"])).hall == group.get("hall")
) # fmt: skip
if usable:
assert isinstance(group, dict) and isinstance(orbits, list)
number = int(group["number"])
params = record["parameters"]
assert isinstance(params, dict)
cell = [params[k] for k in ("a", "b", "c", "alpha", "beta", "gamma")]
letters = _wyckoff_letters(structure, record, len(orbits))
atoms = [(o["element"], int(o["multiplicity"]), letter, list(o["fractional"]))
for o, letter in zip(orbits, letters, strict=True)] # fmt: skip
else:
from vibesketcher.document import symmetry
standard = symmetry.standardize(structure, identifier="bcs", idealize=True,
symprec=symprec) # fmt: skip
report = symmetry.detect(standard, symprec=symprec)
default = find(int(report["number"]))
if default.hall.replace(" ", "") != str(report["hall"]).replace(" ", ""):
shift = origin_shift(find(str(report["hall"])), default)
if shift is None:
raise ValueError(
f"spglib's standard setting of {report['international']} is not related "
"to the BCS default setting by an origin shift"
)
moved = ",".join(str(Fraction(-x).limit_denominator(8)) for x in shift)
standard = transform_setting(standard, f"a,b,c;{moved}", identifier="bcs")
report = symmetry.detect(standard, symprec=symprec)
number = int(report["number"])
cell = list(crystal_parameters(standard))
atoms = []
seen = set()
for site, info in zip(standard.sites, report["sites"], strict=True):
if info["equivalent_to"] in seen:
continue
seen.add(info["equivalent_to"])
multiplicity = sum(
1 for s in report["sites"] if s["equivalent_to"] == info["equivalent_to"]
)
atoms.append((site.symbol, multiplicity, info["wyckoff"], list(site.fractional)))
lines = ["# Space Group ITA number", str(number),
"# Lattice parameters (lengths in angstroems and angles in degrees)",
" ".join(f"{x:.6g}" for x in cell),
"# Number of independent atoms in the asymmetric unit", str(len(atoms)),
"# [atom type] [number] [WP] [x] [y] [z]"] # fmt: skip
counts: Counter[str] = Counter()
for index, (element, multiplicity, letter, frac) in enumerate(atoms, start=1):
counts[element] += 1
position = f"{multiplicity}{letter}" if letter else "-"
coords = " ".join(f"{x % 1.0:.6f}".rstrip("0").rstrip(".") or "0" for x in frac)
lines.append(f"{element}{counts[element]} {index} {position} {coords}")
return "\n".join(lines) + "\n"
def _wyckoff_letters(structure: PeriodicStructure, record: dict[str, Any], count: int) -> list[str]:
"""Letters from a BCS input, else from spglib when installed, else none ("-")."""
given = record.get("bcs")
if isinstance(given, dict):
letters = [re.sub(r"^\d+", "", str(a.get("wyckoff", ""))).replace("-", "")
for a in given.get("atoms", [])] # fmt: skip
if len(letters) == count:
return letters
from vibesketcher.document import symmetry
orbit_of = record.get("orbit_of_site")
if not symmetry.available() or not isinstance(orbit_of, dict):
return [""] * count
report = symmetry.detect(structure)
by_orbit: dict[str, str] = {}
for site, info in zip(structure.sites, report["sites"], strict=True):
by_orbit.setdefault(str(orbit_of.get(site.id)), str(info["wyckoff"]))
orbits = record.get("asymmetric_unit", [])
assert isinstance(orbits, list)
return [by_orbit.get(str(o["site"]), "") for o in orbits]
[docs]
def crystal_parameters(structure: PeriodicStructure) -> tuple[float, ...]:
rows = crystal.lattice_rows(structure)
lengths = np.linalg.norm(rows, axis=1)
def angle(i: int, j: int) -> float:
cosine = float(rows[i] @ rows[j] / (lengths[i] * lengths[j]))
return float(np.degrees(np.arccos(np.clip(cosine, -1, 1))))
return (*(float(x) for x in lengths), angle(1, 2), angle(0, 2), angle(0, 1))
[docs]
def origin_shift(source: Any, target: Any) -> np.ndarray | None:
"""t (multiples of 1/8) with x -> x + t taking ``source``'s operations onto ``target``'s."""
import itertools
wanted = {op.key() for op in target.operations()}
operations = source.operations()
if len(operations) != len(wanted):
return None
for steps in itertools.product(range(8), repeat=3):
t = tuple(Fraction(k, 8) for k in steps)
moved = set()
for op in operations:
r = op.rotation
rt = tuple(sum(r[i][k] * t[k] for k in range(3)) for i in range(3))
moved.add((r, tuple((op.translation[i] + t[i] - rt[i]) % 1 for i in range(3))))
if moved == wanted:
return np.asarray([float(x) for x in t])
return None
# --- (P, p) setting transformations ---------------------------------------------
_TERM = re.compile(r"([+-]?)(\d+(?:/\d+)?|\d*\.\d+)?\*?([abc])(?:/(\d+))?")
[docs]
def parse_transformation(text: str) -> tuple[np.ndarray, np.ndarray]:
"""(P, p) from ``a-b,a+b,c;0,0,1/4`` (ITA Vol. A, Section 1.5)."""
basis, _, origin = text.replace(" ", "").partition(";")
columns = basis.split(",")
if len(columns) != 3 or not all(columns):
raise ValueError(f"transformation {text!r} needs three new basis vectors")
p_matrix = np.zeros((3, 3))
for j, column in enumerate(columns):
position = 0
while position < len(column):
match = _TERM.match(column, position)
if match is None or match.end() == position or not match[3]:
raise ValueError(f"cannot read basis vector {column!r} in {text!r}")
coefficient = float(Fraction(match[2])) if match[2] else 1.0
if match[4]:
coefficient /= int(match[4])
p_matrix["abc".index(match[3]), j] += -coefficient if match[1] == "-" else coefficient
position = match.end()
shift = np.zeros(3)
if origin:
parts = origin.split(",")
if len(parts) != 3:
raise ValueError(f"origin shift in {text!r} needs three components")
shift = np.asarray([float(Fraction(x)) for x in parts])
if np.linalg.det(p_matrix) <= 1e-9:
raise ValueError(f"transformation {text!r} must keep the handedness (det P > 0)")
return p_matrix, shift
[docs]
def transform_setting(
structure: PeriodicStructure, transformation: str, *, identifier: str, tolerance: float = 1e-3
) -> PeriodicStructure:
"""The same crystal in the basis (a,b,c)P with the origin moved to p."""
from vibesketcher.document.symmetry import _rebasis
if sum(structure.pbc) != 3:
raise ValueError("setting transformations apply to 3D crystals")
p_matrix, shift = parse_transformation(transformation)
inverse = np.linalg.inv(p_matrix)
try:
rows, atoms = _rebasis(structure, p_matrix.T, -(inverse @ shift), tolerance)
except ValueError:
raise ValueError(
f"{transformation!r} does not give a cell of this crystal: the new cell would hold "
"a fractional number of atoms, or its vectors are not lattice translations"
) from None
# A cell smaller than the old one needs its vectors to be translations of the crystal.
_check_translations(structure, rows, tolerance)
builder: dict[str, Any] = {
"kind": "setting_transformation",
"source": structure.id,
"transformation": transformation,
"P": p_matrix.tolist(),
"p": shift.tolist(),
"convention": _ITA,
}
return crystal._structure(identifier, rows, atoms, (True, True, True), "unspecified",
builder) # fmt: skip
def _check_translations(structure: PeriodicStructure, rows: np.ndarray, tolerance: float) -> None:
lattice = crystal.lattice_rows(structure)
inverse = np.linalg.inv(lattice)
frac = np.asarray([s.fractional for s in structure.sites], dtype=float)
symbols = [s.symbol for s in structure.sites]
for vector in rows:
step = vector @ inverse
for symbol, f in zip(symbols, frac, strict=True):
moved = f + step
delta = ((frac - moved) + 0.5) % 1 - 0.5
distance = np.linalg.norm(delta @ lattice, axis=1)
if not any(
d < tolerance and s == symbol for d, s in zip(distance, symbols, strict=True)
):
raise ValueError(
"the new basis vectors are not translations of this crystal "
"(a smaller cell needs a centred lattice)"
)
[docs]
def add_to_document(document: Any, structure: PeriodicStructure) -> Any:
from vibesketcher.document.periodic import with_periodic
if any(s.id == structure.id for s in document.periodic_structures):
raise ValueError(f"periodic structure {structure.id!r} already exists")
return with_periodic(document, structure)