Attach stated precision when building structures from CIF and POSCAR

Author: rartino

src/httk/atomistic/asu_recognition.py

@@ -243,7 +243,10 @@ def _group_into_orbits(
             consumed[other] = True
         asu_sites.append(ASUSite(position.letter, parameters, species))
 
-    return ASUStructure(cell, standard, asu_sites, view.species, transform)
+    # The recognized ASU inherits the precision of the structure it came from: nothing
+    # about recognition sharpens the data, and dropping it here would mean the value had
+    # to be guessed again by everything downstream.
+    return ASUStructure(cell, standard, asu_sites, view.species, transform, view.sites.precision)
 
 
 def _lies_in_orbit(point: FracVector, orbit: Sequence[FracVector], cell: Any, tolerance: float) -> bool:

src/httk/atomistic/asu_structure.py

@@ -28,6 +28,7 @@
 from httk.core import FracVector
 
 from . import data
+from ._vector_guards import to_precision
 from .cell import Cell
 from .cell_class_view import CellClassView
 from .cell_like import CellLike
@@ -86,6 +87,7 @@ class ASUStructure:
     _transform: SettingTransform
     _asu_sites: tuple[ASUSite, ...]
     _species: tuple[Species, ...]
+    _coordinate_precision: fractions.Fraction | None
 
     def __init__(
         self,
@@ -94,6 +96,7 @@ def __init__(
         asu_sites: Sequence[ASUSite],
         species: Sequence[SpeciesLike],
         transform: SettingTransform | None = None,
+        coordinate_precision: Any = None,
     ) -> None:
         self._cell = cell if isinstance(cell, Cell) else CellClassView(cell)
         self._spacegroup = spacegroup if isinstance(spacegroup, Spacegroup) else Spacegroup.standard(spacegroup)
@@ -104,6 +107,7 @@ def __init__(
                 f"express the difference as the transform instead"
             )
         self._transform = SettingTransform.identity() if transform is None else transform
+        self._coordinate_precision = to_precision(coordinate_precision)
         self._asu_sites = tuple(asu_sites)
         self._species = tuple(item if isinstance(item, Species) else SpeciesClassView(item) for item in species)
 
@@ -149,6 +153,20 @@ def species(self) -> tuple[Species, ...]:
         """The species referenced by the sites."""
         return self._species
 
+    @property
+    def coordinate_precision(self) -> fractions.Fraction | None:
+        """How precisely the coordinates behind this structure were stated, or ``None``.
+
+        Fractional, and expressed in **this structure's own setting** — the frame the data
+        arrived in — so it needs no transforming on the way to the expanded sites. Recording
+        it here is what lets an asymmetric unit say how good the data behind it was, rather
+        than leaving that to be guessed again downstream.
+
+        It is provenance, never an operating parameter: expansion remains exact and uses no
+        tolerance at all.
+        """
+        return self._coordinate_precision
+
     @property
     def asu(self) -> "ASUStructure":
         """Self, so a view can recognize a backend that already holds an ASU and pass it through."""
@@ -239,7 +257,7 @@ def expand_sites(self) -> Sites:
         The opposite case, a transform onto a larger cell, is covered by
         :meth:`~httk.atomistic.SettingTransform.lattice_cosets`.
         """
-        return Sites(self._expansion[0])
+        return Sites(self._expansion[0], self._coordinate_precision)
 
     def expand_species_at_sites(self) -> tuple[str, ...]:
         """The species name occupying each site produced by :meth:`expand_sites`, in order."""

src/httk/atomistic/asu_structure_view.py

@@ -65,6 +65,7 @@ def __new__(
             asu.asu_sites,
             asu.species,
             asu.transform,
+            asu.coordinate_precision,
         )
         instance._backend = backend
         return instance

src/httk/atomistic/cif_structures.py

@@ -121,7 +121,14 @@ def asu_structure_from_cif(
             )
         asu_sites.append(ASUSite(letter, parameters, name))
 
-    return ASUStructure(cell, standard, asu_sites, list(species_by_name.values()), transform)
+    return ASUStructure(
+        cell,
+        standard,
+        asu_sites,
+        list(species_by_name.values()),
+        transform,
+        data.get("coordinate_precision"),
+    )
 
 
 def cif_setting(data: Mapping[str, Any]) -> Spacegroup:
@@ -182,10 +189,16 @@ def _candidate_settings(data: Mapping[str, Any]) -> list[dict[str, Any]]:
 
 def _cell_from_cif(data: Mapping[str, Any]) -> Cell:
     """The cell, built exactly from the lattice parameters where the file gives them."""
+    precision = data.get("basis_precision")
+    exact = data.get("cell_parameters_exact")
+    if exact is not None and all(value is not None for value in exact):
+        # The text the file wrote, so 5.6402 becomes 56402/10000 rather than the binary
+        # value of float("5.6402").
+        return Cell(CellParams([fractions.Fraction(value) for value in exact]).basis, 1, precision)
     parameters = data.get("cell_parameters")
     if parameters is not None:
-        return Cell(CellParams([fractions.Fraction(str(value)) for value in parameters]).basis)
-    return Cell(data["basis"])
+        return Cell(CellParams([fractions.Fraction(str(value)) for value in parameters]).basis, 1, precision)
+    return Cell(data["basis"], 1, precision)
 
 
 def _exact_positions(data: Mapping[str, Any]) -> list[FracVector]:

src/httk/atomistic/vasp_structures.py

@@ -12,13 +12,16 @@
 """
 
 import fractions
+import math
 from collections.abc import Mapping
 from typing import Any, Callable
 
 from httk.core import SurdScalar, SurdVector, load
 from httk.core.vectors import exactmath
 
+from ._vector_guards import to_surdscalar
 from .cell import Cell
+from .sites import Sites
 from .species import Species
 from .structure import Structure
 
@@ -72,21 +75,24 @@ def structure_from_poscar(data: Mapping[str, Any]) -> Structure:
     scale_str = data.get("scale")
     volume_str = data.get("volume")
     if scale_str is not None:
-        cell = Cell(cell_rows, scale_str)
+        scale: Any = scale_str
     elif volume_str is not None:
         abs_det = abs(_to_fraction(raw_basis.det()))
         if abs_det == 0:
             raise ValueError("Cannot volume-scale a degenerate cell (zero determinant).")
         target_volume = _to_fraction(SurdVector.create(volume_str)._as_scalar())
-        cell = Cell(cell_rows, _cube_root(target_volume / abs_det))
+        scale = _cube_root(target_volume / abs_det)
     else:
-        cell = Cell(cell_rows)
+        scale = 1
+
+    cell = Cell(cell_rows, scale, _basis_precision(data, raw_basis, scale))
 
     if data["cartesian"]:
         # reduced = cart * basis^-1 (row-vector convention); the universal scale cancels.
         reduced: Any = SurdVector.create(data["coords"]) * raw_basis.inv()
     else:
         reduced = data["coords"]
+    sites = Sites(reduced, _coordinate_precision(data, raw_basis))
 
     symbols = data.get("symbols")
     if symbols is None:
@@ -105,7 +111,7 @@ def structure_from_poscar(data: Mapping[str, Any]) -> Structure:
             seen.add(symbol)
             species.append(Species(name=symbol, chemical_symbols=(symbol,), concentration=(1.0,)))
 
-    return Structure(cell, reduced, species, species_at_sites)
+    return Structure(cell, sites, species, species_at_sites)
 
 
 def _structure_from_cif(data: Mapping[str, Any]) -> Structure:
@@ -122,6 +128,55 @@ def _structure_from_cif(data: Mapping[str, Any]) -> Structure:
     return StructureSimpleView(structures[0])
 
 
+def _basis_precision(data: Mapping[str, Any], raw_basis: SurdVector, scale: Any) -> fractions.Fraction | None:
+    """How precisely the POSCAR states its cell, as an absolute length.
+
+    The reader reports the precision of the cell-vector tokens *as written*, before the
+    universal scaling factor multiplies them, so the conversion happens here where the
+    scale is known: ``scale * cell_precision``.
+
+    The scale's **own** written precision is deliberately not folded in, even though the
+    reader reports it. The scaling factor is a definition of units — a choice about whether
+    to carry the lattice constant in the factor or in the rows — and not a quantity measured
+    independently of them. Charging its digits as well would double-count the same
+    measurement, and would give a nonsensical answer in the ordinary case: a POSCAR that
+    writes ``1.0`` would be charged 0.1 of an uncertainty it plainly does not have, and the
+    cell of a 5.64 A structure would be declared good to only half an angstrom.
+
+    For the negative, volume-scaling form the scale is derived from a target volume rather
+    than written; the volume's precision is likewise not propagated through the cube root.
+    """
+    cell_precision = data.get("cell_precision")
+    if cell_precision is None:
+        return None
+    scale_value = abs(_to_fraction(to_surdscalar(scale)))
+    return scale_value * fractions.Fraction(cell_precision)
+
+
+def _coordinate_precision(data: Mapping[str, Any], raw_basis: SurdVector) -> fractions.Fraction | None:
+    """How precisely the POSCAR states its coordinates, in fractional units.
+
+    Direct coordinates are already fractional and pass straight through — and the universal
+    scaling factor cancels for them, exactly as it does for the coordinates themselves.
+
+    Cartesian coordinates are a length, so they are divided by the *shortest* cell edge to
+    become fractional. Shortest because that yields the largest fractional uncertainty,
+    which is the conservative choice; it is an approximation for a strongly oblique cell,
+    where a rigorous bound would need the norm of the inverse basis. The scale cancels here
+    too, since it multiplies the coordinates and the edge lengths alike.
+    """
+    precision = data.get("coordinate_precision")
+    if precision is None:
+        return None
+    if not data["cartesian"]:
+        return fractions.Fraction(precision)
+
+    shortest = min(math.dist((0.0, 0.0, 0.0), row) for row in raw_basis.to_floats())
+    if shortest <= 0:
+        return None
+    return fractions.Fraction(precision) / fractions.Fraction(str(shortest)).limit_denominator(10**12)
+
+
 #: Maps a loaded payload's ``"format"`` tag to the structure builder for it.
 _STRUCTURE_ADAPTERS: dict[str, Callable[[Mapping[str, Any]], Structure]] = {
     "vasp-poscar": structure_from_poscar,

tests/test_precision.py

@@ -9,6 +9,7 @@
 import fractions
 
 import pytest
+from httk.core import FracVector
 
 from httk.atomistic import (
     Cell,
@@ -24,6 +25,20 @@
     same_crystal,
 )
 
+pytest.importorskip("httk.io", reason="the readers live in httk-io")
+
+import io  # noqa: E402
+
+from httk.io.vasp import read_poscar  # noqa: E402
+
+from httk.atomistic import (  # noqa: E402
+    ASUSite,
+    ASUStructure,
+    Spacegroup,
+    recognize_asu,
+    structure_from_poscar,
+)
+
 F = fractions.Fraction
 
 CUBIC = [[5.0, 0, 0], [0, 5.0, 0], [0, 0, 5.0]]
@@ -191,3 +206,77 @@ def test_precision_does_not_affect_equality() -> None:
     vague = Structure(Cell(CUBIC), Sites([[0, 0, 0], [F(1, 2), F(1, 2), F(1, 2)]]), _species(), ["Na", "Na"])
     assert precise == vague
     assert same_crystal(precise, vague)
+
+
+# --- precision arriving from a file ---
+
+
+def _poscar(*, scale: str = "1.0", mode: str = "Direct", coords: str = "0.0000 0.0000 0.0000") -> str:
+    return (
+        "NaCl\n"
+        f"{scale}\n"
+        "5.6400 0.0000 0.0000\n0.0000 5.6400 0.0000\n0.0000 0.0000 5.6400\n"
+        "Na\n1\n"
+        f"{mode}\n{coords}\n"
+    )
+
+
+def test_poscar_direct_coordinates_pass_their_precision_through() -> None:
+    """Direct coordinates are already fractional, and the scale cancels for them."""
+    structure = structure_from_poscar(read_poscar(io.StringIO(_poscar(coords="0.5000 0.5000 0.5000"))))
+    assert structure.coordinate_precision == F(1, 10000)
+    assert structure.basis_precision == F(1, 10000)
+
+
+def test_poscar_cartesian_coordinates_are_converted_to_fractional() -> None:
+    """A Cartesian precision is a length; dividing by the shortest edge makes it fractional."""
+    structure = structure_from_poscar(
+        read_poscar(io.StringIO(_poscar(mode="Cartesian", coords="2.8200 2.8200 2.8200")))
+    )
+    assert float(structure.coordinate_precision) == pytest.approx(1e-4 / 5.64)
+
+
+def test_the_poscar_scale_multiplies_the_basis_precision() -> None:
+    """The cell entries are scaled, so their absolute precision scales with them."""
+    assert structure_from_poscar(read_poscar(io.StringIO(_poscar(scale="1.0")))).basis_precision == F(1, 10000)
+    assert structure_from_poscar(read_poscar(io.StringIO(_poscar(scale="2.0")))).basis_precision == F(1, 5000)
+
+
+def test_the_scales_own_digits_are_not_charged_as_uncertainty() -> None:
+    """The scaling factor defines units; it is not measured apart from the rows it scales.
+
+    Charging its digits would double-count the same measurement, and would declare the cell
+    of a 5.64 A structure good to only half an angstrom because the file wrote ``1.0``.
+    """
+    structure = structure_from_poscar(read_poscar(io.StringIO(_poscar(scale="1.0"))))
+    assert float(structure.basis_precision) == pytest.approx(1e-4)
+
+
+def test_an_asu_structure_carries_and_propagates_its_precision() -> None:
+    """Recorded in the structure's own setting, so it needs no transforming on the way out."""
+    asu = ASUStructure(
+        Cell(CUBIC, 1, BASIS_PRECISION),
+        225,
+        [ASUSite("a", FracVector.create(()), "Na")],
+        _species(),
+        coordinate_precision=COORD_PRECISION,
+    )
+    assert asu.coordinate_precision == COORD_PRECISION
+    assert asu.expand_sites().precision == COORD_PRECISION
+
+    expanded = StructureSimpleView(asu)
+    assert expanded.coordinate_precision == COORD_PRECISION
+    assert expanded.basis_precision == BASIS_PRECISION
+
+
+def test_recognition_carries_the_precision_onto_the_asu() -> None:
+    """Nothing about recognizing symmetry sharpens the data, so the claim is inherited."""
+    structure = Structure(
+        Cell(CUBIC, 1, BASIS_PRECISION),
+        Sites([[0, 0, 0], [F(1, 2), F(1, 2), F(1, 2)]], COORD_PRECISION),
+        _species() + [Species(name="Cl", chemical_symbols=("Cl",), concentration=(1.0,))],
+        ["Na", "Cl"],
+    )
+    recovered = recognize_asu(structure, setting=Spacegroup.standard(221))
+    assert recovered.coordinate_precision == COORD_PRECISION
+    assert recovered.cell.precision == BASIS_PRECISION