diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index b8624e08d..d08a4a0da 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -11,6 +11,7 @@ * Rugnux adds beam-stop holder arms that let part of the beam through to the beam-stop mask, without changing the mask of a sweep that has none. * Rugnux's `--model` rigid-body refinement computes the model's density and solvent mask on all threads, with results unchanged. * Rugnux's `--model` rigid-body refinement computes the model's structure factors from one copy of the model and derives the translation part of its Jacobian exactly, so it runs faster; placements can differ slightly. +* `rugnux --model` also draws a null replicate again when its rigid-body placement ends near an orientation equivalent to the model's, not only when it starts near one. ### 1.0.0-rc.173 diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index dbd77f5a1..8088dbdd5 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/docs/CPU_DATA_ANALYSIS_DECISIONS.md @@ -173,6 +173,13 @@ as fitted, or the comparison would be between a placed model and unplaced nulls be inflated by the placement rather than by the model. Measured on one rotation data set: the crystal's own model +15.0σ, an unrelated protein +1.8σ, and the correct model rigidly rotated 90° −1.0σ. +A replicate is a sample of the null only while it stays away from the model's own solution. A draw +within the rigid body's reach of an orientation equivalent to the model's — under the space group's +rotations or a twin law of the lattice — is drawn again before it is placed, and a replicate the +placement nevertheless carries to within that reach is drawn and placed again afterwards, each +replicate from a generator of its own so the result does not depend on the order the concurrent +replicates finish in. + **R-work carries the decision, not R-free.** Nothing is refined against the working set here — the scale has a handful of parameters (a scale, an anisotropic $B$ and two solvent terms) and the placement six — so R-work carries no optimism, and it has an order of magnitude more reflections than R-free, and so that much more power to separate the two arms. R-free diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index 6f2a20852..d16a53328 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -641,6 +641,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & Fit fit; bool rb_applied = false; double rb_angle_deg = 0, rb_shift_A = 0, r_free_before_rb = 0; + gemmi::Mat33 rb_rotation; // the committed rotation about the centroid, identity if none }; auto place_and_fit = [&](ModelState &ms, const std::vector &obs_in, std::optional fitted, double to_d) -> Placement { @@ -673,6 +674,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & if (commit) { out.rb_applied = true; out.rb_angle_deg = rb.angle_deg; + out.rb_rotation = rb.rotation; out.rb_shift_A = rb.shift_A; out.r_free_before_rb = out.fit.r_free; out.fit = std::move(moved); @@ -857,20 +859,41 @@ ModelValidationResult ValidateAgainstModel(const std::vector & // Each replicate works on its own copy of the model and of its structure factors, and the // replicates share nothing else, so they run concurrently and the real model - which the // maps and the atom-density readout below are still taken from - is never moved. + // + // The check before the fit is only a filter: the rigid body can walk further than its nominal + // reach, and a replicate it carries onto an orientation equivalent to the model's has been + // refined into the model's own solution (measured: 24 deg, R-free 0.56 -> 0.39, which took the + // null's spread from 0.016 to 0.066 and a fitting model below the gate). So where the + // replicate ENDS is checked as well, and one that ended within reach is drawn again and + // placed again - from a generator of its own, so the draws a replicate takes do not depend on + // the order the concurrent replicates finish in. std::vector null_r_work(NULL_REPLICATES), null_margin(NULL_REPLICATES); + std::atomic placed_redraws{0}; auto run_replicate = [&](int i) { - ModelState rep; - rep.st = mdl.st; - SetModelPositions(rep.st.models[0], as_read); - reorient_about_centroid(rep.st.models[0], null_rotation[i]); - compute_model_factors(rep, null_d_min); - std::optional identity_fit; - if (!reindex_ops.empty()) { - IndexingProbe probe = probe_indexing(rep, obs_null); - null_margin[i] = probe.margin; - identity_fit = std::move(probe.identity_fit); // replicates are placed against obs + std::mt19937 rng(NULL_SEED + 1 + i); + gemmi::Mat33 rot = null_rotation[i]; + for (int redraw = 0;; redraw++) { + ModelState rep; + rep.st = mdl.st; + SetModelPositions(rep.st.models[0], as_read); + reorient_about_centroid(rep.st.models[0], rot); + compute_model_factors(rep, null_d_min); + std::optional identity_fit; + if (!reindex_ops.empty()) { + IndexingProbe probe = probe_indexing(rep, obs_null); + null_margin[i] = probe.margin; + identity_fit = std::move(probe.identity_fit); // replicates are placed against obs + } + const Placement placed = place_and_fit(rep, obs_null, std::move(identity_fit), null_d_min); + null_r_work[i] = placed.fit.r_work; + const gemmi::Mat33 ended = placed.rb_rotation.multiply(rot); + if (redraw >= NULL_MAX_REDRAWS || angle_to_nearest_deg(ended, equivalent) >= reach_deg) + break; + ++placed_redraws; + rot = random_rotation(rng); + for (int k = 0; k < NULL_MAX_REDRAWS && angle_to_nearest_deg(rot, equivalent) < reach_deg; k++) + rot = random_rotation(rng); } - null_r_work[i] = place_and_fit(rep, obs_null, std::move(identity_fit), null_d_min).fit.r_work; }; // On threads of their own, not on ParallelFor's pool. A parallel pass reached from a pool // worker runs inline (ParallelFor.h), so there the scale fits and the rigid-body Jacobian of @@ -887,6 +910,10 @@ ModelValidationResult ValidateAgainstModel(const std::vector & })); for (std::future &f : replicate_threads) f.get(); + if (placed_redraws > 0) + logger.Info("Model validation: {} random placement(s) were refined to within the rigid body's reach " + "({:.1f} deg) of an orientation equivalent to the model's, and were drawn and placed again", + placed_redraws.load(), reach_deg); const auto [null_mean, null_sd] = mean_sd(null_r_work); result.fit_tested = true; result.null_replicates = NULL_REPLICATES; diff --git a/rugnux/RigidBodyRefine.cpp b/rugnux/RigidBodyRefine.cpp index 165c22f90..3bd814177 100644 --- a/rugnux/RigidBodyRefine.cpp +++ b/rugnux/RigidBodyRefine.cpp @@ -540,6 +540,15 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, result.b_sol = target.b_sol; const double aa = std::sqrt(q[0] * q[0] + q[1] * q[1] + q[2] * q[2]) / target.RmsRadius(); result.angle_deg = aa * 180.0 / PI; + const double axis_angle[3] = {q[0] / target.RmsRadius(), q[1] / target.RmsRadius(), q[2] / target.RmsRadius()}; + double column[3][3]; + for (int j = 0; j < 3; j++) { + const double e[3] = {j == 0 ? 1.0 : 0.0, j == 1 ? 1.0 : 0.0, j == 2 ? 1.0 : 0.0}; + ceres::AngleAxisRotatePoint(axis_angle, e, column[j]); + } + result.rotation = gemmi::Mat33(column[0][0], column[1][0], column[2][0], + column[0][1], column[1][1], column[2][1], + column[0][2], column[1][2], column[2][2]); result.shift_A = std::sqrt(q[3] * q[3] + q[4] * q[4] + q[5] * q[5]); result.seconds = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); diff --git a/rugnux/RigidBodyRefine.h b/rugnux/RigidBodyRefine.h index 313a50cca..13317df97 100644 --- a/rugnux/RigidBodyRefine.h +++ b/rugnux/RigidBodyRefine.h @@ -21,6 +21,7 @@ class Logger; struct RigidBodyRefineResult { bool converged = false; double angle_deg = 0.0; // magnitude of the rotation about the model centroid + gemmi::Mat33 rotation; // that rotation, Cartesian (identity where nothing moved) double shift_A = 0.0; // magnitude of the translation std::vector zones; // the resolution ladder actually walked, coarsest first int evaluations = 0; // structure-factor evaluations (each one re-fits the scale) diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index e34386a1d..531d55dac 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -1140,3 +1140,53 @@ TEST_CASE("ModelValidation_RigidBodyJacobianMatchesNumericScaleRefit", "[ModelVa CheckRigidBodyJacobian(cryst, far, {0, 1, 2, 3, 4, 5}, 0.975, 0.03); } } + +// The null checks where a replicate ENDED against the orientations equivalent to the model's, from the +// rotation the rigid body reports - which must be the rotation it applied: every atom's offset from the +// centroid after the refinement is that rotation of its offset before. +TEST_CASE("ModelValidation_RigidBodyReportsTheRotationItApplied", "[ModelValidation]") { + Logger logger("ModelValidation_RigidBodyReportsTheRotationItApplied"); + const auto path = WriteTemp("rigid_body_rotation_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); + std::filesystem::remove(path); + 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 3 deg about z through the centroid, so there is a rotation to take back. + 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()); + const double a = 3.0 * PI / 180.0; + const gemmi::Mat33 rz(std::cos(a), -std::sin(a), 0, std::sin(a), std::cos(a), 0, 0, 0, 1); + for (gemmi::Position &p : turned) + p = gemmi::Position(rz.multiply(gemmi::Vec3(p) - centre) + centre); + SetModelPositions(st.models[0], turned); + + const RigidBodyRefineResult result = RefineRigidBody(st.models[0], st.cell, *sg, fobs, 3.0, logger); + REQUIRE(result.angle_deg > 1.0); + const std::vector refined = ModelPositions(st.models[0]); + gemmi::Vec3 moved_centre; + for (const gemmi::Position &p : refined) + moved_centre += p; + moved_centre *= 1.0 / static_cast(refined.size()); + double worst = 0; + for (size_t i = 0; i < refined.size(); i++) { + const gemmi::Vec3 expected = result.rotation.multiply(gemmi::Vec3(turned[i]) - centre); + worst = std::max(worst, (gemmi::Vec3(refined[i]) - moved_centre - expected).length()); + } + CHECK(worst < 1e-6); + const double trace = result.rotation[0][0] + result.rotation[1][1] + result.rotation[2][2]; + CHECK(std::acos(std::clamp((trace - 1) / 2, -1.0, 1.0)) * 180 / PI == + Catch::Approx(result.angle_deg).margin(1e-6)); +}