"""Convert CrystFEL reflection lists to SHELX HKLF 4 format."""
from __future__ import annotations
import math
from pathlib import Path
from typing import Iterable, Iterator, List, Optional, Tuple
__all__ = ["crystfel_to_hklf4", "main"]
def _compute_d_spacing(
h: int, k: int, l: int,
a: float, b: float, c: float,
al: float, be: float, ga: float,
) -> float:
"""
Compute d-spacing in Å for Miller indices (h, k, l) given unit cell
parameters a, b, c (Å) and al, be, ga (degrees). Uses the general
triclinic formula.
"""
ar = math.radians(al)
br = math.radians(be)
gr = math.radians(ga)
ca, cb, cg = math.cos(ar), math.cos(br), math.cos(gr)
sa, sb, sg = math.sin(ar), math.sin(br), math.sin(gr)
v = a * b * c * math.sqrt(
1 - ca**2 - cb**2 - cg**2 + 2 * ca * cb * cg
)
s11 = (b * c * sa) ** 2
s22 = (a * c * sb) ** 2
s33 = (a * b * sg) ** 2
s12 = a * b * c**2 * (ca * cb - cg)
s23 = a**2 * b * c * (cb * cg - ca)
s13 = a * b**2 * c * (ca * cg - cb)
d_star_sq = (
s11 * h**2 + s22 * k**2 + s33 * l**2
+ 2 * s12 * h * k + 2 * s23 * k * l + 2 * s13 * h * l
) / v**2
if d_star_sq <= 0:
return float("inf")
return 1.0 / math.sqrt(d_star_sq)
def _is_table_header(line: str) -> bool:
"""Return True when the line looks like the CrystFEL column header."""
lower = line.lower().strip()
return lower.startswith("h") and "sigma" in lower and "nmeas" in lower
def _last_float(tokens: List[str]) -> Optional[float]:
"""Return the last token that can be parsed as float, if any."""
for tok in reversed(tokens):
try:
return float(tok)
except ValueError:
continue
return None
def _iter_reflections(lines: Iterable[str]) -> Iterator[Tuple[int, int, int, float, float]]:
"""Yield (h, k, l, I, sigI) tuples from a CrystFEL .hkl text stream."""
in_table = False
for raw in lines:
line = raw.strip()
if not line:
continue
lower = line.lower()
if not in_table:
if _is_table_header(lower):
in_table = True
continue
if lower.startswith("end of reflections"):
break
parts = line.split()
if len(parts) < 5:
continue
try:
h, k, l = map(int, parts[0:3])
intensity = float(parts[3])
except ValueError:
continue
sigma_tokens = parts[4:]
# Drop a trailing nmeas column if it looks integer-like.
if sigma_tokens:
try:
int(sigma_tokens[-1])
if len(sigma_tokens) > 1:
sigma_tokens = sigma_tokens[:-1]
except ValueError:
pass
sigma = _last_float(sigma_tokens)
if sigma is None:
continue
yield h, k, l, intensity, sigma
[docs]
def crystfel_to_hklf4(
in_path: str | Path,
out_path: str | Path,
*,
d_min: float | None = None,
cell: tuple[float, float, float, float, float, float] | None = None,
) -> int:
"""
Convert a CrystFEL reflection list ("crystfel.hkl") into a SHELX HKLF 4 file.
The output uses the canonical FORTRAN format ``3I4, 2F8.2, I4``:
- h, k, l as 4-char integers
- intensity and sigma as 8-char floats with 2 decimals
- batch number as 4-char integer (always 0)
The file contains only reflection records followed by a mandatory
zero-terminator record (h=k=l=0). ``HKLF 4`` belongs in the .ins
instruction file, not in the .hkl reflection file.
Parameters
----------
in_path:
Source CrystFEL .hkl file.
out_path:
Destination HKLF 4 path.
d_min:
Optional resolution cutoff in Å. Reflections with d-spacing below
this value are excluded. Requires ``cell`` to be provided.
cell:
Unit cell parameters ``(a, b, c, alpha, beta, gamma)`` in Å and
degrees. Required when ``d_min`` is set.
Returns
-------
int
Number of reflections written.
"""
src = Path(in_path)
dst = Path(out_path)
if not src.is_file():
raise FileNotFoundError(src)
if d_min is not None and cell is None:
raise ValueError("cell parameters are required when d_min is set")
reflections = []
with src.open("r", encoding="utf-8", errors="replace") as fin:
for h, k, l, intensity, sigma in _iter_reflections(fin):
if d_min is not None and cell is not None:
d = _compute_d_spacing(h, k, l, *cell)
if d < d_min:
continue
reflections.append((h, k, l, intensity, sigma))
if not reflections:
raise ValueError("No reflections were parsed from the CrystFEL file.")
max_abs = max(
max(abs(float(intensity)), abs(float(sigma)))
for _h, _k, _l, intensity, sigma in reflections
)
scale = 1.0
if max_abs > 99999.99:
scale = 99999.99 / max_abs
with dst.open("w", encoding="ascii", newline="\n") as fout:
for h, k, l, intensity, sigma in reflections:
fout.write(f"{h:4d}{k:4d}{l:4d}{intensity * scale:8.2f}{sigma * scale:8.2f}{0:4d}\n")
# Mandatory zero-terminator record
fout.write(f"{0:4d}{0:4d}{0:4d}{0.0:8.2f}{0.0:8.2f}{0:4d}\n")
return len(reflections)
[docs]
def main(argv: Optional[List[str]] = None) -> int:
"""Simple CLI wrapper matching the prototype script usage."""
import argparse
parser = argparse.ArgumentParser(
description="Convert CrystFEL .hkl to SHELX HKLF 4",
)
parser.add_argument("input", help="Path to CrystFEL reflection list (crystfel.hkl)")
parser.add_argument("output", help="Destination HKLF 4 file")
parser.add_argument(
"--overwrite",
action="store_true",
help="Allow replacing an existing output file",
)
args = parser.parse_args(argv)
out_path = Path(args.output)
if out_path.exists() and not args.overwrite:
parser.error(f"{out_path} already exists; use --overwrite to replace it")
written = crystfel_to_hklf4(args.input, out_path)
parser.exit(0, f"Wrote {written} reflections to {out_path}\n")
if __name__ == "__main__": # pragma: no cover - convenience entrypoint
raise SystemExit(main())