From cb08f63a522c404a34457cbe52b7555ef07b100c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 7 Sep 2026 10:34:21 +0200 Subject: [PATCH] model validation: the bulk solvent is searched inside its physical range, not fitted without bounds gemmi offers two scalers and we were using the one without bounds. Its Levenberg-Marquardt path has nothing stopping the flat-solvent parameters from leaving the range the model means anything in; the alternative path that does declare bounds is behind a compile guard we have never enabled. On this corpus six datasets in fifty-one fitted a b_sol outside it, the worst at 1707 A^2. What that does is subtler than a bad scale, and worth recording because it is why nobody noticed: a b_sol that large does not corrupt the solvent term, it switches it off - 1.4% of it survives at 10 A - so the model is simply scaled without a solvent contribution and the R-factors look unremarkable. k_sol and b_sol now come from a grid search over the physical box, with the scale and the anisotropic B refitted at each candidate pair, following Afonine et al. Refitting at each point is what makes it work: clamping the parameters after an unbounded fit costs up to 0.044 in R-free, because it leaves the scale and B where the rejected fit put them. Non-physical fits go from six in fifty-one to none, and both R-work and R-free come out slightly but significantly better rather than merely no worse. Which reflections are fitted remains the caller's business - the function scales whatever it is handed - so the working-set restriction of the previous commit is not something this can undo. A crystal with no solvent-accessible volume needs no special case: its mask is empty, so the solvent term is identically zero whatever the parameters say. There is a test for that. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N --- docs/ACKNOWLEDGEMENT.md | 16 +++++ docs/CHANGELOG.md | 1 + docs/CPU_DATA_ANALYSIS.md | 2 + rugnux/CMakeLists.txt | 2 + rugnux/ModelScaling.cpp | 84 ++++++++++++++++++++++ rugnux/ModelScaling.h | 34 +++++++++ rugnux/ModelValidation.cpp | 4 +- tests/CMakeLists.txt | 1 + tests/ModelScalingTest.cpp | 143 +++++++++++++++++++++++++++++++++++++ 9 files changed, 285 insertions(+), 2 deletions(-) create mode 100644 rugnux/ModelScaling.cpp create mode 100644 rugnux/ModelScaling.h create mode 100644 tests/ModelScalingTest.cpp diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 1d39cfe76..1720df60a 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -318,3 +318,19 @@ synchrotron source by R. Kahn, R. Fourme, A. Gadet, J. Janin, C. Dumas and D. An crystallography with synchrotron radiation: photographic data collection and polarization correction" (1982), J. Appl. Cryst. 15, 330-337 [doi:10.1107/S0021889882012060](https://doi.org/10.1107/S0021889882012060). + +**Bulk-solvent correction and overall scaling** — the model's structure factors are put on the +observed scale with an overall factor, an anisotropic B and a flat bulk-solvent term, the flat-mask +model of A. Fokine and A. Urzhumtsev, "Flat bulk-solvent model: obtaining optimal parameters" +(2002), Acta Cryst. D58, 1387-1392 +[doi:10.1107/S0907444902010284](https://doi.org/10.1107/S0907444902010284), which is also the source +of the starting values and of the range those two parameters are physically meaningful over. The +procedure that fits them — a grid search over that range for the solvent pair, with the overall +scale and the anisotropic B refitted at every grid point — follows P. V. Afonine, +R. W. Grosse-Kunstleve and P. D. Adams, "A robust bulk-solvent correction and anisotropic scaling +procedure" (2005), Acta Cryst. D61, 850-855 +[doi:10.1107/S0907444905007894](https://doi.org/10.1107/S0907444905007894). The fit is unweighted, +as in both that procedure and REFMAC5: G. N. Murshudov, P. Skubak, A. A. Lebedev, N. S. Pannu, +R. A. Steiner, R. A. Nicholls, M. D. Winn, F. Long and A. A. Vagin, "REFMAC5 for the refinement of +macromolecular crystal structures" (2011), Acta Cryst. D67, 355-367 +[doi:10.1107/S0907444911001314](https://doi.org/10.1107/S0907444911001314). diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index f85fea482..c94570e18 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -4,6 +4,7 @@ ### 1.0.0-rc.167 * `rugnux --model` fits the model's scale, anisotropic B and bulk-solvent parameters on the working reflections only, so the R-free it reports is measured against a model no free reflection helped scale. +* The bulk-solvent parameters of `rugnux --model` are searched over their physically meaningful range instead of being fitted without bounds, so a model is never scaled with a solvent term that has silently switched itself off. * The rugnux results report opens with a summary - `VERDICT=` (`OK`, `WARNINGS`, `UNUSABLE`, `FAILED`), `VERDICT_TEXT=`, `PATHOLOGY_FLAGS=` with one closed-vocabulary code per condition that warned, and the `WARNING:` lines, which used to close the file - and the sections after it are renumbered 1-5 with no gaps; `REPORT_VERSION` is 12. * `rugnux --developer` writes the full results report - the pipeline-internal keys and the long explanations the default report now leaves out - and `--finalist-ledger` adds the evidence for every space group the search considered, not only the one it adopted. * The results report warns when the merged data carry no usable signal and when too little of reciprocal space was measured inside the fitted resolution, and omits `FITTED_RESOLUTION` where the CC1/2 curve it is fitted on never falls off. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index e9183a44e..6d66d1dcf 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -57,6 +57,8 @@ The methods draw on, and in places reimplement, solutions from: - A. Thorn & G. M. Sheldrick, "ANODE: anomalous and heavy-atom density calculation", *J. Appl. Cryst.* **44** (2011), 1285-1287 (anomalous difference density read at the model's sites). - R. Kahn, R. Fourme, A. Gadet, J. Janin, C. Dumas & D. Andre, "Macromolecular crystallography with synchrotron radiation: photographic data collection and polarization correction", *J. Appl. Cryst.* **15** (1982), 330-337 (the azimuthal polarization factor of §2.2, applied to both the azimuthal profile and the Bragg intensities). - R. J. Read, "Improved Fourier coefficients for maps using phases from partial structures with errors", *Acta Cryst.* **A42** (1986), 140-149 (the sigma_A formalism and the m, D weighting of the map coefficients of §14.4). +- A. Fokine & A. Urzhumtsev, "Flat bulk-solvent model: obtaining optimal parameters", *Acta Cryst.* **D58** (2002), 1387-1392 (the flat bulk-solvent model, its optimal parameters and the range they are physically meaningful over, used when scaling a model to the data in §14). +- P. V. Afonine, R. W. Grosse-Kunstleve & P. D. Adams, "A robust bulk-solvent correction and anisotropic scaling procedure", *Acta Cryst.* **D61** (2005), 850-855 (the grid search over that range that fits k_sol and b_sol, with the overall scale and anisotropic B refitted at each grid point). - K. Shoemake, "Uniform Random Rotations", in *Graphics Gems III*, ed. D. Kirk, Academic Press (1992), 124-132 (the uniform random rotations the model-fit null of §14.5 is built from). - Z. Otwinowski & W. Minor, "Processing of X-ray diffraction data collected in oscillation mode", *Methods Enzymol.* **276** (1997), 307-326 (reweighted, de-biased profile-fit variances). - G. Winter et al., "DIALS: implementation and evaluation of a new integration package", *Acta Cryst.* **D74** (2018), 85-97, and J. Beilsten-Edmands et al., *Acta Cryst.* **D76** (2020), 385-399 (CC1/2 resolution cutoff, merge outlier rejection, scaling error model). diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index 65f0669c1..3484e35e9 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -8,6 +8,8 @@ ADD_LIBRARY(Rugnux STATIC Rugnux.h RugnuxCommandLine.cpp RugnuxCommandLine.h + ModelScaling.cpp + ModelScaling.h ModelValidation.cpp ModelValidation.h RigidBodyRefine.cpp diff --git a/rugnux/ModelScaling.cpp b/rugnux/ModelScaling.cpp new file mode 100644 index 000000000..600b7df3f --- /dev/null +++ b/rugnux/ModelScaling.cpp @@ -0,0 +1,84 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "ModelScaling.h" + +#include +#include + +namespace { + +// R-factor of the current parameters over the fitted reflections. This is what the grid is +// selected on, and it is the quantity the scale exists to make small. +double RFactor(const gemmi::Scaling &scaling) { + double num = 0, den = 0; + for (const auto &p : scaling.points) { + num += std::fabs(p.fobs - scaling.compute_value(p)); + den += p.fobs; + } + return den > 0 ? num / den : 1.0; +} + +} // namespace + +// Following the phenix bulk-solvent and scaling procedure: k_sol and b_sol by a grid search, with +// the overall scale and the anisotropic B refitted at every grid point - Afonine, Grosse-Kunstleve +// & Adams, Acta Cryst. D61, 850-855, 2005, which searches b_sol over 10-80 A^2 in steps of 5. +// The fit is unweighted, as in both phenix and Refmac (Murshudov, Skubak, Lebedev, Pannu, Steiner, +// Nicholls, Winn, Long & Vagin, Acta Cryst. D67, 355-367, 2011, eq. 11). The physical range and the +// starting values are those of Fokine & Urzhumtsev, Acta Cryst. D58, 1387-1392, 2002. +// +// The point of the grid is that k_sol and b_sol cannot leave the physical box: gemmi's own +// fit_parameters() is an unbounded Levenberg-Marquardt, and on this corpus it reached b_sol of +// 1707 A^2 - a solvent term switched off in all but the lowest-resolution shell. Here the solvent +// pair is held fixed at each grid point and only the overall scale and the symmetry-constrained +// anisotropic B are refined, which is the well-conditioned half of the problem and is left to +// gemmi's solver rather than reimplemented. +ModelScaleReport FitModelScale(gemmi::Scaling &scaling, ModelScaleBox box) { + ModelScaleReport report; + report.n_points = static_cast(scaling.points.size()); + if (scaling.points.size() < 20) + return report; + + const bool had_solvent = scaling.use_solvent; + scaling.fix_k_sol = true; // the grid owns the solvent pair; the solver never sees it + scaling.fix_b_sol = true; + + double best_r = -1, best_k_sol = 0.35, best_b_sol = 46.0, best_k_overall = 1.0; + gemmi::SMat33 best_b_star{0, 0, 0, 0, 0, 0}; + + auto try_point = [&](double k_sol, double b_sol) { + scaling.k_sol = k_sol; + scaling.b_sol = b_sol; + scaling.fit_isotropic_b_approximately(); // a fresh starting point for this solvent pair + scaling.fit_parameters(); // k_overall + anisotropic B only + ++report.n_grid; + const double r = RFactor(scaling); + if (best_r < 0 || r < best_r) { + best_r = r; + best_k_sol = k_sol; + best_b_sol = b_sol; + best_k_overall = scaling.k_overall; + best_b_star = scaling.b_star; + } + }; + + // Coarse pass over the whole box, then one refinement pass around the winner. + for (double ks = box.k_lo; ks <= box.k_hi + 1e-9; ks += 0.05) + for (double bs = box.b_lo; bs <= box.b_hi + 1e-9; bs += 10.0) + try_point(ks, bs); + const double k0 = best_k_sol, b0 = best_b_sol; + const double k_hi2 = std::min(box.k_hi, k0 + 0.05); + const double b_hi2 = std::min(box.b_hi, b0 + 10.0); + for (double ks = std::max(box.k_lo, k0 - 0.05); ks <= k_hi2 + 1e-9; ks += 0.025) + for (double bs = std::max(box.b_lo, b0 - 10.0); bs <= b_hi2 + 1e-9; bs += 5.0) + try_point(ks, bs); + + scaling.k_sol = best_k_sol; + scaling.b_sol = best_b_sol; + scaling.k_overall = best_k_overall; + scaling.b_star = best_b_star; + scaling.use_solvent = had_solvent; + report.r_work_fit = best_r; + return report; +} diff --git a/rugnux/ModelScaling.h b/rugnux/ModelScaling.h new file mode 100644 index 000000000..0600097d1 --- /dev/null +++ b/rugnux/ModelScaling.h @@ -0,0 +1,34 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include "gemmi/scaling.hpp" + +// Fit the overall scale, the anisotropic B and the flat bulk solvent of a model to observed +// amplitudes, writing the answer into `scaling` (k_overall, b_star, k_sol, b_sol) so that +// gemmi's own scale_data() can then apply it. +// +// The caller decides which reflections are fitted: FitModelScale reads scaling.points, which +// prepare_points() filled. Pass the WORKING set to keep the free reflections out of the scale +// model - that is a property of the call site, not of this function. +// +// Replaces gemmi's fit_isotropic_b_approximately() + fit_parameters(), which is an unbounded +// Levenberg-Marquardt: nothing there stops k_sol and b_sol reaching values a flat solvent model +// is meaningless at (b_sol of 1707 A^2 was measured on this corpus). Here k_sol and b_sol come +// from a grid search over the physical ranges, so an unphysical answer cannot be represented. +// The box the solvent parameters are searched in. The defaults are the physically reasonable +// range reported by Fokine & Urzhumtsev and adopted as phenix's default grid: k_sol in (0.1, 0.8) +// and b_sol in (10, 80) A^2, of which phenix searches b_sol 10-80 in steps of 5. +struct ModelScaleBox { + double k_lo = 0.10, k_hi = 0.60; + double b_lo = 10.0, b_hi = 80.0; +}; + +struct ModelScaleReport { + int n_points = 0; // reflections the fit actually used + int n_grid = 0; // grid points evaluated + double r_work_fit = 1.0; // R on the fitted reflections, at the chosen solution +}; + +ModelScaleReport FitModelScale(gemmi::Scaling &scaling, ModelScaleBox box = {}); diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index 542459b14..acdf70127 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -2,6 +2,7 @@ // SPDX-License-Identifier: GPL-3.0-only #include "ModelValidation.h" +#include "ModelScaling.h" #include #include @@ -316,8 +317,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & gemmi::Scaling scaling(ucell, sg); scaling.use_solvent = true; scaling.prepare_points(out.fmodel, fobs_work, &ms.fmask); - scaling.fit_isotropic_b_approximately(); - scaling.fit_parameters(); + FitModelScale(scaling); scaling.scale_data(out.fmodel, &ms.fmask); // out.fmodel now holds the scaled, solvent-corrected Fmodel out.k_sol = scaling.k_sol; out.b_sol = scaling.b_sol; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index c140a7610..cf3140099 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -95,6 +95,7 @@ ADD_EXECUTABLE(jfjoch_test XDSPluginTest.cpp MergeScaleTest.cpp AnisotropyAnalysisTest.cpp + ModelScalingTest.cpp TwinningAnalysisTest.cpp TranslationalNCSTest.cpp RfreeFlagsTest.cpp diff --git a/tests/ModelScalingTest.cpp b/tests/ModelScalingTest.cpp new file mode 100644 index 000000000..155ece16a --- /dev/null +++ b/tests/ModelScalingTest.cpp @@ -0,0 +1,143 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include +#include + +#include "../rugnux/ModelScaling.h" +#include "gemmi/scaling.hpp" +#include "gemmi/symmetry.hpp" +#include "gemmi/unitcell.hpp" + +namespace { + +// A synthetic set of scaling points: |Fobs| is generated from a known scale, isotropic B, k_sol +// and b_sol, so the fit has a right answer to find. The "structure factors" are a hash of the +// index, so the set is reproduced bit for bit. +void MakePoints(gemmi::Scaling &scaling, const gemmi::UnitCell &cell, + double k_overall, double b_iso, double k_sol, double b_sol, + double d_min, double d_max) { + scaling.points.clear(); + const int hmax = static_cast(cell.a / d_min) + 1; + const int kmax = static_cast(cell.b / d_min) + 1; + const int lmax = static_cast(cell.c / d_min) + 1; + uint32_t seed = 12345; + auto rnd = [&seed]() { + seed = seed * 1664525u + 1013904223u; + return static_cast((seed >> 8) & 0xffff) / 65535.0; + }; + for (int h = 0; h <= hmax; ++h) + for (int k = 0; k <= kmax; ++k) + for (int l = 0; l <= lmax; ++l) { + if (h == 0 && k == 0 && l == 0) continue; + const gemmi::Miller hkl{{h, k, l}}; + const double d = cell.calculate_d(hkl); + if (d < d_min || d > d_max) continue; + const double stol2 = cell.calculate_stol_sq(hkl); + const std::complex fc(20.0 + 80.0 * rnd(), 0.0); + const std::complex fm(5.0 + 5.0 * rnd(), 0.0); + const std::complex total = fc + k_sol * std::exp(-b_sol * stol2) * fm; + const double fobs = k_overall * std::exp(-b_iso * stol2) * std::abs(total); + gemmi::Scaling::Point p{}; + p.hkl = hkl; + p.stol2 = stol2; + p.fcmol = std::complex(static_cast(fc.real()), 0.f); + p.fmask = std::complex(static_cast(fm.real()), 0.f); + p.fobs = static_cast(fobs); + p.sigma = 1.f; + scaling.points.push_back(p); + } +} + +gemmi::UnitCell Cell() { + gemmi::UnitCell c; + c.set(60.0, 70.0, 80.0, 90.0, 90.0, 90.0); + return c; +} + +} // namespace + +// The fit has to find a solvent pair it was given, not merely land somewhere legal. +TEST_CASE("ModelScaling recovers a known bulk solvent") { + const gemmi::UnitCell cell = Cell(); + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 21 21 21"); + REQUIRE(sg != nullptr); + gemmi::Scaling scaling(cell, sg); + scaling.use_solvent = true; + MakePoints(scaling, cell, /*k_overall=*/0.5, /*b_iso=*/20.0, + /*k_sol=*/0.35, /*b_sol=*/45.0, /*d_min=*/2.0, /*d_max=*/50.0); + REQUIRE(scaling.points.size() > 500); + + const ModelScaleReport report = FitModelScale(scaling); + CHECK(report.n_grid > 0); + // The grid steps are 0.025 in k_sol and 5 in b_sol, so this is one step of tolerance. + CHECK(scaling.k_sol == Catch::Approx(0.35).margin(0.03)); + CHECK(scaling.b_sol == Catch::Approx(45.0).margin(6.0)); + CHECK(report.r_work_fit < 0.05); +} + +// The regression this exists for: gemmi's own fit_parameters() is an unbounded Levenberg-Marquardt +// and, where the data cannot determine the solvent, walks k_sol/b_sol out of the range a flat +// solvent model means anything in (b_sol of 1707 A^2 was measured on real data). Here the solvent +// contribution is made unidentifiable - only high-resolution data, where exp(-b_sol * s^2) is +// indistinguishable from zero for any large b_sol - and the fit must still return a physical pair. +TEST_CASE("ModelScaling stays physical when the solvent is unidentifiable") { + const gemmi::UnitCell cell = Cell(); + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 21 21 21"); + REQUIRE(sg != nullptr); + gemmi::Scaling scaling(cell, sg); + scaling.use_solvent = true; + // Nothing below 3 A: the bulk solvent is a low-resolution feature, so it has almost + // no leverage here. + MakePoints(scaling, cell, /*k_overall=*/0.5, /*b_iso=*/20.0, + /*k_sol=*/0.35, /*b_sol=*/45.0, /*d_min=*/1.2, /*d_max=*/3.0); + REQUIRE(scaling.points.size() > 500); + + FitModelScale(scaling); + CHECK(scaling.k_sol >= 0.15); + CHECK(scaling.k_sol <= 0.50); + CHECK(scaling.b_sol >= 10.0); + CHECK(scaling.b_sol <= 80.0); + + // ... where gemmi's own unbounded fit is free to leave that range. + gemmi::Scaling unbounded(cell, sg); + unbounded.use_solvent = true; + unbounded.points = scaling.points; + unbounded.fit_isotropic_b_approximately(); + unbounded.fit_parameters(); + CHECK(std::isfinite(unbounded.b_sol)); // it converges; it is simply not bounded +} + +// A densely packed crystal with no disordered solvent channels - a small molecule, say - has an +// empty solvent mask, and then Fmask is zero for every reflection. The bulk-solvent contribution +// k_sol * exp(-b_sol s^2) * Fmask is then identically zero WHATEVER k_sol and b_sol come out as, +// so a solvent term fitted where there is no solvent cannot add anything to the model. That is why +// this needs no solvent-content threshold and no switch: the mask already carries the answer. +TEST_CASE("ModelScaling adds nothing when there is no solvent to model") { + const gemmi::UnitCell cell = Cell(); + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 21 21 21"); + REQUIRE(sg != nullptr); + gemmi::Scaling scaling(cell, sg); + scaling.use_solvent = true; + MakePoints(scaling, cell, /*k_overall=*/0.5, /*b_iso=*/20.0, + /*k_sol=*/0.35, /*b_sol=*/45.0, /*d_min=*/2.0, /*d_max=*/50.0); + REQUIRE(scaling.points.size() > 500); + for (auto &p : scaling.points) // an empty mask: no solvent-accessible volume + p.fmask = {0.f, 0.f}; + + FitModelScale(scaling); + // Whatever the search settled on, the model it produces is the solvent-free one. + double worst = 0; + for (const auto &p : scaling.points) + worst = std::max(worst, static_cast( + std::fabs(std::abs(scaling.get_fcalc(p)) - std::abs(p.fcmol)))); + CHECK(worst == Catch::Approx(0.0).margin(1e-6)); + // and the parameters are still reported inside the physical box, not at some arbitrary value + CHECK(scaling.k_sol >= 0.10); + CHECK(scaling.k_sol <= 0.60); + CHECK(scaling.b_sol >= 10.0); + CHECK(scaling.b_sol <= 80.0); +}