diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index b9ab621ee..ed23eb167 100644 --- a/docs/CHANGELOG.md +++ b/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 diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index 84bbb1d45..1eeef5fa5 100644 --- a/rugnux/ModelValidation.cpp +++ b/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(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 &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 & // 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 lattice_ops{gemmi::Op::identity()}; + for (const gemmi::Op &law : ReindexAmbiguityOperators(cell, ambiguity_sg)) + lattice_ops.push_back(law); + std::vector 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 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(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 diff --git a/rugnux/RigidBodyRefine.cpp b/rugnux/RigidBodyRefine.cpp index e60339c57..ccf9c8137 100644 --- a/rugnux/RigidBodyRefine.cpp +++ b/rugnux/RigidBodyRefine.cpp @@ -302,6 +302,24 @@ void SetModelPositions(gemmi::Model &model, const std::vector & // 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 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(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(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 diff --git a/rugnux/RigidBodyRefine.h b/rugnux/RigidBodyRefine.h index a17f3489b..6c15a3b3d 100644 --- a/rugnux/RigidBodyRefine.h +++ b/rugnux/RigidBodyRefine.h @@ -27,6 +27,12 @@ struct RigidBodyRefineResult { std::vector ModelPositions(const gemmi::Model &model); void SetModelPositions(gemmi::Model &model, const std::vector &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