vibe-sketcherDocumentation
Source code for vibesketcher.document.periodic
"""Periodic interchange. Cell rows are angstrom; CIF follows IUCr coreCIF.
Gemmi (optional) parses CIF syntax; unsupported symmetry/magnetism refuse rather
than silently becoming molecular or ordered. Original CIF text remains verbatim.
"""
from __future__ import annotations
import base64
import hashlib
import importlib
import json
import math
from typing import Any
from vibesketcher.document.model import (
_PERIODIC_INTEGER_LIMIT,
Document,
PeriodicBond,
PeriodicSite,
PeriodicStructure,
_periodic_transform,
)
[docs]
def with_periodic(document: Document, structure: PeriodicStructure) -> Document:
data = document.to_json_data()
data["periodic_structures"] = [
s.model_dump(mode="json", exclude_none=True)
for s in document.periodic_structures
if s.id != structure.id
] + [structure.model_dump(mode="json", exclude_none=True)]
return Document.model_validate(data)
def _signature(structure: PeriodicStructure) -> str:
data = structure.model_dump(mode="json", exclude={"source_cif", "source_model_sha256"})
return hashlib.sha256(json.dumps(data, sort_keys=True).encode()).hexdigest()
def _gemmi() -> Any:
try:
return importlib.import_module("gemmi")
except ImportError:
raise ValueError(
"CIF needs the periodic extra: pip install 'vibesketcher[periodic]'"
) from None
[docs]
def from_cif(text: str, *, identifier: str = "crystal1") -> PeriodicStructure:
"""One explicit P1 cell; preserve raw syntax/extra tags or refuse unsupported data."""
from vibesketcher.limits import check_size
check_size(text, "CIF")
gemmi = _gemmi()
try:
parsed = gemmi.cif.read_string(text)
if len(parsed) != 1:
raise ValueError("CIF requires exactly one data block")
block = parsed.sole_block()
tags = [
tag
for item in block
for tag in (list(item.loop.tags) if item.loop else [item.pair[0]] if item.pair else [])
]
if any(
str(tag).lower().startswith(("_atom_site_moment", "_space_group_magn")) for tag in tags
):
raise ValueError("magnetic CIF is not admitted")
for tag in (
"_cell_length_a",
"_cell_length_b",
"_cell_length_c",
"_cell_angle_alpha",
"_cell_angle_beta",
"_cell_angle_gamma",
):
if block.find_value(tag) in (None, "", "?", "."):
raise ValueError(f"CIF requires explicit {tag}")
for tag in ("_atom_site_fract_x", "_atom_site_fract_y", "_atom_site_fract_z"):
values = list(block.find_values(tag))
if not values or any(value in ("?", ".") for value in values):
raise ValueError(f"CIF requires explicit {tag}")
if any(value in ("?", ".") for value in block.find_values("_atom_site_occupancy")):
raise ValueError("CIF occupancy must be explicit when stated")
structure = gemmi.make_small_structure_from_block(block)
group = structure.spacegroup_hm.strip()
ops = list(block.find_values("_space_group_symop_operation_xyz")) or list(
block.find_values("_symmetry_equiv_pos_as_xyz")
)
if group not in ("P 1", "P1", "") or any(
gemmi.cif.as_string(op).replace(" ", "").lower() != "x,y,z" for op in ops
):
raise ValueError("CIF symmetry expansion is not admitted; supply explicit P1 sites")
if not group and not ops:
raise ValueError("CIF must explicitly declare P1 symmetry")
cell = structure.cell
# Gemmi fractional->Cartesian provides the IUCr cell metric convention.
vectors = []
for frac in ((1, 0, 0), (0, 1, 0), (0, 0, 1)):
pos = cell.orthogonalize(gemmi.Fractional(*frac))
vectors.append([pos.x, pos.y, pos.z])
sites = [
PeriodicSite(
id=site.label,
symbol=site.element.name,
fractional=(site.fract.x, site.fract.y, site.fract.z),
occupancy=site.occ,
disorder_group=str(site.disorder_group) if site.disorder_group else None,
)
for site in structure.sites
]
result = PeriodicStructure(
id=identifier,
lattice_angstrom=vectors,
sites=sites,
provenance={"format": "cif", "block": block.name},
)
encoded = block.find_value("_vibesketcher_periodic_json")
if encoded:
# Validate the full typed record, and verify its standard CIF projection.
result = PeriodicStructure.model_validate_json(
base64.b64decode(gemmi.cif.as_string(encoded), validate=True)
)
if len(result.sites) != len(sites) or any(
a.id != b.id
or a.symbol != b.symbol
or a.disorder_group != b.disorder_group
or abs(a.occupancy - b.occupancy) > 1e-10
or any(abs(x - y) > 1e-9 for x, y in zip(a.fractional, b.fractional, strict=True))
for a, b in zip(result.sites, sites, strict=True)
):
raise ValueError("CIF standard sites conflict with periodic metadata")
if any(
abs(x - y) > 1e-8
for row, other in zip(result.lattice_angstrom, vectors, strict=True)
for x, y in zip(row, other, strict=True)
):
# CIF lengths/angles normalize orientation, so compare Gram matrices.
for i in range(3):
for j in range(3):
a = sum(
result.lattice_angstrom[i][k] * result.lattice_angstrom[j][k]
for k in range(3)
)
b = sum(vectors[i][k] * vectors[j][k] for k in range(3))
if abs(a - b) > 1e-7 * max(1, abs(a)):
raise ValueError("CIF metric conflicts with periodic metadata")
result.source_cif = text
result.source_model_sha256 = _signature(result)
return result
except (RuntimeError, TypeError, KeyError, AttributeError) as error:
raise ValueError(f"unsupported or malformed CIF: {error}") from None
[docs]
def to_cif(structure: PeriodicStructure) -> str:
"""Lossless raw readback; normalized P1 plus typed metadata for new structures."""
# Models are mutable: revalidate even an SDK model_copy or edited instance.
structure = PeriodicStructure.model_validate(structure.model_dump())
if structure.source_cif is not None:
if _signature(structure) != structure.source_model_sha256:
raise ValueError(
"CIF source model changed; clear source receipt explicitly before export"
)
reparsed = from_cif(structure.source_cif, identifier=structure.id)
if _signature(reparsed) != _signature(structure):
raise ValueError("CIF source receipt conflicts with the periodic model")
return structure.source_cif
gemmi = _gemmi()
rows = structure.lattice_angstrom
lengths = [math.sqrt(sum(x * x for x in row)) for row in rows]
angles = [
math.degrees(
math.acos(
max(
-1,
min(
1,
sum(x * y for x, y in zip(rows[i], rows[j], strict=True))
/ (lengths[i] * lengths[j]),
),
)
)
)
for i, j in ((1, 2), (0, 2), (0, 1))
]
doc = gemmi.cif.Document()
block = doc.add_new_block(structure.id)
for tag, value in zip(
("a", "b", "c", "alpha", "beta", "gamma"), [*lengths, *angles], strict=True
):
block.set_pair(
"_cell_length_" + tag if tag in ("a", "b", "c") else "_cell_angle_" + tag,
format(value, ".17g"),
)
block.set_pair("_space_group_name_H-M_alt", "'P 1'")
loop = block.init_loop(
"_atom_site_",
["label", "type_symbol", "fract_x", "fract_y", "fract_z", "occupancy", "disorder_group"],
)
for site in structure.sites:
if site.disorder_group and not site.disorder_group.isdecimal():
raise ValueError("normalized CIF requires numeric disorder groups")
loop.add_row(
[
site.id,
site.symbol,
*[format(x, ".17g") for x in site.fractional],
str(site.occupancy),
site.disorder_group or ".",
]
)
encoded = base64.b64encode(
structure.model_dump_json(exclude={"source_cif", "source_model_sha256"}).encode()
).decode()
block.set_pair("_vibesketcher_periodic_json", encoded)
text = str(doc.as_string())
if _signature(from_cif(text)) != _signature(structure):
raise ValueError("normalized CIF failed periodic model readback")
return text
[docs]
def from_qvf_structure(
payload: dict[str, Any], section: dict[str, Any], provenance: dict[str, Any]
) -> PeriodicStructure:
"""Read native QVF cell rows/positions; never infer a molecule from PBC."""
import numpy as np
if len(payload.get("atoms", [])) > 25000:
raise ValueError("periodic atom limit is 25000")
rows = payload.get("lattice_vectors")
if rows is None:
raise ValueError("PBC structure is missing lattice_vectors")
# QVF structure positions and lattice rows are angstrom, spec section 5.1.
sites = []
for index, atom in enumerate(payload["atoms"]):
try:
fractional = np.linalg.solve(np.asarray(rows).T, np.asarray(atom["position"]))
except np.linalg.LinAlgError:
raise ValueError("periodic cell must be nonsingular") from None
sites.append(
PeriodicSite(
id=atom.get("atom_name") or f"s{index + 1}",
symbol=atom["symbol"],
fractional=fractional.tolist(),
occupancy=atom.get("occupancy", 1),
)
)
bonds = []
for pair in payload.get("source_bonds", {}).get("pairs", []):
try:
a, b = sites[pair["i"]], sites[pair["j"]]
if pair["i"] < 0 or pair["j"] < 0:
raise ValueError("negative bond index")
bonds.append(
PeriodicBond(
a=a.id, b=b.id, image=pair.get("image", (0, 0, 0)), order=pair.get("order", 1)
)
)
except (IndexError, KeyError, TypeError) as error:
raise ValueError(f"invalid periodic QVF bond: {error}") from None
return PeriodicStructure(
id=section["id"],
lattice_angstrom=rows,
pbc=payload["pbc"],
sites=sites,
bonds=bonds,
charge=section.get("charge"),
multiplicity=section.get("multiplicity"),
provenance={
"format": "qvf",
"source": provenance,
"structure": payload,
"section": section,
},
)
[docs]
def supercell(
structure: PeriodicStructure, matrix: list[list[int]], *, identifier: str
) -> PeriodicStructure:
"""Integer row-cell transform with explicit parent sites and lattice images.
Fractional row coordinates transform by the inverse matrix (IUCr coreCIF
cell conventions). Coset representatives enumerate Z^3 / rows(matrix).
This changes the physical model; display repetitions are a separate concern.
"""
import itertools
import numpy as np
structure = PeriodicStructure.model_validate(structure.model_dump())
count = _periodic_transform(matrix)
for bond in structure.bonds:
if any(abs(x) > _PERIODIC_INTEGER_LIMIT for x in bond.image):
raise ValueError("supercell bond images exceed the integer arithmetic limit")
transform = np.asarray(matrix, dtype=float)
if count * len(structure.bonds) > 100000:
raise ValueError("supercell exceeds 100000 bonds")
if count * len(structure.sites) > 25000:
raise ValueError("supercell exceeds 25000 sites")
if any(
not structure.pbc[i]
and matrix[i] != [int(i == j) for j in range(3)]
or any(matrix[i][j] and structure.pbc[i] != structure.pbc[j] for j in range(3))
for i in range(3)
):
raise ValueError("supercell must preserve nonperiodic axes")
try:
inverse = np.linalg.inv(transform)
except np.linalg.LinAlgError:
raise ValueError("supercell transform is numerically singular") from None
corners = np.asarray(list(itertools.product((0, 1), repeat=3))) @ transform
lower = np.floor(corners.min(axis=0)).astype(int)
upper = np.ceil(corners.max(axis=0)).astype(int)
volume = math.prod(int(upper[i] - lower[i] + 1) for i in range(3))
if volume > 1_000_000:
raise ValueError("supercell enumeration exceeds resource limit")
images = []
for image in itertools.product(*(range(lower[i], upper[i] + 1) for i in range(3))):
fractional = np.asarray(image) @ inverse
if all(-1e-9 <= x < 1 - 1e-9 for x in fractional):
images.append(image)
if len(images) != count:
raise ValueError("supercell representatives could not be qualified")
def reduce(site: PeriodicSite, image: Any) -> tuple[Any, Any, Any]:
with np.errstate(over="ignore", invalid="ignore"):
frac = (np.asarray(site.fractional) + image) @ inverse
if not np.isfinite(frac).all() or (np.abs(frac) > _PERIODIC_INTEGER_LIMIT).any():
raise ValueError("supercell image bookkeeping exceeds the integer arithmetic limit")
shift = np.floor(frac + 1e-10).astype(int)
canonical = np.asarray(image) - shift @ np.asarray(matrix)
return frac - shift, canonical, shift
sites: list[PeriodicSite] = []
mapped: dict[Any, str] = {}
origins = []
for site in structure.sites:
for image in images:
fractional, canonical, _ = reduce(site, image)
name = f"s{len(sites) + 1}"
sites.append(
site.model_copy(
update={
"id": name,
"fractional": tuple(float(x) for x in fractional),
"parent_site": site.id,
"image": tuple(int(x) for x in canonical),
}
)
)
mapped[(site.id, tuple(canonical))] = name
origins.append((site, canonical, name))
bonds = []
by_id = {site.id: site for site in structure.sites}
outgoing: dict[str, list[PeriodicBond]] = {}
for bond in structure.bonds:
outgoing.setdefault(bond.a, []).append(bond)
for source, image, name in origins:
for bond in outgoing.get(source.id, []):
_, canonical, shift = reduce(by_id[bond.b], image + np.asarray(bond.image))
bonds.append(
PeriodicBond(
a=name,
b=mapped[(bond.b, tuple(canonical))],
image=tuple(int(x) for x in shift),
order=bond.order,
)
)
return PeriodicStructure(
id=identifier,
lattice_angstrom=(transform @ np.asarray(structure.lattice_angstrom)).tolist(),
pbc=structure.pbc,
sites=sites,
bonds=bonds,
parent=structure.id,
supercell_matrix=matrix,
provenance={"operation": "integer supercell"},
charge=None if structure.charge is None else structure.charge * count,
# Total spin is not additive. Do not invent a supercell multiplicity.
multiplicity=None,
)