rugnux --model null: redraw a replicate the rigid body carried onto the model's own orientation

The pre-draw guard keeps a null replicate from STARTING within RigidBodyReachDeg of an orientation
equivalent to the model's (space-group rotations x {identity + twin laws}), but the analytic-Jacobian
rigid body walks further than that reach: on 9yzk one replicate drawn outside the 10.6 deg reach was
walked 23.9 deg onto an equivalent orientation (R-free 0.5578 -> 0.3917), the null SD went
0.016 -> 0.066 and MODEL_FIT flipped ACCEPTED -> REJECTED.

Now where the replicate ENDS is checked too: RefineRigidBody reports the rotation it applied, the
final orientation is that times the drawn one, and a replicate that ended within reach is drawn again
(through the same pre-draw filter) and placed again, capped at NULL_MAX_REDRAWS. The redraws come from
a generator per replicate (NULL_SEED + 1 + i), so they do not depend on the order the concurrent
replicates finish in; the count is logged once per run. The real fit's procedure is unchanged.

Validation, 12 open-arm sets against 20260928-0459_24ae26_r6-pooled: post-placement redraws on 9yzk (1)
and 5epe (1), none elsewhere. 9yzk back to ACCEPTED, +14.15 sigma, null 0.5880 +- 0.0138 (base +12.21,
0.5864 +- 0.0158), p.mtz byte-identical to base. 5epe: MODEL_FIT +41.62 -> +46.39, indexing margin
+40.94 -> +31.97 sigma, same decision. Every other set identical to the previous rigid-bc commit.
-N 1 and the default thread count give identical outputs and logs on 9yzk and 6toc.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C
This commit is contained in:
2026-09-28 12:45:38 +02:00
co-authored by Claude Opus 5.5
parent e0e3e92147
commit 4ce36935b9
6 changed files with 106 additions and 11 deletions
+50
View File
@@ -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<gemmi::ValueSigma<float>> 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<gemmi::Position> turned = ModelPositions(st.models[0]);
gemmi::Vec3 centre;
for (const gemmi::Position &p : turned)
centre += p;
centre *= 1.0 / static_cast<double>(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<gemmi::Position> refined = ModelPositions(st.models[0]);
gemmi::Vec3 moved_centre;
for (const gemmi::Position &p : refined)
moved_centre += p;
moved_centre *= 1.0 / static_cast<double>(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));
}