Fix cif writing issues

Author: rartino

src/httk/atomistic/_writing.py

@@ -1,6 +1,7 @@
 """Private exact neutral-payload serializers used by :func:`httk.core.save`."""
 
 import fractions
+import math
 from collections.abc import Mapping
 from decimal import Decimal, localcontext
 from typing import Any
@@ -86,6 +87,23 @@ def _snap_to_precision(value: fractions.Fraction, precision: Any) -> fractions.F
     return value.limit_denominator(round(1 / precision))
 
 
+def _coordinate_decimals(precision: Any) -> int | None:
+    """Decimal places whose CIF digit-count precision claim (``10**-places``) matches ``precision``.
+
+    A relaxed coordinate is only good to its stated precision. Padding it to the full 16 places
+    would tell the reader (via :func:`httk.core.decimal_precision`, which reads precision off the
+    digit count) that ``0.855`` is exact to ``1e-16``. Emitting ``round(-log10(precision))`` places
+    instead makes the written token claim the precision the structure actually has, so a CIF
+    round-trip recovers it. Returns ``None`` for exact structures (``precision is None``), which
+    keep the full-width padding.
+    """
+    if precision is None or precision <= 0:
+        return None
+    # floor, not round: a written token must never claim finer precision than the data has
+    # (httk reads precision as the coarsest of its numbers), so round the digit count down.
+    return max(1, min(16, math.floor(-math.log10(float(precision)))))
+
+
 def _poscar_token(value: Any) -> str:
     if isinstance(value, fractions.Fraction):
         rational = value
@@ -317,14 +335,20 @@ def append_row(species: Any, index: int, source_label: str, position: FracVector
     # Snap float-noise coordinates onto the short rational their precision resolves (0.355
     # rather than 563492063538/1587301587431). Only for the relaxed unit-cell path (snap=True
     # with a finite precision); the ASU path passes exact Wyckoff reps that must not be rounded.
-    coordinate_precision = structure.coordinate_precision if snap else None
+    # The digit-width claim follows the structure's precision unconditionally (so a relaxed cell
+    # survives repeated CIF round-trips), while value-snapping is only for the relaxed unit-cell
+    # path -- the ASU path passes exact Wyckoff reps that must not be rounded.
+    snap_precision = structure.coordinate_precision if snap else None
     return {
         "format": "cif",
         "approximate": cell_is_approximate,
+        # A finite precision (relaxed data) writes only as many coordinate decimals as it resolves,
+        # so the reader recovers that precision instead of reading full padding as machine-exact.
+        "coordinate_decimals": _coordinate_decimals(structure.coordinate_precision),
         "cell_parameters_exact": cell_parameters,
         "positions_exact": [
             tuple(
-                _render(_snap_to_precision(value, coordinate_precision), field="fractional coordinate")
+                _render(_snap_to_precision(value, snap_precision), field="fractional coordinate")
                 for value in row.to_fractions()
             )
             for row in written_positions
@@ -348,8 +372,10 @@ def _cif_payload_from_structure(obj: Any) -> Mapping[str, object]:
 
     Exact channels retain rational tokens (including ``1/3``); the CIF writer emits
     valid standard decimals and adds ``_httk_*_exact`` companions when those decimals
-    are lossy. Lossy standard values use 16 significant digits; finite coordinate
-    values are padded to 16 decimal places to keep recognition precision stable.
+    are lossy. Lossy standard values use 16 significant digits. An exact structure pads
+    its coordinate decimals to 16 places (reading back as machine-exact); a relaxed one
+    writes only the places its precision resolves (:func:`_coordinate_decimals`), so the
+    digit count round-trips as the precision the data actually has rather than 1e-16.
 
     A relaxed unit-cell structure carries a finite coordinate/cell precision but stores
     coordinates and cell lengths as full-precision rationalizations of floating-point input

src/httk/atomistic/io/cif/cif_writer.py

@@ -301,26 +301,30 @@ def _dual_cif_value(value: object) -> tuple[object, str | None]:
 
 
 def _dual_cif_column(
-    values: Iterable[object], *, pad_coordinates: bool = False
+    values: Iterable[object], *, pad_coordinates: bool = False, decimals: int = _HTTK_CIF_DECIMAL_DIGITS
 ) -> tuple[list[object], list[str | None], bool]:
     standard: list[object] = []
     exact: list[str | None] = []
     for value in values:
         visible, companion = _dual_cif_value(value)
         if pad_coordinates and visible is not None and companion is None:
-            visible = _pad_coordinate_decimal(visible)
+            visible = _pad_coordinate_decimal(visible, decimals)
         standard.append(visible)
         exact.append(companion)
     return standard, exact, any(value is not None for value in exact)
 
 
-def _pad_coordinate_decimal(value: object) -> object:
-    """Give standard coordinate decimals a fixed 16-place precision claim."""
+def _pad_coordinate_decimal(value: object, decimals: int = _HTTK_CIF_DECIMAL_DIGITS) -> object:
+    """Give a standard coordinate decimal a fixed ``decimals``-place precision claim.
+
+    ``decimals`` defaults to 16 (an exact structure's coordinate reads as machine-exact); a relaxed
+    structure passes its precision-derived width so the digit count states the precision it has.
+    """
     text = str(value)
     whole, dot, fraction = text.partition(".")
     if not dot:
-        return text + "." + "0" * _HTTK_CIF_DECIMAL_DIGITS
-    return whole + "." + fraction.ljust(_HTTK_CIF_DECIMAL_DIGITS, "0")
+        return text + "." + "0" * decimals
+    return whole + "." + fraction.ljust(decimals, "0")
 
 
 def _neutral_cif_block(block: Mapping[str, object]) -> dict[str, object]:
@@ -365,8 +369,13 @@ def _neutral_cif_block(block: Mapping[str, object]) -> dict[str, object]:
     raw["loop_atoms"] = loop_atoms
     raw["atom_site_label"] = labels
     raw["atom_site_type_symbol"] = symbols
+    # A relaxed block states how many coordinate decimals its precision resolves; exact blocks
+    # (and reader-round-tripped blocks that omit the key) keep the full 16-place claim.
+    coordinate_decimals = cast("int | None", block.get("coordinate_decimals")) or _HTTK_CIF_DECIMAL_DIGITS
     for index, tag in enumerate(("atom_site_fract_x", "atom_site_fract_y", "atom_site_fract_z")):
-        values, exact_values, has_exact = _dual_cif_column((row[index] for row in positions), pad_coordinates=True)
+        values, exact_values, has_exact = _dual_cif_column(
+            (row[index] for row in positions), pad_coordinates=True, decimals=coordinate_decimals
+        )
         raw[tag] = values
         if has_exact:
             raw[f"httk_{tag}_exact"] = exact_values

tests/test_save.py

@@ -159,9 +159,18 @@ def test_relaxed_cif_snaps_float_noise_to_clean_decimals(tmp_path):
     assert "/" not in "".join(text.splitlines()[12:])  # no fraction tokens in the atom/loop body
     # Clean, snapped values rather than the raw float noise.
     assert "_cell_length_a 5.568" in text
-    assert "0.3550000000000000" in text
-    assert "0.7590000000000000" in text
     assert "0.3549999999999969" not in text
+    # Coordinates are written only to the precision they resolve (~1e-4 -> 4 decimals), not padded
+    # to 16 places, so the digit count does not claim machine precision the data lacks.
+    assert "V V 0.0000 0.3550 0.2500" in text
+    assert "0.3550000000000000" not in text
+    # A CIF round-trip therefore recovers a realistic precision, not 1e-16 -- and stays there on a
+    # second round-trip (the width follows the structure's precision, not just the first-write path).
+    reloaded = load(destination)
+    assert reloaded.coordinate_precision == Fraction(1, 10000)
+    resaved = tmp_path / "relaxed2.cif"
+    save(reloaded, resaved)
+    assert load(resaved).coordinate_precision == Fraction(1, 10000)
 
 
 def test_relaxed_cif_strict_mode_refuses_snapped_cell(tmp_path):