diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md
index 9bedef3c1..698c3255e 100644
--- a/docs/RUGNUX_ADVANCED.md
+++ b/docs/RUGNUX_ADVANCED.md
@@ -353,6 +353,7 @@ Mode — `--mode <name>` (default `mx`):
 | `azint` | Only azimuthal integration (no spot finding/indexing); writes `<prefix>_process.h5` |
 | `scale` | Only re-scale/merge the already-integrated reflections in the input `_process.h5` (no re-integration) |
 | `calibration` | Determine the detector geometry from powder rings; writes `<prefix>.poni` and `<prefix>.json` |
+| `refine` | Refine the `--model` against merged structure factors (an MTZ or a structure-factor mmCIF) and process nothing; see [Refining a model against structure factors](RUGNUX_MODEL_AND_REFERENCE.md#refining-a-model-against-structure-factors-rugnux---mode-refine) (`--refine-score-only`, `--refine-d-min`) |
 
 Calibration (`--mode calibration`):
 
diff --git a/docs/RUGNUX_MODEL_AND_REFERENCE.md b/docs/RUGNUX_MODEL_AND_REFERENCE.md
index ff52d055a..9b6d58b11 100644
--- a/docs/RUGNUX_MODEL_AND_REFERENCE.md
+++ b/docs/RUGNUX_MODEL_AND_REFERENCE.md
@@ -251,6 +251,44 @@ reflections, and a refinement against new data then moves R-free up from the unr
 down. `R_WORK`/`R_FREE` of the model as given are never replaced. On a GPU a refinement takes seconds; the CPU
 engine computes the same and is several times slower.
 
+**How a model is scored.** R-work, R-free and the likelihood target of every stage - the model as given
+(`REFINE_R_WORK_START`) as much as the refined one - come from one model of the structure factors,
+F<sub>model</sub> = k exp(-s<sup>T</sup>Bs/4) (F<sub>calc</sub> + k<sub>sol</sub> exp(-B<sub>sol</sub>s<sup>2</sup>/4) F<sub>mask</sub>):
+the atoms with the dictionary's riding hydrogens and their B made isotropic; a flat bulk solvent masked with the
+REFMAC radii (probe 1.0 A, shrink 0.8 A) on a grid of its own, at most 0.4 A apart, so that the shrink acts at any
+resolution; k<sub>sol</sub> and B<sub>sol</sub> searched over the physical range on the working reflections, and the
+overall scale and anisotropic B fitted to them by robust least squares, so that a few reflections the model cannot
+explain (a beam-stop shadow, an ice ring, a missing part of the model) do not set the scale of the rest. On a set of
+fourteen depositions from 1.1 to 3.6 A this gives REFMAC 5.8's 0-cycle R-work and R-free on the same model and data
+within 0.002 on twelve, and within 0.006 on the other two (`tools/battery/refine_parity.py` runs the comparison;
+`--refine-score-only` prints each term under `REFINE_SCORE_`, see [`--mode refine`](#refining-a-model-against-structure-factors-rugnux---mode-refine)).
+
+## Refining a model against structure factors (`rugnux --mode refine`)
+
+`--mode refine` refines a model against merged structure factors that already exist - processed elsewhere, a
+deposition's, or a rugnux `.mtz` - and processes nothing:
+
+```
+rugnux --mode refine --model model.cif -o refined data.mtz
+```
+
+The input is an MTZ or a structure-factor mmCIF (as the PDB distributes them), gzipped or not, recognised by
+content; the model is PDB or mmCIF, taken as placed in the data's cell and refined in the data's space group. The
+columns are found, not named: amplitudes with their sigmas (F/SIGF, FP/SIGFP, FOBS, ...; calculated amplitudes and
+map coefficients such as FC, FWT, DELFWT are passed over) where the file has them, else intensities (IMEAN/SIGIMEAN,
+I/SIGI, ...) made into amplitudes by French-Wilson, else the mean of the anomalous pairs. The free set is read from
+a FreeR_flag-like column (flag 0 is the test set; where 0 marks the majority, as phenix and CNS write it, the
+minority value is); a file without one gets rugnux's own 5 % set. The refinement is the one `--refine` runs, and
+writes the same files - `<prefix>_refined_model.cif`, `<prefix>_refined_maps.mtz`, `<prefix>_blobs.json` - with
+the `REFINE_` keys in `<prefix>_report.txt`, beside `REFINE_FRENCH_WILSON=` and `REFINE_FREE_FLAGS= FILE|GENERATED`
+saying what had to be made of the input.
+
+| Option | Description |
+| --- | --- |
+| `--model <file>` | The model (required) |
+| `--refine-score-only` | Score the model as given and refine nothing; every term of the score - R-work and R-free overall and in 20 resolution shells, the overall scale, B and anisotropic B, k<sub>sol</sub> and B<sub>sol</sub>, the likelihood over the working and the free reflections and its figure of merit, and the geometry's rmsd and rmsZ - under `REFINE_SCORE_` keys (`--developer` prints them all) |
+| `--refine-d-min <A>` | Leave out reflections finer than this |
+
 ## Free set and model bias
 
 **How the free set is made.** Rugnux draws its own test set (5 %, raised towards about 500 reflections on small data but never above 10 %):
diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt
index e0696c1b3..ff903aef3 100644
--- a/rugnux/CMakeLists.txt
+++ b/rugnux/CMakeLists.txt
@@ -21,6 +21,8 @@ ADD_LIBRARY(JFJochRugnux STATIC
         RugnuxDefaults.h
         RugnuxLive.cpp
         RugnuxLive.h
+        RugnuxRefine.cpp
+        RugnuxRefine.h
         Rugnux.h
         RugnuxCommandLine.cpp
         RugnuxCommandLine.h
diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp
index e25587474..c8bdd3636 100644
--- a/rugnux/ResultReport.cpp
+++ b/rugnux/ResultReport.cpp
@@ -415,25 +415,10 @@ namespace {
         }
         return 0.0;
     }
-}
 
-// ---------------------------------------------------------------------------- the emitter
-// Populates the document. Nothing here writes to a stream: the order of the FILE is the renderer's
-// business, which is what lets the SUMMARY sit above the sections whose warnings it reports.
-ReportDocument BuildReportDocument(const std::string &output_prefix,
-                                   const std::string &input_file,
-                                   const DiffractionExperiment &experiment,
-                                   const ProcessResult &result,
-                                   const RunProvenance &provenance) {
-    ReportDocument doc;
-    const bool rotation = experiment.IsRotationIndexing();
-    const bool merged = result.has_merge_statistics;
-    doc.warnings = result.warnings;
-    // SUMMARY rows filled by the sections that measure them; empty leaves the row out.
-    std::string beam_centre_row, overload_row, overlap_row;
-
-    // ------------------------------------------------------------------ header block
-    {
+    // The header block of every report: the version, the build, the input and how the run was invoked.
+    ReportSection HeaderSection(const std::string &output_prefix, const std::string &input_file,
+                                const RunProvenance &provenance) {
         ReportSection h;
         Add(h, KeyInt("REPORT_VERSION", REPORT_VERSION));
         Add(h, KeyText("RUGNUX_VERSION", jfjoch_version()));
@@ -467,9 +452,136 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
             if (!provenance.gpu_description.empty())
                 Add(h, KeyText("GPU", provenance.gpu_description));
         }
-        doc.sections.push_back(std::move(h));
+        return h;
     }
 
+    // The refinement's keys (Refinement.h): REFINE_ of --refine and of --mode refine, with every term of the
+    // score of the model as given under REFINE_SCORE_.
+    void AddRefinementKeys(ReportSection &s, const RefinementResult &rf, const std::string &output_prefix) {
+        Add(s, Blank());
+        if (!rf.ok) {
+            Add(s, KeyEnum("REFINE", "NOT_PERFORMED", {"PERFORMED", "NOT_PERFORMED"}));
+            Add(s, KeyText("REFINE_REASON", rf.failure_reason));
+        } else {
+            Add(s, KeyEnum("REFINE", "PERFORMED", {"PERFORMED", "NOT_PERFORMED"}));
+            Add(s, KeyText("REFINE_DEVICE", rf.on_gpu ? "GPU" : "CPU", true));
+            Add(s, KeyReal("REFINE_SECONDS", rf.seconds, "{:.1f}", true));
+            Add(s, KeyReal("REFINE_D_MIN", rf.d_min, "{:.2f}"));
+            Add(s, KeyReal("REFINE_R_WORK_START", rf.r_work_start, "{:.4f}"));
+            Add(s, KeyReal("REFINE_R_FREE_START", rf.r_free_start, "{:.4f}"));
+            Add(s, KeyInt("REFINE_RIGID_BODIES", rf.rigid_bodies));
+            Add(s, KeyReal("REFINE_R_WORK_PLACEMENT", rf.r_work_placement, "{:.4f}"));
+            Add(s, KeyReal("REFINE_R_FREE_PLACEMENT", rf.r_free_placement, "{:.4f}"));
+            if (rf.restrained) {
+                Add(s, KeyReal("REFINE_R_WORK", rf.r_work, "{:.4f}"));
+                Add(s, KeyReal("REFINE_R_FREE", rf.r_free, "{:.4f}"));
+                Add(s, KeyReal("REFINE_RMSD_BONDS", rf.geometry.rmsd_bonds, "{:.4f}"));
+                Add(s, KeyReal("REFINE_RMSD_ANGLES", rf.geometry.rmsd_angles, "{:.2f}"));
+                Add(s, KeyReal("REFINE_RMSD_BONDS_START", rf.geometry_start.rmsd_bonds, "{:.4f}", true));
+                Add(s, KeyReal("REFINE_RMSD_ANGLES_START", rf.geometry_start.rmsd_angles, "{:.2f}", true));
+                Add(s, KeyInt("REFINE_CLASHES", rf.geometry.clashes));
+                Add(s, KeyReal("REFINE_XRAY_WEIGHT", rf.weight, "{:.4g}", true));
+                Add(s, KeyInt("REFINE_ATOMS", rf.atoms));
+                Add(s, KeyInt("REFINE_ATOMS_REFINED", rf.refined_atoms));
+                Add(s, KeyInt("REFINE_RESIDUES_FIXED", rf.fixed_residues));
+                if (!rf.fixed_names.empty()) {
+                    std::string names;
+                    for (const std::string &n : rf.fixed_names)
+                        names += (names.empty() ? "" : " ") + n;
+                    Add(s, KeyText("REFINE_RESIDUES_FIXED_NAMES", names));
+                }
+            }
+            Add(s, KeyReal("REFINE_BULK_SOLVENT_K_SOL", rf.k_sol, "{:.3f}", true));
+            Add(s, KeyReal("REFINE_BULK_SOLVENT_B_SOL", rf.b_sol, "{:.1f}", true));
+            int kinds[5] = {0, 0, 0, 0, 0};
+            for (const DensityBlob &b : rf.blobs)
+                kinds[static_cast<int>(b.kind)]++;
+            Add(s, KeyInt("REFINE_BLOBS", static_cast<int64_t>(rf.blobs.size())));
+            Add(s, KeyInt("REFINE_BLOBS_LIGAND", kinds[static_cast<int>(DensityBlob::Kind::Ligand)]));
+            Add(s, KeyInt("REFINE_BLOBS_WATER", kinds[static_cast<int>(DensityBlob::Kind::Water)]));
+            Add(s, KeyInt("REFINE_BLOBS_ION", kinds[static_cast<int>(DensityBlob::Kind::Ion)]));
+            Add(s, KeyInt("REFINE_BLOBS_MODEL_ERROR", kinds[static_cast<int>(DensityBlob::Kind::ModelError)]));
+            std::string top;
+            for (size_t i = 0; i < std::min<size_t>(rf.blobs.size(), 5); i++) {
+                const DensityBlob &b = rf.blobs[i];
+                top += fmt::format("\n  {:<11} peak {:5.1f} sigma, {:6.1f} A^3, at {:7.2f} {:7.2f} {:7.2f}, {:.1f} A from {}",
+                                   DensityBlobKindName(b.kind), b.peak_sigma, b.volume, b.centroid.x, b.centroid.y,
+                                   b.centroid.z, b.nearest_atom_distance, b.nearest_atom);
+            }
+            Add(s, Prose(fmt::format(
+                    " The placed model was refined against these data (--refine): each chain as a rigid body and\n"
+                    " one B per residue, then every atom restrained by the dictionary carried in rugnux, the\n"
+                    " X-ray weight chosen by R-free with the geometry held within bounds. REFINE_ keys describe\n"
+                    " the refined model, which is written as {}_refined_model.cif with its maps; R_WORK and\n"
+                    " R_FREE above remain the model as given. The refinement decided nothing about the\n"
+                    " processing. The strongest unexplained blobs of its difference map ({}_blobs.json):{}",
+                    output_prefix, output_prefix, top.empty() ? std::string(" none") : top)));
+        }
+        if (rf.ok) {
+            const RefinementScore &sc = rf.score_start;
+            Add(s, KeyReal("REFINE_SCORE_R_WORK", sc.r_work, "{:.4f}"));
+            Add(s, KeyReal("REFINE_SCORE_R_FREE", sc.r_free, "{:.4f}"));
+            Add(s, KeyInt("REFINE_SCORE_N_WORK", sc.n_work));
+            Add(s, KeyInt("REFINE_SCORE_N_FREE", sc.n_free));
+            Add(s, KeyReal("REFINE_SCORE_K_OVERALL", sc.k_overall, "{:.5g}", true));
+            Add(s, KeyReal("REFINE_SCORE_B_ISO", sc.b_iso, "{:.2f}", true));
+            Add(s, KeyText("REFINE_SCORE_B_ANISO", fmt::format("{:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f}", sc.b_aniso[0],
+                                                              sc.b_aniso[1], sc.b_aniso[2], sc.b_aniso[3], sc.b_aniso[4],
+                                                              sc.b_aniso[5]), true));
+            Add(s, KeyReal("REFINE_SCORE_K_SOL", sc.k_sol, "{:.4f}", true));
+            Add(s, KeyReal("REFINE_SCORE_B_SOL", sc.b_sol, "{:.2f}", true));
+            Add(s, KeyReal("REFINE_SCORE_MASK_SPACING", sc.mask_spacing, "{:.3f}", true));
+            Add(s, KeyInt("REFINE_SCORE_HYDROGENS", sc.hydrogens, true));
+            Add(s, KeyReal("REFINE_SCORE_NLL_WORK", sc.nll_work, "{:.1f}", true));
+            Add(s, KeyReal("REFINE_SCORE_NLL_FREE", sc.nll_free, "{:.1f}", true));
+            Add(s, KeyReal("REFINE_SCORE_FOM", sc.fom, "{:.4f}", true));
+            Add(s, KeyReal("REFINE_SCORE_RMSD_BONDS", rf.geometry_start.rmsd_bonds, "{:.4f}", true));
+            Add(s, KeyReal("REFINE_SCORE_RMSZ_BONDS", rf.geometry_start.rmsz_bonds, "{:.3f}", true));
+            Add(s, KeyReal("REFINE_SCORE_RMSD_ANGLES", rf.geometry_start.rmsd_angles, "{:.3f}", true));
+            Add(s, KeyReal("REFINE_SCORE_RMSZ_ANGLES", rf.geometry_start.rmsz_angles, "{:.3f}", true));
+            if (!sc.shells.empty()) {
+                ReportEntry e;
+                e.kind = ReportEntry::Kind::Table;
+                e.developer_only = true;
+                e.table.columns = {"D_MAX", "D_MIN", "N_WORK", "N_FREE", "R_WORK", "R_FREE"};
+                e.table.text_header = "   D_MAX   D_MIN  N_WORK  N_FREE  R_WORK  R_FREE";
+                for (const RefinementShell &sh : sc.shells) {
+                    e.table.text_rows.push_back(fmt::format("{:8.2f}{:8.2f}{:8d}{:8d}{:8.4f}{:8.4f}", sh.d_max, sh.d_min,
+                                                            sh.n_work, sh.n_free, sh.r_work, sh.r_free));
+                    std::vector<ReportValue> row;
+                    for (double v : {sh.d_max, sh.d_min, double(sh.n_work), double(sh.n_free), sh.r_work, sh.r_free}) {
+                        ReportValue rv;
+                        rv.type = ReportValue::Type::Real;
+                        rv.real = v;
+                        rv.text = fmt::format("{}", v);
+                        row.push_back(rv);
+                    }
+                    e.table.cells.push_back(std::move(row));
+                }
+                Add(s, std::move(e));
+            }
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------- the emitter
+// Populates the document. Nothing here writes to a stream: the order of the FILE is the renderer's
+// business, which is what lets the SUMMARY sit above the sections whose warnings it reports.
+ReportDocument BuildReportDocument(const std::string &output_prefix,
+                                   const std::string &input_file,
+                                   const DiffractionExperiment &experiment,
+                                   const ProcessResult &result,
+                                   const RunProvenance &provenance) {
+    ReportDocument doc;
+    const bool rotation = experiment.IsRotationIndexing();
+    const bool merged = result.has_merge_statistics;
+    doc.warnings = result.warnings;
+    // SUMMARY rows filled by the sections that measure them; empty leaves the row out.
+    std::string beam_centre_row, overload_row, overlap_row;
+
+    // ------------------------------------------------------------------ header block
+    doc.sections.push_back(HeaderSection(output_prefix, input_file, provenance));
+
     // ------------------------------------------------- 1. DATA SET AND GEOMETRY
     {
         ReportSection s;
@@ -2757,68 +2869,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
         }
         // --refine: the placed model refined (Refinement.h). Keys of their own: R_WORK / R_FREE above stay the
         // unrefined model's.
-        if (result.refinement.has_value()) {
-            const RefinementResult &rf = *result.refinement;
-            Add(s, Blank());
-            if (!rf.ok) {
-                Add(s, KeyEnum("REFINE", "NOT_PERFORMED", {"PERFORMED", "NOT_PERFORMED"}));
-                Add(s, KeyText("REFINE_REASON", rf.failure_reason));
-            } else {
-                Add(s, KeyEnum("REFINE", "PERFORMED", {"PERFORMED", "NOT_PERFORMED"}));
-                Add(s, KeyText("REFINE_DEVICE", rf.on_gpu ? "GPU" : "CPU", true));
-                Add(s, KeyReal("REFINE_SECONDS", rf.seconds, "{:.1f}", true));
-                Add(s, KeyReal("REFINE_D_MIN", rf.d_min, "{:.2f}"));
-                Add(s, KeyReal("REFINE_R_WORK_START", rf.r_work_start, "{:.4f}"));
-                Add(s, KeyReal("REFINE_R_FREE_START", rf.r_free_start, "{:.4f}"));
-                Add(s, KeyInt("REFINE_RIGID_BODIES", rf.rigid_bodies));
-                Add(s, KeyReal("REFINE_R_WORK_PLACEMENT", rf.r_work_placement, "{:.4f}"));
-                Add(s, KeyReal("REFINE_R_FREE_PLACEMENT", rf.r_free_placement, "{:.4f}"));
-                if (rf.restrained) {
-                    Add(s, KeyReal("REFINE_R_WORK", rf.r_work, "{:.4f}"));
-                    Add(s, KeyReal("REFINE_R_FREE", rf.r_free, "{:.4f}"));
-                    Add(s, KeyReal("REFINE_RMSD_BONDS", rf.geometry.rmsd_bonds, "{:.4f}"));
-                    Add(s, KeyReal("REFINE_RMSD_ANGLES", rf.geometry.rmsd_angles, "{:.2f}"));
-                    Add(s, KeyReal("REFINE_RMSD_BONDS_START", rf.geometry_start.rmsd_bonds, "{:.4f}", true));
-                    Add(s, KeyReal("REFINE_RMSD_ANGLES_START", rf.geometry_start.rmsd_angles, "{:.2f}", true));
-                    Add(s, KeyInt("REFINE_CLASHES", rf.geometry.clashes));
-                    Add(s, KeyReal("REFINE_XRAY_WEIGHT", rf.weight, "{:.4g}", true));
-                    Add(s, KeyInt("REFINE_ATOMS", rf.atoms));
-                    Add(s, KeyInt("REFINE_ATOMS_REFINED", rf.refined_atoms));
-                    Add(s, KeyInt("REFINE_RESIDUES_FIXED", rf.fixed_residues));
-                    if (!rf.fixed_names.empty()) {
-                        std::string names;
-                        for (const std::string &n : rf.fixed_names)
-                            names += (names.empty() ? "" : " ") + n;
-                        Add(s, KeyText("REFINE_RESIDUES_FIXED_NAMES", names));
-                    }
-                }
-                Add(s, KeyReal("REFINE_BULK_SOLVENT_K_SOL", rf.k_sol, "{:.3f}", true));
-                Add(s, KeyReal("REFINE_BULK_SOLVENT_B_SOL", rf.b_sol, "{:.1f}", true));
-                int kinds[5] = {0, 0, 0, 0, 0};
-                for (const DensityBlob &b : rf.blobs)
-                    kinds[static_cast<int>(b.kind)]++;
-                Add(s, KeyInt("REFINE_BLOBS", static_cast<int64_t>(rf.blobs.size())));
-                Add(s, KeyInt("REFINE_BLOBS_LIGAND", kinds[static_cast<int>(DensityBlob::Kind::Ligand)]));
-                Add(s, KeyInt("REFINE_BLOBS_WATER", kinds[static_cast<int>(DensityBlob::Kind::Water)]));
-                Add(s, KeyInt("REFINE_BLOBS_ION", kinds[static_cast<int>(DensityBlob::Kind::Ion)]));
-                Add(s, KeyInt("REFINE_BLOBS_MODEL_ERROR", kinds[static_cast<int>(DensityBlob::Kind::ModelError)]));
-                std::string top;
-                for (size_t i = 0; i < std::min<size_t>(rf.blobs.size(), 5); i++) {
-                    const DensityBlob &b = rf.blobs[i];
-                    top += fmt::format("\n  {:<11} peak {:5.1f} sigma, {:6.1f} A^3, at {:7.2f} {:7.2f} {:7.2f}, {:.1f} A from {}",
-                                       DensityBlobKindName(b.kind), b.peak_sigma, b.volume, b.centroid.x, b.centroid.y,
-                                       b.centroid.z, b.nearest_atom_distance, b.nearest_atom);
-                }
-                Add(s, Prose(fmt::format(
-                        " The placed model was refined against these data (--refine): each chain as a rigid body and\n"
-                        " one B per residue, then every atom restrained by the dictionary carried in rugnux, the\n"
-                        " X-ray weight chosen by R-free with the geometry held within bounds. REFINE_ keys describe\n"
-                        " the refined model, which is written as {}_refined_model.cif with its maps; R_WORK and\n"
-                        " R_FREE above remain the model as given. The refinement decided nothing about the\n"
-                        " processing. The strongest unexplained blobs of its difference map ({}_blobs.json):{}",
-                        output_prefix, output_prefix, top.empty() ? std::string(" none") : top)));
-            }
-        }
+        if (result.refinement.has_value())
+            AddRefinementKeys(s, *result.refinement, output_prefix);
         doc.sections.push_back(std::move(s));
     }
 
@@ -3142,6 +3194,42 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
 // -------------------------------------------------------------------------- the text rendering
 // Reads only from the document. Nothing here decides what the report SAYS; it decides only how the
 // text file looks, which is what makes a second rendering a second function rather than a fork.
+// --mode refine: the header, what the input was and what had to be made of it, and the refinement's keys.
+ReportDocument BuildRefineReportDocument(const std::string &output_prefix, const std::string &input_file,
+                                         const std::string &model_file, const RefineOnlyResult &r,
+                                         const RunProvenance &provenance) {
+    ReportDocument doc;
+    doc.sections.push_back(HeaderSection(output_prefix, input_file, provenance));
+    const RefinementResult &rf = r.refinement;
+    ReportSection s;
+    s.title = "1. REFINEMENT";
+    Add(s, KeyText("MODEL_FILE", model_file));
+    Add(s, KeyText("SPACE_GROUP", r.space_group));
+    Add(s, KeyReal("WAVELENGTH", r.wavelength_A, "{:.5f}"));
+    Add(s, KeyInt("REFINE_REFLECTIONS", r.reflections));
+    Add(s, KeyInt("REFINE_FREE_REFLECTIONS", r.free_reflections));
+    Add(s, KeyBool("REFINE_FRENCH_WILSON", r.french_wilson));
+    Add(s, KeyEnum("REFINE_FREE_FLAGS", r.free_flags_generated ? "GENERATED" : "FILE", {"FILE", "GENERATED"}));
+    if (r.d_min_limit > 0)
+        Add(s, KeyReal("REFINE_D_MIN_LIMIT", r.d_min_limit, "{:.3f}"));
+    AddRefinementKeys(s, rf, output_prefix);
+    doc.sections.push_back(std::move(s));
+    doc.verdict = rf.ok ? "OK" : "FAILED";
+    if (!rf.ok)
+        doc.verdict_text = "Refinement did not run: " + rf.failure_reason;
+    else if (r.score_only)
+        doc.verdict_text = fmt::format("Scored as given, nothing refined: R-work {:.4f}, R-free {:.4f}.",
+                                       rf.r_work_start, rf.r_free_start);
+    else
+        doc.verdict_text = fmt::format("Refined: R-work {:.4f}, R-free {:.4f} (as given {:.4f}, {:.4f}).",
+                                       rf.restrained ? rf.r_work : rf.r_work_placement,
+                                       rf.restrained ? rf.r_free : rf.r_free_placement, rf.r_work_start, rf.r_free_start);
+    Add(doc.sections.back(), KeyEnum("VERDICT", doc.verdict, {"OK", "WARNINGS", "UNUSABLE", "FAILED"}));
+    Add(doc.sections.back(), KeyText("VERDICT_TEXT", doc.verdict_text));
+    doc.summary_rows.push_back({"Refinement", doc.verdict_text});
+    return doc;
+}
+
 std::string RenderReportText(const ReportDocument &doc, bool developer) {
     std::ostringstream os;
     // Everything that is not a `KEY= value` line and not blank is written as a `#` comment, so a
diff --git a/rugnux/ResultReport.h b/rugnux/ResultReport.h
index 17705a68c..5078a22cf 100644
--- a/rugnux/ResultReport.h
+++ b/rugnux/ResultReport.h
@@ -9,6 +9,7 @@
 #include "../common/Logger.h"
 #include "ReportDocument.h"
 #include "Rugnux.h"
+#include "RugnuxRefine.h"
 
 // <prefix>_report.txt - what the run DETERMINED, written next to the .mtz/.cif/.hkl on every run
 // that has an output prefix. Modelled on XDS's CORRECT.LP: prose and tables a crystallographer
@@ -65,6 +66,11 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
                                    const ProcessResult &result,
                                    const RunProvenance &provenance = {});
 
+// --mode refine's report: the header, the input as read and the REFINE_ keys.
+ReportDocument BuildRefineReportDocument(const std::string &output_prefix, const std::string &input_file,
+                                         const std::string &model_file, const RefineOnlyResult &result,
+                                         const RunProvenance &provenance = {});
+
 // The text rendering of that document - XDS-style prose, tables and `KEY= value` lines. Reads only
 // from the document, so a second rendering is a second function of the same input.
 std::string RenderReportText(const ReportDocument &doc, bool developer);
diff --git a/rugnux/RugnuxOutputFiles.cpp b/rugnux/RugnuxOutputFiles.cpp
index a2ba9c836..f5650f5c8 100644
--- a/rugnux/RugnuxOutputFiles.cpp
+++ b/rugnux/RugnuxOutputFiles.cpp
@@ -188,3 +188,29 @@ ReportDocument WriteResultReport(const std::string &output_prefix,
     }
     return doc;
 }
+
+ReportDocument WriteRefineOutputs(const RefineOnlyResult &result, const std::string &prefix,
+                                  const std::string &input_file, const std::string &model_file, Logger &logger,
+                                  const RunProvenance &provenance) {
+    ReportDocument doc = BuildRefineReportDocument(prefix, input_file, model_file, result, provenance);
+    if (prefix.empty())
+        return doc;
+    const RefinementResult &rf = result.refinement;
+    if (rf.ok) {
+        const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(result.space_group);
+        const gemmi::UnitCell &c = rf.model->cell;
+        const UnitCell cell{static_cast<float>(c.a), static_cast<float>(c.b), static_cast<float>(c.c),
+                            static_cast<float>(c.alpha), static_cast<float>(c.beta), static_cast<float>(c.gamma)};
+        WritePlacedModel(*rf.model, cell, *sg, prefix + "_refined", logger);
+        WriteModelMaps(*rf.maps, prefix + "_refined", logger);
+        WriteDensityBlobs(rf.blobs, prefix + "_blobs.json", logger);
+    }
+    try {
+        std::ofstream file(prefix + "_report.txt");
+        file.exceptions(std::ios::failbit | std::ios::badbit);
+        file << RenderReportText(doc, provenance.developer);
+    } catch (const std::exception &e) {
+        logger.Warning("Could not write the results report {}_report.txt: {}", prefix, e.what());
+    }
+    return doc;
+}
diff --git a/rugnux/RugnuxOutputFiles.h b/rugnux/RugnuxOutputFiles.h
index 5821bd21c..ee3cab6af 100644
--- a/rugnux/RugnuxOutputFiles.h
+++ b/rugnux/RugnuxOutputFiles.h
@@ -42,6 +42,12 @@ void WriteMergedOutput(const MergedOutput &merged, const std::string &prefix, co
 // A model validation's map coefficients at `prefix`: _maps.mtz.
 void WriteModelMaps(const ModelMaps &maps, const std::string &prefix, Logger &logger);
 
+// --mode refine's files under `prefix`: <prefix>_refined_model.cif, <prefix>_refined_maps.mtz, <prefix>_blobs.json
+// and <prefix>_report.txt; returns the report's document. Writes nothing when the prefix is empty.
+ReportDocument WriteRefineOutputs(const RefineOnlyResult &result, const std::string &prefix,
+                                  const std::string &input_file, const std::string &model_file, Logger &logger,
+                                  const RunProvenance &provenance = {});
+
 // <prefix>_plot.txt (PerImagePlotText).
 void WritePerImagePlot(const std::string &prefix, const PerImagePlot &plot);
 
diff --git a/rugnux/RugnuxRefine.cpp b/rugnux/RugnuxRefine.cpp
new file mode 100644
index 000000000..e70c7611c
--- /dev/null
+++ b/rugnux/RugnuxRefine.cpp
@@ -0,0 +1,54 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "RugnuxRefine.h"
+
+#include <cmath>
+
+#include "../common/UnitCell.h"
+#include "../image_analysis/scale_merge/FrenchWilson.h"
+#include "../image_analysis/scale_merge/RfreeFlags.h"
+
+RefineOnlyResult RunRefineOnly(RugnuxReflections data, const gemmi::Structure &model, double d_min_limit,
+                               const RefinementSettings &settings, Logger &logger) {
+    RefineOnlyResult out;
+    out.d_min_limit = d_min_limit;
+    out.score_only = settings.score_only;
+    const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(data.space_group);
+    if (sg == nullptr) {
+        out.refinement.failure_reason = "unknown space group " + data.space_group;
+        return out;
+    }
+    out.space_group = sg->xhm();
+    out.wavelength_A = data.wavelength_A > 0 ? data.wavelength_A : 1.0;
+    std::vector<MergedReflection> &refl = data.reflections;
+    if (d_min_limit > 0)
+        std::erase_if(refl, [&](const MergedReflection &r) { return r.d < d_min_limit; });
+    if (!data.amplitudes) {
+        FrenchWilsonOptions fw;
+        fw.num_threads = static_cast<int>(settings.nthreads);
+        ApplyFrenchWilson(refl, *sg, fw);
+        out.french_wilson = true;
+    }
+    std::erase_if(refl, [](const MergedReflection &r) { return !std::isfinite(r.F); });
+    if (data.free_column.empty()) {
+        const UnitCell cell{static_cast<float>(data.cell.a), static_cast<float>(data.cell.b),
+                            static_cast<float>(data.cell.c), static_cast<float>(data.cell.alpha),
+                            static_cast<float>(data.cell.beta), static_cast<float>(data.cell.gamma)};
+        AssignRfreeFlags(refl, *sg, 0.05, 500, cell, settings.nthreads);
+        out.free_flags_generated = true;
+    }
+    out.reflections = static_cast<int>(refl.size());
+    for (const MergedReflection &r : refl)
+        out.free_reflections += r.rfree_flag;
+    logger.Info("Refinement input: {} ({}), {} {} in {}, {} free ({}){}", data.name, data.source, refl.size(),
+                data.amplitudes ? "amplitudes " + data.value_column
+                                : "amplitudes by French-Wilson from " + data.value_column,
+                out.space_group, out.free_reflections,
+                out.free_flags_generated ? "drawn here, the file has no free set" : "column " + data.free_column,
+                d_min_limit > 0 ? fmt::format(", to {:.2f} A", d_min_limit) : std::string());
+    if (const gemmi::SpaceGroup *msg = model.find_spacegroup(); msg != nullptr && msg->number != sg->number)
+        logger.Warning("The model names space group {}, the data {}: refined in the data's", msg->xhm(), sg->xhm());
+    out.refinement = RefineModel(model, refl, data.cell, *sg, out.wavelength_A, logger, settings);
+    return out;
+}
diff --git a/rugnux/RugnuxRefine.h b/rugnux/RugnuxRefine.h
new file mode 100644
index 000000000..78e2bf2af
--- /dev/null
+++ b/rugnux/RugnuxRefine.h
@@ -0,0 +1,31 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <string>
+
+#include "gemmi/model.hpp"
+
+#include "../common/Logger.h"
+#include "../image_analysis/structure_refinement/Refinement.h"
+#include "../writer/RugnuxReflectionFile.h"
+
+// rugnux --mode refine: a model refined against merged structure factors it is given, nothing processed. The
+// reflections and the model come in parsed (RugnuxReflectionFile.h, RugnuxModelFile.h); the refinement and what
+// had to be made of the input to run it go out. No file is read or written here.
+struct RefineOnlyResult {
+    RefinementResult refinement;
+    std::string space_group;         // the data's, which the model is refined in
+    double wavelength_A = 0;         // the file's, or 1 A where it states none
+    bool french_wilson = false;      // amplitudes made here from the file's intensities
+    bool free_flags_generated = false;   // the file has no free set: one was drawn here (AssignRfreeFlags)
+    int reflections = 0;             // amplitudes refined against, after the resolution limit
+    int free_reflections = 0;
+    double d_min_limit = 0;          // the limit asked for (0: none)
+    bool score_only = false;         // the model only scored (RefinementSettings::score_only)
+};
+
+// `data` as parsed; `d_min_limit`: leave out reflections finer than this (0: keep all).
+RefineOnlyResult RunRefineOnly(RugnuxReflections data, const gemmi::Structure &model, double d_min_limit,
+                               const RefinementSettings &settings, Logger &logger);
diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp
index 2c27ec893..902484d15 100644
--- a/rugnux/rugnux_cli.cpp
+++ b/rugnux/rugnux_cli.cpp
@@ -77,6 +77,8 @@
 #include "Rugnux.h"
 #include "RugnuxDefaults.h"
 #include "../image_analysis/structure_refinement/ModelValidation.h"
+#include "../writer/RugnuxReflectionFile.h"
+#include "RugnuxRefine.h"
 #include "../writer/WriteModel.h"
 #include "../writer/RugnuxModelFile.h"
 #include "ResultReport.h"
@@ -95,7 +97,7 @@ constexpr double RELATIVE_B_DEFAULT_DEG = 10.0;
 
 // What the run does (--mode). Each of these was once a flag of its own; one option keeps them
 // mutually exclusive by construction.
-enum class RugnuxMode { MX, Azint, Scale, Calibration };
+enum class RugnuxMode { MX, Azint, Scale, Calibration, Refine };
 
 // Azimuthal sectors used by --mode calibration --calibration rings when --azim-phi-bins was not given.
 // The default of 1 is a plain radial profile, in which a ring has been averaged over every direction so
@@ -151,6 +153,13 @@ void print_usage() {
     std::cout << "        azint         Only azimuthal integration (no spot finding/indexing); writes <prefix>_process.h5" << std::endl;
     std::cout << "        scale         Only re-scale/merge the already-integrated reflections in <input> (no re-integration)" << std::endl;
     std::cout << "        calibration   Determine the detector geometry from powder rings; writes <prefix>.poni for pyFAI and <prefix>.json, whose dataset_settings member is a jfjoch_broker dataset_settings body as it stands" << std::endl;
+    std::cout << "        refine        Refine the --model against merged structure factors: <input> is an MTZ or a structure-factor mmCIF (gzipped or not), nothing is processed. As --refine does after a merge, and written alike: <prefix>_refined_model.cif, <prefix>_refined_maps.mtz, <prefix>_blobs.json and REFINE_ keys in <prefix>_report.txt" << std::endl;
+    std::cout << std::endl;
+
+    std::cout << "  Refinement (--mode refine)" << std::endl;
+    std::cout << "   The columns are found, not named: amplitudes with their sigmas (F/SIGF, FP/SIGFP, ...) where the file has them, else intensities (IMEAN/SIGIMEAN, I/SIGI, ...) made into amplitudes by French-Wilson, else the mean of the anomalous pairs; the free set from a FreeR_flag-like column (0 the test set, or the minority value where 0 is the majority, as phenix writes it), drawn here where the file has none. The model is taken as placed in the data's cell, and refined in the data's space group" << std::endl;
+    std::cout << "   --refine-score-only              Score the model as given and refine nothing: R-work/R-free, the scale, the bulk solvent, the likelihood and the geometry, each term under REFINE_SCORE_ (--developer prints them all)" << std::endl;
+    std::cout << "   --refine-d-min <num>             Leave out reflections finer than this (A)" << std::endl;
     std::cout << std::endl;
 
     std::cout << "  Calibration (--mode calibration)" << std::endl;
@@ -388,7 +397,9 @@ enum {
     OPT_DETECTOR_QUARTER_TURNS,
     OPT_FFT_MIN_UNIT_CELL,
     OPT_POLARIZATION,
-    OPT_NO_REFINE_TILT
+    OPT_NO_REFINE_TILT,
+    OPT_REFINE_SCORE_ONLY,
+    OPT_REFINE_D_MIN
 };
 
 static option long_options[] = {
@@ -408,6 +419,8 @@ static option long_options[] = {
     {"crystfel-geometry", required_argument, nullptr, OPT_CRYSTFEL_GEOMETRY},
     {"mr", required_argument, nullptr, OPT_MR},
     {"refine", no_argument, nullptr, OPT_REFINE},
+    {"refine-score-only", no_argument, nullptr, OPT_REFINE_SCORE_ONLY},
+    {"refine-d-min", required_argument, nullptr, OPT_REFINE_D_MIN},
     {"export-unmerged", no_argument, nullptr, OPT_EXPORT_UNMERGED},
     {"no-export-unmerged", no_argument, nullptr, OPT_NO_EXPORT_UNMERGED},
     {"export-unmerged-partials", no_argument, nullptr, OPT_EXPORT_UNMERGED_PARTIALS},
@@ -1041,6 +1054,8 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
     std::string model_pdb;           // --model: PDB to validate merged intensities against (R-free + maps)
     std::string mr_model;            // --mr: --model, placed by molecular replacement straight away
     bool refine_model = false;       // --refine
+    bool refine_score_only = false;  // --refine-score-only
+    double refine_d_min = 0;         // --refine-d-min
     bool export_unmerged = true;            // on by default; --no-export-unmerged turns it off
     bool export_unmerged_partials = false;  // --export-unmerged-partials: the same, one row per image
     std::string dump_observations;   // diagnostic: dump unmerged -P rot3d fulls to this path
@@ -1286,6 +1301,12 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
         case OPT_REFINE:
             refine_model = true;
             break;
+        case OPT_REFINE_SCORE_ONLY:
+            refine_score_only = true;
+            break;
+        case OPT_REFINE_D_MIN:
+            refine_d_min = parse_number_arg<double>(optarg, "--refine-d-min", logger, 0.1, 100.0);
+            break;
         case OPT_EXPORT_UNMERGED:
             export_unmerged = true;
             break;
@@ -1382,8 +1403,9 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
             else if (m == "azint") mode = RugnuxMode::Azint;
             else if (m == "scale") mode = RugnuxMode::Scale;
             else if (m == "calibration") mode = RugnuxMode::Calibration;
+            else if (m == "refine") mode = RugnuxMode::Refine;
             else {
-                logger.Error("Invalid --mode: {} (expected mx|azint|scale|calibration)", m);
+                logger.Error("Invalid --mode: {} (expected mx|azint|scale|calibration|refine)", m);
                 return 1;
             }
             break;
@@ -1753,6 +1775,10 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
         logger.Error("--model and --mr both name the model; give one");
         exit(EXIT_FAILURE);
     }
+    if (mode == RugnuxMode::Refine && model_pdb.empty()) {
+        logger.Error("--mode refine refines a model: give it with --model");
+        exit(EXIT_FAILURE);
+    }
     if (refine_model && mr_model.empty() && model_pdb.empty()) {
         logger.Error("--refine refines a model: give it with --model or --mr");
         exit(EXIT_FAILURE);
@@ -1818,6 +1844,39 @@ static int RunRugnux(int argc, char **argv, bool no_gpu) {
     } else
         logger.Info("Threads: {}{}", nthreads, nthreads_auto ? " (the CPUs this run may use)" : " (-N)");
 
+    // --mode refine: the model against the structure factors given, nothing processed.
+    if (mode == RugnuxMode::Refine) {
+        console->Field("Input", input_file);
+        console->Field("Output", log_file.empty() ? output_prefix + "_*"
+                                                  : fmt::format("{}_*   (full log in {})", output_prefix, log_file));
+        console->Rule('-');
+        console->OnStep("Reading the reflections");
+        RugnuxReflections data;
+        try {
+            data = ReadRugnuxReflectionFile(input_file);
+        } catch (const std::exception &e) {
+            logger.Error("Cannot read the reflections {}: {}", input_file, e.what());
+            exit(EXIT_FAILURE);
+        }
+        console->OnStep(refine_score_only ? "Scoring the model" : "Refining the model");
+        RefinementSettings settings;
+        settings.nthreads = static_cast<size_t>(nthreads);
+        settings.gpu = !no_gpu;
+        settings.score_only = refine_score_only;
+        const RefineOnlyResult refined = RunRefineOnly(std::move(data), model->structure, refine_d_min, settings,
+                                                       logger);
+        if (!refined.refinement.ok)
+            logger.Error("Refinement did not run: {}", refined.refinement.failure_reason);
+        console->OnStep("Writing files");
+        provenance.wall_time_s = std::chrono::duration<double>(
+                std::chrono::steady_clock::now() - invocation_start).count();
+        const ReportDocument doc = WriteRefineOutputs(refined, output_prefix, input_file, model_pdb, logger,
+                                                      provenance);
+        console->Finish();
+        PrintSummary(*console, doc, output_prefix, files_since, provenance.wall_time_s);
+        return refined.refinement.ok ? 0 : 1;
+    }
+
     std::string input_list;
     for (const auto &f: input_files)
         input_list += (input_list.empty() ? "" : ", ") + f;
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 1b39e1607..f94e6495b 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -184,6 +184,7 @@ ADD_EXECUTABLE(jfjoch_test
         ModelMaskGPUTest.cpp
         ModelDensityGPUTest.cpp
         RefinementTest.cpp
+        RefineOnlyTest.cpp
         SearchSpaceGroupTest.cpp
         SearchSpaceGroupTwinTest.cpp
         SyntheticMergedReflections.h
diff --git a/tests/RefineOnlyTest.cpp b/tests/RefineOnlyTest.cpp
new file mode 100644
index 000000000..670971f4a
--- /dev/null
+++ b/tests/RefineOnlyTest.cpp
@@ -0,0 +1,207 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+
+#include <cmath>
+#include <random>
+#include <string>
+#include <vector>
+
+#include <gemmi/dencalc.hpp>
+#include <gemmi/fourier.hpp>
+#include <gemmi/it92.hpp>
+#include <gemmi/mtz.hpp>
+
+#include "../common/Logger.h"
+#include "../rugnux/ResultReport.h"
+#include "../rugnux/RugnuxRefine.h"
+#include "../writer/RugnuxReflectionFile.h"
+
+namespace {
+
+// A small synthetic orthorhombic cell (no real sample) and its reflections h, k, l = 1..6.
+const gemmi::UnitCell test_cell(31.0, 37.0, 43.0, 90.0, 90.0, 90.0);
+
+// An MTZ of `columns` (label, type), every row the hkl followed by values from `value(h, k, l, column)`.
+gemmi::Mtz MakeMtz(const std::vector<std::pair<std::string, char>> &columns,
+                   float (*value)(int h, int k, int l, size_t column)) {
+    gemmi::Mtz mtz(true);
+    mtz.spacegroup = gemmi::find_spacegroup_by_name("P 21 21 21");
+    mtz.set_cell_for_all(test_cell);
+    mtz.add_dataset("test").wavelength = 0.9;
+    for (const auto &[label, type] : columns)
+        mtz.add_column(label, type, -1, -1, false);
+    for (int h = 1; h <= 6; h++)
+        for (int k = 1; k <= 6; k++)
+            for (int l = 1; l <= 6; l++) {
+                mtz.data.insert(mtz.data.end(), {float(h), float(k), float(l)});
+                for (size_t c = 0; c < columns.size(); c++)
+                    mtz.data.push_back(value(h, k, l, c));
+                mtz.nreflections++;
+            }
+    return mtz;
+}
+
+bool InTestSet(int h, int k, int l) { return (h * 7 + k * 13 + l * 17) % 10 == 0; }
+
+}  // namespace
+
+TEST_CASE("RugnuxReflectionsFromMtz takes observed amplitudes over map coefficients", "[refine_only]") {
+    // FWT comes first and is an amplitude with a phase after it; FP/SIGFP is what is observed.
+    const gemmi::Mtz mtz = MakeMtz({{"FWT", 'F'}, {"PHWT", 'P'}, {"FP", 'F'}, {"SIGFP", 'Q'}, {"FreeR_flag", 'I'}},
+                                   [](int h, int k, int l, size_t c) -> float {
+                                       switch (c) {
+                                           case 0: return 999.0f;
+                                           case 1: return 10.0f;
+                                           case 2: return 100.0f + h;
+                                           case 3: return 2.0f;
+                                           default: return InTestSet(h, k, l) ? 0.0f : float(1 + (h + k + l) % 19);
+                                       }
+                                   });
+    const RugnuxReflections r = RugnuxReflectionsFromMtz(mtz, "synthetic.mtz");
+    CHECK(r.amplitudes);
+    CHECK(r.value_column == "FP");
+    CHECK(r.sigma_column == "SIGFP");
+    CHECK(r.free_column == "FreeR_flag");
+    CHECK(r.wavelength_A == Catch::Approx(0.9));
+    REQUIRE(r.reflections.size() == 216);
+    int n_free = 0;
+    for (const MergedReflection &m : r.reflections) {
+        CHECK(m.F == Catch::Approx(100.0 + m.h));
+        CHECK(m.rfree_flag == InTestSet(m.h, m.k, m.l));
+        n_free += m.rfree_flag;
+    }
+    CHECK(r.n_free == n_free);
+}
+
+TEST_CASE("RugnuxReflectionsFromMtz falls back to intensities and reads a 1-marks-free flag", "[refine_only]") {
+    const gemmi::Mtz mtz = MakeMtz({{"IMEAN", 'J'}, {"SIGIMEAN", 'Q'}, {"R-free-flags", 'I'}},
+                                   [](int h, int k, int l, size_t c) -> float {
+                                       if (c == 0) return 50.0f * h;
+                                       if (c == 1) return 5.0f;
+                                       return InTestSet(h, k, l) ? 1.0f : 0.0f;
+                                   });
+    const RugnuxReflections r = RugnuxReflectionsFromMtz(mtz, "synthetic.mtz");
+    CHECK_FALSE(r.amplitudes);
+    CHECK(r.value_column == "IMEAN");
+    CHECK(r.free_column == "R-free-flags");
+    for (const MergedReflection &m : r.reflections) {
+        CHECK(m.I == Catch::Approx(50.0 * m.h));
+        CHECK(std::isnan(m.F));
+        CHECK(m.rfree_flag == InTestSet(m.h, m.k, m.l));
+    }
+}
+
+TEST_CASE("RugnuxReflectionsFromMtz without free flags, and without any value column", "[refine_only]") {
+    const gemmi::Mtz plain = MakeMtz({{"F", 'F'}, {"SIGF", 'Q'}},
+                                     [](int h, int, int, size_t c) -> float { return c == 0 ? 10.0f * h : 1.0f; });
+    const RugnuxReflections r = RugnuxReflectionsFromMtz(plain, "synthetic.mtz");
+    CHECK(r.value_column == "F");
+    CHECK(r.free_column.empty());
+    CHECK(r.n_free == 0);
+
+    const gemmi::Mtz calculated = MakeMtz({{"FC", 'F'}, {"PHIC", 'P'}},
+                                          [](int, int, int, size_t) -> float { return 1.0f; });
+    CHECK_THROWS(RugnuxReflectionsFromMtz(calculated, "synthetic.mtz"));
+}
+
+TEST_CASE("ParseRugnuxReflections reads an MTZ from its bytes", "[refine_only]") {
+    const gemmi::Mtz mtz = MakeMtz({{"F", 'F'}, {"SIGF", 'Q'}, {"FreeR_flag", 'I'}},
+                                   [](int h, int k, int l, size_t c) -> float {
+                                       if (c == 0) return 10.0f * h;
+                                       if (c == 1) return 1.0f;
+                                       return InTestSet(h, k, l) ? 0.0f : 1.0f;
+                                   });
+    std::string bytes;
+    mtz.write_to_string(bytes);
+    const RugnuxReflections r = ParseRugnuxReflections(bytes, "synthetic.mtz");
+    CHECK(r.source == "MTZ");
+    CHECK(r.space_group == "P 21 21 21");
+    CHECK(r.reflections.size() == 216);
+    CHECK(r.free_column == "FreeR_flag");
+}
+
+TEST_CASE("RunRefineOnly scores a model against intensities, with amplitudes and free flags made for it",
+          "[refine_only]") {
+    // Atoms scattered over the cell; the data are their own |F|^2, so the model as given scores well.
+    const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 21 21 21");
+    gemmi::Structure st;
+    st.cell = test_cell;
+    st.spacegroup_hm = sg->xhm();
+    st.setup_cell_images();
+    gemmi::Model model(1);
+    gemmi::Chain chain("A");
+    std::mt19937 rng(7);
+    std::uniform_real_distribution<double> frac(0.0, 1.0);
+    for (int i = 0; i < 120; i++) {
+        gemmi::Residue res;
+        res.name = "HOH";
+        res.seqid = gemmi::SeqId(i + 1, ' ');
+        res.het_flag = 'H';
+        gemmi::Atom a;
+        a.name = "O";
+        a.element = gemmi::El::O;
+        a.pos = test_cell.orthogonalize(gemmi::Fractional(frac(rng), frac(rng), frac(rng)));
+        a.b_iso = 15.0f;
+        a.occ = 1.0f;
+        res.atoms.push_back(a);
+        chain.residues.push_back(res);
+    }
+    model.chains.push_back(chain);
+    st.models.push_back(model);
+
+    gemmi::DensityCalculator<gemmi::IT92<float>, float> dc;
+    dc.d_min = 2.0;
+    dc.grid.unit_cell = test_cell;
+    dc.grid.spacegroup = sg;
+    dc.set_refmac_compatible_blur(st.models[0]);
+    dc.put_model_density_on_grid(st.models[0]);
+    const gemmi::AsuData<std::complex<float>> f =
+            gemmi::transform_map_to_f_phi(dc.grid, true).prepare_asu_data(2.0, dc.blur);
+
+    RugnuxReflections data;
+    data.name = "synthetic.mtz";
+    data.source = "MTZ";
+    data.cell = test_cell;
+    data.space_group = sg->xhm();
+    data.amplitudes = false;
+    for (const auto &hv : f.v) {
+        if (hv.hkl == gemmi::Miller{{0, 0, 0}})
+            continue;
+        const double fc = std::abs(hv.value);
+        MergedReflection m;
+        m.h = hv.hkl[0];
+        m.k = hv.hkl[1];
+        m.l = hv.hkl[2];
+        m.d = static_cast<float>(test_cell.calculate_d(hv.hkl));
+        m.I = static_cast<float>(fc * fc);
+        m.sigma = static_cast<float>(0.02 * fc * fc + 1.0);
+        data.reflections.push_back(m);
+    }
+
+    RefinementSettings settings;
+    settings.score_only = true;
+    settings.gpu = false;
+    settings.nthreads = 2;
+    Logger logger("RefineOnlyTest");
+    const RefineOnlyResult r = RunRefineOnly(data, st, 2.5, settings, logger);
+    REQUIRE(r.refinement.ok);
+    CHECK(r.french_wilson);
+    CHECK(r.free_flags_generated);
+    CHECK(r.free_reflections > 0);
+    CHECK(r.d_min_limit == 2.5);
+    CHECK(r.refinement.d_min >= 2.5);
+    CHECK(r.refinement.r_work_start < 0.15);
+    CHECK(r.refinement.score_start.r_work == Catch::Approx(r.refinement.r_work_start));
+    CHECK(r.refinement.score_start.shells.size() == 20);
+    CHECK_FALSE(r.refinement.restrained);
+
+    const ReportDocument doc = BuildRefineReportDocument("out", "synthetic.mtz", "model.cif", r);
+    CHECK(doc.verdict == "OK");
+    const std::string text = RenderReportText(doc, true);
+    CHECK(text.find("REFINE= PERFORMED") != std::string::npos);
+    CHECK(text.find("REFINE_SCORE_R_WORK=") != std::string::npos);
+    CHECK(text.find("REFINE_FREE_FLAGS= GENERATED") != std::string::npos);
+    CHECK(text.find("REFINE_FRENCH_WILSON= TRUE") != std::string::npos);
+}
diff --git a/tools/battery/refine_parity.py b/tools/battery/refine_parity.py
new file mode 100644
index 000000000..d7efed0fb
--- /dev/null
+++ b/tools/battery/refine_parity.py
@@ -0,0 +1,176 @@
+#!/usr/bin/env python3
+"""refine_parity.py -- rugnux's refinement score against REFMAC's, term by term, on identical inputs.
+
+For one model and one merged MTZ (F, SIGF, FreeR_flag) it scores the model as given, refining nothing, twice:
+
+  * REFMAC 5.8, restrained setup at 0 cycles (`ncyc 0`, `refi bref iso`, riding hydrogens restored as REFMAC
+    does by default), read from its log: R-work / R-free overall and in its 20 resolution bins, the overall
+    scale and B, the anisotropic B, the bulk solvent's k_sol / B_sol, the FOM, -LL / -LLfree and the geometry
+    (rmsd and rmsZ of bonds and angles);
+  * rugnux --mode refine --refine-score-only on the same files, read from its report (REFINE_SCORE_* keys).
+
+The two are printed side by side with their differences, and one JSON line goes to stdout for scripts.
+REFMAC's anisotropic B is printed in the convention exp(-sum B_ij h_i h_j a*_i a*_j), rugnux's in
+exp(-s^T B s / 4); the REFMAC numbers are multiplied by 4 here so that the two columns compare.
+
+Per reflection, REFMAC's model amplitude on the input's scale is recovered from its output MTZ: the output F
+is the input F times REFMAC's overall scale, so FC_ALL_LS * F_in / F_out is the model amplitude R is computed
+with. With --per-reflection the R of each data subset (resolution bins) is recomputed from both programs'
+model amplitudes on the same reflections, rugnux's being FC of its _refined_maps.mtz.
+
+  refine_parity.py --rugnux build/rugnux/rugnux --workdir DIR [--d-min D] [--refmac-extra 'make hydr no'] model data.mtz
+"""
+import argparse
+import json
+import os
+import re
+import subprocess
+
+import gemmi
+import numpy as np
+
+CCP4_SETUP = os.environ.get("CCP4_SETUP", "/opt/xtal/ccp4-9/bin/ccp4.setup-sh")
+
+
+def ccp4_env():
+    out = subprocess.run(["bash", "-c", f"source {CCP4_SETUP} >/dev/null 2>&1; env -0"], capture_output=True).stdout
+    return dict(kv.split("=", 1) for kv in out.decode().split("\0") if "=" in kv)
+
+
+def run_refmac(workdir, model, mtz, d_min, extra):
+    keywords = "labin FP=F SIGFP=SIGF FREE=FreeR_flag\nncyc 0\nweight auto\nrefi bref iso\n" \
+               "make newligand continue\nmake check none\n"
+    if d_min:
+        keywords += f"reso 999 {d_min:.3f}\n"
+    keywords += "".join(line + "\n" for line in extra) + "end\n"
+    log = os.path.join(workdir, "refmac.log")
+    with open(log, "w") as fh:
+        subprocess.run(["nice", "refmac5", "hklin", mtz, "xyzin", model,
+                        "hklout", os.path.join(workdir, "refmac.mtz"), "xyzout", os.path.join(workdir, "refmac.pdb")],
+                       input=keywords.encode(), stdout=fh, stderr=subprocess.STDOUT, cwd=workdir, env=ccp4_env())
+    return parse_refmac(open(log, errors="replace").read())
+
+
+def parse_refmac(text):
+    num = lambda pat: float(re.search(pat, text, re.M).group(1))
+    out = {"r_work": num(r"^Overall R factor\s+=\s+([\d.]+)"), "r_free": num(r"^Free R factor\s+=\s+([\d.]+)")}
+    m = re.search(r"Overall\s+: scale =\s+([-\d.]+), B\s+=\s+([-\d.]+)", text)
+    out["k_overall"], out["b_iso"] = float(m.group(1)), float(m.group(2))
+    m = re.search(r"Partial structure\s+1: scale =\s+([-\d.]+), B\s+=\s+([-\d.]+)", text)
+    if m:
+        out["k_sol"], out["b_sol"] = float(m.group(1)), float(m.group(2))
+    m = re.search(r"B11 =\s*([-\d.]+) B22 =\s*([-\d.]+) B33 =\s*([-\d.]+) B12 =\s*([-\d.]+) B13 =\s*([-\d.]+) "
+                  r"B23 =\s*([-\d.]+)", text)
+    out["b_aniso"] = [4 * float(m.group(i)) for i in range(1, 7)]
+    m = re.search(r"^\s+0\s+([\d.]+)\s+([\d.]+)\s+([\d.]+)\s+([\d.]+)\s+([\d.]+)\s+([\d.]+)\s+([\d.]+)\s+([\d.]+)"
+                  r"\s+([\d.]+)", text, re.M)
+    if m:
+        (out["fom"], out["nll_work"], out["nll_free"], out["rmsd_bonds"], out["rmsz_bonds"], out["rmsd_angles"],
+         out["rmsz_angles"]) = [float(m.group(i)) for i in range(3, 10)]
+    table = text.split("M(4SSQ/LL) NR_used")[1].split("$$")[2]
+    out["shells"] = []
+    for line in table.strip().splitlines():
+        v = line.split()
+        out["shells"].append({"inv_d2": float(v[0]), "n_work": int(v[1]), "r_work": float(v[5]),
+                              "n_free": int(v[7]), "r_free": float(v[10])})
+    return out
+
+
+def run_rugnux(rugnux, workdir, model, mtz, d_min):
+    prefix = os.path.join(workdir, "rugnux")
+    cmd = ["nice", rugnux, "--mode", "refine", "--refine-score-only", "--developer", "--no-gpu", "-N", "8", "-o", prefix,
+           "--model", model]
+    if d_min:
+        cmd += ["--refine-d-min", f"{d_min:.3f}"]
+    subprocess.run(cmd + [mtz], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
+    text = open(prefix + "_report.txt").read()
+    key = lambda k: re.search(rf"^{k}=\s*(.+)$", text, re.M).group(1).strip()
+    out = {k: float(key("REFINE_SCORE_" + k.upper())) for k in
+           ("r_work", "r_free", "k_overall", "b_iso", "k_sol", "b_sol", "nll_work", "nll_free", "fom",
+            "rmsd_bonds", "rmsz_bonds", "rmsd_angles", "rmsz_angles")}
+    out["b_aniso"] = [float(x) for x in key("REFINE_SCORE_B_ANISO").split()]
+    out["mask_spacing"] = float(key("REFINE_SCORE_MASK_SPACING"))
+    return out
+
+
+def per_reflection(workdir, data_mtz):
+    """R per resolution bin of both programs' model amplitudes on the reflections both score."""
+    din = gemmi.read_mtz_file(data_mtz)
+    ref = gemmi.read_mtz_file(os.path.join(workdir, "refmac.mtz"))
+    rug = gemmi.read_mtz_file(os.path.join(workdir, "rugnux_refined_maps.mtz"))
+
+    asu, gops = gemmi.ReciprocalAsu(din.spacegroup), din.spacegroup.operations()
+
+    def table(mtz, labels):
+        a = np.array(mtz, copy=False)
+        cols = {c.label: i for i, c in enumerate(mtz.columns)}
+        out = {}
+        for row in a:
+            h = tuple(asu.to_asu([int(x) for x in row[:3]], gops)[0])
+            out[(h[0] + 512) * 1048576 + (h[1] + 512) * 1024 + h[2] + 512] = row[[cols[l] for l in labels]]
+        return out
+    t_in = table(din, ["F", "FreeR_flag"])
+    t_ref = table(ref, ["F", "FC_ALL_LS"])
+    t_rug = table(rug, ["FC"])
+    rows = []
+    for k, (fo_ref, fc_ref) in t_ref.items():
+        if k in t_in and k in t_rug and fo_ref > 0 and np.isfinite(fc_ref):
+            fo, flag = t_in[k]
+            rows.append((k, fo, flag == 0, fc_ref * fo / fo_ref, t_rug[k][0], 0.0))
+    rows = np.array(rows)
+    if len(rows):
+        h = rows[:, 0].astype(np.int64)
+        hkl = np.stack([h // 1048576 - 512, (h // 1024) % 1024 - 512, h % 1024 - 512], axis=1)
+        rows[:, 5] = [1 / din.cell.calculate_d(list(map(int, x))) ** 2 for x in hkl]
+    return rows
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("model")
+    ap.add_argument("data")
+    ap.add_argument("--rugnux", required=True)
+    ap.add_argument("--workdir", required=True)
+    ap.add_argument("--d-min", type=float, default=None)
+    ap.add_argument("--refmac-extra", action="append", default=[])
+    ap.add_argument("--no-refmac", action="store_true", help="reuse refmac.log in the workdir")
+    args = ap.parse_args()
+    os.makedirs(args.workdir, exist_ok=True)
+    model, data = os.path.abspath(args.model), os.path.abspath(args.data)
+    if args.no_refmac:
+        ref = parse_refmac(open(os.path.join(args.workdir, "refmac.log"), errors="replace").read())
+    else:
+        ref = run_refmac(args.workdir, model, data, args.d_min, args.refmac_extra)
+    rug = run_rugnux(args.rugnux, args.workdir, model, data, args.d_min)
+    print(f"{'term':<14} {'REFMAC':>10} {'rugnux':>10} {'diff':>10}")
+    for k in ("r_work", "r_free", "k_overall", "b_iso", "k_sol", "b_sol", "fom", "nll_work", "nll_free",
+              "rmsd_bonds", "rmsz_bonds", "rmsd_angles", "rmsz_angles"):
+        a, b = ref.get(k), rug.get(k)
+        if a is None or b is None:
+            continue
+        print(f"{k:<14} {a:10.4f} {b:10.4f} {b - a:+10.4f}")
+    print("b_aniso REFMAC(x4) " + " ".join(f"{x:7.2f}" for x in ref["b_aniso"]))
+    print("b_aniso rugnux     " + " ".join(f"{x:7.2f}" for x in rug["b_aniso"]))
+    rows = per_reflection(args.workdir, data)
+    if len(rows):
+        fo, free, fr, fg = rows[:, 1], rows[:, 2] > 0, rows[:, 3], rows[:, 4]
+        print(f"same reflections ({len(rows)}): R-work REFMAC {np.abs(fo - fr)[~free].sum() / fo[~free].sum():.4f} "
+              f"rugnux {np.abs(fo - fg)[~free].sum() / fo[~free].sum():.4f}; R-free REFMAC "
+              f"{np.abs(fo - fr)[free].sum() / fo[free].sum():.4f} rugnux {np.abs(fo - fg)[free].sum() / fo[free].sum():.4f}")
+        # by resolution, on the work set: R of each and the mean model amplitude of rugnux over REFMAC's
+        s2 = rows[:, 5]
+        edges = np.linspace(s2.min(), s2.max() + 1e-12, 11)
+        print(f"{'d_max':>7} {'d_min':>7} {'n':>6} {'R_REFMAC':>9} {'R_rugnux':>9} {'<Fr/Fref>':>10}")
+        for lo, hi in zip(edges[:-1], edges[1:]):
+            sel = (s2 >= lo) & (s2 < hi) & ~free
+            if sel.sum() == 0:
+                continue
+            print(f"{1 / np.sqrt(lo):7.2f} {1 / np.sqrt(hi):7.2f} {sel.sum():6d} "
+                  f"{np.abs(fo - fr)[sel].sum() / fo[sel].sum():9.4f} {np.abs(fo - fg)[sel].sum() / fo[sel].sum():9.4f} "
+                  f"{fg[sel].sum() / fr[sel].sum():10.4f}")
+    ref_short = {k: v for k, v in ref.items() if k != "shells"}
+    print("JSON " + json.dumps({"refmac": ref_short, "rugnux": rug, "refmac_shells": ref["shells"]}))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/writer/CMakeLists.txt b/writer/CMakeLists.txt
index 003a5eaab..7566823e3 100644
--- a/writer/CMakeLists.txt
+++ b/writer/CMakeLists.txt
@@ -38,10 +38,11 @@ TARGET_LINK_LIBRARIES(JFJochModelWriter JFJochLogger gemmi)
 
 # rugnux's merged MTZ, map-coefficient MTZ and placed model, to and from the structured END-message
 # blocks, so that jfjoch_writer writes them with gemmi alone and no analysis library; and the atomic model
-# a run is given, parsed for it.
+# a run is given and the merged reflections --mode refine is given, parsed for it.
 ADD_LIBRARY(JFJochRugnuxFiles STATIC WriteRugnuxFiles.cpp WriteRugnuxFiles.h
         RugnuxResultBlocks.cpp RugnuxResultBlocks.h
-        RugnuxModelFile.cpp RugnuxModelFile.h)
+        RugnuxModelFile.cpp RugnuxModelFile.h
+        RugnuxReflectionFile.cpp RugnuxReflectionFile.h)
 TARGET_LINK_LIBRARIES(JFJochRugnuxFiles JFJochModelWriter JFJochLogger gemmi ZLIB::ZLIB)
 
 TARGET_LINK_LIBRARIES(JFJochWriter JFJochPreview JFJochLogger JFJochHDF5Wrappers CBORStream2FrameSerialize
diff --git a/writer/RugnuxModelFile.cpp b/writer/RugnuxModelFile.cpp
index 287872ff4..e43ef1b45 100644
--- a/writer/RugnuxModelFile.cpp
+++ b/writer/RugnuxModelFile.cpp
@@ -53,10 +53,14 @@ std::string ReadFile(const std::string &path) {
 
 }  // namespace
 
+std::string GunzipIfCompressed(const std::string &bytes) {
+    return bytes.size() >= 2 && bytes[0] == '\x1f' && bytes[1] == '\x8b' ? Gunzip(bytes) : bytes;
+}
+
 RugnuxModel ParseRugnuxModel(const std::string &bytes, const std::string &name, const std::string &pae_json) {
     RugnuxModel m;
     m.name = name;
-    std::string text = bytes.size() >= 2 && bytes[0] == '\x1f' && bytes[1] == '\x8b' ? Gunzip(bytes) : bytes;
+    std::string text = GunzipIfCompressed(bytes);
     gemmi::cif::Document doc;
     m.structure = gemmi::read_structure_from_memory(text.data(), text.size(), name, gemmi::CoorFormat::Detect, &doc);
     if (gemmi::count_atom_sites(m.structure) == 0)
diff --git a/writer/RugnuxModelFile.h b/writer/RugnuxModelFile.h
index 5f242e6b1..57f9c304a 100644
--- a/writer/RugnuxModelFile.h
+++ b/writer/RugnuxModelFile.h
@@ -14,6 +14,9 @@
 // atomic model.
 RugnuxModel ParseRugnuxModel(const std::string &bytes, const std::string &name, const std::string &pae_json = {});
 
+// `bytes` inflated where they are a gzip stream (by the magic number), as they are otherwise.
+std::string GunzipIfCompressed(const std::string &bytes);
+
 // The model file at `path`, with the predicted aligned error the AlphaFold DB puts beside it where it is there
 // (PaePathFor). Named by `path`.
 RugnuxModel ReadRugnuxModelFile(const std::string &path);
diff --git a/writer/RugnuxReflectionFile.cpp b/writer/RugnuxReflectionFile.cpp
new file mode 100644
index 000000000..5fd297502
--- /dev/null
+++ b/writer/RugnuxReflectionFile.cpp
@@ -0,0 +1,247 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "RugnuxReflectionFile.h"
+
+#include <algorithm>
+#include <cctype>
+#include <cmath>
+#include <cstring>
+#include <fstream>
+#include <iterator>
+#include <stdexcept>
+
+#include <gemmi/cif2mtz.hpp>
+#include <gemmi/read_cif.hpp>
+#include <gemmi/refln.hpp>
+
+#include "RugnuxModelFile.h"   // GunzipIfCompressed
+
+namespace {
+
+std::string Lower(std::string s) {
+    std::transform(s.begin(), s.end(), s.begin(), [](unsigned char c) { return std::tolower(c); });
+    return s;
+}
+
+// The sigma of column `c`: the type-Q column named SIG<label>, else the type-Q column right after it.
+const gemmi::Mtz::Column *SigmaOf(const gemmi::Mtz &mtz, const gemmi::Mtz::Column &c, char type) {
+    if (const gemmi::Mtz::Column *s = mtz.column_with_label("SIG" + c.label, nullptr, type))
+        return s;
+    if (c.idx + 1 < mtz.columns.size() && mtz.columns[c.idx + 1].type == type)
+        return &mtz.columns[c.idx + 1];
+    return nullptr;
+}
+
+// Labels of amplitudes that are not observations: calculated structure factors and map coefficients.
+bool CalculatedAmplitude(const std::string &label) {
+    const std::string l = Lower(label);
+    for (const char *p : {"fc", "fwt", "delfwt", "2fofc", "fofc", "f-model", "fmodel", "fcalc", "fpart", "fmask",
+                          "fom", "fdiff", "fano", "dano", "fem"})
+        if (l.starts_with(p))
+            return true;
+    return false;
+}
+
+// A value column of `type` with its sigma: the names in `preferred` first, then any (that `accept` takes).
+template <typename Accept>
+std::pair<const gemmi::Mtz::Column *, const gemmi::Mtz::Column *>
+FindPair(const gemmi::Mtz &mtz, char type, std::initializer_list<const char *> preferred, Accept accept) {
+    for (const char *label : preferred)
+        if (const gemmi::Mtz::Column *c = mtz.column_with_label(label, nullptr, type))
+            if (const gemmi::Mtz::Column *s = SigmaOf(mtz, *c, 'Q'))
+                return {c, s};
+    for (const gemmi::Mtz::Column &c : mtz.columns)
+        if (c.type == type && accept(c.label))
+            if (const gemmi::Mtz::Column *s = SigmaOf(mtz, c, 'Q'))
+                return {&c, s};
+    return {nullptr, nullptr};
+}
+
+// The free-set flag column: FreeR_flag and the other usual names first, then any type-I column with "free" in its name.
+const gemmi::Mtz::Column *FindFree(const gemmi::Mtz &mtz) {
+    for (const char *label : {"FreeR_flag", "FREE", "FREER", "RFREE", "R-free-flags", "FreeRflag", "FREE_FLAG"})
+        if (const gemmi::Mtz::Column *c = mtz.column_with_label(label, nullptr, 'I'))
+            return c;
+    for (const gemmi::Mtz::Column &c : mtz.columns)
+        if (c.type == 'I' && Lower(c.label).find("free") != std::string::npos)
+            return &c;
+    return nullptr;
+}
+
+// The anomalous pair of `value_type` (G or K) with sigmas of `sigma_type` (L or M), named X(+) and X(-).
+bool FindAnomalous(const gemmi::Mtz &mtz, char value_type, char sigma_type, const gemmi::Mtz::Column *out[4]) {
+    for (const gemmi::Mtz::Column &c : mtz.columns) {
+        if (c.type != value_type || !c.label.ends_with("(+)"))
+            continue;
+        const std::string base = c.label.substr(0, c.label.size() - 3);
+        const gemmi::Mtz::Column *m = mtz.column_with_label(base + "(-)", nullptr, value_type);
+        const gemmi::Mtz::Column *sp = SigmaOf(mtz, c, sigma_type);
+        const gemmi::Mtz::Column *sm = m != nullptr ? SigmaOf(mtz, *m, sigma_type) : nullptr;
+        if (m != nullptr && sp != nullptr && sm != nullptr) {
+            out[0] = &c;
+            out[1] = sp;
+            out[2] = m;
+            out[3] = sm;
+            return true;
+        }
+    }
+    return false;
+}
+
+// gemmi's SF-mmCIF -> MTZ conversion with _refln.status preferred to _refln.pdbx_r_free_flag for the free set
+// where status marks any reflection free ('f'), as LoadFCalcFromMtz.cpp reads a deposition.
+std::vector<std::string> SfCifSpec(const gemmi::cif::Loop &loop) {
+    bool status_has_free = false;
+    const int status = loop.find_tag("_refln.status");
+    if (status >= 0)
+        for (size_t i = status; i < loop.values.size(); i += loop.width())
+            if (gemmi::cif::as_string(loop.values[i]) == "f") {
+                status_has_free = true;
+                break;
+            }
+    std::vector<std::string> spec;
+    for (const char **line = gemmi::CifToMtz::default_spec(true); *line != nullptr; ++line) {
+        const std::string l = *line;
+        if ((status_has_free && l.starts_with("pdbx_r_free_flag ")) || (!status_has_free && l.starts_with("status ")))
+            continue;
+        spec.push_back(l);
+    }
+    return spec;
+}
+
+}  // namespace
+
+RugnuxReflections RugnuxReflectionsFromMtz(const gemmi::Mtz &mtz, const std::string &name) {
+    RugnuxReflections out;
+    out.name = name;
+    out.cell = mtz.get_cell();
+    if (mtz.spacegroup == nullptr)
+        throw std::runtime_error(name + " names no space group");
+    out.space_group = mtz.spacegroup->xhm();
+    for (const gemmi::Mtz::Dataset &ds : mtz.datasets)
+        if (ds.wavelength > 0)
+            out.wavelength_A = ds.wavelength;
+
+    auto [value, sigma] = FindPair(mtz, 'F', {"F", "FP", "FOBS", "F-obs", "FMEAN", "F_meas", "Fobs"},
+                                   [](const std::string &l) { return !CalculatedAmplitude(l); });
+    if (value == nullptr) {
+        std::tie(value, sigma) = FindPair(mtz, 'J', {"IMEAN", "I", "IOBS", "I-obs", "Iobs", "I_meas"},
+                                          [](const std::string &) { return true; });
+        out.amplitudes = false;
+    }
+    const gemmi::Mtz::Column *anom[4] = {nullptr, nullptr, nullptr, nullptr};
+    if (value == nullptr) {
+        out.amplitudes = FindAnomalous(mtz, 'G', 'L', anom);
+        if (!out.amplitudes && !FindAnomalous(mtz, 'K', 'M', anom))
+            throw std::runtime_error(name + " holds no amplitudes or intensities with their sigmas");
+        out.value_column = anom[0]->label + "," + anom[2]->label;
+        out.sigma_column = anom[1]->label + "," + anom[3]->label;
+    } else {
+        out.value_column = value->label;
+        out.sigma_column = sigma->label;
+    }
+
+    const gemmi::Mtz::Column *free = FindFree(mtz);
+    std::vector<int> raw_free;
+    const size_t stride = mtz.columns.size();
+    for (size_t r = 0; r < static_cast<size_t>(mtz.nreflections); r++) {
+        const float *row = &mtz.data[r * stride];
+        float v = NAN, s = NAN;
+        if (value != nullptr) {
+            v = row[value->idx];
+            s = row[sigma->idx];
+        } else {
+            // the mean of the two hands where both were measured, the one measured otherwise
+            const float vp = row[anom[0]->idx], sp = row[anom[1]->idx], vm = row[anom[2]->idx], sm = row[anom[3]->idx];
+            if (std::isfinite(vp) && std::isfinite(vm)) {
+                v = 0.5f * (vp + vm);
+                s = 0.5f * std::sqrt(sp * sp + sm * sm);
+            } else if (std::isfinite(vp)) {
+                v = vp;
+                s = sp;
+            } else {
+                v = vm;
+                s = sm;
+            }
+        }
+        if (!std::isfinite(v))
+            continue;
+        MergedReflection m;
+        m.h = static_cast<int32_t>(row[0]);
+        m.k = static_cast<int32_t>(row[1]);
+        m.l = static_cast<int32_t>(row[2]);
+        m.d = static_cast<float>(out.cell.calculate_d({{m.h, m.k, m.l}}));
+        if (out.amplitudes) {
+            m.F = v;
+            m.sigmaF = s;
+        } else {
+            m.I = v;
+            m.sigma = s;
+        }
+        out.reflections.push_back(m);
+        if (free != nullptr) {
+            const float f = row[free->idx];
+            raw_free.push_back(std::isfinite(f) ? static_cast<int>(std::lround(f)) : -1);
+        }
+    }
+    if (out.reflections.empty())
+        throw std::runtime_error(name + " holds no reflection with a value in " + out.value_column);
+
+    if (free != nullptr) {
+        size_t zeros = 0, flagged = 0;
+        for (int f : raw_free) {
+            zeros += f == 0;
+            flagged += f >= 0;
+        }
+        const bool free_is_zero = zeros * 2 <= flagged;
+        for (size_t i = 0; i < out.reflections.size(); i++) {
+            const bool is_free = raw_free[i] >= 0 && (raw_free[i] == 0) == free_is_zero;
+            out.reflections[i].rfree_flag = is_free;
+            out.n_free += is_free;
+        }
+        // a column that marks nothing free is no free set
+        if (out.n_free > 0)
+            out.free_column = free->label;
+    }
+    return out;
+}
+
+RugnuxReflections ParseRugnuxReflections(const std::string &bytes, const std::string &name) {
+    const std::string text = GunzipIfCompressed(bytes);
+    if (text.size() >= 4 && std::memcmp(text.data(), "MTZ ", 4) == 0) {
+        gemmi::Mtz mtz;
+        gemmi::MemoryStream stream(text.data(), text.size());
+        mtz.read_stream(stream, true);
+        RugnuxReflections out = RugnuxReflectionsFromMtz(mtz, name);
+        out.source = "MTZ";
+        return out;
+    }
+    gemmi::cif::Document doc = gemmi::read_cif_from_memory(text.data(), text.size(), name.c_str());
+    std::vector<gemmi::ReflnBlock> blocks = gemmi::as_refln_blocks(std::move(doc.blocks));
+    gemmi::Logger logger;
+    std::string why = name + " holds no merged reflection block";
+    for (gemmi::ReflnBlock &rb : blocks) {
+        if (rb.refln_loop == nullptr)
+            continue;
+        gemmi::CifToMtz converter;
+        converter.spec_lines = SfCifSpec(*rb.refln_loop);
+        try {
+            const gemmi::Mtz mtz = converter.convert_block_to_mtz(rb, logger);
+            RugnuxReflections out = RugnuxReflectionsFromMtz(mtz, name);
+            out.source = "mmCIF block " + rb.block.name;
+            return out;
+        } catch (const std::exception &e) {
+            why = e.what();
+        }
+    }
+    throw std::runtime_error(why);
+}
+
+RugnuxReflections ReadRugnuxReflectionFile(const std::string &path) {
+    std::ifstream f(path, std::ios::binary);
+    if (!f)
+        throw std::runtime_error("cannot open " + path);
+    const std::string bytes{std::istreambuf_iterator<char>(f), std::istreambuf_iterator<char>()};
+    return ParseRugnuxReflections(bytes, path);
+}
diff --git a/writer/RugnuxReflectionFile.h b/writer/RugnuxReflectionFile.h
new file mode 100644
index 000000000..08a2de429
--- /dev/null
+++ b/writer/RugnuxReflectionFile.h
@@ -0,0 +1,50 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <string>
+#include <vector>
+
+#include "gemmi/mtz.hpp"
+#include "gemmi/unitcell.hpp"
+
+#include "../common/Reflection.h"
+
+// Merged reflections a refinement is given (rugnux --mode refine), as their file held them: an MTZ or a
+// structure-factor mmCIF (as the PDB holds them), gzipped or not, the format taken from the content. The
+// columns are found, not named: amplitudes with their sigmas where the file has them, else intensities (the
+// caller turns those into amplitudes), and the free-set flags where there are any. The caller parses; the
+// analysis reads no file.
+struct RugnuxReflections {
+    std::string name;                // what the report calls the file
+    std::string source;              // "MTZ" or "mmCIF block <name>"
+    gemmi::UnitCell cell;
+    std::string space_group;         // extended Hermann-Mauguin symbol
+    double wavelength_A = 0;         // 0: the file does not say
+    // h, k, l and d, with F/sigmaF where amplitudes were read, I/sigma where intensities were; rfree_flag
+    // where the file has free-set flags.
+    std::vector<MergedReflection> reflections;
+    bool amplitudes = true;          // false: intensities, no amplitudes in the file
+    std::string value_column, sigma_column;
+    std::string free_column;         // empty: the file has no free-set flags
+    int n_free = 0;
+};
+
+// The columns of `mtz` (an MTZ, or an mmCIF converted to one), found as RugnuxReflections describes:
+// - amplitudes: a type-F column followed by its type-Q sigma, by the usual names first (F, FP, FOBS, F-obs,
+//   FMEAN, ...), then any such pair that is not a calculated or map coefficient (FC, FWT, DELFWT, F-model ...);
+// - else intensities: a type-J column with its type-Q sigma (IMEAN, I, IOBS, I-obs, ... first); the
+//   anomalous pairs F(+)/F(-) (G/L) or I(+)/I(-) (K/M) averaged where there is nothing else;
+// - free flags: a type-I column named as a free set (FreeR_flag, FREE, R-free-flags, ...: any with "free" in
+//   its name). The test set is flag 0 in the CCP4 convention (0-19 thin shells included); where 0 marks the
+//   majority (1 marks the test set, as phenix and CNS write it), the minority is taken.
+// Throws where there is neither.
+RugnuxReflections RugnuxReflectionsFromMtz(const gemmi::Mtz &mtz, const std::string &name);
+
+// The reflections of a file's bytes, the format from the content (MTZ by its magic, else mmCIF: the first
+// merged block that holds amplitudes or intensities). Throws where the bytes hold none.
+RugnuxReflections ParseRugnuxReflections(const std::string &bytes, const std::string &name);
+
+// The file at `path`, named by it.
+RugnuxReflections ReadRugnuxReflectionFile(const std::string &path);