Merge branch 'model-null-twin' into rc173 (model validation null: no random placement within reach of an equivalent orientation)

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C

docs/CHANGELOG.md

@@ -7,6 +7,7 @@
 * jfjoch_viewer's reference dataset accepts a structure-factor mmCIF as well as an MTZ, and a processing job can name an atomic model to validate the merged data against, as `rugnux --model` does.
 * jfjoch_viewer keeps a separate preferred dataset-info plot for grid scans, where "Spots + background" means the spot count.
 * jfjoch_viewer's dark theme no longer leaves navy buttons, red warnings and chart guide lines at their light-theme colours.
+* `rugnux --model` draws its null's random placements away from every orientation equivalent to the model's under the space group or a twin law, so a random placement can no longer be refined onto the model's own solution and make the crystal's own model read as not fitting.
 * 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.
 
 ### 1.0.0-rc.173

rugnux/ModelValidation.cpp

@@ -91,6 +91,11 @@ constexpr int NULL_REPLICATES = 9;
 // Fixed, so the same data and the same model give the same verdict on every run.
 constexpr unsigned NULL_SEED = 20260902;
 
+// How many draws the null may throw away for lying within the rigid body's reach of the model's own
+// orientation. Only a model small enough for that reach to cover most orientations gets near it, and
+// there the remaining draws are taken as they come rather than the run searching forever.
+constexpr int NULL_MAX_REDRAWS = 1000;
+
 // How far above its own null a fit has to sit before the model is allowed to decide anything. The cut
 // is in sigma of that null and not in R: measured, the R a model that explains nothing reaches moves
 // with the model's atom count and B-factors as much as with the data, so no value of R separates the
@@ -122,6 +127,27 @@ gemmi::Mat33 random_rotation(std::mt19937 &rng) {
             2 * (x * z - y * w),     2 * (y * z + x * w),     1 - 2 * (x * x + y * y)};
 }
 
+// A rotation of the lattice, given on fractional coordinates as a symmetry operator is, in Cartesian
+// coordinates - the frame the null's orientations are drawn in.
+gemmi::Mat33 cartesian_rotation(const gemmi::Op &op, const gemmi::UnitCell &cell) {
+    gemmi::Mat33 f;
+    for (int i = 0; i < 3; i++)
+        for (int j = 0; j < 3; j++)
+            f.a[i][j] = static_cast<double>(op.rot[i][j]) / gemmi::Op::DEN;
+    return cell.orth.mat.multiply(f).multiply(cell.frac.mat);
+}
+
+// The angle, in degrees, of the rotation from `rot` to the nearest of `orientations`.
+double angle_to_nearest_deg(const gemmi::Mat33 &rot, const std::vector<gemmi::Mat33> &orientations) {
+    double best = 180.0;
+    for (const gemmi::Mat33 &o : orientations) {
+        const gemmi::Mat33 d = rot.multiply(o.transpose());
+        const double c = std::clamp((d.a[0][0] + d.a[1][1] + d.a[2][2] - 1.0) / 2.0, -1.0, 1.0);
+        best = std::min(best, std::acos(c) * 180.0 / PI);
+    }
+    return best;
+}
+
 // Everything one fit moves: the model itself, and the structure factors that follow it. The real
 // model has one of these and every null replicate gets a copy of its own, which is what lets the
 // replicates run at the same time without stepping on each other.
@@ -740,11 +766,41 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
         // replicates run at the same time below, and a draw taken on whichever thread reached the
         // generator first would make the answer depend on the scheduling. Replicate i gets rotation i
         // on every run, so the sigma reported is a property of the data and not of the machine.
+        //
+        // A replicate has to start from an orientation the model does NOT have, and there is more than
+        // one orientation to keep away from. The crystal looks the same from every orientation its
+        // lattice's rotations take the model to: the space group's own, and those composed with a twin
+        // law, which is where the model sits against data merged in the other indexing - and the
+        // replicates are fitted to the data as merged. A draw within the rigid body's reach of any of
+        // them is walked onto it and scores like the real model: measured on a cubic crystal with a
+        // twin law, one replicate of nine drawn 7.9 deg from the twin-related orientation was placed
+        // to R-work 0.20 against 0.57 for the other eight, and the spread that one gave the null was
+        // enough for the crystal's own model to read as not fitting. Such a draw is not a sample of
+        // the null, so it is drawn again. Everywhere else the draws are exactly the ones taken before.
+        std::vector<gemmi::Op> lattice_ops{gemmi::Op::identity()};
+        for (const gemmi::Op &law : ReindexAmbiguityOperators(cell, ambiguity_sg))
+            lattice_ops.push_back(law);
+        std::vector<gemmi::Mat33> equivalent;
+        for (const gemmi::Op &law : lattice_ops)
+            for (const gemmi::Op &s : gops.sym_ops)
+                equivalent.push_back(cartesian_rotation(s.combine(law), ucell));
+        const double reach_deg = RigidBodyReachDeg(st.models[0], d_min);
         std::vector<gemmi::Mat33> null_rotation;
         {
             std::mt19937 rng(NULL_SEED);
-            for (int i = 0; i < NULL_REPLICATES; i++)
-                null_rotation.push_back(random_rotation(rng));
+            int redraws = 0;
+            while (static_cast<int>(null_rotation.size()) < NULL_REPLICATES) {
+                const gemmi::Mat33 rot = random_rotation(rng);
+                if (angle_to_nearest_deg(rot, equivalent) < reach_deg && redraws < NULL_MAX_REDRAWS) {
+                    redraws++;
+                    continue;
+                }
+                null_rotation.push_back(rot);
+            }
+            if (redraws > 0)
+                logger.Info("Model validation: {} random placement(s) fell within the rigid body's reach "
+                            "({:.1f} deg) of an orientation equivalent to the model's, and were drawn again",
+                            redraws, reach_deg);
         }
         // 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

rugnux/RigidBodyRefine.cpp

@@ -302,6 +302,24 @@ void SetModelPositions(gemmi::Model &model, const std::vector<gemmi::Position> &
 // the movement to recover is the crystal's, not the molecule's. Splitting it into domains or giving a
 // bound ligand its own six parameters would refine against evidence this data does not separately
 // carry, and the ligand is what the difference map is meant to show rather than model away.
+double RigidBodyReachDeg(const gemmi::Model &model, double d_min) {
+    const std::vector<gemmi::Position> pos = ModelPositions(model);
+    if (pos.empty())
+        return 0.0;
+    gemmi::Position centre;
+    for (const gemmi::Position &p : pos)
+        centre += p;
+    centre *= 1.0 / static_cast<double>(pos.size());
+    double r2 = 0;
+    for (const gemmi::Position &p : pos)
+        r2 += centre.dist_sq(p);
+    const double rms_radius = std::sqrt(r2 / static_cast<double>(pos.size()));
+    // The first zone RefineRigidBody walks, which is d_min itself where d_min is coarser than all
+    // of the ladder.
+    const double first_zone = std::max(LADDER[0], d_min);
+    return first_zone / rms_radius * 180.0 / PI;
+}
+
 //
 // Some of the translation would be a gauge rather than a quantity - the origin is free in all three
 // directions in P1 and along the unique axis in a polar group, and |F| does not change when the whole

rugnux/RigidBodyRefine.h

@@ -27,6 +27,12 @@ struct RigidBodyRefineResult {
 std::vector<gemmi::Position> ModelPositions(const gemmi::Model &model);
 void SetModelPositions(gemmi::Model &model, const std::vector<gemmi::Position> &pos);
 
+// How far, in degrees, a rotation about the centroid can be from the answer and still be walked onto
+// it by RefineRigidBody: the angle that moves the model's rms-radius atom by the resolution of the
+// ladder's first zone. Inside it the first zone sees the rotated model overlap the density it belongs
+// to, and its single broad minimum is that placement.
+double RigidBodyReachDeg(const gemmi::Model &model, double d_min);
+
 // 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