diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index e59cf939..6f54e0d6 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -144,6 +144,10 @@ void WriteMmcifReflections(const std::vector &reflections, out << "_reflns.pdbx_Rrim_I_all " << Fmt(ov.r_meas, 4) << "\n"; out << "_reflns.pdbx_CC_half " << Fmt(ov.cc_half, 4) << "\n"; out << "_reflns.jfjoch_diffrn_ISa " << CifStr(isa) << " # asymptotic I/sigma (Diederichs)\n"; + // Dataset-wide isotropic Wilson B-factor estimate (standard PDBx item), analogous to XDS's + // "WILSON LINE ... B=". Emitted only when the log-linear fit succeeded. + if (std::isfinite(statistics.wilson_b) && statistics.wilson_b > 0.0) + out << "_reflns.B_iso_Wilson_estimate " << Fmt(statistics.wilson_b, 2) << "\n"; // Twinning indicators (no standard mmCIF item; same jfjoch local prefix as ISa above). if (twinning.l_test_pairs > 0) { out << "_reflns.jfjoch_L_test_mean_abs_L " << Fmt(twinning.mean_abs_l, 3) diff --git a/image_analysis/bragg_integration/CalcISigma.cpp b/image_analysis/bragg_integration/CalcISigma.cpp index 144ae3f7..8d0875af 100644 --- a/image_analysis/bragg_integration/CalcISigma.cpp +++ b/image_analysis/bragg_integration/CalcISigma.cpp @@ -1,10 +1,21 @@ // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include +#include +#include + #include "CalcISigma.h" #include "Regression.h" #include "../../common/ResolutionShells.h" +// Upper bound on a physically plausible macromolecular isotropic B-factor (A^2). A fit that lands +// outside (0, WILSON_B_MAX) comes from a bad frame / bad dataset (too few or mis-indexed reflections, +// a degenerate resolution range) rather than real Debye-Waller falloff, so it is rejected as +// indeterminate rather than reported. Radiation-damaged / low-resolution data rarely exceeds ~120 A^2; +// 200 leaves generous head-room while still excluding the hundreds-of-A^2 garbage. +static constexpr float WILSON_B_MAX = 200.0f; + void CalcISigma(DataMessage &msg) { CalcISigma(msg, msg.reflections); } @@ -81,11 +92,80 @@ void CalcWilsonBFactor(DataMessage &msg, if (replace_b && valid_shells > 2) { auto reg_result = regression(shells_mean_one_over_d_square, log_I_mean, valid_shells); + const float b_est = -2.0f * reg_result.slope; - if (reg_result.r_square > 0.3) - msg.b_factor = -2.0f * reg_result.slope; + // Accept only a well-correlated, physically plausible fit. Occasional bad frames (an indexing + // glitch, too few reflections) produce a wildly steep Wilson line and a B of several hundred + // A^2 that pollutes the per-image plot; leaving b_factor unset (rendered as NaN) is better than + // emitting garbage. + if (reg_result.r_square > 0.3 && std::isfinite(b_est) && b_est > 0.0f && b_est < WILSON_B_MAX) + msg.b_factor = b_est; } msg.integration_B_logI = log_I_mean; msg.integration_B_one_over_d_square = shells_mean_one_over_d_square; } + +GlobalWilsonB CalcGlobalWilsonB(const std::vector &merged) { + GlobalWilsonB out; + // A dataset-wide estimate needs enough reflections to average the shell means; below this the + // per-image estimate is the only thing on offer and a global number would be meaningless. + if (merged.size() < 100) + return out; + + float d_min = std::numeric_limits::infinity(), d_max = 0.0f; + for (const auto &r : merged) { + if (std::isfinite(r.I) && r.d > 0.0f) { + d_min = std::min(d_min, r.d); + d_max = std::max(d_max, r.d); + } + } + if (!(d_min < d_max)) + return out; + + // Below ~4 A the Wilson plot is non-linear (bonding/solvent structure), so when the data extend to + // lower resolution than that, restrict the fit to d <= 4 A - the standard Wilson-B convention. If + // the whole dataset is coarser than 4 A, fall back to using all of it. + constexpr double WILSON_LOW_RES_LIMIT_A = 4.0; + const float d_low = (d_max > WILSON_LOW_RES_LIMIT_A && d_min < WILSON_LOW_RES_LIMIT_A) + ? static_cast(WILSON_LOW_RES_LIMIT_A) : d_max; + + const int nshells = 20; + ResolutionShells shells(d_min, d_low, nshells); + std::vector I_sum(nshells, 0.0), sig_sum(nshells, 0.0), count(nshells, 0.0); + for (const auto &r : merged) { + if (!std::isfinite(r.I) || r.d <= 0.0f) + continue; + auto s = shells.GetShell(r.d); // reflections coarser than d_low fall outside -> skipped + if (s) { + I_sum[*s] += r.I; + if (std::isfinite(r.sigma) && r.sigma > 0.0f) + sig_sum[*s] += r.sigma; + ++count[*s]; + } + } + + const auto s2 = shells.GetShellMeanOneOverResSq(); + std::vector x, y; + for (int i = 0; i < nshells; ++i) { + // Skip empty / net-negative shells, and shells past the signal limit (mean I/sigma < 1). The + // latter keeps the fit out of the noise floor: without a resolution cut the weakest high-angle + // shells are background-residual-dominated and flatten the Wilson line, deflating B. This mirrors + // XDS/ctruncate fitting only over the meaningful range and makes the estimate insensitive to how + // far the merged data were carried. + if (count[i] > 0 && I_sum[i] > 0.0 && I_sum[i] > sig_sum[i]) { + x.push_back(s2[i]); + y.push_back(std::log(static_cast(I_sum[i] / count[i]))); + } + } + if (x.size() < 3) + return out; + + const auto reg = regression(x, y, x.size()); + const double b = -2.0 * reg.slope; // ~ exp(-2 B s^2), s^2 = 1/(4 d^2), x = 1/d^2 + out.n_shells = static_cast(x.size()); + out.correlation = std::sqrt(std::clamp(static_cast(reg.r_square), 0.0, 1.0)); + if (std::isfinite(b) && b > 0.0 && b < WILSON_B_MAX) + out.b = b; + return out; +} diff --git a/image_analysis/bragg_integration/CalcISigma.h b/image_analysis/bragg_integration/CalcISigma.h index 21102ff8..6251361c 100644 --- a/image_analysis/bragg_integration/CalcISigma.h +++ b/image_analysis/bragg_integration/CalcISigma.h @@ -3,6 +3,9 @@ #pragma once +#include +#include + #include "../../common/JFJochMessages.h" #include "../../common/Reflection.h" @@ -12,3 +15,14 @@ void CalcWilsonBFactor(DataMessage &msg, bool replace_b = true); void CalcISigma(DataMessage &msg, const std::vector &reflections); void CalcWilsonBFactor(DataMessage &msg, const std::vector &reflections, bool replace_b = true); +// Dataset-wide isotropic Wilson B-factor estimate from a log-linear fit of the shell-averaged merged +// intensity against 1/d^2 (B = -2*slope). Robust because it averages over the whole dataset - the +// per-image CalcWilsonBFactor above is a per-frame estimate and much noisier. Diagnostic only (like +// XDS's "WILSON LINE ... B="), not used in scaling. Returns b = NaN when it cannot be determined. +struct GlobalWilsonB { + double b = std::numeric_limits::quiet_NaN(); // Wilson B-factor estimate (A^2) + double correlation = std::numeric_limits::quiet_NaN(); // |correlation| of the log-linear fit + int n_shells = 0; // resolution shells actually used +}; +GlobalWilsonB CalcGlobalWilsonB(const std::vector &merged); + diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 153c1f0f..7a196d05 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -920,6 +920,9 @@ std::ostream &operator<<(std::ostream &output, const MergeStatistics &in) { output << fmt::format(" {:>8s} ", "Overall"); output << in.overall; output << std::endl; + if (std::isfinite(in.wilson_b) && in.wilson_b > 0.0) + output << fmt::format(" Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f})", + in.wilson_b, in.wilson_b_correlation) << std::endl; output << std::endl; return output; } diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index bc468980..e0968e37 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -51,6 +51,12 @@ struct MergeStatisticsShell { struct MergeStatistics { std::vector shells; MergeStatisticsShell overall; + + // Dataset-wide isotropic Wilson B-factor estimate (A^2) from the log-linear fit of the shell-mean + // merged intensity against 1/d^2 (CalcGlobalWilsonB) - the analogue of XDS's "WILSON LINE ... B=". + // Diagnostic only; not used in scaling. NaN when not determined. + double wilson_b = NAN; + double wilson_b_correlation = NAN; }; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index d53331f0..e2197a3e 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -40,6 +40,7 @@ #include "../image_analysis/scale_merge/TwinningAnalysis.h" #include "../image_analysis/scale_merge/HKLKey.h" #include "../image_analysis/WriteReflections.h" +#include "../image_analysis/bragg_integration/CalcISigma.h" #include "../common/Definitions.h" #include "../common/CorrelationCoefficient.h" #include @@ -1110,6 +1111,17 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { } } + // Dataset-wide Wilson B-factor estimate (like XDS's WILSON LINE B). Diagnostic only - it is not + // fed back into scaling; it just lands in the printed statistics, the mmCIF, and the log. + { + const GlobalWilsonB wilson = CalcGlobalWilsonB(sm.merged); + sm.statistics.wilson_b = wilson.b; + sm.statistics.wilson_b_correlation = wilson.correlation; + if (std::isfinite(wilson.b) && wilson.b > 0.0) + logger.Info("Wilson B-factor estimate: {:.2f} A^2 (correlation {:.3f}, {} shells)", + wilson.b, wilson.correlation, wilson.n_shells); + } + stats_text << sm.statistics; result.merge_statistics_text = stats_text.str(); result.has_merge_statistics = true;