vibe-sketcherDocumentation
Source code for vibesketcher.document.cif_import
"""Experimental CIF import with space-group expansion (#228).
A small pure-Python CIF 1.1 reader: it needs no compiled dependency, so it
works wherever the sketcher does (Gemmi does not install on Python 3.15 yet,
#210). It reads the cell with its standard uncertainties, the symmetry
operations (or, without them, the Hall or Hermann-Mauguin symbol), the atom
sites with occupancies and disorder groups, and expands the asymmetric unit
to the full conventional cell. Partial occupancy is refused unless the caller
chooses how to resolve it; every choice and merge is recorded.
The lossless P1 interchange in :mod:`vibesketcher.document.periodic` is
unchanged; this module is for files from experiments and databases.
Syntax: S. R. Hall, F. H. Allen, I. D. Brown, Acta Cryst. A47, 655 (1991).
"""
from __future__ import annotations
import hashlib
import math
import re
from collections import Counter
from dataclasses import dataclass
from fractions import Fraction
from typing import Any, Literal
import numpy as np
from vibesketcher.document import crystal
from vibesketcher.document.model import PeriodicStructure
from vibesketcher.document.spacegroups import Operation, Setting, compose, find, settings
from vibesketcher.elements import SYMBOLS
Disorder = Literal["refuse", "majority", "group"]
_REFERENCE = "S. R. Hall, F. H. Allen, I. D. Brown, Acta Cryst. A47, 655 (1991)"
# --- tokenizer ---------------------------------------------------------------------
[docs]
@dataclass
class Block:
name: str
items: dict[str, str]
loops: list[dict[str, list[str]]]
[docs]
def value(self, *tags: str) -> str | None:
for tag in tags:
if tag in self.items:
return self.items[tag]
for loop in self.loops:
if tag in loop and len(loop[tag]) == 1:
return loop[tag][0]
return None
[docs]
def column(self, *tags: str) -> tuple[dict[str, list[str]], list[str]] | None:
for tag in tags:
for loop in self.loops:
if tag in loop:
return loop, loop[tag]
if tag in self.items:
return {tag: [self.items[tag]]}, [self.items[tag]]
return None
def _tokens(text: str) -> list[tuple[str, bool]]:
"""(token, quoted) pairs; quoted values are never tags or keywords."""
out: list[tuple[str, bool]] = []
lines = text.replace("\r\n", "\n").replace("\r", "\n").split("\n")
i = 0
while i < len(lines):
line = lines[i]
if line.startswith(";"):
field = [line[1:]]
i += 1
while i < len(lines) and not lines[i].startswith(";"):
field.append(lines[i])
i += 1
if i == len(lines):
raise ValueError("CIF text field (;) is not closed")
out.append(("\n".join(field).strip("\n"), True))
i += 1
continue
pos = 0
while pos < len(line):
ch = line[pos]
if ch.isspace():
pos += 1
elif ch == "#":
break
elif ch in "'\"":
# A quote closes only when followed by whitespace or the line end.
end = pos + 1
while end < len(line) and not (
line[end] == ch and (end + 1 == len(line) or line[end + 1].isspace())
):
end += 1
if end >= len(line):
raise ValueError(f"CIF line {i + 1}: unclosed quote")
out.append((line[pos + 1 : end], True))
pos = end + 1
else:
end = pos
while end < len(line) and not line[end].isspace():
end += 1
out.append((line[pos:end], False))
pos = end
i += 1
return out
[docs]
def parse(text: str) -> list[Block]:
"""The data blocks of a CIF 1.1 text (tags lower-cased)."""
blocks: list[Block] = []
tokens = _tokens(text)
i = 0
current: Block | None = None
while i < len(tokens):
token, quoted = tokens[i]
lower = token.lower()
if not quoted and lower.startswith("data_"):
current = Block(token[5:], {}, [])
blocks.append(current)
i += 1
continue
if current is None:
if not quoted and lower.startswith(("global_", "save_")):
i += 1
continue
raise ValueError("CIF content before the first data_ block")
if not quoted and lower == "loop_":
i += 1
tags = []
while i < len(tokens) and not tokens[i][1] and tokens[i][0].startswith("_"):
tags.append(tokens[i][0].lower())
i += 1
if not tags:
raise ValueError("CIF loop_ without tags")
values = []
while i < len(tokens):
value, q = tokens[i]
if not q and (value.startswith("_") or value.lower() in ("loop_",)
or value.lower().startswith(("data_", "save_"))): # fmt: skip
break
values.append(value)
i += 1
if len(values) % len(tags):
raise ValueError(f"CIF loop with {tags[0]} has a ragged value count")
loop = {t: values[k :: len(tags)] for k, t in enumerate(tags)}
current.loops.append(loop)
continue
if not quoted and token.startswith("_"):
if i + 1 >= len(tokens):
raise ValueError(f"CIF tag {token} has no value")
current.items[lower] = tokens[i + 1][0]
i += 2
continue
raise ValueError(f"unexpected CIF value {token!r}")
if not blocks:
raise ValueError("the text has no CIF data block")
return blocks
# --- values ------------------------------------------------------------------------
_NUMBER = re.compile(r"^([-+]?(?:\d+\.?\d*|\.\d+)(?:[eE][-+]?\d+)?)(?:\((\d+)\))?$")
[docs]
def number(text: str | None, tag: str) -> tuple[float, float | None]:
"""A CIF number and its standard uncertainty, e.g. 4.7602(3) -> (4.7602, 0.0003)."""
if text is None or text in ("?", "."):
raise ValueError(f"CIF requires an explicit {tag}")
match = _NUMBER.match(text.strip())
if match is None:
raise ValueError(f"CIF {tag} {text!r} is not a number")
value = float(match[1])
if match[2] is None:
return value, None
mantissa = match[1].lower().split("e")[0]
decimals = len(mantissa.split(".")[1]) if "." in mantissa else 0
exponent = int(match[1].lower().split("e")[1]) if "e" in match[1].lower() else 0
return value, int(match[2]) * 10.0 ** (exponent - decimals)
def _element(type_symbol: str | None, label: str) -> str:
for text in (type_symbol, label):
if not text or text in ("?", "."):
continue
match = re.match(r"([A-Za-z]{1,2})", text)
if match is None:
continue
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"CIF site {label!r} has no recognisable element")
# --- symmetry ----------------------------------------------------------------------
def _exact(rotation: np.ndarray, translation: np.ndarray) -> Operation:
r = tuple(tuple(int(round(x)) for x in row) for row in rotation)
t = tuple(Fraction(float(x)).limit_denominator(12) % 1 for x in translation)
residual = (np.asarray([float(x) for x in t]) - translation) % 1.0
if not np.allclose(rotation, r) or np.any(np.minimum(residual, 1 - residual) > 1e-4):
raise ValueError("CIF symmetry operation is not a crystallographic operation")
return Operation(r, t) # type: ignore[arg-type]
def _closed(operations: list[Operation]) -> bool:
keys = {op.key() for op in operations}
return all(compose(a, b).key() in keys for a in operations for b in operations)
def _setting_of(operations: list[Operation]) -> Setting | None:
keys = {op.key() for op in operations}
for setting in settings():
theirs = setting.operations()
if len(theirs) == len(keys) and {op.key() for op in theirs} == keys:
return setting
return None
def _symmetry(block: Block, setting_hint: str | None) -> tuple[list[Operation], dict[str, Any]]:
found = block.column("_space_group_symop_operation_xyz", "_symmetry_equiv_pos_as_xyz")
hall = block.value("_space_group_name_hall", "_symmetry_space_group_name_hall")
hm = block.value("_space_group_name_h-m_alt", "_symmetry_space_group_name_h-m")
it_number = block.value("_space_group_it_number", "_symmetry_int_tables_number")
record: dict[str, Any] = {"hall_in_file": hall, "hm_in_file": hm, "number_in_file": it_number}
if found is not None:
operations = []
for text in found[1]:
rotation, translation = crystal.parse_operation(text)
operations.append(_exact(rotation, translation))
unique = list({op.key(): op for op in operations}.values())
if len(unique) != len(operations):
record["duplicate_operations"] = len(operations) - len(unique)
if not _closed(unique):
raise ValueError("the CIF symmetry operations do not form a group")
setting = _setting_of(unique)
record["operations_from"] = "symop loop"
if setting is not None:
claimed = None
for named in (setting_hint, hall, hm, it_number):
if named and named not in ("?", "."):
try:
claimed = find(named)
except ValueError:
continue
break
if claimed is not None and claimed.number != setting.number:
raise ValueError(
f"the CIF operations are those of {setting.symbol} (No. {setting.number}), "
f"but the file names {claimed.symbol} (No. {claimed.number})"
)
record["setting"] = {"number": setting.number, "hm": setting.hm,
"extension": setting.extension, "hall": setting.hall} # fmt: skip
else:
record["setting"] = None
record["warning"] = "the operations are not one of the 564 tabulated settings"
return unique, record
symbol = setting_hint or hall or hm or it_number
if not symbol or symbol in ("?", "."):
raise ValueError("the CIF gives neither symmetry operations nor a space-group symbol")
setting = find(symbol)
if setting_hint is None and hall is None:
# Without operations or a Hall symbol, a symbol with two settings is ambiguous.
alternatives = [s for s in settings() if s.number == setting.number
and s.hm == setting.hm and s.extension != setting.extension] # fmt: skip
if alternatives and not re.search(r":\s*\S+\s*$", symbol):
names = ", ".join(f"{s.hm}:{s.extension}" for s in [setting, *alternatives])
raise ValueError(
f"the CIF names {symbol!r} without operations; it has several settings ({names}). "
"Pass the setting explicitly"
)
record["operations_from"] = "symbol"
record["setting"] = {"number": setting.number, "hm": setting.hm,
"extension": setting.extension, "hall": setting.hall} # fmt: skip
return setting.operations(), record
# --- import ------------------------------------------------------------------------
@dataclass
class _Site:
label: str
element: str
fractional: np.ndarray
esd: list[float | None]
occupancy: float
group: str | None
assembly: str | None
def _sites(block: Block) -> list[_Site]:
found = block.column("_atom_site_fract_x")
if found is None:
if block.column("_atom_site_cartn_x") is not None:
raise ValueError("CIF Cartesian atom sites are not admitted; give fractional ones")
raise ValueError("the CIF has no _atom_site_fract_x loop")
loop = found[0]
def col(tag: str) -> list[str] | None:
return loop.get(tag)
count = len(found[1])
labels = col("_atom_site_label") or [f"site{i + 1}" for i in range(count)]
types = col("_atom_site_type_symbol")
occupancies = col("_atom_site_occupancy")
groups = col("_atom_site_disorder_group")
assemblies = col("_atom_site_disorder_assembly")
flags = col("_atom_site_calc_flag")
out = []
for i in range(count):
if flags and flags[i].lower() == "dum":
continue
coords, esd = [], []
for axis in "xyz":
values = col(f"_atom_site_fract_{axis}")
assert values is not None
value, uncertainty = number(values[i], f"_atom_site_fract_{axis} of {labels[i]}")
coords.append(value)
esd.append(uncertainty)
occupancy = 1.0
if occupancies and occupancies[i] not in ("?", "."):
occupancy = number(occupancies[i], f"occupancy of {labels[i]}")[0]
if not 0 < occupancy <= 1.0 + 1e-6:
raise ValueError(f"CIF site {labels[i]} has occupancy {occupancy}")
group = groups[i] if groups and groups[i] not in ("?", ".") else None
assembly = assemblies[i] if assemblies and assemblies[i] not in ("?", ".") else None
out.append(_Site(labels[i], _element(types[i] if types else None, labels[i]),
np.asarray(coords, dtype=float), esd, min(occupancy, 1.0), group,
assembly)) # fmt: skip
if not out:
raise ValueError("the CIF has no atom sites")
return out
def _resolve_disorder(
sites: list[_Site], policy: Disorder, keep_group: str | None, tolerance: float
) -> tuple[list[_Site], list[dict[str, Any]]]:
partial = [s for s in sites if s.occupancy < 1 - tolerance]
grouped = [s for s in sites if s.group is not None]
if not partial and not grouped:
return sites, []
if policy == "refuse":
names = ", ".join(f"{s.label} ({s.occupancy:g})" for s in (partial or grouped)[:8])
raise ValueError(
f"the CIF has partially occupied or disordered sites: {names}. Choose how to "
"resolve them: disorder='majority' or disorder='group' with a disorder group"
)
decisions: list[dict[str, Any]] = []
kept = list(sites)
if policy == "group":
if keep_group is None:
raise ValueError("disorder='group' needs the disorder group to keep")
available = sorted({s.group for s in grouped if s.group})
if keep_group not in available:
raise ValueError(
f"disorder group {keep_group!r} not in the CIF ({', '.join(available)})"
)
dropped = [s for s in kept if s.group is not None and s.group != keep_group]
kept = [s for s in kept if s.group is None or s.group == keep_group]
decisions.append({"policy": "group", "kept_group": keep_group,
"dropped": [s.label for s in dropped]}) # fmt: skip
else:
# Per assembly, keep the disorder group with the largest mean occupancy.
by_assembly: dict[str | None, dict[str, list[_Site]]] = {}
for s in grouped:
by_assembly.setdefault(s.assembly, {}).setdefault(s.group or "", []).append(s)
for assembly, groups in by_assembly.items():
best = max(groups, key=lambda g: (np.mean([s.occupancy for s in groups[g]]), g))
dropped = [s for g, members in groups.items() if g != best for s in members]
kept = [s for s in kept if s not in dropped]
decisions.append({"policy": "majority", "assembly": assembly, "kept_group": best,
"dropped": [s.label for s in dropped]}) # fmt: skip
# Mixed or vacancy sites that remain: keep the most occupied species per position.
remaining_partial = [s for s in kept if s.occupancy < 1 - tolerance and s.group is None]
if remaining_partial and policy == "group":
names = ", ".join(s.label for s in remaining_partial[:8])
raise ValueError(f"sites outside disorder groups are partially occupied: {names}; "
"use disorder='majority'") # fmt: skip
return kept, decisions
[docs]
def import_cif(
text: str,
*,
identifier: str,
block: str | None = None,
disorder: Disorder = "refuse",
disorder_group: str | None = None,
setting: str | None = None,
merge_tolerance: float = 0.01,
occupancy_tolerance: float = 0.01,
) -> PeriodicStructure:
"""The full conventional cell of one CIF data block (#228).
``disorder`` decides partially occupied and disordered sites:
``"refuse"`` (default), ``"group"`` keeps ``disorder_group`` and drops the
other groups, ``"majority"`` keeps the most occupied group of each disorder
assembly and, on a mixed or partially vacant site, the most occupied
species. Kept sites are written fully occupied; every decision is recorded.
``setting`` names the space-group setting when the file has no operations
and its symbol is ambiguous (e.g. ``"Fd-3m:1"``). Images closer than
``merge_tolerance`` angstrom are merged.
"""
from vibesketcher.limits import check_size
check_size(text, "CIF")
blocks = parse(text)
if block is not None:
chosen = [b for b in blocks if b.name.lower() == block.lower()]
if not chosen:
raise ValueError(
f"CIF data block {block!r} not found ({', '.join(b.name for b in blocks)})"
)
data = chosen[0]
else:
with_sites = [b for b in blocks if b.column("_atom_site_fract_x") is not None]
if len(with_sites) != 1:
names = ", ".join(b.name for b in (with_sites or blocks))
raise ValueError(
f"the CIF has {len(with_sites)} structure blocks ({names}); choose one"
)
data = with_sites[0]
if any(tag.startswith(("_atom_site_moment", "_space_group_magn"))
for tag in [*data.items, *(t for loop in data.loops for t in loop)]): # fmt: skip
raise ValueError("magnetic CIF is not admitted")
cell: dict[str, float] = {}
cell_esd: dict[str, float | None] = {}
for key in ("length_a", "length_b", "length_c", "angle_alpha", "angle_beta", "angle_gamma"):
cell[key], cell_esd[key] = number(data.value(f"_cell_{key}"), f"_cell_{key}")
if min(cell["length_a"], cell["length_b"], cell["length_c"]) <= 0:
raise ValueError("CIF cell lengths must be positive")
lattice = crystal.cell_from_parameters(cell["length_a"], cell["length_b"], cell["length_c"],
cell["angle_alpha"], cell["angle_beta"],
cell["angle_gamma"]) # fmt: skip
operations, symmetry = _symmetry(data, setting)
sites, decisions = _resolve_disorder(_sites(data), disorder, disorder_group,
occupancy_tolerance) # fmt: skip
# Expand every asymmetric site; merge images within the tolerance (minimum image).
expanded: list[tuple[_Site, np.ndarray]] = []
orbits = []
for site in sites:
images: list[np.ndarray] = []
for op in operations:
point = np.asarray(op.apply(tuple(site.fractional)), dtype=float) % 1.0
if not any(_distance(point, other, lattice) < merge_tolerance for other in images):
images.append(point)
orbits.append({"label": site.label, "element": site.element,
"fractional": site.fractional.tolist(), "esd": site.esd,
"occupancy": site.occupancy, "multiplicity": len(images),
"disorder_group": site.group}) # fmt: skip
expanded.extend((site, p) for p in images)
merges: list[dict[str, Any]] = []
final: list[tuple[_Site, np.ndarray]] = []
clash_radius = max(merge_tolerance, 0.5) # closer atoms are disorder or an error
for site, point in expanded:
clash = next((k for k, (_other, q) in enumerate(final)
if _distance(point, q, lattice) < clash_radius), None) # fmt: skip
if clash is None:
final.append((site, point))
continue
other = final[clash][0]
distance = _distance(point, final[clash][1], lattice)
partial = (
site.occupancy < 1 - occupancy_tolerance or other.occupancy < 1 - occupancy_tolerance
)
if disorder == "majority" and partial:
winner = (
other if (other.occupancy, other.label) >= (site.occupancy, site.label) else site
)
loser = site if winner is other else other
final[clash] = (winner, final[clash][1] if winner is other else point)
merges.append({"kept": winner.label, "dropped": loser.label,
"distance_angstrom": round(distance, 4)}) # fmt: skip
continue
if other is site:
raise ValueError(
f"CIF site {site.label} lies {distance:.3f} A from a special position; "
"correct it or raise merge_tolerance"
)
raise ValueError(
f"CIF sites {other.label} and {site.label} generate atoms {distance:.3f} A apart"
+ ("; they look like disorder: choose a disorder policy" if partial else "")
)
if len(final) > 25000:
raise ValueError("the expanded crystal exceeds 25000 sites")
if disorder != "refuse":
merged = {m["kept"] for m in merges}
filled = sorted({s.label for s, _p in final
if s.occupancy < 1 - occupancy_tolerance and s.label not in merged
and s.group is None}) # fmt: skip
if filled:
decisions.append({"policy": disorder, "written_fully_occupied": filled})
counts = Counter(s.element for s, _p in final)
warnings = []
if symmetry.get("warning"):
warnings.append(symmetry["warning"])
formula_sum = data.value("_chemical_formula_sum")
z_text = data.value("_cell_formula_units_z")
if formula_sum and z_text and not decisions and not merges:
expected: dict[str, float] = {}
for element, amount in re.findall(r"([A-Z][a-z]?)\s*([\d.]*)", formula_sum):
expected[element] = expected.get(element, 0.0) + float(amount or 1)
try:
z = number(z_text, "_cell_formula_units_Z")[0]
except ValueError:
z = math.nan
if any(abs(counts[e] - n * z) > 0.05 for e, n in expected.items()) or set(counts) - set(
expected
):
warnings.append(f"the expanded cell {dict(counts)} does not match "
f"_chemical_formula_sum {formula_sum!r} with Z = {z_text}") # fmt: skip
source = {
key: data.value(tag)
for key, tag in (
("doi", "_journal_paper_doi"),
("journal", "_journal_name_full"),
("year", "_journal_year"),
("volume", "_journal_volume"),
("page", "_journal_page_first"),
("cod", "_cod_database_code"),
("icsd", "_database_code_icsd"),
("temperature_k", "_cell_measurement_temperature"),
("diffraction_temperature_k", "_diffrn_ambient_temperature"),
("formula_sum", "_chemical_formula_sum"),
) # fmt: skip
if data.value(tag) not in (None, "?", ".")
}
builder = {
"kind": "cif_import",
"block": data.name,
"sha256": hashlib.sha256(text.encode()).hexdigest(),
"cell": cell,
"cell_esd": cell_esd,
"space_group": symmetry,
"operations": len(operations),
"asymmetric_unit": orbits,
"disorder": decisions,
"merges": merges,
"merge_tolerance": merge_tolerance,
"warnings": warnings,
"source": source,
"reference": _REFERENCE,
}
labels = [s.label for s, _p in final]
ids = crystal._site_ids([s.element for s, _p in final])
structure = crystal._structure(
identifier, lattice, [(s.element, p) for s, p in final], (True, True, True),
"conventional", builder, site_ids=ids,
) # fmt: skip
provenance = dict(structure.provenance)
provenance["format"] = "cif"
provenance["cif_label_of_site"] = dict(zip(ids, labels, strict=True))
return structure.model_copy(update={"provenance": provenance})
def _distance(a: np.ndarray, b: np.ndarray, lattice: np.ndarray) -> float:
delta = a - b
delta -= np.round(delta)
return float(np.linalg.norm(delta @ lattice))
[docs]
def add_to_document(document: Any, text: str, *, identifier: str, **options: Any) -> Any:
"""A copy of ``document`` with the imported crystal added; the ID must be new."""
from vibesketcher.document.periodic import with_periodic
if any(s.id == identifier for s in document.periodic_structures):
raise ValueError(f"periodic structure {identifier!r} already exists")
return with_periodic(document, import_cif(text, identifier=identifier, **options))