rugnux: also write merged reflections as SHELX HKLF-4 (.hkl)

WriteReflections now emits a third format alongside the MTZ and mmCIF: a
SHELX HKLF-4 text file (<prefix>.hkl), the direct input for SHELXC / ANODE /
SHELXD. Fixed FORMAT(3I4,2F8.2), one record per reflection as h k l I sig(I),
ending with the 0 0 0 terminator.

Two things make it a valid substructure-solution input:
- Bijvoet mates are written separately (I(+) at +hkl, I(-) at -hkl) so the
  anomalous differences survive; a reflection with no anomalous split is
  written once as its mean.
- Intensities are put on a common scale so the largest value fits the fixed
  F8.2 field (negative weak intensities would otherwise overflow it). I and
  sigma share the scale, so the anomalous signal is untouched; the absolute
  scale is irrelevant to SHELXC/ANODE, which use only ratios.

The MTZ writer's Bijvoet row-building is factored into a shared BuildMergedRows
helper so the MTZ and HKL outputs agree (MTZ output is byte-equivalent to
before). docs/RUGNUX.md updated to describe the three reflection formats.

Validated end-to-end on lysozyme: SHELXC reads the .hkl and reports real
anomalous signal (<d"/sig> 1.6 at low res), and ANODE locates the Met/Cys S
and a chloride at 6.5-9 sigma from the unmodified file.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

image_analysis/WriteReflections.cpp

@@ -59,6 +59,84 @@ std::string CifStr(const std::string& s) {
     return s;
 }
 
+// One output row per reflection in the standard CCP4 anomalous layout: the merged mean (IMEAN / F)
+// plus the two Bijvoet mates (I(+)/I(-), F(+)/F(-)). The merge keeps the two mates as separate rows
+// (I+ under the ASU representative hkl, I- under -hkl); this collapses them into one row so the MTZ
+// and SHELX writers share one row list. has_anom is set false when NO reflection carries an anomalous
+// split (e.g. the stills path) - callers then omit the +/- columns.
+struct MergedOutRow {
+    int h = 0, k = 0, l = 0;
+    float Imean = NAN, sImean = NAN, Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN;
+    float Fmean = NAN, sFmean = NAN, Fp = NAN, sFp = NAN, Fm = NAN, sFm = NAN;
+    int rfree = 0;
+};
+
+std::vector<MergedOutRow> BuildMergedRows(const std::vector<MergedReflection> &reflections,
+                                          const DiffractionExperiment &experiment,
+                                          bool &has_anom) {
+    std::vector<MergedOutRow> out_rows;
+    has_anom = true;
+
+    if (experiment.GetScalingSettings().GetMergeFriedel()) {
+        // Friedel-merged: IMEAN is the already-merged intensity (r.I). I(+)/I(-) are carried verbatim
+        // from the Bijvoet split the merge kept (rotation always does; scaled non-anomalously), so a weak
+        // anomalous signal is preserved without reprocessing. A reflection with only one mate, or a
+        // centric, gets a missing value (NaN) for the absent hand. When NO reflection has an anomalous
+        // split (e.g. the stills path, which does not compute one) the anomalous columns are omitted.
+        has_anom = std::any_of(reflections.begin(), reflections.end(),
+            [](const MergedReflection& r){ return std::isfinite(r.I_plus) || std::isfinite(r.I_minus); });
+        out_rows.reserve(reflections.size());
+        for (const auto& r : reflections)
+            out_rows.push_back({r.h, r.k, r.l, r.I, r.sigma, r.I_plus, r.sigma_plus, r.I_minus,
+                                r.sigma_minus, r.F, r.sigmaF, r.F_plus, r.sigmaF_plus, r.F_minus,
+                                r.sigmaF_minus, r.rfree_flag ? 1 : 0});
+    } else {
+        // Anomalous: group the two mates by their (shared) Friedel-merged ASU representative, then form
+        // IMEAN / F as their inverse-variance Friedel mean. A single generator gives both the group key
+        // (its hkl, identical for +hkl and -hkl) and which mate this row is (.plus).
+        const HKLKeyGenerator key_gen(false, experiment.GetSpaceGroupNumber().value_or(1));
+        struct AnomRow {
+            int h = 0, k = 0, l = 0;
+            float Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN;
+            float Fp = NAN, sFp = NAN, Fm = NAN, sFm = NAN;
+            int rfree = 0;
+        };
+        std::map<std::tuple<int, int, int>, AnomRow> rows;
+        for (const auto& r : reflections) {
+            const HKLKey key = key_gen(r);
+            AnomRow& row = rows[{key.h, key.k, key.l}];
+            row.h = key.h; row.k = key.k; row.l = key.l;
+            row.rfree = r.rfree_flag ? 1 : 0;
+            if (key.plus) { row.Ip = r.I; row.sIp = r.sigma; row.Fp = r.F; row.sFp = r.sigmaF; }
+            else          { row.Im = r.I; row.sIm = r.sigma; row.Fm = r.F; row.sFm = r.sigmaF; }
+        }
+
+        // Friedel-mean of the two mates by inverse variance (the single mate, if only one was measured).
+        const auto combine = [](float a, float sa, float b, float sb, float& val, float& sig) {
+            const bool ok_a = std::isfinite(a) && sa > 0.0f;
+            const bool ok_b = std::isfinite(b) && sb > 0.0f;
+            if (ok_a && ok_b) {
+                const double wa = 1.0 / (static_cast<double>(sa) * sa);
+                const double wb = 1.0 / (static_cast<double>(sb) * sb);
+                val = static_cast<float>((wa * a + wb * b) / (wa + wb));
+                sig = static_cast<float>(1.0 / std::sqrt(wa + wb));
+            } else if (ok_a) { val = a; sig = sa; }
+            else if (ok_b)   { val = b; sig = sb; }
+            else             { val = NAN; sig = NAN; }
+        };
+
+        out_rows.reserve(rows.size());
+        for (const auto& [hkl, row] : rows) {
+            float i_mean, sig_i_mean, f_mean, sig_f_mean;
+            combine(row.Ip, row.sIp, row.Im, row.sIm, i_mean, sig_i_mean);
+            combine(row.Fp, row.sFp, row.Fm, row.sFm, f_mean, sig_f_mean);
+            out_rows.push_back({row.h, row.k, row.l, i_mean, sig_i_mean, row.Ip, row.sIp, row.Im,
+                                row.sIm, f_mean, sig_f_mean, row.Fp, row.sFp, row.Fm, row.sFm, row.rfree});
+        }
+    }
+    return out_rows;
+}
+
 } // namespace
 
 
@@ -286,77 +364,10 @@ void WriteMtzReflections(const std::vector<MergedReflection> &reflections,
 
     const int dataset_id = ds.id;
 
-    // The merge keeps the two Friedel mates as separate rows (I+ under the ASU representative hkl, I-
-    // under -hkl). Both output modes collapse them into one row per reflection with the standard CCP4
-    // anomalous layout (IMEAN + I(+)/I(-), and the same split for the French-Wilson amplitude), which
-    // aimless / ctruncate / mtz2sca / ANODE read directly. The modes differ only in how each row is
-    // formed, so they build one common row list and share the column/row emission below.
-    struct OutRow {
-        int h = 0, k = 0, l = 0;
-        float Imean = NAN, sImean = NAN, Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN;
-        float Fmean = NAN, sFmean = NAN, Fp = NAN, sFp = NAN, Fm = NAN, sFm = NAN;
-        int rfree = 0;
-    };
-    std::vector<OutRow> out_rows;
+    // One row per reflection in the CCP4 anomalous layout (IMEAN + I(+)/I(-), and the same split for
+    // the French-Wilson amplitude), which aimless / ctruncate / mtz2sca / ANODE read directly.
     bool has_anom = true;
-
-    if (experiment.GetScalingSettings().GetMergeFriedel()) {
-        // Friedel-merged: IMEAN is the already-merged intensity (r.I). I(+)/I(-) are carried verbatim
-        // from the Bijvoet split the merge kept (rotation always does; scaled non-anomalously), so a weak
-        // anomalous signal is preserved without reprocessing. A reflection with only one mate, or a
-        // centric, gets a missing value (NaN) for the absent hand. When NO reflection has an anomalous
-        // split (e.g. the stills path, which does not compute one) the anomalous columns are omitted.
-        has_anom = std::any_of(reflections.begin(), reflections.end(),
-            [](const MergedReflection& r){ return std::isfinite(r.I_plus) || std::isfinite(r.I_minus); });
-        out_rows.reserve(reflections.size());
-        for (const auto& r : reflections)
-            out_rows.push_back({r.h, r.k, r.l, r.I, r.sigma, r.I_plus, r.sigma_plus, r.I_minus,
-                                r.sigma_minus, r.F, r.sigmaF, r.F_plus, r.sigmaF_plus, r.F_minus,
-                                r.sigmaF_minus, r.rfree_flag ? 1 : 0});
-    } else {
-        // Anomalous: group the two mates by their (shared) Friedel-merged ASU representative, then form
-        // IMEAN / F as their inverse-variance Friedel mean. A single generator gives both the group key
-        // (its hkl, identical for +hkl and -hkl) and which mate this row is (.plus).
-        const HKLKeyGenerator key_gen(false, experiment.GetSpaceGroupNumber().value_or(1));
-        struct AnomRow {
-            int h = 0, k = 0, l = 0;
-            float Ip = NAN, sIp = NAN, Im = NAN, sIm = NAN;
-            float Fp = NAN, sFp = NAN, Fm = NAN, sFm = NAN;
-            int rfree = 0;
-        };
-        std::map<std::tuple<int, int, int>, AnomRow> rows;
-        for (const auto& r : reflections) {
-            const HKLKey key = key_gen(r);
-            AnomRow& row = rows[{key.h, key.k, key.l}];
-            row.h = key.h; row.k = key.k; row.l = key.l;
-            row.rfree = r.rfree_flag ? 1 : 0;
-            if (key.plus) { row.Ip = r.I; row.sIp = r.sigma; row.Fp = r.F; row.sFp = r.sigmaF; }
-            else          { row.Im = r.I; row.sIm = r.sigma; row.Fm = r.F; row.sFm = r.sigmaF; }
-        }
-
-        // Friedel-mean of the two mates by inverse variance (the single mate, if only one was measured).
-        const auto combine = [](float a, float sa, float b, float sb, float& val, float& sig) {
-            const bool ok_a = std::isfinite(a) && sa > 0.0f;
-            const bool ok_b = std::isfinite(b) && sb > 0.0f;
-            if (ok_a && ok_b) {
-                const double wa = 1.0 / (static_cast<double>(sa) * sa);
-                const double wb = 1.0 / (static_cast<double>(sb) * sb);
-                val = static_cast<float>((wa * a + wb * b) / (wa + wb));
-                sig = static_cast<float>(1.0 / std::sqrt(wa + wb));
-            } else if (ok_a) { val = a; sig = sa; }
-            else if (ok_b)   { val = b; sig = sb; }
-            else             { val = NAN; sig = NAN; }
-        };
-
-        out_rows.reserve(rows.size());
-        for (const auto& [hkl, row] : rows) {
-            float i_mean, sig_i_mean, f_mean, sig_f_mean;
-            combine(row.Ip, row.sIp, row.Im, row.sIm, i_mean, sig_i_mean);
-            combine(row.Fp, row.sFp, row.Fm, row.sFm, f_mean, sig_f_mean);
-            out_rows.push_back({row.h, row.k, row.l, i_mean, sig_i_mean, row.Ip, row.sIp, row.Im,
-                                row.sIm, f_mean, sig_f_mean, row.Fp, row.sFp, row.Fm, row.sFm, row.rfree});
-        }
-    }
+    const std::vector<MergedOutRow> out_rows = BuildMergedRows(reflections, experiment, has_anom);
 
     mtz.add_column("H", 'H', dataset_id, -1, false);
     mtz.add_column("K", 'H', dataset_id, -1, false);
@@ -406,6 +417,49 @@ void WriteMtzReflections(const std::vector<MergedReflection> &reflections,
     mtz.write_to_file(filename);
 }
 
+void WriteShelxHklReflections(const std::vector<MergedReflection> &reflections,
+                              const DiffractionExperiment &experiment,
+                              const std::string &filename) {
+    bool has_anom = true;
+    const std::vector<MergedOutRow> rows = BuildMergedRows(reflections, experiment, has_anom);
+
+    // SHELX HKLF 4 (SHELXC / ANODE input): fixed FORMAT(3I4,2F8.2), one record per reflection as
+    // h k l I sigma(I). The Bijvoet mates are written separately - I(+) at +hkl, I(-) at -hkl - so the
+    // anomalous differences survive; a reflection with no anomalous split is written once as its mean.
+    // Intensities are put on a common scale so the largest value fits the F8.2 field (the absolute scale
+    // is irrelevant to SHELXC / ANODE, which use only ratios); I and sigma share the scale, so the
+    // anomalous signal is untouched. The file ends with a 0 0 0 terminator record.
+    const auto usable = [](float v, float s) { return std::isfinite(v) && std::isfinite(s) && s > 0.0f; };
+
+    double max_abs = 0.0;
+    for (const auto& r : rows) {
+        if (usable(r.Ip, r.sIp)) max_abs = std::max({max_abs, std::fabs(double(r.Ip)), double(r.sIp)});
+        if (usable(r.Im, r.sIm)) max_abs = std::max({max_abs, std::fabs(double(r.Im)), double(r.sIm)});
+        if (!usable(r.Ip, r.sIp) && !usable(r.Im, r.sIm) && usable(r.Imean, r.sImean))
+            max_abs = std::max({max_abs, std::fabs(double(r.Imean)), double(r.sImean)});
+    }
+    const double scale = (std::isfinite(max_abs) && max_abs > 0.0) ? 9999.0 / max_abs : 1.0;
+
+    std::ofstream out(filename);
+    if (!out)
+        throw std::runtime_error("WriteShelxHklReflections: cannot open " + filename);
+    out << std::fixed << std::setprecision(2);
+    const auto emit = [&out, scale](int h, int k, int l, float I, float sigma) {
+        out << std::setw(4) << h << std::setw(4) << k << std::setw(4) << l
+            << std::setw(8) << scale * I << std::setw(8) << scale * sigma << "\n";
+    };
+    for (const auto& r : rows) {
+        const bool plus = usable(r.Ip, r.sIp);
+        const bool minus = usable(r.Im, r.sIm);
+        if (plus)  emit(r.h, r.k, r.l, r.Ip, r.sIp);
+        if (minus) emit(-r.h, -r.k, -r.l, r.Im, r.sIm);
+        if (!plus && !minus && usable(r.Imean, r.sImean))
+            emit(r.h, r.k, r.l, r.Imean, r.sImean);
+    }
+    emit(0, 0, 0, 0.0f, 0.0f);   // HKLF-4 end-of-data marker
+    out.close();
+}
+
 void WriteReflections(const std::vector<MergedReflection> &reflections,
                       const UnitCell &unitCell,
                       const DiffractionExperiment &experiment,
@@ -413,8 +467,10 @@ void WriteReflections(const std::vector<MergedReflection> &reflections,
                       const std::string &isa,
                       const TwinningAnalysisResult &twinning,
                       const std::string &filename) {
-    // Always write both an MTZ and an mmCIF - each has its uses downstream (MTZ for the CCP4 /
-    // phenix reflection tools, mmCIF for deposition and as the self-describing native format).
+    // Write an MTZ, an mmCIF and a SHELX HKLF-4 .hkl - each has its uses downstream (MTZ for the CCP4 /
+    // phenix reflection tools, mmCIF for deposition and as the self-describing native format, HKLF-4 as
+    // the SHELXC / ANODE substructure-solution input).
     WriteMtzReflections(reflections, unitCell, experiment, filename + ".mtz");
     WriteMmcifReflections(reflections, unitCell, experiment, statistics, isa, twinning, filename + ".cif");
+    WriteShelxHklReflections(reflections, experiment, filename + ".hkl");
 }