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.