From fc6b91555f2d674e6b0c39e81154a66a1014ef09 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 17:13:27 +0200 Subject: [PATCH] Merged mmCIF: deposition keys, R_merge/R_pim/CC*, experiment metadata, scaled unmerged block The .cif now follows the PDBx/mmCIF dictionary's keys: _entry.id, the entry_id of _cell/_symmetry/_reflns, pdbx_ordinal and pdbx_diffrn_id on _reflns and _reflns_shell, diffrn_id on _diffrn_source/_diffrn_detector/ _diffrn_radiation, and _diffrn_radiation_wavelength.id. _software is a loop naming Rugnux for 'data reduction' and for 'data scaling'. The experiment as the input states it: _diffrn.ambient_temp from the sample temperature, the source's general class (from NXmx /entry/source/type, which the HDF5 reader now reads), facility and beamline, the detector make/model in _diffrn_detector.type (where .detector used to carry it) and its class in .detector where the name settles it, and the collection date from /entry/start_time. R_merge and R_pim are formed beside R_meas from the same per-reflection sums (stills and rotation merges) and written with CC* overall and per shell (_reflns.pdbx_Rmerge_I_obs/pdbx_Rpim_I_all/pdbx_CC_star and the _reflns_shell counterparts), as report keys R_MERGE/R_PIM/CC_STAR and as three new columns of the shell tables. R_pim following Weiss (2001), CC* following Karplus & Diederichs (2012). A second data block, data_sample_unmerged, carries the scaled unmerged observations of a rotation merge as _diffrn_refln (measured index, frame, scaled I and sigma, pdbx_scale_value) - the same records as the .hkl, which the wwPDB accepts as scaled unmerged data. ScaledFull now carries the scale and the frame. The unmerged MTZ stays unscaled. Reflection data are unchanged: the merged _refln rows, the .mtz columns and the .hkl are written as before. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- common/Reflection.h | 5 +- docs/ACKNOWLEDGEMENT.md | 5 +- docs/CPU_DATA_ANALYSIS.md | 1 + docs/RUGNUX_ADVANCED.md | 26 +++ docs/RUGNUX_REPORT.md | 7 + image_analysis/WriteReflections.cpp | 210 ++++++++++++++++-- image_analysis/WriteReflections.h | 14 +- image_analysis/scale_merge/Merge.cpp | 29 ++- image_analysis/scale_merge/Merge.h | 12 + .../scale_merge/RotationScaleMerge.cpp | 16 +- reader/HDF5MetadataSource.cpp | 1 + rugnux/ResultReport.cpp | 6 + rugnux/Rugnux.h | 3 + rugnux/RugnuxScaleMerge.cpp | 2 +- rugnux/rugnux_cli.cpp | 7 +- tests/WriteReflectionsTest.cpp | 86 +++++++ 16 files changed, 395 insertions(+), 35 deletions(-) diff --git a/common/Reflection.h b/common/Reflection.h index a6b5ad9c5..cc19f303f 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -68,13 +68,16 @@ 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. +// reading it sees the equivalents and computes Rint itself, and as the scaled unmerged block of the +// mmCIF. 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/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index ca82c3cbc..c2de28fed 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -253,7 +253,10 @@ reporting of I/sigma(I). K. Diederichs and P. A. Karplus, "Improved R-factors for diffraction data analysis in macromolecular crystallography" (1997), Nat. Struct. Biol. 4, 269-275 [doi:10.1038/nsb0497-269](https://doi.org/10.1038/nsb0497-269); P. A. Karplus and K. Diederichs, "Linking crystallographic model and data quality" (2012), Science 336, 1030-1033 -[doi:10.1126/science.1218231](https://doi.org/10.1126/science.1218231); K. Diederichs and +[doi:10.1126/science.1218231](https://doi.org/10.1126/science.1218231); M. S. Weiss, "Global +indicators of X-ray data quality" (2001), J. Appl. Cryst. 34, 130-135 +[doi:10.1107/S0021889800018227](https://doi.org/10.1107/S0021889800018227) (R_pim beside R_meas and +R_merge); K. Diederichs and P. A. Karplus, "Better models by discarding data?" (2013), Acta Cryst. D69, 1215-1222 [doi:10.1107/S0907444913001121](https://doi.org/10.1107/S0907444913001121). diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index bfb7de4f6..b5b39ccce 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -76,6 +76,7 @@ The methods draw on, and in places reimplement, solutions from: - P. R. Evans & G. N. Murshudov, "How good are my data and what is the resolution?", *Acta Cryst.* **D69** (2013), 1204-1214 (AIMLESS: the anisotropic deltaB as the range of the principal components, and diffraction limits from a cone about each principal direction). - G. Assmann, W. Brehm & K. Diederichs, "Identification of rogue datasets in serial crystallography", *J. Appl. Cryst.* **49** (2016), 1021-1028, and G. M. Assmann, M. Wang & K. Diederichs, *Acta Cryst.* **D76** (2020), 636-652 (XDSCC12: sigma-tau CC1/2, delta-CC1/2, the Fisher transformation and the rejection discipline the frame disposition follows). - K. Diederichs & P. A. Karplus, *Nat. Struct. Biol.* **4** (1997), 269-275, and P. A. Karplus & K. Diederichs, *Science* **336** (2012), 1030-1033 (R_meas / R_pim, CC1/2 and CC\*). +- M. S. Weiss, "Global indicators of X-ray data quality", *J. Appl. Cryst.* **34** (2001), 130-135 (R_pim). - IUCr Commission on Crystallographic Nomenclature, "Statistical descriptors in crystallography", *Acta Cryst.* **A45** (1989), 63-75, and *Acta Cryst.* **A51** (1995), 565-569 (uncertainty conventions). (list is not exhaustive; the full citations, with DOIs, are in [ACKNOWLEDGEMENT.md](ACKNOWLEDGEMENT.md)) diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index 2a42120c3..6abcea11e 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -243,6 +243,32 @@ describe the same frame. An older file that was affected by this cannot be repai not recoverable after the fact); such a file now stops with a message naming both cells and the exact `-S`/`-C` override to merge it in its own setting, instead of failing inside the merge. +## 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. + +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 +check reports as unknown and a deposition ignores. + ## Command-line options General: diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index ffaacd498..8f553e584 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -60,6 +60,13 @@ scaled-up counting error gives it, and is in `R_MEAS_WEIGHTED=` at that weight t the number with its noise; a large gap between the two says weak frames are kept at low weight. `MULTIPLICITY=` stays a count of observations. +**`R_MERGE=`** and **`R_PIM=`** are the other two merging R-factors over the same observations as +`R_MEAS=`: R_merge without the multiplicity correction, and the precision-indicating R_pim (Weiss, +2001), which falls as the multiplicity rises. **`CC_STAR=`** is CC\* = sqrt(2 CC1/2 / (1 + CC1/2)) +(Karplus & Diederichs, 2012), the estimated correlation of the merged data with the true signal; it is +absent where CC1/2 is not positive. All three are also columns of the shell tables (`R_merge`, +`R_pim`, `CC*`) and are written per shell into the merged mmCIF. + **`OBSERVATIONS_REJECTED=`** counts the observations the merge's outlier tests removed, and **`OBSERVATIONS_REJECTED_WILSON=`** the part of them the Wilson test took (each listed in the developer report). **`OBSERVATIONS_REJECTED_OVERLOAD=`** is separate from both: the rotation measurements left diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index dd31bae2b..da0fbe3dd 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include #include @@ -68,6 +69,66 @@ std::string CifStr(const std::string& s) { return s; } +std::string Upper(std::string s) { + for (char &c : s) + c = static_cast(std::toupper(static_cast(c))); + return s; +} + +// The general class of the detector, _diffrn_detector.detector in the PDBx enumeration, where its make +// and model name it; empty where they do not. +std::string DetectorClass(const std::string &description) { + const std::string d = Upper(description); + if (d.find("CCD") != std::string::npos) + return "CCD"; + if (d.find("IMAGE PLATE") != std::string::npos) + return "IMAGE PLATE"; + for (const char *pixel : {"EIGER", "PILATUS", "JUNGFRAU"}) + if (d.find(pixel) != std::string::npos) + return "PIXEL"; + return ""; +} + +// The general class of the source, _diffrn_source.source in the PDBx enumeration, from the NXmx +// /entry/source/type the file states; empty where it states none or one PDBx has no class for. +std::string SourceClass(const std::string &nxmx_type) { + const std::string t = Upper(nxmx_type); + if (t.find("SYNCHROTRON") != std::string::npos) return "SYNCHROTRON"; + if (t.find("FREE-ELECTRON") != std::string::npos || t.find("FREE ELECTRON") != std::string::npos) + return "FREE ELECTRON LASER"; + if (t.find("ROTATING ANODE") != std::string::npos) return "ROTATING ANODE"; + if (t.find("FIXED TUBE") != std::string::npos) return "SEALED TUBE"; + if (t.find("METAL JET") != std::string::npos) return "LIQUID ANODE"; + return ""; +} + +// The date part of an ISO-8601 time stamp (YYYY-MM-DD...), which is all the PDBx collection date +// takes; empty for anything else. +std::string IsoDate(const std::string ×tamp) { + if (timestamp.size() < 10 || timestamp[4] != '-' || timestamp[7] != '-') + return ""; + return timestamp.substr(0, 10); +} + +// _entry, _cell and _symmetry: what every data block of the file restates, so each one stands alone. +void WriteBlockHeader(std::ostream &out, const std::string &entry_id, const UnitCell &unitCell, + const gemmi::SpaceGroup &sg) { + out << "_entry.id " << entry_id << "\n"; + out << "#\n"; + out << "_cell.entry_id " << entry_id << "\n"; + out << "_cell.length_a " << Fmt(unitCell.a, 3) << "\n"; + out << "_cell.length_b " << Fmt(unitCell.b, 3) << "\n"; + out << "_cell.length_c " << Fmt(unitCell.c, 3) << "\n"; + out << "_cell.angle_alpha " << Fmt(unitCell.alpha, 2) << "\n"; + out << "_cell.angle_beta " << Fmt(unitCell.beta, 2) << "\n"; + out << "_cell.angle_gamma " << Fmt(unitCell.gamma, 2) << "\n"; + out << "#\n"; + out << "_symmetry.entry_id " << entry_id << "\n"; + out << "_symmetry.space_group_name_H-M " << CifStr(sg.hm) << "\n"; + out << "_symmetry.Int_Tables_number " << sg.number << "\n"; + out << "#\n"; +} + // 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 @@ -162,6 +223,8 @@ 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) { @@ -171,9 +234,16 @@ void WriteMmcifReflections(const std::vector &reflections, out << std::fixed; + // 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. + const std::string entry_id = "sample"; + const gemmi::SpaceGroup &sg = experiment.GetSpaceGroupOrP1(); + // ---------- data block ---------- out << "data_sample" << "\n"; out << "#\n"; + WriteBlockHeader(out, entry_id, unitCell, sg); // ---------- _audit ---------- out << "_audit.revision_id 1\n"; @@ -182,35 +252,54 @@ void WriteMmcifReflections(const std::vector &reflections, out << "#\n"; // ---------- _software ---------- - out << "_software.name 'Rugnux'\n"; - out << "_software.version " << CifStr(jfjoch_version()) << "\n"; - out << "_software.classification 'data reduction'\n"; + // Rugnux integrated the images (data reduction) and scaled and merged them (data scaling), the + // two rows the wwPDB asks a data-processing program to be named under. + out << "loop_\n"; + out << "_software.pdbx_ordinal\n"; + out << "_software.name\n"; + out << "_software.version\n"; + out << "_software.classification\n"; + out << "1 Rugnux " << CifStr(jfjoch_version()) << " 'data reduction'\n"; + out << "2 Rugnux " << CifStr(jfjoch_version()) << " 'data scaling'\n"; out << "#\n"; - // ---------- _cell ---------- - out << "_cell.length_a " << Fmt(unitCell.a, 3) << "\n"; - out << "_cell.length_b " << Fmt(unitCell.b, 3) << "\n"; - out << "_cell.length_c " << Fmt(unitCell.c, 3) << "\n"; - out << "_cell.angle_alpha " << Fmt(unitCell.alpha, 2) << "\n"; - out << "_cell.angle_beta " << Fmt(unitCell.beta, 2) << "\n"; - out << "_cell.angle_gamma " << Fmt(unitCell.gamma, 2) << "\n"; - - const gemmi::SpaceGroup &sg = experiment.GetSpaceGroupOrP1(); - - // ---------- _symmetry ---------- - out << "_symmetry.space_group_name_H-M " << CifStr(sg.hm) << "\n"; - out << "_symmetry.Int_Tables_number " << sg.number << "\n"; + // ---------- _exptl_crystal / _diffrn ---------- + out << "_exptl_crystal.id 1\n"; + out << "#\n"; + out << "_diffrn.id 1\n"; + out << "_diffrn.crystal_id 1\n"; + if (const auto t = experiment.GetSampleTemperature_K()) + out << "_diffrn.ambient_temp " << Fmt(*t, 1) << "\n"; out << "#\n"; - // ---------- _diffrn_source / _diffrn_detector ---------- + // ---------- _diffrn_radiation / _diffrn_source / _diffrn_detector ---------- + out << "_diffrn_radiation.diffrn_id 1\n"; + out << "_diffrn_radiation.wavelength_id 1\n"; + out << "_diffrn_radiation.pdbx_diffrn_protocol 'SINGLE WAVELENGTH'\n"; + out << "#\n"; + out << "_diffrn_radiation_wavelength.id 1\n"; + out << "_diffrn_radiation_wavelength.wavelength " << Fmt(experiment.GetWavelength_A(), 5) << "\n"; + out << "#\n"; + + // The source and the detector as the input file names them: the facility and the beamline (NXmx + // /entry/source/name and /entry/instrument/name) and the detector's make and model. The general + // classes PDBx enumerates are written only where those names settle them. + out << "_diffrn_source.diffrn_id 1\n"; + if (const std::string source = SourceClass(experiment.GetSourceType()); !source.empty()) + out << "_diffrn_source.source " << CifStr(source) << "\n"; if (!experiment.GetSourceName().empty()) out << "_diffrn_source.pdbx_synchrotron_site " << CifStr(experiment.GetSourceName()) << "\n"; - if (!experiment.GetInstrumentName().empty()) out << "_diffrn_source.pdbx_synchrotron_beamline " << CifStr(experiment.GetInstrumentName()) << "\n"; + out << "_diffrn_source.pdbx_wavelength_list " << Fmt(experiment.GetWavelength_A(), 5) << "\n"; + out << "#\n"; - out << "_diffrn_radiation_wavelength.wavelength " << Fmt(experiment.GetWavelength_A(), 5) << "\n"; - out << "_diffrn_detector.detector " << CifStr(experiment.GetDetectorDescription()) << "\n"; + out << "_diffrn_detector.diffrn_id 1\n"; + if (const std::string cls = DetectorClass(experiment.GetDetectorDescription()); !cls.empty()) + out << "_diffrn_detector.detector " << CifStr(cls) << "\n"; + out << "_diffrn_detector.type " << CifStr(experiment.GetDetectorDescription()) << "\n"; + if (const std::string date = IsoDate(provenance.collection_date); !date.empty()) + out << "_diffrn_detector.pdbx_collection_date " << date << "\n"; // Detector geometry actually used for integration (refined, when geometry refinement ran - the rotation // two-pass or the stills global refinement update it on experiment_ before the written pass). jfjoch_ // local-data-name items: feedback of the distance / beam centre the data was reduced with. @@ -236,6 +325,9 @@ void WriteMmcifReflections(const std::vector &reflections, // Anomalous signal-to-noise (SigAno) is written only when an anomalous split was made, so a // non-anomalous merge keeps its previous stats block / shell-loop columns unchanged. const bool has_anom = std::isfinite(ov.abs_diff_over_sigma_anomalous); + out << "_reflns.entry_id " << entry_id << "\n"; + out << "_reflns.pdbx_ordinal 1\n"; + out << "_reflns.pdbx_diffrn_id 1\n"; out << "_reflns.d_resolution_high " << Fmt(ov.d_min, 2) << "\n"; out << "_reflns.d_resolution_low " << Fmt(ov.d_max, 2) << "\n"; out << "_reflns.number_obs " << ov.unique_reflections << "\n"; @@ -243,8 +335,11 @@ void WriteMmcifReflections(const std::vector &reflections, out << "_reflns.pdbx_redundancy " << Fmt(mult(ov), 2) << "\n"; out << "_reflns.percent_possible_obs " << Fmt(compl_pct(ov), 1) << "\n"; out << "_reflns.pdbx_netI_over_sigmaI " << Fmt(ov.mean_i_over_sigma, 2) << "\n"; + out << "_reflns.pdbx_Rmerge_I_obs " << Fmt(ov.r_merge, 4) << "\n"; out << "_reflns.pdbx_Rrim_I_all " << Fmt(ov.r_meas, 4) << "\n"; + out << "_reflns.pdbx_Rpim_I_all " << Fmt(ov.r_pim, 4) << "\n"; out << "_reflns.pdbx_CC_half " << Fmt(ov.cc_half, 4) << "\n"; + out << "_reflns.pdbx_CC_star " << Fmt(CCStar(ov.cc_half), 4) << "\n"; if (has_anom) out << "_reflns.pdbx_absDiff_over_sigma_anomalous " << Fmt(ov.abs_diff_over_sigma_anomalous, 3) << " # SigAno = <|dano|>/\n"; @@ -360,6 +455,8 @@ void WriteMmcifReflections(const std::vector &reflections, } out << "loop_\n"; + out << "_reflns_shell.pdbx_ordinal\n"; + out << "_reflns_shell.pdbx_diffrn_id\n"; out << "_reflns_shell.d_res_high\n"; out << "_reflns_shell.d_res_low\n"; out << "_reflns_shell.number_measured_obs\n"; @@ -367,17 +464,23 @@ void WriteMmcifReflections(const std::vector &reflections, out << "_reflns_shell.pdbx_redundancy\n"; out << "_reflns_shell.percent_possible_obs\n"; out << "_reflns_shell.meanI_over_sigI_obs\n"; + out << "_reflns_shell.Rmerge_I_obs\n"; out << "_reflns_shell.pdbx_Rrim_I_all\n"; + out << "_reflns_shell.pdbx_Rpim_I_all\n"; out << "_reflns_shell.pdbx_CC_half\n"; + out << "_reflns_shell.pdbx_CC_star\n"; if (has_anom) out << "_reflns_shell.pdbx_absDiff_over_sigma_anomalous\n"; + int ordinal = 0; for (const auto &s : statistics.shells) { if (s.unique_reflections == 0) continue; - out << Fmt(s.d_min, 2) << " " << Fmt(s.d_max, 2) << " " + out << ++ordinal << " 1 " + << Fmt(s.d_min, 2) << " " << Fmt(s.d_max, 2) << " " << s.total_observations << " " << s.unique_reflections << " " << Fmt(mult(s), 2) << " " << Fmt(compl_pct(s), 1) << " " - << Fmt(s.mean_i_over_sigma, 2) << " " << Fmt(s.r_meas, 4) << " " << Fmt(s.cc_half, 4); + << Fmt(s.mean_i_over_sigma, 2) << " " << Fmt(s.r_merge, 4) << " " << Fmt(s.r_meas, 4) << " " + << Fmt(s.r_pim, 4) << " " << Fmt(s.cc_half, 4) << " " << Fmt(CCStar(s.cc_half), 4); if (has_anom) out << " " << Fmt(s.abs_diff_over_sigma_anomalous, 3); out << "\n"; @@ -454,6 +557,67 @@ 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(); } @@ -1072,6 +1236,7 @@ void WriteReflections(const std::vector &reflections, const MergeStatistics &statistics, const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, + const OutputProvenance &provenance, const std::string &filename, size_t nthreads, const std::vector *scaled_fulls) { @@ -1081,6 +1246,7 @@ void WriteReflections(const std::vector &reflections, // fulls where the merge provides them (rotation), the merged reflections otherwise (stills). WriteMtzReflections(reflections, unitCell, experiment, filename + ".mtz"); WriteMmcifReflections(reflections, unitCell, experiment, statistics, error_model, twinning, + scaled_fulls ? *scaled_fulls : std::vector{}, 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 697cbe87c..bdac24f65 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -27,13 +27,24 @@ struct ErrorModelReport { std::string b; }; -// nthreads: workers for the per-reflection row formatting, which is the bulk of the file. +// What the reflection files record about where they came from, beyond the experiment itself. +// Empty strings are left out. +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. 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); @@ -92,6 +103,7 @@ void WriteReflections(const std::vector &reflections, const MergeStatistics &statistics, const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, + const OutputProvenance &provenance, const std::string &filename, size_t nthreads, const std::vector *scaled_fulls = nullptr); diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 4c63aa15c..af996e7d6 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -669,15 +669,23 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m } // R_meas per shell: sum over reflections of sqrt(n/(n-1)) * sum_i|I_i - |, over sum of I_i. + // R_merge and R_pim beside it, with the factor 1 and sqrt(1/(n-1)) - R_pim following Weiss (2001) + // J. Appl. Cryst. 34, 130-135. std::vector rmeas_num(n_shells, 0.0), rmeas_den(n_shells, 0.0); - double rmeas_num_all = 0.0, rmeas_den_all = 0.0; + std::vector rmerge_num(n_shells, 0.0), rpim_num(n_shells, 0.0); + double rmeas_num_all = 0.0, rmeas_den_all = 0.0, rmerge_num_all = 0.0, rpim_num_all = 0.0; for (const auto &[key, ra]: rmeas_obs) { if (ra.n < 2 || ra.shell < 0 || ra.shell >= n_shells) continue; const double factor = std::sqrt(static_cast(ra.n) / (ra.n - 1)); + const double pim_factor = std::sqrt(1.0 / (ra.n - 1)); rmeas_num[ra.shell] += factor * ra.sum_abs_dev; + rmerge_num[ra.shell] += ra.sum_abs_dev; + rpim_num[ra.shell] += pim_factor * ra.sum_abs_dev; rmeas_den[ra.shell] += ra.sum_I; rmeas_num_all += factor * ra.sum_abs_dev; + rmerge_num_all += ra.sum_abs_dev; + rpim_num_all += pim_factor * ra.sum_abs_dev; rmeas_den_all += ra.sum_I; } @@ -701,6 +709,8 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m ss.cc_half = sa.cc_half.GetCC(); ss.cc_ref = sa.cc_ref.GetCC(); ss.r_meas = rmeas_den[s] > 0.0 ? rmeas_num[s] / rmeas_den[s] : NAN; + ss.r_merge = rmeas_den[s] > 0.0 ? rmerge_num[s] / rmeas_den[s] : NAN; + ss.r_pim = rmeas_den[s] > 0.0 ? rpim_num[s] / rmeas_den[s] : NAN; } auto &overall = out.overall; @@ -728,6 +738,8 @@ MergeStatistics MergeOnTheFly::MergeStats(const std::vector &m overall.cc_half = cc_half_overall.GetCC(); overall.cc_ref = cc_ref_overall.GetCC(); overall.r_meas = rmeas_den_all > 0.0 ? rmeas_num_all / rmeas_den_all : NAN; + overall.r_merge = rmeas_den_all > 0.0 ? rmerge_num_all / rmeas_den_all : NAN; + overall.r_pim = rmeas_den_all > 0.0 ? rpim_num_all / rmeas_den_all : NAN; return out; } @@ -798,15 +810,18 @@ std::ostream &operator<<(std::ostream &output, const MergeStatisticsShell &in) { double multiplicity = in.unique_reflections > 0 ? static_cast(in.total_observations) / in.unique_reflections : NAN; - output << fmt::format("{:8d} {:8d} {:8d} {:>8s} {:>7s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}", + output << fmt::format("{:8d} {:8d} {:8d} {:>8s} {:>7s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}", in.total_observations, in.unique_reflections, in.possible_unique_reflections, Measured(completeness, "%"), Measured(multiplicity), Measured(in.mean_i_over_sigma), + Measured(in.r_merge*100.0, "%"), Measured(in.r_meas*100.0, "%"), + Measured(in.r_pim*100.0, "%"), Measured(in.cc_half*100.0, "%"), + Measured(CCStar(in.cc_half)*100.0, "%"), Measured(in.cc_ref*100.0, "%"), Measured(in.cc_anom*100.0, "%"), Measured(in.abs_diff_over_sigma_anomalous, "", 2)); @@ -814,14 +829,14 @@ std::ostream &operator<<(std::ostream &output, const MergeStatisticsShell &in) { } std::string ShellTableHeader() { - return fmt::format(" {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>7s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}", - "d_min", "N_obs", "N_uniq", "N_possib", "Compl", "Mult", "", "R_meas", "CC1/2", "CCref", - "CCanom", "SigAno"); + return fmt::format(" {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>7s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s} {:>8s}", + "d_min", "N_obs", "N_uniq", "N_possib", "Compl", "Mult", "", "R_merge", "R_meas", "R_pim", + "CC1/2", "CC*", "CCref", "CCanom", "SigAno"); } std::string ShellTableRule() { - return fmt::format(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->7s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}", - "", "", "", "", "", "", "", "", "", "", "", ""); + return fmt::format(" {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->7s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s} {:->8s}", + "", "", "", "", "", "", "", "", "", "", "", "", "", "", ""); } std::ostream &operator<<(std::ostream &output, const MergeStatistics &in) { diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 6119ab318..9ab944876 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -39,6 +39,12 @@ struct MergeStatisticsShell { // error model (WeightedRmeasTerms); with equal sigmas it is r_meas. A frame kept at a small weight // then no longer sets the number with its noise. NaN where it is not computed (stills, per hand). double r_meas_weighted = NAN; + // The other two merging R-factors, over the same observations and means as r_meas and differing + // from it only in the per-reflection factor on sum_i|I_i-|: 1 for R_merge, sqrt(1/(n-1)) for + // the precision-indicating R_pim (Weiss (2001) J. Appl. Cryst. 34, 130-135). Reported for + // deposition, which asks for all three. + double r_merge = NAN; + double r_pim = NAN; // Anomalous signal-to-noise (XDS "SigAno" / mmCIF pdbx_absDiff_over_sigma_anomalous): // <|I(+)-I(-)|> / over acentric reflections measured in both hands. NaN when @@ -207,6 +213,12 @@ struct MergeStatistics { }; +// CC* = sqrt(2 CC1/2 / (1 + CC1/2)), the correlation of the merged data with the noise-free signal +// (Karplus & Diederichs (2012) Science 336, 1030-1033). NaN where CC1/2 is not positive. +inline double CCStar(double cc_half) { + return cc_half > 0.0 ? std::sqrt(2.0 * cc_half / (1.0 + cc_half)) : NAN; +} + std::ostream &operator<<(std::ostream &output, const MergeStatisticsShell &in); std::ostream &operator<<(std::ostream &output, const MergeStatistics &in); // The header and the rule of the table operator<< prints, for a table with the same columns rendered diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 57e35e070..8e6c9c39b 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -5409,7 +5409,8 @@ 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)}; + 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}; } }); } @@ -5639,6 +5640,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool double rmeas_num_all = 0.0, rmeas_den_all = 0.0; std::vector urmeas_num(n_shells, 0.0), urmeas_den(n_shells, 0.0); double urmeas_num_all = 0.0, urmeas_den_all = 0.0; + std::vector rmerge_num(n_shells, 0.0), rpim_num(n_shells, 0.0); + double rmerge_num_all = 0.0, rpim_num_all = 0.0; // R_meas and the observation counts behind the reported multiplicity. Report-only, and the last // full pass over the observations there is - so, like the correction surfaces, it is not made on a // merge whose numbers nobody reads. @@ -5733,6 +5736,11 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool const double factor = std::sqrt(static_cast(r.n) / (r.n - 1)); urmeas_num[r.shell] += factor * r.sum_abs_dev; urmeas_den[r.shell] += r.sum_I; urmeas_num_all += factor * r.sum_abs_dev; urmeas_den_all += r.sum_I; + // R_merge and R_pim on the same sums, unweighted as r_meas is: factor 1, and + // sqrt(1/(n-1)) following Weiss (2001) J. Appl. Cryst. 34, 130-135. + const double pim_factor = std::sqrt(1.0 / (r.n - 1)); + rmerge_num[r.shell] += r.sum_abs_dev; rmerge_num_all += r.sum_abs_dev; + rpim_num[r.shell] += pim_factor * r.sum_abs_dev; rpim_num_all += pim_factor * r.sum_abs_dev; } } @@ -6003,6 +6011,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool ss.cc_ref = NAN; ss.r_meas = urmeas_den[s] > 0.0 ? urmeas_num[s] / urmeas_den[s] : NAN; ss.r_meas_weighted = rmeas_den[s] > 0.0 ? rmeas_num[s] / rmeas_den[s] : NAN; + ss.r_merge = urmeas_den[s] > 0.0 ? rmerge_num[s] / urmeas_den[s] : NAN; + ss.r_pim = urmeas_den[s] > 0.0 ? rpim_num[s] / urmeas_den[s] : NAN; ss.abs_diff_over_sigma_anomalous = sig_den[s] > 0.0 ? sig_num[s] / sig_den[s] : NAN; ss.cc_anom = cc_anom_shell[s].GetCC(); } @@ -6019,6 +6029,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool overall.cc_ref = NAN; overall.r_meas = urmeas_den_all > 0.0 ? urmeas_num_all / urmeas_den_all : NAN; overall.r_meas_weighted = rmeas_den_all > 0.0 ? rmeas_num_all / rmeas_den_all : NAN; + overall.r_merge = urmeas_den_all > 0.0 ? rmerge_num_all / urmeas_den_all : NAN; + overall.r_pim = urmeas_den_all > 0.0 ? rpim_num_all / urmeas_den_all : NAN; overall.abs_diff_over_sigma_anomalous = sig_den_all > 0.0 ? sig_num_all / sig_den_all : NAN; overall.cc_anom = cc_anom_overall.GetCC(); @@ -6039,6 +6051,8 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool hs.cc_half = ha[sh].cc_half.GetCC(); hs.r_meas = ha[sh].rmeas_den > 0.0 ? ha[sh].rmeas_num / ha[sh].rmeas_den : NAN; hs.r_meas_weighted = NAN; + hs.r_merge = NAN; // not formed per hand + hs.r_pim = NAN; out.per_hand.overall.total_observations += ha[sh].total_obs; out.per_hand.overall.unique_reflections += ha[sh].unique; out.per_hand.overall.possible_unique_reflections += ha[sh].possible; diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index 0444a668f..12ba9a265 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -855,6 +855,7 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen InstrumentMetadata metadata; metadata.InstrumentName(master_file->GetString("/entry/instrument/name")); metadata.SourceName(master_file->GetString("/entry/source/name")); + metadata.SourceType(master_file->GetString("/entry/source/type")); dataset->experiment.ImportInstrumentMetadata(metadata); // The rotation axis is whatever the file calls it. The name is free-form throughout the API, diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 7a752dae1..2ac62f675 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -995,8 +995,14 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyReal("R_MEAS", o.r_meas, "{:.4f}")); if (std::isfinite(o.r_meas_weighted)) Add(s, KeyReal("R_MEAS_WEIGHTED", o.r_meas_weighted, "{:.4f}")); + if (std::isfinite(o.r_merge)) + Add(s, KeyReal("R_MERGE", o.r_merge, "{:.4f}")); + if (std::isfinite(o.r_pim)) + Add(s, KeyReal("R_PIM", o.r_pim, "{:.4f}")); if (std::isfinite(o.cc_half)) Add(s, KeyReal("CC_HALF", o.cc_half, "{:.4f}")); + if (std::isfinite(CCStar(o.cc_half))) + Add(s, KeyReal("CC_STAR", CCStar(o.cc_half), "{:.4f}")); if (std::isfinite(o.abs_diff_over_sigma_anomalous)) Add(s, KeyReal("SIGANO", o.abs_diff_over_sigma_anomalous, "{:.3f}")); if (std::isfinite(o.cc_anom)) diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 4df80b6b8..51b284640 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -40,6 +40,7 @@ #include "../image_analysis/geom_refinement/BeamCenterFromBackground.h" // BeamCenterEstimate #include "../image_analysis/rotation_indexer/RotationIndexer.h" // RotationIndexerResult #include "../image_analysis/geom_refinement/PowderCalibration.h" // CalibrationMethod, CalibrationResult +#include "../image_analysis/WriteReflections.h" // OutputProvenance class JFJochReader; class JFJochHDF5Reader; @@ -70,6 +71,8 @@ struct ProcessConfig { // Output prefix for the _process.h5 (and scaled reflections). Empty => process without writing. std::string output_prefix; + // What the reflection files record about where they came from (WriteReflections). + OutputProvenance provenance; // Write the per-image _process.h5. Defaults true; rugnux turns it off when merging // (the .mtz/.cif is the wanted output and the h5 is large) unless --write-process-h5 is given. diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 823d15df6..e00a0f95a 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -3477,7 +3477,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { f.h = h[0]; f.k = h[1]; f.l = h[2]; } WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics, - em_report, result.twinning, config_.output_prefix, + em_report, result.twinning, config_.provenance, config_.output_prefix, static_cast(std::max(1, config_.nthreads)), &sm.scaled_fulls); writing_mark.Stop(); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 9e74ef787..003f59bde 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1664,6 +1664,9 @@ static int RunRugnux(int argc, char **argv) { } logger.Info("Loaded dataset from {}", input_file); + // Recorded in the reflection files. + OutputProvenance output_provenance; + output_provenance.collection_date = dataset->arm_date; { // The header as the file states it, before anything is measured: what a question about a run // is answered from first. @@ -2229,7 +2232,8 @@ static int RunRugnux(int argc, char **argv) { error_model_a > 0 ? fmt::format("{:.3f}", error_model_a) : std::string(), error_model_b > 0 ? fmt::format("{:.4e}", error_model_b) : std::string()}; WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, merged_statistics, - em_report, twinning, output_prefix, static_cast(nthreads), &scaled_fulls); + em_report, twinning, output_provenance, output_prefix, static_cast(nthreads), + &scaled_fulls); } // --mode scale re-merges stored reflections, so it determines a merging result and gets the @@ -2986,6 +2990,7 @@ static int RunRugnux(int argc, char **argv) { // Rugnux; the experiment above carries all algorithm settings). ProcessConfig config; config.mode = ProcessMode::FullAnalysis; + config.provenance = output_provenance; config.start_image = start_image; config.end_image = end_image; config.stride = image_stride; diff --git a/tests/WriteReflectionsTest.cpp b/tests/WriteReflectionsTest.cpp index 7f3e070ab..9c95b7c07 100644 --- a/tests/WriteReflectionsTest.cpp +++ b/tests/WriteReflectionsTest.cpp @@ -9,10 +9,14 @@ #include #include +#include +#include #include "../common/DiffractionExperiment.h" #include "../image_analysis/WriteReflections.h" #include "../image_analysis/IntegrationOutcome.h" +#include "../image_analysis/scale_merge/Merge.h" +#include "../image_analysis/scale_merge/TwinningAnalysis.h" #include #include "SyntheticMergedReflections.h" @@ -217,3 +221,85 @@ TEST_CASE("SHELX .hkl of unmerged fulls: fixed 3I4,2F8.2, scaled to fit, 0 0 0 t CHECK(lines[2] == " 1 2 3 -2.50 1.00"); CHECK(lines[3] == " 0 0 0 0.00 0.00"); } + +TEST_CASE("CC* from CC1/2", "[write_reflections][portable]") { + CHECK(CCStar(1.0) == Catch::Approx(1.0)); + CHECK(CCStar(0.5) == Catch::Approx(std::sqrt(2.0 / 3.0))); + CHECK(std::isnan(CCStar(0.0))); + CHECK(std::isnan(CCStar(-0.1))); + CHECK(std::isnan(CCStar(NAN))); +} + +TEST_CASE("Merged mmCIF: PDBx keys, merging R-factors and the scaled unmerged block", + "[write_reflections][portable]") { + jfjoch_test::SyntheticMergeParams params; + params.true_space_group = "P 43 21 2"; + params.twin_supergroup = "P 43 21 2"; + params.d_min_A = 5.0; + const auto reflections = jfjoch_test::GenerateSyntheticMerged(params); + REQUIRE(!reflections.empty()); + + auto experiment = TestExperiment(); + experiment.SampleTemperature_K(100.0f); + + MergeStatistics statistics; + MergeStatisticsShell shell; + shell.d_min = 5.0f; shell.d_max = 30.0f; + shell.total_observations = 400; shell.unique_reflections = 100; shell.possible_unique_reflections = 110; + shell.mean_i_over_sigma = 10.0; + shell.cc_half = 0.98; + shell.r_merge = 0.05; shell.r_meas = 0.06; shell.r_pim = 0.03; + 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"; + + 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); + gemmi::cif::Document doc = gemmi::read_cif_gz(path); + std::filesystem::remove(path); + + REQUIRE(doc.blocks.size() == 2); + gemmi::cif::Block &merged = doc.blocks[0]; + CHECK(merged.name == "sample"); + CHECK(*merged.find_value("_reflns.pdbx_ordinal") == "1"); + CHECK(*merged.find_value("_reflns.pdbx_diffrn_id") == "1"); + CHECK(*merged.find_value("_reflns.pdbx_Rmerge_I_obs") == "0.0500"); + CHECK(*merged.find_value("_reflns.pdbx_Rpim_I_all") == "0.0300"); + CHECK(gemmi::cif::as_number(*merged.find_value("_reflns.pdbx_CC_star")) + == Catch::Approx(std::sqrt(2 * 0.98 / 1.98)).margin(1e-4)); + CHECK(*merged.find_value("_diffrn.ambient_temp") == "100.0"); + CHECK(*merged.find_value("_diffrn_detector.pdbx_collection_date") == "2026-07-13"); + + gemmi::cif::Table software = merged.find("_software.", {"pdbx_ordinal", "classification"}); + REQUIRE(software.length() == 2); + CHECK(gemmi::cif::as_string(software[0][1]) == "data reduction"); + CHECK(gemmi::cif::as_string(software[1][1]) == "data scaling"); + + gemmi::cif::Table shells = merged.find("_reflns_shell.", {"pdbx_ordinal", "Rmerge_I_obs", + "pdbx_Rpim_I_all", "pdbx_CC_star"}); + REQUIRE(shells.length() == 1); + 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)); +}