From e5c96cd41b49e89c890dccff3bb29eb1d6eb8826 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 7 Sep 2026 15:27:30 +0200 Subject: [PATCH] model validation: the placement target carries a bulk solvent that means something The rigid-body target refits the scale at every evaluation, deliberately - held fixed, the target would measure the scale as much as the placement, and the body would translate to repair a scale error instead of moving to where the density is. That refit was gemmi's unbounded fit, the one already replaced for the reported scale, and here it was worse: measured over a corpus of deposited models, 40% of the evaluations that decide where the model goes came out with a b_sol outside 10-80 A^2, on 57 of 89 datasets, ranging from -8072 to +1721. A negative b_sol is a solvent term that GROWS with resolution. One crystal ran its entire committed placement between 290 and 430 A^2, and that placement went into the reported maps. The bulk solvent is fitted once per zone instead, inside the same physical box the reported fit searches, and then held while the overall scale and the anisotropic B keep following the body. That split is the point: k_sol and b_sol describe the crystal's disordered solvent, not the fit of one placement, and measured across a whole zone they drift by a single grid step. Fitting them at every evaluation costs three times the wall clock, makes the scaler 82% of the run, and moves discontinuously under a forward difference - which is poison for a numerical Jacobian. This costs 4%, leaves no evaluation outside the box, and lands the body within 0.37 degrees of the expensive version, against 4.33 degrees for what it replaces. Placements change on fourteen crystals in eighty-nine. R-free is a wash on the mean; the step buys more total R-free from fewer commits, and the gain sits where the runaway actually bit. Two things found while auditing the file and left as they were, because they are right: the refinement sees working reflections only, end to end, and the gauge that removes the origin-free directions of a polar group carries more than half the whole refined translation in ninety-four of a hundred and thirty-nine polar zones. It had no test. It has one now. An observation with no calculated amplitude gets a zero residual, which drops it from the target rather than scoring it as a perfect fit, and is now counted and reported - a large count says the model's reflection conditions do not match the data's, which is a statement about the model. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N --- docs/CHANGELOG.md | 1 + rugnux/RigidBodyRefine.cpp | 48 +++++++++++++-- rugnux/RigidBodyRefine.h | 6 +- tests/ModelValidationTest.cpp | 108 ++++++++++++++++++++++++++++++++-- 4 files changed, 149 insertions(+), 14 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index afae86fda..c30d50fdc 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -6,6 +6,7 @@ * `rugnux --model` reports CC(model, data) - the correlation of the merged intensities with the placed, scaled model - by resolution shell, on the same shells as CC1/2, with the reflection count and a significance for each. * `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 rigid-body placement of `rugnux --model` uses the same bounded bulk solvent as the reported fit, so a model is no longer placed against a target carrying a solvent term with no physical meaning. * 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. * `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/rugnux/RigidBodyRefine.cpp b/rugnux/RigidBodyRefine.cpp index a2a725c55..e9a34551f 100644 --- a/rugnux/RigidBodyRefine.cpp +++ b/rugnux/RigidBodyRefine.cpp @@ -18,6 +18,7 @@ #include "gemmi/scaling.hpp" // Scaling (bulk solvent + anisotropic B) #include "gemmi/solmask.hpp" // SolventMasker +#include "ModelScaling.h" // FitModelScale #include "../common/JFJochMath.h" // PI #include "../common/Logger.h" @@ -82,11 +83,14 @@ public: for (const auto &hv : fobs_.v) sum += hv.value.value; f_mean_ = fobs_.v.empty() ? 1.0 : sum / static_cast(fobs_.v.size()); + solvent_fitted_ = false; } size_t NumObservations() const { return fobs_.v.size(); } double JacobianStep() const { return JACOBIAN_STEP_FRACTION * d_min_; } int evaluations = 0; + int unmatched = 0; // zone observations with no calculated amplitude to compare against + double k_sol = 0, b_sol = 0; // the bulk solvent the zone's target was evaluated with bool Residuals(const double q[6], double *residuals) { ++evaluations; @@ -116,24 +120,51 @@ public: // Re-fitted at every evaluation: with the scale held at the starting placement's value the // target would measure the scale as much as the placement, and the body would translate to // repair a scale error instead of moving where the density is. + // + // The bulk solvent is not part of that scale. k_sol and b_sol describe the disordered solvent + // of the crystal rather than the fit of one placement, so they are fitted once per zone - by + // FitModelScale, inside the same physical box the reported fit is searched in - and then held + // while the overall scale and the anisotropic B follow the body. Leaving them free at every + // evaluation, which is what gemmi's unbounded Levenberg-Marquardt did here, puts a solvent + // term with no physical meaning inside the target that decides where the model goes: measured + // over a corpus of deposited models, 40% of the evaluations came out with b_sol outside + // 10-80 A^2, some of them negative, which is a solvent that GROWS with resolution. gemmi::Scaling scaling(cell_, &sg_); scaling.use_solvent = true; scaling.prepare_points(fcalc, fobs_, &fmask); if (scaling.points.empty()) return false; + if (!solvent_fitted_) { + FitModelScale(scaling); + k_sol = scaling.k_sol; + b_sol = scaling.b_sol; + solvent_fitted_ = true; + } + scaling.k_sol = k_sol; + scaling.b_sol = b_sol; + scaling.fix_k_sol = true; + scaling.fix_b_sol = true; scaling.fit_isotropic_b_approximately(); scaling.fit_parameters(); scaling.scale_data(fcalc, &fmask); - // Both are sorted and in the same ASU, so one merge pass matches them. + // Both are sorted and in the same ASU, so one merge pass matches them. An observation with no + // calculated amplitude gets residual 0, which drops it from the target rather than scoring it + // as a perfect fit: its Jacobian row below comes out zero as well, and the set that matches is + // fixed by the cell, the group and the zone, so it does not move as the body does. It is + // counted and reported because a large count is a statement about the model rather than about + // this refinement - a group whose reflection conditions the data do not obey leaves half of + // them with nothing to compare against. auto c = fcalc.v.begin(); + unmatched = 0; for (size_t i = 0; i < fobs_.v.size(); ++i) { const gemmi::Miller &h = fobs_.v[i].hkl; while (c != fcalc.v.end() && c->hkl < h) ++c; - residuals[i] = (c != fcalc.v.end() && c->hkl == h) - ? (fobs_.v[i].value.value - std::abs(c->value)) / f_mean_ - : 0.0; + const bool matched = c != fcalc.v.end() && c->hkl == h; + if (!matched) + ++unmatched; + residuals[i] = matched ? (fobs_.v[i].value.value - std::abs(c->value)) / f_mean_ : 0.0; } return true; } @@ -147,6 +178,7 @@ private: gemmi::AsuData> fobs_; double d_min_ = 0; double f_mean_ = 1; + bool solvent_fitted_ = false; }; // Ceres' own numeric differentiation steps by |x| * relative_step_size, which is zero at the start of @@ -299,8 +331,10 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, any_zone_solved = any_zone_solved || summary.IsSolutionUsable(); const gemmi::Vec3 along = gauge.multiply(gemmi::Vec3(q[3], q[4], q[5])); q[3] -= along.x; q[4] -= along.y; q[5] -= along.z; - logger.Debug("Rigid body zone {:.1f} A: {} reflections, {} iterations, {} evaluations, {:.2f} s, " - "rotation {:.3f} deg, translation {:.3f} A", zone, zone_obs.v.size(), + logger.Debug("Rigid body zone {:.1f} A: {} reflections ({} without a calculated amplitude), " + "solvent k_sol {:.2f} b_sol {:.0f} A^2, {} iterations, {} evaluations, {:.2f} s, " + "rotation {:.3f} deg, translation {:.3f} A", zone, zone_obs.v.size(), ev.unmatched, + ev.k_sol, ev.b_sol, summary.iterations.empty() ? 0 : summary.iterations.size() - 1, ev.evaluations - evaluations_before, std::chrono::duration(std::chrono::steady_clock::now() - zone_t0).count(), @@ -311,6 +345,8 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, placement.Apply(q, base, model); // Ceres left the model at a Jacobian probe; put it at the answer result.evaluations = ev.evaluations; result.converged = any_zone_solved; + result.k_sol = ev.k_sol; + result.b_sol = ev.b_sol; const double aa = std::sqrt(q[0] * q[0] + q[1] * q[1] + q[2] * q[2]) / placement.rms_radius; result.angle_deg = aa * 180.0 / PI; result.shift_A = std::sqrt(q[3] * q[3] + q[4] * q[4] + q[5] * q[5]); diff --git a/rugnux/RigidBodyRefine.h b/rugnux/RigidBodyRefine.h index 3bd51283f..2f57fdc5d 100644 --- a/rugnux/RigidBodyRefine.h +++ b/rugnux/RigidBodyRefine.h @@ -19,6 +19,7 @@ struct RigidBodyRefineResult { std::vector zones; // the resolution ladder actually walked, coarsest first int evaluations = 0; // structure-factor evaluations (each one re-fits the scale) double seconds = 0.0; + double k_sol = 0.0, b_sol = 0.0; // the bulk solvent the last zone's target was scored through }; // Every atom position of `model`, in the order the model iterates them. @@ -28,8 +29,9 @@ void SetModelPositions(gemmi::Model &model, const std::vector & // Refine the placement of `model` as one rigid body against the observed amplitudes: an angle-axis // rotation about the model's own centroid followed by a translation, six parameters, over a // coarse-to-fine resolution ladder. Each evaluation recomputes Fcalc and the bulk-solvent mask for -// the moved model and re-fits the scale (k_overall, anisotropic B, k_sol, b_sol), so the target -// measures the placement and not the scale. +// the moved model and re-fits the scale (k_overall and the anisotropic B), so the target measures +// the placement and not the scale. The two flat-solvent constants are fitted once per zone, inside +// their physical range, and then held: they describe the crystal rather than the placement. // // `fobs` should be the working set only - the free reflections are what the caller decides on. The // model is left MOVED (whether or not the refinement helped): the caller scores it and puts it back diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index d638cb287..ed7282c2e 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -53,12 +53,16 @@ ATOM 2 C . CB GLY A 1 12.000 14.000 16.000 1.00 20.00 "ATOM 2 CB GLY A 1 12.000 14.000 16.000 1.00 20.00 C\n" "END\n"; - // A synthetic "protein": carbons filling one asymmetric unit of a small P2(1)2(1)2(1) cell. The - // space group matters - in P1 the origin is free in all three directions, so |F| does not change - // when the whole content is translated and there is no translation to recover. No specimen is - // involved; the positions come out of a fixed seed. - std::string ClusterPdb() { - std::string pdb = "CRYST1 30.000 34.000 38.000 90.00 90.00 90.00 P 21 21 21 4\n"; + // A synthetic "protein": carbons filling one asymmetric unit of a small cell. The space group + // matters - in P1 the origin is free in all three directions, so |F| does not change when the + // whole content is translated and there is no translation to recover. No specimen is involved; + // the positions come out of a fixed seed. + const char *kCryst = "CRYST1 30.000 34.000 38.000 90.00 90.00 90.00 P 21 21 21 4\n"; + // The same cell in a polar group: b is the unique axis, so the origin is free along y alone. + const char *kPolarCryst = "CRYST1 30.000 34.000 38.000 90.00 90.00 90.00 P 1 2 1 2\n"; + + std::string ClusterPdb(const char *cryst = kCryst) { + std::string pdb = cryst; std::mt19937 rng(20260902); std::uniform_real_distribution x(2, 14), y(2, 16), z(2, 18); char line[96]; @@ -185,6 +189,98 @@ TEST_CASE("ModelValidation_RigidBodyRecoversASmallShift", "[ModelValidation]") { std::filesystem::remove(path); } +// A polar space group leaves the origin free along one direction: moving the whole cell content +// along it multiplies every structure factor by a phase and changes no amplitude, so the data +// cannot say where the body sits along it and the refinement must not pretend otherwise. The check +// is closed - the "observed" amplitudes are the model's own - and the model is displaced in all +// three directions at once, so the same run says both what is recovered and what is left alone. +TEST_CASE("ModelValidation_RigidBodyLeavesThePolarDirectionAlone", "[ModelValidation]") { + Logger logger("ModelValidation_RigidBodyLeavesThePolarDirectionAlone"); + + const auto path = WriteTemp("rigid_body_polar_test.pdb", ClusterPdb(kPolarCryst).c_str()); + gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect); + const gemmi::SpaceGroup *sg = st.find_spacegroup(); + REQUIRE(sg != nullptr); + st.setup_cell_images(); + + const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger); + REQUIRE_FALSE(ref.empty()); + gemmi::AsuData> fobs; + fobs.unit_cell_ = st.cell; + fobs.spacegroup_ = sg; + for (const auto &r : ref) + fobs.v.push_back({{{r.h, r.k, r.l}}, {std::sqrt(r.I), 1.0f}}); + fobs.ensure_sorted(); + + // b is the unique axis of P 1 2 1, so y is the free direction and x and z are determined. + const std::vector original = ModelPositions(st.models[0]); + std::vector displaced; + for (const gemmi::Position &p : original) + displaced.emplace_back(p.x + 0.35, p.y + 0.50, p.z - 0.30); + SetModelPositions(st.models[0], displaced); + + const RigidBodyRefineResult result = + RefineRigidBody(st.models[0], st.cell, *sg, fobs, 3.0, logger); + CHECK(result.converged); + + const std::vector refined = ModelPositions(st.models[0]); + gemmi::Vec3 left; + for (size_t i = 0; i < original.size(); i++) + left += gemmi::Vec3(refined[i]) - gemmi::Vec3(original[i]); + left *= 1.0 / static_cast(original.size()); + logger.Info("Rigid-body polar test: left over ({:.3f}, {:.3f}, {:.3f}) A", left.x, left.y, left.z); + CHECK(std::fabs(left.x) < 0.10); + CHECK(std::fabs(left.z) < 0.10); + // Along b nothing was refined away, because there is nothing there to refine. + CHECK(left.y == Catch::Approx(0.50).margin(0.02)); + + std::filesystem::remove(path); +} + +// The bulk solvent the placement is scored through has to stay inside the range a flat solvent model +// means anything in. gemmi's own scaler is an unbounded Levenberg-Marquardt and reaches b_sol of +// hundreds or thousands of A^2, which does not corrupt a reported number here but distorts the target +// that decides where the model goes, at every one of the hundreds of evaluations. +TEST_CASE("ModelValidation_RigidBodySolventStaysPhysical", "[ModelValidation]") { + Logger logger("ModelValidation_RigidBodySolventStaysPhysical"); + + const auto path = WriteTemp("rigid_body_solvent_test.pdb", ClusterPdb().c_str()); + gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect); + const gemmi::SpaceGroup *sg = st.find_spacegroup(); + REQUIRE(sg != nullptr); + st.setup_cell_images(); + + const auto ref = ModelReferenceIntensities(path, {}, {}, 3.0, logger); + REQUIRE_FALSE(ref.empty()); + gemmi::AsuData> fobs; + fobs.unit_cell_ = st.cell; + fobs.spacegroup_ = sg; + for (const auto &r : ref) + fobs.v.push_back({{{r.h, r.k, r.l}}, {std::sqrt(r.I), 1.0f}}); + fobs.ensure_sorted(); + + // Turned right around, so the model explains nothing: that is where an unbounded solvent fit + // has nothing to hold it and runs away. + std::vector turned = ModelPositions(st.models[0]); + gemmi::Vec3 centre; + for (const gemmi::Position &p : turned) + centre += p; + centre *= 1.0 / static_cast(turned.size()); + for (gemmi::Position &p : turned) + p = gemmi::Position(centre - (gemmi::Vec3(p) - centre)); + SetModelPositions(st.models[0], turned); + + const RigidBodyRefineResult result = + RefineRigidBody(st.models[0], st.cell, *sg, fobs, 3.0, logger); + logger.Info("Rigid-body solvent test: k_sol {:.3f}, b_sol {:.1f} A^2", result.k_sol, result.b_sol); + CHECK(result.k_sol >= 0.10); + CHECK(result.k_sol <= 0.60); + CHECK(result.b_sol >= 10.0); + CHECK(result.b_sol <= 80.0); + + std::filesystem::remove(path); +} + // sigma_A is what says how much of the model to believe, so the two ends of its range are what the // weighting has to get right: a model that explains the data completely, and one that explains none // of it.