"""DFTB+ ``.gen`` geometry -> CIF text: the export writer (M10).
Pure stdlib and pure computation -- no science imports, no AiiDA, no I/O, no
Rich. The ``cif`` command feeds this the ``geom.out.gen`` text that DFTB+
relaxed to (fetched from the user's AiiDA repository, or cached in
``calc_result.json``) and writes the returned document to disk; keeping the
conversion here means every cell/coordinate rule is unit-testable in the
science-free dev venv, and no CIF is ever produced by a dependency we do not
own.
The written document is deliberately plain: space group ``P 1`` (a relaxed
DFTB+ cell carries no symmetry claim -- the atoms have moved), one
``_atom_site_*`` loop in fractional coordinates, and whatever provenance the
caller passes as ``_mcli_*`` items (material id, WorkChain/CalcJob pk,
energy, ...). It round-trips through the extract runner's own CIF reader,
and VESTA/ASE read it directly.
Two deliberate non-transformations, so an exported CIF is exactly what the
calculation produced:
* coordinates are **not** wrapped into ``[0, 1)`` -- an atom that relaxed to
``0.999998`` stays there rather than jumping across the cell;
* the ``.gen`` origin row is ignored (as ASE does) -- lattice vectors alone
define the cell.
``extract_runner.py`` holds the *inverse* helpers (CIF -> cell vectors); it
is a self-contained runner staged to compute backends and may never import
this module, so the two directions live apart by design.
"""
from __future__ import annotations
import math
import re
from collections.abc import Mapping, Sequence
from dataclasses import dataclass
from typing import Any
#: A cartesian or fractional triple.
Vector = tuple[float, float, float]
#: Three lattice vectors (rows), in angstrom.
Lattice = tuple[Vector, Vector, Vector]
#: Provenance keys the caller may attach; they become ``_mcli_<key>`` items.
_EXTRA_KEY_RE = re.compile(r"[A-Za-z0-9_]+")
#: Below this, a float is float noise and serializes as an honest zero
#: (``cos(90 deg)`` is ~6e-17, and ``-0.00000000`` in a CIF helps nobody).
_ZERO_SNAP = 5e-9
[docs]
@dataclass(frozen=True)
class Geometry:
"""One periodic structure, ready to serialize.
Attributes
----------
species : tuple of str
One element symbol per site, already expanded from the ``.gen`` species
table.
fractional : tuple
One ``(x, y, z)`` fractional triple per site, in the same order as
``species``.
lattice : tuple
The three cell vectors (rows), in angstrom.
"""
species: tuple[str, ...]
fractional: tuple[Vector, ...]
lattice: Lattice
# ---------------------------------------------------------------------------
# .gen parsing
# ---------------------------------------------------------------------------
def _floats(tokens: Sequence[str], row: str) -> Vector:
"""Read three floats from a ``.gen`` row.
Parameters
----------
tokens
The row's whitespace-split tokens; the first three are used.
row
The original row text, quoted back in the error message.
Returns
-------
tuple
The ``(x, y, z)`` triple.
Raises
------
ValueError
If any of the three tokens is not a number. The message names the
offending row -- a geometry that cannot be read must never be guessed
at.
"""
try:
x, y, z = (float(token) for token in tokens[:3])
except ValueError:
raise ValueError(f"gen row has non-numeric coordinates: {row!r}") from None
return (x, y, z)
[docs]
def parse_gen(text: str) -> Geometry:
"""Parse DFTB+ ``.gen`` text into a :class:`Geometry`.
Parameters
----------
text
The file's contents. Blank lines and ``#`` comments are ignored.
Returns
-------
Geometry
The structure, always in fractional coordinates.
Raises
------
ValueError
For a cluster geometry (``C``), which has no cell to write; for an
unknown geometry type; and for any malformed header, species row, atom
row, or lattice block. Every rejection names the offending row -- a
geometry that cannot be read must never be guessed at.
Notes
-----
Both periodic forms are supported: ``F`` (fractional -- what our HSD
template and DFTB+'s ``geom.out.gen`` produce) and ``S`` (cartesian,
converted here against the lattice).
"""
lines = [
stripped
for stripped in (line.strip() for line in text.splitlines())
if stripped and not stripped.startswith("#")
]
if len(lines) < 3:
raise ValueError("gen text is too short to hold a header, a species row and atoms")
header = lines[0].split()
if len(header) < 2:
raise ValueError(f"gen header must read '<natoms> <C|S|F>', got {lines[0]!r}")
try:
natoms = int(header[0])
except ValueError:
raise ValueError(f"gen header has a non-integer atom count: {lines[0]!r}") from None
if natoms < 1:
raise ValueError(f"gen header declares {natoms} atoms")
kind = header[1].upper()
if kind == "C":
raise ValueError(
"gen geometry is a cluster ('C') and carries no lattice; a CIF needs a periodic cell"
)
if kind not in {"S", "F"}:
raise ValueError(f"unknown gen geometry type {header[1]!r} (expected C, S or F)")
species = lines[1].split()
if not all(symbol.isalpha() for symbol in species) or not species:
raise ValueError(f"gen species row does not look like element symbols: {lines[1]!r}")
atom_rows = lines[2 : 2 + natoms]
if len(atom_rows) < natoms:
raise ValueError(
f"gen header declares {natoms} atoms but only {len(atom_rows)} atom row(s) follow"
)
symbols: list[str] = []
coordinates: list[Vector] = []
for row in atom_rows:
tokens = row.split()
if len(tokens) < 5:
raise ValueError(f"gen atom row must read '<index> <species> <x> <y> <z>': {row!r}")
try:
species_index = int(tokens[1])
except ValueError:
raise ValueError(f"gen atom row has a non-integer species index: {row!r}") from None
if not 1 <= species_index <= len(species):
raise ValueError(
f"gen atom row names species {species_index}, outside 1..{len(species)}: {row!r}"
)
symbols.append(species[species_index - 1])
coordinates.append(_floats(tokens[2:5], row))
tail = lines[2 + natoms :]
if len(tail) < 4:
raise ValueError(
"periodic gen geometry needs an origin row and 3 lattice vectors after the atoms, "
f"found {len(tail)} row(s)"
)
vectors = tuple(_floats(row.split(), row) for row in tail[1:4])
lattice: Lattice = (vectors[0], vectors[1], vectors[2])
if kind == "S":
coordinates = [cart_to_frac(position, lattice) for position in coordinates]
return Geometry(tuple(symbols), tuple(coordinates), lattice)
# ---------------------------------------------------------------------------
# Cell geometry (the inverse of extract_runner's cell_to_lattice / frac_to_cart)
# ---------------------------------------------------------------------------
[docs]
def cell_volume(lattice: Lattice) -> float:
"""Compute the cell volume.
Parameters
----------
lattice
The three cell vectors, in angstrom.
Returns
-------
float
Volume in cubic angstrom -- the absolute determinant, so a
left-handed cell reports a positive volume.
Examples
--------
>>> cell_volume(((4.0, 0.0, 0.0), (0.0, 4.0, 0.0), (0.0, 0.0, 4.0)))
64.0
"""
(ax, ay, az), (bx, by, bz), (cx, cy, cz) = lattice
return abs(ax * (by * cz - bz * cy) - ay * (bx * cz - bz * cx) + az * (bx * cy - by * cx))
[docs]
def cell_parameters(lattice: Lattice) -> tuple[float, float, float, float, float, float]:
"""Convert lattice vectors to the six CIF cell parameters.
Parameters
----------
lattice
The three cell vectors, in angstrom.
Returns
-------
tuple
``(a, b, c, alpha, beta, gamma)`` with lengths in angstrom and angles
in degrees. Following the crystallographic convention the CIF tags
name, ``alpha`` is the angle between ``b`` and ``c``, ``beta`` between
``a`` and ``c``, and ``gamma`` between ``a`` and ``b``.
Raises
------
ValueError
If any vector has zero length -- it has no direction, so the honest
answer is a refusal rather than a NaN angle.
Examples
--------
>>> cell_parameters(((3.0, 0.0, 0.0), (0.0, 3.0, 0.0), (0.0, 0.0, 5.0)))
(3.0, 3.0, 5.0, 90.0, 90.0, 90.0)
"""
lengths = [math.sqrt(sum(component * component for component in vector)) for vector in lattice]
for index, length in enumerate(lengths):
if length <= 0.0:
raise ValueError(f"lattice vector {index + 1} has zero length; the cell is degenerate")
def angle(first: int, second: int) -> float:
dot = sum(u * v for u, v in zip(lattice[first], lattice[second], strict=True))
cosine = min(1.0, max(-1.0, dot / (lengths[first] * lengths[second])))
return math.degrees(math.acos(cosine))
a, b, c = lengths
return (a, b, c, angle(1, 2), angle(0, 2), angle(0, 1))
[docs]
def cart_to_frac(position: Vector, lattice: Lattice) -> Vector:
"""Convert cartesian coordinates to fractional ones.
The inverse of ``extract_runner.frac_to_cart``.
Parameters
----------
position
An ``(x, y, z)`` triple in angstrom.
lattice
The three cell vectors to express it against.
Returns
-------
tuple
The fractional triple. Not wrapped into ``[0, 1)`` -- see the module
docstring.
Raises
------
ValueError
If the lattice vectors are linearly dependent, so the cell has no
volume and cannot be inverted.
"""
inverse = _invert3(lattice)
x, y, z = (sum(position[row] * inverse[row][column] for row in range(3)) for column in range(3))
return (x, y, z)
def _invert3(matrix: Lattice) -> Lattice:
"""Invert a 3x3 matrix by adjugate over determinant.
Parameters
----------
matrix
The rows to invert.
Returns
-------
tuple
The inverse, in the same row form.
Raises
------
ValueError
If the determinant is effectively zero -- a singular cell has no
volume.
"""
(a, b, c), (d, e, f), (g, h, i) = matrix
adjugate = (
(e * i - f * h, c * h - b * i, b * f - c * e),
(f * g - d * i, a * i - c * g, c * d - a * f),
(d * h - e * g, b * g - a * h, a * e - b * d),
)
determinant = a * adjugate[0][0] + b * adjugate[1][0] + c * adjugate[2][0]
if abs(determinant) < 1e-12:
raise ValueError("lattice vectors are linearly dependent; the cell has no volume")
return (
(adjugate[0][0] / determinant, adjugate[0][1] / determinant, adjugate[0][2] / determinant),
(adjugate[1][0] / determinant, adjugate[1][1] / determinant, adjugate[1][2] / determinant),
(adjugate[2][0] / determinant, adjugate[2][1] / determinant, adjugate[2][2] / determinant),
)
# ---------------------------------------------------------------------------
# CIF writing
# ---------------------------------------------------------------------------
def _snap(value: float) -> float:
return 0.0 if abs(value) < _ZERO_SNAP else value
def _fmt(value: float) -> str:
return f"{_snap(value):.8f}"
def _cif_value(value: Any) -> str:
"""Render one Python value as a CIF data value.
Parameters
----------
value
A string, number, bool, or ``None``.
Returns
-------
str
``?`` (CIF's honest unknown) for ``None`` and for NaN/infinity, a bare
literal for numbers, and quoted text otherwise. Single quotes inside a
string become double quotes -- a CIF value must not carry its own
delimiter, and provenance text is never worth failing an export over.
"""
if value is None:
return "?"
if isinstance(value, bool):
return "'true'" if value else "'false'"
if isinstance(value, int):
return str(value)
if isinstance(value, float):
return "?" if math.isnan(value) or math.isinf(value) else _fmt(value)
return "'" + str(value).replace("'", '"').replace("\n", " ") + "'"
[docs]
def geometry_to_cif(
geometry: Geometry,
*,
data_name: str,
extra: Mapping[str, Any] | None = None,
) -> str:
"""Serialize a :class:`Geometry` as a P1 CIF document.
Parameters
----------
geometry
The structure to write.
data_name
Names the data block -- a material id. CIF forbids whitespace here, and
an empty name too.
extra
Provenance items written as ``_mcli_<key>`` tags in the given order.
Keys must match ``[A-Za-z0-9_]+``; values may be strings, numbers,
bools, or ``None``.
Returns
-------
str
The complete CIF document, newline-terminated.
Raises
------
ValueError
For an empty or whitespace-bearing ``data_name``; a geometry whose
species and coordinate counts disagree or that holds no sites; a
degenerate cell; or a provenance key with illegal characters.
"""
if not data_name or any(character.isspace() for character in data_name):
raise ValueError(
f"CIF data block name must be non-empty and whitespace-free: {data_name!r}"
)
if len(geometry.species) != len(geometry.fractional):
raise ValueError(
f"geometry has {len(geometry.species)} species but "
f"{len(geometry.fractional)} coordinate triple(s)"
)
if not geometry.species:
raise ValueError("geometry holds no atom sites")
a, b, c, alpha, beta, gamma = cell_parameters(geometry.lattice)
lines = [
"# Generated by MatMelt -- DFTB+ relaxed geometry, P1 (no symmetry claimed).",
f"data_{data_name}",
f"_cell_length_a {_fmt(a)}",
f"_cell_length_b {_fmt(b)}",
f"_cell_length_c {_fmt(c)}",
f"_cell_angle_alpha {_fmt(alpha)}",
f"_cell_angle_beta {_fmt(beta)}",
f"_cell_angle_gamma {_fmt(gamma)}",
f"_cell_volume {_fmt(cell_volume(geometry.lattice))}",
f"_symmetry_space_group_name_H-M {_cif_value('P 1')}",
"_symmetry_Int_Tables_number 1",
f"_chemical_formula_sum {_cif_value(formula_sum(geometry.species))}",
]
for key, value in (extra or {}).items():
if not _EXTRA_KEY_RE.fullmatch(key):
raise ValueError(f"CIF provenance key must be [A-Za-z0-9_]+, got {key!r}")
lines.append(f"_mcli_{key}".ljust(33) + _cif_value(value))
lines += [
"",
"loop_",
"_symmetry_equiv_pos_as_xyz",
" 'x, y, z'",
"",
"loop_",
"_atom_site_type_symbol",
"_atom_site_label",
"_atom_site_occupancy",
"_atom_site_fract_x",
"_atom_site_fract_y",
"_atom_site_fract_z",
]
seen: dict[str, int] = {}
for symbol, (x, y, z) in zip(geometry.species, geometry.fractional, strict=True):
seen[symbol] = seen.get(symbol, 0) + 1
label = f"{symbol}{seen[symbol]}"
lines.append(f"{symbol:<4} {label:<6} 1.0000 {_fmt(x)} {_fmt(y)} {_fmt(z)}")
return "\n".join(lines) + "\n"
[docs]
def gen_to_cif(
gen_text: str,
*,
data_name: str,
extra: Mapping[str, Any] | None = None,
) -> str:
"""Parse DFTB+ ``.gen`` text and serialize it as a P1 CIF document.
The one call the ``cif`` command makes per material.
Parameters
----------
gen_text
The ``geom.out.gen`` contents DFTB+ relaxed to.
data_name
Names the CIF data block -- a material id.
extra
Provenance items, written as ``_mcli_<key>`` tags.
Returns
-------
str
The complete CIF document.
Raises
------
ValueError
From ``parse_gen`` for unreadable geometry, or from
``geometry_to_cif`` for an unwritable one.
"""
return geometry_to_cif(parse_gen(gen_text), data_name=data_name, extra=extra)