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 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N

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).

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.

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).

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

rugnux/ModelScaling.cpp

@@ -0,0 +1,84 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "ModelScaling.h"
+
+#include <algorithm>
+#include <cmath>
+
+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<float> &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<float> &scaling, ModelScaleBox box) {
+    ModelScaleReport report;
+    report.n_points = static_cast<int>(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<double> 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;
+}

rugnux/ModelScaling.h

@@ -0,0 +1,34 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// 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<float> &scaling, ModelScaleBox box = {});

rugnux/ModelValidation.cpp

@@ -2,6 +2,7 @@
 // SPDX-License-Identifier: GPL-3.0-only
 
 #include "ModelValidation.h"
+#include "ModelScaling.h"
 
 #include <algorithm>
 #include <cmath>
@@ -316,8 +317,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
         gemmi::Scaling<float> 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;

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

tests/ModelScalingTest.cpp

@@ -0,0 +1,143 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+
+#include <cmath>
+#include <algorithm>
+#include <complex>
+
+#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<float> &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<int>(cell.a / d_min) + 1;
+    const int kmax = static_cast<int>(cell.b / d_min) + 1;
+    const int lmax = static_cast<int>(cell.c / d_min) + 1;
+    uint32_t seed = 12345;
+    auto rnd = [&seed]() {
+        seed = seed * 1664525u + 1013904223u;
+        return static_cast<double>((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<double> fc(20.0 + 80.0 * rnd(), 0.0);
+                const std::complex<double> fm(5.0 + 5.0 * rnd(), 0.0);
+                const std::complex<double> 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<float>::Point p{};
+                p.hkl = hkl;
+                p.stol2 = stol2;
+                p.fcmol = std::complex<float>(static_cast<float>(fc.real()), 0.f);
+                p.fmask = std::complex<float>(static_cast<float>(fm.real()), 0.f);
+                p.fobs = static_cast<float>(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<float> 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<float> 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<float> 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<float> 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<double>(
+                            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);
+}