From 4120b0daa8c53aaf8fc3eb5acb9aae1ca5dcb606 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 18:51:56 +0200 Subject: [PATCH] Merged mmCIF: drop the scaled unmerged block The unmerged observations already leave a run as the unmerged MTZ (and the scaled fulls as the .hkl), and deposition does not come straight from rugnux, so the _diffrn_refln block made the .cif up to ten times larger for no consumer. The .cif is one data block again; the keys, software rows, metadata, R_merge/R_pim/CC* and provenance stay. ScaledFull loses the scale and frame it carried only for that block. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- common/Reflection.h | 5 +- docs/RUGNUX_ADVANCED.md | 30 +++------ image_analysis/WriteReflections.cpp | 67 +------------------ image_analysis/WriteReflections.h | 7 +- .../scale_merge/RotationScaleMerge.cpp | 3 +- tests/WriteReflectionsTest.cpp | 24 +------ 6 files changed, 21 insertions(+), 115 deletions(-) diff --git a/common/Reflection.h b/common/Reflection.h index cc19f303f..a6b5ad9c5 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -68,16 +68,13 @@ struct Reflection { // measurement, every correction and the per-frame scale applied, its sigma under the error model the // merge weighted it with - but NOT averaged with its symmetry mates. h k l are the index it was // measured at, not its ASU representative. Written as the SHELX HKLF 4 file, so that the program -// reading it sees the equivalents and computes Rint itself, and as the scaled unmerged block of the -// mmCIF. +// reading it sees the equivalents and computes Rint itself. struct ScaledFull { int32_t h = 0; int32_t k = 0; int32_t l = 0; float I = NAN; float sigma = NAN; - float scale = NAN; // the factor the scaling applied: I is the integrated intensity times this - float image_number = NAN; // the frame the full peaks on (0-based, fractional) }; struct MergedReflection { diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index 75e2a69a2..72a683ba6 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -245,25 +245,17 @@ exact `-S`/`-C` override to merge it in its own setting, instead of failing insi ## The merged mmCIF (`.cif`) -The `.cif` is written for deposition and follows the PDBx/mmCIF dictionary. It holds two data -blocks: - -- **`data_sample`** — the merged reflections (`_refln`: I, F and the Bijvoet mates, with the free set - in `_refln.status`) and what describes them: the cell and space group; `_software` with Rugnux - named twice, once for `data reduction` and once for `data scaling`; the experiment as the input - file states it (`_diffrn.ambient_temp` from the sample temperature, the wavelength, the source's - facility and beamline and its general class, the detector's make and model in - `_diffrn_detector.type` and its class in `_diffrn_detector.detector`, and the collection date from - the NXmx start time); and the merging statistics overall (`_reflns`) and per shell - (`_reflns_shell`): completeness, multiplicity, , R_merge, R_meas, R_pim, CC1/2 and CC\*, - and SigAno where an anomalous split was made. A quantity the input does not state is left out - rather than guessed. -- **`data_sample_unmerged`** (rotation) — the scaled **unmerged** observations the merge averaged, - as `_diffrn_refln`: the index each was measured at, the frame it peaks on (`pdbx_image_id`, - numbered from 1 as the unmerged MTZ's `BATCH`), the scaled intensity and its σ under the error - model, and the scale factor applied (`pdbx_scale_value`). These are the same records as the `.hkl`. - The wwPDB accepts scaled unmerged data beside the merged set. The `_unmerged.mtz` is a different - thing and stays unscaled, for the programs that fit a scale model of their own. +The `.cif` follows the PDBx/mmCIF dictionary and holds one data block, `data_sample`: the merged +reflections (`_refln`: I, F and the Bijvoet mates, with the free set in `_refln.status`) and what +describes them - the cell and space group; `_software` with Rugnux named twice, once for +`data reduction` and once for `data scaling`; the experiment as the input file states it +(`_diffrn.ambient_temp` from the sample temperature, the wavelength, the source's facility and +beamline and its general class, the detector's make and model in `_diffrn_detector.type` and its +class in `_diffrn_detector.detector`, and the collection date from the NXmx start time); and the +merging statistics overall (`_reflns`) and per shell (`_reflns_shell`): completeness, multiplicity, +, R_merge, R_meas, R_pim, CC1/2 and CC\*, and SigAno where an anomalous split was made. A +quantity the input does not state is left out rather than guessed. The unmerged observations are in +the `_unmerged.mtz` (unscaled) and the `.hkl` (scaled), not in the `.cif`. Items PDBx does not define — ISa, the error model, the twinning and anisotropy diagnostics, the refined distance and beam centre — are written under the `jfjoch_` local prefix, which a dictionary diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 179c137e5..7b49f877a 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -140,7 +140,7 @@ void AddProvenanceHistory(gemmi::Mtz &mtz, const OutputProvenance &provenance) { add("Collected: ", provenance.collection_date); } -// _entry, _cell and _symmetry: what every data block of the file restates, so each one stands alone. +// _entry, _cell and _symmetry. void WriteBlockHeader(std::ostream &out, const std::string &entry_id, const UnitCell &unitCell, const gemmi::SpaceGroup &sg) { out << "_entry.id " << entry_id << "\n"; @@ -253,7 +253,6 @@ void WriteMmcifReflections(const std::vector &reflections, const MergeStatistics &statistics, const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, - const std::vector &scaled_fulls, const OutputProvenance &provenance, const std::string &filename, size_t nthreads) { @@ -266,7 +265,7 @@ void WriteMmcifReflections(const std::vector &reflections, // The keys PDBx asks for are written with the one value this file has for each: one entry, one // crystal, one diffraction experiment (diffrn_id 1, which _pdbx_reflns_twin below also points at) - // and one wavelength. The merged data come first, as the wwPDB expects of a structure-factor file. + // and one wavelength. const std::string entry_id = "sample"; const gemmi::SpaceGroup &sg = experiment.GetSpaceGroupOrP1(); @@ -603,66 +602,6 @@ void WriteMmcifReflections(const std::vector &reflections, out << "#\n"; out << "# End of reflections\n"; - // ---------- the scaled unmerged observations, as a second data block ---------- - // The observations the merge above averaged - each full with the scale and the error-model sigma - // the merge gave it, at the index it was measured at - which is what the wwPDB accepts as scaled - // unmerged data. The same records as the .hkl: an observation without a usable sigma is in - // neither. pdbx_image_id is the frame the full peaks on, numbered from 1 as the unmerged MTZ's - // BATCH is. (The unmerged MTZ itself stays unscaled, for the programs that fit their own scales.) - if (!scaled_fulls.empty()) { - out << "data_sample_unmerged\n"; - out << "#\n"; - WriteBlockHeader(out, entry_id, unitCell, sg); - out << "_exptl_crystal.id 1\n"; - out << "#\n"; - out << "_diffrn.id 1\n"; - out << "_diffrn.crystal_id 1\n"; - out << "#\n"; - out << "_diffrn_radiation_wavelength.id 1\n"; - out << "_diffrn_radiation_wavelength.wavelength " << Fmt(experiment.GetWavelength_A(), 5) << "\n"; - out << "#\n"; - out << "loop_\n"; - out << "_diffrn_refln.diffrn_id\n"; - out << "_diffrn_refln.id\n"; - out << "_diffrn_refln.pdbx_image_id\n"; - out << "_diffrn_refln.index_h\n"; - out << "_diffrn_refln.index_k\n"; - out << "_diffrn_refln.index_l\n"; - out << "_diffrn_refln.intensity_net\n"; - out << "_diffrn_refln.intensity_sigma\n"; - out << "_diffrn_refln.pdbx_scale_value\n"; - - // Numbered in order over the usable fulls, so the id of each is known before the rows are - // formatted on several workers. - std::vector id(scaled_fulls.size() + 1, 0); - for (size_t i = 0; i < scaled_fulls.size(); ++i) { - const auto &f = scaled_fulls[i]; - id[i + 1] = id[i] + ((std::isfinite(f.I) && std::isfinite(f.sigma) && f.sigma > 0.0f) ? 1 : 0); - } - const size_t nrow = scaled_fulls.size(); - const size_t nw = std::max(nthreads, 1); - const int nch = static_cast(ThreadsForWork(nrow, nw, 4096)); - std::vector block(nch); - ParallelChunks(nch, nw, [&](int tlo, int thi) { - for (int t = tlo; t < thi; ++t) { - const size_t lo = nrow * t / nch, hi = nrow * (t + 1) / nch; - std::string &s = block[t]; - s.reserve((hi - lo) * 64); - for (size_t i = lo; i < hi; ++i) { - if (id[i + 1] == id[i]) - continue; - const auto &f = scaled_fulls[i]; - const long frame = std::isfinite(f.image_number) ? 1 + std::lround(f.image_number) : 0; - s += "1 " + std::to_string(id[i + 1]) + " " + (frame > 0 ? std::to_string(frame) : "?") + " " - + std::to_string(f.h) + " " + std::to_string(f.k) + " " + std::to_string(f.l) + " " - + Fmt(f.I, 3) + " " + Fmt(f.sigma, 3) + " " + Fmt(f.scale, 5) + "\n"; - } - } - }); - for (const std::string &s : block) - out.write(s.data(), static_cast(s.size())); - out << "#\n"; - } out.close(); } @@ -1296,7 +1235,7 @@ void WriteReflections(const std::vector &reflections, // fulls where the merge provides them (rotation), the merged reflections otherwise (stills). WriteMtzReflections(reflections, unitCell, experiment, provenance, filename + ".mtz"); WriteMmcifReflections(reflections, unitCell, experiment, statistics, error_model, twinning, - scaled_fulls ? *scaled_fulls : std::vector{}, provenance, + provenance, filename + ".cif", nthreads); if (scaled_fulls && !scaled_fulls->empty()) WriteShelxHklReflections(*scaled_fulls, filename + ".hkl", nthreads); diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 6fc02f042..559b3e662 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -37,17 +37,14 @@ struct OutputProvenance { std::string collection_date; // as the input states it (NXmx /entry/start_time), ISO-8601 }; -// The merged reflections, their statistics and the experiment, as one PDBx/mmCIF data block; and, when -// scaled_fulls holds any, the same observations unmerged and scaled (_diffrn_refln) as a second block - -// the scaled unmerged data the wwPDB accepts beside the merged set. nthreads: workers for the -// per-reflection row formatting, which is the bulk of the file. +// The merged reflections, their statistics and the experiment, as one PDBx/mmCIF data block. +// nthreads: workers for the per-reflection row formatting, which is the bulk of the file. void WriteMmcifReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, const MergeStatistics &statistics, const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, - const std::vector &scaled_fulls, const OutputProvenance &provenance, const std::string &filename, size_t nthreads); diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 8e6c9c39b..144e85d5a 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -5409,8 +5409,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool for (int j = lo; j < hi; ++j) { const Obs &o = fulls[kept[j]]; const float I_corr = o.I * o.corr; - result.scaled_fulls[j] = {o.h, o.k, o.l, I_corr, corrected_sigma(o, I_corr, o.sigma * o.corr), - o.corr, o.image_number}; + result.scaled_fulls[j] = {o.h, o.k, o.l, I_corr, corrected_sigma(o, I_corr, o.sigma * o.corr)}; } }); } diff --git a/tests/WriteReflectionsTest.cpp b/tests/WriteReflectionsTest.cpp index c66b4a33f..848a3d02c 100644 --- a/tests/WriteReflectionsTest.cpp +++ b/tests/WriteReflectionsTest.cpp @@ -232,7 +232,7 @@ TEST_CASE("CC* from CC1/2", "[write_reflections][portable]") { CHECK(std::isnan(CCStar(NAN))); } -TEST_CASE("Merged mmCIF: PDBx keys, merging R-factors and the scaled unmerged block", +TEST_CASE("Merged mmCIF: PDBx keys, merging R-factors and provenance", "[write_reflections][portable]") { jfjoch_test::SyntheticMergeParams params; params.true_space_group = "P 43 21 2"; @@ -254,12 +254,6 @@ TEST_CASE("Merged mmCIF: PDBx keys, merging R-factors and the scaled unmerged bl statistics.shells = {shell}; statistics.overall = shell; - // Three fulls, one without a usable sigma - that one is left out, as from the .hkl. - const std::vector fulls{ - {1, 2, 3, 200.0f, 20.0f, 1.25f, 0.4f}, - {-1, -2, -3, 100.0f, 12.0f, 0.8f, 9.6f}, - {4, 5, 6, 10.0f, 0.0f, 1.0f, 3.0f}, - }; OutputProvenance provenance; provenance.collection_date = "2026-07-13T13:48:52.413Z"; provenance.command_line = "rugnux -o 'my run' /data/x_master.h5"; @@ -268,11 +262,11 @@ TEST_CASE("Merged mmCIF: PDBx keys, merging R-factors and the scaled unmerged bl const auto path = (std::filesystem::temp_directory_path() / "rugnux_merged_test.cif").string(); WriteMmcifReflections(reflections, TETRAGONAL_CELL, experiment, statistics, ErrorModelReport{}, - TwinningAnalysisResult{}, fulls, provenance, path, 2); + TwinningAnalysisResult{}, provenance, path, 2); gemmi::cif::Document doc = gemmi::read_cif_gz(path); std::filesystem::remove(path); - REQUIRE(doc.blocks.size() == 2); + REQUIRE(doc.blocks.size() == 1); gemmi::cif::Block &merged = doc.blocks[0]; CHECK(merged.name == "sample"); CHECK(*merged.find_value("_reflns.pdbx_ordinal") == "1"); @@ -301,18 +295,6 @@ TEST_CASE("Merged mmCIF: PDBx keys, merging R-factors and the scaled unmerged bl CHECK(shells[0][1] == "0.0500"); CHECK(shells[0][2] == "0.0300"); - gemmi::cif::Block &unmerged = doc.blocks[1]; - CHECK(unmerged.name == "sample_unmerged"); - gemmi::cif::Table obs = unmerged.find("_diffrn_refln.", {"id", "pdbx_image_id", "index_h", - "intensity_net", "intensity_sigma", - "pdbx_scale_value"}); - REQUIRE(obs.length() == 2); - CHECK(obs[0][0] == "1"); - CHECK(obs[0][1] == "1"); // peaks on frame 0.4 -> image 1 - CHECK(obs[1][1] == "11"); // 9.6 -> image 11 - CHECK(obs[1][2] == "-1"); - CHECK(gemmi::cif::as_number(obs[0][3]) == Catch::Approx(200.0)); - CHECK(gemmi::cif::as_number(obs[1][5]) == Catch::Approx(0.8)); } TEST_CASE("MTZ history: command line and input, in records of at most 80 characters",