diff --git a/common/Reflection.h b/common/Reflection.h index b5ae95f65..bcb614e40 100644 --- a/common/Reflection.h +++ b/common/Reflection.h @@ -60,6 +60,19 @@ struct Reflection { bool clipped = false; }; +// One full reflection of a rotation sweep as the merge used it - its partials summed into one +// 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. +struct ScaledFull { + int32_t h = 0; + int32_t k = 0; + int32_t l = 0; + float I = NAN; + float sigma = NAN; +}; + struct MergedReflection { int32_t h = 0; int32_t k = 0; diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 9a34391b6..c487269d3 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,7 @@ ### 1.0.0-rc.174 +* Rugnux: `.hkl` now holds unmerged, scaled full reflections (SHELX HKLF 4) on rotation data, so SHELXL computes Rint itself. * Rugnux reads Rigaku d*TREK SMV images (Saturn CCD), including detector 2theta, image orientation and encoded pixel overflows. ### 1.0.0-rc.173 diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index f2e6c5d82..a1a23340f 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -57,7 +57,7 @@ in both. A rotation run that merges — the default — leaves seven files next ``` myrun.mtz merged intensities + French-Wilson amplitudes, for CCP4 / phenix myrun.cif the same, as mmCIF - the self-describing format, and what to deposit -myrun.hkl the same, as SHELX HKLF 4 - feed this to SHELXC / SHELXD / ANODE +myrun.hkl the scaled observations unmerged, as SHELX HKLF 4 - for SHELXL, SHELXC / ANODE myrun_unmerged.mtz every observation before scaling, for pointless / aimless / careless - the largest file of the run (--no-export-unmerged skips it) myrun_P1.mtz the same observations merged in P1, so a wrong space-group call can be @@ -85,7 +85,7 @@ A few things worth knowing before reaching for more flags: (see [Rotation data](RUGNUX_TUTORIAL.md#rotation-data)). `-S` takes a Hermann-Mauguin symbol (`P43212`) or a space-group number (`96`), whichever is to hand. - **Anomalous data are there without `-A`.** A rotation merge always keeps the Bijvoet split: a - default run's `.mtz` carries `I(+)`/`I(-)` beside `IMEAN`, and its `.hkl` the ±hkl mates — + default run's `.mtz` carries `I(+)`/`I(-)` beside `IMEAN`, and its `.hkl` every observation at the index it was measured at — `FRIEDELS_LAW= TRUE` in the report says how the *statistics* were counted, not that the signal was averaged away. What `-A` changes is the counting basis and the error model: each hand becomes a merged observation of its own, so multiplicity, completeness and ⟨I/σ⟩ are counted anomalously and diff --git a/docs/RUGNUX_INTEGRATION.md b/docs/RUGNUX_INTEGRATION.md index c98c769ce..06f3596bf 100644 --- a/docs/RUGNUX_INTEGRATION.md +++ b/docs/RUGNUX_INTEGRATION.md @@ -38,12 +38,17 @@ ignores them, and one that does can find them without guessing. > whole-range value. There is no version marker inside the file, so a number taken from an older > `.cif` is not comparable with one taken from a newer one. -**SHELX HKLF 4** (`.hkl`). Fixed-format `3I4,2F8.2` — `h k l I σ(I)`, one record per -reflection, terminated by a `0 0 0` record — which is what **SHELXC**, **SHELXD** and **ANODE** -expect. Two properties worth knowing before using it: +**SHELX HKLF 4** (`.hkl`). Fixed-format `3I4,2F8.2` — `h k l I σ(I)`, terminated by a +`0 0 0` record — which is what **SHELXL**, **SHELXC**, **SHELXD** and **ANODE** expect. Two +properties worth knowing before using it: -- **Bijvoet mates are written separately**, `I(+)` at `+hkl` and `I(-)` at `-hkl`, so the anomalous - differences survive into SHELXC; a reflection with no anomalous split is written once, as its mean. +- **On rotation data it is unmerged**: one record per full reflection (its partials summed), with the + per-frame scale and every correction applied and the σ(I) the merge weighted it with, but not + averaged with its symmetry equivalents, and at the index it was measured at — the chemical + crystallographer's convention, so SHELXL computes Rint and Rsigma itself and Friedel mates keep + their own records. Observations the merge rejected as outliers, or that lie beyond its resolution + cut, are left out. There is no batch number column (in HKLF 4 that selects a BASF scale factor). + A **stills** run writes the merged reflections instead, Bijvoet mates at `+hkl` and `-hkl`. - **Intensities are rescaled** by a single global factor so the largest value fits the `F8.2` field. `I` and `σ(I)` share that factor, so every ratio — and therefore the anomalous signal — is untouched, but the absolute scale is not meaningful. This matters only if you intend to compare @@ -371,12 +376,12 @@ should have produced rather than trust the exit status. **Nothing has to be switched on to get the anomalous signal.** A default rotation merge keeps the Bijvoet split, whether or not `-A` was given: `myrun.mtz` carries `I(+)`/`I(-)` and `F(+)`/`F(-)` -beside the means, and `myrun.hkl` writes each mate as its own record, `I(+)` at `+hkl` and `I(-)` at -`-hkl`. `-A` changes what the merging statistics are counted over, not whether the signal is in the -file. The one case with no anomalous columns at all is a **stills** run, which computes no Bijvoet -split; there `myrun.hkl` holds means only and there is nothing for SHELXC to work with. Unmerged -data are not wanted anywhere in this chain either, so a run with `--no-export-unmerged` is not -missing a file SHELX needs. +beside the means, and `myrun.hkl` holds every observation unmerged at the index it was measured at, so +SHELXC sees both hands. `-A` changes what the merging statistics are counted over, not whether the +signal is in the file. The one case with no anomalous signal at all is a **stills** run, which +computes no Bijvoet split; there `myrun.hkl` holds means only and there is nothing for SHELXC to work +with. `myrun_unmerged.mtz` is not part of this chain (it is unscaled), so a run with +`--no-export-unmerged` is not missing a file SHELX needs. **HKLF 4 carries no metadata**, so the cell and the space group have to be repeated on the SHELXC command — take them from `UNIT_CELL_CONSTANTS` and `SPACE_GROUP_NAME` in section 2 of the diff --git a/docs/RUGNUX_TUTORIAL.md b/docs/RUGNUX_TUTORIAL.md index 90972a914..f7664c38e 100644 --- a/docs/RUGNUX_TUTORIAL.md +++ b/docs/RUGNUX_TUTORIAL.md @@ -205,10 +205,12 @@ goniometer axis but you want per-frame stills processing anyway, add `--force-st - `.cif` — mmCIF, for deposition and as the self-describing native format (also carries the merging statistics, ISa, twinning and radiation-damage indicators). - `.hkl` — SHELX **HKLF 4** text (`h k l I σ(I)`, fixed `3I4,2F8.2`), the direct input for - **SHELXC / ANODE / SHELXD**. Bijvoet mates are written separately (`I(+)` at `+hkl`, `I(-)` at - `-hkl`) so the anomalous signal is preserved; intensities are put on a common scale so the largest - value fits the fixed-width field (the absolute scale is irrelevant to SHELXC/ANODE), and the file - ends with the `0 0 0` terminator record. + **SHELXL** and **SHELXC / ANODE / SHELXD**. On rotation data it holds the scaled full reflections + **unmerged** — every correction applied, symmetry equivalents not averaged, each at the index it was + measured at — so SHELXL reports Rint and Rsigma itself and the anomalous signal is all there; + stills runs write the merged reflections (Bijvoet mates at `+hkl` and `-hkl`). Intensities are put + on a common scale so the largest value fits the fixed-width field, and the file ends with the + `0 0 0` terminator record. All three carry the **refined unit cell** (from rotation indexing) and the **space group determined from systematic absences** (constrained to the indexed lattice symmetry). diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 798467ec4..281cf5629 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -532,55 +532,35 @@ void WriteMtzReflections(const std::vector &reflections, mtz.write_to_file(filename); } -void WriteShelxHklReflections(const std::vector &reflections, - const DiffractionExperiment &experiment, - const std::string &filename, - size_t nthreads) { - bool has_anom = true; - const std::vector rows = BuildMergedRows(reflections, experiment, has_anom); - - // SHELX HKLF 4 (SHELXC / ANODE input): fixed FORMAT(3I4,2F8.2), one record per reflection as - // h k l I sigma(I). The Bijvoet mates are written separately - I(+) at +hkl, I(-) at -hkl - so the - // anomalous differences survive; a reflection with no anomalous split is written once as its mean. - // Intensities are put on a common scale so the largest value fits the F8.2 field (the absolute scale - // is irrelevant to SHELXC / ANODE, which use only ratios); I and sigma share the scale, so the - // anomalous signal is untouched. The file ends with a 0 0 0 terminator record. - const auto usable = [](float v, float s) { return std::isfinite(v) && std::isfinite(s) && s > 0.0f; }; - - double max_abs = 0.0; - for (const auto& r : rows) { - if (usable(r.Ip, r.sIp)) max_abs = std::max({max_abs, std::fabs(double(r.Ip)), double(r.sIp)}); - if (usable(r.Im, r.sIm)) max_abs = std::max({max_abs, std::fabs(double(r.Im)), double(r.sIm)}); - if (!usable(r.Ip, r.sIp) && !usable(r.Im, r.sIm) && usable(r.Imean, r.sImean)) - max_abs = std::max({max_abs, std::fabs(double(r.Imean)), double(r.sImean)}); - } - const double scale = (std::isfinite(max_abs) && max_abs > 0.0) ? 9999.0 / max_abs : 1.0; - - std::ofstream out(filename); - if (!out) - throw std::runtime_error("WriteShelxHklReflections: cannot open " + filename); - // Up to two records per reflection, five formatted numbers each. Built in parallel into per-worker - // blocks and handed to the file in order, exactly as the mmCIF rows are; "%.2f" right-aligned in - // the fixed field is what `fixed` + `setprecision(2)` + `setw` made the stream write. - const auto column = [](std::string &s, const std::string &v, size_t width) { - if (v.size() < width) s.append(width - v.size(), ' '); - s.append(v); - }; - const auto num2 = [](double v) { - char b[64]; - const int n = std::snprintf(b, sizeof b, "%.2f", v); - return std::string(b, static_cast(std::clamp(n, 0, static_cast(sizeof b) - 1))); - }; - const auto emit = [&column, &num2, scale](std::string &s, int h, int k, int l, float I, float sigma) { - column(s, std::to_string(h), 4); - column(s, std::to_string(k), 4); - column(s, std::to_string(l), 4); - column(s, num2(scale * I), 8); - column(s, num2(scale * sigma), 8); +namespace { + // The SHELX HKLF 4 record, FORMAT(3I4,2F8.2): h k l I sigma(I). "%.2f" right-aligned in the fixed + // field is what `fixed` + `setprecision(2)` + `setw` made the stream write. + void AppendHklf4(std::string &s, int h, int k, int l, double I, double sigma) { + const auto column = [&s](const std::string &v, size_t width) { + if (v.size() < width) s.append(width - v.size(), ' '); + s.append(v); + }; + const auto num2 = [](double v) { + char b[64]; + const int n = std::snprintf(b, sizeof b, "%.2f", v); + return std::string(b, static_cast(std::clamp(n, 0, static_cast(sizeof b) - 1))); + }; + column(std::to_string(h), 4); + column(std::to_string(k), 4); + column(std::to_string(l), 4); + column(num2(I), 8); + column(num2(sigma), 8); s.push_back('\n'); - }; - { - const size_t nrow = rows.size(); + } + + // An HKLF 4 file of nrow source rows, each appending its records with append_row(i, s). Built in + // parallel into per-worker blocks and handed to the file in order, exactly as the mmCIF rows are, + // and closed with the 0 0 0 end-of-data record. + template + void WriteHklf4File(const std::string &filename, size_t nrow, size_t nthreads, AppendRow append_row) { + std::ofstream out(filename); + if (!out) + throw std::runtime_error("WriteShelxHklReflections: cannot open " + filename); const size_t nw = std::max(nthreads, 1); const int nch = static_cast(ThreadsForWork(nrow, nw, 4096)); std::vector block(nch); @@ -589,24 +569,71 @@ void WriteShelxHklReflections(const std::vector &reflections, const size_t lo = nrow * t / nch, hi = nrow * (t + 1) / nch; std::string &s = block[t]; s.reserve((hi - lo) * 2 * 29); - for (size_t i = lo; i < hi; ++i) { - const auto &r = rows[i]; - const bool plus = usable(r.Ip, r.sIp); - const bool minus = usable(r.Im, r.sIm); - if (plus) emit(s, r.h, r.k, r.l, r.Ip, r.sIp); - if (minus) emit(s, -r.h, -r.k, -r.l, r.Im, r.sIm); - if (!plus && !minus && usable(r.Imean, r.sImean)) - emit(s, r.h, r.k, r.l, r.Imean, r.sImean); - } + for (size_t i = lo; i < hi; ++i) + append_row(i, s); } }); for (const std::string &s : block) out.write(s.data(), static_cast(s.size())); + std::string tail; + AppendHklf4(tail, 0, 0, 0, 0.0, 0.0); + out.write(tail.data(), static_cast(tail.size())); } - std::string tail; - emit(tail, 0, 0, 0, 0.0f, 0.0f); // HKLF-4 end-of-data marker - out.write(tail.data(), static_cast(tail.size())); - out.close(); + + // Intensities are put on a common scale so the largest value fits the F8.2 field (the absolute scale + // is irrelevant to every SHELX program, which refines or ignores it); I and sigma share the scale. + double Hklf4Scale(double max_abs) { + return (std::isfinite(max_abs) && max_abs > 0.0) ? 9999.0 / max_abs : 1.0; + } + + bool UsableForHkl(float v, float s) { return std::isfinite(v) && std::isfinite(s) && s > 0.0f; } +} + +void WriteShelxHklReflections(const std::vector &reflections, + const DiffractionExperiment &experiment, + const std::string &filename, + size_t nthreads) { + bool has_anom = true; + const std::vector rows = BuildMergedRows(reflections, experiment, has_anom); + + // One record per merged reflection. The Bijvoet mates are written separately - I(+) at +hkl, I(-) + // at -hkl - so the anomalous differences survive; a reflection with no anomalous split is written + // once as its mean. + double max_abs = 0.0; + for (const auto& r : rows) { + if (UsableForHkl(r.Ip, r.sIp)) max_abs = std::max({max_abs, std::fabs(double(r.Ip)), double(r.sIp)}); + if (UsableForHkl(r.Im, r.sIm)) max_abs = std::max({max_abs, std::fabs(double(r.Im)), double(r.sIm)}); + if (!UsableForHkl(r.Ip, r.sIp) && !UsableForHkl(r.Im, r.sIm) && UsableForHkl(r.Imean, r.sImean)) + max_abs = std::max({max_abs, std::fabs(double(r.Imean)), double(r.sImean)}); + } + const double scale = Hklf4Scale(max_abs); + WriteHklf4File(filename, rows.size(), nthreads, [&](size_t i, std::string &s) { + const auto &r = rows[i]; + const bool plus = UsableForHkl(r.Ip, r.sIp); + const bool minus = UsableForHkl(r.Im, r.sIm); + if (plus) AppendHklf4(s, r.h, r.k, r.l, scale * r.Ip, scale * r.sIp); + if (minus) AppendHklf4(s, -r.h, -r.k, -r.l, scale * r.Im, scale * r.sIm); + if (!plus && !minus && UsableForHkl(r.Imean, r.sImean)) + AppendHklf4(s, r.h, r.k, r.l, scale * r.Imean, scale * r.sImean); + }); +} + +void WriteShelxHklReflections(const std::vector &fulls, + const std::string &filename, + size_t nthreads) { + // One record per full, at the index it was measured at: SHELXL / SHELXC average the equivalents + // themselves and report Rint and Rsigma from them, which a merged file hides. No batch number - in + // HKLF 4 that column selects a BASF scale factor in SHELXL, which these fulls do not want. + double max_abs = 0.0; + for (const auto &f : fulls) + if (UsableForHkl(f.I, f.sigma)) + max_abs = std::max({max_abs, std::fabs(double(f.I)), double(f.sigma)}); + const double scale = Hklf4Scale(max_abs); + WriteHklf4File(filename, fulls.size(), nthreads, [&](size_t i, std::string &s) { + const auto &f = fulls[i]; + if (UsableForHkl(f.I, f.sigma)) + AppendHklf4(s, f.h, f.k, f.l, scale * f.I, scale * f.sigma); + }); } namespace { @@ -1043,12 +1070,17 @@ void WriteReflections(const std::vector &reflections, const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, const std::string &filename, - size_t nthreads) { + size_t nthreads, + const std::vector *scaled_fulls) { // Write an MTZ, an mmCIF and a SHELX HKLF-4 .hkl - each has its uses downstream (MTZ for the CCP4 / - // phenix reflection tools, mmCIF for deposition and as the self-describing native format, HKLF-4 as - // the SHELXC / ANODE substructure-solution input). + // phenix reflection tools, mmCIF for deposition and as the self-describing native format, HKLF-4 for + // SHELXL refinement and SHELXC / ANODE substructure solution). The .hkl holds the unmerged scaled + // 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, filename + ".cif", nthreads); - WriteShelxHklReflections(reflections, experiment, filename + ".hkl", nthreads); + if (scaled_fulls && !scaled_fulls->empty()) + WriteShelxHklReflections(*scaled_fulls, filename + ".hkl", nthreads); + else + WriteShelxHklReflections(reflections, experiment, filename + ".hkl", nthreads); } diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 9d1ff5421..697cbe87c 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -42,12 +42,17 @@ void WriteMtzReflections(const std::vector &reflections, const DiffractionExperiment &experiment, const std::string &filename); -// SHELX HKLF-4 text file (h k l I sigma(I), Bijvoet mates separate) for SHELXC / ANODE. -// nthreads: workers for the per-reflection row formatting, as for the mmCIF. +// SHELX HKLF-4 text file, FORMAT(3I4,2F8.2) h k l I sigma(I), scaled so the largest value fits the +// field, ended by a 0 0 0 record. nthreads: workers for the row formatting, as for the mmCIF. +// Merged reflections, Bijvoet mates as separate records (the stills output). void WriteShelxHklReflections(const std::vector &reflections, const DiffractionExperiment &experiment, const std::string &filename, size_t nthreads); +// Unmerged scaled fulls, one record each at its measured index (the rotation output). +void WriteShelxHklReflections(const std::vector &fulls, + const std::string &filename, + size_t nthreads); // Unmerged observations in the column and batch-header layout POINTLESS writes: aimless, pointless, // careless and iotbx.merging_statistics all read that layout. H K L are the ASU indices and M/ISYM @@ -88,4 +93,5 @@ void WriteReflections(const std::vector &reflections, const ErrorModelReport &error_model, const TwinningAnalysisResult &twinning, const std::string &filename, - size_t nthreads); \ No newline at end of file + size_t nthreads, + const std::vector *scaled_fulls = nullptr); diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 81dabcd22..12f4dba1d 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -4885,6 +4885,23 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool }); } + // The same fulls unmerged, for the SHELX file: those the final merge above accumulated - usable, + // not rejected by either outlier test, inside the cut (judged on the group's d, as the merged + // reflections are) - each with the sigma the merge weighted it by. + if (export_scaled_fulls && full_stats && !for_search) { + result.scaled_fulls.clear(); + for (int i = 0; i < n_full; ++i) { + const int g = mf.group[i]; + if (g < 0 || rejected_obs[i] || !std::isfinite(merged_I[g])) + continue; + if (effective_d_min && acc[g].d < *effective_d_min) + continue; + const Obs &o = fulls[i]; + const float I_corr = o.I * o.corr; + result.scaled_fulls.push_back({o.h, o.k, o.l, I_corr, corrected_sigma(o, I_corr, o.sigma * o.corr)}); + } + } + // Asymptotic I/sigma. ISa is by definition the I -> infinity limit of the signal-to-noise, i.e. the // reproducibility of the strongest reflections (Diederichs, Acta Cryst. D66 (2010), 733-740). The // (a, b) fit above spans the whole intensity range, and a mild excess of scatter at intermediate diff --git a/image_analysis/scale_merge/RotationScaleMerge.h b/image_analysis/scale_merge/RotationScaleMerge.h index 6b5a81216..45ce1f6fa 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.h +++ b/image_analysis/scale_merge/RotationScaleMerge.h @@ -84,6 +84,9 @@ public: int scaling_iterations_partials = 0; int scaling_iterations_fulls = 0; bool scaling_converged = true; + // The fulls the merge kept, unmerged (see ScaledFull): what passed the outlier tests, inside the + // resolution cut. Filled only when SetExportScaledFulls asked for it, and never on a search merge. + std::vector scaled_fulls; }; // experiment: read live (its space group is changed by the caller between Run() calls). @@ -124,6 +127,10 @@ public: // unmerged MTZ describe the merge that was written. void SetWriteBackPerFrameScale(bool on) { write_back_per_frame_scale = on; } + // Whether Run() also hands back the merge's fulls unmerged (Result::scaled_fulls). Off by default: + // it is a copy as long as the fulls, wanted only by a caller that writes them. + void SetExportScaledFulls(bool on) { export_scaled_fulls = on; } + // Override the high-resolution cut for the next Run() - used to gate the de-novo P1 search pass at // >= 1 without cutting the final in-symmetry merge. Reset to the manual limit afterwards. void SetDMinLimit(std::optional d_min_A) { d_min_limit = d_min_A; } @@ -419,6 +426,7 @@ private: bool gpu_active_ = false; #endif bool write_back_per_frame_scale = true; // see SetWriteBackPerFrameScale + bool export_scaled_fulls = false; // see SetExportScaledFulls // --- helpers (each a flat pass; see the .cpp) --- // Turn the per-frame mean background under the reflections (accumulated by the ingest fill loop) into diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index aed1270d0..974404ae1 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -6449,6 +6449,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // The reference-range table and its ISa (--report-resolution); see ProcessResult. std::optional reference_statistics; double reference_isa = 0.0; + // The merge's fulls unmerged, for the .hkl (rotation, written merges only - see below). + std::vector scaled_fulls; }; // The reference path computes each image's G once (per-image scaling against the // reference); the scaling loop below is skipped, so G is stable across the two passes. @@ -6830,6 +6832,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // merge of its own - only the first merge, which the two-pass quality guard reads. // already_run: the engine's result for this merge, when it was made earlier (see the P1 // cross-check); it is then only reported, exactly as if it had been made here. + bool export_scaled_fulls = true; // off around the P1 cross-check, see scale_and_merge auto scale_and_merge = [&](const std::string &label, bool for_search, bool measure_cc = false, std::optional already_run = std::nullopt) @@ -6839,6 +6842,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // The geometry pre-pass merges only to choose the space group and to give the quality // guard something to judge the second pass against; its reflections are never written, // so the parts of the merge that only fill in an output file are skipped there. + // Only a merge that may be written hands its fulls back: they are a copy as long as the + // fulls, and the search merges, the pre-pass and the P1 cross-check never write a .hkl. + rsm->SetExportScaledFulls(!for_search && !geometry_prepass && write_files && config_.write_merged + && export_scaled_fulls); auto r = already_run ? std::move(*already_run) : rsm->Run(for_search, /*full_stats=*/!geometry_prepass, /*measure_cc_before_corrections=*/measure_cc); @@ -6855,7 +6862,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.search_scaling_converged = r.scaling_converged; return ScaleMergeResult{std::move(r.merged), std::move(r.statistics), r.cc_half_before_corrections, - std::move(r.reference_statistics), r.reference_isa}; + std::move(r.reference_statistics), r.reference_isa, + std::move(r.scaled_fulls)}; } // Stills (rotation goes through RotationScaleMerge above): self-scale each image against the // running merge with ScaleOnTheFly (fixed partiality), then merge directly. This runs even @@ -9022,9 +9030,16 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b ? fmt::format("{:.2f}", result.error_model_isa_asymptotic) : std::string(), result.error_model_a > 0 ? fmt::format("{:.3f}", result.error_model_a) : std::string(), result.error_model_b > 0 ? fmt::format("{:.4e}", result.error_model_b) : std::string()}; + // The fulls carry the index they were measured at in the scaler's frame; every relabelling + // the merged reflections went through since (merge_to_written) takes them to the written one. + if (!(merge_to_written == gemmi::Op::identity())) + for (auto &f : sm.scaled_fulls) { + const gemmi::Op::Miller h = merge_to_written.apply_to_hkl({{f.h, f.k, f.l}}); + 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, - static_cast(std::max(1, config_.nthreads))); + static_cast(std::max(1, config_.nthreads)), &sm.scaled_fulls); // P1 cross-check dataset. The group the files above are written in was chosen by the // search, and if that choice is wrong nothing in them says so - every statistic was @@ -9114,7 +9129,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b experiment_.SpaceGroupNumber(1); if (!p1_merged_early) rsm->SetWriteBackPerFrameScale(false); + export_scaled_fulls = false; auto p1 = scale_and_merge("P1 cross-check", false, false, std::move(p1_merged_early)); + export_scaled_fulls = true; rsm->SetWriteBackPerFrameScale(true); // The scaler still holds the observations in the indexing they were merged in, so every // relabelling since (the written setting, the model's indexing) is applied to this merge diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 4eb457f6b..3a8364894 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -1781,6 +1781,7 @@ static int RunRugnux(int argc, char **argv) { const auto scale_start = std::chrono::steady_clock::now(); std::vector merged_reflections; + std::vector scaled_fulls; // the .hkl: rotation merges only MergeStatistics merged_statistics; double error_model_isa = 0.0; double error_model_isa_asymptotic = 0.0; @@ -1806,11 +1807,13 @@ static int RunRugnux(int argc, char **argv) { RotationScaleMerge rsm(experiment, reflections, experiment.GetUnitCell(), scaling_iter, nthreads, logger); rsm.Ingest(); + rsm.SetExportScaledFulls(!output_prefix.empty()); // No pass to compare against here - --mode scale merges the stored reflections once - so // decline the extra merge the pre-correction CC1/2 would cost. auto r = rsm.Run(false, /*full_stats=*/true, /*measure_cc_before_corrections=*/false); merged_reflections = std::move(r.merged); merged_statistics = std::move(r.statistics); + scaled_fulls = std::move(r.scaled_fulls); error_model_isa = r.isa; error_model_isa_asymptotic = r.isa_asymptotic; error_model_isa_resolved = r.isa_resolved; @@ -2030,7 +2033,7 @@ 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)); + em_report, twinning, output_prefix, static_cast(nthreads), &scaled_fulls); } // --mode scale re-merges stored reflections, so it determines a merging result and gets the diff --git a/tests/WriteReflectionsTest.cpp b/tests/WriteReflectionsTest.cpp index 24a532447..83ed12ab7 100644 --- a/tests/WriteReflectionsTest.cpp +++ b/tests/WriteReflectionsTest.cpp @@ -4,6 +4,7 @@ #include #include +#include #include #include @@ -185,3 +186,31 @@ TEST_CASE("Unmerged MTZ: built on several workers, the same rows as on one", "[w CHECK(parallel.batches.back().number == serial.batches.back().number); } } + +TEST_CASE("SHELX .hkl of unmerged fulls: fixed 3I4,2F8.2, scaled to fit, 0 0 0 terminator", + "[write_reflections][portable]") { + // Two fulls of the same reflection and a Friedel mate: written as measured, nothing averaged. The + // largest |I| or sigma is put at 9999.00; an unusable full (sigma not positive) is left out. + const std::vector fulls{ + {1, 2, 3, 200000.0f, 2000.0f}, + {-1, -2, -3, 100000.0f, 1200.0f}, + {1, 2, 3, -50.0f, 20.0f}, + {4, 5, 6, 10.0f, 0.0f}, + }; + const auto path = (std::filesystem::temp_directory_path() / "jfjoch_scaled_fulls_test.hkl").string(); + WriteShelxHklReflections(fulls, path, 2); + + std::ifstream in(path); + std::vector lines; + for (std::string line; std::getline(in, line);) + lines.push_back(line); + std::filesystem::remove(path); + + REQUIRE(lines.size() == 4); + for (const auto &line : lines) + CHECK(line.size() == 28); + CHECK(lines[0] == " 1 2 3 9999.00 99.99"); + CHECK(lines[1] == " -1 -2 -3 4999.50 59.99"); + CHECK(lines[2] == " 1 2 3 -2.50 1.00"); + CHECK(lines[3] == " 0 0 0 0.00 0.00"); +}